DESY 21–013

DO–TH 21/04

TTP21–004

SAGEX–21–03

February 2021

The QED Initial State Corrections to

the Forward-Backward Asymmetry of

$e^{+}e^{-}\rightarrow\gamma^{*}/{Z^{0}}^{*}$ to Higher Orders

J. Blümleina, A. De Freitasa, and K. Schönwaldb

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

b Institut für Theoretische Teilchenphysik,
Karlsruher Institut für Technologie (KIT) D-76128 Karlsruhe, Germany

###### Abstract

The QED initial state corrections are calculated to the forward–backward
asymmetry for $e^{+}e^{-}\rightarrow\gamma^{*}/{Z^{0}}^{*}$ in the leading
logarithmic approximation to $O(\alpha^{6}L^{6})$ extending the known
corrections up to $O(\alpha^{2}L^{2})$ in analytic form. We use the method of
massive on-shell operator matrix elements and present the radiators both in
Mellin-$N$ and momentum fraction $z$-space. Numerical results are presented
for various energies around the $Z$-peak by also including energy cuts. These
corrections are of relevance for the precision measurements at the FCC$\\_$ee.

Dedicated to the Memory of Tini Veltman, who made it possible to probe the
Standard Model at high precision.

##  1 Introduction

The measurement of the forward–backward asymmetry for the process
$e^{+}e^{-}\rightarrow\gamma^{*}/Z^{*}$ provides an excellent possibility to
determine the running fine structure constant $\alpha_{\rm QED}(s)=4\pi a(s)$
near $s=M_{Z}^{2}$ at high precision. At the planned future $e^{+}e^{-}$
facilities which operate at high energy and at large luminosity, like the ILC,
CLIC [1, 2, 3, 4], the FCC_ee [5], and also muon colliders [6], if operating
in the vicinity of the $Z$-peak, one will obtain highly precise data. Any
measurement based on these data needs a theoretical description of even higher
precision [7]. Recently, higher order inclusive corrections were calculated
for $e^{+}e^{-}\rightarrow\gamma^{*}/Z^{*}$ to $O(\alpha_{s}^{6}L^{5})$ up to
the first three orders in $L$ at the respective order in $\alpha$ in Refs. [8,
9, 10, 11] confirming the results to $O(\alpha^{2})$ in Ref. [12] and
correcting Refs. [13, 14], where $L=\ln(s/m_{e}^{2})$. This will lead to a
change of the analysis codes TOPAZ [15, 16] and ZFITTER [17] and may require a
re-analysis of the data taken at LEP [18].

The first order QED initial state radiative (ISR) corrections to the
forward–backward asymmetry have been mutually calculated, see [19, 23, 16, 20,
21, 22, 24]. Furthermore, the initial-final state interference and final state
corrections are known at this order, cf. [16] for a survey. Furthermore,
electroweak [25, 26] and QCD corrections [27] have also been calculated.
Starting at $O(\alpha^{2})$ also the contributions to the leading order series
of $O(\alpha^{k}L^{k}),\leavevmode\nobreak\ k\geq 1$, with
$L=\ln(s/m_{e}^{2})$111The leading order contributions to the direct terms
were obtain in [28]. receive besides the inclusive contribution another one,
related to the angular structure, being present in all sub-leading terms as
well. Yet the $O(\alpha^{k}L^{k})$ terms are universal, since they do not
depend on the process-dependent Wilson coefficients, cf. [11]. The inclusive
terms were computed in [11] and include besides the anomalous dimensions and
massive OMEs also the inclusive massless Wilson coefficients up to
$O(\alpha^{2})$ [29, 30].

In the present paper we will calculate the angular dependent leading
logarithmic contributions to the radiators to $O(\alpha^{6}L^{6})$ as a first
specific contribution which emerges for the forward–backward asymmetry. To 2nd
order, these corrections were obtained in [23]. The corresponding radiators
can be represented by iterated integrals over the alphabet of the harmonic
polylogarithms [31] and cyclotomic harmonic polylogarithms [32] for cyclotomy
c = 4. We also determine efficient representations for these quantities
allowing a fast numerical analysis and present the corresponding corrections
for the forward–backward asymmetry in the vicinity of the $Z$-resonance. These
corrections, unlike the inclusive ones, do not lead to distribution–valued
radiators.

The paper is organized as follows. In Section 2 we calculate the higher order
QED initial state corrections to the forward–backward asymmetry. Here we use
the packages Sigma [33, 34] and HarmonicSums [35, 32, 36, 37, 31, 38]. The
leading logarithmic radiators of $H_{FB}(z)$ to $O(\alpha^{6}L^{6})$ of the
angular-dependent terms are presented in Section 3. We also derive the
expansions of the radiators in the regions $z\rightarrow 0$ and $z\rightarrow
1$. Numerical and phenomenological results are presented in Section 4. In
Appendix A we give a brief account of the aforementioned radiators in Mellin
$N$-space. The radiators are given in computer-readable form in an attachment
to this paper.

##  2 The Forward-Backward Asymmetry around the $Z$-Peak

The forward-backward asymmetry is formed out of the partial cross sections
integrating over the angle $\theta$ for the forward and backward hemispheres
separately,

$\displaystyle\sigma_{F}=2\pi\int\limits_{0}^{1}{\rm d}\cos(\theta)\frac{{\rm
d}\sigma}{{\rm d}\Omega},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \sigma_{B}=2\pi\int\limits_{-1}^{0}{\rm
d}\cos(\theta)\frac{{\rm d}\sigma}{{\rm d}\Omega}\,.$ (1)

The angle $\theta$ is defined between the incoming electron $e^{-}$ and the
outgoing muon $\mu^{-}$ from $\gamma^{*}/Z^{*}$ decay. The forward-backward
asymmetry is defined by

$\displaystyle A_{\text{FB}}(s)$ $\displaystyle=$
$\displaystyle\frac{\sigma_{F}(s)-\sigma_{B}(s)}{\sigma_{T}(s)},$ (2)

with $\sigma_{T}(s)=\sigma_{F}(s)+\sigma_{B}(s)$. At Born level this reduces
to [39]

$\displaystyle\sigma_{FB}^{(0)}(s)$ $\displaystyle=$
$\displaystyle\sigma_{F}^{(0)}(s)-\sigma_{B}^{(0)}(s)=\frac{\pi\alpha^{2}}{s}N_{C,f}\left(1-\frac{4m_{f}^{2}}{s}\right)G_{3}(s)\,,$
(3) $\displaystyle\sigma_{T}^{(0)}(s)$ $\displaystyle=$
$\displaystyle\sigma_{F}^{(0)}(s)+\sigma_{B}^{(0)}(s)=\frac{4\pi\alpha^{2}}{3s}N_{C,f}\sqrt{1-\frac{4m_{f}^{2}}{s}}\left[\left(1+\frac{2m_{f}^{2}}{s}\right)G_{1}(s)-6\frac{m_{f}^{2}}{s}G_{2}(s)\right],$
(4)

with $m_{f}$ the final state fermion mass, $m_{f}\equiv m_{\mu}$. $N_{C,f}$ is
the number of colors of the final state fermion, with $N_{C,f}=1$ in the
present case. $s$ is the cms energy, and the effective couplings
$G_{i}(s)|_{i=1...3}$ read

$\displaystyle G_{1}(s)$ $\displaystyle=$ $\displaystyle
G_{1,1}+G_{1,2}+G_{1,3}=Q_{e}^{2}Q_{f}^{2}+2Q_{e}Q_{f}v_{e}v_{f}{\sf
Re}[\chi_{Z}(s)]+(v_{e}^{2}+a_{e}^{2})(v_{f}^{2}+a_{f}^{2})|\chi_{Z}(s)|^{2},$
(5) $\displaystyle G_{2}(s)$ $\displaystyle=$
$\displaystyle(v_{e}^{2}+a_{e}^{2})a_{f}^{2}|\chi_{Z}(s)|^{2},$ (6)
$\displaystyle G_{3}(s)$ $\displaystyle=$ $\displaystyle
G_{3,1}+G_{3,2}=2Q_{e}Q_{f}a_{e}a_{f}{\sf
Re}[\chi_{Z}(s)]+4v_{e}v_{f}a_{e}a_{f}|\chi_{Z}(s)|^{2}.$ (7)

For later use we define

$\displaystyle\sigma_{T}^{\gamma\gamma}$ $\displaystyle=$ $\displaystyle
F_{1}\leavevmode\nobreak\ G_{1,1},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $ $\displaystyle\sigma_{T}^{\gamma
Z}$ $\displaystyle=$
$\displaystyle\sigma_{T}^{Z\gamma}=\frac{1}{2}F_{1}\leavevmode\nobreak\
G_{1,2},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ $ $\displaystyle\sigma_{T}^{ZZ}$ $\displaystyle=$
$\displaystyle F_{1}\leavevmode\nobreak\ G_{1,3}-6F_{3}\leavevmode\nobreak\
G_{2}$ (8) $\displaystyle\sigma_{FB}^{\gamma Z}$ $\displaystyle=$
$\displaystyle\sigma_{FB}^{Z\gamma}=\frac{1}{2}F_{2}\leavevmode\nobreak\
G_{3,1},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ $ $\displaystyle\sigma_{FB}^{ZZ}$ $\displaystyle=$
$\displaystyle F_{2}\leavevmode\nobreak\ G_{3,2}.$ (9)

with

$\displaystyle F_{1}$ $\displaystyle=$
$\displaystyle\frac{4\pi\alpha^{2}}{3s}N_{C,f}\sqrt{1-\frac{4m_{f}^{2}}{s}}\left(1+\frac{2m_{f}^{2}}{s}\right),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
F_{2}=\frac{\pi\alpha^{2}}{s}N_{C,f}\left(1-\frac{4m_{f}^{2}}{s}\right),$
$\displaystyle F_{3}$ $\displaystyle=$
$\displaystyle\frac{4\pi\alpha^{2}}{3s}N_{C,f}\left(1-\frac{4m_{f}^{2}}{s}\right)\frac{m_{f}^{2}}{s}.$
(10)

The reduced $Z$–propagator is given by

$\displaystyle\chi_{Z}(s)=\frac{s}{s-M_{Z}^{2}+iM_{Z}\Gamma_{Z}},$ (11)

where $M_{Z}$ and $\Gamma_{Z}$ are the mass and the width of the $Z$ boson.
$Q_{e,f}$ are the electromagnetic charges of the electron $(Q_{e}=-1)$ and the
final state fermion, respectively, and the electro–weak couplings $v_{i}$ and
$a_{i}$ read

$\displaystyle v_{e}$
$\displaystyle=\frac{1}{\sin\theta_{w}\cos\theta_{w}}\left[I^{3}_{w,e}-2Q_{e}\sin^{2}\theta_{w}\right],\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $
$\displaystyle a_{e}$
$\displaystyle=\frac{1}{\sin\theta_{w}\cos\theta_{w}}I^{3}_{w,e},$ (12)
$\displaystyle v_{f}$
$\displaystyle=\frac{1}{\sin\theta_{w}\cos\theta_{w}}\left[I^{3}_{w,f}-2Q_{f}\sin^{2}\theta_{w}\right],\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $
$\displaystyle a_{f}$
$\displaystyle=\frac{1}{\sin\theta_{w}\cos\theta_{w}}I^{3}_{w,f}\leavevmode\nobreak\
.$ (13)

$\theta_{w}$ denotes the weak mixing angle, and $I^{3}_{w,i}=\pm 1/2$ the
third component of the weak isospin for up and down particles, respectively.

When accounting for initial-state-radiation, one obtains the following
representation of the $A_{FB}$ introducing the radiators $H_{e}^{LL}$ and
$H_{FB}^{LL}$ [23], using the notation in [19],

$\displaystyle A_{\text{FB}}(s)$ $\displaystyle=$
$\displaystyle\frac{1}{\sigma_{T}(s)}\int\limits_{z_{0}}^{1}{\rm
d}z\,\frac{4z}{(1+z)^{2}}\tilde{H}_{e}^{LL}(z)\sigma_{FB}^{(0)}(zs),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\sigma_{T}(s)=\int\limits_{z_{0}}^{1}{\rm
d}z\,H_{e}(z)\sigma_{T}^{(0)}(zs)\,,$

with

$\displaystyle\tilde{H}_{e}^{LL}(z)=\left[H_{e}^{LL}(z)+H_{FB}^{LL}(z)\right].$
(15)

The normalization factor $\sigma_{T}(s)$ will be calculated considering all
corrections derived in Ref. [11]. We will consider different values for the
threshold $z_{0}$. Like in Refs. [8, 9, 10] $z_{0}$ is chosen as
$z_{0}=4m_{\tau}^{2}/s$, with $m_{\tau}$ the mass of the $\tau$ lepton.
Another choice will be $z_{0}=0.99$ or $0.999$ in accordance with [40]. The
variable $z$ is given by $z=s^{\prime}/s$, where $s^{\prime}$ denotes the
virtuality of the gauge boson $\gamma^{*}$ or $Z^{*}$. Furthermore, $H_{e}(z)$
is the radiator in the inclusive case, given in Ref. [11] and $H_{FB}^{LL}(z)$
denotes the leading-log radiator in the angular dependent case up to
$O(\alpha^{6}L^{6})$,

$\displaystyle H_{FB}^{LL}(z)=\int\limits_{0}^{1}{\rm
d}x_{1}\int\limits_{0}^{1}{\rm
d}x_{2}\left(\frac{(1+z)^{2}}{(x_{1}+x_{2})^{2}}-1\right)\Gamma_{ee}^{LL}(x_{2})\Gamma_{ee}^{LL}(x_{1})\delta(x_{1}x_{2}-z).$
(16)

Here $\Gamma_{ee}^{LL}(x)$ is the leading-log operator matrix element. The
product $\Gamma_{ee}^{LL}(x_{1})\Gamma_{ee}^{LL}(x_{2})$ is consistently
expanded up to $O(\alpha^{6}L^{6})$. The latter radiator does only contain
universal contributions. The radiators obey the expansion

$\displaystyle H_{i}(z)$ $\displaystyle=$
$\displaystyle\delta_{ie}\delta(1-z)+\sum_{k=1}^{\infty}a^{k}\sum_{l=0}^{k}L^{l}H^{(k,l)}_{i}(z),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ i=e,FB,$ (17) $\displaystyle
H_{i}^{LL}(z)$ $\displaystyle=$
$\displaystyle\delta_{ie}\delta(1-z)+\sum_{k=1}^{\infty}(aL)^{k}H^{(k),LL}_{i}(z),$
(18)

where $a=\alpha/(4\pi)$ and $\alpha$ denotes the fine structure constant.

We proceed in the following way in order to evaluate Eq. (16) analytically.
First, we calculate the Mellin transform

$\displaystyle\mathcal{M}[H_{FB}^{LL}(z)](n)$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{1}{\rm
d}zz^{n}H_{FB}^{LL}(z)=\int\limits_{0}^{1}{\rm d}x_{1}\int\limits_{0}^{1}{\rm
d}x_{2}x_{1}^{n}x_{2}^{n}\left(\frac{(1+x_{1}x_{2})^{2}}{(x_{1}+x_{2})^{2}}-1\right)\Gamma_{ee}^{LL}(x_{2})\Gamma_{ee}^{LL}(x_{1})\leavevmode\nobreak\
.$

Since this integral is not suited to be integrated with the package
HarmonicSums directly, we compute the generating function

$\displaystyle\mathcal{G}[H_{FB}^{LL}(z)](t)$ $\displaystyle=$
$\displaystyle\sum\limits_{n=0}^{\infty}t^{n}\mathcal{M}[H_{FB}^{LL}(z)](n)$
(20) $\displaystyle=$ $\displaystyle\int\limits_{0}^{1}{\rm
d}x_{1}\int\limits_{0}^{1}{\rm
d}x_{2}\frac{1}{1-tx_{1}x_{2}}\left(\frac{(1+x_{1}x_{2})^{2}}{(x_{1}+x_{2})^{2}}-1\right)\Gamma_{ee}^{LL}(x_{2})\Gamma_{ee}^{LL}(x_{1})\leavevmode\nobreak\
,$

which resums the Mellin-kernel into a denominator which can be easily
integrated over. After the integration over $x_{1}$ and $x_{2}$ we are left
with generalized iterated integrals evaluated at argument $z=1$ which contain
the parameter $t$ in their letters. We use differential equations to pull this
parameter into the argument. This is straightforward since the limit $t\to 0$
always exists and is easily expressed in terms of known constants. Afterwards
we can use the HarmonicSums command GetMoment to get the $n$-space expression
and GeneralInvMellin to arrive at the final result. The radiators are
consistently expanded in $(aL)$ to $O((aL)^{6})$ and are expressed using the
variable $\sqrt{z}$ rather than $z$ to obtain a unique representation
concerning the contributing iterated integrals.

At $O(\alpha)$ the complete QED initial state corrections [21, 22, 16] are
known and are given by

$\displaystyle A_{\text{FB}}^{(1)}(s)$ $\displaystyle=$
$\displaystyle\frac{1}{\sigma_{T}(s)}\int\limits_{z_{0}}^{1}{\rm
d}z\,\frac{4z}{(1+z)^{2}}a(s)[\tilde{H}_{e}^{(1),LL}(z)L+\tilde{H}_{e}^{(1,0)}(z)]\sigma_{T}^{(0)}(zs)\,,$
(21)

and

$\displaystyle{H}_{e}^{(1),LL}(z)$ $\displaystyle=$ $\displaystyle
4\left[\frac{1+z^{2}}{1-z}\right]_{+}$ (22)
$\displaystyle\tilde{H}_{e}^{(1,0)}(z)$ $\displaystyle=$ $\displaystyle
4\left[-\left[\frac{1+z^{2}}{1-z}\right]_{+}+\left(2\zeta_{2}-\frac{1}{2}\right)\delta(1-z)+\frac{1+z^{2}}{1-z}\left[2\ln(1+z)-2\ln(2)-\ln(z)\right]\right].$

##  3 The leading-log radiators $H_{FB}^{LL}$

We obtain the following leading-log radiators $H_{FB}^{(k),LL}$ :

$\displaystyle H_{FB}^{(1),LL}(z)$ $\displaystyle=$ $\displaystyle 0$ (24)
$\displaystyle H_{FB}^{(2),LL}(z)$ $\displaystyle=$
$\displaystyle\frac{2(1-z)(1+z)^{2}}{z}+2\pi\frac{(1-z)^{2}}{\sqrt{z}}-8(1+z)H_{0}-8(1-z)^{2}\frac{H_{\\{4,0\\}}}{\sqrt{z}}$
(25) $\displaystyle H_{FB}^{(3),LL}(z)$ $\displaystyle=$
$\displaystyle-\frac{16(1-z)\big{(}4+11z+4z^{2}\big{)}}{3z}-\pi\biggl{[}\frac{4\big{(}2-3z-2z^{2}-3z^{3}+2z^{4}\big{)}}{3z^{3/2}}$
(26)
$\displaystyle+\frac{4(1-z)(1+5z)}{\sqrt{z}}H_{0}+\frac{16(1-z)^{2}}{\sqrt{z}}H_{\\{4,1\\}}\biggr{]}+\biggl{[}\frac{4(1+z)\big{(}5-18z-19z^{2}\big{)}}{3z}$
$\displaystyle-\frac{16(1-z)(1-7z)}{\sqrt{z}}H_{\\{4,0\\}}-96(1+z)H_{\\{4,1\\}}\biggr{]}H_{0}-8(1+z)H_{0}^{2}$
$\displaystyle+\biggl{[}\frac{16(1-z)(1+z)^{2}}{z}-\frac{64(1-z)^{2}}{\sqrt{z}}H_{\\{4,0\\}}\biggr{]}H_{1}+\biggl{[}\frac{16\big{(}2-3z-2z^{2}-3z^{3}+2z^{4}\big{)}}{3z^{3/2}}$
$\displaystyle+\frac{64(1-z)^{2}}{\sqrt{z}}H_{\\{4,1\\}}\biggr{]}H_{\\{4,0\\}}+\biggl{[}-\frac{16(1-z)(1+z)^{2}}{z}+64(1+z)H_{0}$
$\displaystyle+\frac{64(1-z)^{2}}{\sqrt{z}}H_{\\{4,0\\}}\biggr{]}H_{-1}-64(1+z)H_{0,1}+\frac{32(1-z)^{2}}{\sqrt{z}}H_{0,\\{4,0\\}}+96(1+z)H_{0,\\{4,1\\}}$
$\displaystyle+\frac{64(1-z)^{2}}{\sqrt{z}}H_{1,\\{4,0\\}}-\frac{64(1-z)^{2}}{\sqrt{z}}H_{\\{4,0\\},\\{4,1\\}}-64(1+z)H_{-1,0}$
$\displaystyle-\frac{64(1-z)^{2}}{\sqrt{z}}H_{-1,\\{4,0\\}}+20(1+z)\zeta_{2},$

with $H_{\vec{w}}\equiv H_{\vec{w}}(\sqrt{z})$ and
$\zeta_{k},k\in\mathbb{N},k\geq 2$ are the values of Riemann’s $\zeta$
function at integer argument. The iterated integrals are given by

$\displaystyle
H_{b,\vec{a}}(z)=\int_{0}^{z}dyf_{b}(y)H_{\vec{a}}(y),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ H_{\emptyset}=1.$ (27)

The indices $b,a_{i}$ refer to the letters in the following alphabet

$\displaystyle\mathfrak{A}=\left\\{f_{0}=\frac{1}{z},f_{1}=\frac{1}{1-z},f_{-1}=\frac{1}{1+z},f_{\\{4,0\\}}=\frac{1}{1+z^{2}},f_{\\{4,1\\}}=\frac{z}{1+z^{2}}\right\\},$
(28)

where the last two letters in (28) are cyclotomic. The function
$H_{FB}^{(2),LL}(z)$ has been calculated in [23] in an equivalent
representation. The radiators $H_{FB}^{(4-6),LL}(z)$ are too voluminous to be
presented here and are given in computer-readable form in the attachment.

We finally illustrate the structure of the radiators expanding them for their
leading small and a few large $z$ terms. One obtains

$\displaystyle zH_{FB}^{(2),LL}(z)$ $\displaystyle\simeq$ $\displaystyle
2+O(\sqrt{z})$ (29) $\displaystyle zH_{FB}^{(3),LL}(z)$ $\displaystyle\simeq$
$\displaystyle-\frac{8\pi}{3\sqrt{z}}+\frac{10}{3}\ln(z)-\frac{32}{3}+O(\sqrt{z}\ln(z))$
(30) $\displaystyle zH_{FB}^{(4),LL}(z)$ $\displaystyle\simeq$
$\displaystyle\frac{32}{9z}-\frac{11}{6}\ln^{2}(z)+\frac{10}{27}\ln(z)-38-43\zeta_{2}$
(31)
$\displaystyle+\frac{\pi}{\sqrt{z}}\biggl{(}\frac{4}{3}\ln(z)+\frac{328}{27}\biggr{)}+O(\sqrt{z}\ln^{2}(z))$
$\displaystyle zH_{FB}^{(5),LL}(z)$ $\displaystyle\simeq$
$\displaystyle-\frac{512}{81z}+\frac{5}{36}\ln^{3}(z)+\frac{97}{18}\ln^{2}(z)+\biggl{[}-\frac{299}{27}+\frac{215\zeta_{2}}{18}\biggr{]}\ln(z)+\frac{72283}{243}$
(32)
$\displaystyle+\frac{2551}{27}\zeta_{2}+\frac{320}{9}\zeta_{3}+\frac{\pi}{\sqrt{z}}\biggl{(}\frac{2}{3}\ln^{2}(z)+\frac{92}{27}\ln(z)+\frac{460}{81}$
$\displaystyle+\frac{4}{9}\zeta_{2}\biggr{)}+O(\sqrt{z}\ln^{3}(z))$
$\displaystyle zH_{FB}^{(6),LL}(z)$ $\displaystyle\simeq$
$\displaystyle-\frac{10240}{729z}-\frac{256}{81z}\ln(z)+\frac{97}{720}\ln^{4}(z)+\frac{79}{162}\ln^{3}(z)+\biggl{(}\frac{17059}{2430}-\frac{19}{180}\zeta_{2}\biggr{)}\ln^{2}(z)$
(33)
$\displaystyle+\biggl{(}-\frac{159109}{7290}-\frac{2825}{81}\zeta_{2}+\frac{136}{9}\zeta_{3}\biggr{)}\ln(z)-\frac{2459744}{10935}+\frac{87173}{1215}\zeta_{2}-\frac{53216}{405}\zeta_{3}$
$\displaystyle+\frac{5659}{300}\zeta_{2}^{2}+\frac{\pi}{\sqrt{z}}\biggl{(}-\frac{1}{18}\ln^{3}(z)-\frac{32}{9}\ln^{2}(z)+\biggl{[}-\frac{39212}{1215}-\frac{73}{9}\zeta_{2}\biggr{]}\ln(z)$
$\displaystyle-\frac{401806}{3645}-\frac{15904}{405}\zeta_{2}-\frac{128}{9}\zeta_{3}\biggr{)}+O(\sqrt{z}\ln^{4}(z)).$

For the large $z$ representation we set $u=1-\sqrt{z}$ and obtain

$\displaystyle H_{FB}^{(2),LL}(z)$ $\displaystyle\simeq$ $\displaystyle
32u-16u^{2}+\frac{112}{3}u^{3}+\frac{56}{3}u^{4}+O(u^{5})$ (34) $\displaystyle
H_{FB}^{(3),LL}(z)$ $\displaystyle\simeq$
$\displaystyle\frac{64}{3}u-\frac{416}{3}u^{2}+\frac{1408}{9}u^{3}-\frac{736}{9}u^{4}+\Biggl{(}-256u+128u^{2}-\frac{896}{3}u^{3}-\frac{448}{3}u^{4}\Biggr{)}$
(35) $\displaystyle\times\ln(2u)+O(u^{5})$ $\displaystyle H_{FB}^{(4),LL}(z)$
$\displaystyle\simeq$
$\displaystyle\frac{27328}{27}u-\frac{13664}{27}u^{2}+\frac{162080}{81}u^{3}+\frac{36112}{81}u^{4}+\biggl{(}1024u-512u^{2}$
(36)
$\displaystyle+\frac{3584}{3}u^{3}+\frac{1792}{3}u^{4}\biggr{)}\left(\ln^{2}(2u)-\zeta_{2}\right)+\biggl{(}1024-\frac{3328}{3}u+\frac{2176}{3}u^{2}\biggr{)}$
$\displaystyle\times u^{2}\ln(2u)+O(u^{5})$ $\displaystyle H_{FB}^{(5),LL}(z)$
$\displaystyle\simeq$
$\displaystyle\frac{387712}{81}u-\frac{523328}{81}u^{2}+\frac{2954144}{243}u^{3}-\frac{328688}{243}u^{4}+\biggl{(}\frac{2048}{3}+\frac{11264}{3}u$
(37)
$\displaystyle-\frac{34304}{9}u^{2}+\frac{28928}{9}u^{3}\biggr{)}u\left(\zeta_{2}-\ln^{2}(2u)\right)+\biggl{(}-\frac{8192}{3}+\frac{4096}{3}u-\frac{28672}{9}u^{2}$
$\displaystyle-\frac{14336}{9}u^{3}\biggr{)}u\left(2\zeta_{3}-3\zeta_{2}\ln(2u)+\ln^{3}(2u)\right)+\biggl{(}-\frac{73216}{9}+\frac{30464}{9}u$
$\displaystyle-\frac{419840}{27}u^{2}-\frac{108544}{27}u^{3}\biggr{)}u\ln(2u)+O(u^{5})$
$\displaystyle H_{FB}^{(6),LL}(z)$ $\displaystyle\simeq$
$\displaystyle\frac{153977536}{3645}u-\frac{136846688}{3645}u^{2}+\frac{1140995264}{10935}u^{3}+\frac{60397216}{10935}u^{4}+\biggl{(}-\frac{894976}{27}$
(38)
$\displaystyle+\frac{300032}{27}u-\frac{24477184}{405}u^{2}-\frac{7326464}{405}u^{3}\biggr{)}u\left(\zeta_{2}-\ln^{2}(2u)\right)+\biggl{(}\frac{16384}{5}-\frac{8192}{5}u$
$\displaystyle+\frac{57344}{15}u^{2}+\frac{28672}{15}u^{3}\biggr{)}u\left(\zeta_{2}^{2}+\frac{5}{3}\ln^{4}(2u)\right)+\biggl{(}\frac{32768}{9}+\frac{81920}{9}u-\frac{1132544}{135}u^{2}$
$\displaystyle+\frac{1276928}{135}u^{3}\biggr{)}u\left(2\zeta_{3}-3\zeta_{2}\ln(2u)+\ln^{3}(2u)\right)+\biggl{(}-\frac{886784}{27}+\frac{1338368}{27}u$
$\displaystyle-\frac{35939456}{405}u^{2}+\frac{5414336}{405}u^{3}-\biggl{[}-\frac{131072}{3}+\frac{65536}{3}u-\frac{458752}{9}u^{2}$
$\displaystyle-\frac{229376}{9}u^{3}\biggr{]}\zeta_{3}\biggr{)}u\ln(2u)+\biggl{(}-32768+16384u-\frac{114688}{3}u^{2}-\frac{57344}{3}u^{3}\biggr{)}$
$\displaystyle\times u\zeta_{2}\ln^{2}(2u)+O(u^{5}).$

These terms are the beginning of the respective series used in the numerical
representation to which we turn now.

##  4 Numerical results

For the numerical evaluation of the new radiators we calculated 60 terms of
the expansions around $\sqrt{z}=0$ and $\sqrt{z}=1$ using HarmonicSums. The
constants associated with cyclotomic harmonic polylogarithms for cyclotomy
c=4, which are needed for the expansion around $\sqrt{z}=1$ are known at least
up to weight w = 5 [32]. However, all constants except of $\pi$ cancel in the
final results. We switch between these two expansions at $z=1/4$, where the
absolute difference of the respective approximation is always smaller than
$10^{-12}$. The usual harmonic polylogarithms are calculated using the
implementations in [42, 43]. In the numerical illustration we refer to the
electroweak parameters given in [41].

In Figure 1 we illustrate the effect of the initial–state QED corrections on
the forward–backward asymmetry for a wider range of cms energies around the
$Z$-peak with $\sqrt{s}\in[85,96]\leavevmode\nobreak\ {\rm GeV}$ to provide an
overall impression on the size of the radiative corrections, which will be
detailed below.

Let us define

$\displaystyle\Delta A_{FB}=1-\frac{A_{FB}^{(l)}}{A_{FB}^{(0)}},$ (39)

Figure 1: $A_{FB}$ and its initial state QED corrections as a function of
$\sqrt{s}$. Black (dashed) the Born approximation, blue (dotted) the
$O(\alpha)$ improved approximation, red (full) also including the leading-log
improvement up to $O(\alpha^{6})$ for $s^{\prime}/s\geq 4m_{\tau}^{2}/s$.
Figure 2: $\Delta A_{FB}$ in % as defined in Eq. (39) as a function of
$\sqrt{s}$. Black (dashed) the $O(\alpha)$ improved approximation, blue
(dotted) the $O(\alpha^{2}L^{2})$ improved approximation, red (full) also
including the leading-log improvement up to $O(\alpha^{6})$ for
$s^{\prime}/s\geq 4m_{\tau}^{2}/s$.

where $l$ denotes the order to which the initial state radiation is taken into
account. $l=0$ corresponds to the Born approximation. In Figure 2 we
illustrate the relative size of the ISR corrections to the Born cross section.

In Table 1 we illustrate the effect of the different orders of the QED initial
state corrections of the forward–backward asymmetry at the $Z$-peak and two
more values of $s=s_{\pm}$ [40]. The ISR corrections are not monotonic order
by order, even leading to a sign change in the corrected values of $A_{FB}$ at
$s=M_{Z}^{2}$. Here the ISR corrections are largest.

If we compare the known one-loop corrected values of $A_{FB}$ with the highest
radiative correction calculated in this paper we obtain further corrections of
$-3$ % for $s_{-}$ and $-1\%$ for $s_{+}$. These values supersede previous
estimates, based on assumptions, in [40].

| $A_{FB}(s_{-})$ | $A_{FB}(M_{Z}^{2})$ | $A_{FB}(s_{+})$
---|---|---|---
$\mathcal{O}(\alpha^{0})$ | $-0.3564803$ | $0.0225199$ | $0.2052045$
$\mathcal{O}(\alpha L^{1})$ | $-0.2945381$ | $-0.0094232$ | $0.1579347$
$\mathcal{O}(\alpha L^{0})$ | $-0.2994478$ | $-0.0079610$ | $0.1611962$
$\mathcal{O}(\alpha^{2}L^{2})$ | $-0.3088363$ | $0.0014514$ | $0.1616887$
$\mathcal{O}(\alpha^{3}L^{3})$ | $-0.3080578$ | $0.0000198$ | $0.1627252$
$\mathcal{O}(\alpha^{4}L^{4})$ | $-0.3080976$ | $0.0001587$ | $0.1625835$
$\mathcal{O}(\alpha^{5}L^{5})$ | $-0.3080960$ | $0.0001495$ | $0.1625911$
$\mathcal{O}(\alpha^{6}L^{6})$ | $-0.3080960$ | $0.0001499$ | $0.1625911$
Table 1: $A_{FB}$ evaluated at $s_{-}=(87.9\,{\rm GeV})^{2}$, $M_{Z}^{2}$ and
$s_{+}=(94.3\,{\rm GeV})^{2}$ for the cut $z>4m_{\tau}^{2}/s$.

If cuts of $s^{\prime}$ are applied the ISR corrections turn out to be
smaller, as illustrated in Tables 2 and 3. As noted in [19] also an
approximative treatment of the forward–backward asymmetry

$\displaystyle A_{\text{FB}}(s)$ $\displaystyle\approx$
$\displaystyle\frac{1}{\sigma_{T}(s)}\int\limits_{z_{0}}^{1}{\rm
d}z\,\frac{4z}{(1+z)^{2}}H_{e}(z)\sigma_{FB}^{(0)}(zs)\,,$ (40)

has been used in the literature [24] close to the $Z$-peak, where $H_{e}(z)$
denotes the inclusive

| $A_{FB}(s_{-})$ | $A_{FB}(M_{Z}^{2})$ | $A_{FB}(s_{+})$
---|---|---|---
$\mathcal{O}(\alpha^{0})$ | $-0.2714847$ | $0.01545540$ | $0.2539839$
$\mathcal{O}(\alpha L^{1})$ | $-0.2913715$ | $0.01244012$ | $0.2787800$
$\mathcal{O}(\alpha L^{0})$ | $-0.2923065$ | $0.01266498$ | $0.2794039$
$\mathcal{O}(\alpha^{2}L^{2})$ | $-0.2922610$ | $0.01489549$ | $0.2760952$
$\mathcal{O}(\alpha^{3}L^{3})$ | $-0.2917702$ | $0.01430419$ | $0.2764488$
$\mathcal{O}(\alpha^{4}L^{4})$ | $-0.2918380$ | $0.01439578$ | $0.2763865$
$\mathcal{O}(\alpha^{5}L^{5})$ | $-0.2918313$ | $0.01438608$ | $0.2763934$
$\mathcal{O}(\alpha^{6}L^{6})$ | $-0.2918318$ | $0.01438681$ | $0.2763929$
Table 2: $A_{FB}$ evaluated at $s_{-}=(87.9\,{\rm GeV})^{2}$, $M_{Z}^{2}$ and
$s_{+}=(94.3\,{\rm GeV})^{2}$ for the cut $z>0.99$.

radiator appearing in case of $\sigma_{T}(s)$ [11]. In the sub–leading
corrections to the radiators [11], Wilson coefficients appear, which are
different for $\sigma_{T}$ and $\sigma_{FB}$. Therefore the radiator for
$\sigma_{T}$ cannot be used in its sub–leading terms for this reason. The
numerical results in Table 4–6 show, that contributions due to the angular
radiators are indeed suppressed if measuring at the $Z$-peak w.r.t. to the
inclusive ones also in higher orders. This holds quite irrespectively of the
cuts in $s^{\prime}$.

| $A_{FB}(s_{-})$ | $A_{FB}(M_{Z}^{2})$ | $A_{FB}(s_{+})$
---|---|---|---
$\mathcal{O}(\alpha^{0})$ | $-0.2683019$ | $0.01538044$ | $0.2576633$
$\mathcal{O}(\alpha L^{1})$ | $-0.2905494$ | $0.01623701$ | $0.2797643$
$\mathcal{O}(\alpha L^{0})$ | $-0.2913658$ | $0.01630129$ | $0.2805232$
$\mathcal{O}(\alpha^{2}L^{2})$ | $-0.2899912$ | $0.01655798$ | $0.2787005$
$\mathcal{O}(\alpha^{3}L^{3})$ | $-0.2898133$ | $0.01641903$ | $0.2787137$
$\mathcal{O}(\alpha^{4}L^{4})$ | $-0.2898401$ | $0.01645212$ | $0.2786941$
$\mathcal{O}(\alpha^{5}L^{5})$ | $-0.2898362$ | $0.01644636$ | $0.2786984$
$\mathcal{O}(\alpha^{6}L^{6})$ | $-0.2898368$ | $0.01644713$ | $0.2786978$
Table 3: $A_{FB}$ evaluated at $s_{-}=(87.9\,{\rm GeV})^{2}$, $M_{Z}^{2}$ and $s_{+}=(94.3\,{\rm GeV})^{2}$ for the cut $z>0.999$. | all $(s_{-})$ | ang. $(s_{-})$ | all $(M_{Z}^{2})$ | ang. $(M_{Z}^{2})$ | all $(s_{+})$ | ang. $(s_{+})$
---|---|---|---|---|---|---
$\mathcal{O}(\alpha^{2}L^{2})$ | $-9\cdot 10^{-3}$ | $-9\cdot 10^{-5}$ | $9\cdot 10^{-3}$ | $-1\cdot 10^{-5}$ | $5\cdot 10^{-4}$ | $-6\cdot 10^{-5}$
$\mathcal{O}(\alpha^{3}L^{3})$ | $8\cdot 10^{-4}$ | $-4\cdot 10^{-6}$ | $-1\cdot 10^{-3}$ | $-8\cdot 10^{-7}$ | $1\cdot 10^{-3}$ | $-3\cdot 10^{-6}$
$\mathcal{O}(\alpha^{4}L^{4})$ | $-4\cdot 10^{-5}$ | $-3\cdot 10^{-7}$ | $1\cdot 10^{-4}$ | $-7\cdot 10^{-8}$ | $-1\cdot 10^{-4}$ | $-7\cdot 10^{-8}$
$\mathcal{O}(\alpha^{5}L^{5})$ | $2\cdot 10^{-6}$ | $-1\cdot 10^{-8}$ | $-9\cdot 10^{-6}$ | $-4\cdot 10^{-9}$ | $8\cdot 10^{-6}$ | $1\cdot 10^{-8}$
$\mathcal{O}(\alpha^{6}L^{6})$ | $-6\cdot 10^{-8}$ | $-6\cdot 10^{-10}$ | $4\cdot 10^{-7}$ | $-1\cdot 10^{-10}$ | $-3\cdot 10^{-8}$ | $2\cdot 10^{-9}$
Table 4: Contribution to $A_{FB}$ at different orders of $\alpha L$ evaluated at $s_{-}=(87.9\,{\rm GeV})^{2}$, $M_{Z}^{2}$ and $s_{+}=(94.3\,{\rm GeV})^{2}$ for the cut $z>4m_{\tau}^{2}/s$. all: full contributions, ang.: only the angular dependent contributions. | all $(s_{-})$ | ang. $(s_{-})$ | all $(M_{Z}^{2})$ | ang. $(M_{Z}^{2})$ | all $(s_{+})$ | ang. $(s_{+})$
---|---|---|---|---|---|---
$\mathcal{O}(\alpha^{2}L^{2})$ | $5\cdot 10^{-5}$ | $-5\cdot 10^{-8}$ | $2\cdot 10^{-3}$ | $-2\cdot 10^{-9}$ | $-3\cdot 10^{-3}$ | $5\cdot 10^{-8}$
$\mathcal{O}(\alpha^{3}L^{3})$ | $5\cdot 10^{-4}$ | $-3\cdot 10^{-8}$ | $-6\cdot 10^{-4}$ | $-1\cdot 10^{-9}$ | $4\cdot 10^{-4}$ | $3\cdot 10^{-8}$
$\mathcal{O}(\alpha^{4}L^{4})$ | $-7\cdot 10^{-5}$ | $-9\cdot 10^{-9}$ | $9\cdot 10^{-5}$ | $-3\cdot 10^{-10}$ | $-6\cdot 10^{-5}$ | $9\cdot 10^{-9}$
$\mathcal{O}(\alpha^{5}L^{5})$ | $7\cdot 10^{-6}$ | $-2\cdot 10^{-9}$ | $-1\cdot 10^{-5}$ | $-4\cdot 10^{-11}$ | $7\cdot 10^{-6}$ | $2\cdot 10^{-9}$
$\mathcal{O}(\alpha^{6}L^{6})$ | $-5\cdot 10^{-7}$ | $-2\cdot 10^{-10}$ | $7\cdot 10^{-7}$ | $-4\cdot 10^{-12}$ | $-5\cdot 10^{-7}$ | $2\cdot 10^{-10}$
Table 5: Contribution to $A_{FB}$ at different orders of $\alpha L$ evaluated
at $s_{-}=(87.9\,{\rm GeV})^{2}$, $M_{Z}^{2}$ and $s_{+}=(94.3\,{\rm
GeV})^{2}$ for the cut $z>0.99$. all: full contributions, ang.: only the
angular dependent contributions.

For completeness we also quantify the final state and initial–final
interference contributions at $O(\alpha)$, cf. e.g. [19, 16, 21, 22] for the
respective numerator, $\sigma_{FB}$, and denominator, $\sigma_{T}$
contributions to $A_{FB}=\sigma_{FB}/\sigma_{T}$. The inclusive final state
correction to $\sigma_{FB}$ and $\sigma_{T}$ imply the factors

$\displaystyle 1+\delta_{FB}^{F},\leavevmode\nobreak\ \leavevmode\nobreak\
\text{with}\leavevmode\nobreak\ \leavevmode\nobreak\
\delta_{FB}^{F}=0\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
1+\delta_{T}^{F},\leavevmode\nobreak\ \leavevmode\nobreak\
\text{with}\leavevmode\nobreak\ \leavevmode\nobreak\ \delta_{T}^{F}=3a$ (41)

for the $\mu^{+}\mu^{-}$ final state. The multiplicative numerical correction
to $A_{FB}$ is $0.998149$.

The initial–final interference term [16, 22] has a more involved
representation compared to the final state corrections, since there are
different correction terms to the different pieces of the angular averaged
Born cross section. In particular these terms are important because of
logarithmically large terms of $\ln(1-s^{\prime}_{\rm min}/s)$ in the region
$s^{\prime}_{\rm min}$ close to $s$ as chosen in [40]. The correction term
$1+\delta_{FB}^{FI}$ modifying $A_{FB}$ to $O(\alpha)$ is given by

$\displaystyle\delta_{FB}^{FI}$ $\displaystyle=$
$\displaystyle-\frac{1}{\sigma_{T}^{\rm Born}}\Biggl{\\{}\sigma_{FB}^{\gamma
Z}{\sf Re}[R_{T}^{\gamma\gamma}+R_{T}^{ZZ}]+\sigma_{FB}^{ZZ}{\sf
Re}[R_{T}^{ZZ}]\Biggr{\\}}$ (42)
$\displaystyle+\frac{1}{A_{FB}}\Biggl{\\{}{\sf
Re}[R_{FB}^{\gamma\gamma}]\frac{\sigma_{T}^{\gamma\gamma}}{\sigma_{T}^{\rm
Born}}+{\sf Re}[R_{FB}^{\gamma\gamma}+R_{FB}^{ZZ}]\frac{\sigma_{T}^{\gamma
Z}}{\sigma_{T}^{\rm Born}}+{\sf
Re}[R_{FB}^{ZZ}]\frac{\sigma_{T}^{ZZ}}{\sigma_{T}^{\rm Born}}\Biggr{\\}},$

with $R_{A}^{ij},\leavevmode\nobreak\ \leavevmode\nobreak\
A=T,FB,\leavevmode\nobreak\ \leavevmode\nobreak\ i,j=\gamma,Z$, the
corresponding radiators [16, 22, 26] and

$\displaystyle\delta_{FB}^{FI}$ $\displaystyle=$
$\displaystyle-0.0821\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ \leavevmode\nobreak\
\sqrt{s}=87.9\leavevmode\nobreak\ {\rm GeV}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
s^{\prime}/s>0.99$ (43) $\displaystyle\delta_{FB}^{FI}$ $\displaystyle=$
$\displaystyle+0.0928\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ \leavevmode\nobreak\
\sqrt{s}=94.3\leavevmode\nobreak\ {\rm GeV}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
s^{\prime}/s>0.99$ (44) $\displaystyle\delta_{FB}^{FI}$ $\displaystyle=$
$\displaystyle-0.1954\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ \leavevmode\nobreak\
\sqrt{s}=87.9\leavevmode\nobreak\ {\rm GeV}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
s^{\prime}/s>0.999$ (45) $\displaystyle\delta_{FB}^{FI}$ $\displaystyle=$
$\displaystyle+0.2223\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ \leavevmode\nobreak\
\sqrt{s}=94.3\leavevmode\nobreak\ {\rm GeV}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
s^{\prime}/s>0.999$ (46)

The tight cuts in (43–46) imply a correction term of up to $\pm 9\%$ for
$z>0.99$ and up to $\pm 22\%$ for $z>0.999$.

| all $(s_{-})$ | ang. $(s_{-})$ | all $(M_{Z}^{2})$ | ang. $(M_{Z}^{2})$ | all $(s_{+})$ | ang. $(s_{+})$
---|---|---|---|---|---|---
$\mathcal{O}(\alpha^{2}L^{2})$ | $1\cdot 10^{-3}$ | $-5\cdot 10^{-10}$ | $3\cdot 10^{-4}$ | $2\cdot 10^{-11}$ | $-2\cdot 10^{-3}$ | $5\cdot 10^{-10}$
$\mathcal{O}(\alpha^{3}L^{3})$ | $2\cdot 10^{-4}$ | $-4\cdot 10^{-10}$ | $-1\cdot 10^{-4}$ | $2\cdot 10^{-11}$ | $1\cdot 10^{-5}$ | $4\cdot 10^{-10}$
$\mathcal{O}(\alpha^{4}L^{4})$ | $-3\cdot 10^{-5}$ | $-2\cdot 10^{-10}$ | $3\cdot 10^{-5}$ | $9\cdot 10^{-12}$ | $-2\cdot 10^{-5}$ | $2\cdot 10^{-10}$
$\mathcal{O}(\alpha^{5}L^{5})$ | $4\cdot 10^{-6}$ | $-5\cdot 10^{-11}$ | $-6\cdot 10^{-6}$ | $3\cdot 10^{-12}$ | $4\cdot 10^{-6}$ | $5\cdot 10^{-11}$
$\mathcal{O}(\alpha^{6}L^{6})$ | $-5\cdot 10^{-7}$ | $-1\cdot 10^{-11}$ | $8\cdot 10^{-7}$ | $5\cdot 10^{-13}$ | $-6\cdot 10^{-7}$ | $1\cdot 10^{-11}$
Table 6: Contribution to $A_{FB}$ at different orders of $\alpha L$ evaluated
at $s_{-}=(87.9\,{\rm GeV})^{2}$, $M_{Z}^{2}$ and $s_{+}=(94.3\,{\rm
GeV})^{2}$ for the cut $z>0.999$. all: full contributions, ang.: only the
angular dependent contributions.

##  Appendix A $N$ space results

In the following we are presenting the first radiators in Mellin $N$-space.

$\displaystyle H_{FB}^{(1),LL}(N)$ $\displaystyle=$ $\displaystyle 0$ (47)
$\displaystyle H_{FB}^{(2),LL}(N)$ $\displaystyle=$
$\displaystyle\frac{8\big{(}3N^{2}+3N-1\big{)}P_{1}}{(N-1)N^{2}(N+1)^{2}(N+2)(2N-1)(2N+3)}$
(48)
$\displaystyle-\frac{32\big{(}4N^{2}+4N-1\big{)}(-1)^{N}}{(2N-1)(2N+1)(2N+3)}[S_{-1}+\ln(2)]$
$\displaystyle H_{FB}^{(3),LL}(N)$ $\displaystyle=$
$\displaystyle\frac{64\big{(}3N^{2}+3N-1\big{)}P_{2}}{(N-1)N^{2}(N+1)^{2}(N+2)(2N-1)(2N+3)}S_{1}$
(49)
$\displaystyle-\frac{4P_{7}}{3(N-1)^{2}N^{3}(N+1)^{3}(N+2)^{2}(2N-3)(2N-1)^{2}(2N+3)^{2}(2N+5)}$
$\displaystyle+(-1)^{N}\Biggl{\\{}\frac{8P_{4}}{N(N+1)(2N-1)(2N+1)(2N+3)}S_{-2}$
$\displaystyle-\frac{256\big{(}4N^{2}+4N-1\big{)}}{(2N-1)(2N+1)(2N+3)}\left[S_{-1,1}-\frac{1}{2}\ln^{2}(2)+\sum_{i=1}^{N}\frac{\ln(2)+S_{-1}\big{(}i\big{)}}{1+2i}\right]$
$\displaystyle-\frac{128P_{5}}{3(2N-3)(2N-1)^{2}(2N+1)^{2}(2N+3)^{2}(5+2N)}S_{-1}$
$\displaystyle+\frac{4P_{3}}{N(N+1)(2N-1)(2N+1)(2N+3)}\zeta_{2}$
$\displaystyle-\frac{128P_{6}}{3(2N-3)(2N-1)^{2}(2N+1)^{2}(2N+3)^{2}(2N+5)}\ln(2)\Biggr{\\}},$

with the polynomials

$\displaystyle P_{1}$ $\displaystyle=$ $\displaystyle
4N^{4}+8N^{3}-N^{2}-5N-3$ (50) $\displaystyle P_{2}$ $\displaystyle=$
$\displaystyle 4N^{4}+8N^{3}-N^{2}-5N-3,$ (51) $\displaystyle P_{3}$
$\displaystyle=$ $\displaystyle 32N^{4}+16N^{3}-48N^{2}-14N+9,$ (52)
$\displaystyle P_{4}$ $\displaystyle=$ $\displaystyle
64N^{4}+80N^{3}-24N^{2}-22N+9,$ (53) $\displaystyle P_{5}$ $\displaystyle=$
$\displaystyle
128N^{7}+256N^{6}-320N^{5}-144N^{4}+1128N^{3}+104N^{2}-936N+27,$ (54)
$\displaystyle P_{6}$ $\displaystyle=$ $\displaystyle
896N^{7}+2944N^{6}-512N^{5}-7344N^{4}-2760N^{3}+2816N^{2}+36N-243,$ (55)
$\displaystyle P_{7}$ $\displaystyle=$ $\displaystyle
5120N^{14}+35840N^{13}+73856N^{12}-22784N^{11}-252848N^{10}-201200N^{9}$ (56)
$\displaystyle+193136N^{8}+264800N^{7}-43255N^{6}-113445N^{5}+7512N^{4}+13587N^{3}$
$\displaystyle+1539N^{2}+8262N-1620.$

Since the expressions for the higher order radiators are voluminous, they are
given in a computer-readable file attached to this paper. Here the functions
$S_{\vec{a}}\equiv S_{\vec{a}}(N)$ denote the harmonic sums [36, 37]

$\displaystyle S_{b,\vec{a}}(N)=\sum_{k=1}^{N}\frac{({\rm
sign}(b))^{k}}{k^{|b|}}S_{\vec{a}}(k),\leavevmode\nobreak\
\leavevmode\nobreak\ S_{\emptyset}=1,\leavevmode\nobreak\ \leavevmode\nobreak\
a,b_{i}\in\mathbb{Z}\backslash\\{0\\}.$ (57)

Despite the fact that already at 2nd order in $z$–space cyclotomic harmonic
polylogarithms occur, cyclotomic harmonic sums are only contributing from 3rd
order onward. As usual, the $N$–space expressions turn out to be structurally
simpler than the $z$–space expressions. However, there are evanescent poles at
half integer arguments $N=1/2,3/2,...$. As has been seen in the small $z$
expansion of the radiators, they are all tractable for $N>1$. The pure singlet
kernels are introducing a rightmost pole at $N=1$ and poles left to this are
allowed.

One may numerically represent the higher order corrections also starting from
Mellin $N$ space, as described in Refs. [45, 38]. For this asymptotic
representations and the recurrences of the respective expressions need to be
known, which can be easily obtained using the package HarmonicSums. The
recurrences are implied by the quantities itself in terms of the hierarchic
sum–structures.

Acknowledgments.
We would like to thank W. Hollik, P. Janot, P. Marquard, G. Passarino, and C.
Schneider for discussions. This project has received funding from the European
Union’s Horizon 2020 research and innovation programme under the Marie
Sklodowska-Curie grant agreement No. 764850, SAGEX, COST action CA16201:
Unraveling new physics at the LHC through the precision frontier.

## References

  * [1] E. Accomando et al. Phys. Rept. 299 (1998) 1–78 [hep-ph/9705442];
J.A. Aguilar-Saavedra et al. hep-ph/0106315;
International Linear Collider Reference Design Report, ILC-REPORT-2007-001,
Eds. J. Brau, Y. Okada, and N. Walker; Vol. 1–4.
G. Aarons et al. [ ILC Collaboration ], International Linear Collider
Reference Design Report Volume 2: Physics At The ILC, [arXiv:0709.1893 [hep-
ph]];
http://www.linearcollider.org/ILC

  * [2] H. Aihara et al. [Linear Collider Collaboration], The International Linear Collider. A Global Project, arXiv:1901.09829 [hep-ex].
  * [3] J. Mnich, The International Linear Collider: Prospects and Possible Timelines, arXiv: 1901.10206 [hep-ex].
  * [4] S. van der Meer, The CLIC Project and the Design for an $e^{+}e^{-}$ Collider, CLIC-NOTE-68, (1988);
R.W. Assmann et al., CLIC Study Team, A 3 TeV $e^{+}e^{-}$ Linear Collider
Based on CLIC Technology, CERN 2000-008;
E. Accomando et al. [CLIC Physics Working Group], Physics at the CLIC multi-
TeV linear collider, arXiv:hep-ph/0412251;
P. Roloff et al. [CLIC and CLICdp Collaborations], The Compact Linear e+e-
Collider (CLIC): Physics Potential, arXiv:1812.07986 [hep-ex].

  * [5] A. Abada et al. [FCC Collaboration], Future Circular Collider : 2 The Lepton Collider (FCC$\\_$ee), CERN-ACC-2018-0057, Eur. Phys. J. Special Topics 228 (2019) no.2, 261–623;
http://tlep.web.cern.ch/

  * [6] J.P. Delahaye et al., Muon Colliders, arXiv:1901.06150 [physics.acc-ph].
  * [7] M.J.G. Veltman, Rev. Mod. Phys. 72 (2000) 341–349
  * [8] J. Blümlein, A. De Freitas, C.G. Raab and K. Schönwald, Phys. Lett. B 791 (2019) 206–209 [arXiv:1901.08018 [hep-ph]].
  * [9] J. Blümlein, A. De Freitas, C.G. Raab and K. Schönwald, Phys. Lett. B 801 (2020) 135196 [arXiv:1910.05759 [hep-ph]].
  * [10] J. Blümlein, A. De Freitas, C. Raab and K. Schönwald, Nucl. Phys. B 956 (2020) 115055 [arXiv:2003.14289 [hep-ph]].
  * [11] J. Ablinger, J. Blümlein, A. De Freitas and K. Schönwald, Nucl. Phys. B 955 (2020) 115045 [arXiv:2004.04287 [hep-ph]].
  * [12] J. Blümlein, A. De Freitas and W.L. van Neerven, Nucl. Phys. B 855 (2012) 508–569 [arXiv:1107.4638 [hep-ph]].
  * [13] F.A. Berends, W.L. van Neerven and G.J.H. Burgers, Nucl. Phys. B 297 (1988) 429–478. Erratum: [Nucl. Phys. B 304 (1988) 921–922].
  * [14] B.A. Kniehl, M. Krawczyk, J.H. Kühn and R.G. Stuart, Phys. Lett. B 209 (1988) 337–342.
  * [15] G. Montagna, O. Nicrosini, F. Piccinini and G. Passarino, Comput. Phys. Commun. 117 (1999) 278–289 [arXiv:hep-ph/9804211 [hep-ph]].
  * [16] D.Y. Bardin, M. Grünewald and G. Passarino, Precision calculation project report, arXiv:hep-ph/9902452 [hep-ph];
D.Y. Bardin and G. Passarino, The standard model in the making: Precision
study of the electroweak interactions, International series of monographs on
physics 104 (Calendron Press, Oxford, 1999).

  * [17] D.Y. Bardin, P. Christova, M. Jack, L. Kalinovskaya, A. Olchevski, S. Riemann and T. Riemann, Comput. Phys. Commun. 133 (2001) 229–395 [arXiv:hep-ph/9908433 [hep-ph]];
A.B. Arbuzov, M. Awramik, M. Czakon, A. Freitas, M.W. Grünewald, K. Mönig, S.
Riemann and T. Riemann, Comput. Phys. Commun. 174 (2006) 728–758 [arXiv:hep-
ph/0507146 [hep-ph]].

  * [18] S. Schael et al. Phys. Rept. 427 (2006) 257–454 [hep-ex/0509008].
  * [19] M. Böhm, W. Hollik, D.Y. Bardin, W. Beenakker, F.A. Berends, M.S. Bilenky, G. Burgers, J.E. Campagne, A. Djouadi, O. Fedorenko, et al., Forward-backward asymmetries, CERN-TH-5536-89.
  * [20] M. Greco and A. F. Grillo, Lett. Nuovo Cim. 15 (1976) 174–178;
M. Greco, G. Pancheri-Srivastava and Y. Srivastava, Nucl. Phys. B 171 (1980)
118–140 [Erratum: Nucl. Phys. B 197 (1982) 543–454];
O. Nicrosini and L. Trentadue, Z. Phys. C 39 (1988) 479–486;
J. E. Campagne and R. Zitoun, Phys. Lett. B 222 (1989) 497–500;
O. Nicrosini and L. trentadue, Structure functions techniques in $e^{+}e^{-}$
collisions in : Radiative corrections for $e^{+}e^{-}$ Collisions, (Springer,
Berlin, 1989), Ed. J.H. Kühn, 25–54;
W. Beenakker, F.A. Berends and S.C. van der Marck, Phys. Lett. B 251 (1990)
299–304;
S. Jadach, B.F.L. Ward and Z. Was, Phys. Lett. B 257 (1991) 213–218;
S. Jadach and S. Yost, Phys. Rev. D 100 (2019) no.1, 013002 [arXiv:1801.08611
[hep-ph]].

  * [21] D.Y. Bardin, M.S. Bilenky, A. Chizhov, A. Sazonov, Y. Sedykh, T. Riemann and M. Sachwitz, Phys. Lett. B 229 (1989) 405–408.
  * [22] D.Y. Bardin, M.S. Bilenky, A. Sazonov, Y. Sedykh, T. Riemann and M. Sachwitz, Phys. Lett. B 255 (1991) 290–296 [arXiv:hep-ph/9801209 [hep-ph]].
  * [23] W. Beenakker, F.A. Berends and W.L. van Neerven, Applications of renormalization group methods to radiative corrections, in : Radiative corrections for $e^{+}e^{-}$ Collisions, (Springer, Berlin, 1989), Ed. J.H. Kühn, 3–24. Print-89-0445 (Leiden).
  * [24] T. Riemann and Z. Was, Mod. Phys. Lett. A 4 (1989) 2487–2491;
Z. Was and S. Jadach, Phys. Rev. D 41 (1990) 1425–1437;
G. Montagna, O. Nicrosini and L. Trentadue, Phys. Lett. B 231 (1989) 492–496

  * [25] B.W. Lynn and R.G. Stuart, Nucl. Phys. B 253 (1985) 216–230;
D.Y. Bardin, M. S. Bilenky, G. Mitselmakher, T. Riemann and M. Sachwitz, Z.
Phys. C 44 (1989) 493–502.;
W.J.P. Beenakker, Electroweak corrections: techniques and applications, PhD
Thesis, (Leiden University, 1989);
D.Y. Bardin, M.S. Bilenky, A. Chizhov, A. Sazonov, O. Fedorenko, T. Riemann
and M. Sachwitz, Nucl. Phys. B 351 (1991) 1–48 [arXiv:hep-ph/9801208 [hep-
ph]];
A. Leike and T. Riemann, Z. Phys. C 51 (1991) 113–118;
W. Beenakker, S.C. van der Marck and W. Hollik, Nucl. Phys. B 365 (1991)
24–78;
W. Beenakker and W. Hollik, Phys. Lett. B 269 (1991) 425–431;
W. Beenakker, A. Denner and A. Kraft, Nucl. Phys. B 410 (1993) 219–244;
W. Beenakker and G. Passarino, Phys. Lett. B 425 (1998) 199–207 [arXiv:hep-
ph/9710376 [hep-ph]];
G. Duckeck and W. Hollik, Electroweak precision tests at LEP, Springer Tracts
Mod. Phys. 162 (2000) 1–161;
J.H. Kühn, A.A. Penin and V.A. Smirnov, Eur. Phys. J. C 17 (2000) 97–105
[arXiv:hep-ph/9912503 [hep-ph]];
B. Feucht, J.H. Kühn, A.A. Penin and V.A. Smirnov, Phys. Rev. Lett. 93 (2004)
101802 [arXiv:hep-ph/0404082 [hep-ph]];
G. Bell, J.H. Kühn and J. Rittinger, Eur. Phys. J. C 70 (2010) 659—671
[arXiv:1004.4117 [hep-ph]];
A.V. Manohar and M. Trott, Phys. Lett. B 711 (2012) 313–316 [arXiv:1201.3926
[hep-ph]]; I. Dubovyk, A. Freitas, J. Gluza, T. Riemann and J. Usovitsch, JHEP
08 (2019) 113 [arXiv:1906.08815 [hep-ph]].

  * [26] D.Y. Bardin, M.S. Bilenky, A. Chizhov, A. Sazonov, O. Fedorenko, T. Riemann and M. Sachwitz, Nucl. Phys. B 351 (1991) 1–48 [arXiv:hep-ph/9801208 [hep-ph]].
  * [27] J. Jersak, E. Laermann and P.M. Zerwas, Phys. Lett. B 98 (1981) 363–366; Phys. Rev. D 25 (1982) 1218–1228 [Erratum: Phys. Rev. D 36 (1987) 310];
A. Djouadi, Z. Phys. C 39 (1988) 561–568;
A.B. Arbuzov, D.Y. Bardin and A. Leike, Mod. Phys. Lett. A 7 (1992) 2029–2038
[Erratum: Mod. Phys. Lett. A 9 (1994) 1515];
J.B. Stav and H.A. Olsen, Phys. Rev. D 52 (1995) 1359–1368;
G. Altarelli and B. Lampe, Nucl. Phys. B 391 (1993) 3–22;
P. Nason and B.R. Webber, Nucl. Phys. B 421 (1994) 473–517 [Erratum: Nucl.
Phys. B 480 (1996) 755];
P.J. Rijken and W.L. van Neerven, Phys. Lett. B 386 (1996) 422–428 [arXiv:hep-
ph/9604436 [hep-ph]]; Nucl. Phys. B 487 (1997) 233-282 [arXiv:hep-ph/9609377
[hep-ph]]; Phys. Lett. B 392 (1997) 207–215 [arXiv:hep-ph/9609379 [hep-ph]];
V. Ravindran and W.L. van Neerven, Phys. Lett. B 445 (1998) 214–222
[arXiv:hep-ph/9809411 [hep-ph]];
W. Bernreuther, R. Bonciani, T. Gehrmann, R. Heinesch, T. Leineweber, P.
Mastrolia and E. Remiddi, Nucl. Phys. B 750 (2006) 83–107; [arXiv:hep-
ph/0604031 [hep-ph]];
J.A. Aguilar-Saavedra, W. Bernreuther and Z.G. Si, Phys. Rev. D 86 (2012)
115020 [arXiv:1209.6352 [hep-ph]];
J. Gao and H.X. Zhu, Phys. Rev. Lett. 113 (2014) no.26, 262001
[arXiv:1410.3165 [hep-ph]];
W. Bernreuther, L. Chen, O. Dekkers, T. Gehrmann and D. Heisler, JHEP 01
(2017) 053 [arXiv:1611.07942 [hep-ph]];
D. d’Enterria and C. Yan, Revised QCD effects on the Z $\to b\bar{b}$ forward-
backward asymmetry, [arXiv:2011.00530 [hep-ph]];
A.B. Arbuzov, T.V. Kopylova and G.A. Seilkhanova, Mod. Phys. Lett. A 35 (2020)
no.25, 2050210 [arXiv:2003.14054 [hep-ph]].

  * [28] M. Skrzypek, Acta Phys. Polon. B23 (1992) 135–172;
M. Jezabek, Z. Phys. C56 (1992) 285–288;
M. Przybycien, Acta Phys. Polon. B 24 (1993) 1105–1114 [hep-th/9511029];
J. Blümlein, S. Riemersma and A. Vogt, Eur. Phys. J. C 1 (1998) 255–259 [hep-
ph/9611214];
A.B. Arbuzov, Phys. Lett. B 470 (1999) 252–258 [hep-ph/9908361];
A.B. Arbuzov, JHEP 0107 (2001) 043 [hep-ph/9907500];
J. Blümlein and H. Kawamura, Nucl. Phys. B 708 (2005) 467–510 [hep-
ph/0409289]; Eur. Phys. J. C 51 (2007) 317–333 [arXiv:hep-ph/0701019].

  * [29] R. Hamberg, W.L. van Neerven and T. Matsuura, Nucl. Phys. B 359 (1991) 343–405 Erratum: [Nucl. Phys. B 644 (2002) 403–404].
  * [30] R.V. Harlander and W.B. Kilgore, Phys. Rev. Lett. 88 (2002) 201801 [hep-ph/0201206].
  * [31] E. Remiddi and J.A.M. Vermaseren, Int. J. Mod. Phys. A15 (2000) 725–754 [hep-ph/9905237].
  * [32] J. Ablinger, J. Blümlein and C. Schneider, J. Math. Phys. 52 (2011) 102301 [arXiv:1105.6063 [math-ph]].
  * [33] C. Schneider, Sém. Lothar. Combin. 56 (2007) 1–36 article B56b.
  * [34] C. Schneider, Simplifying Multiple Sums in Difference Fields, in: Computer Algebra in Quantum Field Theory: Integration, Summation and Special Functions Texts and Monographs in Symbolic Computation eds. C. Schneider and J. Blümlein (Springer, Wien, 2013) 325–360 [arXiv:1304.4134 [cs.SC]].
  * [35] J. Ablinger, J. Blümlein and C. Schneider, J. Math. Phys. 54 (2013) 082301 [arXiv:1302.0378 [math-ph]];
J. Ablinger, J. Blümlein, C.G. Raab and C. Schneider, J. Math. Phys. 55 (2014)
112301 [arXiv:1407.1822 [hep-th]];
J. Blümlein, Comput. Phys. Commun. 159 (2004) 19–54 [arXiv:hep-ph/0311046
[hep-ph]];
J. Ablinger, J. Blümlein and C. Schneider, J. Phys. Conf. Ser. 523 (2014)
012060 [arXiv:1310.5645 [math-ph]];
J. Ablinger, PoS (LL2014) 019 [arXiv:1407.6180[cs.SC]]; A Computer Algebra
Toolbox for Harmonic Sums Related to Particle Physics, Diploma Thesis, JKU
Linz, 2009, arXiv:1011.1176[math-ph]; Computer Algebra Algorithms for Special
Functions in Particle Physics, Ph.D. Thesis, Linz U. (2012)
arXiv:1305.0687[math-ph]; PoS (LL2016) 067; Experimental Mathematics 26 (2017)
[arXiv:1507.01703 [math.CO]]; PoS (RADCOR2017) 001 [arXiv:1801.01039 [cs.SC]];
arXiv:1902.11001 [math.CO]; PoS (LL2018) 063.

  * [36] J. Vermaseren, Int. J. Mod. Phys. A14 (1999) 2037–2076 [hep-ph/9806280].
  * [37] J. Blümlein and S. Kurth, Phys. Rev. D60 (1999) 014018 [hep-ph/9810241].
  * [38] J. Blümlein, Comput. Phys. Commun. 180 (2009) 2218–2249 [arXiv:0901.3106 [hep-ph]].
  * [39] M. Böhm, A. Denner, and H. Joos, Gauge Theories of the Strong and Electroweak Interaction, (B.G. Teubner, Stuttgart, 2001);
A. Denner and S. Dittmaier, Electroweak Radiative Corrections for Collider
Physics, arXiv:1912.06823 [hep-ph].

  * [40] P. Janot, JHEP 02 (2016) 053 [Erratum: JHEP 11 (2017) 164] [arXiv:1512.05544 [hep-ph]].
  * [41] P.A. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. (2020) 083C01.
  * [42] T. Gehrmann and E. Remiddi, Comput. Phys. Commun. 141 (2001) 296–312 [hep-ph/0107173].
  * [43] J. Ablinger, J. Blümlein, M. Round and C. Schneider, Comput. Phys. Commun. 240 (2019) 189–201 [arXiv:1809.07084 [hep-ph]].
  * [44] J. Blümlein and H. Kawamura, Phys. Lett. B 553 (2003) 242–250 [hep-ph/0211191];
A.B. Arbuzov, Phys. Part. Nucl. 50 (2019) no.6, 721–825.

  * [45] J. Blümlein, Comput. Phys. Commun. 133 (2000) 76–104 [hep-ph/0003100];
J. Blümlein and S.O. Moch, Phys. Lett. B 614 (2005) 53–61 [hep-ph/0503188].


# Methods for the Numerical Analysis of Boundary Value Problem of Partial
Differential Equations Based on Kolmogorov Superposition Theorem

Korney Tomashchuk

###### Abstract

This research introduces a new method for the transition from partial to
ordinary differential equations that is based on the Kolmogorov superposition
theorem. In this paper, we discuss the numerical implementation of the
Kolmogorov theorem and propose an approach that allows us to apply the theorem
to represent partial derivatives of multivariate function as a combination of
ordinary derivatives of univariate functions. We tested the method by running
a numerical experiment with the Poisson equation. As a result, we managed to
get a system of ordinary differential equations whose solution coincides with
a solution of the initial partial differential equation.

Keywords: Kolmogorov superposition theorem, partial differential equations,
boundary value problem of ordinary differential equations.

## 1 Introduction

Partial differential equations are a powerful instrument of mathematics that
allows to model various processes of different complexities. PDEs111Partial
differential equations are applied to many real-life problems. For example,
this family of equations may be used for modeling blood flow in the human
venous system or for price evaluation of some financial instrument. However,
very often partial differential equations are very hard to solve. Some of them
we even can not solve at all, so we use modern methods of machine learning to
approximate the solution. Obviously, bifurcation analysis of partial
differential equations is even more complex task to solve and mathematicians
have not managed to develop a general procedure for how to do it.

Ordinary differential equation is a much simpler class of equations. The
number of processes that can be described by them is limited, but we can solve
any type of them using methods of numerical integration and analyse behaviour
of a dynamic system generated by them using bifurcation analysis.

Thus in this research, we worked on the derivation of the method for the
transition from partial differential equations to equivalent ordinary
differential equations, as this will allow us to apply all the existing
methods of solution and analysis of ODEs222Ordinary differential equation to
PDEs.

The main research aim of our work is to develop a method, based on the
Kolmogorov superposition theorem, for the transition from partial to ordinary
differential equations, such that we will be able not only to solve PDE’s
using method of simple numerical integration, but we will be able to apply
existing methods of bifurcation analysis to partial differential equations.

## 2 Problem Statement

The problem can be formulated as follows. Consider a boundary problem of
partial differential equations that is given by (2.1) and corresponding
boundary conditions (2.2)

$F\Big{(}x_{1},\dots,x_{n},f(x_{1},\dots,x_{n}),\frac{\partial
f(x_{1},\dots,x_{n})}{\partial
x_{i_{1}}},\frac{\partial^{2}f(x_{1},\dots,x_{n})}{\partial x_{i_{1}}\partial
x_{i_{2}}},\dots,\frac{\partial^{N}f(x_{1},\dots,x_{n})}{\partial
x_{i_{1}}\dots\partial x_{i_{N}}}\Big{)}=0$ (2.1)

$\begin{split}\varphi_{k_{1}}\Big{(}x_{1},\dots,x_{n},f(x_{1},\dots,x_{n}),\frac{\partial
f(x_{1},\dots,x_{n})}{\partial
x_{i_{1}}},&\frac{\partial^{2}f(x_{1},\dots,x_{n})}{\partial x_{i_{1}}\partial
x_{i_{2}}},\\\ &\dots,\frac{\partial^{N-1}f(x_{1},\dots,x_{n})}{\partial
x_{i_{1}}\dots\partial
x_{i_{N}}}\Big{)}\Bigg{|}_{x_{t_{1}}=x_{t_{1}}^{min}}=0\\\
\varphi_{k_{2}}\Big{(}x_{1},\dots,x_{n},f(x_{1},\dots,x_{n}),\frac{\partial
f(x_{1},\dots,x_{n})}{\partial
x_{i_{1}}},&\frac{\partial^{2}f(x_{1},\dots,x_{n})}{\partial x_{i_{1}}\partial
x_{i_{2}}},\\\ &\dots,\frac{\partial^{N-1}f(x_{1},\dots,x_{n})}{\partial
x_{i_{1}}\dots\partial
x_{i_{N}}}\Big{)}\Bigg{|}_{x_{t_{2}}=x_{t_{2}}^{max}}=0\\\ \end{split}$ (2.2)

where
$k_{1},k_{2}\in\mathbb{N},t_{1},t_{2}\in\\{1,\dots,n\\},i_{j}\in\\{1,\dots,n\\}$,
$j=1,\dots,N$. We seek to find ordinary differential equation (2.3)

$\tilde{F}\Big{(}x,y(x),\frac{dy(x)}{dx},\frac{d^{2}y(x)}{dx^{2}},\dots,\frac{d^{\tilde{N}}y(x)}{dx^{\tilde{N}}}\Big{)}=0,$
(2.3)

such that solution of the equation coincides with the solution of initial
problem (2.1, 2.2). To solve this kind of problems we developed a method based
on KST333Kolmogorov superposition theorem that we present in the sections 4-7.

## 3 Related Works

### 3.1 Method for transition from partial to ordinary differential equations

As partial differential equations have ben used for mathematical modeling in
various fields for a long time, the idea of reducing them to ordinary
differential equations have been studied by many mathematicians. In a review
paper [2] group of researchers give a detailed review of existing methods for
reduction of PDEs to ODEs. They carry various experiments in order to compare
and outline pros and cons of each method. Two other mathematicians in their
work [6] write about superiority of the Kontorovich method on other existing
methods supporting their opinions with several numerical experiments. However,
existing methods have their own disadvantages. Some of them use approximations
that lead to uncertainties. Other methods that use orthogonal basis, such as
Fourier based methods, break bifurcation picture of the equations making
conclusions of bifurcation theory irrelevant. Our method is based on a strict
equality that is known as Kolmorogov superposition theorem.

### 3.2 Kolmogorov superposition theorem

The method that we propose in this paper relies heavily on Sprecher’s
formulation of the Kolmogorov superposition theory that he introduced in the
article [7]. This formulation is more convenient in terms of numerical
realisation. In the same article he proposes algorithms for calculating inner
and outer functions that are required in order to apply the theorem on
practice.

Despite that Sprecher’s formulation of the theorem was correct, the algorithm
for calculating internal functions was wrong. For certain values of parameters
the monotonicity property of these functions is violated. Another
Mathematician Mario Koppen who also dedicated part of his work to study the
Kolmogorov theorem, spotted this mistake and corrected. Proper algorithm for
constructing inner functions was introduced in the paper [5].

Jonas Actor in his PhD thesis [1] on the study of the properties of the
Kolmogorov theorem, give a review of existing formulations of Kolmogorov
theorem, and proposes a new type of improved inner functions that are
Lipschitz continuous. Initial functions that were proposed by David Sprecher
have very steep slope on some segments. This can lead to computational
problems. In contrast, Lipschitz continuous functions have more controlled
slope.

T. Hedberb also studied properties of the Komogorov superposition theorem. In
the work [3] he derived his own formulation of the theorem that does not use
strictly defined inner functions. Instead, he have showed his reformulation is
correct for quasi-all functions if only they satisfy certain conditions
described in [3].

## 4 Numerical implementation of Kolmogorov superposition theorem

### 4.1 Theorem formulation review

The method for transition from partial to ordinary differential equations is
based on Sprecher’s reformulation of Kolmogorov theorem. We start with
providing a quick overview of the theorem. Consider a function of $n$
variables $f(x_{1},\dots,x_{n})$. As stated in the paper [7],
$f(x_{1},\dots,x_{n})$ can be written as a sum of functions of one variable
(4.1)

$f(x_{1},x_{2},\dots,x_{n})=\sum_{q=0}^{2n}\Phi_{q}\Big{[}\sum_{p=1}^{n}\alpha_{p}\psi(x_{p}+aq)\Big{]},$
(4.1)

where constants $a,\alpha_{p},p=1,\dots,n$ and inner function $\psi(x)$ are
the same for different $f(x_{1},\dots,x_{n})$. Moreover, function $\psi(x)$ is
a specific function that can be calculated numerically using algorithm of
Mario Koppen described in [5]. Initial algorithm proposed by David A. Sprecher
contained a mistake for some combination of algorithm parameters. Further, we
provide a quick review of a correct algorithm.

### 4.2 Numerical implementation of inner $\psi(x)$ function

Let’s introduce notation $D$ for the set of terminating rational numbers that
are defined according to equation (4.2).

$d_{k}=\sum_{j=1}^{k}\frac{i_{j}}{\gamma^{k}},i_{j}=0,\dots,\gamma-1,$ (4.2)

where $k,\gamma$ are parameters of the algorithm such that $\gamma\geq 2n+2$
and $k\in\mathbb{N}$. Next, we define constants $a,\alpha_{p}$ according to
the formulas (4.3).

$\begin{split}a&=\frac{1}{\gamma(\gamma-1)}\\\
\alpha_{p}&=\sum_{j=1}^{\infty}\gamma^{-(p-1)\frac{n^{p}-1}{n-1}}\end{split}$
(4.3)

Then, the inner function $\psi(x)$ is recursively defined by the formula (4.4)

$\psi^{k}(d_{k})=\begin{cases}d_{k},\text{ for }k=1\\\
\psi^{k-1}(d_{k}-\frac{i_{k}}{\gamma_{k}})+\frac{i_{k}}{\gamma^{(n^{k}-1)/(n-1)}}\text{
for }k>1,i_{k}<\gamma-1\\\
\frac{1}{2}\Big{(}\psi_{k}(d_{k}-\frac{1}{\gamma^{k}})+\psi_{k-1}(d_{k}+\frac{1}{\gamma_{k}})\Big{)}\text{
for }k>1,i_{k}=\gamma-1\end{cases}$ (4.4)

where $i_{k}$ are decimal integers of rational number $d_{k}=0.i_{1}i_{2}\dots
i_{k}$. In most of the problems the function $f(x_{1},\dots,x_{n})$ is defined
beyond the unit cube. For $x>1$, function $\psi(x)$ is calculated using
formula (4.5).

$\psi(x)=\psi(x-[x])+[x],$ (4.5)

where $[x]$ is the integer part of $x$.

According to the work [5] of Mario Koppen $\gamma=10$ is the best value of
parameter for $n=2$. Graphs for function $\psi(x)$ for $k=2,3,4$ and
$\gamma=10$ are depicted in the Figure 1. Function $\psi(x)$ is defined only
at certain points $d_{k}$ but for the better representation we use linear
interpolation.

Figure 1: $\psi(x)$ for $k=2,3,4$ and $\gamma=10$

As we are going to apply Kolmogorov theorem to the PDE problems, we should be
able to define derivatives of inner functions. $\psi(x)$ cannot be
analytically differentiated. Thus, we consider difference analogue of the
derivative (4.6).

$\psi^{\prime}(x)=\frac{\psi(x+\Delta)-\psi(x)}{\Delta},$ (4.6)

where $\Delta$ is an increment that depends on the choice of the parameters
$\gamma$ and $k$. As $\psi(x)$ is defined only at a finite number of points
$d_{k}$, we take $\Delta=\gamma^{-k}$ equal to the size of the step between
two neighbour points from the set $D_{k}$. Graphs of the first two derivates
of the function $\psi(x)$ with parameters $k=4,\gamma=10$ are depicted in the
Figure 2.

Figure 2: Graphs of $\psi(x),\psi^{\prime}(x),\psi^{\prime\prime}(x)$ with
parameters $k=4,\gamma=10$.

## 5 Representation of a multivariate function as a sum of functions of one
variable

Transition from partial to ordinary differential equations consists of two
main parts. First, we should represent function of several variable as a
combination of functions of one variable. Second, substitute partial
derivatives with corresponding ordinary. To complete the first stage, we use
Kolmogorov theorem and the fact that $a=\frac{1}{\gamma(\gamma-1)}$ is a small
parameter. For recommended value of $\gamma=10$ for $n=2$, $a$ is equal to
$\frac{1}{90}$. As there are no upper bound for a value of $\gamma$, we can
consider taking $\gamma\rightarrow\infty$, as a result $a\rightarrow 0$. Thus,
we can consider taking Taylor series of (4.1) in the vicinity $a=0$ (see
equation (5.1)).

$\begin{split}\sum_{q=0}^{2n}\Phi_{q}\Big{(}\sum_{p=0}^{n}&\alpha_{p}\psi(x_{p}+qa)\Big{)}=\sum_{m=0}^{M}\frac{a^{m}}{m!}\frac{d^{m}}{da^{m}}\Bigg{(}\sum_{q=0}^{2n}\Phi_{q}\Big{(}\sum_{i=p}^{n}\alpha_{p}\psi(x_{p}+aq)\Big{)}\Bigg{)}\Bigg{|}_{a=0}\end{split}$
(5.1)

In order to get representation of Taylor series with $M$ terms of KST theorem,
we use Bell’s polynomials and Faà di Bruno’s formula. Definition of both are
provided below.

### 5.1 Bell’s polynomial

By definition, Bell’s polynomials are defined by the formula

$B_{n,k}(x_{1},\dots,x_{n-k+1})=\sum\bigg{(}\frac{n!}{j_{1}!j_{2}!\dots
j_{n-k+1}!}\prod_{i=1}^{n-k+1}\Big{(}\frac{x_{i}}{i!}\Big{)}^{j_{i}}\bigg{)},$
(5.2)

where outer sum is taken over all non-negative numbers
$j_{1},j_{2},\dots,j_{n-k+1}$ that satisfy two conditions (5.3, 5.4).

$\sum_{i=1}^{n-k+1}j_{i}=k$ (5.3) $\sum_{i=1}^{n-k+1}ij_{i}=n$ (5.4)

### 5.2 Faà di Bruno’s formula

Faà di Bruno’s formula is given by the equation (5.5)

$\frac{d^{m}}{dx^{m}}f(g(x))=\sum\bigg{(}\frac{m!}{j_{1}!j_{2}!\dots
j_{m}!}\cdot
f^{(j_{1}+\dots+j_{m})}\big{(}g(x)\big{)}\prod_{i=1}^{m}\Big{(}\frac{g^{(i)}(x)}{i!}\Big{)}^{j_{i}}\bigg{)},$
(5.5)

where external sum is taken over all nonnegative integers that satisfy second
condition of Bell’s polynomial (5.4). It is used to define high-order
derivatives of complex functions. Using Bell’s polynomials, we can rewrite
formula (5.5) in the following form.

$\frac{d^{m}}{dx^{m}}f(g(x))=\sum_{k=0}^{m}f^{(k)}\big{(}g(x)\big{)}B_{m,k}\big{(}g^{\prime}(x),g^{\prime\prime}(x),\dots,g^{(m-k-1)}(x)\big{)}$
(5.6)

### 5.3 Taylor series of Kolmogorov superposition theorem

We can substitute formula (5.6) to achieve compact representation of Taylor
series (5.1) as shown in the equation (5.7)

$\begin{split}&\frac{d^{m}}{da^{m}}\Bigg{(}\sum_{q=0}^{2n}\Phi_{q}\Big{(}\sum_{p=1}^{n}\alpha_{p}\psi(x_{p}+aq)\Big{)}\Bigg{)}\Bigg{|}_{a=0}=\sum_{k=0}^{m}\sum_{q=0}^{2n}\Phi_{q}^{(k)}\Big{(}\sum_{p=1}^{n}\alpha_{p}\psi(x_{p})\Big{)}\cdot\\\
&\cdot
B_{m,k}\big{(}q\sum_{p=1}^{n}\alpha_{p}\psi^{\prime}(x_{p}),q^{2}\sum_{p=1}^{n}\alpha_{p}\psi^{\prime\prime}(x_{p}),\dots,q^{(m-k-1)}\sum_{p=1}^{n}\alpha_{p}\psi^{(m-k-1)}(x_{p})\big{)}\end{split}$
(5.7)

Using definition of Bell’s polynomial and applying sequence of simple
transformations, we get more convinient representation (5.8).

$\begin{split}&B_{m,k}\big{(}q\sum_{p=1}^{n}\alpha_{p}\psi^{\prime}(x_{p}),q^{2}\sum_{p=1}^{n}\alpha_{p}\psi^{\prime\prime}(x_{p}),\dots,q^{(m-k-1)}\sum_{p=1}^{n}\alpha_{p}\psi^{(m-k-1)}(x_{p})\big{)}=\\\
&=\sum\bigg{(}\frac{m!}{j_{1}!j_{2}!\dots
j_{m-k+1}!}\prod_{i=1}^{m-k+1}\Big{(}\frac{q^{i}\sum_{p=1}^{n}\alpha_{p}\psi^{(i)}(x_{p})}{i!}\Big{)}^{j_{i}}\bigg{)}=q^{m}\tilde{B}_{m,k}(x_{1},x_{2}\dots,x_{n})\end{split}$
(5.8)

where

$\tilde{B}_{m,k}(x_{1},x_{2}\dots,x_{n})=\sum\bigg{(}\frac{m!}{j_{1}!j_{2}!\dots
j_{m-k+1}!}\prod_{i=1}^{m-k+1}\Big{(}\frac{\sum_{p=1}^{n}\alpha_{p}\psi^{(i)}(x_{p})}{i!}\Big{)}^{j_{i}}\bigg{)}$
(5.9)

Then, Taylor series of KST representation of function of $n$ variables with
arbitrary number of terms ($M$) can be represented with the equation (5.10)

$\begin{split}\sum_{q=0}^{2n}\Phi_{q}\Big{(}\sum_{p=0}^{n}\alpha_{p}\psi(x_{p}+qa)\Big{)}=\sum_{m=0}^{M}\sum_{k=0}^{m}\tilde{B}_{m,k}(x_{1},\dots,x_{n})\sum_{q=0}^{2n}\frac{a^{m}q^{m}}{m!}\Phi_{q}^{(k)}\Big{(}\sum_{p=1}^{n}\alpha_{p}\psi(x_{p})\Big{)}\end{split}$
(5.10)

Finally, by making a substitution $z=\sum_{p=1}^{n}\alpha_{p}\psi(x_{p})$, we
derive representation of multivariate $f(x_{1},x_{2},\dots,x_{n})$ as a sum of
univariate functions $\Phi_{q}(z)$ (see equation (5.11)).

$f(x_{1},x_{2},\dots,x_{n})=\sum_{m=0}^{M}\sum_{k=0}^{m}\tilde{B}_{m,k}(x_{1},x_{2}\dots,x_{n})\sum_{q=0}^{2n}\frac{a^{m}q^{m}}{m!}\Phi_{q}^{(k)}(z),$
(5.11)

## 6 Method for transition from partial to ordinary differential equations

Method for transition from partial to ordinary differential equations that we
derived can be applied only to variational formulation of PDE problem. Let’s
consider functional $\Xi$

$\Xi=\int_{x}L\Big{(}x_{1},\dots,x_{n},f(x_{1},\dots,x_{n}),\frac{\partial
f(x_{1},\dots,x_{n})}{\partial
x_{i_{1}}},\frac{\partial^{2}f(x_{1},\dots,x_{n})}{\partial x_{i_{1}}\partial
x_{i_{2}}},\dots,\frac{\partial^{N}f(x_{1},\dots,x_{n})}{\partial
x_{i_{1}}\dots\partial x_{i_{N}}}\Big{)}dx$ (6.1)

such that solution of variational problem (6.2) coincide with the solution of
initial PDE problem (2.1, 2.2).

$\Xi\longrightarrow\underset{x}{extr}$ (6.2)

In order to get rid of partial derivatives in functional, we use obtained
formula (5.11). First, we change variable of integrations. We parameterise
$x_{2},\dots,x_{n}$ and $x_{1}$ express through other parametrised variables
as shown in (6.3).

$\begin{split}x_{2}&=\tilde{x}_{2}\\\ \dots&\dots\dots{\\\
}x_{n}&=\tilde{x}_{n}\\\
x_{1}&=\psi^{-1}\Big{(}\frac{z-\sum_{p=2}^{n}\alpha_{p}\psi(x_{p})}{\alpha_{1}}\Big{)}\end{split}$
(6.3)

Boundaries of integration change according to formulas (6.4)

$\begin{split}&x_{1}^{min}\mapsto
z^{min}=\alpha_{1}\psi(x_{1}^{min})+\alpha_{2}\psi(\tilde{x}_{2})+\dots+\alpha_{n}\psi(\tilde{x}_{n})\\\
&x_{1}^{max}\mapsto
z^{max}=\alpha_{1}\psi(x_{1}^{max})+\alpha_{2}\psi(\tilde{x}_{2})+\dots+\alpha_{n}\psi(\tilde{x}_{n})\\\
\end{split}$ (6.4)

and integration factor according to determinant of Jacobian matrix (6.5).
Determinant of a diagonal matrix equal to product of diagonal elements.

$\begin{split}&J=\begin{bmatrix}z^{\prime}_{x_{1}}&z^{\prime}_{x_{2}}&\dots&z^{\prime}_{x_{n}}\\\
{x^{\prime}_{2}}_{x_{1}}&{x^{\prime}_{2}}_{x_{2}}&\dots&{x^{\prime}_{2}}_{x_{n}}\\\
\dots&\dots&\dots&\dots\\\
{x^{\prime}_{n}}_{x_{1}}&{x^{\prime}_{n}}_{x_{2}}&\dots&{x^{\prime}_{n}}_{x_{n}}\end{bmatrix}=\begin{bmatrix}\alpha_{1}\psi^{\prime}(x_{1})&\alpha_{2}\psi^{\prime}(x_{2})&\dots&\alpha_{n}\psi^{\prime}(x_{n})\\\
0&1&\dots&0\\\ \dots&\dots&\dots&\dots\\\ 0&0&\dots&1\end{bmatrix}\\\
&\implies
dx_{1}=\frac{dz}{\alpha_{1}\psi^{\prime}\circ\psi^{-1}\Big{(}\frac{z-\sum_{p=2}^{n}\alpha_{p}\psi(\tilde{x}_{p})}{\alpha_{1}}\Big{)}}\end{split}$
(6.5)

Secondly, we redefine partial derivatives of $f(x_{1},\dots,x_{n})$ as shown
in (6.6)

$\begin{split}&\frac{\partial^{N}f(x_{1},\dots,x_{n})}{\partial
x_{i_{1}}\dots\partial
x_{i_{N}}}=\sum_{m=0}^{M}\sum_{k=0}^{m}\Bigg{[}\sum_{j=0}^{N}\bigg{[}\frac{\partial^{N}}{\partial
x_{i_{1}}\dots\partial
x_{i_{N}}}\Big{(}\tilde{B}_{m,k}(x)\sum_{q=0}^{2n}\frac{a^{m}q^{m}}{m!}\Phi_{q}^{(k)}(z)\Big{)}\bigg{]}\Bigg{]}\end{split}$
(6.6)

where $\tilde{B}_{m,k}$ and $z$ are known from KST theorem. Thus, in terms of
unknown functions $\Phi_{q}(z)$ we get expression that contains only ordinary
derivatives. Applying described procedure to the functional (6.1), we derive
equivalent functional (6.7) that does not contain partial derivatives.

$\begin{split}\tilde{\Xi}=\int_{\tilde{x}}\int_{z^{min}}^{z^{max}}L\Big{(}z,\tilde{x}_{2},\dots,\tilde{x}_{n},\Phi_{0}(z),\dots,\Phi_{2n}(z),&\\\
\frac{d\Phi_{0}(z)}{dz},\dots,\frac{d\Phi_{2n}(z)}{dz}&,\dots,\frac{d^{N}\Phi_{0}(z)}{dz^{N}},\dots,\frac{d^{N}\Phi_{2n}(z)}{dz^{N}}\Big{)}dzd\tilde{x}\end{split}$
(6.7)

Then, initial PDE problem is reduced to a much simpler variational problem
with ordinary derivatives (6.8). This problem can be solved using standard
methods of calculus of variations.

$\begin{split}\forall\\{\tilde{x}_{2},\dots,\tilde{x}_{n}\\},\>\tilde{x}_{2}\in[\tilde{x}_{2}^{min},&\tilde{x}_{2}^{max}],\dots,\tilde{x}_{n}\in[\tilde{x}_{n}^{min},\tilde{x}_{n}^{max}]:\\\
&\tilde{\Xi}\rightarrow\underset{z}{extr}\end{split}$ (6.8)

## 7 Numerical experiment with Poisson equation

### 7.1 Equivalent variational formulation

In order to check our method, we decided to apply it to the partial
differential equation with known solution which is the Poisson equation. Let’s
consider the following formulation (7.1) with boundary conditions (7.2).

$\nabla^{2}u(x_{1},x_{2})=\sin(\pi x_{1})\sin(\pi x_{2})$ (7.1)
$u(0,x_{2})=0,u(1,x_{2})=0,u(x_{1},0)=0,u(x_{1},1)=0$ (7.2)

Solution for problem (7.1, 7.2) is given by equation (7.3).

$u^{*}(x_{1},x_{2})=\frac{\sin(\pi x_{1})\sin(\pi x_{2})}{-2\pi^{2}}$ (7.3)

One can easily check that (7.3) satisfies both equation and boundary
conditions. As shown in Section 6, we should first find equivalent variational
formulation. For PDE problem (7.1, 7.2) corresponding problem of calculus of
variations is given by (7.4). The whole derivation procedure is shown in the
appendix A.

$\begin{split}\int_{0}^{1}\int_{0}^{1}\bigg{[}&-\Big{(}\frac{\partial{u(x_{1},x_{2})}}{\partial{x_{1}}}\Big{)}^{2}-\Big{(}\frac{\partial{u(x_{1},x_{2})}}{\partial{x_{2}}}\Big{)}^{2}-2\sin(\pi
x_{1})\sin(\pi x_{2})u(x_{1},x_{2})-\\\
&-2\frac{\partial^{2}u(x_{1},x_{2})}{\partial
x_{1}^{2}}u(x_{1},x_{2})-2\frac{\partial^{2}u(x_{1},x_{2})}{\partial
x_{2}^{2}}u(x_{1},x_{2})\bigg{]}dx_{1}dx_{2}\longrightarrow extr\end{split}$
(7.4)

### 7.2 ODE boundary problem of Poisson equation

In order to derive corresponding system of ordinary differential equations, we
apply the procedure that we described in section 6. Let’s consider the
simplest case when we remove all terms of Taylor series except the first one.
Then, according to (5.11) for $n=2,M=0$ function $u(x_{1},x_{2})$ can be
written as (7.5)

$u(x_{1},x_{2})=\sum_{q=0}^{4}\Phi_{q}(z),$ (7.5)

where $z=\alpha_{1}\psi(x_{1})+\alpha_{2}\psi(x_{2})$. Further, following
algorithm described above we apply substitutions (7.6)

$\begin{split}&x_{2}=\tilde{x}_{2}\\\
&x_{1}=\psi^{-1}\Big{(}\frac{z-\alpha_{2}\psi(\tilde{x}_{2})}{\alpha_{1}}\Big{)}\end{split}$
(7.6)

and change factor and boundaries of integration (7.7, 7.8).

$\begin{split}x_{1}&=0\rightarrow z^{min}=\alpha_{2}\psi(\tilde{x}_{2})\\\
x_{1}&=1\rightarrow z^{max}=\alpha_{1}+\alpha_{2}\psi(\tilde{x}_{2})\\\
\end{split}$ (7.7)
$\begin{split}J&=\begin{bmatrix}z^{\prime}_{x_{1}}&z^{\prime}_{x_{2}}\\\
{x_{2}}^{\prime}_{x_{1}}&{x_{2}}^{\prime}_{x_{2}}\end{bmatrix}=\begin{bmatrix}\alpha_{1}\psi^{\prime}(x_{1})&\alpha_{2}\psi^{\prime}(x_{2})\\\
0&1\end{bmatrix}\\\ &\implies
dx_{1}=\frac{dz}{\alpha_{1}\psi^{\prime}(x_{1})}=\frac{dz}{\alpha_{1}\psi^{\prime}\psi^{-1}\big{(}\frac{z-\alpha_{2}\psi(\tilde{x}_{2})}{\alpha_{1}}\big{)}}\end{split}$
(7.8)

Next step is to define partial derivatives up to second through ordinary
derivatives of functions $\Phi_{q}(z)$ (see equation (7.9)).

$\begin{split}&\frac{\partial u(x_{1},x_{2})}{\partial
x_{i}}=\alpha_{i}\psi^{\prime}(x_{i})\sum_{q=0}^{4}\Phi_{q}^{\prime}(z)\\\
&\frac{\partial^{2}u(x_{1},x_{2})}{\partial
x_{i}^{2}}=\alpha_{i}\psi^{\prime\prime}(x_{i})\sum_{q=0}^{4}\Phi_{q}^{\prime}(z)+\alpha_{i}^{2}\psi^{\prime
2}(x_{i})\sum_{q=0}^{4}\Phi_{q}^{\prime\prime}(z)\\\ \end{split}$ (7.9)

As a result, using all above substitutions, we reformulate PDE Poisson problem
as a problem of finding extremum of a functional that contains only ordinary
derivatives (7.10)

$\begin{split}\int_{0}^{1}&\int_{z_{min}}^{z_{max}}\frac{1}{\alpha_{1}\psi^{\prime}\psi^{-1}\big{(}\frac{z-\alpha_{2}\psi(\tilde{x}_{2})}{\alpha_{1}}\big{)}}\bigg{[}-\Big{(}\alpha_{1}\psi^{\prime}\psi^{-1}\big{(}\frac{z-\alpha_{2}\psi(\tilde{x}_{2})}{\alpha_{1}}\big{)}\sum_{q=0}^{4}\Phi_{q}^{\prime}(z)\Big{)}^{2}-\\\
&-\Big{(}\alpha_{2}\psi^{\prime}(\tilde{x}_{2})\sum_{q=0}^{4}\Phi_{q}^{\prime}(z)\Big{)}^{2}+2\sum_{q=0}^{4}\Phi_{q}(z)\cdot\Big{[}f\Big{(}\psi^{-1}\big{(}\frac{z-\alpha_{2}\psi(\tilde{x}_{2})}{\alpha_{1}}\big{)},\tilde{x}_{2}\Big{)}-\\\
&-\Big{(}\alpha_{1}\psi^{\prime\prime}\psi^{-1}\big{(}\frac{z-\alpha_{2}\psi(\tilde{x}_{2})}{\alpha_{1}}\big{)}\sum_{q=0}^{4}\Phi_{q}^{\prime}(z)+\alpha_{1}^{2}\psi^{\prime
2}\psi^{-1}\big{(}\frac{z-\alpha_{2}\psi(\tilde{x}_{2})}{\alpha_{1}}\big{)}\sum_{q=0}^{4}\Phi_{q}^{\prime\prime}(z)\Big{)}-\\\
&-2\Big{(}\alpha_{2}\psi^{\prime\prime}(\tilde{x}_{2})\sum_{q=0}^{4}\Phi_{q}^{\prime}(z)+\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})\sum_{q=0}^{4}\Phi_{q}^{\prime\prime}(z)\Big{)}\bigg{]}dzdx_{2}\longrightarrow
extr\end{split}$ (7.10)

In order to solve problem (7.10) we should take variational derivatives with
respect to each independent variable and equate to zero. First, we should
identify the number of independent variables by variating functional with
respect to each of the functions $\Phi_{q}(z)$ and identify rank of the matrix
of coefficients of the highest derivatives. To make formulas more readable, we
use notation
$\psi^{-1}\big{(}\frac{z-\alpha_{2}\psi(\tilde{x}_{2})}{\alpha_{1}}\big{)}=\tilde{x}_{1}$.
Then, variational derivative with respect to $\Phi_{i}(z)$ is equal to (7.11).

$\begin{split}\int_{0}^{1}&\int_{z_{min}}^{z_{max}}\bigg{[}\frac{f(\tilde{x}_{1},\tilde{x}_{2})}{\alpha_{1}\psi^{\prime}(\tilde{x}_{1})}-\frac{\alpha_{1}^{2}\psi^{\prime
2}(\tilde{x}_{1})+\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})}{\alpha_{1}\psi^{\prime}(\tilde{x}_{1})}\sum_{q=0}^{4}\Phi_{q}^{\prime\prime}(z)-\\\
&-\frac{\alpha_{1}^{2}\psi^{\prime
2}(\tilde{x}_{1})\psi^{\prime\prime}(\tilde{x}_{1})-\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})\psi^{\prime\prime}(\tilde{x}_{1})}{\alpha_{1}^{2}\psi^{\prime
3}(\tilde{x}_{1})}\sum_{q=0}^{4}\Phi_{q}^{\prime}(z)-\\\
&-\Big{(}\frac{\alpha_{1}\alpha_{2}\psi^{\prime
2}(\tilde{x}_{1})\psi^{\prime\prime}(\tilde{x}_{1})\psi^{\prime\prime}(\tilde{x}_{2})}{\alpha_{1}^{3}\psi^{\prime
5}(\tilde{x}_{1})}\\\ \ &-\frac{\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})\big{(}3\psi^{\prime\prime}(\tilde{x}_{1})-\psi(\tilde{x}_{1})\psi^{\prime\prime\prime}(\tilde{x}_{1})\big{)}}{\alpha_{1}^{3}\psi^{\prime
5}(\tilde{x}_{1})}\Big{)}\sum_{q=0}^{4}\Phi_{q}(z)\bigg{]}\delta{\Phi_{i}(z)}dzd\tilde{x}_{2}-\\\
-&\int_{0}^{1}\Bigg{[}\bigg{[}\frac{\alpha_{2}\Big{(}\alpha_{2}\psi^{\prime
2}(\tilde{x}_{2})\psi^{\prime\prime}(\tilde{x}_{1})+\alpha_{1}\psi^{\prime
2}(\tilde{x}_{1})\psi^{\prime\prime}(\tilde{x}_{2})\Big{)}}{\alpha_{1}^{2}\psi^{\prime
3}(\tilde{x}_{1})}\bigg{]}\sum_{q=0}^{4}\Phi_{q}(z)\delta{\Phi_{i}(z)}\Bigg{|}_{z_{min}}^{z_{max}}+\\\
&+\bigg{[}\alpha_{1}\psi^{\prime}(\tilde{x}_{1})+\frac{\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})}{\alpha_{1}\psi^{\prime}(\tilde{x}_{1})}\bigg{]}\sum_{q=0}^{4}\Phi_{q}(z)\delta{\Phi_{i}^{\prime}(z)}\Bigg{|}_{z_{min}}^{z_{max}}\Bigg{]}d\tilde{x}_{2}=0\end{split}$
(7.11)

We can see that coefficient of $\Phi_{i}^{\prime\prime}(z)$ does not depend on
the variable of variation. This implies that rank of the matrix of
coefficients is 1 and that there is only one independent variable. We define
this variable according to formula (7.12)

$U(z)=\sum_{q=0}^{4}\Phi_{q}(z)$ (7.12)

Making this substitution and variating functional with respect to $U(z)$, we
get (7.13)

$\begin{split}\int_{0}^{1}&\int_{z_{min}}^{z_{max}}\bigg{[}\frac{f(\tilde{x}_{1},\tilde{x}_{2})}{\alpha_{1}\psi^{\prime}(\tilde{x}_{1})}-\frac{\alpha_{1}^{2}\psi^{\prime
2}(\tilde{x}_{1})+\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})}{\alpha_{1}\psi^{\prime}(\tilde{x}_{1})}U^{\prime\prime}(z)-\\\
&-\frac{\alpha_{1}^{2}\psi^{\prime
2}(\tilde{x}_{1})\psi^{\prime\prime}(\tilde{x}_{1})-\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})\psi^{\prime\prime}(\tilde{x}_{1})}{\alpha_{1}^{2}\psi^{\prime
3}(\tilde{x}_{1})}U^{\prime}(z)\\\ &-\frac{\alpha_{1}\alpha_{2}\psi^{\prime
2}(\tilde{x}_{1})\psi^{\prime\prime}(\tilde{x}_{1})\psi^{\prime\prime}(\tilde{x}_{2})+\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})\big{(}3\psi^{\prime\prime}(\tilde{x}_{1})-\psi(\tilde{x}_{1})\psi^{\prime\prime\prime}(\tilde{x}_{1})\big{)}}{\alpha_{1}^{3}\psi^{\prime
5}(\tilde{x}_{1})}U(z)\bigg{]}\delta{\Phi_{i}(z)}dzd\tilde{x}_{2}-\\\
-&\int_{0}^{1}\Bigg{[}\bigg{[}\frac{\alpha_{2}\Big{(}\alpha_{2}\psi^{\prime
2}(\tilde{x}_{2})\psi^{\prime\prime}(\tilde{x}_{1})+\alpha_{1}\psi^{\prime
2}(\tilde{x}_{1})\psi^{\prime\prime}(\tilde{x}_{2})\Big{)}}{\alpha_{1}^{2}\psi^{\prime
3}(\tilde{x}_{1})}\bigg{]}U(z)\delta{U(z)}\Bigg{|}_{z_{min}}^{z_{max}}+\\\
&+\bigg{[}\alpha_{1}\psi^{\prime}(\tilde{x}_{1})+\frac{\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})}{\alpha_{1}\psi^{\prime}(\tilde{x}_{1})}\bigg{]}U(z)\delta{U^{\prime}(z)}\Bigg{|}_{z_{min}}^{z_{max}}\Bigg{]}d\tilde{x}_{2}=0\end{split}$
(7.13)

Variation of functional is equal to zero if and only if integrand and boundary
conditions are equal to zero. Then, we get ordinary differential equation of
second degree (7.14)

$\begin{split}&\frac{\alpha_{1}^{2}\psi^{\prime
2}(\tilde{x}_{1})+\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})}{\alpha_{1}\psi^{\prime}(\tilde{x}_{1})}U^{\prime\prime}(z)+\frac{\alpha_{1}^{2}\psi^{\prime
2}(\tilde{x}_{1})\psi^{\prime\prime}(\tilde{x}_{1})-\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})\psi^{\prime\prime}(\tilde{x}_{1})}{\alpha_{1}^{2}\psi^{\prime
3}(\tilde{x}_{1})}U^{\prime}(z)+\\\ &+\frac{\alpha_{1}\alpha_{2}\psi^{\prime
2}(\tilde{x}_{1})\psi^{\prime\prime}(\tilde{x}_{1})\psi^{\prime\prime}(\tilde{x}_{2})+\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})\big{(}3\psi^{\prime\prime}(\tilde{x}_{1})-\psi(\tilde{x}_{1})\psi^{\prime\prime\prime}(\tilde{x}_{1})\big{)}}{\alpha_{1}^{3}\psi^{\prime
5}(\tilde{x}_{1})}U(z)-\\\
&-\frac{f(\tilde{x}_{1},\tilde{x}_{2})}{\alpha_{1}\psi^{\prime}(\tilde{x}_{1})}=0\end{split}$
(7.14)

with boundary conditions (7.15).

$\begin{split}&\bigg{[}\frac{\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})\psi^{\prime\prime}(0)+\alpha_{1}\alpha_{2}\psi^{\prime
2}(0)\psi^{\prime\prime}(\tilde{x}_{2})}{\alpha_{1}^{2}\psi^{\prime
3}(0)}+\alpha_{1}\psi^{\prime}(0)+\frac{\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})}{\alpha_{1}\psi^{\prime}(0)}\bigg{]}U(\alpha_{2}\psi(\tilde{x}_{2}))=0\\\
&\\\ &\bigg{[}\frac{\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})\psi^{\prime\prime}(1)+\alpha_{1}\alpha_{2}\psi^{\prime
2}(1)\psi^{\prime\prime}(\tilde{x}_{2})}{\alpha_{1}^{2}\psi^{\prime
3}(1)}+\alpha_{1}\psi^{\prime}(1)+\frac{\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})}{\alpha_{1}\psi^{\prime}(1)}\bigg{]}U(\alpha_{1}+\alpha_{2}\psi(\tilde{x}_{2}))=0\\\
\end{split}$ (7.15)

Making substitution $U^{\prime}(z)=W(z)$, we get boundary problem of ordinary
differential equations of the first degree (7.16) that can be solved using
standard methods of numerical integration.

$\begin{cases}U^{\prime}(z)=W(z)\\\ \\\ \begin{split}W^{\prime}(z)&=\\\
&=\frac{f(\tilde{x}_{1},\tilde{x}_{2})}{\alpha_{1}^{2}\psi^{\prime
2}(\tilde{x}_{1})+\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})}-\frac{\alpha_{1}^{2}\psi^{\prime
2}(\tilde{x}_{1})\psi^{\prime\prime}(\tilde{x}_{1})-\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})\psi^{\prime\prime}(\tilde{x}_{1})}{\alpha_{1}\psi^{\prime
2}(\tilde{x}_{1})\big{(}\alpha_{1}^{2}\psi^{\prime
2}(\tilde{x}_{1})+\alpha_{2}^{2}\psi^{\prime 2}(\tilde{x}_{2})\big{)}}W(z)-\\\
&-\frac{\alpha_{1}\alpha_{2}\psi^{\prime
2}(\tilde{x}_{1})\psi^{\prime\prime}(\tilde{x}_{1})\psi^{\prime\prime}(\tilde{x}_{2})+\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})\big{(}3\psi^{\prime\prime}(\tilde{x}_{1})-\psi(\tilde{x}_{1})\psi^{\prime\prime\prime}(\tilde{x}_{1})\big{)}}{\alpha_{1}^{2}\psi^{\prime
2}(\tilde{x}_{1})\big{(}\alpha_{1}^{2}\psi^{\prime
2}(\tilde{x}_{1})+\alpha_{2}^{2}\psi^{\prime
2}(\tilde{x}_{2})\big{)}}U(z)\end{split}\end{cases}$ (7.16)

### 7.3 Results of simulation

Chosen Poisson equation is two dimensional ($n=2$), then as stated in KST
$\gamma\geq 6$. We choose $\gamma=10$ and other parameters are defined as
follows: $a=\frac{1}{90}$, $\alpha_{1}=1$, $\alpha_{2}=0.10100010000000001$.
For sake of simplicity, we choose $k=1$, as for this value of $k$, inner
function $\psi(x)$ is a simple identity function. Using linear interpolation,
we get $\psi(x)=x$ as shown in Figure 3.

Figure 3: $\psi(x)$ for $k=1,\gamma=10$. Orange dots are points in which
function $\psi(x)$ id defined. Black line is the result of linear
interpolation.

To solve boundary problem (7.16, 7.15), we use Newton-Raphson method as by
introducing slight modifications it allows us to carry bifurcation analysis
for boundary value problems of ordinary differential equations. The results of
simulation are shown in the Figure 4.

Figure 4: Solutions of Poisson equaton. Blue graph is a solution of system
(7.16, 7.15). Orange graph is a theoretical solution of Poisson equation (7.1,
7.2).

We can see that even using only first term in Taylor series and the simplest
variant of inner function $\psi(x)$, we managed to reduce partial differential
equation to the system of ordinary differential equation with the same
solution.

## 8 Conclusion

In this paper, we introduced a new method for the transition from partial to
ordinary differential equations. In comparison with existing methods, ours
does not use approximation. Instead, it is based on the strict equality -
Kolmogorov superposition theorem. We tested the proposed method on the Poisson
equation and compared the solution of the resulting ODE system with the
solution of the initial partial differential equation. In the zero
approximation solutions coincided.

Further, we plan to test the proposed method on non-linear partial
differential equations with known solution, such as Karman equation. Secondly,
we want to apply this method together with bifurcation analysis to make sure
that the method does not influence the behavior of the dynamic system
generated by the equation.

## References

  * [1] J. Actor.  Computation for the Kolmogorov Superposition Theorem.  PhD thesis, 2018.
  * [2] J. Awrejcewicz, V. Krysko-Jr, L. Kalutsky, M. Zhigalov, and V. Krysko.  Review of the methods of transition from partial to ordinary differential equations: From macro-to nano-structural dynamics.  Archives of Computational Methods in Engineering, pages 1–33, 2021\.
  * [3] T. Hedberg.  The kolmogorov superposition theorem, appendix ii to hs shapiro, topics in approximation theory.  Lecture notes in Math, 187:267–275, 1971.
  * [4] S. I. Khavinson.  Best approximation by linear superpositions (approximate nomography), volume 159.  American Mathematical Soc., 1997.
  * [5] M. Köppen.  On the training of a kolmogorov network.  In International Conference on Artificial Neural Networks, pages 474–479. Springer, 2002.
  * [6] P. Singhatanadgid and T. Singhanart.  The kantorovich method applied to bending, buckling, vibration, and 3d stress analyses of plates: A literature review.  Mechanics of advanced materials and structures, 26(2):170–188, 2019\.
  * [7] D. Sprecher.  A numerical implementation of kolmogorov’s superpositions.  Neural Networks, 9(5):765–772, 1996.
  * [8] D. A. Sprecher.  A numerical implementation of kolmogorov’s superpositions ii.  Neural Networks, 10(3):447–457, 1997.
  * [9] D. A. Sprecher.  From Algebra to Computational Algorithms: Kolmogorov and Hilbert’s Problem 13.  Docent Press, 2017.
  * [10] A. G. Vitushkin.  On hilbert’s thirteenth problem and related questions.  Russian Mathematical Surveys, 59(1):11, 2004.

## Appendix

### A. Variational formulation of Poisson equation

In order to find variational formulation of Poisson equation we start with
functional (A.1).

$\int_{0}^{1}\int_{0}^{1}\Bigg{[}\Big{(}\frac{\partial{u(x_{1},x_{2})}}{\partial{x_{1}}}\Big{)}^{2}+\Big{(}\frac{\partial{u(x_{1},x_{2})}}{\partial{x_{2}}}\Big{)}^{2}-2\sin(\pi
x_{1})\sin(\pi
x_{2})u(x_{1},x_{2})\Bigg{]}dx_{1}dx_{2}\rightarrow\underset{x_{1},x_{2}}{extr}$
(A.1)

To check whether it satisfies partial differential equation (7.1) with
boundary conditions (7.2), we solve the given problem by taking variational
derivative with respect to $u(x_{1},x_{2})$ and equate it to zero.

$\begin{split}\int_{0}^{1}\int_{0}^{1}\Bigg{[}2\frac{\partial{u(x_{1},x_{2})}}{\partial{x_{1}}}\delta{u^{\prime}_{x_{1}}(x_{1},x_{2})}&+2\frac{\partial{u(x_{1},x_{2})}}{\partial{x_{2}}}\delta{u^{\prime}_{x_{2}}(x_{1},x_{2})}-\\\
&-2\sin(\pi x_{1})\sin(\pi
x_{2})\delta{u(x_{1},x_{2})}\Bigg{]}dx_{1}dx_{2}=0\end{split}$ (A.2)

Next, we get rid of derivatives in variations using integration by parts.

$\begin{split}\int_{0}^{1}\int_{0}^{1}\Bigg{[}&2\nabla^{2}{u(x_{1},x_{2})}\delta{u}-2\sin(\pi
x_{1})\sin(\pi x_{2})\delta{u(x_{1},x_{2})}\Bigg{]}dx_{1}dx_{2}+\\\
&+2\int_{0}^{1}\frac{\partial u(x_{1},x_{2})}{\partial
x_{2}}\delta{u(x_{1},x_{2})}\Bigg{|}_{0}^{1}dx_{1}+2\int_{0}^{1}\frac{\partial
u(x_{1},x_{2})}{\partial
x_{1}}\delta{u(x_{1},x_{2})}\Bigg{|}_{0}^{1}dx_{2}=0\end{split}$ (A.3)

Integrand satisfies initial PDE (7.1) but natural boundary conditions do not
match with (7.2). Let’s add boundary terms that do not affect integrand part
but change the boundary conditions so that they coincide with the PDE initial
conditions (7.2). As a result, we get functional (A.4).

$\begin{split}\int_{0}^{1}\int_{0}^{1}\Bigg{[}&\Big{(}\frac{\partial{u(x_{1},x_{2})}}{\partial{x_{1}}}\Big{)}^{2}+\Big{(}\frac{\partial{u(x_{1},x_{2})}}{\partial{x_{2}}}\Big{)}^{2}-2\sin(\pi
x_{1})\sin(\pi x_{2})u(x_{1},x_{2})\Bigg{]}dx_{1}dx_{2}-\\\
&-\int_{0}^{1}2\frac{\partial u(x_{1},x_{2})}{\partial
x_{1}}u(x_{1},x_{2})\Bigg{|}_{0}^{1}dx_{2}-\int_{0}^{1}2\frac{\partial
u(x_{1},x_{2})}{\partial
x_{2}}u(x_{1},x_{2})\Bigg{|}_{0}^{1}dx_{1}\end{split}$ (A.4)

In order to apply our method for transition from partial to ordinary
differential equations, we drag boundary terms under the integral. Finally, we
get equivalent variational formulation for Poisson equation.

$\begin{split}\int_{0}^{1}\int_{0}^{1}\bigg{[}&-\Big{(}\frac{\partial{u(x_{1},x_{2})}}{\partial{x_{1}}}\Big{)}^{2}-\Big{(}\frac{\partial{u(x_{1},x_{2})}}{\partial{x_{2}}}\Big{)}^{2}-2\sin(\pi
x_{1})\sin(\pi x_{2})u(x_{1},x_{2})-\\\
&-2\frac{\partial^{2}u(x_{1},x_{2})}{\partial
x_{1}^{2}}u(x_{1},x_{2})-2\frac{\partial^{2}u(x_{1},x_{2})}{\partial
x_{2}^{2}}u(x_{1},x_{2})\bigg{]}dx_{1}dx_{2}\longrightarrow extr\end{split}$
(A.5)