# Automatic computation of Feynman integrals containing linear propagators via
auxiliary mass flow

Zhi-Feng Liu<EMAIL_ADDRESS>School of Physics and State Key Laboratory
of Nuclear Physics and Technology, Peking University, Beijing 100871, China
Yan-Qing Ma<EMAIL_ADDRESS>School of Physics and State Key Laboratory of
Nuclear Physics and Technology, Peking University, Beijing 100871, China
Center for High Energy Physics, Peking University, Beijing 100871, China
Collaborative Innovation Center of Quantum Matter, Beijing 100871, China

###### Abstract

We proposed a recipe to systematically calculate Feynman integrals containing
linear propagators using the auxiliary mass flow method. The key of the recipe
is to introduce an quadratic term for each linear propagator and then using
differential equations to get rid of their effects. As an application, we
calculated all master integrals of vacuum integrals containing a gauge link up
to 4 loops, and we checked the results by nontrivial dimensional recurrence
relations.

##  I Introduction

As particle physics enters the era of precision, it is necessary to compute
perturbative quantum field theory to high order in coupling constants. To this
end, the evaluation of Feynman integrals is an unavoidable task. The main
strategy at present is to decompose all Feynman integrals in a problem to a
small set of bases, called master integrals (MIs), using integration-by-parts
(IBP) reduction Chetyrkin:1981qh ; Laporta:2000dsw ; Gluza:2010ws ;
Schabinger:2011dz ; vonManteuffel:2012np ; Lee:2013mka ; vonManteuffel:2014ixa
; Larsen:2015ped ; Peraro:2016wsq ; Mastrolia:2018uzb ; Liu:2018dmc ;
Guan:2019bcx ; Klappert:2019emp ; Peraro:2019svx ; Frellesvig:2019kgj ;
Wang:2019mnn ; Smirnov:2019qkx ; Klappert:2020nbg ; Boehm:2020ijp ;
Basat:2021xnn ; Heller:2021qkz ; Bendle:2021ueg , and then calculate these
MIs.

Many methods exist to calculate MIs Hepp:1966eg ; Roth:1996pd ; Binoth:2000ps
; Heinrich:2008si ; Smirnov:2015mct ; Borowka:2015mxa ; Borowka:2017idc ;
Boos:1990rg ; Smirnov:1999gc ; Tausk:1999vh ; Czakon:2005rk ; Smirnov:2009up ;
Gluza:2007rt ; Laporta:2000dsw ; Lee:2009dh ; Kotikov:1990kg ; Kotikov:1991pm
; Remiddi:1997ny ; Gehrmann:1999as ; Argeri:2007up ; MullerStach:2012mp ;
Henn:2013pwa ; Henn:2014qga ; Moriello:2019yhu ; Hidding:2020ytt ;
Catani:2008xa ; Rodrigo:2008fp ; Bierenbaum:2010cy ; Bierenbaum:2012th ;
Tomboulis:2017rvd ; Runkel:2019yrs ; Capatti:2019ypt ; Aguilera-
Verdugo:2020set ; Song:2021vru ; Dubovyk:2022frj ; Liu:2017jxz ; Liu:2020kpc ;
Liu:2021wks ; AMFlow . Among them, the auxiliary mass flow (AMF) method
Liu:2017jxz ; Liu:2020kpc ; Liu:2021wks ; AMFlow , a special case of
differential equations (DEs) method, calculates MIs by setting up and solving
DEs of MIs with respect to an auxiliary mass term $\eta$ (denoted as
$\eta$-DEs). This is systematic and efficient, as far as $\eta$-DEs can be set
up by using reduction strategy. Equipped with an iterative strategy
Liu:2021wks  and the block-triangular reduction form Liu:2018dmc ;
Guan:2019bcx , the AMF method becomes so powerful that many Feynman integrals
in cutting-edge problems, which are very challenging for other methods, can be
calculated (see Fig. 4 in Ref. Liu:2021wks ). At the same time, the AMF method
has already been used in many physical processes Zhang:2018mlo ; Zhang:2020atv
; Yang:2020msy ; Bronnum-Hansen:2020mzk ; Bronnum-Hansen:2021olh ; Wu:2021tzo
; Bronnum-Hansen:2021pqc ; Baranowski:2021gxe .

The AMF method is systematic only if boundary conditions of the $\eta$-DEs at
$\eta\to\infty$ can be systematically obtained. If all inverse propagators
defining Feynman integrals are quadratic in loop momenta, which kind of
Feynman integrals will be called quadratic integrals, then regions at
$\eta\to\infty$ can be systematically identified Beneke:1997zp ;
Smirnov:1999bza  and contribution of each region can be calculated within the
AMF framework LoopDrop . However, life becomes not that easy when some inverse
propagators defining Feynman integrals are linear in loop momenta, which kind
of Feynman integrals will be called linear integrals. In this case,
identification of regions at $\eta\to\infty$ is not that systematic and
calculation of each region becomes nontrivial Zhang:2018mlo ; Zhang:2020atv .
As linear integrals show up frequently in region expansion Beneke:1997zp  and
in effective field theories, it will be good if they can be calculated
systematically by the powerful AMF method.

In this paper, by introducing an auxiliary quadratic term for each linear
propagator, we develop a recipe to calculate any linear integral
systematically within the AMF framework. In the rest of the paper, we first
present the general idea of our recipe and then give some examples to
illustrate its details. As an application, we calculate all master integrals
of 4-loop vacuum integrals containing a gauge link for the first time, which
are useful for studying parton distribution functions Ji:2013dva ; Ji:2014gla
; Radyushkin:2017cyf ; Ma:2014jla ; Ma:2017pxb . All results in this paper
have been checked by using dimensional recurrence relations Tarasov:1996br ;
Tarasov:2000sf ; Lee:2009dh .

##  II The recipe

Consider a general Feynman integral

${I}_{\vec{\nu}}=\prod_{j=1}^{L}\int\frac{\mathrm{d}^{D}\ell_{j}}{\mathrm{i}\pi^{D/2}}\prod_{a=1}^{N}\mathcal{D}_{a}^{-\nu_{a}},$
(1)

where $L$ is the number of loops, $N$ is the number of propagators and
irreducible scalar products, $\mathcal{D}_{a}$ are inverse propagators, and
$\vec{\nu}=(\nu_{1},\nu_{2},\cdots,\nu_{N})$ are integers. We introduce a
corresponding auxiliary quadratic integral ${\widetilde{I}}_{\vec{\nu}}(x)$
defined by

${{\widetilde{I}}}_{\vec{\nu}}(x)=\prod_{j=1}^{L}\int\frac{\mathrm{d}^{D}\ell_{j}}{\mathrm{i}\pi^{D/2}}\prod_{a=1}^{N}\widetilde{\mathcal{D}}_{a}^{-\nu_{a}}.$
(2)

If $\mathcal{D}_{a}$ is an irreducible scalar product or is quadratic in loop
momenta, then $\widetilde{\mathcal{D}}_{a}={\mathcal{D}}_{a}$. Otherwise, it
can be generally expressed as 111If $\mathcal{D}_{a}$ has a negative
infinitesimal imaginary part, we can simply substitute it by
$-\mathcal{D}_{a}$, leaving an overall factor $-1$.

$\displaystyle\mathcal{D}_{a}$
$\displaystyle=p_{a}\cdot\left(\sum_{j=1}^{L}c_{aj}\ell_{j}\right)+\Delta_{a}+\mathrm{i}0^{+}$
$\displaystyle\equiv p_{a}\cdot{\ell}_{a}+\Delta_{a}+\mathrm{i}0^{+},$ (3)

where $p_{a}$ is a fixed Lorentz vector, $\Delta_{a}$ and $c_{aj}$ are fixed
Lorentz scalars, and we then define

$\displaystyle\widetilde{\mathcal{D}}_{a}$
$\displaystyle=x\left({\ell}_{a}^{~{}2}+\mathrm{i}0^{+}\right)+p_{a}\cdot{\ell}_{a}+\Delta_{a}$
$\displaystyle=x\left[\left({\ell}_{a}+\frac{p_{a}}{2x}\right)^{2}+\frac{\Delta_{a}}{x}-\frac{p_{a}^{2}}{4x^{2}}+\mathrm{i}0^{+}\right].$
(4)

Clearly, the original integral can be obtained by

${I}_{\vec{\nu}}=\lim_{x\to 0^{+}}{{\widetilde{I}}}_{\vec{\nu}}(x),$ (5)

which correctly recovers the Feynman prescription.

Denoting MIs of the obtained auxiliary quadratic integrals as
$\overrightarrow{{\widetilde{J}}}(x)$, by using IBP reduction, we can get the
following reduction relation 222With IBP reduction, we in fact only need to
calculate MIs of original linear integrals. Then we can alway choose
${{\widetilde{I}}}_{\vec{\nu}}(x)$ as one of the MIs of auxiliary quadratic
integrals.,

$\displaystyle{{\widetilde{I}}}_{\vec{\nu}}(x)=\sum_{i=1}^{n}c_{\vec{\nu};i}(x){{\widetilde{J}}}_{i}(x),$
(6)

valid for any value of $x$, where $n$ is the number of MIs and
$c_{\vec{\nu};i}(x)$ are rational functions of $x$. And
${{\widetilde{I}}}_{\vec{\nu}}(0^{+})$ can be achieved once
$\overrightarrow{{\widetilde{J}}}(x)$ are known. Again using IBP, we can set
up a system of DEs for MIs with respect to $x$,

$\displaystyle\frac{\partial}{\partial
x}\overrightarrow{{\widetilde{J}}}(x)=A(x){\widetilde{J}}_{(}x),$ (7)

where $A(x)$ is a $n\times n$ matrix with elements being rational functions of
$x$. Because Feynman integrals have Feynman parametric representation, then
$\overrightarrow{{\widetilde{J}}}(x)$ have the following asymptotic expansion
(see e.g. Ref. Henn:2014qga )

${{\widetilde{J}}}_{k}(x)=\sum_{s}\sum_{i=0}^{n_{s}}x^{s}\ln^{i}(x)\sum_{j=0}^{\infty}{{\widetilde{J}}}_{k}^{s\,i\,j}x^{j}\,,$
(8)

where $s$ is a linear function of space-time dimension $D=4-2\epsilon$,
$n_{s}$ is an integer determined by $s$, and ${{\widetilde{J}}}_{k}^{s\,i\,j}$
are functions of $\epsilon$. Taking the limit $x\to 0^{+}$ is then straight
forward 333We assume that dimensional regularization can regularize all
singularities of $x$. Therefore, all terms where $s$ depends on $\epsilon$ can
be ignored as $x\to 0^{+}$. Furthermore, if $s$ is independent of $\epsilon$,
it can be only non-negative integers and $n_{s}=0$.

$\displaystyle{{\widetilde{J}}}_{k}(0^{+})$
$\displaystyle={{\widetilde{J}}}_{k}^{0\,0\,0}\,,$ (9)
$\displaystyle{{\widetilde{J}}}_{k}^{\prime}(0^{+})$
$\displaystyle={{\widetilde{J}}}_{k}^{0\,0\,1}\,,$ (10)

and so on. Then based on Eq. (6), we can get any desired integral
${I}_{\vec{\nu}}={{\widetilde{I}}}_{\vec{\nu}}(0^{+})$.

Taking advantage of the DEs (7), there are only $n$ independent
${{\widetilde{J}}}_{k}^{s\,i\,j}$, which can be fixed by matching the
expansion (8) with values at a given point
$\overrightarrow{{\widetilde{J}}}(x_{0})$, as far as $x_{0}$ is chosen
sufficiently small so that the above expansion is convergent at $x=x_{0}$.
Therefore, once $\overrightarrow{{\widetilde{J}}}(x_{0})$ is known, the
asymptotic expansion of $\overrightarrow{{\widetilde{J}}}(x)$ around $x=0$ can
be achieved and thus ${I}_{\vec{\nu}}$ can be obtained.

For a given regular point $x=x_{0}$, $\overrightarrow{{\widetilde{J}}}(x_{0})$
are quadratic integrals that can be calculated by using the AMF method
Liu:2017jxz ; Liu:2020kpc ; Liu:2021wks ; AMFlow . A convenient choice can be
replacing $\mathrm{i}0^{+}$ by $\mathrm{i}\eta$ in all modified propagators,
i.e., propagators which are linear in loop momenta when $x\to 0^{+}$. The
advantage of this choice is that the number of MIs does not increase after
introducing $\mathrm{i}\eta$. As $\eta\to\infty$, the newly defined integrals
in AMF method are essentially single-mass vacuum integrals Liu:2021wks , which
can be easily calculated within the AMF framework LoopDrop . With this
boundary condition, we can obtain $\overrightarrow{{\widetilde{J}}}(x_{0})$ by
flowing $\eta$ from $\infty$ to $0^{+}$ using $\eta$-DEs.

As numerically solving DEs is very efficient, our method is thus systematic
and efficient as far as DEs can be set up.

##  III A pedagogical example

To elaborate how to use the above proposed recipe to calculate linear
integrals, we take the simplest linear integral with one linear propagator and
one quadratic propagator as an example,

$\displaystyle{I}_{(1,1)}=\int\frac{\mathrm{d}^{D}\ell}{\mathrm{i}\pi^{D/2}}\frac{1}{(\ell^{2}+\mathrm{i}0^{+})(\ell\cdot
n+q\cdot n+\mathrm{i}0^{+})},$ (11)

where $\ell$ is the loop momentum, and $q$ and $n$ are external momenta with
$n^{2}=-1$. Clearly, this integral only depends on the scale $q\cdot n$.
Without loss of generality, we choose $q\cdot n=1$. This Feynman integral is
relevant to the Feynman diagram shown in Fig. 1, where double line stands for
gauge link which contributes linear propagators. We call this a one-loop
bubble integral because there is no external legs of particles.

Figure 1: One-loop vacuum integral containing a a gauge link. Double lines
denote the gauge link.

The corresponding auxiliary quadratic integral is defined as

$\displaystyle{\widetilde{I}}_{(1,1)}(x)=\int\frac{\mathrm{d}^{D}\ell}{\mathrm{i}\pi^{D/2}}\frac{1}{\widetilde{\mathcal{D}}_{1}\widetilde{\mathcal{D}}_{2}},$
(12)

with

$\displaystyle\widetilde{\mathcal{D}}_{1}=$
$\displaystyle\ell^{2}+\mathrm{i}0^{+},$ (13)
$\displaystyle\widetilde{\mathcal{D}}_{2}=$ $\displaystyle
x\left[\left(\ell+\frac{n}{2x}\right)^{2}+\frac{1+4x}{4x^{2}}+\mathrm{i}0^{+}\right].$
(14)

For finite $x$, there are two MIs which can be defined as

$\displaystyle{\widetilde{J}}_{1}(x)$
$\displaystyle=\int\frac{\mathrm{d}^{D}\ell}{\mathrm{i}\pi^{D/2}}\frac{1}{\widetilde{\mathcal{D}}_{2}},$
(15) $\displaystyle{\widetilde{J}}_{2}(x)$
$\displaystyle=\int\frac{\mathrm{d}^{D}\ell}{\mathrm{i}\pi^{D/2}}\frac{x}{\widetilde{\mathcal{D}}_{1}\widetilde{\mathcal{D}}_{2}},$
(16)

with ${\widetilde{I}}_{(1,1)}(x)=\frac{{\widetilde{J}}_{2}(x)}{x}$.

Using IBP, we get the system of $x$-DEs

$\frac{\partial}{\partial
x}\overrightarrow{{\widetilde{J}}}(x)=\left(\begin{array}[]{cc}\frac{2(-1+\epsilon)(1+2x)}{x(1+4x)}&0\\\
-\frac{4(-1+\epsilon)x}{1+4x}&\frac{1}{x}\end{array}\right)\overrightarrow{{\widetilde{J}}}(x).$
(17)

Besides $x=\infty$, there are two additional singularities $x=0$ and
$x=-\frac{1}{4}$, which correspond to infinity mass term and zero mass term in
$\widetilde{\mathcal{D}}_{2}$, respectively. Therefore, we can choose boundary
conditions at, e.g., $x_{0}=\frac{1}{16}$, which point can be well estimated
by asymptotic expansions around $x=0$. Based on DEs (17), we get the following
asymptotic expansions

$\displaystyle{\widetilde{J}}_{1}(x)$
$\displaystyle=x^{-2+2\epsilon}\sum_{i=0}^{\infty}a_{i}x^{i},$ (18)
$\displaystyle{\widetilde{J}}_{2}(x)$
$\displaystyle=c_{0}x+x^{2\epsilon}\sum_{i=0}^{\infty}b_{i}x^{i},$ (19)

which gives

$\displaystyle{\widetilde{I}}_{(1,1)}(x)$
$\displaystyle=\frac{{\widetilde{J}}_{2}(x)}{x}=c_{0}+x^{2\epsilon-1}\sum_{i=0}^{\infty}b_{i}x^{i}.$
(20)

As $x\to 0^{+}$, we find ${\widetilde{J}}_{1}(0^{+})=0$ and
${\widetilde{I}}_{(1,1)}(0^{+})=c_{0}$. The DEs (17) also provide recurrence
relations of coefficients

$\displaystyle a_{i}$ $\displaystyle=\frac{8-4\epsilon-4i}{i}a_{i-1},$ (21)
$\displaystyle b_{i}$
$\displaystyle=-\frac{4(-1+\epsilon)}{-1+2\epsilon+i}a_{i}-\frac{4(-2+2\epsilon+i)}{-1+2\epsilon+i}b_{i-1},$
(22) $\displaystyle b_{0}$
$\displaystyle=-\frac{4(-1+\epsilon)}{-1+2\epsilon}a_{0},$ (23)

where $i\geq 1$. Therefore, only $a_{0}$ and $c_{0}$ are free coefficients.

To fix $a_{0}$ and $c_{0}$, we need to know
$\overrightarrow{{\widetilde{J}}}(\frac{1}{16})$. Using the package AMFlow
AMFlow , we get

$\displaystyle{{\widetilde{J}}_{1}}(\frac{1}{16})=-80.00000000000000\epsilon^{-1}$
$\displaystyle+(316.7393839660332-251.3274122871835\mathrm{i})$
$\displaystyle-(338.0371701002434-995.0661217702465\mathrm{i})\epsilon$
$\displaystyle-(311.267338694263+1888.809135035203\mathrm{i})\epsilon^{2}$
$\displaystyle+(1510.395112926408+2295.761140388504\mathrm{i})\epsilon^{3}$
$\displaystyle-(2656.967479760546+2012.921788373378\mathrm{i})\epsilon^{4}$
$\displaystyle+(3337.360161271413+1359.608268441597\mathrm{i})\epsilon^{5}$
$\displaystyle-(3544.782677804504+737.597502301953\mathrm{i})\epsilon^{6}$
$\displaystyle+(3503.444773310217+330.827170478017\mathrm{i})\epsilon^{7}$
$\displaystyle-(3393.220492080171+125.304259848717\mathrm{i})\epsilon^{8},$
(24)
$\displaystyle{{\widetilde{J}}_{2}}(\frac{1}{16})=1.000000000000000\epsilon^{-1}$
$\displaystyle+(-3.361601777683940+3.141592653589793\mathrm{i})$
$\displaystyle+(2.540630265434342-10.560783449066254\mathrm{i})\epsilon$
$\displaystyle+(4.39414268443412+18.31705093747636\mathrm{i})\epsilon^{2}$
$\displaystyle-(15.75459001722795+20.93897855970841\mathrm{i})\epsilon^{3}$
$\displaystyle+(25.55137360766245+17.56661536372871\mathrm{i})\epsilon^{4}$
$\displaystyle-(31.12100936594129+11.47483328701037\mathrm{i})\epsilon^{5}$
$\displaystyle+(32.43597076695672+6.06440829051716\mathrm{i})\epsilon^{6}$
$\displaystyle-(32.01028557435537+2.66344640964930\mathrm{i})\epsilon^{7}$
$\displaystyle+(30.90646894444592+0.99159094290134\mathrm{i})\epsilon^{8},$
(25)

where we only keep 10 orders in $\epsilon$ expansion with 16-digit precision,
although higher order and more digits can be obtained easily. We will keep the
same accuracy in all examples in this paper. Calculation of
$\overrightarrow{{\widetilde{J}}}(\frac{1}{16})$ in this example is certainly
a very simple problem, but for general complicated problems we can still use
AMF method to calculate them. By matching the asymptotic expansions of
$\overrightarrow{{\widetilde{J}}}(x)$ at $x=0$ with
$\overrightarrow{{\widetilde{J}}}(\frac{1}{16})$, we can determine $a_{0}$ and
$c_{0}$. Eventually we obtain

$\displaystyle{I}_{(1,1)}=c_{0}=-2.000000000000000\epsilon^{-1}$
$\displaystyle+(-0.072979947957153-6.283185307179586\mathrm{i})$
$\displaystyle-(2.356397451352731+0.229273268361557\mathrm{i})\epsilon$
$\displaystyle+(0.19031305777140-28.07369204230733\mathrm{i})\epsilon^{2}$
$\displaystyle-(2.342969999912346+0.156392715314289\mathrm{i})\epsilon^{3}$
$\displaystyle+(0.29757525976318-113.32027706384633\mathrm{i})\epsilon^{4}$
$\displaystyle-(2.334468092829592+0.075946730925850\mathrm{i})\epsilon^{5}$
$\displaystyle+(0.3248443395722-453.6966444901993\mathrm{i})\epsilon^{6}$
$\displaystyle-(2.333380642220398+0.034360215639063\mathrm{i})\epsilon^{7}$
$\displaystyle+(0.331273935869-1814.969761821040\mathrm{i})\epsilon^{8}.$ (26)

##  IV Vacuum integrals containing a gauge link up to 4 loops

Now we apply our recipe to calculate vacuum integrals containing a gauge link
up to 4 loops. These results are useful for studying parton distribution
functions Ji:2013dva ; Ji:2014gla ; Radyushkin:2017cyf ; Ma:2014jla ;
Ma:2017pxb . On the one hand they can serve as boundary conditions for 4-loop
parton correlation functions calculation, and on the other hand they are
needed to renormalize 3-loop parton correlation functions Braun:2018brg ;
Li:2020xml . In general, for a $L$-loop vacuum integral with a gauge link,
there must be at least $L$ quadratic propagators and at most $L$ linear
propagators. Or else, it is either reducible or scaleless. With this
restriction, up to 4 loops we find that there are 64 nonzero and irreducible
sectors in total which give 65 MIs. The 64 nonzero and irreducible sectors are
displayed in the the ancillary file together with results of all the 65 MIs.

Figure 2: The most complex integral family for 4-loop vacuum integrals
containing a gauge link. Double lines denote the gauge link.

We take the most complicated integral family shown in Fig. 2 as an example to
explain our calculation details. All other families can be calculated
similarly and final results are provided in the attached ancillary file.

A general integral in this family can be expressed as

${I}_{\vec{\nu}}=\left(\prod_{j=1}^{4}\int\frac{\mathrm{d}^{D}\ell_{j}}{\mathrm{i}\pi^{D/2}}\right)\frac{\mathcal{D}_{12}^{-\nu_{12}}\mathcal{D}_{13}^{-\nu_{13}}\mathcal{D}_{14}^{-\nu_{14}}}{\mathcal{D}_{1}^{\nu_{1}}\mathcal{D}_{2}^{\nu_{2}}\cdots\mathcal{D}_{11}^{\nu_{11}}},$
(27)

where ${\mathcal{D}_{1},\mathcal{D}_{2},\cdots,\mathcal{D}_{8}}$ are quadratic
denominators, ${\mathcal{D}_{9},\mathcal{D}_{10},\mathcal{D}_{11}}$ are linear
denominators, and ${\mathcal{D}_{12},\mathcal{D}_{13},\mathcal{D}_{14}}$ are
irreducible numerators, with explicit form

$\displaystyle\mathcal{D}_{1}=\ell_{1}^{2},\,\mathcal{D}_{2}=\ell_{2}^{2},\,\mathcal{D}_{3}=\ell_{3}^{2},\,\mathcal{D}_{4}=\ell_{4}^{2},\,\mathcal{D}_{5}=(\ell_{2}+\ell_{4})^{2},$
$\displaystyle\mathcal{D}_{6}=(\ell_{1}+\ell_{2}+\ell_{3})^{2},\,\mathcal{D}_{7}=(\ell_{1}+\ell_{2}+\ell_{4})^{2},$
$\displaystyle\mathcal{D}_{8}=(\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4})^{2},\,\mathcal{D}_{9}=(\ell_{1}+q)\cdot
n,$ $\displaystyle\mathcal{D}_{10}=(\ell_{1}+\ell_{2}+q)\cdot
n,\,\mathcal{D}_{11}=(\ell_{1}+\ell_{2}+\ell_{3}+q)\cdot n,$
$\displaystyle\mathcal{D}_{12}=(\ell_{1}+\ell_{4})^{2},\,\mathcal{D}_{13}=(\ell_{2}+\ell_{3})^{2},\,\mathcal{D}_{14}=\ell_{4}\cdot
n,$

where Feynman prescription $+\mathrm{i}0^{+}$ is omitted. This family has 32
MIs in total. Similar to one-loop case, the single scale $q\cdot n$ can be
extracted as an overall factor, and thus we can simply set $q\cdot n=1$.

Based on our recipe, we define a corresponding auxiliary family of integrals

${\widetilde{I}}_{\vec{\nu}}(x)=\left(\prod_{j=1}^{4}\int\frac{\mathrm{d}^{D}\ell_{j}}{\mathrm{i}\pi^{D/2}}\right)\frac{\widetilde{\mathcal{D}}_{12}^{-\nu_{12}}\widetilde{\mathcal{D}}_{13}^{-\nu_{13}}\widetilde{\mathcal{D}}_{14}^{-\nu_{14}}}{\widetilde{\mathcal{D}}_{1}^{\nu_{1}}\widetilde{\mathcal{D}}_{2}^{\nu_{2}}\cdots\widetilde{\mathcal{D}}_{11}^{\nu_{11}}},$
(28)

where

$\displaystyle\widetilde{\mathcal{D}}_{9}$
$\displaystyle={x}\left[{(\ell_{1}+\frac{n}{2x})^{2}+\frac{1+4x}{4x^{2}}}+\mathrm{i}0^{+}\right],$
$\displaystyle\widetilde{\mathcal{D}}_{10}$
$\displaystyle={x}\left[{(\ell_{1}+\ell_{2}+\frac{n}{2x})^{2}+\frac{1+4x}{4x^{2}}}+\mathrm{i}0^{+}\right],$
$\displaystyle\widetilde{\mathcal{D}}_{11}$
$\displaystyle={x}\left[{(\ell_{1}+\ell_{2}+\ell_{3}+\frac{n}{2x})^{2}+\frac{1+4x}{4x^{2}}}+\mathrm{i}0^{+}\right],$

and $\widetilde{\mathcal{D}}_{i}={\mathcal{D}}_{i}$ for $i\neq 9,10,11$. After
changing linear propagators to quadratic propagators, the number of MIs in the
auxiliary family increases to 120. We use
$\overrightarrow{{\widetilde{J}}}(x)$ to denote these MIs. Each
${\widetilde{J}}_{i}(x)$ is defined by a specific list of exponents
$\vec{\nu}^{i}$ via the following relation

$\displaystyle{\widetilde{J}}_{i}(x)=x^{\nu_{9}^{i}+\nu_{10}^{i}+\nu_{11}^{i}}{\widetilde{I}}_{\vec{\nu}^{i}}(x),$
(29)

where, for convenience, the factor $x^{\nu_{9}^{i}+\nu_{10}^{i}+\nu_{11}^{i}}$
is introduced to cancel the overall factor $x$ in
$\widetilde{\mathcal{D}}_{9},\,\widetilde{\mathcal{D}}_{10},\,\widetilde{\mathcal{D}}_{11}$.

For a regular point, say $x_{0}=\frac{1}{16}$, we calculate the corresponding
MIs $\overrightarrow{{\widetilde{J}}}(\frac{1}{16})$ using the package
AMFlowAMFlow  with the mass mode at the first step of the iteration. The
calculation is automatic and can be decomposed to the following steps. First,
$\mathrm{i}0^{+}$ in $\widetilde{\mathcal{D}}_{9/10/11}$ is replaced by
$\mathrm{i}\eta$. Then, $\eta$-DEs are set up and boundary conditions at
$\eta\to\infty$ are calculated Liu:2021wks ; LoopDrop , where
KiraKlappert:2020nbg  is employed to do IBP reduction. And finally, the
$\eta$-DEs are solved to realize the flow of $\eta$ from $\infty$ to $0^{+}$.

In this problem, $x$-DEs of $\overrightarrow{{\widetilde{J}}}(x)$ can be
derived from $\eta$-DEs during calculating
$\overrightarrow{{\widetilde{J}}}(\frac{1}{16})$, details of which are given
in the Appendix A. With boundary conditions at $x_{0}=\frac{1}{16}$, the
$x$-DEs are solved (similar to solve $\eta$-DEs) by using AMFlowAMFlow  to
obtain results at $x=0^{+}$, which gives final results of our desired linear
integrals. Following is the result of the top-sector corner integral

$\displaystyle{I}_{(1,\cdots,1,0,0,0)}=2.467401100272340\epsilon^{-4}$
$\displaystyle-(2.10725949249774-31.00627668029982\mathrm{i})\epsilon^{-3}$
$\displaystyle-(142.4061631115384+26.4806037633531\mathrm{i})\epsilon^{-2}$
$\displaystyle+(258.9104506294486-157.4236379069742\mathrm{i})\epsilon^{-1}$
$\displaystyle+(-1638.226426723546+1859.681566826831\mathrm{i})$
$\displaystyle-(4750.006669884407+11690.896001211223\mathrm{i})\epsilon$
$\displaystyle+(22176.46421510774+23525.15655373790\mathrm{i})\epsilon^{2}$
$\displaystyle-(120064.4744857791+165534.0450504650\mathrm{i})\epsilon^{3}$
$\displaystyle+(855055.4310121035+539390.2035592209\mathrm{i})\epsilon^{4}$
$\displaystyle-(3796471.438494197+1678278.567906508\mathrm{i})\epsilon^{5}.$
(30)

Note that higher precision and higher orders in $\epsilon$ can also be
achieved easily.

##  V Integrals with external legs

Now let us consider more general linear integrals containing both gauge links
and external legs. Take 4-loop integrals containing one external leg shown in
Fig. 3 as an example. The top-sector corner integral can be expressed as

${I}_{(1,1,1,1,1,1)}(q\cdot n,p\cdot
n)=\left(\prod_{j=1}^{4}\int\frac{\mathrm{d}^{D}\ell_{j}}{\mathrm{i}\pi^{D/2}}\right)\frac{1}{\mathcal{D}_{1}\mathcal{D}_{2}\mathcal{D}_{3}\mathcal{D}_{4}\mathcal{D}_{5}\mathcal{D}_{6}},$
(31)

with

$\displaystyle\mathcal{D}_{1}=\ell_{1}^{2},\,\mathcal{D}_{2}=\ell_{2}^{2},\,\mathcal{D}_{3}=\ell_{3}^{2},\,\mathcal{D}_{4}=\ell_{4}^{2},\,\mathcal{D}_{5}=(\ell_{1}+q)\cdot
n,$
$\displaystyle\mathcal{D}_{6}=(\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4}+q-p)\cdot
n,$ (32)

where $p$ is a light-like external momentum, and Feynman prescription
$+\mathrm{i}0^{+}$ for each inverse propagator is omitted.

Figure 3: 4-loop Feynman diagram with a gauge link and a external leg.

To calculate integrals in this example, we define auxiliary quadratic
integrals by using the following inverse propagators

$\displaystyle\widetilde{\mathcal{D}}_{5}$
$\displaystyle={x}\left[{(\ell_{1}+\frac{n}{2x})^{2}+\frac{1}{4x^{2}}}+\frac{q\cdot
n}{x}+\mathrm{i}0^{+}\right],$ $\displaystyle\widetilde{\mathcal{D}}_{6}$
$\displaystyle={x}\bigg{[}(\ell_{1}+\ell_{2}+\ell_{3}+\ell_{4}+\frac{n}{2x})^{2}$
$\displaystyle\quad\quad+\frac{1}{4x^{2}}+\frac{(q-p)\cdot
n}{x}+\mathrm{i}0^{+}\bigg{]},$

and $\widetilde{\mathcal{D}}_{i}={\mathcal{D}}_{i}$ for $i\leq 4$. After
changing linear propagators to quadratic propagators, the number of MIs
increases from 2 to 6. Using AMFlow AMFlow , we calculate auxiliary quadratic
integrals at a fixed point $x=x_{0}$, and then solve $x$-DEs to obtain
integrals at $x=0^{+}$. Without loss of generality, for the phase-space point
$q\cdot n=1$, $p\cdot n=\frac{1}{2}$, we obtain

$\displaystyle{I}_{(1,1,1,1,1,1)}(1,\frac{1}{2})=-0.001388888888888889\epsilon^{-2}$
$\displaystyle-(0.01673049985543654+0.01745329251994330\mathrm{i})\epsilon^{-1}$
$\displaystyle+(-0.2181793358910845-0.2102416617468981\mathrm{i})$
$\displaystyle-(1.816593593006501+3.660426889433229\mathrm{i})\epsilon$
$\displaystyle-(15.81900673709547+33.89466583909418\mathrm{i})\epsilon^{2}$
$\displaystyle-(110.3704374791621+401.1363831252734\mathrm{i})\epsilon^{3}$
$\displaystyle-(810.083299387532+3287.604829797105\mathrm{i})\epsilon^{4}$
$\displaystyle-(5155.23653311233+33468.66058680208\mathrm{i})\epsilon^{5}$
$\displaystyle-(34800.6654379141+258370.0763552148\mathrm{i})\epsilon^{6}$
$\displaystyle-(211074.015357447+2454130.498362980\mathrm{i})\epsilon^{7}.$
(33)

When there are more gauge links or more external legs, the calculation process
is similar.

##  VI Self-consistency check by dimensional recurrence relations

For any integral ${I}_{\vec{\nu}}^{D}$, where we explicitly write down the
dependence of spacetime dimension $D$, there is a parametric representation
(see, e.g., Bitoun:2017nre )

$\displaystyle{I}_{\vec{\nu}}^{D}$
$\displaystyle=\left(\prod_{j=1}^{L}\int\frac{\mathrm{d}^{D}\ell_{j}}{\mathrm{i}\pi^{D/2}}\right)\prod_{i=1}^{N}\mathcal{D}_{i}^{-\nu_{i}}$
$\displaystyle=\left(\prod_{i=1}^{N}\int_{0}^{\infty}\frac{x_{i}^{\nu_{i}-1}\mathrm{d}x_{i}}{\Gamma(\nu_{i})}\right)\frac{e^{-F/U}}{U^{D/2}},$
(34)

and thus

${I}_{\vec{\nu}}^{D-2}=\left(\prod_{i=1}^{N}\int_{0}^{\infty}\frac{x_{i}^{\nu_{i}-1}\mathrm{d}x_{i}}{\Gamma(\nu_{i})}\right)\frac{e^{-F/U}}{U^{D/2}}U.$
(35)

As $U$ is a homogeneous polynomial of $x_{i}$, the above expression can be
translated into a linear combination of integrals at dimension $D$ with
different $\vec{\nu}$, which can be reduced to MIs by using IBP identities.
Supposing $\vec{J}^{D}$ is a set of MIs, the above procedure results in
dimensional recurrence relations Tarasov:1996br ; Tarasov:2000sf

$\vec{J}^{D-2}=C^{D}\vec{J}^{D},$ (36)

where $C^{D}$ is a matrix rationally depending on $D$.

Up to 3 loops, vacuum integrals containing a gauge link have been calculated
in Ref. Li:2020xml  by solving dimensional recurrence relations using the
method in Ref. Lee:2009dh . The obtained results agree with ours. For 4-loop
integrals containing a gauge link, we find it very hard to solve the
corresponding dimensional recurrence relations. However, the dimensional
recurrence relations can provide a highly nontrivial self-consistency check
for our results. We calculate MIs at two different space-time dimensions
$D_{0}$ and $D_{0}-2$, and we find that the obtained results satisfy the
relations Eq. (36) within precision.

##  VII Summary

In summary, we develop a recipe to calculate linear integrals using the AMF
method Liu:2017jxz ; Liu:2020kpc ; Liu:2021wks . For any given linear
integral, our recipe is to introduce an auxiliary quadratic term for each
linear propagator, and the obtained auxiliary quadratic integral can be
calculated systematically using the AMF method. Taking result of the auxiliary
quadratic integral calculated at fixed auxiliary quadratic terms as boundary
condition and using differential equations to push the auxiliary quadratic
terms to zero, effects of auxiliary quadratic terms will die out eventually
and we get result of the target linear integral. This recipe of calculating
linear integrals is very systematic and will be implemented in the package
AMFlow AMFlow .

As linear integrals show up frequently in region expansion and in effective
field theories, our recipe will be useful in phenomenological studies. As the
first application, we have calculated all MIs of vacuum integrals containing a
gauge link up to 4 loops, which are useful to study parton distribution
functions. Our results have been checked by nontrivial dimensional recurrence
relations.

##  VIII Acknowledgments

We thank K.T. Chao, X. Liu, Z.Y. Li, X. Li, X. Guan, W.H. Wu, C. Meng, R.H. Wu
for many useful communications and discussions. The work is supported in part
by the National Natural Science Foundation of China (Grants No. 11875071, No.
11975029), the National Key Research and Development Program of China under
Contracts No. 2020YFA0406400, and the High-performance Computing Platform of
Peking University.

##  Appendix A Derivation of relation between $x$-DEs and $\eta$-DEs for
vacuum integrals

For vacuum integrals containing a gauge link, $x$-DEs can be obtained from
$\eta$-DEs. We take the family in Fig. 1 as an example to explain this, and
the final relation is independent of the specific family. MIs after
introducing $x$ are defined by
$\overrightarrow{{\widetilde{J}}}(x)=\left(\widetilde{J}_{1}(x),\widetilde{J}_{2}(x)\right)^{T}$,
with

$\displaystyle\widetilde{J}_{1}(x)=\int\frac{\mathrm{d}^{D}\ell}{\mathrm{i}\pi^{D/2}}\frac{1}{(\ell+\frac{n}{2x})^{2}+\frac{1+4x}{4x^{2}}},$
(37)
$\displaystyle\widetilde{J}_{2}(x)=\int\frac{\mathrm{d}^{D}\ell}{\mathrm{i}\pi^{D/2}}\frac{1}{\ell^{2}\left((\ell+\frac{n}{2x})^{2}+\frac{1+4x}{4x^{2}}\right)},$
(38)

where Feynman prescription $+\mathrm{i}0^{+}$ in each denominator is omitted.
The corresponding $x$-DEs are denoted as

$\displaystyle\frac{\partial}{\partial
x}\overrightarrow{{\widetilde{J}}}(x)=A(x)\overrightarrow{{\widetilde{J}}}(x).$
(39)

When calculating $\overrightarrow{{\widetilde{J}}}(x_{0})$ using AMF method,
we define modified integrals
$\overrightarrow{{\widehat{J}}}(\eta)=\left(\widehat{J}_{1}(\eta),\widehat{J}_{2}(\eta)\right)^{T}$,
with

$\displaystyle\widehat{J}_{1}(\eta)=\int\frac{\mathrm{d}^{D}\ell}{\mathrm{i}\pi^{D/2}}\frac{1}{(\ell+\frac{n}{2x_{0}})^{2}+\frac{1+4x_{0}}{4x_{0}^{2}}-\eta},$
(40)
$\displaystyle\widehat{J}_{2}(\eta)=\int\frac{\mathrm{d}^{D}\ell}{\mathrm{i}\pi^{D/2}}\frac{1}{\ell^{2}\left((\ell+\frac{n}{2x_{0}})^{2}+\frac{1+4x_{0}}{4x_{0}^{2}}-\eta\right)}.$
(41)

Note that in the above we introducing $-\eta$ instead of $\mathrm{i}\eta$ for
convenience, which is only a matter of convention. The corresponding
$\eta$-DEs are denoted as

$\displaystyle\frac{\partial}{\partial\eta}\overrightarrow{{\widehat{J}}}(\eta)=M(\eta)\overrightarrow{{\widehat{J}}}(\eta).$
(42)

We will express $A(x)$ by $M(\eta)$.

Let us define
$\overrightarrow{K}(p^{2},m^{2})=\left(K_{1}(p^{2},m^{2}),K_{2}(p^{2},m^{2})\right)^{T}$
with

$\displaystyle
K_{1}(p^{2},m^{2})=\int\frac{\mathrm{d}^{D}\ell}{\mathrm{i}\pi^{D/2}}\frac{1}{(\ell+p)^{2}+m^{2}},$
(43) $\displaystyle
K_{2}(p^{2},m^{2})=\int\frac{\mathrm{d}^{D}\ell}{\mathrm{i}\pi^{D/2}}\frac{1}{\ell^{2}\left((\ell+p)^{2}+m^{2}\right)}.$
(44)

Because Feynman integrals have definite mass dimension, there is a rescaling
relation

$\displaystyle\overrightarrow{K}((1+\delta)p^{2},(1+\delta)m^{2})=(1+\delta)^{\alpha}\overrightarrow{K}(p^{2},m^{2}),$
(45)

where $\alpha$ is a diagonal matrix with diagonal elements denoting mass
dimension. In this example it is given by

$\alpha=\left(\begin{array}[]{cc}\frac{D}{2}-1&0\\\
0&\frac{D}{2}-2\end{array}\right).$ (46)

Expansion of Eq. (45) to $O(\delta)$ gives

$\displaystyle\left(p^{2}\frac{\partial}{\partial
p^{2}}+m^{2}\frac{\partial}{\partial
m^{2}}-\alpha\right)\overrightarrow{K}(p^{2},m^{2})=0.$ (47)

For convenience, we denote

$\displaystyle y$ $\displaystyle=-\frac{1}{4x^{2}},\,z=\frac{1+4x}{4x^{2}},$
$\displaystyle y_{0}$
$\displaystyle=-\frac{1}{4x_{0}^{2}},\,z_{0}=\frac{1+4x_{0}}{4x_{0}^{2}}.$

Then integrals are related by

$\displaystyle\overrightarrow{{\widetilde{J}}}(x)$
$\displaystyle=\overrightarrow{K}(y,z),$ (48)
$\displaystyle\overrightarrow{{\widehat{J}}}(\eta)$
$\displaystyle=\overrightarrow{K}(y_{0},z_{0}-\eta).$ (49)

Based on $\eta$-DEs (42) we get

$\displaystyle\frac{\partial}{\partial\eta}\overrightarrow{K}(y_{0},z_{0}-\eta)$
$\displaystyle=M(\eta)\overrightarrow{K}(y_{0},z_{0}-\eta).$ (50)

By choosing $\eta=z_{0}-z$ in the above equation, we get

$\displaystyle\frac{\partial}{\partial z}\overrightarrow{K}(y_{0},z)$
$\displaystyle=-M(z_{0}-z)\overrightarrow{K}(y_{0},z),$ (51)

and thus

$\displaystyle\frac{\partial}{\partial z}\overrightarrow{K}(y,z)$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
z}\left[\left(\frac{y}{y_{0}}\right)^{\alpha}\overrightarrow{K}(y_{0},\frac{y_{0}}{y}z)\right]$
$\displaystyle=$
$\displaystyle-\left(\frac{y}{y_{0}}\right)^{\alpha-1}M(z_{0}-\frac{y_{0}}{y}z)\overrightarrow{K}(y_{0},\frac{y_{0}}{y}z)$
$\displaystyle=$
$\displaystyle-\left(\frac{y}{y_{0}}\right)^{\alpha-1}M(z_{0}-\frac{y_{0}}{y}z)\left(\frac{y}{y_{0}}\right)^{-\alpha}\overrightarrow{K}(y,z).$
(52)

where Eq. (45) has been used. Therefore we have

$\displaystyle\frac{\partial}{\partial
x}\overrightarrow{{\widetilde{J}}}(x)=\left(\frac{\partial y}{\partial
x}\frac{\partial}{\partial y}+\frac{\partial z}{\partial
x}\frac{\partial}{\partial z}\right)\overrightarrow{K}(y,z)$
$\displaystyle=\left[\frac{\partial y}{\partial
x}\left(\frac{\alpha}{y}-\frac{z}{y}\frac{\partial}{\partial
z}\right)+\frac{\partial z}{\partial x}\frac{\partial}{\partial
z}\right]\overrightarrow{K}(y,z)$ $\displaystyle=\frac{\partial y}{\partial
x}\frac{\alpha}{y}\overrightarrow{{\widetilde{J}}}(x)-\left(\frac{y}{y_{0}}\right)^{\alpha-1}M(z_{0}-\frac{y_{0}}{y}z)\left(\frac{y}{y_{0}}\right)^{-\alpha}$
$\displaystyle\hskip 56.9055pt\times\left[-\frac{z}{y}\frac{\partial
y}{\partial x}+\frac{\partial z}{\partial
x}\right]\overrightarrow{{\widetilde{J}}}(x),$ (53)

where the Eq. (47) has been used. The above result means

$\displaystyle A(x)=$ $\displaystyle\frac{1}{y}\frac{\partial y}{\partial
x}\alpha+\frac{z^{2}}{y}\left(\frac{\partial}{\partial x}\frac{y}{z}\right)$
$\displaystyle\times\left(\frac{y}{y_{0}}\right)^{\alpha-1}M\left(z_{0}-y_{0}\frac{z}{y}\right)\left(\frac{y}{y_{0}}\right)^{-\alpha},$
(54)

which connects $x$-DEs to $\eta$-DEs for any vacuum integrals containing a
gauge link. It is useful because $M$ can be more easily obtained.

## References

  * (1) K. G. Chetyrkin and F. V. Tkachov, Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops, Nucl. Phys. B192 (1981) 159–204 [InSPIRE].
  * (2) S. Laporta, High precision calculation of multiloop Feynman integrals by difference equations, Int. J. Mod. Phys. A 15 (2000) 5087–5159 [hep-ph/0102033] [InSPIRE].
  * (3) J. Gluza, K. Kajda, and D. A. Kosower, Towards a Basis for Planar Two-Loop Integrals, Phys. Rev. D83 (2011) 045012 [arXiv:1009.0472] [InSPIRE].
  * (4) R. M. Schabinger, A New Algorithm For The Generation Of Unitarity-Compatible Integration By Parts Relations, JHEP 01 (2012) 077 [arXiv:1111.4220] [InSPIRE].
  * (5) A. von Manteuffel and C. Studerus, Reduze 2 - Distributed Feynman Integral Reduction, [arXiv:1201.4330] [InSPIRE].
  * (6) R. N. Lee, LiteRed 1.4: a powerful tool for reduction of multiloop integrals, J. Phys. Conf. Ser. 523 (2014) 012059 [arXiv:1310.1145] [InSPIRE].
  * (7) A. von Manteuffel and R. M. Schabinger, A novel approach to integration by parts reduction, Phys. Lett. B744 (2015) 101–104 [arXiv:1406.4513] [InSPIRE].
  * (8) K. J. Larsen and Y. Zhang, Integration-by-parts reductions from unitarity cuts and algebraic geometry, Phys. Rev. D93 (2016) 041701 [arXiv:1511.01071] [InSPIRE].
  * (9) T. Peraro, Scattering amplitudes over finite fields and multivariate functional reconstruction, JHEP 12 (2016) 030 [arXiv:1608.01902] [InSPIRE].
  * (10) P. Mastrolia and S. Mizera, Feynman Integrals and Intersection Theory, JHEP 02 (2019) 139 [arXiv:1810.03818] [InSPIRE].
  * (11) X. Liu and Y.-Q. Ma, Determining arbitrary Feynman integrals by vacuum integrals, Phys. Rev. D 99 (2019) 071501 [arXiv:1801.10523] [InSPIRE].
  * (12) X. Guan, X. Liu, and Y.-Q. Ma, Complete reduction of two-loop five-light-parton scattering amplitudes, Chin. Phys. C 44 (2020) 9 [arXiv:1912.09294] [InSPIRE].
  * (13) J. Klappert and F. Lange, Reconstructing rational functions with FireFly, Comput. Phys. Commun. 247 (2020) 106951 [arXiv:1904.00009] [InSPIRE].
  * (14) T. Peraro, FiniteFlow: multivariate functional reconstruction using finite fields and dataflow graphs, JHEP 07 (2019) 031 [arXiv:1905.08019] [InSPIRE].
  * (15) H. Frellesvig, F. Gasparotto, S. Laporta, M. K. Mandal, P. Mastrolia, L. Mattiazzi, and S. Mizera, Decomposition of Feynman Integrals on the Maximal Cut by Intersection Numbers, JHEP 05 (2019) 153 [arXiv:1901.11510] [InSPIRE].
  * (16) Y. Wang, Z. Li, and N. Ul Basat, Direct reduction of multiloop multiscale scattering amplitudes, Phys. Rev. D 101 (2020) 076023 [arXiv:1901.09390] [InSPIRE].
  * (17) A. V. Smirnov and F. S. Chuharev, FIRE6: Feynman Integral REduction with Modular Arithmetic, Comput. Phys. Commun. 247 (2020) 106877 [arXiv:1901.07808] [InSPIRE].
  * (18) J. Klappert, F. Lange, P. Maierhöfer, and J. Usovitsch, Integral reduction with Kira 2.0 and finite field methods, Comput. Phys. Commun. 266 (2021) 108024 [arXiv:2008.06494] [InSPIRE].
  * (19) J. Boehm, M. Wittmann, Z. Wu, Y. Xu, and Y. Zhang, IBP reduction coefficients made simple, JHEP 12 (2020) 054 [arXiv:2008.13194] [InSPIRE].
  * (20) N. u. Basat, Z. Li, and Y. Wang, Reduction of the planar double-box diagram for single-top production via auxiliary mass flow, Phys. Rev. D 104 (2021) 056020 [arXiv:2102.08225] [InSPIRE].
  * (21) M. Heller and A. von Manteuffel, MultivariateApart: Generalized partial fractions, Comput. Phys. Commun. 271 (2022) 108174 [arXiv:2101.08283] [InSPIRE].
  * (22) D. Bendle, J. Boehm, M. Heymann, R. Ma, M. Rahn, L. Ristau, M. Wittmann, Z. Wu, and Y. Zhang, Two-loop five-point integration-by-parts relations in a usable form, [arXiv:2104.06866] [InSPIRE].
  * (23) K. Hepp, Proof of the Bogolyubov-Parasiuk theorem on renormalization, Commun. Math. Phys. 2 (1966) 301–326 [InSPIRE].
  * (24) M. Roth and A. Denner, High-energy approximation of one loop Feynman integrals, Nucl. Phys. B 479 (1996) 495–514 [hep-ph/9605420] [InSPIRE].
  * (25) T. Binoth and G. Heinrich, An automatized algorithm to compute infrared divergent multiloop integrals, Nucl. Phys. B585 (2000) 741–759 [hep-ph/0004013] [InSPIRE].
  * (26) G. Heinrich, Sector Decomposition, Int. J. Mod. Phys. A23 (2008) 1457–1486 [arXiv:0803.4177] [InSPIRE].
  * (27) A. V. Smirnov, FIESTA4: Optimized Feynman integral calculations with GPU support, Comput. Phys. Commun. 204 (2016) 189–199 [arXiv:1511.03614] [InSPIRE].
  * (28) S. Borowka, G. Heinrich, S. P. Jones, M. Kerner, J. Schlenk, and T. Zirke, SecDec-3.0: numerical evaluation of multi-scale integrals beyond one loop, Comput. Phys. Commun. 196 (2015) 470–491 [arXiv:1502.06595] [InSPIRE].
  * (29) S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner, J. Schlenk, and T. Zirke, pySecDec: a toolbox for the numerical evaluation of multi-scale integrals, Comput. Phys. Commun. 222 (2018) 313–326 [arXiv:1703.09692] [InSPIRE].
  * (30) E. E. Boos and A. I. Davydychev, A Method of evaluating massive Feynman integrals, Theor. Math. Phys. 89 (1991) 1052–1063 [InSPIRE]. [Teor. Mat. Fiz.89,56(1991)].
  * (31) V. A. Smirnov, Analytical result for dimensionally regularized massless on shell double box, Phys. Lett. B460 (1999) 397–404 [hep-ph/9905323] [InSPIRE].
  * (32) J. B. Tausk, Nonplanar massless two loop Feynman diagrams with four on-shell legs, Phys. Lett. B 469 (1999) 225–234 [hep-ph/9909506] [InSPIRE].
  * (33) M. Czakon, Automatized analytic continuation of Mellin-Barnes integrals, Comput. Phys. Commun. 175 (2006) 559–571 [hep-ph/0511200] [InSPIRE].
  * (34) A. V. Smirnov and V. A. Smirnov, On the Resolution of Singularities of Multiple Mellin-Barnes Integrals, Eur. Phys. J. C 62 (2009) 445–449 [arXiv:0901.0386] [InSPIRE].
  * (35) J. Gluza, K. Kajda, and T. Riemann, AMBRE: A Mathematica package for the construction of Mellin-Barnes representations for Feynman integrals, Comput. Phys. Commun. 177 (2007) 879–893 [arXiv:0704.2423] [InSPIRE].
  * (36) R. N. Lee, Space-time dimensionality D as complex variable: Calculating loop integrals using dimensional recurrence relation and analytical properties with respect to D, Nucl. Phys. B830 (2010) 474–492 [arXiv:0911.0252] [InSPIRE].
  * (37) A. V. Kotikov, Differential equations method: New technique for massive Feynman diagrams calculation, Phys. Lett. B254 (1991) 158–164 [InSPIRE].
  * (38) A. V. Kotikov, Differential equation method: The Calculation of N point Feynman diagrams, Phys. Lett. B 267 (1991) 123–127 [InSPIRE]. [Erratum: Phys.Lett.B 295, 409–409 (1992)].
  * (39) E. Remiddi, Differential equations for Feynman graph amplitudes, Nuovo Cim. A110 (1997) 1435–1452 [hep-th/9711188] [InSPIRE].
  * (40) T. Gehrmann and E. Remiddi, Differential equations for two loop four point functions, Nucl. Phys. B580 (2000) 485–518 [hep-ph/9912329] [InSPIRE].
  * (41) M. Argeri and P. Mastrolia, Feynman Diagrams and Differential Equations, Int. J. Mod. Phys. A 22 (2007) 4375–4436 [arXiv:0707.4037] [InSPIRE].
  * (42) S. Müller-Stach, S. Weinzierl, and R. Zayadeh, Picard-Fuchs equations for Feynman integrals, Commun. Math. Phys. 326 (2014) 237–249 [arXiv:1212.4389] [InSPIRE].
  * (43) J. M. Henn, Multiloop integrals in dimensional regularization made simple, Phys. Rev. Lett. 110 (2013) 251601 [arXiv:1304.1806] [InSPIRE].
  * (44) J. M. Henn, Lectures on differential equations for Feynman integrals, J. Phys. A48 (2015) 153001 [arXiv:1412.2296] [InSPIRE].
  * (45) F. Moriello, Generalised power series expansions for the elliptic planar families of Higgs + jet production at two loops, JHEP 01 (2020) 150 [arXiv:1907.13234] [InSPIRE].
  * (46) M. Hidding, DiffExp, a Mathematica package for computing Feynman integrals in terms of one-dimensional series expansions, Comput. Phys. Commun. 269 (2021) 108125 [arXiv:2006.05510] [InSPIRE].
  * (47) S. Catani, T. Gleisberg, F. Krauss, G. Rodrigo, and J.-C. Winter, From loops to trees by-passing Feynman’s theorem, JHEP 09 (2008) 065 [arXiv:0804.3170] [InSPIRE].
  * (48) G. Rodrigo, S. Catani, T. Gleisberg, F. Krauss, and J.-C. Winter, From multileg loops to trees (by-passing Feynman’s Tree Theorem), Nucl. Phys. B Proc. Suppl. 183 (2008) 262–267 [arXiv:0807.0531] [InSPIRE].
  * (49) I. Bierenbaum, S. Catani, P. Draggiotis, and G. Rodrigo, A Tree-Loop Duality Relation at Two Loops and Beyond, JHEP 10 (2010) 073 [arXiv:1007.0194] [InSPIRE].
  * (50) I. Bierenbaum, S. Buchta, P. Draggiotis, I. Malamos, and G. Rodrigo, Tree-Loop Duality Relation beyond simple poles, JHEP 03 (2013) 025 [arXiv:1211.5048] [InSPIRE].
  * (51) E. T. Tomboulis, Causality and Unitarity via the Tree-Loop Duality Relation, JHEP 05 (2017) 148 [arXiv:1701.07052] [InSPIRE].
  * (52) R. Runkel, Z. Szőr, J. P. Vesga, and S. Weinzierl, Causality and loop-tree duality at higher loops, Phys. Rev. Lett. 122 (2019) 111603 [arXiv:1902.02135] [InSPIRE]. [Erratum: Phys.Rev.Lett. 123, 059902 (2019)].
  * (53) Z. Capatti, V. Hirschi, D. Kermanschah, and B. Ruijl, Loop-Tree Duality for Multiloop Numerical Integration, Phys. Rev. Lett. 123 (2019) 151602 [arXiv:1906.06138] [InSPIRE].
  * (54) J. J. Aguilera-Verdugo, F. Driencourt-Mangin, R. J. Hernández-Pinto, J. Plenter, S. Ramirez-Uribe, A. E. Renteria Olivo, G. Rodrigo, G. F. R. Sborlini, W. J. Torres Bobadilla, and S. Tracz, Open Loop Amplitudes and Causality to All Orders and Powers from the Loop-Tree Duality, Phys. Rev. Lett. 124 (2020) 211602 [arXiv:2001.03564] [InSPIRE].
  * (55) Q. Song and A. Freitas, On the evaluation of two-loop electroweak box diagrams for $e^{+}e^{-}\to HZ$ production, JHEP 04 (2021) 179 [arXiv:2101.00308] [InSPIRE].
  * (56) I. Dubovyk, A. Freitas, J. Gluza, K. Grzanka, M. Hidding, and J. Usovitsch, Evaluation of multi-loop multi-scale Feynman integrals for precision physics, [arXiv:2201.02576] [InSPIRE].
  * (57) X. Liu, Y.-Q. Ma, and C.-Y. Wang, A Systematic and Efficient Method to Compute Multi-loop Master Integrals, Phys. Lett. B779 (2018) 353–357 [arXiv:1711.09572] [InSPIRE].
  * (58) X. Liu, Y.-Q. Ma, W. Tao, and P. Zhang, Calculation of Feynman loop integration and phase-space integration via auxiliary mass flow, Chin. Phys. C 45 (2021) 013115 [arXiv:2009.07987] [InSPIRE].
  * (59) X. Liu and Y.-Q. Ma, Multiloop corrections for collider processes using auxiliary mass flow, [arXiv:2107.01864] [InSPIRE].
  * (60) X. Liu and Y.-Q. Ma, AMFlow: a Mathematica Package for Feynman integrals computation via Auxiliary Mass Flow, (2022) .
  * (61) P. Zhang, C.-Y. Wang, X. Liu, Y.-Q. Ma, C. Meng, and K.-T. Chao, Semi-analytical calculation of gluon fragmentation into1S${}_{0}^{[1,8]}$ quarkonia at next-to-leading order, JHEP 04 (2019) 116 [arXiv:1810.07656] [InSPIRE].
  * (62) P. Zhang, C. Meng, Y.-Q. Ma, and K.-T. Chao, Gluon fragmentation into ${}^{3}{P}_{J}^{\left[1,8\right]}$ quark pair and test of NRQCD factorization at two-loop level, JHEP 08 (2021) 111 [arXiv:2011.04905] [InSPIRE].
  * (63) Q. Yang, R.-Y. Zhang, M.-M. Long, S.-M. Wang, W.-G. Ma, J.-W. Zhu, and Y. Jiang, QCD Corrections to $e^{+}e^{-}\rightarrow H^{\pm}W^{\mp}$ in Type-I THDM at Electron Positron Colliders, Chin. Phys. C 44 (2020) 093101 [arXiv:2005.11010] [InSPIRE].
  * (64) C. Brønnum-Hansen and C.-Y. Wang, Contribution of third generation quarks to two-loop helicity amplitudes for W boson pair production in gluon fusion, JHEP 01 (2021) 170 [arXiv:2009.03742] [InSPIRE].
  * (65) C. Brønnum-Hansen and C.-Y. Wang, Top quark contribution to two-loop helicity amplitudes for $Z$ boson pair production in gluon fusion, JHEP 05 (2021) 244 [arXiv:2101.12095] [InSPIRE].
  * (66) R.-H. Wu, Y.-S. Zuo, C. Meng, Y.-Q. Ma, and K.-T. Chao, NLO effects for $\Omega$ QQQ baryons in QCD Sum Rules, Chin. Phys. C 45 (2021) 093103 [arXiv:2104.07384] [InSPIRE].
  * (67) C. Brønnum-Hansen, K. Melnikov, J. Quarroz, and C.-Y. Wang, On non-factorisable contributions to t-channel single-top production, JHEP 11 (2021) 130 [arXiv:2108.09222] [InSPIRE].
  * (68) D. Baranowski, M. Delto, K. Melnikov, and C.-Y. Wang, On phase-space integrals with Heaviside functions, [arXiv:2111.13594] [InSPIRE].
  * (69) M. Beneke and V. A. Smirnov, Asymptotic expansion of Feynman integrals near threshold, Nucl.Phys. B522 (1998) 321–344 [hep-ph/9711391] [InSPIRE].
  * (70) V. A. Smirnov, Problems of the strategy of regions, Phys. Lett. B 465 (1999) 226–234 [hep-ph/9907471] [InSPIRE].
  * (71) Z.-F. Liu and Y.-Q. Ma, Feynman integrals are completely determined by linear algebra, (2022) .
  * (72) X. Ji, Parton Physics on a Euclidean Lattice, Phys.Rev.Lett. 110 (2013) 262002 [arXiv:1305.1539] [InSPIRE].
  * (73) X. Ji, Parton Physics from Large-Momentum Effective Field Theory, Sci.China Phys.Mech.Astron. 57 (2014) 1407–1412 [arXiv:1404.6680] [InSPIRE].
  * (74) A. Radyushkin, Quasi-parton distribution functions, momentum distributions, and pseudo-parton distribution functions, Phys. Rev. D 96 (2017) 034025 [arXiv:1705.01488] [InSPIRE].
  * (75) Y.-Q. Ma and J.-W. Qiu, Extracting Parton Distribution Functions from Lattice QCD Calculations, Phys. Rev. D 98 (2018) 074021 [arXiv:1404.6860] [InSPIRE].
  * (76) Y.-Q. Ma and J.-W. Qiu, Exploring Partonic Structure of Hadrons Using ab initio Lattice QCD Calculations, Phys. Rev. Lett. 120 (2018) 022003 [arXiv:1709.03018] [InSPIRE].
  * (77) O. Tarasov, Connection between Feynman integrals having different values of the space-time dimension, Phys.Rev. D54 (1996) 6479–6490 [hep-th/9606018] [InSPIRE].
  * (78) O. V. Tarasov, Application and explicit solution of recurrence relations with respect to space-time dimension, Nucl. Phys. B Proc. Suppl. 89 (2000) 237–245 [hep-ph/0102271] [InSPIRE].
  * (79) V. M. Braun, A. Vladimirov, and J.-H. Zhang, Power corrections and renormalons in parton quasidistributions, Phys. Rev. D 99 (2019) 014013 [arXiv:1810.00048] [InSPIRE].
  * (80) Z.-Y. Li, Y.-Q. Ma, and J.-W. Qiu, Extraction of Next-to-Next-to-Leading-Order Parton Distribution Functions from Lattice QCD Calculations, Phys. Rev. Lett. 126 (2021) 072001 [arXiv:2006.12370] [InSPIRE].
  * (81) T. Bitoun, C. Bogner, R. P. Klausen, and E. Panzer, Feynman integral relations from parametric annihilators, Lett. Math. Phys. 109 (2019) 497–564 [arXiv:1712.09215] [InSPIRE].


# A Simple Architecture for Enterprise Large Language Model Applications based
on Role based security and Clearance Levels using Retrieval-Augmented
Generation or Mixture of Experts

Atilla Özgür1 Atilla Özgür is partially supported by Grant DFG UY-56/5-1
Yılmaz Uygun 1

( 1School of Business, Social and Decision Sciences, Constructor University
Bremen, Bremen, Germany
{aoezguer<EMAIL_ADDRESS>
###### Abstract

This study proposes a simple architecture for Enterprise application for Large
Language Models (LLMs) for role based security and NATO clearance levels. Our
proposal aims to address the limitations of current LLMs in handling security
and information access. The proposed architecture could be used while
utilizing Retrieval-Augmented Generation (RAG) and fine tuning of Mixture of
experts models (MoE). It could be used only with RAG, or only with MoE or with
both of them. Using roles and security clearance level of the user, documents
in RAG and experts in MoE are filtered. This way information leakage is
prevented.

Keywords: Large language models, Mixture of Experts, role-based security,
clearance levels

##  1 Introduction

According to World Economic Forum [9], Venture Capital investments in the area
of artificial intelligence are about $290 billion in the last 5 years. With
the introduction of ChatGPT by OpenAI [21], interest in the large language
models (LLMs) has exploded.

Even though ChatGPT provides an application programming interface (API) to the
developers, it is a closed model. Developers could only utilize provided API.
To answer ChatGPT from OpenAI, large language model called LLaMA [27] is
introduced by Meta (Facebook). With the advent of LLaMA, the era of open
source Large Language Models started. Open source developers trained their own
models in different domains using LLaMA and published these models. Due to
success of this, other companies followed and published their own models in
open source fashion. For example, Microsoft published Phi [1] and Google
published gemini [10] open source LLMs. A more complete list could be found in
[30].

Additionally, a lot of different client libraries for these open source models
are introduced. With these client libraries, it is easier than ever to produce
applications which utilizes LLMs in their portfolio. Security of the these
applications are very important. But in academia, instead of general security,
most of the time security for LLM itself are the focus. Security of LLMs are
reviewed in a lot of different articles [31, 3, 33, 5, 18, 13, 8, 11, 28, 19,
29, 22, 7, 14, 17]. But these articles are mostly about the attack on LLMs
itself. Even articles about privacy is mostly concerned from general usage of
LLMs. According to the best of our knowledge, how to secure LLMs in
enterprises from the view point of security clearance levels and role based
security is not investigated.

A recent lawsuit by New York Times [12] shows that is is possible to reproduce
exact trained documents using appropriate prompts. Normally, this is only a
copyright and privacy problem but in the context of military or enterprise
applications, such a problem becomes a big security problem. Think about the
following example: A custom LLM application developed for NATO military
documents. If an user of this NATO LLM application could reproduce a Secret
NATO document when his security clearance is only for Confidential documents,
this reproduction of document would be a very big security problem.

This study proposes a simple role based security architecture for custom LLM
applications. This approach could be used while utilizing Retrieval-Augmented
Generation (RAG) and fine tuning of Mixture of experts models (MoE). It could
be used only with RAG, or only with MoE or with both of them. While using RAG,
only documents which user’s roles have access to will be returned from RAG.
While using MoE, only experts which user’s roles have access will be
consulted. This way, information leakage is prevented and security of the
application will be ensured.

##  2 Background Information

To be able to understand proposed model more easily, some background
information is needed. In this section, this background information about NATO
clearance levels, role based security, introduction to large language models,
Retrieval-Augmented Generation (RAG) and Mixture of experts are given.

###  2.1 NATO Clearance Levels

The North Atlantic Treaty Organization (NATO) is a military alliance. NATO
consists of 32 member states and is established after World War II. NATO like
many international organizations deals with sensitive information. This
information can range from battle plans to diplomatic communications and
intelligence reports. Leaking this information could have serious consequences
[4].

NATO has following four security classifications:

  1. 1.

Cosmic top secret

  2. 2.

Secret

  3. 3.

Confidential

  4. 4.

Restricted

Even though it is not counted among the classifications, not classified
category is also used to show that a document or an information could be
shared outside NATO but their rights belong to NATO [24].

Clearance levels or security classifications are used for following purposes:
Protect Classified Information, Minimize Risks and Ensure Trust Between Allies
Classified information is categorized according to its sensitivity. For
example size information of F-35 fighter jet could be not classified while
radar sensitive painting information of F-35 could be Cosmic top secret. In
short, NATO clearance levels are a security categorization designed to prevent
sensitive information.

###  2.2 Role-Based Access Control (RBAC)

Role-Based Access Control (RBAC) is a security approach that manages access to
resources within a system. Idea of RBAC has been around since the start of
multi user computers [26]. RBAC resolves around permissions, roles and user.

Figure 1: Example Role Based Access

Permissions control access to resources or abilities of users. An example for
permission could be log on to system. Roles are used as containers for
permissions. For example a human resources (HR) role could be used to contain
permission related to personnel management tasks. In most systems, a role
could also contain other roles. Figure 1 shows that Normal User role belongs
to other roles. Similar to Roles, Users also could belong to more than one
role.

When security clearance levels are used, most of the time, these clearance
levels names are extended together with normal roles. For example, Operator
role will be extended as Operator not classified, Operator restricted and so
on.

###  2.3 Introduction to large language models

Large Language Models are supervised learning models. That is they are trained
with using known input and outputs. LLMs are trained using a lot of text from
internet to repeatedly predict next words using previous words. See Table 1
for an example sentence: My favorite food is a Döner with spicy sauce.

Table 1: LLM Text Generation Example Input A |  Output
---|---
My favorite food is a |  Döner
My favorite food is a Döner |  with
My favorite food is a Döner with |  spicy
My favorite food is a Döner with spicy |  sauce

####  2.3.1 Training of large language models

For training LLMs, terabytes of text from internet are used. Basic workflow
for training LLM are something like below:

  1. 1.

Download 10+TB of text

  2. 2.

Get a cluster of 6k+ GPUs

  3. 3.

Train your neural network, pay $2M, wait for 12 days

  4. 4.

Obtain base model

This training could be seen as something like compressing the training
dataset. Since LLMs could only produce values that are in their training data,
for specialized tasks fine tuning of LLMs should be done. Fine tuning should
be done with quality data specially prepared for the task. Basic workflow for
fine tuning LLM is something like below:

  1. 1.

create curated dataset of quality instructions

  2. 2.

fine tune base model on this data, wait 1 day

  3. 3.

Obtain assistant model

  4. 4.

Test your model if necessary go to step 1

  5. 5.

Deploy on your servers

  6. 6.

Monitor, go to step 1

Since fine tuning uses less data, it costs less and could be repeated more.
When using Mixture of experts, first fine tuning would be slower since we will
need to fine tune multiple expert LLMs. But subsequent fine tuning for experts
will be faster since whenever new documents come only the relevant experts
will be re-fine tuned.

###  2.4 Retrieval-Augmented Generation (RAG)

As we have explained in section 2.3, LLMs are trained with massive amounts of
data. As all machine learning models, LLMs are also depended on statistical
patterns in their training data. For example, an LLM model, which is trained
in 2023, will not be able to answer questions about Euro 2024 (European
Football Championship).

Retrieval-Augmented Generation (RAG) introduced by Facebook researchers[16]
address these limitation by connecting LLMs to update data sources. These
sources could be news articles, company internal knowledge base or databases
like wikipedia.

Figure 2: Retrieval-Augmented Generation (RAG) workflow

RAG workflow could be seen in Figure 2

When a new prompt comes to LLM system, similar documents to this new prompt
are searched in databases. Most of the time a vector database is used for fast
response times. Then, LLM uses this context enriched prompt to give answers.

RAG makes LLM outputs more reliable using factual databases. Like previous
example of Euro 2024, RAG enables LLMs to use latest information if their
training data is older. RAG could also be adapted to specific domains using
relevant databases.

In our proposed role-based access control architecture, every document in RAG
databases will have allowed Roles information. When similar documents are
searched for RAG prompts, only documents which asking user have access to will
be searched. This way, information leakage will be prevented.

###  2.5 Mixture of experts

Jacobs et al [15] introduced mixture-of-experts (MoE) model in 1991. The MoE
models utilizes individual neural networks as experts instead of a single
neural network. The idea of mixture of experts are very similar to ensemble
learning [23]. Ensemble learning combines multiple classifiers in different
fashion. According to Polikar [23], among the first examples of Ensemble
Learning [6] was in 1979 by Dasarathy and Sheela. Mixture of experts name is
more popular in deep learning literature while ensemble learning is more
popular in general machine learning literature.

In mixture of experts, different neural networks to solve the problem are
called experts. Mixture of expert model uses a router or coordinator neural
network to decide which experts to utilize [20]. Router assigns a ”gating
score” to each network to indicate how relevant the input to the expert is. A
softmax function is used to transform the gating scores to a probability
distribution. When doing the predictions only experts with highest gating
scores are activated. This router idea is very similar to voting in the
ensemble classifiers. Voting classifiers are used in a lot of different
domains like intrusion detection [32, 2, 25, 35, 34]

In our role based access control architecture, gating scores of some experts
will be zero according to the roles of the caller. For example, a normal user
is calling the large language model with a human resource prompt. Only the
experts which normal user has access to will be activated. Other experts like
human resources expert will not be called at all.

##  3 Proposed Architecture based on Role based security and Clearance Levels

###  3.1 User, Role, Security Clearance to Documents Mapping

Most enterprises map their users to one or more roles. This mapping could be
stored in directory services like active directory (Microsoft) or Enterprise
Resource planning databases. For our LLM application, this user to role and
user to security clearance levels mappings should be accessible from its
programming interface. Additionally, role to documents and security clearance
levels to documents should also be accessible. This mapping information could
be stored in LLM application’s own database or a web service could be provided
to the application. Basic mapping could be seen in Figure 3.

Figure 3: Entity Relationship Diagram for User to Role and Clearance Level to
Document Mapping

Every User has more than one role, every User has exactly one security
clearance. Every Document has exactly one security clearance. Every Role has
access to zero or more documents. It is implied that if a user does not have
access to necessary clearance level, he will not be able to access the
document.

###  3.2 Training LLM using Roles and Security Clearance Levels

Proposed architecture presumes using local open source LLMs due to security
requirements. But if strict local security is not required then commercial
LLMs like OpenAI ChatGPT with Retrieval-Augmented Generation (RAG) could be
used. When commercial LLMs are used, training step may not be necessary
according to requirements or only not confidential documents could be used for
training.

For local open source LLMs, Mixture of Experts model should be used. For this
case, it is advised to train multiple experts for every role multiplied by
five, that is four clearance levels plus not classified. For example, let’s
assume that the LLM application has four roles: HR, Accounting, Normal User,
IT as shown in Figure 1. Then, we should train $4*5=20$ experts. This experts
will be named like HR Not Classified, HR Restricted, HR Confidential, HR
Secret, HR Cosmic top secret. See full names in Table 2.

Table 2: Roles and Clearance Levels for Experts |  HR |  Accounting |  Normal User |  IT
---|---|---|---|---
Not classified |  HR Not classified |  Accounting Not classified |  Normal User Not classified |  IT Not classified
Restricted |  HR Restricted |  Accounting Restricted |  Restricted |  IT Restricted
Confidential |  HR Confidential |  Accounting Confidential |  Normal User Confidential |  IT Confidential
Secret |  HR Secret |  Accounting Secret |  Normal User Secret |  IT Secret
Cosmic top secret |  HR Cosmic top secret |  Accounting Cosmic top secret |  Normal User Cosmic top secret |  IT Cosmic top secret

###  3.3 Inference using Retrieval-Augmented Generation (RAG)

While doing the inference using Retrieval-Augmented Generation (RAG), role
information and security clearance levels should be used. Most of the vector
databases allow use of filters. Using filters, only documents a user has
access to should be returned.

Figure 4:  Sequence Diagram for Role/Clearance level based access to LLM only
using RAG

As could be seen in Figure 4, this architecture does not presume mixture of
experts or local open source model; therefore, this role based model could
also be used with commercial LLMs like ChatGPT.

###  3.4 Inference using Mixture of experts models (MoE)

In this inference model, local open source LLM should be used. Experts in MoE
are trained using only documents for which role and security clearance level
has an access to. For performance reasons, filtering in the router part of MoE
would be useful. If router does not ask answers from experts for whom user has
no security clearance or necessary role, the application will perform faster.
Full workflow could be seen in Figure 5.

Figure 5:  Sequence Diagram for Role/Clearance level based access to LLM only
using MoE

###  3.5 Inference using both RAG and MoE

Doing inference using both Retrieval-Augmented Generation (RAG) and Mixture of
experts models (MoE) could be seen in Figure 6. Basically, this workflow is
combination of previous two workflows.

Figure 6: Role based access to LLM Sequence Diagram

##  4 Conclusion

A simple role based security architecture for Large Language Model (LLM)
applications are proposed. Background information necessary for the
understanding of the architecture given in multiple sections. Proposed
approach is usable with Retrieval-Augmented Generation (RAG) and/or Mixture of
experts models (MoE). Normally usage of local open source LLM models are
assumed but RAG version is also suitable for commercial LLMs. Sequence
diagrams for workflows using RAG, MoE and RAG+MoE are given.

## References

  * [1] Marah Abdin and et al.  Phi-3 technical report: A highly capable language model locally on your phone, 2024.
  * [2] Muhammad Ali, Mansoor-ul Haque, Muhammad Hanif Durad, Anila Usman, Syed Muhammad Mohsin, Hana Mujlid, and Carsten Maple.  Effective network intrusion detection using stacking-based ensemble approach.  International Journal of Information Security, 22(6):1781–1798, July 2023.
  * [3] Arijit Ghosh Chowdhury, Md Mofijul Islam, Vaibhav Kumar, Faysal Hossain Shezan, Vaibhav Kumar, Vinija Jain, and Aman Chadha.  Breaking down the defenses: A comparative survey of attacks on large language models, 2024.
  * [4] NATO Security Comittee.  A short guide to the handling of nato classified information.  Technical report, NATO, 1958.
  * [5] Badhan Chandra Das, M. Hadi Amini, and Yanzhao Wu.  Security and privacy challenges of large language models: A survey, 2024\.
  * [6] B.V. Dasarathy and B.V. Sheela.  A composite classifier system design: Concepts and methodology.  Proceedings of the IEEE, 67(5):708–713, 1979.
  * [7] Erik Derner, Kristina Batistič, Jan Zahálka, and Robert Babuška.  A security risk taxonomy for large language models, 2023.
  * [8] Mohamed Amine Ferrag, Fatima Alwahedi, Ammar Battah, Bilel Cherif, Abdechakour Mechri, and Norbert Tihanyi.  Generative ai and large language models for cyber security: All insights you need, 2024.
  * [9] World Economic Forum.  How venture capital is investing in ai in the top five global economies — and shaping the ai ecosystem, May 2024.
  * [10] Gemma Team and et al.  Gemma: Open models based on gemini research and technology, 2024.
  * [11] Takeshi Goto, Kensuke Ono, and Akira Morita.  Evaluating the cybersecurity robustness of commercial llms against adversarial prompts: A promptbench analysis.  TechRxiv, March 2024.
  * [12] Michael M. Grynbaum and Ryan Mac.  The times sues openai and microsoft over a.i. use of copyrighted work, 2023.
  * [13] Mohammed Hassanin and Nour Moustafa.  A comprehensive overview of large language models (llms) for cyber defences: Opportunities and directions, 2024.
  * [14] Umar Iqbal, Tadayoshi Kohno, and Franziska Roesner.  Llm platform security: Applying a systematic evaluation framework to openai’s chatgpt plugins, 2023.
  * [15] Robert A. Jacobs, Michael I. Jordan, Steven J. Nowlan, and Geoffrey E. Hinton.  Adaptive mixtures of local experts.  Neural Computation, 3(1):79–87, February 1991.
  * [16] Patrick Lewis, Ethan Perez, Aleksandra Piktus, Fabio Petroni, Vladimir Karpukhin, Naman Goyal, Heinrich Küttler, Mike Lewis, Wen-tau Yih, Tim Rocktäschel, Sebastian Riedel, and Douwe Kiela.  Retrieval-augmented generation for knowledge-intensive nlp tasks, 2020\.
  * [17] Jiaju Lin, Haoran Zhao, Aochi Zhang, Yiting Wu, Huqiuyue Ping, and Qin Chen.  Agentsims: An open-source sandbox for large language model evaluation, 2023.
  * [18] Zefang Liu.  A review of advancements and applications of pre-trained language models in cybersecurity.  In 2024 12th International Symposium on Digital Forensics and Security (ISDFS). IEEE, April 2024.
  * [19] Sichun Luo, Wei Shao, Yuxuan Yao, Jian Xu, Mingyang Liu, Qintong Li, Bowei He, Maolin Wang, Guanzhi Deng, Hanxu Hou, Xinyi Zhang, and Linqi Song.  Privacy in llm-based recommendation: Recent advances and future directions, 2024.
  * [20] Saeed Masoudnia and Reza Ebrahimpour.  Mixture of experts: a literature survey.  Artificial Intelligence Review, 42(2):275–293, May 2012.
  * [21] OpenAI.  Chatgpt (mar 14 version) [large language model], 2022.
  * [22] Rahul Pankajakshan, Sumitra Biswal, Yuvaraj Govindarajulu, and Gilad Gressel.  Mapping llm security landscapes: A comprehensive stakeholder risk assessment proposal, 2024.
  * [23] R. Polikar.  Ensemble learning.  Scholarpedia, 4(1):2776, 2009.  revision #186077.
  * [24] Alasdair Roberts.  Nato, secrecy, and the right to information.  East EuropeanConstitutional Review, 12 2002.
  * [25] Yakub Kayode Saheed and Sanjay Misra.  A voting gray wolf optimizer-based ensemble learning models for intrusion detection in the internet of things.  International Journal of Information Security, 23(3):1557–1581, January 2024.
  * [26] Ravi S. Sandhu.  Role-based Access Control, pages 237–286.  Elsevier, 1998.
  * [27] Hugo Touvron, Thibaut Lavril, Gautier Izacard, Xavier Martinet, Marie-Anne Lachaux, Timothée Lacroix, Baptiste Rozière, Naman Goyal, Eric Hambro, Faisal Azhar, Aurelien Rodriguez, Armand Joulin, Edouard Grave, and Guillaume Lample.  Llama: Open and efficient foundation language models, 2023.
  * [28] Shang Wang, Tianqing Zhu, Bo Liu, Ming Ding, Xu Guo, Dayong Ye, Wanlei Zhou, and Philip S. Yu.  Unique security and privacy threats of large language model: A comprehensive survey, 2024.
  * [29] Fangzhou Wu, Ning Zhang, Somesh Jha, Patrick McDaniel, and Chaowei Xiao.  A new era in llm security: Exploring security concerns in real-world llm-based systems, 2024.
  * [30] Eugene Yan.  Open large language models list, 2024.
  * [31] Yifan Yao, Jinhao Duan, Kaidi Xu, Yuanfang Cai, Zhibo Sun, and Yue Zhang.  A survey on large language model (llm) security and privacy: The good, the bad, and the ugly.  High-Confidence Computing, 4(2):100211, June 2024.
  * [32] Hayette Zeghida, Mehdi Boulaiche, and Ramdane Chikh.  Securing mqtt protocol for iot environment using ids based on ensemble learning.  International Journal of Information Security, 22(4):1075–1086, March 2023.
  * [33] Jie Zhang, Haoyu Bu, Hui Wen, Yu Chen, Lun Li, and Hongsong Zhu.  When llms meet cybersecurity: A systematic literature review, 2024.
  * [34] Atilla Özgür and Hamit Erdem.  Feature selection and multiple classifier fusion using genetic algorithms in intrusion detection systems.  Journal of the Faculty of Engineering and Architecture of Gazi University, 33:0 – 0, 2018.
  * [35] Atilla Özgür, Fatih Nar, and Hamit Erdem.  Sparsity-driven weighted ensemble classifier.  International Journal of Computational Intelligence Systems, 11:962–978, 2018.