DESY 21–096

DO–TH 21/21

SAGEX–21–13–E

June 2021

The N3LO Scheme-invariant QCD Evolution of the Non-

singlet Structure Functions $F^{\rm NS}_{2}(x,Q^{2})$ and $g_{1}^{\rm
NS}(x,Q^{2})$

J. Blümlein and M. Saragnese

Deutsches Elektronen–Synchrotron, DESY,
Platanenallee 6, D-15738 Zeuthen, Germany

###### Abstract

We present the scheme-invariant unpolarized and polarized flavor non–singlet
evolution equation to N3LO for the structure functions $F_{2}(x,Q^{2})$ and
$g_{1}(x,Q^{2})$ including the charm- and bottom quark effects in the
asymptotic representation. The corresponding evolution is based on the
experimental measurement of the non–singlet structure functions at a starting
scale $Q_{0}^{2}$. In this way the evolution does only depend on the strong
coupling constant $\alpha_{s}(M_{Z})$ or the QCD scale $\Lambda_{\rm QCD}$ and
the charm and bottom quark masses $m_{c}$ and $m_{b}$ and provides one of the
cleanest ways to measure the strong coupling constant in future high
luminosity deep–inelastic scattering experiments. The yet unknown parts of the
4–loop anomalous dimensions introduce only a marginal error in this analysis.

##  1 Introduction

The measurement of the strong coupling constant $\alpha_{s}(M_{Z})$ from
precision data on deep–inelastic scattering is one of the cleanest ways to
obtain this fundamental parameter of the Standard Model. The present world
data allow measurements of $\delta\alpha_{s}(M_{Z})/\alpha_{s}(M_{Z})$ at the
level of $\sim O(1\%)$, [1, 3, 2]. In the standard analyses one is fitting the
scaling violations of deep–inelastic structure functions and obtains
$\alpha_{s}(M_{Z})$ together with the parameters of the parton distribution
functions, observing the correlation of all parameters, cf. [4]. Compared to
this, the measurement of $\alpha_{s}(M_{Z})$ from $R_{e^{+}e^{-}}$, cf. [5],
does not need to fit a larger amount of further parameters. The different
precision determinations of $\alpha_{s}(M_{Z})$ using different processes do
not yet agree and further precision measurements are therefore needed.

For deep–inelastic scattering a direct determination of $\alpha_{s}(M_{Z})$ is
also possible in the case the starting distributions of the evolution are
measured experimentally. In the flavor non–singlet case the corresponding
non–singlet structure functions at a scale $Q_{0}^{2}$ provide the input and
the data at $Q^{2}>Q_{0}^{2}$ are fitted to a physical evolution function,
which only depends on $\alpha_{s}(M_{Z})$, provided that the input for the
values of the heavy quark masses $m_{c}$ and $m_{b}$ is known at high
precision, as also required in the measurement of $\alpha_{s}(M_{Z})$ from
$R_{e^{+}e^{-}}$. Here $Q^{2}$ denotes the virtuality of the exchanged gauge
boson. Therefore scheme–invariant evolution equations allow for a direct
determination of the strong coupling constant using measured one dimensional
input distributions which are here functions of the Bjorken variable $x$. An
analysis of this kind has been performed in Ref. [6] in the unpolarized case,
effectively up to N3LO, however, only considering massless quarks with some
further assumptions, see also Ref. [7]. The sensitivity of the structure
function on $\alpha_{s}(M_{Z})$ or $\Lambda_{\rm QCD}$ is traditionally
illustrated by the slope $\partial\ln(F_{i}(x,Q^{2}))/\partial\ln(Q^{2})$
determined experimentally [8, 9, 10].

In the present paper we complete the formalism by the single- and double–mass
heavy flavor corrections for the Wilson coefficients up to three–loop order
and provide numerical results for these effects, using the results given in
Refs. [11, 12, 13, 14, 15, 16, 17]111A numeric implementation of the heavy
flavor Wilson coefficients to two–loop order has been given in [18]., as well
as the massless contributions. Scheme–invariant evolution is usually performed
in Mellin $N$ space for technical reasons, because it is simpler in analytic
form compared to the corresponding $z$-space representation. In the polarized
case the reference to factorization scheme–invariant quantities will also
solve the associated $\gamma_{5}$ problem, which usually arises performing the
loop calculation in $D=4+\varepsilon$ space–time dimensions. The strong
coupling constant is used in the $\overline{\sf MS}$ scheme. A main goal of
the present study is to outline the different evolution effects for high
precision measurements.

The measurement of scaling violations concerns that of massless parton
evolution and is there implied by the ultraviolet singularity of the
respective twist–2 local operators. Referring to the scale evolution using
observables necessarily requests to map to the massless partons, which implies
in $N$ space the algebraic use of the Wilson coefficients, already at the hard
scale $Q^{2}$, as an obstacle, which are different for the various hard
processes. For the latter ones it has to be guaranteed that only twist–2 terms
contribute, cf. [20, 19, 7], by putting the cuts $W^{2}>15\leavevmode\nobreak\
\mathrm{GeV}^{2},Q^{2}>10\leavevmode\nobreak\ \mathrm{GeV}^{2}$, where $W$ is
the hadronic mass. The natural scale of evolution in deep–inelastic scattering
is $Q^{2}$.

The evolution equation is solved analytically and one final numerical contour
integral around the singularities of the corresponding problem delivers the
evolved flavor non–singlet structure function in $x$ space. We first consider
the flavor decomposition to provide realistic input distributions and derive
the corresponding evolution operators for the structure function. Here we also
discuss the remainder systematics of the yet not completely known four–loop
non–singlet splitting functions $P^{\rm NS,(3),\pm}$. Finally we are providing
numerical results in the unpolarized and polarized case. Precision analyses of
this kind can be performed in the physics programme at the future collider
experiments at the EIC [21] and LHeC [22].

##  2 Flavor Decomposition

Since the evolution for the three different flavor non–singlet distributions
and the singlet evolution are different, the experimental input at the
starting scale $Q_{0}^{2}$ has to be projected correspondingly by combining
different deep–inelastic structure functions. One may refer to the following
well–known decomposition, see also [23]. Let

$\displaystyle v_{k^{2}-1}^{\pm}$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{k}(q_{l}\pm\bar{q}_{l})-k(q_{k}\pm\bar{q}_{k}),$ (1)

with $q_{i}$ the quark distributions and

$\displaystyle v_{1}^{\pm}$ $\displaystyle=$ $\displaystyle 0$ (2)
$\displaystyle v_{3}^{\pm}$ $\displaystyle=$
$\displaystyle(u\pm\bar{u})-(d\pm\bar{d})$ (3) $\displaystyle v_{8}^{\pm}$
$\displaystyle=$ $\displaystyle(u\pm\bar{u})+(d\pm\bar{d})-2(s\pm\bar{s}),$
(4)

one has

$\displaystyle q_{i}+\bar{q}_{i}$ $\displaystyle=$
$\displaystyle\frac{1}{N_{F}}\Sigma-\frac{1}{i}v^{+}_{i^{2}-1}+\sum_{l=i+1}^{N_{F}}\frac{1}{l(l-1)}v^{+}_{l^{2}-1},$
(5) $\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{N_{F}}(q_{l}+\bar{q}_{l}).$ (6)

The nucleon structure functions for pure photon exchange in the case of three
light flavors $(u,d,s)$ are then given at leading order (LO) by

$\displaystyle F_{2}^{p}$ $\displaystyle=$ $\displaystyle
x\left[\frac{2}{9}\Sigma+\frac{1}{6}v_{3}^{+}+\frac{1}{18}v_{8}^{+}\right]$
(7) $\displaystyle F_{2}^{d}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[F_{2}^{p}+F_{2}^{n}\right]=x\left[\frac{2}{9}\Sigma+\frac{1}{18}v_{8}^{+}\right].$
(8)

A synonymous decomposition applies to the structure functions $xg_{1}^{p}$ and
$xg_{1}^{d}$. To project onto the singlet distribution $\Sigma$ directly one
usually needs charged current structure functions, as measured in neutrino
experiments [24] or at facilities like the planned LHeC project [22] or the
EIC [21] in the unpolarized and polarized case, at high luminosity, with [25]

$\displaystyle\frac{1}{2}\left[W_{2}^{p,+}+W_{2}^{p,-}\right]$
$\displaystyle=$ $\displaystyle x\Sigma$ (9)

above the charm threshold [26]. Here the index $\pm$ denotes the exchange of a
$W^{+}$ or a $W^{-}$ boson, respectively. Otherwise additional information on
sea–quarks is necessary. The flavor non–singlet combinations we are
considering are given by

$\displaystyle F_{2}^{\rm NS}$ $\displaystyle=F_{2}^{p}-F_{2}^{d}$
$\displaystyle=\frac{1}{6}xC_{q}^{\rm NS,+}\otimes v_{3}^{+}$ (10)
$\displaystyle xg_{1}^{\rm NS}$ $\displaystyle=xg_{1}^{p}-xg_{1}^{d}$
$\displaystyle=\frac{1}{6}x\Delta C_{q}^{\rm NS,+}\otimes\Delta v_{3}^{+},$
(11)

where the Wilson coefficients are $C_{q}$ and $\Delta C_{q}$ contain light and
heavy flavor contributions, and $\otimes$ denotes the Mellin convolution

$\displaystyle[A\otimes
B](x)=\int_{0}^{1}dx_{1}\int_{0}^{1}dx_{2}\delta(x-x_{1}x_{2})A(x_{1})B(x_{2}).$
(12)

The Mellin transform

$\displaystyle{\rm\bf M}[f(x)](N)=\int_{0}^{1}dxx^{N-1}f(x)$ (13)

turns Eq. (12) into the product

$\displaystyle{\rm\bf M}[(A\otimes B)(x)](N)={\rm\bf M}[A(x))](N)\cdot{\rm\bf
M}[B(x))](N).$ (14)

Before forming the structure function difference in (10, 11) one has to unfold
the nuclear corrections of the deuteron structure functions. The lowest Mellin
moment of (11) is given by the polarized Bjorken sum rule [27].

##  3 The Non–singlet Evolution

The evolution operator of scheme–invariant flavor non–singlet evolution,
$E_{\rm NS}$, is obtained as follows.222Here and below we will work in Mellin
$N$ space.

The evolution equation for the non–singlet structure functions can be written
as

$\displaystyle\frac{d}{d\ln(Q^{2})}\ln\left[F^{\rm NS}(Q^{2})\right]$
$\displaystyle=$ $\displaystyle\frac{d}{d\ln(Q^{2})}\ln\left[C^{\rm
NS}(Q^{2})\right]+\frac{d}{d\ln(Q^{2})}\ln\left[q^{\rm NS}(Q^{2})\right].$
(15)

Its solution is given by

$\displaystyle F^{\rm NS}(Q^{2})$ $\displaystyle=$ $\displaystyle E_{\rm
NS}(Q^{2},Q_{0}^{2})\cdot F^{\rm NS}(Q^{2}_{0}).$ (16)

The Wilson coefficient is given by

$\displaystyle C(Q^{2})$ $\displaystyle=$ $\displaystyle
1+\sum_{k=1}^{\infty}a_{s}^{k}(Q^{2})C_{k},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ C_{k}=c_{k}+h_{k}\left(L_{c},L_{b}\right).$ (17)

Here $c_{k}$ denote the expansion coefficients of the massless Wilson
coefficients and $h_{k}$ of the massive Wilson coefficient, with

$\displaystyle
L_{c}=\ln\left(\frac{Q^{2}}{m_{c}^{2}}\right),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
L_{b}=\ln\left(\frac{Q^{2}}{m_{b}^{2}}\right)$ (18)

and $m_{c,b}$ are the on–shell charm and bottom quark masses.

In the non–singlet case the heavy flavor corrections contribute form
$O(a_{s}^{2})$ onward. One has

$\displaystyle h_{1}$ $\displaystyle=$ $\displaystyle 0$ (19) $\displaystyle
h_{2}$ $\displaystyle=$ $\displaystyle\hat{h}_{2}(L_{c})+\hat{h}_{2}(L_{b})$
(20) $\displaystyle h_{3}$ $\displaystyle=$
$\displaystyle\hat{h}_{3}(L_{c})+\hat{h}_{3}(L_{b})+\hat{\hat{h}}_{3}(L_{c},L_{b}),$
(21)

where $\hat{h}_{i}$ denote the single mass and $\hat{\hat{h}}_{3}$ the double
mass contributions.

One may rewrite the differential operator

$\displaystyle\frac{d}{d\ln(Q^{2})}=\frac{da_{s}(Q^{2})}{d\ln(Q^{2})}\cdot\frac{d}{da_{s}(Q^{2})}$
(22)

with

$\displaystyle\frac{da_{s}}{d\ln(Q^{2})}=-\sum_{k=0}^{\infty}\beta_{k}a_{s}^{k+2}.$
(23)

The evolution equation for the non–singlet quark density is given by

$\displaystyle\frac{d}{da_{s}}\ln\left[q^{\rm
NS}(Q^{2})\right]=-\frac{1}{2}\frac{\sum_{k=0}^{\infty}a_{s}^{k+1}P_{k}}{\sum_{k=0}^{\infty}\beta_{k}a_{s}^{k+2}},$
(24)

where $\beta_{k}$ are expansion coefficients of the QCD–$\beta$ function and
$P_{k}$ are the splitting functions. The anomalous dimensions are related to
the splitting functions333Our normalizations are such that a factor of two has
to be applied to those given in [28, 29]. by

$\displaystyle\gamma_{ij}^{(k)}(N)=-\int_{0}^{1}dxx^{N-1}P_{ij}^{(k)}(x).$
(25)

The solution of Eq. (15) to N3LO reads

$\displaystyle E_{\rm NS}(Q^{2},Q_{0}^{2})$ $\displaystyle=$
$\displaystyle\left(\frac{a}{a_{0}}\right)^{-\frac{P_{0}}{2\beta_{0}}}\Biggl{\\{}1+\frac{a-a_{0}}{2\beta_{0}^{2}}\biggl{\\{}\Bigl{[}1+a^{2}C_{2}(Q^{2})-a_{0}^{2}C_{2}(Q_{0}^{2})\Bigr{]}\bigl{(}2\beta_{0}^{2}C_{1}-\beta_{0}P_{1}+\beta_{1}P_{0}\bigr{)}$
(26)
$\displaystyle-\frac{\bigl{(}a^{2}-a_{0}^{2}\bigr{)}}{4\beta_{0}^{3}}\bigl{(}2\beta_{0}^{2}C_{1}-\beta_{0}P_{1}+\beta_{1}P_{0}\bigr{)}\Bigl{[}2\beta_{0}^{3}C_{1}^{2}+\beta_{0}^{2}P_{2}-\beta_{0}\beta_{1}P_{1}+\bigl{(}\beta_{1}^{2}-\beta_{0}\beta_{2}\bigr{)}P_{0}\Bigr{]}$
$\displaystyle+\frac{\bigl{(}a^{2}+aa_{0}+a_{0}^{2}\bigr{)}}{3\beta_{0}^{2}}\Bigl{[}2\beta_{0}^{4}C_{1}^{3}-\beta_{0}^{3}P_{3}+\beta_{0}^{2}\beta_{1}P_{2}+\bigl{(}\beta_{0}^{2}\beta_{2}-\beta_{0}\beta_{1}^{2}\bigr{)}P_{1}$
$\displaystyle+\bigl{(}\beta_{0}^{2}\beta_{3}-2\beta_{0}\beta_{1}\beta_{2}+\beta_{1}^{3}\bigr{)}P_{0}\Bigr{]}+\frac{a-a_{0}}{4\beta_{0}^{2}}\bigl{(}2\beta_{0}^{2}C_{1}-\beta_{0}P_{1}+\beta_{1}P_{0}\bigr{)}^{2}$
$\displaystyle+\frac{(a-a_{0})^{2}}{24\beta_{0}^{4}}\bigl{(}2\beta_{0}^{2}C_{1}-\beta_{0}P_{1}+\beta_{1}P_{0}\bigr{)}^{3}-\frac{a+a_{0}}{2\beta_{0}}\Bigl{[}2\beta_{0}^{3}C_{1}^{2}+\beta_{0}^{2}P_{2}-\beta_{0}\beta_{1}P_{1}$
$\displaystyle+P_{0}\bigl{(}\beta_{1}^{2}-\beta_{0}\beta_{2}\bigr{)}\Bigr{]}\biggr{\\}}+a^{2}C_{2}(Q^{2})-a_{0}^{2}C_{2}(Q_{0}^{2})-C_{1}\Bigl{[}a^{3}C_{2}(Q^{2})-a_{0}^{3}C_{2}(Q_{0}^{2})\Bigr{]}$
$\displaystyle+a^{3}C_{3}(Q^{2})-a_{0}^{3}C_{3}(Q_{0}^{2})\Biggr{\\}}$

and $a=a_{s}(Q^{2}),\leavevmode\nobreak\ a_{0}=a_{s}(Q^{2}_{0})$ and
$P_{i}(N)\equiv P_{i}$ denotes the Mellin transform of $P_{i}(z)$. The heavy
quark contributions to the Wilson coefficients are given by [30, 15, 16, 31]

$\displaystyle\hat{h}_{2}^{(Q)}$ $\displaystyle=$
$\displaystyle-\frac{\beta_{0,Q}}{4}P_{qq}^{(0)}\ln^{2}\left(\frac{Q^{2}}{m^{2}}\right)+\frac{1}{2}\hat{P}_{qq}^{(1),\rm
NS}\ln\left(\frac{Q^{2}}{m^{2}}\right)+a_{qq}^{(2),\rm
NS}+\frac{\beta_{0,Q}}{4}\zeta_{2}P_{qq}^{(0)}+\hat{C_{q}}^{(2),\rm NS}$ (27)
$\displaystyle\hat{h}_{3}^{(Q)}$ $\displaystyle=$
$\displaystyle-\frac{1}{6}P_{qq}^{(0)}\beta_{0,Q}\left(\beta_{0}+2\beta_{0,Q}\right)\ln^{3}\left(\frac{Q^{2}}{m^{2}}\right)+\frac{1}{4}\Biggl{[}-2P_{qq}^{(1),\rm
NS}\beta_{0,Q}+2\hat{P}_{qq}^{(1),\rm NS}\left(\beta_{0}+\beta_{0,Q}\right)$
(28)
$\displaystyle-\beta_{1,Q}P_{qq}^{(0)}\Biggr{]}\ln^{2}\left(\frac{Q^{2}}{m^{2}}\right)-\frac{1}{2}\Biggl{[}-\hat{P}_{qq}^{(2),\rm
NS}-\left(4a_{qq,Q}^{(2),\rm
NS}+\zeta_{2}\beta_{0,Q}P_{qq}^{(0)}\right)\left(\beta_{0}+\beta_{0,Q}\right)$
$\displaystyle-
P_{qq}^{(0)}\beta_{1,Q}^{(1)}\Biggr{]}\ln\left(\frac{Q^{2}}{m^{2}}\right)+4\bar{a}_{qq,Q}^{(2),\rm
NS}(\beta_{0}+\beta_{0,Q})+P_{qq}^{(0)}\beta_{1,Q}^{(2)}+\frac{1}{6}P_{qq}^{(0)}\beta_{0}\beta_{0,Q}\zeta_{3}$
$\displaystyle+\frac{1}{4}P_{qq}^{(1),\rm NS}\beta_{0,Q}\zeta_{2}-2\delta
m_{1}^{(1)}\beta_{0,Q}P_{qq}^{(0)}-\delta m_{1}^{(0)}\hat{P}_{qq}^{(1),\rm
NS}+2\delta m_{1}^{(-1)}a_{qq,Q}^{(2),\rm NS}+a_{qq,Q}^{(3),\rm NS}$
$\displaystyle+\Biggl{[}-\frac{\beta_{0,Q}}{4}P_{qq}^{(0)}\ln^{2}\left(\frac{Q^{2}}{m^{2}}\right)+\frac{1}{2}\hat{P}_{qq}^{(1),\rm
NS}\ln\left(\frac{Q^{2}}{m^{2}}\right)+a_{qq}^{(2),\rm
NS}+\frac{\beta_{0,Q}}{4}\zeta_{2}P_{qq}^{(0)}\Biggr{]}C_{q}^{(1),\rm NS}$
$\displaystyle+\hat{C}_{q}^{(3),\rm NS}.$

The two–mass three–loop contributions [32] read

$\displaystyle\hat{\hat{h}}_{3}^{\rm NS}$ $\displaystyle=$ $\displaystyle
P_{qq}^{(0)}\beta_{0,Q}^{2}\left[\frac{2}{3}\left(L_{c}^{3}+L_{b}^{3}\right)+\frac{1}{2}\left(L_{c}^{2}L_{b}+L_{c}L_{b}^{2}\right)\right]-\beta_{0,Q}\hat{P}_{qq}^{(1),\rm
NS}\left(L_{c}^{2}+L_{b}^{2}\right)-\left[4a_{qq,Q}^{(2),\rm
NS}\beta_{0,Q}\right.$ (29)
$\displaystyle\left.-\frac{1}{2}\beta_{0,Q}^{2}P_{qq}^{(0)}\zeta_{2}\right](L_{c}+L_{b})+8\bar{a}_{qq,Q}^{(2),\rm
NS}\beta_{0,Q}+\tilde{a}^{(3),\rm NS}_{qq,Q}(m_{c},m_{b},Q^{2}).$

The two–mass term is the same in the unpolarized and polarized case. In the
r.h.s. of Eqs. (27–29) we define

$\displaystyle\hat{f}(x,N_{F})=f(x,N_{F}+1)-f(x,N_{F}).$ (30)

In the $\overline{\sf MS}$ scheme the iterative solution for $a_{s}(Q^{2})$ is
[33]

$\displaystyle a_{s}(Q^{2})$ $\displaystyle=$
$\displaystyle\frac{1}{\beta_{0}L}-\frac{\beta_{1}}{\beta_{0}^{3}L^{2}}\ln(L)+\frac{1}{\beta_{0}^{3}L^{3}}\Biggl{[}\frac{\beta_{1}^{2}}{\beta_{0}^{2}}(\ln^{2}(L)-\ln(L)-1)+\frac{\beta_{2}}{\beta_{0}}\Biggr{]}$
(31)
$\displaystyle+\frac{1}{\beta_{0}^{4}L^{4}}\Biggl{[}\frac{\beta_{1}^{3}}{\beta_{0}^{3}}\Biggl{(}-\ln^{3}(L)+\frac{5}{2}\ln^{2}(L)+2\ln(L)-\frac{1}{2}\Biggr{)}-3\frac{\beta_{1}\beta_{2}}{\beta_{0}^{2}}\ln(L)+\frac{\beta_{3}}{2\beta_{0}}\Biggr{]},$

with $L=\ln(Q^{2}/\Lambda_{\rm QCD}^{2})$. Here the integration constant for
solving (23) is chosen by $(\beta_{1}/\beta_{0}^{2})\ln(\beta_{0})$ [34]. The
expansion coefficients of the $\beta$-function to N3LO were calculated in
[35]. The flavor matching conditions were given in [33]. The expansion
coefficients of the renormalized mass were given in [37, 36]. The constant and
$O(\varepsilon)$ parts of the massive unrenormalized operator matrix elements
at $O(a_{s}^{k})$ are denoted by $a_{ij}^{(k)}$ and $\bar{a}_{ij}^{(k)}$,
respectively, cf. [11, 12, 38, 39].

The three–loop massless Wilson coefficients are expressed by effective
representations. Otherwise we use the analytic Mellin space representations.
After algebraic reduction [40], they depend on 32 harmonic sums [41, 42] up to
weight w=6 only, which is particular to the flavor non–singlet case to
three–loop orders. Other massive Wilson coefficients have a more involved
structure [43, 44, 46, 45]. It is useful to represent harmonic sums with an
alternating index by their Mellin transform, to eliminate spurious $(-1)^{N}$
terms [42] before performing the analytic continuation from even or odd
integers to the complex plane. The singularities of the problem are located at
the non–positive integers, around which the contour integral is performed, cf.
[47, 48, 49, 50]. For the analytic continuation we follow Refs. [52, 51].

In the case of the structure function $F_{2}$ the non–singlet anomalous
dimensions are $P^{{\rm NS},+}$ and the expansion coefficients of the Wilson
coefficient up to $c_{3}$ were calculated in [53]. In the case of the
structure function $g_{1}$ the anomalous dimensions are $P^{{\rm NS},-}$, cf.
[28], and the Wilson coefficient has been given in [54].444Lower orders of the
non–singlet anomalous dimensions and the massless non–singlet Wilson
coefficients have been calculated in [55, 56, 57, 58, 59, 60, 61, 62, 63, 15,
64] and [65, 66, 67, 68, 63] and in the polarized case in [69, 70, 71, 72, 73,
74] and [76, 77, 78, 75]. Below we will also use the combination

$\displaystyle F_{2}^{\rm h}(N,Q^{2})=\left[E_{\rm NS}-\left.E_{\rm
NS}\right|_{\rm h=0}\right]F_{2}(N,Q_{0}^{2}).$ (32)

The four-loop non–singlet anomalous dimensions are not yet completely known as
a function of $N$. However, a series of moments has been calculated in [79,
80, 81, 82, 83, 84]. We follow the earlier investigation in Ref. [6] and
compare these moments with the Padé-approximation

$\displaystyle P_{qq}^{3,\pm,\rm NS}(N)\approx\frac{P_{qq}^{2,\pm,\rm
NS}(N)^{2}}{P_{qq}^{1,\pm,\rm NS}(N)},$ (33)

with the exact moments in [79, 80, 81, 82, 83, 84] in Table 1. Furthermore,
the leading $N_{F}$ terms for the even moments have been predicted in [85].

From the 2nd moment, which agrees within $21\%$, cf. also [6], the accuracy
improves to 2.2 % for the known even moments at $N=16$ and to 2.6% for the odd
moments at $N=15$. For the first moment for $\gamma^{-,\rm NS}$ the Padé-
approximation delivers even the correct result, after using the l’Hospital
rule.

$N$ | $\delta\gamma^{+,\rm NS}$ | $N$ | $\delta\gamma^{-,\rm NS}$
---|---|---|---
2 | 0.208822541 | 1 | 0.0
4 | 0.123728742 | 3 | 0.147102092
6 | 0.087155544 | 5 | 0.101634935
8 | 0.064949195 | 7 | 0.074593595
10 | 0.049680399 | 9 | 0.056598595
12 | 0.038394815 | 11 | 0.043633919
14 | 0.029638565 | 13 | 0.033767853
16 | 0.022602035 | 15 | 0.025956941
Table 1:  The relative error comparing the exact moments of the four–loop
anomalous dimensions, $\gamma^{(3),\pm,\rm NS}$, with the Padé approximation
(33).

The leading small $x$ terms for $P_{3}^{\rm NS,+}$ and $P_{2}^{\rm NS,-}$ have
been given in [86] after correcting some misprints there in Ref. [87], see
also [88]. However, yet unknown sub–leading terms, having been studied at the
known lower orders, are numerically dominant over the leading small $x$
corrections and several of these sub–leading terms are needed, cf. [87, 49],
to be calculated to obtain reliable quantitative prediction in this region of
$x$. Table 1 shows that (33) provides an excellent model. The earlier
assumption of an error of 100% in [6] on this relation has been very
conservative and still led to an error in $\Lambda_{\rm QCD}$ of
$2\leavevmode\nobreak\ {\rm MeV}$ only, which could now be improved in
principle. Yet the experimental accuracy at this level cannot be reached at
present, since the current experimental error amounts to $\delta\Lambda_{\rm
QCD}=26\leavevmode\nobreak\ {\rm MeV}$ [6].

##  4 Numerical Results

The measured input distributions contain correlated errors, which are
parameterized by the respective fits. cf. [6, 89, 19]. Their evolution,
including the corresponding correlation matrix, has to be performed to provide
the error prediction of the structure function at every fixed value of
$\Lambda_{\rm QCD},m_{c}$ and $m_{b}$ respectively. Here one may use the
corresponding formulae for the correlated error propagation given in Ref. [6]
and extensions of them, which are straightforward. Within future data analyses
one will probably import the values of the heavy quark masses from the world
data analyses. For the charm quark mass it has already been shown that its
value obtained from deep–inelastic scattering data fully agrees with other
precision measurements [90]. In our illustrations we will use the values
$m_{c}=1.59\leavevmode\nobreak\ \mathrm{GeV}$ [90] and
$m_{b}=4.78\leavevmode\nobreak\ \mathrm{GeV}$ [91].

For the input distribution in the unpolarized case we refer to the one of Ref.
[6]

$\displaystyle F_{2}^{\rm NS}(x,Q_{0}^{2})$ $\displaystyle=$ $\displaystyle
C^{\rm NS}(x,Q_{0}^{2})\otimes xq^{\rm NS}(x,Q_{0}^{2})$ (34) $\displaystyle
xq^{\rm NS}(x,Q_{0}^{2})$ $\displaystyle=$
$\displaystyle\frac{1}{3}\biggl{[}0.262\
x^{0.298}(1-x)^{4.032}(1+6.042\sqrt{x}+35.49x)$ (35)
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ -1.085\
x^{0.5}(1-x)^{5.921}(1-3.618\sqrt{x}+16.41x)\biggr{]}$

at $Q_{0}^{2}=4\leavevmode\nobreak\ \mathrm{GeV}^{2}$. In the polarized case
we use the fit to $xg_{1}^{\rm NS}(x,Q_{0}^{2})$ of the structure function of
Ref. [89] at $Q_{0}^{2}=10\leavevmode\nobreak\ \mathrm{GeV}^{2}$,

$\displaystyle xg_{1}^{\rm NS}(x,Q_{0}^{2})$ $\displaystyle=$ $\displaystyle
2.4312\cdot 10^{-5}x^{-0.38573}(1-x)^{2.69522}$ (36)
$\displaystyle\times(-1+279.624\sqrt{x}+1239.76x+7053.93x^{3/2}+1866.208x^{2})$

Figure 1: Left: The structure function $F_{2}^{\rm NS}$ at N3LO. Right: The
structure function $xg_{1}^{\rm NS}$ at N3LO. Full lines:
$Q^{2}=10\leavevmode\nobreak\ \mathrm{GeV}^{2}$; dashed lines:
$100\leavevmode\nobreak\ \mathrm{GeV}^{2}$; dash-dotted lines:
$1000\leavevmode\nobreak\ \mathrm{GeV}^{2}$; dotted lines:
$10000\leavevmode\nobreak\ \mathrm{GeV}^{2}$.

In the numeric illustration our reference starting scale will be chosen to be
$Q_{0}^{2}=10\leavevmode\nobreak\ \mathrm{GeV}^{2}$ both in the unpolarized
and polarized case. In Figures 1 we show the scheme–invariant evolution of the
non–singlet structure functions $F_{2}^{\rm NS}$ and $xg_{1}^{\rm NS}$ in the
kinematic region $Q^{2}\in[10,10^{4}]\leavevmode\nobreak\ \mathrm{GeV}^{2}$.
In Figures 2 we expand the representation for the region of larger values of
$x$. Both structure functions show a falling behaviour both towards small and
large values of $x$. To improve the measurement of the strong coupling
constant future measurements of these structure functions are necessary at the
1% level.

Figure 2: The same as Figure 1 but with expanded large $x$ region.

Figure 3: Left: The relative contribution of $F_{2}^{\rm NS}$ in the evolution
from $Q^{2}=10\leavevmode\nobreak\ \mathrm{GeV}^{2}$ to
$10000\leavevmode\nobreak\ \mathrm{GeV}^{2}$. Right: The same for the
structure function $xg_{1}^{\rm NS}$.

In Figures 3 we illustrate the relative effect of the scale evolution in
$Q^{2}$ both for $F_{2}^{\rm NS}$ and $xg_{1}^{\rm NS}$ comparing to the
starting scale $Q_{0}^{2}$. There is a fixpoint at $x_{0}\sim 0.1$ with a
positive evolution with $Q^{2}$ at lower $x$ and a negative evolution at
larger $x$.

Figure 4: Left: The relative contributions from LO (dotted lines), NLO (dashed
lines) and NNLO (full lines) to the structure function $F_{2}^{\rm NS}$ at
N3LO at $Q^{2}=100\leavevmode\nobreak\ \mathrm{GeV}^{2}$ as an example. Right:
The same for the structure function $xg_{1}^{\rm NS}$.

In Figures 4 we show the ratio of the results obtained at leading order (LO),
next-to-leading order (NLO), and next-to-next-to leading order (NNLO) to the
N3LO results at $Q^{2}=100\leavevmode\nobreak\ \mathrm{GeV}^{2}$. This
illustrates the convergence of the perturbative series of the corrections and
the necessity to include N3LO corrections at an accuracy of the data in the 1%
region.

Figure 5: Left: The relative contribution of the heavy flavor contributions
due to $c$ and $b$ quarks to the structure function $F_{2}^{\rm NS}$ at N3LO;
dashed lines: $100\leavevmode\nobreak\ \mathrm{GeV}^{2}$; dashed-dotted lines:
$1000\leavevmode\nobreak\ \mathrm{GeV}^{2}$; dotted lines:
$10000\leavevmode\nobreak\ \mathrm{GeV}^{2}$. Right: The same for the
structure function $xg_{1}^{\rm NS}$ at N3LO.

In Figures 5 we illustrate the relative size of the heavy flavor parts for the
same region in $Q^{2}$ in the unpolarized and polarized cases. In the
important region $x\leq 0.4$ the heavy flavor corrections reach the size of
$\sim 1\%$. Therefore they are of importance for precision analyses. Larger
relative effects are found at larger values of $x$. There, however, the
structure functions and their heavy flavor content drop rapidly.

In Figures 6 we illustrate the effect of the half difference if putting
$P_{qq}^{3,\pm,\rm NS}=2{P_{qq}^{2,\pm,\rm NS}}^{2}/P_{qq}^{2,\pm,\rm NS}$ and
$P_{qq}^{3,\pm,\rm NS}=0$ for both $F_{2}^{\rm NS}$ and $xg_{1}^{\rm NS}$.
This rescaled correction is in the sub–percent range. Moreover, the impact on
$\Lambda_{\rm QCD}$ comes from the slope in $Q^{2}$ which is seen to be rather
small.

Figure 6: The effect of the variation of $P_{qq}^{3,\pm,\rm NS}$ around the
value in Eq. (33) by $\pm 100\%$. Dashed lines: $Q^{2}=100\leavevmode\nobreak\
\mathrm{GeV}^{2}$, Dash-dotted lines: $Q^{2}=1000\leavevmode\nobreak\
\mathrm{GeV}^{2}$; Dotted lines: $Q^{2}=10000\leavevmode\nobreak\
\mathrm{GeV}^{2}$. Left: $F_{2}^{\rm NS}$; Right: $xg_{1}^{\rm NS}$.

##  5 Conclusions

We have calculated the evolution of flavor non–singlet structure functions in
the unpolarized and polarized case, which become available in high luminosity
experiments performed both on proton and deuteron data with comparable
statistics. The evolution has been performed in Mellin $N$ space using our
code QCDEVO. Referring to the structure functions directly has the advantage
that the input distributions are measured. The evolution does only depend on
the QCD scale $\Lambda_{\rm QCD}$ and the heavy quark masses $m_{c}$ and
$m_{b}$, the latter of which have been measured very precisely already [91].
Therefore, this method allows for a one parameter fit of $\Lambda_{\rm QCD}$
only, and provides a method which is systematically very stable. The input,
after corrections for the deuteron wave function, is given by the difference
of two structure functions which both can be measured at very high statistics
at future facilities [21, 22]. The heavy flavor corrections amount to a
contribution of $\sim 1\%$ in the region $x<0.4$.

Let us finally mention that one may also consider the scheme–invariant singlet
evolution [92, 93, 62, 65, 94, 95, 99, 96, 97, 98, 100] in addition to the one
in the flavor non–singlet case, provided the corresponding flavor
decomposition can be carried out. This will usually require more than just the
deep–inelastic data, because of the sea–quark combinations. A possible way
would consist in referring to the Drell–Yan cross section here. Already this
part complicates the analysis. Furthermore, one needs to measure the slope
$\partial F_{i}/\partial\ln(Q^{2})$ at the input scale $Q_{0}^{2}$ together
with $F_{i}(Q_{0}^{2})$ in an uncorrelated way. A second possibility would
consist in using the pair $F_{2}^{\rm S}(Q_{0}^{2}),F_{L}^{\rm S}(Q_{0}^{2})$.
However, the measurement of the structure function $F_{L}^{\rm S}$ is much
more difficult than that of $F_{2}^{\rm S}$ and in the past not enough
statistics has been collected in this case, cf. [101, 102, 103], see, however,
[104]. To perform the necessary flavor decomposition, one needs to perform all
these measurements both at proton and deuteron targets. By scheme–invariant
evolution the fitting problem of input distributions does not exist and no
further optimization by neural network techniques is necessary since the input
distribution is fully determined by its measurement.

Acknowledgment. We would like to thank K. Chetyrkin, K. Schönwald, and A. Vogt
for discussions. This project has received funding from the European Union’s
Horizon 2020 research and innovation programme under the Marie Skłodowska-
Curie grant agreement No. 764850, SAGEX.

## References

  * [1] S. Bethke et al., Workshop on Precision Measurements of $\alpha_{s}$, arXiv:1110.0016 [hep-ph].
  * [2] S. Moch et al., High precision fundamental constants at the TeV scale, arXiv:1405.4781 [hep-ph].
  * [3] S. Alekhin, J. Blümlein and S. O. Moch, Mod. Phys. Lett. A 31 (2016) no.25, 1630023.
  * [4] A. Accardi et al., Eur. Phys. J. C 76 (2016) no.8, 471 [arXiv:1603.08906 [hep-ph]];
S. Alekhin, J. Blümlein, S. Moch and R. Placakyte, Phys. Rev. D 96 (2017)
no.1, 014011 [arXiv:1701.05838 [hep-ph]].

  * [5] P.A. Baikov, K.G. Chetyrkin, J.H. Kühn and J. Rittinger, Phys. Rev. Lett. 108 (2012) 222003 [arXiv: 1201.5804 [hep-ph]].
  * [6] J. Blümlein, H. Böttcher and A. Guffanti, Nucl. Phys. B 774 (2007) 182–207 [hep-ph/0607200].
  * [7] J. Blümlein and H. Böttcher, Proc. DIS 2012, Bonn, 237–241 (DESY, Hamburg, 2012) [arXiv:1207.3170 [hep-ph]].
  * [8] F. Eisele, Rept. Prog. Phys. 49 (1986) 233–339.
  * [9] A.C. Benvenuti et al. [BCDMS Collaboration], Phys. Lett. B 223 (1989) 490–496.
  * [10] J. Blümlein, M. Klein, G. Ingelman and R. Rückl, Z. Phys. C 45 (1990) 501–513.
  * [11] M. Buza, Y. Matiounine, J. Smith, R. Migneron and W.L. van Neerven, Nucl. Phys. B 472 (1996) 611–658 [hep-ph/9601302].
  * [12] I. Bierenbaum, J. Blümlein and S. Klein, Nucl. Phys. B 780 (2007) 40–75 [hep-ph/0703285].
  * [13] M. Buza, Y. Matiounine, J. Smith and W.L. van Neerven, Nucl. Phys. B 485 (1997) 420–456 [hep-ph/9608342].
  * [14] A. Behring, I. Bierenbaum, J. Blümlein, A. De Freitas, S. Klein and F. Wißbrock, Eur. Phys. J. C 74 (2014) no.9, 3033 [arXiv:1403.6356 [hep-ph]].
  * [15] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round, C. Schneider, and F. Wißbrock, Nucl. Phys. B 886 (2014) 733–823 [arXiv:1406.4654 [hep-ph]].
  * [16] A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C. Schneider, Nucl. Phys. B 897 (2015) 612–644 [arXiv:1504.08217 [hep-ph]].
  * [17] J. Blümlein, A. De Freitas, M. Saragnese, C. Schneider and K. Schönwald, The Logarithmic Contributions to the Polarized $O(\alpha^{3}_{s})$ Asymptotic Massive Wilson Coefficients and Operator Matrix Elements in Deeply Inelastic Scattering, Phys. Rev. D (2021) in print [arXiv:2105.09572 [hep-ph]].
  * [18] S.I. Alekhin and J. Blümlein, Phys. Lett. B 594 (2004) 299–307 [hep-ph/0404034].
  * [19] S. Alekhin, J. Blümlein and S. Moch, Phys. Rev. D 86 (2012) 054009 [arXiv:1202.2281 [hep-ph]].
  * [20] J. Blümlein and H. Böttcher, Phys. Lett. B 662 (2008) 336–340 [arXiv:0802.0408 [hep-ph]].
  * [21] D. Boer et al., Gluons and the quark sea at high energies: Distributions, polarization, tomography, arXiv:1108.1713 [nucl-th].
  * [22] J.L. Abelleira Fernandez et al. [LHeC Study Group], J. Phys. G 39 (2012) 075001 [arXiv:1206.2913 [physics.acc-ph]].
  * [23] A. Vogt, Comput. Phys. Commun. 170 (2005) 65–92 [hep-ph/0408244].
  * [24] M. Boscolo, J.P. Delahaye and M. Palmer, The future prospects of muon colliders and neutrino factories, [arXiv:1808.01858 [physics.acc-ph]].
  * [25] J. Blümlein, Prog. Part. Nucl. Phys. 69 (2013) 28–84 [arXiv:1208.6087 [hep-ph]].
  * [26] N. Schmitz, Neutrinophysik, (Teunber, Stuttgart, 1997).
  * [27] J. D. Bjorken, Phys. Rev. D 1 (1970) 1376–1379.
  * [28] S. Moch, J.A.M. Vermaseren and A. Vogt, Nucl. Phys. B 688 (2004) 101–134 [hep-ph/0403192].
  * [29] A. Vogt, S. Moch and J.A.M. Vermaseren, Nucl. Phys. B 691 (2004) 129–181 [hep-ph/0404111].
  * [30] J. Blümlein, G. Falcioni and A. De Freitas, Nucl. Phys. B 910 (2016) 568–617 [arXiv:1605.05541 [hep-ph]].
  * [31] I. Bierenbaum, J. Blümlein and S. Klein, Nucl. Phys. B 820 (2009) 417–482 [arXiv:0904.3563 [hep-ph]].
  * [32] J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, C. Schneider and F. Wißbrock, Nucl. Phys. B 921 (2017) 585–688 [arXiv:1705.07030 [hep-ph]].
  * [33] K.G. Chetyrkin, B.A. Kniehl and M. Steinhauser, Phys. Rev. Lett. 79 (1997) 2184–2187 [hep-ph/9706430].
  * [34] W.A. Bardeen, A.J. Buras, D.W. Duke and T. Muta, Phys. Rev. D 18 (1978) 3998–4017.
  * [35] D.J. Gross and F. Wilczek, Phys. Rev. Lett. 30 (1973) 1343–1346;
H.D. Politzer, Phys. Rev. Lett. 30 (1973) 1346–1349;
W.E. Caswell, Phys. Rev. Lett. 33 (1974) 244–246;
D.R.T. Jones, Nucl. Phys. B 75 (1974) 531–538;
O.V. Tarasov, A.A. Vladimirov and A.Y. Zharkov, Phys. Lett. B 93 (1980)
429–432;
S.A. Larin and J.A.M. Vermaseren, Phys. Lett. B 303 (1993) 334–336 [hep-
ph/9302208];
T. van Ritbergen, J.A.M. Vermaseren and S.A. Larin, Phys. Lett. B 400 (1997)
379–384 [hep-ph/9701390];
M. Czakon, Nucl. Phys. B 710 (2005) 485–498 [hep-ph/0411261];
P.A. Baikov, K.G. Chetyrkin and J.H. Kühn, Phys. Rev. Lett. 118 (2017) no.8,
082002 [arXiv:1606.08659 [hep-ph]];
F. Herzog, B. Ruijl, T. Ueda, J.A.M. Vermaseren and A. Vogt, JHEP 02 (2017)
090 [arXiv:1701.01404 [hep-ph]];
T. Luthe, A. Maier, P. Marquard and Y. Schröder, JHEP 10 (2017) 166
[arXiv:1709.07718 [hep-ph]];
K.G. Chetyrkin, G. Falcioni, F. Herzog and J.A.M. Vermaseren, JHEP 10 (2017)
179 [arXiv:1709.08541 [hep-ph]].

  * [36] D.J. Broadhurst, N. Gray and K. Schilcher, Z. Phys. C 52 (1991) 111–122.
  * [37] N. Gray, D.J. Broadhurst, W. Grafe and K. Schilcher, Z. Phys. C 48 (1990) 673–680.
  * [38] I. Bierenbaum, J. Blümlein, S. Klein and C. Schneider, Nucl. Phys. B 803 (2008) 1–41 [arXiv:0803.0273 [hep-ph]].
  * [39] I. Bierenbaum, J. Blümlein, A. De Freitas, S. Klein, and K. Schönwald, The $O(\alpha_{s}^{2})$ Polarized Heavy Flavor Corrections to Deep-Inelastic Scattering at $Q^{2}\gg m^{2}$, DESY 15–004.
  * [40] J. Blümlein, Comput. Phys. Commun. 159 (2004) 19–54 [hep-ph/0311046].
  * [41] J.A.M. Vermaseren, Int. J. Mod. Phys. A 14 (1999) 2037–2076 [hep-ph/9806280].
  * [42] J. Blümlein and S. Kurth, Phys. Rev. D 60 (1999) 014018 [hep-ph/9810241].
  * [43] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C. Schneider, Nucl. Phys. B 890 (2014) 48–151 [arXiv:1409.1135 [hep-ph]].
  * [44] J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round and C. Schneider, Nucl. Phys. B 885 (2014) 280–317 [arXiv:1405.4259 [hep-ph]].
  * [45] J. Blümlein, J. Ablinger, A. Behring, A. De Freitas, A. von Manteuffel, C. Schneider and C. Schneider, PoS (QCDEV2017) 031 [arXiv:1711.07957 [hep-ph]].
  * [46] J. Ablinger et al., DESY 15–112.
  * [47] M. Diemoz, F. Ferroni, E. Longo and G. Martinelli, Z. Phys. C 39 (1988) 21–37.
  * [48] M. Glück, E. Reya and A. Vogt, Z. Phys. C 48 (1990) 471–482.
  * [49] J. Blümlein and A. Vogt, Phys. Rev. D 58 (1998) 014020 [hep-ph/9712546].
  * [50] R.K. Ellis, Z. Kunszt and E.M. Levin, Nucl. Phys. B 420 (1994) 517–549, Erratum: [Nucl. Phys. B 433 (1995) 498].
  * [51] J. Blümlein, Comput. Phys. Commun. 180 (2009) 2218–2249 [arXiv:0901.3106 [hep-ph]].
  * [52] J. Blümlein, Clay Math. Proc. 12 (2010) 167–188 [arXiv:0901.0837 [math-ph]].
  * [53] J.A.M. Vermaseren, A. Vogt and S. Moch, Nucl. Phys. B 724 (2005) 3–182 [hep-ph/0504242].
  * [54] S. Moch, J.A.M. Vermaseren and A. Vogt, Nucl. Phys. B 813 (2009) 220–258 [arXiv:0812.4168 [hep-ph]].
  * [55] D.J. Gross and F. Wilczek, Phys. Rev. D 8 (1973) 3633–3652.
  * [56] H. Georgi and H.D. Politzer, Phys. Rev. D 9 (1974) 416–420.
  * [57] E.G. Floratos, D.A. Ross and C.T. Sachrajda, Nucl. Phys. B 129 (1977) 66–88 [Erratum-ibid. B 139 (1978) 545–546].
  * [58] A. Gonzalez-Arroyo, C. Lopez and F.J. Yndurain, Nucl. Phys. B 153 (1979) 161–186.
  * [59] A. Gonzalez-Arroyo, C. Lopez and F.J. Yndurain, Nucl. Phys. B 159 (1979) 512–527.
  * [60] G. Curci, W. Furmanski and R. Petronzio, Nucl. Phys. B 175 (1980) 27–92.
  * [61] E.G. Floratos, R. Lacaze and C. Kounnas, Phys. Lett. 98B (1981) 89–95.
  * [62] E.G. Floratos, C. Kounnas and R. Lacaze, Nucl. Phys. B 192 (1981) 417–462.
  * [63] S. Moch and J.A.M. Vermaseren, Nucl. Phys. B 573 (2000) 853–907 [hep-ph/9912355].
  * [64] A. Behring, J. Blümlein, A. De Freitas, A. Goedicke, S. Klein, A. von Manteuffel, C. Schneider and K. Schönwald, Nucl. Phys. B 948 (2019) 114753 [arXiv:1908.03779 [hep-ph]].
  * [65] W. Furmanski and R. Petronzio, Z. Phys. C 11 (1982) 293–314.
  * [66] W.L. van Neerven and E.B. Zijlstra, Phys. Lett. B 272 (1991) 127–133.
  * [67] E. B. Zijlstra and W. L. van Neerven, Phys. Lett. B 297 (1992) 377–384.
  * [68] E. B. Zijlstra and W. L. van Neerven, Nucl. Phys. B 383 (1992) 525–574.
  * [69] K. Sasaki, Prog. Theor. Phys. 54 (1975) 1816–1827.
  * [70] M.A. Ahmed and G.G. Ross, Phys. Lett. 56B (1975) 385–390.
  * [71] G. Altarelli and G. Parisi, Nucl. Phys. B 126 (1977) 298–318.
  * [72] R. Mertig and W. L. van Neerven, Z. Phys. C 70 (1996) 637–654 [hep-ph/9506451].
  * [73] W. Vogelsang, Phys. Rev. D 54 (1996) 2023–2029 [hep-ph/9512218].
  * [74] W. Vogelsang, Nucl. Phys. B 475 (1996) 47–72 [hep-ph/9603366].
  * [75] E.B. Zijlstra and W.L. van Neerven, Nucl. Phys. B 417 (1994) 61–100 Errata: [Nucl. Phys. B 426 (1994) 245]; [Nucl. Phys. B 773 (2007) 105-106]; [Nucl. Phys. B 501 (1997) 599].
  * [76] J. Kodaira, S. Matsuda, T. Muta, K. Sasaki and T. Uematsu, Phys. Rev. D 20 (1979) 627–629.
  * [77] J. Kodaira, S. Matsuda, K. Sasaki and T. Uematsu, Nucl. Phys. B 159 (1979) 99–124.
  * [78] I. Antoniadis and C. Kounnas, Phys. Rev. D 24 (1981) 505–525.
  * [79] P.A. Baikov and K.G. Chetyrkin, Nucl. Phys. Proc. Suppl. 160 (2006) 76–79.
  * [80] P.A. Baikov, K.G. Chetyrkin and J.H. Kühn, Nucl. Part. Phys. Proc. 261-262 (2015) 3–18 [arXiv: 1501.06739 [hep-ph]].
  * [81] V.N. Velizhanin, Four loop anomalous dimension of the third and fourth moments of the non-singlet twist-2 operator in QCD, arXiv:1411.1331 [hep-ph].
  * [82] B. Ruijl, T. Ueda, J.A.M. Vermaseren, J. Davies and A. Vogt, PoS (LL2016) 071 [arXiv:1605.08408 [hep-ph]].
  * [83] J. Davies, A. Vogt, B. Ruijl, T. Ueda and J.A.M. Vermaseren, Nucl. Phys. B 915 (2017) 335–362 [arXiv:1610.07477 [hep-ph]].
  * [84] S. Moch, B. Ruijl, T. Ueda, J.A.M. Vermaseren and A. Vogt, JHEP 1710 (2017) 041 [arXiv:1707.08315 [hep-ph]].
  * [85] J.A. Gracey, Phys. Lett. B 322 (1994) 141–146 [hep-ph/9401214].
  * [86] R. Kirschner and L.N. Lipatov, Nucl. Phys. B 213 (1983) 122–148.
  * [87] J. Blümlein and A. Vogt, Phys. Lett. B 370 (1996) 149–155 [hep-ph/9510410].
  * [88] J. Bartels, B.I. Ermolaev and M.G. Ryskin, Z. Phys. C 70 (1996) 273–280 [hep-ph/9507271].
  * [89] J. Blümlein and H. Böttcher, Nucl. Phys. B 841 (2010) 205–230 [arXiv:1005.3113 [hep-ph]].
  * [90] S. Alekhin, J. Blümlein, K. Daum, K. Lipka and S. Moch, Phys. Lett. B 720 (2013) 172–176 [arXiv: 1212.2355 [hep-ph]].
  * [91] P.A. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020) and 2021 update.
  * [92] W.A. Bardeen and A.J. Buras, Phys. Lett. B 86 (1979) 61–66, Erratum: [Phys. Lett. B 90 (1980) 485].
  * [93] A.J. Buras, Rev. Mod. Phys. 52 (1980) 199–276.
  * [94] M. Glück and E. Reya, Phys. Rev. D 25 (1982) 1211–1217.
  * [95] G. Grunberg, Phys. Rev. D 29 (1984) 2315–2338.
  * [96] S. Catani, Z. Phys. C 75 (1997) 665–678 [hep-ph/9609263].
  * [97] R.S. Thorne, Nucl. Phys. B 512 (1998) 323–392 [hep-ph/9710541].
  * [98] J. Blümlein, V. Ravindran and W.L. van Neerven, Nucl. Phys. B 586 (2000) 349–381 [hep-ph/0004172].
  * [99] F.J. Yndurain The Theory of Quark and Gluon Interactions, (Springer, Berlin, 1999), 3rd Edition.
  * [100] J. Blümlein and A. Guffanti, Nucl. Phys. Proc. Suppl. 152 (2006) 87–91 [hep-ph/0411110].
  * [101] C. Adloff et al. [H1], Phys. Lett. B 393 (1997), 452-464 [arXiv:hep-ex/9611017 [hep-ex]].
  * [102] V. Andreev et al. [H1], Eur. Phys. J. C 74 (2014) no.4, 2814 [arXiv:1312.4821 [hep-ex]].
  * [103] S. Chekanov et al. [ZEUS], Phys. Lett. B 682 (2009), 8-22 [arXiv:0904.1092 [hep-ex]].
  * [104] P. Agostini et al. [LHeC and FCC-he Study Group], The Large Hadron-Electron Collider at the HL-LHC, arXiv:2007.14491 [hep-ex].


# A Complex-LASSO Approach for Localizing Forced Oscillations in Power Systems
††thanks: This material is based upon work supported by the National Science
Foundation under Grant No. OAC-1934766 and PSERC project S-87.

Rajasekhar Anguluri, Nima Taghipourbazargani, Oliver Kosut, and Lalitha Sankar
School of Electrical, Computer, and Energy Engineering
Arizona State University, Tempe, AZ 85281 USA
<EMAIL_ADDRESS>
###### Abstract

We study the problem of localizing multiple sources of forced oscillations
(FOs) and estimating their characteristics, such as frequency, phase, and
amplitude, using noisy PMU measurements. For each source location, we model
the input oscillation as a sum of unknown sinusoidal terms. This allows us to
obtain a linear relationship between measurements and the inputs at the
unknown sinusoids’ frequencies in the frequency domain. We determine these
frequencies by thresholding the empirical spectrum of the noisy measurements.
Assuming sparsity in the number of FOs’ locations and the number of sinusoids
at each location, we cast the location recovery problem as an
$\ell_{1}$-regularized least squares problem in the complex domain—i.e.,
complex-LASSO (linear shrinkage and selection operator). We numerically solve
this optimization problem using the complex-valued coordinate descent method,
and show its efficiency on the IEEE 68-bus, 16 machine and WECC 179-bus,
29-machine systems.

###### Index Terms:

Forced oscillations, complex-LASSO, sparsity, sampled data system, PMU
measurements.

##  I Introduction

Detecting and localizing forced oscillations (FOs) is crucial for safety and
reliability of power systems. Detection helps determine if the manifested
oscillations in PMU measurements are FOs; whereas localization helps to
mitigate the oscillations either by disabling the sources (e.g., malfunctioned
controllers and cyclic loads) or by injecting certain counteracting signals[1,
2]. Compared to localization, detecting oscillations is easier and has been
well studied in power systems [3, 4].

We focus on localizing the sources of FO, which amounts to finding $m^{*}$
true sources from $m$ possible sources. A brute-force search requires
searching across ${m\choose m^{*}}$ configurations, which is computationally
intractable for large $m$. One way to tackle this challenge is to find source
locations by minimizing performance measures, such as system-theoretic norms
(e.g., $\mathcal{H}_{2}$ and $\mathcal{H}_{\infty}$) and information-theoretic
based measures [5].

Our approach aims to leverage the fact that the FO sources are sparse relative
to the number of possible sources ($m^{*}\ll m$) [6, 1]. In fact, we find the
source locations and the associated input parameters using an $\ell_{1}$-norm
regularized optimization problem. In several inverse problems with sparsity
constraints—finding unknown sparse parameters using a few
measurements—$\ell_{1}$-norm regularization has shown to accurately recover
the sparsity pattern of the unknown parameter than the standard
$\ell_{2}$-norm regularization [7]. Sparsity has been long recognized in power
system applications, especially in state and topology estimation problems;
however, sparsity methods are less explored for localizing FO sources [8].

By working in the frequency-domain, we encode sparsity of the number of
locations and the number of sinusoids jointly in a single unknown vector
$\mathbf{U}_{K}$, where $\mathbf{U}_{K}$ satisfies
$\mathbf{Y}_{K}=\mathbf{G}_{K}\mathbf{U}_{K}$ (see Section II for more
details). This kind of joint sparsity might not be possible in the time-
domain, especially, for arbitrary forced inputs [9].

We summarize our FO location recovery method below:

  1. 1.

Using the discrete time Fourier transform (DTFT), we obtain the complex-valued
linear model $\mathbf{Y}_{K}=\mathbf{G}_{K}\mathbf{U}_{K}$ that relates
measurements and the unknown oscillatory inputs at frequencies at which these
inputs oscillate. Here, $K$ is the number of frequencies at all locations.

  2. 2.

We then use the fast Fourier transform (FFT) to determine these frequencies by
thresholding the spectrum of measurements collected over a finite time
interval.

  3. 3.

For the model in step 1, we solve an $\ell_{1}$ regularized least squares
(henceforth, LASSO) problem in the complex domain via coordinated descent
method [10] to infer: (i) the number of true locations and the number of
sinusoids at any given location, and (ii) the frequency, phase, and amplitude
of each sinusoid.

Our recovery algorithm is simple as it requires thresholding the spectrum and
solving a convex optimization problem. We validate its performance on two
benchmark systems.

Related literature: In [11, 12], sources were identified using spectral
properties of the transfer functions between sources and measurements. Along
these lines, in [3], sources are localized using ratios between signal content
at the harmonics and the fundamental frequency. By representing generators’
frequency response as an effective admittance matrix and oscillations as
current injections, [13] identifies sources by comparing predicted against
measured current spectrum, and [14] uses a Bayesian approach. Finally, in [8],
a sparse plus low rank decomposition of Hankel measurement matrices has been
suggested to localize the sources.

However, unlike our method, algorithms in the above studies localize sources
(mostly single) but do not jointly estimate the input. Further, with the
exception of works in [8] and [14], sparsity has not been explored in a
systematic way. The work closest to ours is [14], which considers an
$\ell_{1}$-norm regularized optimization framework in a Bayesian setting.
However, [14] models FOs as current injections and generators as frequency
response functions. Instead, we model FOs as exogenous inputs, which is
applicable to any dynamical system, not just power systems. Moreover, we
consider practical aspects, such as spectral leakage in the frequency domain.

##  II System dynamics under forced oscillations

Consider the following linearized multi-machine dynamics excited by external
unknown inputs in state-space form:

$\displaystyle\begin{split}\dot{\mathbf{x}}(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t),\quad\mathbf{y}(t)=\mathbf{C}\mathbf{x}(t),\quad\forall
t\geq 0.\end{split}$ (1)

The state vector $\mathbf{x}(t)\in\mathbb{R}^{n}$ includes machine internal
dynamics and variables associated with closed-loop controllers. We model the
input $\mathbf{u}(t)\in\mathbb{R}^{m}$ as the forced oscillation vector:

$\displaystyle\mathbf{u}(t)\triangleq\begin{bmatrix}u_{1}(t)\\\ u_{2}(t)\\\
\vdots\\\
u_{m}(t)\end{bmatrix}=\begin{bmatrix}\sum_{l=1}^{M_{1}}a_{1,l}\sin(\omega_{1,l}t+\phi_{1,l})\\\
\sum_{l=1}^{M_{2}}a_{2,l}\sin(\omega_{2,l}t+\phi_{2,l})\\\ \vdots\\\
\sum_{l=1}^{M_{m}}a_{m,l}\sin(\omega_{m,l}t+\phi_{m,l})\end{bmatrix}$ (2)

where $a_{r,l}\geq 0$, $\omega_{r,l}=2\pi f_{r,l}\geq 0$, and $\phi_{r,l}$ are
the amplitude, angular frequency, and phase of the $l^{th}$ sinusoid term at
the $r^{th}$ location, and $r\in\\{1,\ldots,m\\}$. Finally,
$\mathbf{y}(t)\in\mathbb{R}^{p}$ is the measurement obtained from $p$ sensors.
For simplicity, we build the theory using noise-free measurements, and later
consider measurement noise in simulations.

By modeling forced input as a sum of multiple weighted sinusoids, we allow for
different oscillatory waveforms. Thus, our approach generalizes the approaches
recovering a single sinusoidal input, which are well studied in the literature
[6, 1]. Let
$\mathbf{a}_{r}(t)=[a_{r,1}\sin(\omega_{r,1}t+\phi_{r,1})\,\cdots\,a_{r,M_{r}}\sin(\omega_{r,M_{r}}t+\phi_{r,M_{r}})]$
be the vector of sinusoids at the $r^{th}$ location. We make the following
assumption on $\mathbf{u}(t)$ in (1) and $\mathbf{a}_{r}(t)$.

###### Assumption 1

Both $\mathbf{u}(t)$ in (1) and $\mathbf{a}_{r}(t)$ defined above are sparse;
that is, $\|\mathbf{u}(t)\|_{0}\ll m$ and $\|\mathbf{a}_{r}(t)\|_{0}\ll
M_{r}$, where $\|\mathbf{z}\|_{0}$ counts the number of non-zero entries in
the vector $\mathbf{z}$.

As PMUs record measurements at discrete time instants, we sample (1) with the
sampling period $T$ (e.g., $F\\!=\\!1/T\\!=\\!30$–$60$ Hz for PMUs). Let
$\mathbf{A}_{d}=\exp(\mathbf{A}T)$ and
$\mathbf{B}_{d}=\int_{0}^{T}\exp(\mathbf{A}(T-s))ds\,\mathbf{B}$. Let
$k=0,1,\ldots$ and define $\mathbf{x}[k]\triangleq\mathbf{x}(kT)$,
$\mathbf{u}[k]\triangleq\mathbf{u}(kT)$, and
$\mathbf{y}[k]\triangleq\mathbf{y}(kT)$. Suppose that $\mathbf{u}(t)$ is a
piecewise constant111Frequencies of real FOs are smaller than the sampling
frequency of PMUs. during $kT\leq t\leq(k+1)T$. Then

$\displaystyle\begin{split}\mathbf{x}[k+1]\\!=\\!\mathbf{A}_{d}\mathbf{x}[k]+\mathbf{B}_{d}\mathbf{u}[k],\quad\mathbf{y}[k]\\!=\\!\mathbf{C}\mathbf{x}[k]\end{split}$
(3)

describe the discrete-time sampled model of (1). As our focus is on
oscillations triggered by inputs, we assume $\mathbf{x}[0]=\mathbf{0}$.

Define the matrix-valued transfer function
$\mathbf{H}[z]\\!=\\!\mathbf{C}(z\mathbf{I}-\mathbf{A}_{d})^{-1}\mathbf{B}_{d}\in\mathbb{C}^{p\times
m}$, where $z\in\mathbb{C}$. Then, from the standard linear system analysis,
we conclude that $\mathbf{Y}[z]\\!=\\!\mathbf{H}[z]\mathbf{U}[z]$, where
$\mathbf{Y}[z]\in\mathbb{C}^{p\times 1}$ and
$\mathbf{U}[z]\in\mathbb{C}^{m\times 1}$ are the $Z$-transforms of
$\mathbf{y}[k]$ and $\mathbf{u}[k]$. Finally, at $z=\exp(j\Omega)$, where
$\Omega\in(0,2\pi)$ and $j^{2}=-1$, we have the DTFT representation of (3):

$\displaystyle\mathbf{Y}[e^{j\Omega}]=\mathbf{H}[e^{j\Omega}]\mathbf{U}[e^{j\Omega}].$
(4)

The benefit of working in the Fourier (or frequency) domain is that
$\mathbf{U}[e^{j\Omega}]$ can be expressed as a sum of weighted Dirac delta
functions. As a result, $\mathbf{Y}[e^{j\Omega}]$ can be expanded in terms of
basis functions that encode both the source location and frequency of the
input sinuosoids. Using this observation, we later show that the source
localization problem can be cast as a solution to a simple $\ell_{1}$
regularized least-squares problem.

Since $\mathbf{u}(t)$ in (2) is a sum of sinusoids, it follows that

$\displaystyle\mathbf{U}[e^{j\Omega}]\\!=\\!j\pi\begin{bmatrix}\sum_{l=1}^{M_{1}}a_{1,l}\,e^{j\phi_{1,l}}[\delta(\Omega\\!+\\!\widetilde{\omega}_{1,l})\\!-\\!\delta(\Omega\\!-\\!\widetilde{\omega}_{1,l})]\\\
\sum_{l=1}^{M_{2}}a_{2,l}\,e^{j\phi_{2,l}}[\delta(\Omega\\!+\\!\widetilde{\omega}_{2,l})\\!-\\!\delta(\Omega\\!-\\!\widetilde{\omega}_{2,l})]\\\
\vdots\\\
\sum_{l=1}^{M_{m}}a_{m,l}\,e^{j\phi_{m,l}}[\delta(\Omega\\!+\\!\widetilde{\omega}_{m,l})\\!-\\!\delta(\Omega\\!-\\!\widetilde{\omega}_{m,l})]\end{bmatrix}$
(5)

where $\delta(x)$ is the Dirac delta function and
$\widetilde{\omega}_{k,l}=\omega_{k,l}T$. Note that $\mathbf{U}[e^{j\Omega}]$
is non-zero only for
$\Omega\in\\{\pm\widetilde{\omega}_{1,1},\ldots,\pm\widetilde{\omega}_{m,M_{m}}\\}$.
Suppose that for $\Omega$ in this set,
$\mathbf{H}[e^{j\Omega}]\neq\mathbf{0}$. Then, from (4), we have
$\mathbf{Y}[e^{j\Omega}]\neq\mathbf{0}$. If this is not the case, we cannot
recover the sinusoid oscillating with $\Omega$ using measurements. Finally, we
have $\mathbf{Y}[e^{j\Omega}]=\mathbf{Y}^{*}[e^{-j\Omega}]$. Thus, we focus
only on $\widetilde{\omega}_{k,l}\geq 0$.

With a slight abuse of notation, let
$\\{\widetilde{\omega}_{1},\ldots,\widetilde{\omega}_{K}\\}$ be the set of
frequencies where
$\mathbf{Y}[e^{j\widetilde{\omega}_{l}}]=\mathbf{H}[e^{j\widetilde{\omega}_{l}}]\mathbf{U}[e^{j\widetilde{\omega}_{l}}]\neq\mathbf{0}$,
where both $\widetilde{\omega}_{l}$ and $K$ are unknown. Consider the
following model:

$\displaystyle\underbrace{\begin{bmatrix}\mathbf{Y}[e^{j\widetilde{\omega}_{1}}]\\\
\mathbf{Y}[e^{j\widetilde{\omega}_{2}}]\\\ \vdots\\\
\mathbf{Y}[e^{j\widetilde{\omega}_{K}}]\end{bmatrix}}_{\triangleq\mathbf{Y}_{K}}\\!=\\!\underbrace{\begin{bmatrix}\mathbf{H}[e^{j\widetilde{\omega}_{1}}]&\\\
&\mathbf{H}[e^{j\widetilde{\omega}_{2}}]\\\ &&\ddots\\\
&&&\mathbf{H}[e^{j\widetilde{\omega}_{K}}]\end{bmatrix}}_{\triangleq\mathbf{H}_{K}}\underbrace{\begin{bmatrix}\mathbf{U}[e^{j\widetilde{\omega}_{1}}]\\\
\mathbf{U}[e^{j\widetilde{\omega}_{2}}]\\\ \vdots\\\
\mathbf{U}[e^{j\widetilde{\omega}_{K}}]\end{bmatrix}}_{\triangleq\mathbf{U}_{K}}$
(6)

where $\mathbf{Y}_{K}\in\mathbb{C}^{pK\times 1}$,
$\mathbf{H}_{K}\in\mathbb{C}^{pK\times mK}$, and
$\mathbf{U}_{K}\in\mathbb{C}^{mK\times 1}$. In light of Assumption 1, note
that $K\ll[M_{1}+\ldots+M_{m}]$. Further, the non-zero components of
$\mathbf{U}[e^{j\widetilde{\omega}_{l}}]\in\mathbb{C}^{m}$ correspond to
locations with the inputs containing sinusoids of frequency
$\widetilde{\omega}_{l}$. Thus, from the sparsity pattern of $\mathbf{U}_{K}$
alone we can determine source locations and their sinusoidal frequencies (see
Fig. 1). The values of $\mathbf{U}_{K}$ provide phase and amplitude
information.

The model in (6) captures the relationship between measurements and the input
locations at multiple frequencies; thus, for a sinusoid input, i.e.,
$\mathbf{U}[e^{j\widetilde{\omega}_{l}}]$ is a delta function,
$\mathbf{Y}[e^{j\widetilde{\omega}_{l}}]$ is infinite for any $l$. This is a
natural aspect of taking DTFTs of infinite time-length sinusoids. However, in
practice, we only use measurements for a finite time interval and use FFTs to
compute $\mathbf{Y}[e^{j\widetilde{\omega}_{l}}]$, which is consequently
finite.

Figure 1: (a) shows the discrete-time power system model with $m=4$ possible
locations but only locations $1$ and $4$ are injecting oscillations. Each of
these inputs is a sum of two sinusoids with different parameters except for
${\widetilde{\omega}}_{11}={\widetilde{\omega}}_{42}$. In Fig. (b), input
frequencies are found by computing the spectrum of the two sensory
measurements. Since ${\widetilde{\omega}}_{11}={\widetilde{\omega}}_{42}$, we
see only three dominant peaks in the spectrum, i.e., $K=3$, and the smaller
peaks are due to spectral leakage. In Fig. (c), we visualize linear model in
(6). For any $\widetilde{\omega}_{k,l}$, the four entries (square blocks) in
$\mathbf{U}[e^{j\widetilde{\omega}_{k,l}}]$ correspond to four possible
locations. Locations having sinusoid oscillating with
$\widetilde{\omega}_{k,l}$ are highlighted in color. The top input
$\mathbf{U}[e^{j\widetilde{\omega}_{1,1}}]$ has two non-zero entries
(highlighted in brown and green) indicating that locations $1$ and $4$ have a
sinusoid with $\widetilde{\omega}_{1,1}=\widetilde{\omega}_{4,2}$.

##  III Complex-LASSO for Source Localization

In practice, the total number of sensors (p) could be less than the possible
number of sources (m). For example, PMUs can measure only bus level quantities
but not the internal signals of control devices. Thus, model in (6) is an
under-determined system and we cannot obtain $\mathbf{U}_{K}$ by means of
ordinary least squares. As a result, for the model in (6), we consider the
$\ell_{1}$-regularized optimization problem to both localize the sources and
identify the frequencies of the sinusoids at each source:

$\displaystyle\widehat{\mathbf{U}}_{K}\\!=\\!\operatorname*{arg\,min}_{\begin{subarray}{c}\mathbf{U}_{K}\in\mathbb{C}^{mK\times
1}\end{subarray}}\left\\{\frac{1}{2}\left\lVert\mathbf{Y}_{K}-\mathbf{H}_{K}\mathbf{U}_{K}\right\rVert_{2}^{2}\\!+\\!\lambda\|\mathbf{U}_{K}\|_{1}\right\\}$
(7)

where $\lambda\geq 0$ is the tuning parameter. For the vector $\mathbf{z}$,
the $\ell_{1}$-norm is $\|\mathbf{z}\|_{1}=\sum_{i}|z_{i}|$, where
$|z_{i}|=\sqrt{\text{Re}(z_{i})^{2}+\text{Im}(z_{i})^{2}}$. The regularization
term $\|\mathbf{U}_{K}\|_{1}$ promotes sparsity in $\widehat{\mathbf{U}}_{K}$.
We henceforth refer to the problem in (7) as the complex-LASSO.

Akin to $\mathbf{U}_{K}$ in (6), define
$\widehat{\mathbf{U}}_{K}^{\mathsf{T}}=[\widehat{\mathbf{U}}[e^{j\widetilde{\omega}_{1}}]^{\mathsf{T}}\ldots\widehat{\mathbf{U}}[e^{j\widetilde{\omega}_{K}}]^{\mathsf{T}}]$.
The block diagonal form of $\mathbf{H}_{K}$ and additive property of the
$\ell_{1}$-norm allow us to compute
$\widehat{\mathbf{U}}[e^{j\widetilde{\omega}_{l}}]$ in a distributed fashion:

$\displaystyle\widehat{\mathbf{U}}[e^{j\widetilde{\omega}_{l}}]\\!=\\!\operatorname*{arg\,min}_{\begin{subarray}{c}\mathbf{U}\in\mathbb{C}^{m\times
1}\end{subarray}}\left\\{\frac{1}{2}\left\lVert\mathbf{Y}[e^{j\widetilde{\omega}_{l}}]-\mathbf{H}[e^{j\widetilde{\omega}_{l}}]\mathbf{U}\right\rVert_{2}^{2}\\!+\\!\lambda\|\mathbf{U}\|_{1}\right\\}.$

The above optimization problem is computationally convenient and extremely
useful to quickly determine the sources injecting oscillations with a
particular frequency of interest.

To solve (7), we need $\mathbf{H}_{K}$, $\mathbf{Y}_{K}$, and $K$. The matrix
$\mathbf{H}_{K}$ is computed using the power system matrices
$(\mathbf{A}_{d},\mathbf{B}_{d},\mathbf{C})$. If not available, one can use
empirically determined transfer functions. To compute $K$ and
$\mathbf{Y}_{K}$, we obtain the vector-valued $N$-point DFT222In simulations,
for computational speedup, we use the FFT. of
$\mathbf{y}[L],\ldots,\mathbf{y}[N-1+L]$ using

$\displaystyle\widetilde{\mathbf{Y}}[q]\triangleq\frac{2}{N}\sum_{k=L}^{N-1+L}\mathbf{y}[k]e^{-j\frac{2\pi
lk}{N}}\quad(q=0,1,\ldots,N-1)$ (8)

where we choose $L\gg 0$ such that $\mathbf{y}[k]$, for all $k\geq L$, is in
steady state. For any $q$, recall that
$\|\widetilde{\mathbf{Y}}[q]\|_{\infty}=\max_{i}|\widetilde{\mathbf{Y}}^{(i)}[q]|$,
where $\widetilde{\mathbf{Y}}^{(i)}[q]$ is the $i^{th}$ scalar in
$\widetilde{\mathbf{Y}}^{(i)}[q]\in\mathbb{C}^{p}$.

The input (angular) frequencies can be determined using
$\widetilde{\omega}_{l}=2\pi l/N$, where $l$ satisfies
$\|\widetilde{\mathbf{Y}}[l]\|_{\infty}\neq 0$. However, due to the sensor
noise and spectral leakage, $\|\widetilde{\mathbf{Y}}[l]\|_{\infty}$ could be
non-zero even when $\widetilde{\omega}_{l}$ is not the true input frequency.
The noise can be attenuated by filtering the measurements. Instead, we reduce
spectral leakage by multiplying measurements with the Hamming window [15]. For
these processed measurements, let
$S=\\{l:\|\widetilde{\mathbf{Y}}[l]\|_{\infty}>\tau\\}$, where $\tau>0$ is the
user-defined threshold. Then $K=|S|$, where $|\cdot|$ is the cardinality of
the set $S$, gives us the total number of input frequencies. Finally, we
replace $\mathbf{Y}[e^{j\widetilde{\omega}_{l}}]$ in $\mathbf{Y}_{K}$ (given
by (6)) with $\widetilde{\mathbf{Y}}[l]$, and evaluate the corresponding
$\mathbf{H}[e^{j\widetilde{\omega}_{l}}]$ at $\widetilde{\omega}_{l}=2\pi
l/N$, where $l\in S$.

With $\mathbf{H}_{K}$ and $\mathbf{Y}_{K}$ at our disposal, we now can solve
(7) using coordinate descent method [10] in the complex domain. We summarize
our source recovery method in Algorithm 1, which recovers locations, and input
parameters in discrete-time domain. We obtain the continuous-time frequency as
$f_{l}=\widetilde{\omega}_{l}/(2\pi T)$. Instead, we recover amplitudes by
$|\widehat{\mathbf{U}}_{K}|/\pi$, where $|\cdot|$ is the complex magnitude,
applied for entry in $\widehat{\mathbf{U}}_{K}$.

Input: $\tau>0$, $\lambda\geq 0$, $\\{\mathbf{Y}[k]\\}_{k=L}^{N-1-L}$, and
$(\mathbf{A}_{d},\mathbf{B}_{d},\mathbf{C})$.

Step 1: Compute $\widetilde{\mathbf{Y}}[q]$ using (8) and process them using
the Hamming window.

Step 2: Set $S=\\{l:\|\widetilde{\mathbf{Y}}[l]\|_{\infty}>\tau\\}$ and
$K=|S|$.

Step 3: Let
$\mathbf{Y}_{K}=[\widetilde{\mathbf{Y}}[l_{1}]\ldots\widetilde{\mathbf{Y}}[l_{K}]]$.
Evaluate $\mathbf{H}[e^{j\widetilde{\omega}_{l_{i}}}]$ in (6) using
$\widetilde{\omega}_{l_{i}}=2\pi l_{i}/N$, and $l_{i}\in S$

Step 4: Solve complex-LASSO problem in (7) using the coordinate descent method
[10].

Return: $\widetilde{\omega}_{l}$;
$\widehat{\mathbf{U}}[e^{j\widetilde{\omega}_{l}}]$; and source locations:
indices of non-zero entries in
$\widehat{\mathbf{U}}[e^{j\widetilde{\omega}_{l}}]$, $l=\\{1,\ldots,K\\}$.

Algorithm 1 Source localization via complex-LASSO

##  IV Simulation Results

We apply Algorithm 1 to recover FO locations and oscillatory input parameters
in two benchmark power systems. We add Gaussian noise with the SNR of 10 dB to
the measurable quantities (see below). We let $1/T=F=30$ and $N=600$ .

Tuning parameter ($\lambda$) selection: For $\lambda=0$, the solution to (7)
is given by the least squares solution. For
$\lambda\geq\lambda_{\max}\triangleq\|\mathbf{Y}_{K}^{T}\mathbf{H}_{K}\|_{\infty}$,
$\widehat{\mathbf{U}}_{K}=\mathbf{0}$ (a fully sparse vector) [10, Section
2.5]. Thus, we set $\lambda=\alpha\lambda_{\text{max}}$ and choose
$\alpha\in[0,1]$ by performing a sensitivity analysis with respect to the true
positive rate (TPR)—ratio of the number of correctly identified non-zero
entries in $\widehat{\mathbf{U}}_{K}$ to the number of non-zero entries in
$\mathbf{U}_{K}$—and the false positive rate (FPR)—ratio of the number of
incorrectly identified non-zero entries in $\widehat{\mathbf{U}}_{K}$ to the
number of non-zero entries in $\mathbf{U}_{K}$.  TPR and FPR depend on the
noisy measurements and so does $\alpha$. Thus to get a reliable and
trustworthy $\alpha$, we perform sensitivity analysis for 20 realizations of
$\\{y[k]\\}_{k=L}^{N-1+L}$. We then pick an $\alpha$ that yields
$\text{TPR}=1$ and $\text{FPR}=0$ for as many realizations as possible.

TABLE I: True and estimated values of the amplitude, frequency and phase of
the forced oscillations inputs

System under study |  FO Location | Input Signal
---|---|---
Amplitude (p.u.) | Frequency (Hz) | Phase (rad)
Parameter | True | Estimate ($\sigma$)* | Parameter | True | Estimate | Parameter | True | Estimate ($\sigma$)*
68-bus system (near PMUs) | 4 | $a_{4,1}$ | $0.01$ | $0.0126~{}(0.003)$ | $f_{4,1}$ | $2$ | $2$ | $\phi_{4,1}$ | $0.3$ | $0.3302~{}(0.024)$
$a_{4,2}$ | $0.01$ | $0.0081~{}(0.0018)$ | $f_{4,2}$ | $2.5$ | $2.5$ | $\phi_{4,2}$ | $0.4$ | $0.4401~{}(0.012)$
10 | $a_{10,1}$ | $0.01$ | $0.0087~{}(0.005)$ | $f_{10,1}$ | $1.5$ | $1.5$ | $\phi_{10,1}$ | $0.1$ | $0.1326~{}(0.01)$
$a_{10,2}$ | $0.01$ | $0.0093~{}(0.006)$ | $f_{10,2}$ | $1$ | $1$ | $\phi_{10,2}$ | $0.3$ | $0.2804~{}(0.024)$
13 | $a_{13,1}$ | $0.02$ | $0.0169~{}(0.004)$ | $f_{13,1}$ | $3.5$ | $3.5$ | $\phi_{13,1}$ | $0.1$ | $0.0893~{}(0.007)$
$a_{13,2}$ | $0.02$ | $0.0170~{}(0.002)$ | $f_{13,2}$ | $0.8$ | $0.8$ | $\phi_{13,2}$ | $0.2$ | $0.1771~{}(0.0052)$
WECC-179 bus system | 5 | $a_{5,1}$ | $0.02$ | $0.0244~{}(0.002)$ | $f_{5,1}$ | $1$ | $1$ | $\phi_{5,1}$ | $0.3$ | $0.2993~{}(0.0133)$
$a_{5,2}$ | $0.01$ | $0.0117~{}(0.005)$ | $f_{5,2}$ | $0.8$ | $0.8$ | $\phi_{5,2}$ | $0.4$ | $0.3405~{}(0.0305)$
14 | $a_{14,1}$ | $0.03$ | $0.0214~{}(0.008)$ | $f_{14,1}$ | $0.7$ | $0.7$ | $\phi_{14,1}$ | $0.2$ | $0.2134~{}(0.0068)$
$a_{14,2}$ | $0.02$ | $0.0170~{}(0.001)$ | $f_{14,2}$ | $1.5$ | $1.5$ | $\phi_{14,2}$ | $0.3$ | $0.2613~{}(0.0048)$
27 | $a_{27,1}$ | $0.04$ | $0.0350~{}(0.007)$ | $f_{27,1}$ | $2$ | $2$ | $\phi_{27,1}$ | $0.1$ | $0.1076~{}(0.0098)$
$a_{27,2}$ | $0.01$ | $0.0094~{}(0.003)$ | $f_{27,2}$ | $1.2$ | $1.2$ | $\phi_{27,2}$ | $0.2$ | $0.2567~{}(0.0457)$

  * *

Standard deviation of estimated parameters over 20 realizations of noisy
measurements.

Figure 2: One-sided power spectrum of measurements for different case studies.
In each plot, we overlay the spectrum associated with different measurement
buses. The black dashed line represents the thresholding parameter $\tau$,
which equals $0.2\times 10^{-3}$ in Fig. (a) and (b), and $0.4\times 10^{-3}$
in Fig. (c). Theoretically, if inputs are observable, the spectrum of each PMU
measurement should contain dominant peaks at all input frequencies. However,
as evident in all plots, not even for one PMU, we can see dominant peaks at
all input frequencies. This can be attributed to the measurement noise and
spectral leakage. However, by considering all PMUs, we see sharp dominant
peaks at all input frequencies. Interestingly, for IEEE-68 bus case, the
magnitude of spectral peaks is smaller when PMUs are away from the true source
locations. Thus it makes sense to consider measurements from nearby PMUs if we
have prior coarse knowledge of source locations. Figure 3: Sensitivity
analysis of the scaling parameter $\alpha$ in terms of TPR and FPR for
different cases. The error bars in the figures illustrate the maximum and
minimum of twenty values (as described in the main text) of TPR and FPR.
Instead, the solid line is the average value Interestingly, from panels (a)
and (b), we note that the PMUs in the vicinity of the FO sources resulted in a
wider range of $\alpha$ values that yield a better TPR ($100\%$) and FPR
($0\%$).

###  IV-A Case 1: IEEE 68 bus, 16 machine system [16]

Figure 4: IEEE NETS/NYPS 16 machine 68 bus system [16]. Among $m=16$
excitation control inputs (green circles), only $m^{*}=3$ locations (red
colored bus no.s: $4$, $10$, $13$) are excited by FO inputs.

Each generator is represented by ten states: the rotor angle, angular
frequency, damper winding flux leakages, and states corresponding to
excitation systems. We obtain the state space matrices using the power system
toolbox (PST) [17]. See Fig. 4 and Table I for the inputs’ description. We
consider scenarios where PMUs are near and far away from the sources. For the
first, PMUs are at buses 68, 31, and 17 (highlighted with a blue sensor icon
in Fig. 4). For the latter, we place PMUs at buses 23, 41, and 64 (highlighted
with a red sensor icon in Fig. 4). For all realizations, we correctly
determined the input frequencies by thresholding the power spectrum of the
voltage magnitude measurements (near and far case). Fig. 2 shows the spectrum
of an arbitrary measurement realization.

The sensitivity analysis of the parameter $\alpha$ in terms of TPR and FPR for
the two scenarios are shown in Fig. 3 (a) and (b), respectively. Due to space
limitations, we focus only on nearby PMUs scenario. For $0.08\leq\alpha\leq
0.14$, we were able to accurately find true FO locations in most of the
realizations. Finally, for $\alpha=0.14$, we report average and standard
deviation of estimated input parameters in Table I.

###  IV-B Case 2: Reduced WECC 179 bus 29 machine system [18]

Each generator is modeled as a second order classical model and has two
states: rotor angle and angular frequency. We use the small signal analysis
tool (SSAT) to extract state space matrices of the WECC model described in
[18]. We add FO input signals (parameters reported in Table I) to the
mechanical torque of three generators at buses, labeled 5, 14, and 27 in [18,
Fig. 1]. The PMU buses are 4, 18, and 47. Our measurements consist of the
rotor angle of generators [12]. Fig. 3 (c) illustrates the sensitivity
analysis of the parameter $\alpha$ in terms of TPR and FPR. Note that for
$\alpha=0.11$, we found the correct source locations.

##  V Conclusions

This paper studies the effectiveness of complex-LASSO for (i) localizing
sources of multiple forced oscillatory inputs and (ii) estimating the input
parameters. We show that sparsity in the number of locations and sparsity in
the number of sinusoids at a location are clearly manifested in the frequency
domain. This observation led us to cast the localization problem as an
$\ell_{1}$-regularized least squares problem, which we solve numerically via a
complex-valued coordinate descent method. Our localization-estimation
Algorithm 1 is simple and has a potential to be integrated into real-time grid
operations. For simplicity, we present our results assuming the knowledge of
dynamic models; However, our approach is general and can even work with
empirically determined transfer functions.

## References

  * [1] M. Ghorbaniparvar.  Survey on forced oscillations in power system.  J. of Modern Power Systems and Clean Energy, 5(5):671–682, 2017\.
  * [2] D. J. Trudnowski and R. Guttromson.  A strategy for forced oscillation suppression.  IEEE Trans. Power Syst., 35(6):4699–4708, 2020.
  * [3] S. Roy, W. Ju, N. Nayak, and B. Lesieutre.  Localizing power-grid forced oscillations based on harmonic analysis of synchrophasor data.  In 55th Annual Conference on Information Sci. and Sys., pages 1–5, 2021.
  * [4] H. Ye et al.  Analysis and detection of forced oscillation in power system.  IEEE Trans. Power Syst., 32(2):1149–1160, 2017.
  * [5] K. Manohar, J. N. Kutz, and S. L. Brunton.  Optimal sensor and actuator selection using balanced model reduction.  IEEE Trans. Autom. Control, pages 1–1, 2021.
  * [6] B. Wang and K. Sun.  Location methods of oscillation sources in power systems: a survey.  J. of Modern Power Sys. and Clean Energy, 5(2):151–159, 2017.
  * [7] J. A. Tropp and S. J. Wright.  Computational methods for sparse solution of linear inverse problems.  Proceed. of the IEEE, 98(6):948–958, 2010.
  * [8] T. Huang, N. M. Freris, P. R. Kumar, and L. Xie.  A synchrophasor data-driven method for forced oscillation localization under resonance conditions.  IEEE Trans. on Power Systems, 35(5):3927–3939, 2020.
  * [9] R. Anguluri, O. Kosut, and L. Sankar.  Localization and estimation of unknown forced inputs: A group LASSO approach.  IEEE Transactions on Control of Network Systems.  Submitted.
  * [10] J. H. Friedman, T. Hastie, and R. Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  Journal of Statistical Software, 33(1):1–22, 2010.
  * [11] N. Zhou, M. Ghorbaniparvar, and S. Akhlaghi.  Locating sources of forced oscillations using transfer functions.  In 2017 IEEE Power and Energy Conference at Illinois (PECI), pages 1–8. 2017.
  * [12] M. Luan, S. Li, D. Gan, and D. Wu.  Frequency domain approaches to locate forced oscillation source to control device.  International Journal of Electrical Power & Energy Systems, 117:105704, 2020.
  * [13] S. C. Chevalier, P. Vorobev, and K. Turitsyn.  Using effective generator impedance for forced oscillation source location.  IEEE Trans. Power Syst., 33(6):6264–6277, 2018.
  * [14] S. C. Chevalier, P. Vorobev, and K. Turitsyn.  A bayesian approach to forced oscillation source location given uncertain generator parameters.  IEEE Trans. Power Syst., 34(2):1641–1649, 2019.
  * [15] P. Stoica and R. Moses.  Introduction to Spectral Analysis.  Upper Saddle River, NJ: Prentice-Hall, 1997.
  * [16] A. K. Singh et al.  Report on the 68-bus, 16-machine, 5-area system.  IEEE PES Task Force on Benchmark Systems for Stability Controls. Ver, 3, 2013.
  * [17] J.H. Chow and K.W. Cheung.  A toolbox for power system dynamics and control engineering education and research.  IEEE Transactions on Power Systems, 7(4):1559–1564, 1992.
  * [18] S. Maslennikov et al.  A test cases library for methods locating the sources of sustained oscillations.  In 2016 IEEE Power and Energy Society General Meeting (PESGM), pages 1–5, 2016.