Differential equations, recurrence relations, and quadratic constraints for $L$-loop two-point massive tadpoles and propagators

TL;DR

This paper develops a comprehensive differential equation and recurrence relation framework for L-loop two-point massive tadpole integrals, expressing solutions via Lauricella functions and deriving quadratic constraints, advancing analytical methods in multiloop Feynman integrals.

Contribution

The paper introduces a Pfaffian differential system for L-loop tadpole integrals and their relation to Lauricella functions, providing new tools for multiloop integral analysis.

Findings

Derived differential equations with Pfaffian form for L-loop integrals.

Expressed solutions in terms of Lauricella functions $F_C^{(L)}$.

Obtained quadratic constraints for solutions near integer dimensions.

Abstract

We consider $L$-loop two-point tadpole (watermelon) integral with arbitrary masses, regularized both dimensionally and analytically. We derive differential equation system and recurrence relations (shifts of dimension and denominator powers). Since the $L$-loop sunrise integral corresponds to the $(L+1)$-loop watermelon integral with one cut line, our results are equally applicable to the former. The obtained differential system has a Pfaffian form and is linear in dimension and analytic regularization parameters. In general case, the solutions of this system can be expressed in terms of the Lauricella functions $F_C^{(L)}$ with generic parameters. Therefore, as a by-product, we obtain, to our knowledge for the first time, the Pfaffian system for $F_C^{(L)}$ for arbitrary $L$. The obtained system has no apparent singularities. Near odd dimension and integer denominator powers the system can be easily transformed into canonical form. Using the symmetry properties of the matrix in the right-hand side of the differential system, we obtain quadratic constraints for the expansion of solutions near integer dimension and denominator powers. In particular, we obtain quadratic constraints for Bessel moments similar to those discovered by Broadhurst and Roberts.

Equations

$$\mathcal{T}(m_0,\dots ,m_L)=2^{D-1}i{\pi }^{D/2}\Gamma \left(\frac{D}{2}\right)\int \delta \left({\sum }_{l=0}^L{p}_l\right)\prod _{l=0}^L\left[\frac{d^D{p}_l}{i{\pi }^{D/2}}\frac{m_l\Gamma \left({\alpha }_l\right)}{{\left(m_l^2-p_l^2-i0\right)}^{{\alpha }_l}}\right].$$

$$\mathcal{S}(m_0,\dots ,m_L|q^2)=2^{D-1}i{\pi }^{D/2}\Gamma \left(\frac{D}{2}\right)\int \delta \left(q+{\sum }_{l=0}^L{p}_k\right)\prod _{l=0}^L\left[\frac{d^D{p}_l}{i{\pi }^{D/2}}\frac{m_k\Gamma \left({\alpha }_l\right)}{{\left(m_l^2-p_l^2-i0\right)}^{{\alpha }_l}}\right].$$

$$\begin{gathered} (m_0{)}^{D-1}\mathcal{S}(m_1,\dots ,m_L|m_0^2)=2^{D-2}i{\pi }^{D/2-1}{\Gamma }^2\left(\frac{D}{2}\right) \\ \times \int \delta \left({\sum }_{l=0}^{L+1}{p}_k\right)\frac{d^D{p}_0}{i{\pi }^{D/2}}\left[\frac{m_0}{m_0^2-p_0^2-i0}-\frac{m_0}{m_0^2-p_0^2+i0}\right] \\ \times \prod _{l=1}^L\left[\frac{d^D{p}_l}{i{\pi }^{D/2}}\frac{m_k\Gamma \left({\alpha }_l\right)}{{\left(m_l^2-p_l^2-i0\right)}^{{\alpha }_l}}\right], \end{gathered}$$

$$\mathcal{T}=2^{D-1}{\pi }^{D/2}\Gamma \left(\frac{D}{2}\right)\int \delta \left({\sum }_{l=0}^L{p}_k\right)\prod _{l=0}^L\left[\frac{d^D{p}_l}{{\pi }^{D/2}}\frac{m_k\Gamma \left({\alpha }_l\right)}{{\left(p_l^2+m_l^2\right)}^{{\alpha }_l}}\right],$$

$$\mathcal{T}=\underset{0}{\overset{\infty }{\int }}dx{x}^{D-1}\prod _{l=0}^L\int \frac{d^{D}p}{{\pi }^{D/2}}\frac{m_k\Gamma ({\alpha }_k)e^{ipx}}{(p^2+m_k^2{)}^{{\alpha }_k}}=\underset{0}{\overset{\infty }{\int }}dx{x}^{D-1}\prod _{l=0}^L{P}_0({\mu }_k,m_k,x).$$

$$\begin{gathered} P_0(\mu ,m,x)=m\Gamma (\alpha )\int \frac{d^{D}p}{{\pi }^{D/2}}\frac{e^{ipx}}{(p^2+m^2{)}^{\alpha }}=m\underset{0}{\overset{\infty }{\int }}d\lambda {\lambda }^{\alpha -1}\int \frac{d^{D}p}{{\pi }^{D/2}}{e}^{ipx-\lambda \left(p^2+m^2\right)} \\ =m\underset{0}{\overset{\infty }{\int }}d\lambda {\lambda }^{\frac{1-\mu }{2}-1}{e}^{-\lambda m^2-\frac{x^2}{4\lambda }}=2m{\left(\frac{x}{2m}\right)}^{\frac{1-\mu }{2}}{K}_{\frac{\mu -1}{2}}\left(mx\right). \end{gathered}$$

$$T_{\mathbf{a}}\left(D,\mathbf{\mu },\mathbf{m}\right)=\underset{0}{\overset{\infty }{\int }}dx{x}^{D-1}{t}_{\mathbf{a}},$$

$$t_{\mathbf{a}}=\prod _{l=0}^L{P}_{a_l}({\mu }_l,m_l,x),$$

$$\begin{gathered} P_1(\mu ,m,x)=\frac{1}{m}{\partial }_x{P}_0(\mu ,m,x)=-\frac{x}{2}\underset{0}{\overset{\infty }{\int }}d\lambda {\lambda }^{-\frac{1+\mu }{2}-1}{e}^{-\lambda m^2-\frac{x^2}{4\lambda }} \\ =-2m{\left(\frac{x}{2m}\right)}^{\frac{1-\mu }{2}}{K}_{\frac{\mu +1}{2}}\left(mx\right). \end{gathered}$$

$$\begin{gathered} T_{\mathbf{a}}(D,\mathbf{\mu })=\prod _l{\left(\frac{-1}{2m_l}\right)}^{a_l}\mathcal{T}(D\to D+|\mathbf{a}|,{\mu }_l\to {\mu }_l+2a_l) \\ =\prod _l{\left(\frac{-1}{2m_l}\right)}^{a_l}\mathcal{T}(D\to D+|\mathbf{a}|,{\alpha }_l\to {\alpha }_l+|\mathbf{a}|/2-a_l), \end{gathered}$$

$$m{\partial }_m{P}_0$$

$$m{\partial }_m{P}_1$$

$$m{\partial }_{m}P=(mx\sigma +\mu \bar{n})P,$$

$$x{\partial }_{x}P=(m{\partial }_m-\mu )P=(mx\sigma -\mu n)P,$$

$$\underset{0}{\overset{\infty }{\int }}dx{x}^{D-1}\left(xM-W\right)t=0.$$

$$M=\sum _{l=0}^L{m}_l{\sigma }_l,W=\sum _{l=0}^L{\mu }_l{n}_l-D,$$

$$\underset{0}{\overset{\infty }{\int }}dx{x}^{D}t=M^{-1}WT$$

$$\sum _l{\eta }_l{m}_l=0$$

$${\partial }_{m_l}T=A_{l}T,A_l={\sigma }_l{M}^{-1}W+\frac{{\mu }_l{\bar{n}}_l}{m_l}.$$

$$dT=AT,$$

$$A=\sum _{l=0}^L{A}_{l}d{m}_l=dM{M}^{-1}W+\sum _{l=0}^L{\mu }_l{\bar{n}}_l\frac{d{m}_l}{m_l}.$$

$$\mathcal{P}={\left(-1\right)}^{{\sum }_i{n}_i}=\prod _i(1-2n_i).$$

$$P_0\left(\mu -2,m,x\right)$$

$$P_1\left(\mu -2,m,x\right)$$

$$T(\mathbf{\mu }-2{\mathbf{e}}_l)=R_l\left(D,\mathbf{\mu }\right)T(\mathbf{\mu }),$$

$$R_l\left(D,\mathbf{\mu }\right)=-\frac{1}{2m_l}\left(A_l-\frac{{\bar{n}}_l}{m_l}\right)=-\frac{1}{2m_l}\left({\sigma }_l{M}^{-1}W+\frac{\left({\mu }_l-1\right){\bar{n}}_l}{m_l}\right)$$

$$R(D,\mathbf{\mu })=M^{-1}W,$$

$$T(D+1,\mathbf{\mu })=R(D,\mathbf{\mu })T(D,\mathbf{\mu }).$$

$$T\left(D-2,\mathbf{\mu }-2\mathbf{e}\right)=\mathfrak{R}(D,\mathbf{\mu })T(D,\mathbf{\mu }),$$

$$\mathfrak{R}(D,\mathbf{\mu })=R^{-1}(D-2,\mathbf{\mu }-2\mathbf{e})R^{-1}(D-1,\mathbf{\mu }-2\mathbf{e})\prod _{l=0}^L{R}_l\left(D,\mathbf{\mu }-2{\sum }_{j=l+1}^L{\mathbf{e}}_j\right)$$

$$K_{\frac{{\mu }_l\pm 1}{2}}\left(m_{l}x\right)\to \mp \frac{{\sigma }_l\pi }{2\cos \frac{\pi {\mu }_l}{2}}{I}_{{\sigma }_l\frac{{\mu }_l\pm 1}{2}}\left(m_{l}x\right),$$

$$V^{(\mathbf{\rho })}=\underset{0}{\overset{\infty }{\int }}dx{x}^{D-1}P({\mu }_0,m_0,x)\otimes \underset{l=1}{\overset{L}{⨂}}{Q}^{({\rho }_l)}({\mu }_l,m_l,x),$$

Full text

1 Introduction

Multiloop massive sunrise and watermelon integrals (see Fig. 1) are ubiquitous in various physical applications. Being seemingly simple, they, in fact, bring with themselves a lot of complexity in any calculation. The reason is that their $\epsilon$-expansion can not be expressed in terms of the standard functional basis, viz., multiple polylogarithms Goncharov1998 ; Weinzierl:2002hv . Already starting from two loops one observes complicated iterated integrals involving modular forms and/or elliptic curves Adams2017b ; Broedel2018c . For $L$ loops, the $\epsilon$ expansion of massive sunrise integrals is likely to involve periods of certain $(L-1)$-dimensional Calabi-Yau varieties (cf. Refs. arXiv0407327v1.math ; vanhove2019feynman ).

For the general case of different masses the closed expressions for these integrals are known since Ref. Berends:1993ee (see also Ref. Kalmykov:2016lxx ) in terms of Lauricella functions $F_{C}^{(L)}$ defined via multiple hypergeometric series. Although these functions have a long history of investigation, starting from lauricella1893sulle , and much is known about their properties, still many questions remain unanswered. In particular, the functions which appear in coefficients of their expansion in the Laurent series in parameters are not well investigated. Meanwhile, these coefficients are the most important from the point of view of multiloop calculations. From the purely mathematical point of view, among the four flavors of Lauricella functions $F_{A}$, $F_{B}$, $F_{C}$, and $F_{D}$, the functions $F_{C}^{(L)}$ are the only ones for which the differential system in Pfaffian form is not yet known for $L>2$.

In the present paper we derive the Pfaffian system of the differential equations for the $L$-loop watermelon and sunrise graphs with arbitrary masses regularized both dimensionally and analytically. Similar to Refs. Berends:1993ee ; Kalmykov:2016lxx , we find the fundamental system of solutions in terms of the Lauricella functions $F_{C}^{(L)}$ with generic parameters and arguments. In this way, our results give a Pfaffian system for Lauricella $F_{C}^{(L)}$ for all $L$, to the best of our knowledge, for the first time111For the case $L=2$ a Pfaffian system was found in Ref. kato1988pfaffian ..

In addition to the differential system, we also derive the operators shifting dimension or the powers of denominators. From the point of view of $F_{C}^{(L)}$ functions these operators provide a certain form of the contiguity relations. Near $D=1$ we reduce the differential system to canonical form which allows one to obtain the coefficients of expansion in dimension and analytic regulators in terms of the Goncharov’s polylogarithms.

Last but not least, we observe that the matrix in the right-hand side of the obtained differential system has special features which result in the existence of the bilinear relations for the solutions of the systems that differ in the sign of the right-hand side (cf. Ref. Goto2013 ). Combined with the shift operators, these relations allow us to obtain the infinite set of the quadratic relations for the coefficients of expansion in dimension and analytic regulators near any integer point, including the most interesting case of even $D$.

2 Differential equations in Pfaffian form

We consider the $L$-loop watermelon integral, Fig.1, with different masses and powers of denominators (exponents),

[TABLE]

Here we imply the Lorentzian signature, so that $p^{2}=p_{0}^{2}-\boldsymbol{p}^{2}$. Note that we have chosen a somewhat unconventional overall normalization factor $2^{D-1}\Gamma\left(\tfrac{D}{2}\right)\prod_{l}m_{l}\Gamma\left(\alpha_{l}\right)$ in order to simplify some expressions below.

Before moving any further let us discuss the applicability of our present consideration to the sunrise integral

[TABLE]

The obvious relation between $\mathcal{S}$ and $\mathcal{T}$ is $\mathcal{S}|_{q^{2}=0}=\mathcal{T}$. However, we can also represent the $(L-1)$-loop sunrise integral as an $L$-loop watermelon integral with one cut line:

[TABLE]

where $m_{0}=\sqrt{q^{2}}$. The two terms in square brackets on the last line correspond to two different deformations of the integration contour over $p_{0}^{0}$. As the differential system we are looking for is not aware of the integration contours, it is clear that $m_{0}^{D-1}\mathcal{S}(m_{1},\ldots,m_{L}|m_{0}^{2})$ satisfies the same differential equations as $\mathcal{T}(m_{0},\ldots,m_{L})$ with $\alpha_{0}=1$. Therefore, in the rest of the paper we will concentrate on the watermelon integral $\mathcal{T}$ only.

After the Wick rotation we can rewrite Eq. (1) as

[TABLE]

where now $p^{2}=p_{0}^{2}+\boldsymbol{p}^{2}$. Representing the function $\delta(p)$ as $\intop\frac{d^{D}x}{(2\pi)^{D}}e^{ipx}$, we obtain

[TABLE]

As we shall see immediately, the integral $\int\frac{d^{D}p}{\pi^{D/2}}\frac{m\Gamma(\alpha)e^{ipx}}{(p^{2}+m^{2})^{\alpha}}$ depends on $D$ and $\alpha$ only via the combination $\mu=D-2\alpha+1$ and we have, therefore, introduced

[TABLE]

Here $K_{\nu}(y)$ is the Macdonald function.

In order to derive the differential equations for this integral, we introduce the set of $2^{L+1}$ functions labeled by a binary number (a string of zeros and ones) $\boldsymbol{a}=a_{0}a_{1}\ldots a_{L}$,

[TABLE]

where

[TABLE]

$P_{0}$ is defined in eq. (6), and

[TABLE]

The original integral (4) is recovered as $\mathcal{T}=T_{\boldsymbol{0}}\left(D,\boldsymbol{\mu},\boldsymbol{m}\right)$. The rationale behind introducing auxiliary $2^{L+1}-1$ functions is that it turns out to be possible to derive a closed system of linear differential equations for $T_{\boldsymbol{a}}$. Note that all components of $T_{\boldsymbol{a}}$ can be expressed via $\mathcal{T}$ with shifted dimension and indices:

[TABLE]

where $|\boldsymbol{a}|=\sum_{l}a_{l}$. It is convenient to treat the $2^{L+1}$ quantities $V_{\boldsymbol{a}}$ as components of the vector $V$ in $\underbrace{\mathbb{C}^{2}\otimes\ldots\otimes\mathbb{C}^{2}}_{L+1\text{ factors}}$, where $V_{\boldsymbol{a}}$ can be either $T_{\boldsymbol{a}}$ or $t_{\boldsymbol{a}}$. In what follows, we assume that the factors in $\mathbb{C}^{2}\otimes\ldots\otimes\mathbb{C}^{2}$ are numbered starting from zero, so that, e.g., $L$-th factor is the right-most one. Besides, to lighten notation we often omit the arguments of functions.

In order to derive the differential system for $T_{\boldsymbol{a}}$ we differentiate $P_{0,1}(\mu,m,x)$ with respect to $m$ and use the differential equation for $K_{\nu}(y)$. We obtain

[TABLE]

In matrix notations we may rewrite this as

[TABLE]

where $P=\begin{pmatrix}P_{0}\\ P_{1}\end{pmatrix}$, $\sigma=\sigma_{x}=\begin{pmatrix}0&1\\ 1&0\end{pmatrix}$, and $\bar{n}=1-n=\frac{1}{2}(1+\sigma_{z})=\begin{pmatrix}1&0\\ 0&0\end{pmatrix}$.

Using these formulae, we see that the expression for the derivative of $T$, in addition to $T$ multiplied by a matrix, contains the integrals of the form $\intop_{0}^{\infty}dx\,x^{D}\,t$. Note the additional factor $x$ in the integrand as compared to the definition (7). In order to get rid of this factor, we consider the identity $\int dx\partial_{x}(x^{D}t)=0$. Explicitly differentiating and using the relation

[TABLE]

which, again, follows from the differential equation for $K_{\nu}$, we obtain

[TABLE]

Here

[TABLE]

and the operators $\sigma_{l}$ and $n_{l}$ act in $l$-th factor of $\mathbb{C}^{2}\otimes\ldots\otimes\mathbb{C}^{2}$ as $\sigma$ and $n$, respectively. Therefore,

[TABLE]

whenever $M^{-1}$ exists. Since the operators $\sigma_{l}$ commute with each other and their eigenvalues are $\eta_{l}=\pm 1$, the eigenvalues of $M$ are $\sum_{l}\eta_{l}m_{l}$. Therefore, $M$ is not invertible iff

[TABLE]

for some choice of signs $\eta_{l}$.

Thus, using Eqs. (12) and (16) we obtain the following expression for the derivative of $T$ with respect to $m_{k}$

[TABLE]

Equivalently, the above equations can be written as

[TABLE]

Remarkably, the differential form $A$ is closed ($dA=0$) which can be easily proved if one takes into account the pairwise commutativity of $\sigma_{l}$. Then the integrability condition $dA=A\wedge A$ requires that $A\wedge A=0$, see Appendix A for the derivation. The system (19) equipped with the condition $dA=0$ is said to be in Pfaffian form. We observe that all singularities of the differential form $A$ are located on the $L+1$ hyperplanes defined by $m_{l}=0$ and $2^{L}$ hyperplanes defined by Eq. (17) (cf. Ref. Goto2013 ). The latter hyperplanes correspond to vanishing of the eigenvalues of the matrix $M$.

Let us note that the differential system (19) in fact splits into two separate systems, each for $2^{L}$ functions $T_{\boldsymbol{a}}$ with $|\boldsymbol{a}|$ being either even or odd. This due to the commutativity of $A$ with the parity operator

[TABLE]

Of course, we are mostly interested in the subsystem for even components which involves our original tadpole integral $\mathcal{T}$.

3 Shifting exponents and dimension

Let us now obtain the recurrence relations in $\alpha_{k}$ and $D$. For the former we use the identities

[TABLE]

Then, using Eq. (18), we have

[TABLE]

where

[TABLE]

and $\boldsymbol{e}_{l}=(\ldots,0,\underset{l\text{-th}}{1},0,\ldots)$ is the vector with $l$-th component equal to $1$ and all other components equal to zero.

Let us now derive the operator which shifts the dimension. Note that the identity (16) defines the operator

[TABLE]

which shifts $D$ by $+1$ at fixed $\boldsymbol{\mu}$, i.e.

[TABLE]

In order to obtain the operator which shifts $D$ at fixed $\boldsymbol{\alpha}$, we should first use the operators $R_{l}$ to shift $\mu\to\mu-2$ and then shift the dimension by $-2$ at fixed $\mu$ using the operator $R^{-1}$. We obtain

[TABLE]

where

[TABLE]

is the operator which shifts $D$ by $-2$ at fixed $\boldsymbol{\alpha}$. Here and below $\boldsymbol{e}=\sum_{j=0}^{L}\boldsymbol{e}_{j}=(1,1,\ldots,1)$ and the product of matrices corresponds to the multiplication from the left to the right, so that $\prod_{l=0}^{L}R_{l}=R_{0}\ldots R_{L}$.

4 Basis of solution space

While the function $T$ as defined in Eq. (7) is a specific solution of the differential system (19) and recurrence relations (24) and (27), it is natural to consider the linear space of all solutions of these equations. In the present section we fix a basis of $2^{L+1}$ functions in this space. In order to do this, we note that the derivation of the differential equations and recurrence relations in Sections 2 and 3 can be explicitly repeated if some of the Macdonald functions $K_{\nu}$ are replaced with $\frac{I_{\pm\nu}}{\sin\pi\nu}$. More specifically, let us make the following replacement in the definition of $t_{\boldsymbol{a}}$, Eq. (8):

[TABLE]

for all $l\in S$, where $S$ is some subset of $E=\{0,1,\ldots,L\}$. Here each $\sigma_{l}$ can be either $+1$ or $-1$. Defined in this way, the function $T_{\boldsymbol{a}}$ formally satisfies the differential and difference equations (19), (24), and (27). There may be however a problem with the convergence of the integral over $x$ due to the exponential growth of the modified Bessel functions. For all masses being positive, the convergence condition is $\sum_{l\in S}m_{l}<\sum_{l\in E\backslash S}m_{l}$. In particular, $E\backslash S$ should not be empty. Therefore, we take $E\backslash S=\{0\}$, i.e. $S=\{1,\ldots,L\},\,$ and define the set of functions

[TABLE]

where $\boldsymbol{\rho}=\rho_{1}\rho_{2}\ldots\rho_{L}$ is a binary number. We assume that $m_{l}>0$ and $m_{0}>\sum_{l=1}^{L}m_{l}$, so that the integrand decays exponentially at large $x$. Using the identity

[TABLE]

we can express $T$, Eq. (7), via $V^{(\boldsymbol{\rho})}$ as

[TABLE]

The set of functions (31) contains $2^{L}$ functions numbered by the vector $\boldsymbol{\rho}$. Meanwhile, the solution space of the differential system (19) is $2^{L+1}$-dimensional. In order to obtain the whole set of solutions, we use the symmetry of the differential system related to the operator $\mathcal{P}$, Eq. (21). Therefore, we can define functions

[TABLE]

which are the solutions of the differential system. Note that we have introduced the additional factor $(-1)^{|\boldsymbol{\rho}|}$ for further convenience. In terms of the components we have

[TABLE]

so that the solutions with $\varrho_{0}=|\boldsymbol{\rho}|\ (\mathrm{mod}\ 2)$ have nonzero components with even $|\boldsymbol{a}|$.

Together with the different choices of $\boldsymbol{\rho}$, this gives us $2^{L+1}$ solutions of the differential system (19). It is also easy to see that $\mathcal{P}R_{l}=R_{l}\mathcal{P}$ and, therefore, $V^{(\varrho_{0},\boldsymbol{\rho})}$ is also a solution of the $\mu$-shifting recurrences (24). As to the the operator $R$, Eq. (27), we have $\mathcal{P}R=-R\mathcal{P}$ and therefore

[TABLE]

where we used the notation $\bar{x}=1-x$. Thus, the functions (35) satisfy the relations shifting $D$ by $2$,

[TABLE]

which follow from Eq. (27)222Note, that we can also construct $2^{L+1}$ functions which satisfy (27), e.g., by taking $V^{(\boldsymbol{\rho})}=V^{(0\boldsymbol{\rho})}+V^{(1\boldsymbol{\rho})}$ and $(V^{(0\boldsymbol{\rho})}-V^{(1\boldsymbol{\rho})})g(D)$, where $g$ is any antiperiodic function of $D$, i.e., $g(D+1)=-g(D)$..

Note that our basis functions $V^{(\varrho_{0},\boldsymbol{\rho})}$ can be expressed via Lauricella function $F_{C}^{(L)}$ (cf. Berends:1993ee ; Kalmykov:2016lxx ). This function is defined via the hypergeometric series

[TABLE]

where $(a)_{k}=\Gamma(a+k)/\Gamma(a)$ is the Pochhammer symbol. In order to do this we use the expansion of the modified Bessel function $I_{\nu}$:

[TABLE]

Substituting this expansion in Eq. (31) and taking the integral over $x$ with the help of the identity

[TABLE]

we obtain

[TABLE]

where

[TABLE]

The functions $V^{(\varrho_{0},\boldsymbol{\rho})}$ can be obtained from Eq. (42) with the help of Eq. (36).

Note that

[TABLE]

where $x_{l}=(m_{l}/m_{0})^{2}$ and the parameters

[TABLE]

are generic for generic $\boldsymbol{\mu}$ and $D$. Therefore, $V^{({1\ldots 1})}_{{0\ldots 0}}$ is proportional to $F_{C}^{(L)}$ with generic parameters. Recall that the system (19) splits into two subsystems, for even and odd components. Thus, we see that the subsystem of $2^{L}$ equations for even components gives a Pfaffian system for Lauricella $F_{C}^{(L)}$.

5 Bilinear relations and expansion near integer $D$ and $\boldsymbol{\mu}$.

The fundamental matrix of the differential system (19) can be written in the form of path-ordered exponent

[TABLE]

where ${\mathcal{C}}$ is some path in $L+1$-dimensional space of mass parameters starting from some fixed point $\boldsymbol{m}_{0}$ and ending in $\boldsymbol{m}$333Note that the matrix $A$ has singularities when either $\det M=0$ or $m_{l}=0$. Naturally, the path $\mathcal{C}$ should not intersect the singular locus of $A$, and, therefore, $U$ depends not only on the end points $\boldsymbol{m}_{0}$ and $\boldsymbol{m}$, but also on the equivalence class of $\mathcal{C}$.. Note that $A$ is linear homogeneous in $D$ and $\boldsymbol{\mu}$ which, for the expansion near $D=\boldsymbol{\mu}=0$, corresponds to a certain generalization of the canonical form of Ref. Henn2013 to the case of analytical regularization. Thanks to this property, the expansion around the point $D=\boldsymbol{\mu}=0$ has the form

[TABLE]

where the coefficients $C_{n,\boldsymbol{l}}(\boldsymbol{m})$ are the iterated integrals expressed in terms of the Goncharov’s polylogarithms Goncharov1998 ; Weinzierl:2002hv .

Let us note that the recurrence relations obtained in the previous section allow us to express via the same polylogarithmic functions also the expansions near any integer $D$ and even $\boldsymbol{\mu}$. In particular, we can do it for the expansion near any odd integer value of $D$ and integer exponents $\boldsymbol{\alpha}$.

In the points where at least one of $\mu_{l}$ is an odd number, the recurrence relations can not be used to express the corresponding expansion via polylogarithms. However, as we shall see soon, they can be used to obtain the quadratic relations for the expansion coefficients, see Fig. 2. These relations are closely related to the ones described in Ref. Lee2018 in the case when matrix $A$ is symmetric and proportional to $\epsilon$ or $\epsilon+1/2$.

In order to obtain these relations, we first note that the matrix $A(D,\boldsymbol{\mu})$ has the following symmetry

[TABLE]

Let $F_{\pm}$ be the two solutions of the equations

[TABLE]

Then, from Eq. (48) it follows that the bilinear form $F^{T}_{-}WF_{+}$ is independent of $\boldsymbol{m}$:

[TABLE]

Now we note that $A(-D,-\boldsymbol{\mu})=-A(D,\boldsymbol{\mu})$ and, therefore, $F_{-}(D,\boldsymbol{\mu})$ obeys the same equation as $F_{+}(-D,-\boldsymbol{\mu})$. So, given two solutions $T_{1,2}(D,\boldsymbol{\mu},\boldsymbol{m})$ of Eq. (19), we can construct the conserved bilinear form

[TABLE]

where $C_{12}(D,\boldsymbol{\mu})$ in the right-hand side denotes some quantity which, in general, depends on $D$ and $\boldsymbol{\mu}$, but not on $\boldsymbol{m}$.

Now let $T_{1,2}$ also satisfy the recurrence relations (24) and (27). Suppose that we are interested in the expansion near $D=D^{\star}$ and $\boldsymbol{\mu}=\boldsymbol{\mu}^{\star}$, where $2D^{\star}$ and $\mu_{l}^{\star}$ are integer. Then we can write

[TABLE]

where $C_{12}(D,\boldsymbol{\mu})$ is the same function as in Eq. (51) and

[TABLE]

Using the commutation relations derived in the Appendix A one can show that $B$ has the symmetry

[TABLE]

When expanded in $\epsilon$ and $\boldsymbol{\tau}$, the identity (52) provides a lot of quadratic constraints for the expansion coefficients of the solutions near $D=D^{\star}$ and $\boldsymbol{\mu}=\boldsymbol{\mu}^{\star}$. Equations (51) and (52) are important results of the present paper.

Let us now specialize Eq. (51) to the case $T_{1}=V^{(\varrho_{0},\boldsymbol{\rho})}$, $T_{2}=V^{(\tilde{\varrho}_{0},\boldsymbol{\tilde{\rho}})}$. In order to calculate the constant in the right-hand side, we consider the limit $m_{l}\to 0$ ($l=1,\ldots,L$) at fixed $m_{0}$. In this limit the main contribution comes from the region $m_{0}x\sim 1$ and from Eq. (42) we see that

[TABLE]

where $F$ is a function of $({m_{l}}/{m_{0}})^{2}$, analytic and nonzero at the origin. The value of $F$ at the origin can be obtained by replacing $F_{C}^{(L)}\to 1$ in Eq. (42). Since we know that the right-hand side of Eq. (51) does not depend on masses, it can be nonzero only if $\boldsymbol{\tilde{\rho}}=\boldsymbol{\rho}$. Indeed, the power $m_{l}^{\mu_{l}(\tilde{\rho}_{l}-{\rho}_{l})}$ can not be canceled if $\tilde{\rho}_{l}\neq\rho_{l}$. Besides, the exponent $1-\delta_{a_{l}\rho_{l}}$ in Eq. (55) shows that only the contribution of the components with $a_{l}=\rho_{l}$ ($l=1,2,\ldots,L$) may contribute to the limit $m_{l}\to 0$. Taking the integral over $x$, we have

[TABLE]

Using this asymptotics in the left-hand side of Eq. (51), we fix the right-hand side of this identity. We find

[TABLE]

Recall that $\bar{\varrho}_{0}=1-\varrho_{0}$. Let us note that in Ref. Goto2013 some bilinear relations between Lauricella functions $F_{C}$ have been derived using different methods. It would be interesting to compare them with Eq. (57).

Summing over $\varrho_{0}$ and $\tilde{\varrho}_{0}$ and using the identity

[TABLE]

we obtain

[TABLE]

Finally, summing over $\boldsymbol{\rho}$ and $\boldsymbol{\tilde{\rho}}$ and using Eq. (34) we have

[TABLE]

Note that we can obtain the bilinear relation separately for $\mathcal{P}$-even part $T_{+}=\frac{1+\mathcal{P}}{2}T$. In order to do this, in addition to Eq. (58), we also use the identity

[TABLE]

We obtain

[TABLE]

6 Special cases when the system acquires a block-triangular form

6.1 Removing analytic regularization

Let us first consider the special case $\mu_{l}=\mu=D-1$, which corresponds to the removal of the analytical regularization. The differential system (19) in this case should have a block-triangular form, corresponding to the decoupling of the trivial “clover-leaf” integrals arising from the contraction of one of $L+1$ lines. In order to see this structure, it is convenient to pass to the new functions

[TABLE]

The differential system for these functions has the form

[TABLE]

where $N=\sum_{l=0}^{L}n_{l}$. Note that Eq. (64) corresponds to $(\epsilon+1/2)$-form discussed in Ref. Lee2018 . In particular, this form allows one to write down the $\epsilon$-expansion of the solution near $D=1$ in terms of Goncharov polylogarithms (see Appendix B for details). Let us observe that the equations for each of the components

[TABLE]

decouple and acquire the form

[TABLE]

with the general solution being

[TABLE]

These are just the equations for the clover-leaf integrals. One can check explicitly that $Y_{10\ldots 0},\,Y_{01\ldots 0},\ldots,Y_{0\ldots 01}$ are indeed evaluated to the corresponding expressions for the clover-leaf integrals after substitution of Eq. (7) into Eq. (63). The simplest way to check this is to represent $Y(\mu)$ as

[TABLE]

and to note that the corresponding elements of $W$ are vanishing in this limit. Then, only the pole part of the corresponding components of $T(D-1,\mu)$ have to be evaluated (note that $W$ is a diagonal matrix).

Let us consider now the homogeneous differential system obtained by putting components (66) to zero. This is the system for the maximally cut tadpole diagram. It has the same form (64), where now the action of $H$ is restricted to the subspace with

[TABLE]

Note that in this subspace the operator $N-1$ is invertible. From now on to the end of this subsection, we assume that all operators are restricted to this subspace and that all vector functions belong to it. Then, using the property

[TABLE]

and following the same path as in the previous section, we obtain the bilinear constraint

[TABLE]

Let us now determine the basis of solutions in the subspace constrained by Eq. (70). We define

[TABLE]

One can check explicitly that $2^{L+1}-L-1$ such functions with $\varrho_{0}\neq|\boldsymbol{\rho}|$ satisfy the constraints (70) and, therefore, form the basis we are looking for444In particular, when $|\boldsymbol{\rho}|\geqslant 2$, the function Y^{(\boldsymbol{\rho})}=M\,V^{(\boldsymbol{\rho})}\big{|}_{\begin{subarray}{c}D=\mu+1\\ \mu_{l}=\mu\end{subarray}}=Y^{(0,\boldsymbol{\rho})}+Y^{(1,\boldsymbol{\rho})} is the homogeneous solution..

For two basis functions, $Y^{(\tilde{\varrho}_{0},\boldsymbol{\tilde{\rho}})\intercal}$ and $Y^{(\varrho_{0},\boldsymbol{\rho})}$, using Eqs. (57) and (69), we have

[TABLE]

Note that the right-hand side of this identity only makes sense when $\varrho_{0}\neq|\boldsymbol{\rho}|$, otherwise the argument of $\sin$ function in the denominator becomes zero.

In order to obtain the quadratic relations for the coefficients of $\epsilon$-expansion near $\mu=\mu^{\star}$ for some integer $\mu^{\star}$, we have to substitute $\mu=\mu^{\star}-2\epsilon$ in Eq. (74) and use the relation

[TABLE]

where \mathfrak{R}(\mu)=\mathfrak{R}(D,\boldsymbol{\mu})\big{|}_{\begin{subarray}{c}D=\mu\\ \mu_{l}=\mu\end{subarray}} is the operator of Eq. (29) restricted to the subspace (70).

Note that the quadratic relations obtained by substituting Eq. (75) into Eq. (74) are valid only for the solutions of homogeneous equations. A natural question arises: whether one can obtain similar relations also for generic solutions of inhomogeneous equations? Within our approach, the answer is negative. The bilinear relations (51) are valid for any two solutions $T_{1}$ and $T_{2}$ even if we put $D-1=\mu_{l}=\mu$. However, in order to obtain the quadratic relations near any integer point $\mu=\mu^{\star}\geqslant 0$, we have to shift the first argument in the left factor $T_{1}^{\intercal}$ at least twice by the operator $R^{-1}$. Then the second operator $R^{-1}(-\mu,-\mu)$ does not make sense (since $R(-\mu,-\mu)=-\mu M^{-1}(N-1)$ is not invertible).

6.2 Groups of identical lines

Let us now discuss the case when some lines share the same $\mu$ and $m$. Then we can write

[TABLE]

where $r_{k}$ is the number of lines with mass $m_{k}$ and parameter $\mu_{k}$, so that $\sum_{k}r_{k}=L+1$. There is a symmetry group $\cal G$ generated by the permutations of lines sharing the same mass and parameter. This symmetry, in particular, holds for the matrix $A$ in the right-hand side of the differential system, which now can be written as

[TABLE]

where

[TABLE]

Here $u_{k}=\sum_{s=0}^{k-1}r_{s}$, and the sequence $u_{k},u_{k}+1,\ldots,u_{k}+r_{k}-1$ enumerates the lines in $k$-th group. It is remarkable that $A$ is expressed solely via $S_{xk}$ and $S_{zk}$ which are among the standard generators of the Lie algebra ${\mathfrak{sl}}(2,\mathbb{C})$. We can express via the same generators other operators which we have used previously, in particular $\mathcal{P}=i^{L+1}\prod_{k=0}^{K}e^{-i\pi S_{zk}}$.

Note that there is a subtle point here that we want to stress. When all $r_{k}$ are even, the matrix $M$ obviously has zero eigenvalue and, therefore, is not invertible. Thus the matrix $A$ is ill-defined. So, from now on we assume that at least one $r_{k}$ is odd, i.e., at least one group contains odd number of lines.

The symmetry of the matrix $A$ leads to the splitting of the system into separate subsystems, each corresponding to a specific irreducible representation of $\cal G$. Since the tadpole integral (76) is invariant under the action of this group, we are mostly interested in the subspace spanned by the tensors $T_{\boldsymbol{a}}$ symmetric with respect to all permutations of indices from this symmetry group.

Alternatively, we can reduce the number of auxiliary functions from $2^{L+1}=2^{\sum_{k}r_{k}}$ to $\prod_{k}(r_{k}+1)$ from the very beginning. Namely, we can introduce the functions

[TABLE]

where $\boldsymbol{a}=a_{0}\ldots a_{K}$ and $a_{k}$ runs from [math] to $r_{k}$, and ${n\choose k}=\frac{n!}{k!(n-k)!}$ is the binomial coefficient.

As before, it is convenient to treat the quantities $T_{\boldsymbol{a}}$ as components of the vector $T$ in $\mathbb{C}^{r_{0}+1}\otimes\ldots\otimes\mathbb{C}^{r_{K}+1}$. Then the operators $S_{xk}$ and $S_{zk}$ act on $k$-th factor as usual spin operators, i.e., the matrices with nonzero elements being

[TABLE]

Again, thanks to the parity operator the differential system consisting of $\prod_{k}(r_{k}+1)$ equations splits into two subsystems corresponding to $\mathcal{P}=+1$ and $\mathcal{P}=-1$, the first involving components with $|\boldsymbol{a}|=\sum_{k}a_{k}$ even, the second with $|\boldsymbol{a}|$ odd. Note that since we assume that at least one $r_{k}$ is odd, the numbers of even and odd components coincide and are equal to $\frac{1}{2}\prod_{k}(r_{k}+1)$.

Let us now discuss the parameter and dimension shifting relations. Here we have to take into account that shifting separately each parameter is not compatible with the permutation symmetry. Fortunately, for the dimension shift at fixed $\alpha_{k}$ it is sufficient to have possibility to shift all parameters simultaneously. Therefore, we need to define the action of

[TABLE]

on the symmetric subspace. After some transformations we obtain

[TABLE]

where $S_{-k}=S_{xk}-iS_{yk}=S_{xk}+[S_{xk},S_{zk}]$ which is the matrix with nonzero elements being

[TABLE]

Let us now construct the basis of solutions. We remind that we restrict ourselves to the case when at least one of $r_{k}$ is odd. We assume that $r_{0}$ is odd. Similarly to Section 4, we want to replace in Eq. (76) some $P_{0}$ with $Q_{0}^{(0)}$ or $Q_{0}^{(1)}$ defined in Eq. (32). Assuming that $m_{0}>\sum_{k=1}^{K}r_{k}m_{k}$, for convergence of the integral over $x$ we keep $[P_{0}(\mu_{0},m_{0},x)]^{\lceil r_{0}/2\rceil}$ and replace all other $P_{0}$ with $Q_{0}^{(0)}$ or $Q_{0}^{(1)}$. Thus we have

[TABLE]

where $\rho_{0}=0,1,\ldots,\lfloor r_{0}/2\rfloor$, $\rho_{k}=0,1,\ldots,r_{k}$ ($k>0$), and $\boldsymbol{\rho}=\rho_{1},\ldots\rho_{K}$. Note that since $r_{0}$ is odd, we have $\lceil r_{0}/2\rceil=(r_{0}+1)/2$, $\lfloor r_{0}/2\rfloor=(r_{0}-1)/2$. Then there are $(\lfloor r_{0}/2\rfloor+1)\prod_{k=1}^{K}(r_{k}+1)=\tfrac{1}{2}\prod_{k=0}^{K}(r_{k}+1)$ functions $V_{\boldsymbol{0}}^{(\boldsymbol{\rho})}$. This is exactly the correct dimension of $\mathcal{P}$-even part of solution space. If we modify similarly other components of $T_{\boldsymbol{a}}$ from Eq. (80), we immediately see that we have exactly $\tfrac{1}{2}\prod_{k=0}^{K}(r_{k}+1)$ basis vectors also in $\mathcal{P}$-odd part of solution space, which also coincides with the correct dimension of this subspace. Note that the components other than $T_{\boldsymbol{0}}$ involve proper symmetrization within each group of equivalent lines. For the sake of completeness we present here the result:

[TABLE]

Again, we have the quadratic relations of the form

[TABLE]

If $r_{0}=1$, the function $C(\rho_{0},\tilde{\rho}_{0},\boldsymbol{\rho}|D,\boldsymbol{\mu})$ in the right-hand side can be easily deduced from Eq. (59):

[TABLE]

However, it is not quite clear how to fix this function for $r_{0}>1$, and we leave this question for further investigations. The reason why the case $r_{0}=1$ is simple for us is because, while constructing the basis of solutions, we always assumed that one mass ($m_{0}$) is larger than the sum of all others. Therefore we did not have obligations to consider the questions of analytical continuation of the obtained solutions across (or around) the singular hypersurfaces defined by Eq. (17). However, if $r_{0}>1$, we would have to do it.

Finally, we note that if we remove the analytical regularization, the $\mathcal{P}$-even subsystem contains $K+1$ decoupled equations and the number of its homogeneous solutions becomes $\frac{1}{2}\prod_{k}(r_{k}+1)-K-1$.

6.3 Removing dimensional regularization at $D=2$

Let us now discuss some important features of the obtained results at $D=\mu+1=2$. Below we assume that $L>1$. The matrix $A$ in this case has the form

[TABLE]

The second term is diagonal and it is clear that the equations for the components with $|\boldsymbol{a}|\neq 2$ form a closed subsystem. On the other hand, if we pass to $Y$, Eq. (63), we obtain the system

[TABLE]

with $H$ defined in Eq. (65). In this form, the equations for $Y_{\boldsymbol{a}}$ components with $|\boldsymbol{a}|=1$ decouple. Let us consider the linear transformation of column of functions $T_{\boldsymbol{a}}$ which replaces $L+1$ entries $T_{110\ldots 0},\,T_{1010\ldots 0},\ldots,\,T_{10\ldots 10},\,T_{10\ldots 01}$, and $T_{0110\ldots 0}$ with $Y_{100\ldots 0},\,Y_{010\ldots 0},\ldots,Y_{0\ldots 01}$. It is easy to check that this transformation is non-degenerate. With the new column of functions $\tilde{T}$ the differential system acquires the block-triangular form depicted in Fig. 3.

Let us observe that the component $T_{\boldsymbol{0}}$ enters the closed system of $2^{L}-(L+1)L/2$ equations determined by the block $B_{2}$. Therefore, in a certain sense, we find that the number of master integrals at $D=2$ is equal to $2^{L}-(L+1)L/2$. This is in agreement with known cases $L\leqslant 4$, Adams2013 ; TancrediPC2019 .

Let us comment about the number of master integrals at $D=2$ when there are groups of identical masses. In the notations of previous subsection, we obtain that the number of master integrals is equal to $\frac{1}{2}\prod_{k}(r_{k}+1)-(K+1)(K+2)/2+\sum_{k}\delta_{r_{k},1}$. In particular, for $L$-loop sunrise graph with equal masses we have $L+2-1-1=L$ master integrals.

Basis of solutions.

Constructing the basis of solution space at $D=2$ appears to be extremely tricky and deserves a dedicated consideration elsewhere. Here we will only explain the complications that arise on the way.

First, one might be tempted to pass to the basis of functions $V^{[\boldsymbol{\upsilon}]}$ defined as

[TABLE]

where $\boldsymbol{\upsilon}=\upsilon_{1}\upsilon_{2}\ldots\upsilon_{L}$ is a binary number. The functions of this basis have the following limit of $\boldsymbol{a}=\boldsymbol{0}$ component:

[TABLE]

However, other components of $V^{[\boldsymbol{\upsilon}]}$ diverge in the limit $\epsilon\to 0$ unless $\boldsymbol{\upsilon}=1\ldots 1$ because fo the singular asymptotics of $K_{1}$ function at small argument.

In order to get rid of the divergent overall factor in Eq. (31) at $\mu\to 1$, let us define the functions

[TABLE]

Then if we naively take limit $\epsilon\to 0$ under the integral sign, we can come to a wrong conclusion that all $U^{(\boldsymbol{\rho})}$ tend to one and the same expression independent of $\boldsymbol{\rho}$:

[TABLE]

However, Eq. (95) is correct only for $U^{(\boldsymbol{0})}$. The reason why it breaks down for $\boldsymbol{\rho}\neq\boldsymbol{0}$ is that the expansion of $I_{-1+\epsilon}(mx)$ in $x$ starts from the singular term $\Gamma[\epsilon]^{-1}(mx/2)^{\epsilon-1}$. Although this term is suppressed in $\epsilon$, it may survive after the integration over $x$ due to its singular nature. More precisely, consider the identities

[TABLE]

They show that when $\alpha$ is close to $-n$ ($n=0,1,\ldots$), the integral over $x$ gives poles in $\epsilon$. A thorough investigation allows us to establish that nonzero corrections to the naive limit (95) appear only in the components with $|\boldsymbol{a}|=2$. Moreover, we can explicitly calculate these corrections using Eqs. (6.3):

[TABLE]

Here we imply that $\boldsymbol{\rho}\neq\boldsymbol{0}$ and that the first term in square brackets should be omitted when $|\boldsymbol{\rho}|=1$. We can construct the solutions which have nonzero $\epsilon\to 0$ limit exactly in one component with $|\boldsymbol{a}|=2$. Let us adopt the notation $\boldsymbol{1}_{l_{1},l_{2},\ldots}$ for the vector $\boldsymbol{\rho}$ with $\rho_{l_{1}}=\rho_{l_{2}}=\ldots=1$ and other components equal to zero. Then we define

[TABLE]

with the property $[U(b,c)]_{\boldsymbol{a}}=\delta_{\boldsymbol{a},\boldsymbol{1}_{b,c}}\frac{\prod_{l=0}^{L}(2m_{l})}{m_{b}m_{c}}$ (here $0\leqslant b<c\leqslant L$). These ${L+1\choose 2}=L(L+1)/2$ solutions are obviously the solutions of the homogeneous equations for the components with $N=2$. Note that any difference $U^{\left(\boldsymbol{\rho}\right)}-U^{\left(\boldsymbol{0}\right)}$ tends to a linear combination of $U(b,c)$. Therefore, we construct the combinations that vanish in the limit $\epsilon\to 0$,

[TABLE]

for any $\boldsymbol{\rho}$ with $|\boldsymbol{\rho}|\geqslant 3$. The idea is that we can now define the solutions

[TABLE]

In particular, if $L=3$, we obtain that

[TABLE]

has a finite limit $\epsilon\to 0$, while each term in the square brackets contains divergences in components with $|\boldsymbol{a}|=2$. Unfortunately, at $L>3$ the $2^{L}-L(L+1)/2-1$ solutions $J^{\left(\boldsymbol{\rho}\right)}$ appear to be not linearly independent and, in order to proceed to the full basis of solutions, we have to determine on the next step the linear combinations of $J^{\left(\boldsymbol{\rho}\right)}/\epsilon$ which vanish in the limit $\epsilon\to 0$. The combinatorics appears to be quite involved and we reserve the full consideration for our future work.

6.4 Equal mass sunrise integral in $D=2$

The case of equal mass $L$-loop sunrise integral in $D=2$ has probably received the most attention in the literature, see, in particular, Refs. Muller-StachWeinzierlZayadeh2012 ; Vanhove:2014wqa ; BlochKerrVanhove2014 ; Primo2017 ; Adams2017b ; Broedel2018b .

So, we consider the two groups of lines, one containing a single ($r_{0}=1$) line with $m_{0}=\sqrt{q^{2}}$, the other consisting of $r_{1}=L+1$ lines with mass $m_{1}=m$. It is convenient to put $q^{2}=1$.

Similar to Eq. (93), we define the solutions $V^{[\upsilon]}$ ($\upsilon=0,1,\ldots,L+1$) as follows

[TABLE]

where $V^{(\boldsymbol{\rho})}=V^{(\rho_{0}\rho_{1})}$ are defined in Eq. (87). Note that when $\upsilon\geqslant 2$ the function $V^{[\upsilon]}$ is the solution of the homogeneous system for the ‘cut’ sunrise integral. This is because it is a linear combination of $V^{(\boldsymbol{\rho})}$ with $|\boldsymbol{\rho}|>1$ which are all the homogeneous solutions, see footnote 4 on page 4. Again, $00$ components of $V^{[\upsilon]}$ all have a finite limit at $D=2$:

[TABLE]

As we explained in the previous section, the limit of components with $\boldsymbol{a}\neq 00$ of each $V^{[\upsilon]}$ can contain divergences. However, we have some freedom in definition which might help us to get rid of these divergences. Namely, we can construct

[TABLE]

where $C_{\upsilon\omega}$ are some coefficients with the asymptotics $C_{\upsilon\omega}(\epsilon)=O(\epsilon^{\omega})$. These new functions $\widetilde{V}^{[\upsilon]}$ have the same properties as ${V}^{[\upsilon]}$. Namely, they have the same limit of $00$ component as ${V}^{[\upsilon]}$ and are the solutions of the homogeneous system for the cut sunrise graph. Now the idea is that we can fix the coefficients $C_{\upsilon\omega}(\epsilon)$ in such a way as to eliminate all divergences in the limit $\epsilon\to 0$. We were able to check this for a few low-loop cases555Note that this is not the case for the sunrise integral with different masses..

Instead of running into complications connected with the basis $V^{(0\rho_{1})}$, we can pass to the basis of momenta spanned by the functions $\mathrm{IKM}[\{0,1\}_{1},\{\upsilon,L+1-\upsilon\}_{m},s]$ defined as

[TABLE]

Similar functions appear in Refs. Broadhurst:2016myo ; Zhou:2017jnm ; Broadhurst:2018tey . For any specified $L$ it is easy to construct the matrix $\mathcal{M}$ for transition to moments basis using the matrix $R$ which shifts the power of $x$ in the integrand. Namely, the $n$-th row of the matrix $\mathcal{M}$ has the form

[TABLE]

where $[\ldots]_{0}$ denotes the zeroth row of the matrix in the brackets. In Appendix C we present a few examples of the quadratic relations obtained in this way. These quadratic relations resemble those discovered by Broadhurst and Roberts in Ref. Broadhurst:2018tey . The difference is that our relations explicitly depend on one parameter $q^{2}/m^{2}$ while the relations discovered in Ref. Broadhurst:2018tey are obtained for the symmetric point $q^{2}=m^{2}$. However, it is not possible to obtain the latter from the former by simply putting $q^{2}=m^{2}$. One obstacle is that in the symmetric point (when $q^{2}=m^{2}$) the number of independent moments drops by about a half. Another related problem is that our identities tend to develop singularities at $q^{2}=m^{2}$.

Alternatively, we can start from Eq. (88) for $K=0$, use the shift operators $R$, Eq. (27), and $\mathcal{R}_{k}$, Eq. (83), and the matrix $\mathcal{M}$ above. In this way we’ve been able to reproduce the matrix $D_{N}$ in the left-hand side of Eq. (5.1) of Ref. Broadhurst:2018tey up to at least $10$ loops. However, in order to obtain the right-hand side of the same equation, we need the generalization of (89) to the case $r_{0}>1$.

7 Conclusion

In the present paper we considered the $L$-loop watermelon and sunrise integrals regularized both dimensionally and analytically. We have constructed a special set of “master integrals” for which we have managed to derive a Pfaffian differential system.

The basis of solutions of this system is expressed in terms of the Lauricella function $F_{C}$ with $L$ independent arguments and generic parameters (cf. Berends:1993ee ; Kalmykov:2016lxx ). Thus, we obtain the Pfaffian differential system equivalent to the Lauricella $F_{C}^{(L)}$ second-order differential system. To the best of our knowledge, the Pfaffian differential system for $F_{C}^{(L)}$ was previously known only for $L\leqslant 2$. It worth to note that the obtained system does not have apparent singularities.

Using the symmetry properties of the matrix in the right-hand side of the differential system and following the path similar to that in Ref. Lee2018 , we obtain the bilinear relations between the solutions of the system for opposite signs of $D$ and $\boldsymbol{\mu}$.

We derive the operators which shift dimension and/or exponents of the propagators by $\pm 1$. These operators, in particular, allow one to obtain the canonical form of the differential system near any odd integer value of $D$ and integer exponents $\boldsymbol{\alpha}$. Therefore, the expansion of the solutions in the dimensional and analytic regularization parameters in this case can be written in terms of polylogarithmic functions. Thanks to the shift operators, the obtained bilinear relations lead to the quadratic relations between the coefficients of expansion of the solutions near any integer point in $D,\boldsymbol{\mu}$ space.

We have considered the integrals with integer exponents of the denominators. We observe that the quadratic relations for the solutions of the whole differential system do not exist in this case. Instead, we have derived the quadratic relations for the solutions of the homogeneous differential system.

A Pfaffian differential system can be easily restricted to any subvariety in the space of variables. In particular, we have considered the case when the masses and exponents of some lines coincide. In this case the matrix in the right-hand side of the differential system can be naturally rewritten in terms of the generators of $\mathfrak{sl}(2,\mathbb{C})$ acting in the $(r_{k}+1)$-dimensional irreducible representation ($r_{k}$ is the number of identical lines in the $k$-th group).

Finally, we have considered the differential system for the tadpole integral in two dimensions ($D=2,\,\alpha_{l}=1$). The system acquires a block-triangular form, with the nontrivial block of size $2^{L}-(L+1)L/2$. Using the quadratic relations for the maximally cut sunrise integral with equal masses in $D=2$ we obtain for any given $L$ the quadratic relations for the moments of the product of the Bessel and Macdonald functions. However, we were not able to derive the closed formula for these relations for general $L$. The found relations, in a sense, generalize the quadratic relations obtained by Broadhurst and Roberts Broadhurst:2018tey to the case of arbitrary incoming momentum. However, it is not easy to obtain the latter from the former.

The consideration of the present paper, while satisfactorily describing the generic case, can not be considered as exhaustive when it concerns special cases. In particular, we don’t fully understand how to construct the basis of solution space and the quadratic relations for the case $D=2$. We also did not consider the special configurations of masses when $M$ is not invertible. In particular, it happens in physically relevant case when each group contains even number of identical lines. We underline that once the number of loops $L$ is fixed, one can rely on a brute-force approach to tackle the above issues. The principal problem is to solve them for general $L$.

Acknowledgements.

We thank Mikhail Kalmykov for many useful discussions and David Broadhurst for valuable comments on the manuscript. R. Lee would like to express special thanks to the Mainz Institute for Theoretical Physics (MITP) of the Cluster of Excellence PRISMA+ (Project ID 39083149) for its hospitality and support. R. Lee is supported by the grant of the “Basis” foundation for theoretical physics.

Appendix A Compatibility conditions

Let us examine the compatibility conditions of Eqs. (19), (24), and (27). First, it is easy to see that $dA=0$, thanks to the identity $[\sigma_{l},\sigma_{k}]=0$. Then we have to check that $[A_{l},A_{k}]=0$, or equivalently, $A\wedge A=0$. We have

[TABLE]

Other conditions have the form

[TABLE]

These identities can be checked along the same lines as Eq. (108).

Let us also write down two useful identities which allow one to check the symmetry property (54) of the matrix $B$ in the right-hand side of the quadratic relation (52). They read

[TABLE]

Appendix B Expansion near $D=1$

We will consider $\epsilon$-expansion of $Y$ for $D=1-2\epsilon$. Then the $\epsilon$-expansion of $T$ can be obtained with the help of Eq. (63). Let us rescale $m_{k}\to m_{k}x$ ($k=1,\ldots,L$) and consider the differential equation with respect to $x$:

[TABLE]

where the set (or alphabet) $\Lambda$ contains [math] and all letters of the form

[TABLE]

where $L$-tuple $\{\eta_{1},\ldots,\eta_{L}\}$ has elements $\pm 1$, and

[TABLE]

Then we choose the fundamental matrix of solutions of Eq. (116) as

[TABLE]

This definition corresponds to the small-$x$ asymptotics

[TABLE]

The $\epsilon$-expansion of the fundamental matrix has the form

[TABLE]

where the sum $\sum_{a_{1}\ldots a_{n}}$ runs over all words of length $n$ with letters from the alphabet $\Lambda$, and $G(a_{1},\ldots,a_{n}|x)$ denotes the Goncharov’s polylogarithms Goncharov1998 ; Weinzierl:2002hv defined recursively by

[TABLE]

Let us represent the specific solution $Y_{a_{0}\boldsymbol{\rho}}$ as

[TABLE]

where $C_{a_{0}\boldsymbol{\rho}}^{(\varrho_{0},\boldsymbol{\rho})}$ are the integration constants. In order to find $C_{a_{0}\boldsymbol{\rho}}^{(\varrho_{0},\boldsymbol{\rho})}$, we consider the asymptotics $x\to 0$. Using Eqs. (56) and (69) we have

[TABLE]

Appendix C Examples of quadratic relations for $\mathrm{IKM}$ functions.

Let us present a few examples of the quadratic relations for $\mathrm{IKM}$ functions related to the cut sunrise integral. For $L$-loop sunrise integral the overall number of Bessel and Macdonald functions in the integrand is equal to $L+2$. We use the $\mathrm{IKM}$ functions as defined in Eq. (106).

For $L=2$ we obtain the following nontrivial quadratic relation:

[TABLE]

For $L\geqslant 3$ we obtain similar relations which are too lengthy to be presented here. These relations, however, simplify on the pseudo-thresholds.

For $L=3$ we have

[TABLE]

[TABLE]

[TABLE]

For $L=4$ we have

[TABLE]

with

[TABLE]

At $m=1/3$ we have one relation

[TABLE]

We note that we could have easily proceeded to more loops.