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

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.


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 -expansion can not be expressed in terms of the standard functional basis, viz., multiple polylogarithms [1,2]. Already starting from two loops one observes complicated iterated integrals involving modular forms and/or elliptic curves [3,4]. For L loops, the expansion of massive sunrise integrals is likely to involve periods of certain (L − 1)-dimensional Calabi-Yau varieties (cf. Refs. [5,6]).
For the general case of different masses the closed expressions for these integrals are known since Ref. [7] (see also Ref. [8]) in terms of Lauricella functions F (L) C defined via multiple hypergeometric series. Although these functions have a long history of investigation, starting from [9], 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 (L) C 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. [7,8], we find the fundamental system of solutions in terms of the Lauricella functions F (L) C with generic parameters and arguments. In this way, our results give a Pfaffian system for Lauricella F (L) C for all L, to the best of our knowledge, for the first time 1 .
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 (L) C 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. [11]). 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. Here we imply the Lorentzian signature, so that p 2 = p 2 0 − p 2 . Note that we have chosen a somewhat unconventional overall normalization factor 2 D−1 Γ D 2 l m l Γ (α l ) in order to simplify some expressions below.

Differential equations in Pfaffian form
Before moving any further let us discuss the applicability of our present consideration to the sunrise integral S(m 0 , . . . , m L |q 2 ) = 2 (2.2) 1 For the case L = 2 a Pfaffian system was found in Ref. [10].  The obvious relation between S and T is S| q 2 =0 = T . However, we can also represent the (L − 1)-loop sunrise integral as an L-loop watermelon integral with one cut line: where m 0 = 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 D−1 0 S(m 1 , . . . , m L |m 2 0 ) satisfies the same differential equations as T (m 0 , . . . , m L ) with α 0 = 1. Therefore, in the rest of the paper we will concentrate on the watermelon integral T only.
After the Wick rotation we can rewrite Eq. (2.1) as where now p 2 = p 2 0 + p 2 . Representing the function δ(p) as As we shall see immediately, the integral d D p π D/2 mΓ(α)e ipx (p 2 +m 2 ) α depends on D and α only via the combination µ = D − 2α + 1 and we have, therefore, introduced Here K ν (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) a = a 0 a 1 . . . a L , where

8)
P 0 is defined in eq. (2.6), and The original integral (2.4) is recovered as T = T 0 (D, µ, m). 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 a . Note that all components of T a can be expressed via T with shifted dimension and indices: where |a| = l a l . It is convenient to treat the 2 L+1 quantities V a as components of the vector V in C 2 ⊗ . . . ⊗ C 2

L+1 factors
, where V a can be either T a or t a . In what follows, we assume that the factors in C 2 ⊗. . .⊗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 a we differentiate P 0,1 (µ, m, x) with respect to m and use the differential equation for K ν (y). We obtain In matrix notations we may rewrite this as m∂ m P = (mxσ + µn)P , (2.12) where P = P 0 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 ∞ 0 dx x D t. Note the additional factor x in the integrand as compared to the definition (2.7). In order to get rid of this factor, we consider the identity dx∂ x (x D t) = 0. Explicitly differentiating and using the relation x∂ x P = (m∂ m − µ)P = (mxσ − µn)P , (2.13) which, again, follows from the differential equation for K ν , we obtain (2.14) Here 15) and the operators σ l and n l act in l-th factor of C 2 ⊗ . . . ⊗ C 2 as σ and n, respectively. Therefore, for some choice of signs η l . Thus, using Eqs. (2.12) and (2.16) we obtain the following expression for the derivative of T with respect to m k Equivalently, the above equations can be written as Remarkably, the differential form A is closed (dA = 0) which can be easily proved if one takes into account the pairwise commutativity of σ l . Then the integrability condition dA = A ∧ A requires that A ∧ A = 0, see Appendix A for the derivation. The system (2.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. (2.17) (cf. Ref. [11]). The latter hyperplanes correspond to vanishing of the eigenvalues of the matrix M .
Let us note that the differential system (2.19) in fact splits into two separate systems, each for 2 L functions T a with |a| being either even or odd. This due to the commutativity of A with the parity operator (2.21) Of course, we are mostly interested in the subsystem for even components which involves our original tadpole integral T .

Shifting exponents and dimension
Let us now obtain the recurrence relations in α k and D. For the former we use the identities Then, using Eq. (2.18), we have In order to obtain the operator which shifts D at fixed α, we should first use the operators R l to shift µ → µ − 2 and then shift the dimension by −2 at fixed µ using the operator where is the operator which shifts D by −2 at fixed α. Here and below e = L j=0 e j = (1, 1, . . . , 1) and the product of matrices corresponds to the multiplication from the left to the right, so that L l=0 R l = R 0 . . . R L .

Basis of solution space
While the function T as defined in Eq. (2.7) is a specific solution of the differential system (2.19) and recurrence relations (3.3) and (3.6), 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 ν are replaced with I ±ν sin πν . More specifically, let us make the following replacement in the definition of t a , Eq. (2.8): where ρ = ρ 1 ρ 2 . . . ρ L is a binary number. We assume that m l > 0 and m 0 > L l=1 m l , so that the integrand decays exponentially at large x. Using the identity The set of functions (4.2) contains 2 L functions numbered by the vector ρ. Meanwhile, the solution space of the differential system (2.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 P, Eq. (2.21). Therefore, we can define functions which are the solutions of the differential system. Note that we have introduced the additional factor (−1) |ρ| for further convenience. In terms of the components we have so that the solutions with 0 = |ρ| (mod 2) have nonzero components with even |a|.
Together with the different choices of ρ, this gives us 2 L+1 solutions of the differential system (2.19). It is also easy to see that PR l = R l P and, therefore, V ( 0 ,ρ) is also a solution of the µ-shifting recurrences (3.3). As to the the operator R, Eq. (3.6), we have PR = −RP and therefore where we used the notationx = 1−x. Thus, the functions (4.6) satisfy the relations shifting D by 2, which follow from Eq. (3.6) 2 . Note that our basis functions V ( 0 ,ρ) can be expressed via Lauricella function F (L) C (cf. [7,8]). This function is defined via the hypergeometric series where (a) k = Γ(a + k)/Γ(a) is the Pochhammer symbol. In order to do this we use the expansion of the modified Bessel function I ν : .
Substituting this expansion in Eq. (4.2) and taking the integral over x with the help of the identity where (4.14) The functions V ( 0 ,ρ) can be obtained from Eq. (4.13) with the help of Eq. (4.7). Note that where x l = (m l /m 0 ) 2 and the parameters are generic for generic µ and D. Therefore, V with generic parameters. Recall that the system (2.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 (L) C .

Bilinear relations and expansion near integer D and µ.
The fundamental matrix of the differential system (2.19) can be written in the form of path-ordered exponent where C is some path in L + 1-dimensional space of mass parameters starting from some fixed point m 0 and ending in m 3 . Note that A is linear homogeneous in D and µ which, for the expansion near D = µ = 0, corresponds to a certain generalization of the canonical form of Ref. [12] to the case of analytical regularization. Thanks to this property, the expansion around the point D = µ = 0 has the form where the coefficients C n,l (m) are the iterated integrals expressed in terms of the Goncharov's polylogarithms [1,2]. 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 µ. In particular, we can do it for the expansion near any odd integer value of D and integer exponents α.
In the points where at least one of µ 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. [13] in the case when matrix A is symmetric and proportional to or + 1/2. In order to obtain these relations, we first note that the matrix A(D, µ) has the following symmetry Let F ± be the two solutions of the equations Then, from Eq. (5.3) it follows that the bilinear form F T − W F + is independent of m: Now we note that A(−D, −µ) = −A(D, µ) and, therefore, F − (D, µ) obeys the same equation as F + (−D, −µ). So, given two solutions T 1,2 (D, µ, m) of Eq. (2.19), we can construct the conserved bilinear form where C 12 (D, µ) in the right-hand side denotes some quantity which, in general, depends on D and µ, but not on m. Now let T 1,2 also satisfy the recurrence relations (3.3) and (3.6). Suppose that we are interested in the expansion near D = D and µ = µ , where 2D and µ l are integer. Then we can write where C 12 (D, µ) is the same function as in Eq. (5.6) and Using the commutation relations derived in the Appendix A one can show that B has the symmetry When expanded in and τ , the identity (5.7) provides a lot of quadratic constraints for the expansion coefficients of the solutions near D = D and µ = µ . Equations (5.6) and (5.7) are important results of the present paper.
Let us now specialize Eq. (5.6) to the case . In order to calculate the constant in the right-hand side, we consider the limit m l → 0 (l = 1, . . . , L) at fixed m 0 . In this limit the main contribution comes from the region m 0 x ∼ 1 and from Eq. (4.13) we see that 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 . Since we know that the right-hand side of Eq. (5.6) does not depend on masses, it can be nonzero only ifρ = ρ. Indeed, the power m µ l (ρ l −ρ l ) l can not be canceled ifρ l = ρ l . Besides, the exponent 1−δ a l ρ l in Eq. (5.10) shows that only the contribution of the components with a l = ρ l (l = 1, 2, . . . , L) may contribute to the limit m l → 0. Taking the integral over x, we have Using this asymptotics in the left-hand side of Eq. (5.6), we fix the right-hand side of this identity. We find Recall that¯ 0 = 1 − 0 . Let us note that in Ref. [11] some bilinear relations between Lauricella functions F C have been derived using different methods. It would be interesting to compare them with Eq. (5.12).
Summing over 0 and˜ 0 and using the identity we obtain (5.14) Finally, summing over ρ andρ and using Eq. (4.5) we have (5.15) Note that we can obtain the bilinear relation separately for P-even part T + = 1+P 2 T . In order to do this, in addition to Eq. (5.13), we also use the identity (5.16) We obtain 6 Special cases when the system acquires a block-triangular form

Removing analytic regularization
Let us first consider the special case µ l = µ = D − 1, which corresponds to the removal of the analytical regularization. The differential system (2.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 D=µ+1 µ l =µ (6.1) The differential system for these functions has the form dY = µHY , (6.2) where N = L l=0 n l . Note that Eq. (6.2) corresponds to ( + 1/2)-form discussed in Ref. [13]. In particular, this form allows one to write down the -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 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, µ) have to be evaluated (note that W is a diagonal matrix). Let us consider now the homogeneous differential system obtained by putting components (6.4) to zero. This is the system for the maximally cut tadpole diagram. It has the same form (6.2), where now the action of H is restricted to the subspace with 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 and following the same path as in the previous section, we obtain the bilinear constraint Let us now determine the basis of solutions in the subspace constrained by Eq. (6.8). We define One can check explicitly that 2 L+1 −L−1 such functions with 0 = |ρ| satisfy the constraints (6.8) and, therefore, form the basis we are looking for 4 . For two basis functions, Y (˜ 0 ,ρ) and Y ( 0 ,ρ) , using Eqs. (5.12) and (6.7), we have . (6.12) Note that the right-hand side of this identity only makes sense when 0 = |ρ|, otherwise the argument of sin function in the denominator becomes zero. In order to obtain the quadratic relations for the coefficients of -expansion near µ = µ for some integer µ , we have to substitute µ = µ − 2 in Eq. (6.12) and use the relation is the operator of Eq. (3.8) restricted to the subspace (6.8).
Note that the quadratic relations obtained by substituting Eq. (6.13) into Eq. (6.12) 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 (5.6) are valid for any two solutions T 1 and T 2 even if we put D − 1 = µ l = µ. However, in order to obtain the quadratic relations near any integer point µ = µ 0, we have to shift the first argument in the left factor T 1 at least twice by the operator R −1 . Then the second operator R −1 (−µ, −µ) does not make sense (since R(−µ, −µ) = −µM −1 (N − 1) is not invertible).

Groups of identical lines
Let us now discuss the case when some lines share the same µ and m. Then we can write where r k is the number of lines with mass m k and parameter µ k , so that k r k = L + 1.
There is a symmetry group 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 where S xk = 1 2 Here u k = k−1 s=0 r s , and the sequence u k , u k + 1, . . . , 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 sl(2, C). We can express via the same generators other operators which we have used previously, in particular P = i L+1 K k=0 e −iπ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 G. Since the tadpole integral (6.14) is invariant under the action of this group, we are mostly interested in the subspace spanned by the tensors T 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 k r k to k (r k + 1) from the very beginning. Namely, we can introduce the functions where a = a 0 . . . a K and a k runs from 0 to r k , and n k = n! k!(n−k)! is the binomial coefficient. As before, it is convenient to treat the quantities T a as components of the vector T in C r 0 +1 ⊗ . . . ⊗ 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 (S xk ) l,l−1 = (S xk ) l−1,l = 1 2 l(r k + 1 − l) (l = 1, . . . r k ) , (6.19) (S zk ) ll = r k /2 − l (l = 0, . . . r k ) . (6.20) Again, thanks to the parity operator the differential system consisting of k (r k + 1) equations splits into two subsystems corresponding to P = +1 and P = −1, the first involving components with |a| = k a k even, the second with |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 1 2 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 α k it is sufficient to have possibility to shift all parameters simultaneously. Therefore, we need to define the action of on the symmetric subspace. After some transformations we obtain which is the matrix with nonzero elements being 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. (6.14) some P 0 with Q where ρ 0 = 0, 1, . . . , r 0 /2 , ρ k = 0, 1, . . . , r k (k > 0), and ρ = ρ 1 , . . . ρ K . Note that since r 0 is odd, we have r 0 /2 = (r 0 + 1)/2, r 0 /2 = (r 0 − 1)/2. Then there are 0 . This is exactly the correct dimension of P-even part of solution space. If we modify similarly other components of T a from Eq. (6.18), we immediately see that we have exactly 1 2 K k=0 (r k + 1) basis vectors also in P-odd part of solution space, which also coincides with the correct dimension of this subspace. Note that the components other than T 0 involve proper symmetrization within each group of equivalent lines. For the sake of completeness we present here the result: (6.26) If r 0 = 1, the function C(ρ 0 ,ρ 0 , ρ|D, µ) in the right-hand side can be easily deduced from Eq. (5.14): 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. (2.17). However, if r 0 > 1, we would have to do it. Finally, we note that if we remove the analytical regularization, the P-even subsystem contains K + 1 decoupled equations and the number of its homogeneous solutions becomes 1 2 k (r k + 1) − K − 1.

Removing dimensional regularization at D = 2
Let us now discuss some important features of the obtained results at D = µ+1 = 2. Below we assume that L > 1. The matrix A in this case has the form The second term is diagonal and it is clear that the equations for the components with |a| = 2 form a closed subsystem. On the other hand, if we pass to Y , Eq. (6.1), we obtain the system dY = HY It is easy to check that this transformation is non-degenerate. With the new column of func-tionsT the differential system acquires the block-triangular form depicted in Fig. 3. Let us observe that the component T 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 4, [14,15]. 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 1 2 k (r k + 1) − (K + 1)(K + 2)/2 + k δ 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 [υ] defined as  where υ = υ 1 υ 2 . . . υ L is a binary number. The functions of this basis have the following limit of a = 0 component: However, other components of V [υ] diverge in the limit → 0 unless υ = 1 . . . 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. (4.2) at µ → 1, let us define the functions (6.32) Then if we naively take limit → 0 under the integral sign, we can come to a wrong conclusion that all U (ρ) tend to one and the same expression independent of ρ: However, Eq. (6.33) is correct only for U (0) . The reason why it breaks down for ρ = 0 is that the expansion of I −1+ (mx) in x starts from the singular term Although this term is suppressed in , it may survive after the integration over x due to its singular nature. More precisely, consider the identities They show that when α is close to −n (n = 0, 1, . . .), the integral over x gives poles in . A thorough investigation allows us to establish that nonzero corrections to the naive limit (6.33) appear only in the components with |a| = 2. Moreover, we can explicitly calculate these corrections using Eqs. (6.34): Here we imply that ρ = 0 and that the first term in square brackets should be omitted when |ρ| = 1. We can construct the solutions which have nonzero → 0 limit exactly in one component with |a| = 2. Let us adopt the notation 1 l 1 ,l 2 ,... for the vector ρ with ρ l 1 = ρ l 2 = . . . = 1 and other components equal to zero. Then we define . These L+1 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 (ρ) −U (0) tends to a linear combination of U (b, c). Therefore, we construct the combinations that vanish in the limit → 0, for any ρ with |ρ| 3. The idea is that we can now define the solutions In particular, if L = 3, we obtain that has a finite limit → 0, while each term in the square brackets contains divergences in components with |a| = 2. Unfortunately, at L > 3 the 2 L − L(L + 1)/2 − 1 solutions J (ρ) 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 (ρ) / which vanish in the limit → 0. The combinatorics appears to be quite involved and we reserve the full consideration for our future work.
So, we consider the two groups of lines, one containing a single (r 0 = 1) line with m 0 = 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. (6.31), we define the solutions V [υ] (υ = 0, 1, . . . , L + 1) as follows where V (ρ) = V (ρ 0 ρ 1 ) are defined in Eq. (6.25). Note that when υ 2 the function V [υ] is the solution of the homogeneous system for the 'cut' sunrise integral. This is because it is a linear combination of V (ρ) with |ρ| > 1 which are all the homogeneous solutions, see footnote 4 on page 14. Again, 00 components of V [υ] all have a finite limit at D = 2: (6.42) As we explained in the previous section, the limit of components with a = 00 of each V [υ] can contain divergences. However, we have some freedom in definition which might help us to get rid of these divergences. Namely, we can construct where C υω are some coefficients with the asymptotics C υω ( ) = O( ω ). These new functions V [υ] have the same properties as V [υ] . Namely, they have the same limit of 00 component as V [υ] and are the solutions of the homogeneous system for the cut sunrise graph. Now the idea is that we can fix the coefficients C υω ( ) in such a way as to eliminate all divergences in the limit → 0. We were able to check this for a few low-loop cases 5 . Instead of running into complications connected with the basis V (0ρ 1 ) , we can pass to the basis of momenta spanned by the functions IKM[{0, 1} 1 , {υ, L + 1 − υ} m , s] defined as (6.44) Similar functions appear in Refs. [21][22][23]. For any specified L it is easy to construct the matrix 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 M has the form where [. . .] 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. [23]. The difference is that our relations explicitly depend on one parameter q 2 /m 2 while the relations discovered in Ref. [23] 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. (6.26) for K = 0, use the shift operators R, Eq. (3.6), and R k , Eq. (6.21), and the matrix 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. [23] up to at least 10 loops. However, in order to obtain the right-hand side of the same equation, we need the generalization of (6.27) to the case r 0 > 1.

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. [7,8]). Thus, we obtain the Pfaffian differential system equivalent to the Lauricella F (L) C second-order differential system. To the best of our knowledge, the Pfaffian differential system for F (L) C was previously known only for L 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. [13], we obtain the bilinear relations between the solutions of the system for opposite signs of D and µ.
We derive the operators which shift dimension and/or exponents of the propagators by ±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 α. 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, µ 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 sl(2, 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, α 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 [23] 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.
Other conditions have the form where the sum a 1 ...an runs over all words of length n with letters from the alphabet Λ, and G(a 1 , . . . , a n |x) denotes the Goncharov's polylogarithms [1,2] defined recursively by G(a 1 , a 2 , . . . , a n |x) = We note that we could have easily proceeded to more loops.