Yangian Bootstrap for Conformal Feynman Integrals

We explore the idea to bootstrap Feynman integrals using integrability. In particular, we put the recently discovered Yangian symmetry of conformal Feynman integrals to work. As a prototypical example we demonstrate that the D-dimensional box integral with generic propagator powers is completely fixed by its symmetries to be a particular linear combination of Appell hypergeometric functions. In this context the Bloch-Wigner function arises as a special Yangian invariant in 4D. The bootstrap procedure for the box integral is naturally structured in algorithmic form. We then discuss the Yangian constraints for the six-point double box integral as well as for the related hexagon. For the latter we argue that the constraints are solved by a set of generalized Lauricella functions and we comment on complications in identifying the integral as a certain linear combination of these. Finally, we elaborate on the close relation to the Mellin-Barnes technique and argue that it generates Yangian invariants as sums of residues.


I. INTRODUCTION
Theoretical predictions for particle phenomenology strongly depend on our understanding of Feynman integrals. When the number of loops and legs increases, computations quickly become intractable. Facing these problems, theorists are challenged to identify new methods to evaluate these integrals and to unveil their deeper mathematical structure. Recently a new infinite dimensional Yangian symmetry was identified for a large class of so-called fishnet Feynman graphs [1,2]. In the present paper we explore this connection between Feynman integrals and the theory of integrable models, which play a crucial role for developing analytical methods in all areas of physics. Notably, these scalar fishnet integrals furnish some of the most important building blocks of quantum field theory at any loop order. Their integrability properties can be understood through their interpretation as correlation functions of an integrable bi-scalar fishnet model, which represents an elegant reduction of deformed N = 4 super Yang-Mills theory [3]. Via this relation the integrability features of the AdS/CFT correspondence find their way to phenomenologically relevant building blocks of generic quantum field theories. Moreover, this makes connection to an alternative interpretation of fishnet integrals in terms of integrable vertex models, which was discovered by Zamolodchikov almost fourty years ago [4].
In the present paper we investigate the constraining power of the Yangian for conformal Feynman integrals. In particular, we discuss the respective constraints for the first two non-trivial cases of fishnet graphs in four dimensions. These are the completely off-shell one-loop box and the two-loop double box integral [5]: (1) Conformal symmetry allows to write these integrals as functions of 2 and 9 cross ratios, respectively. The case of the double box integral is particularly interesting since it has not been solved so far. Recently, a lot of progress was made on understanding the 7-cross-ratio limit of this integral, for which the two middle legs in momentum space (dotted lines) are put on shell [6,7]. In this limit the integral is known to be described by elliptic functions [8]. Due to its interesting relation to the double box [9] as outlined below, we also discuss the 9-variable hexagon integral:Î Also for this integral results are only known in a threeparticle on-shell case resulting in a function of 6 variables [10,11]. The Yangian symmetry employed in this paper provides the algebraic foundation of rational integrable models. Traditionally it appears as a symmetry of integrable S-matrices in two dimensions, where it typically fixes the scattering matrix completely, cf. [12]. One may thus expect similarly strong constraints for the above box, double box and hexagon integrals.
In the following we show that indeed the Yangian can be used to fix the D-dimensional box integral with generic arXiv:1912.05561v1 [hep-th] 11 Dec 2019 propagator powers. We then discuss the analogous constraints for the double box and the related hexagon integral. These constraints are formulated as systems of differential equations in the conformal cross ratios for the respective Feynman integrals. For the hexagon we argue that the Yangian constraints are solved by a large set of generalized (Srivastava-Daoust) Lauricella series in 9 variables, whose exact domains of convergence remain unclear. We discuss a recursive strategy to fix overall constants of the considered integrals by relating them to the star-triangle equation in coincidence limits of external points. Finally, we comment on the close relation of this bootstrap approach to the Mellin-Barnes technique and the common convergence issues faced for the considered six-point integrals. We close with an extended outlook pointing at various promising future directions.
Integrals built from such vertices are conformal, i.e. they transform covariantly under the differential generators J A of the conformal Lie algreba so(1, D + 1), whose densities read [13] P µ = −i∂ µ , Here the conformal dimension ∆ has to reflect the weight of the respective integral, e.g. for the above one-loop integral (3) one sets ∆ j = a j for j = 1, . . . , n. Due to their conformal symmetry, these integrals can be written in the form I n = V n φ(u 1 , . . . , u N ).
Here the prefactor V n carries the conformal weight of the integral while the variables u j denote the conformal cross ratios whose number N depends on the number n of external points. For n = 3 it is not possible to construct conformal cross ratios and hence, the above function φ is constant. This is reflected in the well known star-triangle or uniqueness relation, which holds for the conformality condition a + b + c = D: Here we have defined a = D/2 − a as well as see e.g. [14], with Γ x = Γ(x) denoting the Gamma function. At four points, the function φ in (5) becomes nontrivial due to the presence of two non-trivial conformal invariants z andz defined by Hence, conformal symmetry is no longer sufficient to completely fix the four-point function and it is natural to ask how one can further constrain the conformal fourpoint integralÎ 4 =φ(z,z)/x 2 13 x 2 24 . As noted in the introduction, the class of fishnet graphs was recently shown to feature an infinite dimensional extension of the conformal Lie algebra [1,2]. The so-called Yangian algebra is generated by the level-zero generators given in (4) and the level-one generators taking the form Here f A BC denote the dual structure constants of the conformal algebra so(1, D + 1) and the so-called evaluation parameters s j are numbers associated to each Feynman graph, cf. [1] and section H. For instance, the levelone momentum generator reads (10) Notably, these level-one generators act non-locally on the external legs of the Feynman graphs, i.e. they have a non-trivial coproduct [15]. The resulting invariance equations can be translated into partial differential equations (PDEs) in the conformal cross ratios. More explicitly, the application of the level-one generator to the integral I n yields Here, the coefficients PDE jk denote differential operators depending only on the cross ratios and we can employ conformal transformations in order to vary their prefactors independently. As shown in appendix A the above Yangian invariance condition requires that [16] PDE jk φ = 0, at least as long as we have no more than six external points. Notably, this makes connection between the Yangian symmetry of fishnet Feynman graphs and systems of differential equations for Feynman integrals which have been studied in various contexts. We will exploit these Yangian differential equations in the following.

III. BLOCH-WIGNER FUNCTION FROM YANGIAN SYMMETRY
In order to illustrate the constraining power of the Yangian algebra we start by considering the one-loop box integral for the special case of propagator weights a j = 1 in D = 4 dimensions: The Yangian invariance of the box integral translates into a system of two partial differential equations where the differential operators D j are given by Clearly, a solution to the first differential equation will not automatically solve the second equation and vice versa. It is thus natural to ask which boundary conditions lead to simultaneous solutions. We consider boundary conditions on the line z =z, which is the natural boundary of the kinematic space described by the x i , cf. appendix A for a more detailed discussion. Expanding the equations around this boundary, we find that the combination of both differential equations constrains the boundary conditions to four possible functions. In order to find the complete solution of the above system, we introduce the coordinates z 1 and z 2 as the real and imaginary part of z = z 1 + iz 2 , and we expand the equations around a generic point a = (z −z)/2i. Moreover, we introduce the function ψ(z 1 , z 2 ) = z 2φ (z 1 , z 2 ) and expand around the line z 2 = a: Note that this form is completely general (for generic a), since we do not need to consider the possibility that φ diverges for generic a. The above differential equations (14) now translate into differential equations for the coefficient functions f a,n . In particular, the full solution for the above box integral can be obtained from f a,0 via the relation ψ(z 1 , a) = f a,0 (z 1 ). Here f a,0 is essentially found by solving an ordinary third-order differential equation, which can be done straightforwardly in Mathematica. The integration constants appearing in the solution of these ordinary differential equations can be fixed, e.g. by requiring that ∂ a f a,0 (z 1 ) = f a,1 (z 1 ).  (6) is fixed by the levelzero Yangian symmetry, i.e. by invariance under the conformal Lie algebra generators (4). Similarly, the four-point integral (13) is fixed by invariance under the level-one Yangian generators (9) supplemented by permutation symmetries.
In agreement with the boundary conditions, the full solution to the Yangian constraints obtained in this way has four free parameters c j : Here we have defined Obviously, the box integral is invariant under permutations of any of its external legs. This results in functional relations forφ(z,z): Imposing these permutation symmetries onφ uniquely fixes the solution to be given by the well known Bloch-Wigner function f 4 given in (23) divided by an overall factor z −z:φ Below we will also demonstrate that the overall constant is fixed by the star-triangle integral (6) and takes the value c 4 = π 2 . This is in agreement with the results in the literature [17].
In conclusion, the four-point box integral is completely fixed by its symmetries. Note that we did not assume any boundary conditions, nor did we use an ansatz to obtain the solution. The situation resembles the star-triangle relation at three points, which is fixed by the level zero of the Yangian, i.e. by the conformal Lie algebra symmetry, see table I.

IV. PARAMETRIC BOX IN D DIMENSIONS
Next we would like to understand more generic Yangian-invariant four-point functions. A natural ex-tension is to generalize the above four-point box to D spacetime dimensions and to introduce generic propagator powers: Conformal symmetry requires that a + b + c + d = D and that the scaling dimensions entering (4) take values Note that using the star-triangle relation (6), this integral can be mapped (modulo an external propagator) to a two-loop integral with two connected three-point stars: Here the parameter e is fixed through the constraint a + b + e = D/2. Note that the propagators on the right hand side sum up to D at each of the two integration vertices. For D = 6 this integral is the natural fourpoint Yangian invariant composed of three-point vertices and with propagator weights To investigate the Yangian invariance of the above Ddimensional integral with generic propagator powers we write and note that the evaluation parameters for the Yangian generators are given in equation (H1). For conciseness we introduce the Euler operators θ j = v j ∂ vj with and the shorthand θ jk = θ j +θ k . The Yangian constraints then translate into the following parametric differential equations: Here greek parameters are given in terms of latin propagator powers and the spacetime dimension D: Importantly, equations (34) can be identified with the system of partial differential equations defining the Appell hypergeometric function F 4 of two variables u and v [18]: Here the Pochhammer symbol is given by the ratio of Gamma functions (λ, k) = Γ λ+k /Γ λ . In agreement with our findings on the special case in the previous section, it is well known that the space of solutions to the above PDEs is spanned by four functions [19]: The final steps for fixing this integral will be outlined in detail in section V. For completeness let us already note that we can employ the permutation symmetries of the box integral to completely fix the solution up to an overall constant N 4 : The overall constant can be fixed by comparison with the star-triangle integral in a coincidence limit of two external points: If we send one of the external points of the above fourpoint invariant to infinity via a conformal transformation, this result perfectly agrees with the triangle integral computed by Boos and Davydychev [20].
As already pointed out in the classic reference [19], the limit a, b, c, d → 1 of unit propagator powers in D = 4 is subtle, since in this limit the above four solutions (37)- (40) coincide. Moreover, their coefficients in (41) diverge. Careful investigation shows that the solution is given by, cf. [21]: Here h 1 (u, v) = F 4 (1, 1, 1, 1, u, v) and using the notation f ,α := ∂ α f and f ,αβ := ∂ α ∂ β f we have defined This result indeed reproduces the Bloch-Wigner function (26) as found above.

V. BOOTSTRAPPING THE BOX
In this section we demonstrate explicitly how to bootstrap the box integral with generic propagator powers from scratch. This is particularly instructive in view of the more involved examples considered in the subsequent sections. In order to solve the Yangian differential equations, we make a power series ansatz and translate the PDEs into the following set of recurrence relations for the coefficient functions g αβγγ m,n : These are straightforwardly solved using Mathematica and the solution can be brought to the following form [22], which is of course only determined up to an overall constant: We will refer to this expression as the fundamental solution. Note that in order to show that g mn solves the above recurrence equations, it is not necessary to assume that m, n are integers. We have hence found a formal solution of the PDEs for every x, y ∈ [0, 1): The solution with x = y = 0 corresponds to the solution g 1 given in (37), i.e. to the unshifted Appell function F 4 : Hence, G 1 inherits its convergence properties from F 4 , i.e. for x = y = 0 the power series (48) converges if Here, the sum in (48) effectively only extends over m, n ∈ N since the above solution (47) implies that g αβγγ m,−n = g αβγγ −m,n = 0 ∀m, n ∈ N.
In fact, this can also be observed directly by inspecting the recurrence relations (45,46). Note that if we move x or y slightly away from zero, the sum in equation (48) will extend over all of Z and diverge. We assume that a convergent series is only obtained if x and y are chosen in such a way that the series terminates at a lower or upper bound for both m and n. To achieve this we can identify all zeros of the solution (47), generalizing (51) for (x, y) = (0, 0). This limits us to the following 12 choices for (x, y): Hence, we have 12 solutions of the Yangian PDEs, which are of the form (48) and for which the series terminates.
Anticipating their interpretation we have already split these into three categories.
Let us see how this basis of solutions is related to the four functions g j=1,2,3,4 given in equations (37)-(40) of the previous section IV. Using the identity (G2) for Gamma functions, we immediately see that G 1 corresponds to the previous solution g 1 . For the case (x, y) = (1 − γ, 0) we have Now, note that which follows directly from the properties of the fundamental solution (47). We have thus found that which is related to g 2 by a constant factor that can be obtained from equation (49). In a similar fashion, we find the relations Modulo overall constants we have thus established the correspondence G j ↔ g j for j = 1, 2, 3, 4, i.e. we have identified the solutions in Region I with the four solutions discussed in section IV. So what is the meaning of the remaining eight solutions? For the case (x, y) = (0, −α), we note that which implies Comparing with the convergence condition (50), we see that in this case the series expansion is hence convergent if |u/v| + |1/v| < 1. Note that we have not found a new solution beyond the four we already encountered. The additional solutions correspond to analytic continuations of the four original series to different regions of kinematical space, see figure 1: Region II: Region III: In order to see the relation to the original solutions explicitly, we can employ the functional identity (G3) for the Appell function F 4 given in appendix G. This yields the relation Similar relations can be established for the remaining choices of (x, y) as well. We have thus derived from scratch that the solution to the Yangian PDEs is described by a linear combination of four series converging around u = v = 0: Here, we have suppressed the dependence on the cross ratios u, v as well as the parameters α, β, γ, γ for the functions G j and the coefficients c j . a. Permutation Symmetries. As anticipated in section IV, the coefficients can be constrained by employing the invariance of the integral I 4 under simultaneous permutations of the external points x j and the associated propagator powers a, b, c, d entering the solutions through the relations (35). In order to derive the consequences of permutation invariance in a compact form, we consider the generators σ x 1 = (1234) and σ x 2 = (12) of the permutation group S 4 , which act on the external legs of the Feynman diagram. These generators act on the cross ratios and parameters α, β, γ, γ as We recall the relation between the function φ(u, v) and the integral and note the functional relations that follow from invariance under σ 1 and σ 2 , respectively: The invariance under σ 1 allows to express the coefficients c 2 , c 3 and c 4 in the ansatz (64) in terms of c 1 : The additional invariance under σ 2 implies functional relations for c 1 as a function of the parameters α, β, γ, γ . The simplest way to state these relations is to note that the function is invariant under both the actions of σ 1 and σ 2 on the parameters. As a function of the parameters a, b, c, d, N 4 is hence invariant under all permutations of its arguments. This allows us to express all coefficients appearing in our ansatz (64) in terms of the coefficient N αβγγ 4 . b. Overall Constant. The above requirement does not determine the coefficient N αβγγ 4 uniquely and we employ the coincidence limit 2 → 1 of external points of the Feynman diagram in order to fix it. Applying this limit to the box integral, we find (70) implies that (u, v) → (0, 1) and we can write We thus read off that On the basis functions G j appearing in our ansatz (64), this limit acts as (we assume γ < 1) Here we have employed the identity (G4) for the Gauß hypergeometric function to end up with expressions in terms of Gamma functions. Hence, we have obtained the relation Note that the latter form makes the permutation symmetry manifest. We have thus bootstrapped the Ddimensional box integral (27) with generic propagator powers and obtained the result where for j = 1, . . . , 4 we have defined Here (x j , y j ) label the four shifts in Region I of (52) and we remind the reader that the propagator powers a, b, c, d are related to the greek parameters via (35).
It may be useful to summarize the algorithmic steps that allowed us to bootstrap the above box integral:

VI. SIX-POINT DOUBLE BOX
Being a member of the class of fishnet Feynman graphs discussed in [1,2], also the double box integral is invariant under the conformal Yangian algebra: (75) In this case conformal symmetry dictates that I 3,3 is of the formÎ with a conformally invariant functionφ 3,3 of nine cross ratios that we define as in appendix A, cf. [23]. The partial differential equations arising from the Yangian levelone symmetry read with 6 of the 15 differential operators PDE jk given by Above, for compactness, the Euler operators θ j = u j ∂ uj are packaged into Moreover, the remaining 9 Yangian differential operators PDE jk of the set (12), which annihilate φ 3,3 , are obtained from the following permutations of cross ratio labels: Notably, these permutations of the cross ratios u i leave φ 3,3 invariant: of the six external legs, respectively. An important point to note is that these permutations not only leave the integral invariant, but also the level-one momentum generator (10). Therefore, the full invariance equation (11) stays invariant under this permutation, which makes it easy to identify pairs of differential equations that are related by the corresponding functional identity. Further functional identities generalizing the invariance under the permutations above for the box integral, are listed in appendix D. Similar differential equations as given in (78) can be written down for the double box integral with generic propagator powers where we have stripped off a prefactor and conformality requires a + b + c + = D and d + e + f + = D.
Solving these differential equations in nine variables is obviously a much more involved task than for the twovariable box function. Moreover, the double box integral has less permutation symmetries than the totally symmetric cross integral. It is thus reasonable to approach this problem from a more symmetric direction and to consider a simpler situation.

VII. HEXAGON
The double box integral in D dimensions is related to the (D + 2)-dimensional hexagon via the following simple differential equation relating the respective conformally invariant functions [9] ∂ u8 φ 3,3 (u 1 , . . . , u 9 , D) = − π D/2−1 Γ φ 6 (u 1 , . . . , u 9 , D + 2), (86) which holds true for D/2 − = 1 and with φ 3,3 and φ 6 as defined in appendix C. In fact, using the expressions provided in appendix C it can be shown that the following slightly stronger equation holds true b c d e Note that here the Feynman diagrams do not represent the full integrals but rather the conformally invariant functions φ 3,3 (u 1 , . . . , u 9 ) on the left hand side and φ 6 (u 1 , . . . , u 7 , u , u 9 ) on the right hand side, respectively. The above relation implies that the hexagon integral obeys similar differential equations as the double box. In fact, we can give an argument independent of the double box, which shows that the hexagon is Yangian invariant in three and six spacetime dimensions. Firstly, in three dimensions, the hexagon is simply the fundamental Yangian-invariant vertex, similar to the box integral in four dimensions [2]. In six dimensions, Yangian-invariance follows from the following two observations: i) In [2] it was noted that six-dimensional Feynman graphs with propagator weights 2 and built from three-point vertices, are Yangian invariant (similar results hold for deformed propagator powers), ii) Using the star-triangle relation (6), the hexagon multiplied by external propagators on the left hand side can be related to a three-point graph shown on the right hand side: (89) Here we have a + b + c + d + e + f = D and the prefactor Moreover, we have redefined the cross ratios employed in the conformal parametrization above according to w 5 = u 2 u 3 u 9 , w 6 = u 2 u 3 u 4 u 5 u 9 , w 7 = u 8 u 9 , w 8 = u 1 u 2 u 3 u 8 u 9 , w 9 = u 6 u 8 u 9 . (91) These turn out to be convenient in order to write the fundamental solution to the Yangian recurrence equations in the form of a Taylor series.
Having established the Yangian symmetry and the conformal parametrization of the hexagon integral, we employ the evaluation parameters given in equation (H3) and apply the Yangian level-one generator P µ to the above expression. This yields an invariance equation of the form (11), from which we read off the 15 partial differential equations collected in appendix E. We then employ a series ansatz in terms of the cross ratios (91): Here, for convenience we have set The recurrence equations for f n1...n9 , which follow from imposing the Yangian PDEs on the above series ansatz, are listed in appendix F. Notably, these equations appear too complicated to be solved by elementary means. However, a fundamental solution to these recurrences can be obtained from the Feynman parameter representation of the hexagon integral given in (C3). Taylor-expanding this representation in the cross ratios and integrating order by order yields the expression where we use the shorthands M 1 = α + 9 k=1 n k , (95) M 2 = β 1 + n 1 + n 5 + n 8 + n 9 , M 3 = β 2 + n 2 + n 6 + n 7 + n 8 , M 4 = β 3 + n 3 + n 4 + n 5 + n 6 , M 5 = γ 1 + n 1 + n 2 + n 3 + n 5 + n 6 + n 8 , M 6 = γ 2 + n 4 + n 5 + n 6 + n 7 + n 8 + n 9 , and the greek parameters encode the propagator powers of (89): This is the analogue of the function G 1 given in (49) Here the coefficient c 1 is given by with the shorthand ρ x = Γ x Γ 1−x . As before, additional solutions of the Yangian invariance conditions can be obtained by summing over a shifted lattice with base point (x 1 , . . . , x 9 ): Obviously we have H 1 = H 0...0 . Restricting to base points for which the series terminates in all nine parameters, we find 2530 possible sets (x 1 , . . . , x 9 ), which compares to the situation of the box integral as follows: Box Hexagon Variables 2 9 Series 12 2530 As for the case of the box integral, we expect that these Yangian invariants are series representations converging in different domains, but are linked by functional relations similar to equation (G3). One may thus expect that the total number of Yangian invariants is lower than the number 2530 of series representations found above. In the algorithm outlined at the end of section V, the next step would be to classify the above sets (x 1 , . . . , x 9 ), i.e. the zeros of the fundamental solution h n1...n9 , by their kinematic region. For the box integral this can most efficiently be done by employing the shift identities listed in appendix B. However, it is not clear that similar shift identities exist for the given fundamental solution of the hexagon. This obscures the identification of a linear combination of the above series that represents the full hexagon integral. Moreover, a full analysis of the domain of convergence of all series representations seems to require a significant improvement in the current understanding of the properties of the above generalized Lauricella functions. We thus leave further steps into these directions for future work.
As argued above, the double box integral is even more involved than the hexagon considered in this section. This makes it clear that gaining full control over the hexagon bootstrap is a natural prerequisite for further investigations of the double box discussed in section VI.

VIII. RECURSIVE STRUCTURE AND OVERALL CONSTANTS
As demonstrated in the previous section VII, Yangian symmetry does not fix the considered six-point integrals completely. This underlines the need for further constraints required to eventually bootstrap these integrals. As argued in section V for the case of the box integral, also the considered six-point integrals can be recursively related to the star-triangle relation which can thus be used to e.g. fix their overall constants: While the four-point situation was already discussed in section V (see (71)), let us explain the six-point cases in more detail.
Note that we can similarly take coincidence limits of different external points leading to further equations which constrain the coefficients of Yangian invariant functions and in particular the overall constant of the integral. b. Double Box. The case of the double box integral is slightly more involved. Consider the conformal double box with parameters obeying a + b + c + = D and + e + f + g = D: We now take the coincidence limit 2 → 1 and 5 → 4 of the external points such that lim 5→4 2→1 and we use the star-triangle relation (6) on the first integral to find the box integral lim 5→4 2→1 Note that the sum of propagators in the remaining box integral gives + d + e + f = D, i.e. the integral has conformal Yangian symmetry. We can take a further coincidence limit 4 → 3 and use again the star-triangle relation to find lim 5→4 2→1 4→3 For the cross ratios the above consecutive triple coincidence limit corresponds to Hence, the limit on right hand side of (84) can be written as lim 5→4 2→1 4→3 We thus read off that lim wj →ŵj φ 3,3 (w 1 , . . . , w 9 ) = X a+b,c, X D/2−c,D/2+c−f,f .

(113)
Again, this is only the result of one possible coincidence limit and we can obtain further constraints by taking other limits.

IX. YANGIAN INVARIANTS AND MELLIN-BARNES INTEGRALS
Integrability is very constraining and if properly understood one can expect that it completely fixes physical observables through the underlying symmetry constraints. In [1,2] it was shown that certain Feynman graphs provide a means to obtain an infinite class of Yangian invariants. In the previous section we have shown that these Yangian invariants have a fine structure, i.e. there are more elementary Yangian invariants whose linear combination is selected by imposing further symmetries (e.g. permutation invariance) of the considered Feynman integrals. So what is the construction principle underlying these more elementary Yangian invariants and what is the most natural way to fix their linear combination? In order to get more insights into this, it is useful to compare the above construction to the Mellin-Barnes technique for obtaining certain Feynman integrals.
Let us discuss the box integral in more detail. Its Mellin-Barnes representation is most conveniently obtained by iteratively applying the rule to the Feynman parameter representation in appendix C and integrating out the Feynman parameter integrals.
Here the contour C extends from −i∞ to +i∞ and is chosen such that it separates the two pole series of the Gamma functions. If the intersection between the areas characterized by Re(−z) > 0 and Re(λ + z) > 0 is non-empty, the contour C can be taken to be a straight line c + iR with c being a real number such that both Gamma functions have positive real part. Carrying out the above procedure yields the following Mellin-Barnes representation for the box integral where and with N 4 as defined in equation (42). Here, κ 4 labels a point in R 2 and is again chosen such that all Gamma functions have positive real part. In figure 2 the latter region is depicted as a meshed (blue) triangle in the center of the plot and κ 4 corresponds to the orange dot, which can be moved within the fundamental triangle without changing the value of the integral. The standard method to compute integrals of the form (114) and (115) is to use Jordan's lemma in conjunction with the residue theorem [26], leading to series representations, which under favorable circumstances can be summed. To apply Jordan's lemma, one first needs to analyze in which domains of the integration space the integrand is a decreasing function. An important quantity in this context is the vector Θ which for denominator-free integrands of the form is defined as Evaluating this quantity for the integrals (114) and (115) shows that both have Θ = 0 thus corresponding to the so-called degenerate case [27]. In this case, the Gamma functions are essentially balanced and integration contours can typically be closed in multiple ways leading to series representations which are valid in different kinematic regimes. For example, the integration contour C in equation (114) can be closed via the left or right halfplane resulting in series representations which converge for |A/B| > 1 and |A/B| < 1, respectively. Similarly, we will see that multiple series representations of the box integral coexist which are analytic continuations of each other. To obtain these, we proceed by analyzing the singularity structure of the Mellin-Barnes integrand (116). The Gamma functions have poles for non-positive integer values of their arguments. In the Re(z 1 )-Re(z 2 )-plane these poles correspond to singular lines, see figure 2. For example, Γ −z1 has poles for z 1 = m with m being a non-negative integer and these are depicted as vertical solid orange lines. Similarly, the Gamma functions involving only z 2 lead to the green horizontal lines while those depending on the linear combination of both z variables correspond to the diagonal lines. An interesting point to note is that due to the special structure of the Mellin-Barnes integrand (116), finding the zeros of the fundamental solution (47) is essentially the same thing as finding the zeroth representatives of all (infinite) families of singular lines. Let us now turn to the question of how to express the Mellin-Barnes integral (115) as a sum of residues. The residues need to be computed at points where singular lines intersect but we have not yet explained which subgroups of poles should be summed to obtain a valid series representation of the original integral. However, for two-dimensional integrals there exists a simple graphical procedure which can be used to find all of these regions [27,28]. The first step consists of drawing an arbitrary cone R with vertex κ 4 in the Re(z 1 )-Re(z 2 )-plane. The cone is called compatible with the families of singular lines if each line intersects at most one side of the cone R. In figure 2 we have drawn the three cones R I,II,III (boundaries of the red areas) that are compatible with the six families of singular lines. Once such a compatible cone is found, the integral can be expressed as where the summation ranges over all intersection points that lie inside the compatible cone R i . As an example, let us compute the representation that results from summing all residues inside cone R I . Obviously, there are four families of poles in this cone, which can be parametrized as in complete agreement with table (52). Computing residues at positions z * 1 yields where (λ, k) is the Pochhammer symbol as defined in (36). The residues at positions z * 1 obviously correspond to the fundamental solution (47). Since calculating residues in cone R I is essentially trivial, we leave the computation of the other residues to the reader and merely state the final result for the Mellin-Barnes integral (116) with g i (u, v) and N 4 as defined in section IV. Note that the above result is in complete agreement with equation (41) and the Mellin-Barnes result for the triangle integral of [20]. Summing the residues in the other two cones is most conveniently done by performing the change of variables z 1 = z 1 and z 2 = z 1 + z 2 and yields similar expressions but with Appell functions depending on (u/v, 1/v) and (v/u, 1/u), respectively. The expressions obtained by summing over residues in cone R II and R III exactly agree with those obtained by applying the F 4 -identity (G3) to equation (122), thus showing that all three expressions are indeed analytic continuations of each other which converge in different kinematic regions. The above arguments now also make it clear why we chose to label the cones in exactly the same way as we labeled kinematic regions in section V: the three cones are in one-toone correspondence with the three kinematic regions in figure 1.
Finally, let us note that all four families of poles { z * 1 , z * 2 , z * 3 , z * 4 } in cone R I individually lead to a Yangian invariant quantity. This statement follows immediately from the discussion in section V and does also apply to the eight remaining families of poles. This shows that in order to obtain a Yangian invariant, one merely needs to sum over all residues originating from the same type of intersection of singular lines, see figure 2. Summing over all residues inside a given cone is apparently not required for Yangian symmetry.
Having discussed the box integral in great detail, let us now turn to the D-dimensional hexagon integral (89). Applying nine times the Mellin-Barnes identity to the Feynman parameter representation (C3) yields where and with the cross ratios w i as defined in appendix A. In complete analogy with the two-dimensional case, the vector κ 6 ∈ R 9 is defined such that all Gamma functions have positive real part.
Evaluating the vector Θ as defined in equation (118) shows that the above integral is degenerate as well, i.e. Θ = 0, so that presumably multiple series representations coexist. Since the integration space is nine-dimensional, the graphical method outlined above can no longer be applied and one needs to rely on purely algebraic methods to find all compatible cones. This, however, is left for future work. Instead, we will content ourselves with computing the analog of the fundamental solution (121). For this, we only need to find the residues at the intersection points z * 1 = (n 1 , n 2 , n 3 , n 4 , n 5 , n 6 , n 7 , n 8 , n 9 ). We obtain where the M j=1,...,6 were defined in (95) and again we used Pochhammer notation to emphasize the hypergeometric nature of the residues. More precisely, summing over all the residues z * 1 yields a multivariate hypergeometric series of Srivastava-Daoust type, see [24,25]. Picking other sets of residues leads to similar expressions with some linear combinations of n i 's replaced by others, all of them representing individual Yangian invariants.
Let us finish this section with a remark on the consideration of D-dimensional integrals with deformed propagator powers. While the deformation naively just adds another layer of complexity, it actually turns out to be a blessing in the context of Mellin-Barnes integrals as it disentangles different sets of poles which would otherwise overlap. For the box integral the latter statement becomes transparent by comparing the result (122) to the result for the undeformed box integral (43). In case of the Hexagon integral it even seems that the deformation is what makes the residue theorem applicable in the first place since in the undeformed case there exist singular points in which more than nine singular hyperplanes intersect, thus making the residues a priori no longer welldefined.

X. CONCLUSIONS AND OUTLOOK
In this paper we have demonstrated that it is possible to bootstrap Feynman integrals using their Yangian symmetry. In the case of the 2-variable box integral we have shown in full detail that the Yangian constrains the functional form of the integral to a space spanned by four Appell hypergeometric functions F 4 . Their linear combination is fixed through the integral's permutation symmetries, and the overall constant is determined by relating the integral to the star-triangle relation in a coincidence limit of external points. Hence, we have completely bootstrapped the D-dimensional conformal box integral with generic propagator powers. For the much harder 9-variable hexagon and double box integrals, we have discussed the analogous Yangian constraints, which in each case can be translated into a system of 15 differential equations in the conformal cross ratios. For the hexagon PDEs we have argued that these constraints are solved by a set of 2530 generalized Lauricella series. Due to this large number and the poor understanding of the convergence properties of these solutions, it was not possible to identify a linear combination of these series that corresponds to the integral.
These investigations suggest plenty of directions that require further understanding. Firstly, the discussion in section IX illustrates the close connection to the Mellin-Barnes technique for the computation of Feynman integrals, which in turn can be understood through the close connection between Mellin-Barnes integrals and hypergeometric functions. In our context, the Mellin-Barnes integrand is closely related to the fundamental solution of the Yangian recurrence equations. In fact, both approaches share similar problems for integrals with a larger number of variables. These are to identify the correct linear combination of series solutions and to understand their convergence properties. Already in the 2-variable case a proper convergence analysis is laborious, see e.g. [28] for the explicit discussion of the convergence of 2variable Mellin-Barnes integrals. This underlines the importance of getting better control over the mathematical properties of the often poorly understood multi-variable generalizations of hypergeometric functions. A serious convergence analysis for the 9-variable case seems indeed very hard.
Let us point out that in this paper we have observed that the Yangian differential equations for Feynman integrals can be formulated for generic spacetime dimension D, whereas the symmetry found in [2] was phrased in different but fixed spacetime dimensions D = 3, 4, 6. While here the approach with generic D emerged naturally, the case most interesting for phenomenological applications is D = 4 with unit propagator powers. It is thus natural to ask whether working directly in this limit, the considered integrals can be bootstrapped more easily. For the case of the box integral we have seen in section III that indeed in this limit the Yangian con-straints yield the solution by elementary means. This solution, however, seems less algorithmic than the bootstrap for generic propagator powers and thus less simple to generalize. The subtleties of the limit of unit propagator powers discussed in section IV show that it has various advantages to work with the deformed integral. Nevertheless it is clearly interesting to further investigate the unit-propagator bootstrap, e.g. by studying the Yangian constraints on the symbol of the respective function, similar to the approach of [11] for the hexagon integral with three massless and three massive corners.
In addition to the Yangian constraints, here we used the permutation symmetries of the box integral as well as a coincidence limit of two external legs to fix it. Aesthetically it would be more pleasing to fix an integral by integrability (alias Yangian symmetry) alone. That this is indeed in reach is suggested by the recursive structure described in section VIII. Similar to the conformal symmetry of scattering amplitudes in N = 4 super Yang-Mills theory [29], it may be possible to include this structure into the representation of the Yangian on Feynman integrals. Certainly in the case of on-shell legs, the conformal differential equations acquire inhomogeneities, and the resulting equations have been shown to yield powerful tools for the computation of Feynman integrals [30,31].
For the double box, the case with unit propagator powers and on-shell legs is known to be described by elliptic functions, whose theory in the context of Feynman integrals is still under construction [6][7][8]32]. It would be very interesting to apply the Yangian PDEs studied in this paper to an ansatz for this integral, once such an ansatz becomes available. Moreover, explicitly relating the above elliptic formalism and the hypergeometric building blocks that naturally emerge in the context of the Yangian PDEs should be a worthwhile goal.
The relation of Yangian symmetry to the PDEs (12) shows that the roots of the Yangian symmetry lie in systems of partial differential equations in the conformal cross ratios. Notably, there are various strategies to write down differential equations for Feynman integrals. In particular, the more formal approach of the recent papers [33][34][35] using the Gelfand-Kapranov-Zelevinsky systems seems closely related to ours. It would be interesting to study the systematics behind this relation and to see in how far the resulting systems of differential equations agree.
The main tool of the present paper is the Yangian Hopf algebra acting on certain Feynman graphs. Recently, similar algebraic structures were found in the context of other classes of Feynman integrals, in particular also in the context of Appell and Lauricella hypergeometric functions, see e.g. [36,37]. It should be enlightening to investigate the parallels in these approaches and to understand whether both algebraic structures coincide or coexist.
Curiously, integrability also enters the scene of conformal correlation functions from a different direction. In [38] and several follow-up works it has been shown that conformal blocks can be understood as eigenfunctions of an integrable Calogero-Sutherland Hamiltonian. There the eigenvalue equation is obtained from the conformal Casimir equation known to hold for the conformal blocks. Understanding the connection between that approach and the integrability properties of conformal correlators employed in the present paper should be instructive. A natural starting point is the box integral considered here, which can be interpreted as a correlation function in the fishnet theory of [3].
While the present paper deals with the constraints for scalar Feynman integrals, also Feynman integrals including fermions can be shown to obey a Yangian symmetry [2]. The respective diagrams are again interpreted as correlators in a generalized fishnet model, see also [39,40]. An obvious task is thus to bootstrap the simplest examples of fermionic Feynman integrals and to see how far this approach reaches for those cases.
Certainly, it is an interesting question on its own to understand the constraining power of integrability in the context of four-dimensional high energy physics. However, a more ambitious goal of this program is to develop efficient integrability methods for the computation of Feynman integrals and to understand the deeper mathematical structures underlying quantum field theory. Here it will be important to further extend the traditional integrability toolbox to the situations at hand. For the case of the yet unknown six-point integrals discussed in this paper, the present status report furnishes the groundwork for further progress into this direction.

ACKNOWLEDGMENTS
We thank Lance Dixon, Jürg Kramer, Julian Miczajka and Cristian Vergu for helpful discussions. Moreover, we are grateful to Julian Miczajka for comments and suggestions on the manuscript. DM and FL thank the Pauli Center and ETH Zürich for support and hospitality. The work of HM has been supported by the grant no. 615203 from the European Research Council under the FP7 and by the Swiss National Science Foundation through the NCCR SwissMAP. DM was supported by DFF-FNU through grant number DFF-FNU 4002-00037. The work of FL is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)-Projektnummer 363895012.

Appendix A: Details on Conformal Cross Ratios
In order to understand the degrees of freedom that remain after imposing conformal invariance, it is helpful to consider the conformal compactification of our underlying spacetime, on which the conformal group acts linearly [41]. In the case of Euclidean four-dimensional space, the conformal compactification can be realized by the light cone in R (1,5) , modulo the identification X ∼ λX. We can map this space to R 4 and vice versa employing the identifications We consider six points {X i } and constrain their coordinates as far as possible using SO(1, 5) transformations, i.e. conformal symmetry. For the first two points, it is clear that using only SO (5)  We can then determine the stabilizer of [X 1 ] by requiring that M ∈ so(1, 2) only acts as a scaling on y = (1, 0, 1) and exponentiation of the elements spanning this space. In this way, we find that we can reach the form which leaves us with SO(3) as the stabilizer of these three points. It is then straight-forward to constrain the following three points to x µ 2 = (z 1 , z 2 , 0, 0), (A6) x µ 5 = (z 3 , z 4 , z 5 , 0), (A7) x µ 6 = (z 6 , z 7 , z 8 , z 9 ), (A8) the points on the Dirac light cone (A1) follow from the relation (A2). It becomes clear from the above construction that in the case of four points, there are two degrees of freedom (compared to the 16−15 = 1 one could expect based on the dimension of the conformal group), since a stabilizer group SO(2) remains. We also note that the range of the coordinates z i is clear, since we can always pick the respective points in R 4 they represent. However, by performing rotations with angle π in R 4 , we can always enforce that In terms of the vectors [X i ] on the Dirac cone, the conformal cross ratios are given by It is helpful, to express the z-variables in terms of one set of independent cross ratios. This facilitates to check whether any other set of cross ratios is independent (by expressing it in terms of the first set) and is also a handy tool in order to derive the relations between two given sets of cross ratios. For this purpose, we consider the following set of cross ratios v 1 = u 1234 = z 2 1 + z 2 2 , v 2 = u 1432 = (z 1 − 1) 2 + z 2 2 , v 3 = u 1435 = (z 3 − 1) 2 + z 2 4 + z 2 5 , v 4 = u 1534 = z 2 3 + z 2 4 + z 2 5 , v 5 = u 1234 u 1425 = (z 1 − z 3 ) 2 + (z 2 − z 4 ) 2 + z 2 5 , (A11) v 6 = u 1436 = (z 6 − 1) 2 + z 2 7 + z 2 8 + z 2 9 , v 7 = u 1634 = z 2 6 + z 2 7 + z 2 8 + z 2 9 , v 8 = u 1234 u 1426 = (z 1 − z 6 ) 2 + (z 2 − z 7 ) 2 + z 2 8 + z 2 9 , v 9 = u 1534 u 1456 = (z 3 − z 6 ) 2 + (z 4 − z 7 ) 2 + (z 5 − z 8 ) 2 + z 2 9 .
For the first two of the above cross ratios, we also employ the notation with z = z 1 + iz 2 andz its complex conjugate. For these, we note the relations from which we read off that u, v are restricted to the domain Solving these forz i shows that these cross ratios are restricted by the relation covering the opposite domain of the Euclidean cross ratios and overlapping only along the line 4u = (1+u−v) 2 .