Yangian Bootstrap for Massive Feynman Integrals

We extend the study of the recently discovered Yangian symmetry of massive Feynman integrals and its relation to massive momentum space conformal symmetry. After proving the symmetry statements in detail at one and two loop orders, we employ the conformal and Yangian constraints to bootstrap various one-loop examples of massive Feynman integrals. In particular, we explore the interplay between Yangian symmetry and hypergeometric expressions of the considered integrals. Based on these examples we conjecture single series representations for all dual conformal one-loop integrals in D spacetime dimensions with generic massive propagators.


Introduction and Summary
Conformal symmetry plays an important role in theoretical physics. On the one hand it represents an (approximate) symmetry of many interesting models. On the other hand, it furnishes a powerful tool that puts strong constraints on a theory's observables and leads to intriguing mathematical structures. In this paper we explore an extension of conformal symmetry into two directions: its applicability to situations with masses as well as its embedding into an infinite dimensional Yangian symmetry. While there is a clear phenomenological motivation for going beyond the realm of massless particles, the extension to a conformal Yangian brings us to the theory of integrable models where one may expect that physical quantities of interest are fixed completely by the underlying symmetry.
While conformal symmetry can be studied on different levels of a given model, here we are interested in its impact on the elementary building blocks of quantum field theory, i.e. on Feynman integrals. The focus of the present paper lies on the question of how to bootstrap massive Feynman integrals by using conformal symmetry or its Yangian extension. In fact, we will discuss two instances of conformal symmetry, i.e. an 'ordinary' conformal symmetry that is here naturally formulated in momentum space and dual conformal symmetry acting on dual region momenta. In both cases we will discuss representations of the symmetry that also act on the particles' masses, which can be interpreted as extra-dimensional components of the coordinate vectors. The Yangian algebra employed here is then understood as the closure of these two conformal algebras, cf. [1][2][3].
The dual conformal symmetry of certain Feynman integrals has a long history in the massless as well as in the massive situation, see e.g. [4][5][6], and it strongly reduces the number of variables a function of interest depends on. Among the dual conformal Feynman integrals, infinite classes of diagrams of fishnet structure feature an even larger Yangian symmetry as was shown in [2,7] for the massless case and a massive version was recently found in [3]. The Yangian algebra is well known to underly rational integrable models [8][9][10], and it includes the dual conformal symmetry at the zeroth level of its infinite set of generators. In the case of two-dimensional field theories, this nonlocal symmetry typically fixes the scattering matrix completely [11]. The distinguished role of fishnet-type Feynman integrals can be understood from the fact that their conformal Yangian symmetry is inherited from planar N = 4 super Yang-Mills (SYM) theory via a particular double scaling limit of its gammadeformation. This double scaling limit yields a massless fishnet theory [12], whose correlators or scattering amplitudes are in one-to-one correspondence with individual Feynman graphs of fishnet structure. Similarly, a massive fishnet theory can be obtained from a double-scaling limit of N = 4 SYM theory on the Coulomb branch, allowing to identify massive Feynman integrals with Yangian invariant scattering amplitudes [13]. Also the massive version of the Yangian can be understood as the closure of massive dual conformal symmetry and a novel massive extension of ordinary conformal symmetry [3]. In this paper we will study the constraints of the Yangian and its (dual) conformal sub-algebras.
The idea to bootstrap Feynman integrals using their Yangian symmetry was first discussed in [14] for the examples of the massless box, hexagon and double box integrals. While the 2-variable box integral was shown to be completely fixed by its symmetries, in this first approach it was not possible to fix the linear combination of formal Yangian invariant building blocks for the 9-variable hexagon and double box integrals. Here an important step was recently made in [15], where this linear combination was determined using a multivariable extension of the Mellin-Barnes techniques, cf. [16]. In order to refine the algorithmic approach towards Feynman integrals from Yangian symmetry it would be desirable to study examples that interpolate between the above 2-and 9-variable cases. Here the recent extension of Yangian symmetry to Feynman integrals with massive propagators comes in handy, since switching on individual masses allows us to slowly increase the number of variables. In fact, initial examples of massive integrals were obtained from Yangian symmetry in [3].
In the present paper we expand on the details of the massive Yangian and conformal symmetry presented in [3]. First we prove the symmetry statements at one-and two-loop orders in detail. We then elaborate on the relation between Yangian symmetry and the massive extension of momentum space conformal symmetry. We explicitly discuss the implications of this momentum space symmetry on a few examples in Section 6. In the massless limit, the resulting constraints are precisely those that have recently been studied in the context of the momentum space conformal bootstrap (see e.g. [17][18][19][20][21][22][23]) with applications in cosmology or condensed matter physics.
We then systematically apply the bootstrap approach to massive Feynman integrals with generic propagator powers at one loop order. Since the treatment strongly depends on the choice of variables, it is useful to discuss different examples in detail, even if simpler cases can (in principle) be obtained as limits of more complex cases. In particular, we will discuss representations of the considered integrals in terms of different hypergeometric functions. The results are summarized in the following Based on the intuition gained with these examples, we finally conjecture two different series representations for the most generic massive n-point integrals at one-loop order, with propagator powers a j that obey the dual conformal constraint, i.e. they add up to the spacetime dimension D = j a j : The two conjectured series representations correspond to expressing the integral in terms of two different sets of variables The A-series represents an n-point generalization of Gauß' hypergeometric function 2 F 1 and Srivastava's triple hypergeometric series H C for 2 and 3 points, respectively. The B-series closely resembles a representation given by Aomoto [24]. 1 Moroever we note that the vtype variables are distinguished since in the case of unit propagator powers a j = 1, elegant polylogarithmic expressions for this class of integrals are known up to five points [25]. Interestingly, this family of all-mass n-gon integrals has been found to have beautiful relations to geometry, see e.g. [25][26][27]. Our analysis suggests that this beauty is closely connected to the underlying Yangian symmetry which essentially fixes these integrals completely. Note the hint at integrability in the 1998 paper [26] by Davydychev and Delbourgo, where the simplification of the integral representation for the constraint j a j = D was already called a "generalization of the so-called uniqueness formula for massless triangle diagrams". 2 In the non-dual-conformal case of unconstrained propagator powers, hypergeometric representations for these one-loop integrals were obtained in [29].
Throughout this paper we use the following notation for non-dual-conformal and dual conformal one-loop n-point integrals with propagator weights a j , respectively, which is also reflected in the subsection titles (cf. the above table of contents): In the dual conformal case denoted by •, the propagator weights obey j a j = D. We close the paper with an outlook in Section 10.

Massive Dual Conformal Symmetry
In this section we discuss the massive dual conformal symmetry of Feynman integrals. Integrals with this symmetry are particularly interesting in the context of the present paper since only these are invariant under the whole tower of Yangian generators and thus maximally constrained. The distinguished feature of these integrals is that the powers of propagators entering into an integration vertex obey the dual conformal constraint n j=1 a j = D. (2.1) After having discussed the case of one-loop integrals in large detail, we will comment on generalizations to higher loop orders.
One-Loop Integrals and Dual Conformal Transformations. We begin by discussing the dual conformal symmetry of Feynman integrals with massive external legs. As a simple example, we consider the one-loop n-gon integral with arbitrary propagator powers, .

(2.2)
The above integral is immediately invariant under four-dimensional Poincaré transformations. In order to have an object that is additionally scale and conformally invariant, we introduce a prefactor V n , Any appropriate prefactor will lead to scale invariance of the combined object under simultaneous rescalings of the x-coordinates and the masses. If the propagator weights satisfy the constraint (2.1), the function φ n is also invariant under the (D + 1)-dimensional inversion Here, we note that the indexμ runs from 1 to D + 1 and the additional component of the vectors x j is given by x D+1 j = m j . To denote an index running from 1 to D, we employ the unhatted version µ. Moreover, we use the abbreviation (2.5) For the action of the inversion map, we note that where we have applied the ordinary D-dimensional inversion to x 0 , which can be achieved by using an appropriate substitution under the integral. Correspondingly, the integral transforms as If the conformal constraint (2.1) is satisfied, the integral thus transforms by a factor. For the function φ n to be an invariant, we hence need to construct a prefactor which transforms as There are several possible choices for such a prefactor, in particular since we can obtain the respective scalings from combinations ofx 2 ij or m i . A simple choice of prefactor satisfying the above constraint is given by however, other choices of prefactors can be more convenient and we will also employ different ones below. The combination of the above inversion with D-dimensional translations yields the special conformal transformations in D + 1 dimensions, albeit with the extra-dimensional component c D+1 set to zero, since I n is not invariant under the respective translation. These transformations are generated by the conformal generator Next we note that the invariance of φ n under the above generator can be translated to an invariance statement for I n with an adapted generator, which is easy to see using the explicit prefactor given in (2.9) but holds for any prefactor satisfying (2.8). The integral I n is hence invariant under the dual conformal generators if we set the weights ∆ j equal to the propagator powers a j . That is, for J a denoting one of the above generators we have J a I n = 0. (2.14) We remind the reader that the indexμ runs from 1 to D + 1 while µ runs from 1 to D. The generators can hence also be understood as massless generators in D + 1 dimensions. Note however that in order to have invariance, we have to restrict to indices µ, ν, e.g. the integral is not invariant under translations in the mass dimension. The generators given above satisfy the conformal algebra On a massless leg j, the same representation applies with m j ≡ x D+1 j = 0. Summing up, we have found that the one-loop graph (2.2) is invariant under the dual conformal generators (2.13) provided that the propagator weights satisfy the constraint a j = D.
Higher Loop Integrals. The above invariance statement carries over to higher loop graphs if we demand that at each vertex, the joining propagator weights sum up to the spacetime dimension, Here, the variables a j denote the weights of external propagators whereas the b k correspond to internal propagators. In order to see this, consider a multi-loop integral of the form where V i andṼ i denote the set of external or internal points connected to y i . The internal propagators need to be massless in order to have dual conformal symmetry, since we are not integrating over the (D + 1)-component of the internal points, i.e. the mass. After carrying out the inversion given in (2.6), we pick up a factor of Given the above constraint, this factor cancels at each vertex.
Conformal Variables. Due to its invariance, the function φ n can be parametrized in terms of conformally invariant variables, simplifying its functional form considerably. In the massless case, the natural variables are conformal cross ratios. In the massive case, there are (at least) three natural kinds of massive conformal variables: Here we use the abbreviationx Sometimes it is useful to multiply these variables by overall constants, to add constants to them, or to consider the inverse of these variables, which may lead to a more natural form of the resulting differential equations. Clearly, the v k ij and the w kl ij are not independent of the u ij , since Therefore, in the general n-point case with all masses non-vanishing, one would try to find an independent set among the n(n − 1)/2 different u ij . However, the other cross ratios become important as soon as we consider special cases where some of the masses are set to zero. Setting k masses to zero reduces the number of degrees of freedom by k, but it leads to the vanishing of k i=1 (n − i) of the u ij . Therefore, one may need to extend the set of non-vanishing u ij by an independent subset of the v k ij and w kl ij . We note that we can and will use the freedom to select a set of independent cross ratios in order to simplify the form of the Yangian PDEs.
Independent Variables. In the general case it can be difficult to make sure that a given set of the above variables is indeed independent, and which combinations of values can be reached by choosing appropriate xμ j . In order to answer such questions systematically, we can employ the construction of the Dirac cone [30], which is also used in the construction of the conformal compactification of Minkowski or Euclidean space. To this end, we map xμ to a (D + 3)-dimensional lightlike vector with components We can map X back to x via (2.23) Note that the latter mapping is invariant under a rescaling of X and hence we consider equivalence classes of lightlike vectors X or the light-cone in projective space.
The main advantage of this approach is that conformal transformations of x correspond to linear mappings of X, which makes them much easier to treat. In our concrete case we act with transformations belonging to SO(1, D + 1), embedded in such a way that it acts trivially on the (D + 1)-component of X, which corresponds to the mass component of the spacetime vector x.
In the following we consider a configuration of 4 points xμ i and successively exhaust our freedom to employ SO(1, D + 1) transformations in order to reach a set of fixed configurations. This approach gives a different parametrization of the conformally invariant degrees of freedom of the configuration. It has the advantage that the variables we obtain are independent by construction and their range is clear. We can then check if a given set of (generalized) conformal cross ratios of the form (2.19) is indeed independent by expressing them in terms of the new variables.
We employ the notation such that the (massive) conformal symmetry SO(1, D + 1) keeps the last element of the above vector fixed. We consider the case of at least one of the external legs being massive and assume without loss of generality that m 1 = 0. The case of all masses vanishing was discussed in detail in [14]. Next, note that the vector is timelike (since we are leaving out a nonvanishing spatial component of a lightlike vector), and we can hence find a transformation in SO(1, D + 1), such that We note that this corresponds to and we have effectively used the freedom to scale our variables to set m 1 = 1 and the remaining masses are effectively measured in units of m 1 , which we make explicit in the following by using the notationm i = m i /m 1 . The stabilizer of the above vector X 1 in SO(1, D + 1) is given by the obvious SO(D + 1) and fixing the following vectors is a straight-forward exercise leading to the configuration  = 6  2  1  1  1  1  3  3  3  3  3  4  6  6  6  6  5  10  10  10  10  6  14  15  15  15  7 18 20 21 21 Clearly, after fixing n points, we have a stabilizer of SO(D + 2 − n), provided that n ≤ D + 1. The number of independent, conformally invariant variables is thus given by see also Table 1 for the case of few particles. The above derivation assumes that at least one of the masses is non-vanishing. We can conclude that, as long as one non-vanishing mass remains, we only need to subtract one for every constraint such as masses being equal or vanishing in order to find the corresponding number of degrees of freedom. In fact, for n ≥ 3, this procedure remains valid for the case of all masses vanishing, cf. Appendix A in [14].
Example: 3 Points, m 1 m 2 m 3 . As a simple example, we consider the case of three external points and three distinct, non-zero masses. We take the generalized conformal cross ratios to be From the configurations obtained above, we note (setting D = 4 for the moment), x 1 = (0, 0, 0, 0, 1), x 2 = (0, 0, 0, 0,m 2 ), The cross ratios are thus given by These expressions can in principle be employed to find out what values the triple (u, v, w) can take by solving for m i .
Here, f a bc denote the dual structure constants of the above massive, dual-conformal algebra, i.e. the spacetime indices are summed from 1 to D + 1. Introducing a free parameter y, the level-one momentum generator reads The other Yangian level-one and extra generators are listed in (A.1) and (A.2). Setting y = 1 corresponds to the choice of considering the whole algebra so(1, D + 2) for the construction of the level-one Yangian generators. Leaving out the contribution to the summation from µ = D + 1 corresponds to setting y = 0. It is interesting to note that at the one-loop level P µ is a symmetry for any value of y, as we will see below.
Let us pause here for a moment and introduce some additional notation. We denote a generator acting on the sites l through r of an n-site object as We will drop the subscripts if they are clear from context. For a 2-point Yangian level-one generator acting on legs j and k, we introduce the notation

The Symmetry at One Loop
We show that the above level-one generator is a symmetry of generic scalar n-point Feynman integrals at one-loop order with massive propagators, The propagator powers a j and the spacetime dimension D are arbitrary, and we use the notation x µ jk = x µ j − x µ k . These integrals are invariant under all permutations of the external legs which are accompanied by the respective permutations of the propagator weights a j and the masses m j . It turns out that already the integrand is invariant, which implies that J a I n = 0. (3.8) Before we go on to prove the above result, let us discuss the implications of the permutation symmetry of the n-gon integral I n . Concretely, we consider the Yangian level-one generator with the evaluation parameters s (n) j given by for the one-loop integral we consider. We note that we can split up the generator as n−1 J a n−1 + s (n) n J a n . (3.11) Here, the restriction to the bi-local part corresponds to leaving out the local contribution governed by the evaluation parameters. Next, we consider the permutation P = (n − 1, n), along with the respective permutations of the weights a j , and note that P −1 J a (1,n),n P = J a (1,n−2),n + 1 2 f a bc J c (1,n−2) J b (n−1,n) − J a (n−1,n),n bi-local +s (n) n−1 J a n +s (n) n J a n−1 , where s (n) n−1 = 1 2 a n − 1 2 a (1,n−2) , s (n) n = − 1 2 a n−1 − 1 2 a (1,n−2) , (3.13) s (n) n−1 = 1 2 a n−1 − 1 2 a (1,n−2) ,s (n) n = − 1 2 a n − 1 2 a (1,n−2) . (3.14) Due to the permutation symmetry of the integral I n , the transformed generator is a symmetry as well, and hence so is the difference of the generators J a (1,n) − P −1 J a (1,n) P = f a bc J c n−1 J b n − a n−1 J a n + a n J a n−1 = 2 J a n−1,n . The n-point invariance thus implies a two-point invariance which, due to the full permutation invariance of the integral, holds for any given pair of points. Conversely, it is easy to see that the evaluation parameters are designed in such a way that Two-Point Level-One Invariance. Finally, let us prove the two-point level-one invariance of the integrand (3.7). We begin by considering the level-one momentum density at y = 0 Plugging the level-zero densities (2.13) into the above expression yields where Since the integrand of (3.7) is factorized and the level-one density (3.17) only contains derivatives with respect to points j and k, it suffices to consider the action of generator density on a product of two propagators. Using the abbreviationx 2 0j = x 2 0j + m 2 j we find Carrying out the contractions yields and consequently Checking the two-point invariance of the the integrand (3.7) under the remaining level-one and level-one extra generators is a straightforward exercise. For y = 0 and a dual-conformal integrand (2.1) there is of course no need for further checks as the algebra relations already guarantee the invariance under the remaining level-one generators. However, let us stress that level-one invariance as well as level-one extra invariance of the generic scalar n-point Feynman integrals (3.7) hold no matter whether the conformal constraint (2.1) is satisfied or not. This observation will later on allow us to also consider non-dual-conformal integrals which do not have full level-zero symmetry.
In summary, we have found that one-loop integrals (and in fact already the integrands) are invariant under Yangian level-one generators,  Internal Mass. The last finding can be puzzling, since e.g. the generator P µ extra involves the densities of the generator L µD+1 , which mixes the mass with the other dimensions and is hence itself not a symmetry of the Feynman integral. We can illustrate the significance of the contribution P µ extra by allowing the internal mass m 0 to be non-vanishing. We are thus considering the Feynman integral (3.25) Note that this integral arises from the one with vanishing m 0 by shifting the external masses, We thus find that the level-one generators act on the shifted Feynman integral as The above conjugation is evaluated as where, [A, B] (n) denotes the n-fold commutator These commutators are easy to evaluate by employing the commutation relations (2.15) and (3.1). Specifying to the level-one momentum generator, we note that P µ only commutes with P D+1 if we set y = 1 in (3.3) and thus include the contribution P µ extra summing over the mass component as well. We conclude that for non-vanishing m 0 the operator P is only a symmetry if y = 1.

The Symmetry at Two Loops
The invariance of two-loop graphs can be derived from the invariance of the constituting one-loop graphs. Concretely, we consider a two-loop graph that arises from joining two one-loop graphs with l + 1 and r + 1 legs, respectively, thus having n = l + r legs in total: We have noted in the above discussion that at one-loop order the diagram need not be invariant under the underlying dual conformal symmetry in order to have level-one invariance. For the two-loop discussion, however, level-zero invariance is more critical and we will set y = 0 in the following, since the Lorentz generator L D+1µ is not a symmetry of the Feynman diagrams we consider, and thus the extra generators J extra will not generate separate symmetries at two loops, cf. (3.24).
Dual Conformal Case. We would like to show that there is a set of 2-loop evaluation parameters s (2,n) j such that the above generator (3.9) becomes a symmetry of this diagram. Here we denote the -loop evaluation parameters by s ( ,n) j . To this end, note that we can split up a generic level-one generator as follows: The terms in the last line are due to the differences of the evaluation parameters for the generators acting on the diagrams containing 1 or 2 loops, respectively. When acting on the level-zero invariant I n , we note that where the dual Coxeter number c arises from the contraction Consequently, we have Here, we have used that I (2) n is invariant under the partial level-one generators acting on the first l and last r legs, respectively, since already the integrands of the constituent one-loop graphs are invariant.
We can then take the above equation as a definition of the evaluation parameters for the Yangian level-one generators at the two-loop level. Concretely, this gives The dual Coxeter number can be inferred by commuting the constituents of the level-one momentum generator (3.4). This gives and consequently Non-Dual-Conformal Case. At the one-loop level, we noticed that all level-one generators are symmetries even if the conformal constraints (2.16) are not satisfied. This finding partially persists at the two-loop level, where we restrict ourselves to the level-one momentum generator, again setting y = 0. For this generator, many aspects of the above discussion are still valid. The only difference is that I (2) n is no longer annihilated by the dilatation generator D. Instead, we note that which only vanishes if the conformal constraints are satisfied. We can then modify (3.34) to yield from which we read off the evaluation parameters We note that we can always adapt the evaluation parameters by employing a shift s k → s k + C, which acts as on the level-one generator. Since P µ is a symmetry generator itself, we can omit the last term when discussing the respective level-one generator. In this way, we obtain the evaluation parameters These evaluation parameters can be stated in terms of the following rule: Pick an arbitrary leg as leg one and set Then move clockwise around the diagram and update the next parameter as if legs k and k + 1 are attached to the same vertex, or as if the vertices of legs k and k + 1 are connected by an internal propagator with weight b. This rule was also stated in [3]. We note that it applies to all cases we discussed above and we have hence omitted the explicit reference to the loop number.  Summary and Higher Loops. The above findings at one-and two-loop orders are summarized in Table 2. Here the statements at higher loop orders reflect the following conjecture formulated in [3]: Feynman graphs cut along a closed contour from one of the three regular tilings of the plane with massless internal propagators have full Yangian symmetry in the dual conformal case, or only P symmetry in the non-dual-conformal case, respectively; external propagators can be massive or massless. This conjecture is motivated by the fact that in the massless limit, these are the classes of Feynman diagrams that are known to enjoy Yangian symmetry [2]. It is further supported by numerical evidence obtained as follows.
The level-one momentum generator P was applied to the Feynman parametrization of examples of Feynman graphs in the respective categories. The resulting expression was then integrated numerically in Mathematica and it was compared whether the given uncertainty neighborhood of the result includes zero. It would certainly be desirable to find an analytic proof of this conjecture on massive higher loop integrals similar to the one in the massless case [2], or to extend the numerical tests with more advanced numerical techniques.

Two-Point Yangian Invariants
Let us discuss some subtleties specific to two points. Consider a two-point invariant under the level-one symmetry which obeys J a I 2 = 0, (3.49) where at two points we have The one-loop integrals are actually invariant under the level-zero symmetry if the dual conformal constraint a 1 +a 2 = D holds, i.e. they obey (2.14) This implies that we can write From (3.10) we can read off the one-loop evaluation parameters s Hence, for the dual Coxeter number c = D as given in (3.38) this equation is trivial if the dual conformal constraint a 1 + a 2 = D holds and thus yields no non-trival constraints that could help to determine the integral.
Nevertheless, there are two ways to obtain non-trivial constraints on a two-point invariant from the Yangian level-one generators: • Firstly, at one loop order we can employ the extra symmetry J extra , see (3.24) and the two-point examples in Section 8.
• Secondly, we can give up on (parts of) the level-zero symmetry, e.g. the special conformal symmetry K µ , which amounts to relaxing the dual conformal condition j a j = D.
Examples of how the resulting constraints can be used are discussed in Section 7.

Recursions from P D+1
As seen above, there is an extra contribution to the level-one generators containing the D + 1-components of the generator densities, in particular the D + 1 momentum density Recursions from P D+1 . While these generator densities feature in the level-one symmetry of one-loop diagrams, they do not combine to form level-zero symmetry generators. It is well known that the resulting differential equations in the mass variables are very helpful in solving Feynman integrals [31]. As an example, acting on the one-loop integral (2.2), .., a n ] = 2ia j m j I n [a 1 , ..., a j + 1, ..., a n ], (3.59) i.e. the set of all diagrams transforms covariantly under the action of the D + 1 momentum generator density. Hence, after the kinematic dependence of a diagram has been reduced to a linear combination of basis functions using the level-zero and level-one symmetries, these covariance equations can be used to constrain the propagator weight dependence of the expansion coefficients. Since there is an independent covariance equation for every mass m j , these constraints are most powerful for integrals with all propagators massive, where the dependence on all propagator powers a k is fixed. Otherwise, only the dependence on the propagator powers of massive propagators are constrained. Below, we use these identities to fix the propagator weight dependence of the non-dualconformal all-mass two-point one-loop integral, such that the only remaining freedom is given by numerical prefactors. This information can then be transported to all other two-point cases by taking massless and conformal limits.
From Conformal to Non-Conformal. Another important application of the D + 1 momentum densities is to generate non-dual-conformal integrals from dual-conformal ones. As we argued in Section 2, Feynman diagrams are dual conformal if for any vertex the sum of all propagator weights gives the spacetime dimension D. Acting on the respective Feynman integral with a D + 1 momentum density raises the propagator weight of a leg and breaks the conformal condition at the corresponding vertex. Repeated action allows to extract the integral at an arbitrary integer distance to conformality from the knowledge of the dual-conformal integral.
As an example, consider the case of the two-point one-loop integral with one nonvanishing mass. In section 8.1 and section 7.3 we calculate both the conformal as well as the non-conformal version of this integral respectively. They are given by and respectively. The precise relation between the two results is given by (3.62) From the non-dual-conformal result, we find To be concrete, taking n = 1, we have in agreement with the derivative of the conformal integral. Hence, knowledge of the dual-conformal case is actually enough to derive the results of many non-dual-conformal cases, at least the ones with integer deviations from conformality. These cases are also the ones important for phenomenological applications. Still, having the results of arbitrarily many non-dual-conformal integrals does not necessarily allow to derive the continuous dependence on propagator weights and spacetime-dimension.

Yangian Bootstrap in a Nutshell
In this section we discuss the Yangian bootstrap algorithm introduced in [14] and extend it to the massive situation. After a brief explanation of how to extract the Yangian PDEs from the invariance equations, we illustrate the algorithm by means of a simple example. Similar steps can be taken for the integrals discussed in the subsequent sections, but below we will refrain from giving all details and rather highlight some key elements.

Extracting Yangian PDEs
Above we have argued that certain classes of (massive) Feynman integrals are invariant under the Yangian algebra. In order to evaluate the constraints from Yangian symmetry on a given integral, it is useful to translate the Yangian invariance equations into differential equations for the function φ n that depends only on a reduced set of variables which takes into account Poincaré, scaling or full (dual) conformal symmetry. Doing this requires two steps: First, we bring the PDEs to a canonical form, e.g. for the level-one momentum generator we have Here the form of the vectors a µ jk depends on the particular choice variables. As a second step, we may exploit the spacetime symmetry (e.g. conformal symmetry) of the integral to argue that the vectors a µ jk are in fact independent, 4 which implies the following set of partial differential equations: Similarly we can obtain a separate set of one-loop PDEs, which follow from the extra symmetries, see (3.24). Determining the respective integral then boils down to solving the resulting set of PDEs and to using additional input (e.g. permutation symmetry in the external legs) to identify a unique invariant that corresponds to the integral under study. It is instructive to discuss an example.
Example: 3 Points, m 1 m 2 m 3 (Dual Conformal). As an example consider the conformal three-point integral with three non-vanishing masses and the kinematic variables given in (2.31), This integral is fully bootstrapped in Section 8.5. We set I 3 = V 3 φ 3 (u, v, w) and choose the prefactor V 3 that carries the scaling weight according to Now we would like to evaluate the level-one momentum invariance which takes the form Here the symbols PDE jk represent differential operators in the conformal varibles u, v and w.
In the following we will argue that the vectors are independent such that we can read off three partial differential equations in the cross ratios: 5 In order to see the independence of the above vectors, we employ the configurations derived in Section 2, Since for the case of equal masses the number of variables is reduced, we assume that the masses are different and ordered according to Now we can investigate, whether the vectors a µ jk given in (4.7) are independent after conformal transformations. In the three-variable case, the invariance equation (4.6) takes the form Setting the parameter c µ of the special conformal transformation to zero, yields On the other hand, taking the c 1 , c 4 -derivative and then setting c µ to zero gives Hence, for m 2 > 1, we can conclude that PDE 13 and PDE 23 independently annihilate φ 3 . Finally, taking the c 4 -derivative and then setting c µ to 0 yields (after using that PDE 13 and Hence, also the differential operator PDE 12 annihilates φ 3 .
5 For explicit PDEs for the considered integral see Section 8.5.

Algorithmic Solution of PDEs
In [14] an algorithmic solution of the Yangian PDEs for massless Feynman integrals was presented. The same steps can also be applied to massive Feynman integrals: 1. Translate the symmetry PDEs (4.2) into recurrence equations.
2. Find a fundamental solution of the recurrence equations.
3. Find all zeros of the fundamental solution yielding a solution basis.
4. Classify these zeros by their kinematic region. 5. In a given kinematic region, use further constraints to fix the linear combination of basis functions.
For conciseness we will not go through all these steps for all of the examples discussed in the subsequent sections. We explicitly illustrate the above algorithm on the following pedagogical example.

2 Points
Consider the two-point integral where we do not impose any (dual conformal) constraint on the parameters a 1 , a 2 and D.
We choose to write The integral is trivially annihilated by P and L but not by K and D since we do not have level-zero symmetry under K and D, cf. Section 3.3 . Hence, P and P extra do not yield any non-trivial PDEs. However, invariance under K and D give rise to the same equation, namely to Euler's hypergeometric differential equation where we have introduced the parameters To convert this Yangian PDE into a recurrence equation we make the series ansatz Here, a priori the sum runs over all integer numbers and we even allow for a non-integer shift x ∈ [0, 1) of the (here one-dimensional) summation lattice. Then the hypergeometric PDE translates into the following recurrence equation for the coefficient functions g n : 0 = (n + α)(n + β)g n − (n + 1)(n + γ)g n+1 . (4.20) Modulo an overall constant this equation has a unique fundamental solution where we make use of the Pochhammer symbol While the above representation g αβγ n of the fundamental solution in terms of Pochhammer symbols may be better behaved in certain limits, we can alternatively express the above series via f αβγ where for integer n we can write Here the overall constant C depends on the parameters α, β, γ and the base point x. Hence, for every x ∈ [0, 1) the following series represents a formal solution of the hypergeometric PDE: The above series G αβγ x (u) terminates for the following four choices of x, which we take as a condition for convergence: These four choices for x correspond to the zeros of the fundamental solution (4.23), and anticipating their interpretation we have classified them according to two different kinematic regions. Moreover, the fundamental solution (4.23), which is symmetric in the first two parameters, obeys the following shift identities (plus the identities with α and β exchanged): (4.28) Note that these identities and similar identities for the more complicated functions considered in the rest of the paper get additional prefactors when formulated in terms of the g αβγ n of (4.23).
Evaluating the series defined by (4.25) explicitly we find the following expressions in terms of Gauß' hypergeometric function 2 F 1 : The above algorithmic procedure generalizes to the more involved examples discussed in the subsequent sections. In this simple case, however, the two functions G I,1 and G I,2 in Region I are also found when solving the Yangian PDE (4.17) directly in Mathematica: Finally, we can use the P D+1 recursion from Section 3.4 to determine the a 1 -dependence of the prefactors c j . Here, (3.59) reads In this particular case, this constraint could be solved using identities of 2 F 1 , but in anticipation of the more complicated cases in the following sections, we instead expand both sides of (4.34) in a series and compare term by term, Solving the resulting difference equations for c 1 and c 2 , we find where e 1 and e 2 are constant in a 1 but may still depend on a 2 and D.
To fix the full dependence of the prefactors, one needs to start from the more symmetric two-mass integral given in (7.45), which we bootstrap in Section 7.5. Taking the m 2 → 0 limit, the answer agrees with the a 1 dependence we derived here and the full result given in [32] modulo a conventional overall factor: Conformal Limit. Due to distinguished role of integrals with dual conformal symmetry, let us now consider the dual conformal limit in which we expect to find the below result (8.2). Setting D = a 1 + a 2 + 2 , we have Accordingly, in the limit → 0 (4.33) takes the form where A ± = (a 1 ± a 2 )/2 + 1. Since c 1 given in (4.37) is of order , we conclude lim D→a 1 +a 2 in agreement with (8.2).

A Change of Perspective: Massive Conformal Symmetry in Momentum Space
In this section we look at Yangian symmetry in region momentum space from the momentum space point of view. We begin by deriving the momentum space representation of the dual level-one momentum generator, thereby introducing a novel massive generalization of momentum space conformal symmetry, cf. [3]. After verifying the algebraic consistency of this representation, we explore the idea to bootstrap Feynman integrals in momentum space by utilizing the newly gained insights.
Translating Generators. A natural question is whether the massive dual conformal Yangian symmetry can be understood as the closure of the massive dual conformal Lie algebra symmetry and an ordinary (massive) conformal symmetry, similar to the situation in the massless case [1]. To address this question, we focus on the level-one momentum generator in dual space where and rewrite it in terms of momentum variables. The latter are related to the dual variables in the following way: 3) The mass variables, on the contrary, stay untouched. Note that momentum conservation implies the identification x µ n+1 = x µ 1 which will be implicitly assumed henceforth. Inverting the above relation yields which shows that the dual variables are in fact region momenta, thus explaining the arbitrary reference point. In order to rewrite the generator (5.1) as a generator in momentum space, we furthermore need to express the derivatives in equation (5.1) in terms of derivatives in momentum space. To this end, we apply the chain rule and obtain Substituting the expressions (5.4) and (5.5) into equation (5.1) and simplifying the result yields Note that momentum conservation always allows us to eliminate one momentum in favor of the remaining n − 1 momenta. The massive n-gon integrals for example can be expressed as so that the n-th momentum does not appear. This choice turns out to be very convenient since it allows us to to drop all terms containing a derivative with respect to p n in equation (5.6). Finally, we recall that level-one momentum invariance requires the scaling dimensions ∆ i to be equal to a i while the evaluation parameters s j have to be chosen according to (3.10).
For the combination appearing in equation (5.6) this implies Algebra. With the above results at our disposal, we can now define the massive representation of the conformal algebra in momentum space as A straightforward computation of the commutation relations yields that these generators indeed satisfy the same algebra relations as the massless generators, i.e.
where the action on an n-fold tensor product space reads This confirms the statement that (5.10) is merely a different representation.
Consistency Check: One Loop Invariance. In order to check the consistency of equation (5.6), let us verify that the right-hand side indeed annihilates the massive n-gons (5.8).
For simplicity, we set y = 0. Denoting the integrand of (5.8) as i n we find withK µ k denoting the massless special conformal generator acting on the loop momentum. Due to the fact thatK µ k is itself a total derivative, the above expression vanishes when integrated over.

Momentum Space Conformal Bootstrap
In this section we explore the idea to bootstrap Feynman integrals in momentum space rather than in dual space. The motivation to do so is twofold. First, this idea seems natural as the ubiquitous (dual) level-one momentum generator has been identified as the special conformal generator in momentum space which is local and thus easier to handle. Second, such an analysis bridges the gap between the Yangian bootstrap approach and the study of conformal constraints in momentum space pursued e.g. in [18,20].

3 Points, 000: Appell F 4
In order to start with an example that is in fact completely fixed by momentum space conformal symmetry, we begin by considering the non-dual-conformal massless star integral with three external points While its dual conformal cousin with a 1 + a 2 + a 3 = D is uniquely fixed by the star-triangle relation, see e.g. [33], no such statement holds for the non-dual-conformal version of the integral. We therefore employ the conformal momentum space symmetry from above to constrain the function. To do so, we first express the star integral in terms of momenta by using equation (5.3) and obtain Next, we employ the scaling equation to justify the following ansatz: Eliminating p 3 from the ansatz by using momentum conservation and acting on it withK µ m=0 and∆ i as specified in (5.9),i.e.∆ i = a i + a i+1 , yields where Here, the Greek parameters are related to the propagator powers and dimension in the following way: Since p 1 and p 2 can be freely varied, both equations have to be fulfilled independently. We recognize these partial differntial equations as the system defining the Appell function F 4 , see also [18,20] for similar discussions of the conformal momentum space constraints. This comes as no surprise since the triangle integral has been shown to evaluate to a linear combination of four F 4 functions more than 30 years ago [32]. Furthermore, the triangle can be obained from the box integral by sending one of the external points to infinity [5]. The latter was recently computed in [14] by utilizing the Yangian bootstrap approach and we can use the exact same techniques to reproduce the result stated in [32]. Here, we only give a brief summary of the necessary steps. For more details the reader is referred to [14].
In order to solve the partial differential equations from above, we make a power series ansatz Acting with the PDEs (6.6) on this ansatz yields recurrence relations for the coefficients f αβγγ kn which can straightforwardly be solved, for example, by using Mathematica Note that this expression leads to a (formal) solution of the PDEs for any value of x and y in (6.9). However, for generic values the series will most likely be divergent for any value of u and v because the sum extends over a shifted Z-lattice. Only if x and y are chosen in such a way that the series terminates at a lower or upper bound for both k and n the series has a chance of being convergent. A careful investigation of the zeros of fundamental solution (6.10) shows that there are 12 choices for (x, y) for which the power series terminates and converges. However, only four of them converge in a region around the origin in the u-v-plane which is the region that we want to focus on here. These are In the final step, we employ the permutation symmetries of the triangle 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 relation (6.2): 14) The result (6.13) can also be shown to be in full agreement with the Feynman parametrization of the function φ(u, v) reading

3 Points, m 1 00
Let us now consider the same integral with one massive leg: Again, we utilize the scaling equation to justify an ansatz of the following form: Here, the signs have been chosen for later convenience. Eliminating p 3 from the ansatz using momentum conservation and acting on the integral withK µ m for∆ i = a i + a i+1 yields K µ m I m 1 00 where Since the number of independent momenta has not increased compared to the massless case, we again find two partial differential equations which need to be satisfied independently. However, as the mass introduces an additional degree of freedom, the function (6.17) now depends on three scale invariant variables instead of two. Hence, the number of PDEs is not sufficient to fully constrain the function. To make matters worse, the Pμ extra symmetry equation is trivially fulfilled and does therefore not yield any additional constraints.
To be more explicit, we make the series ansatz The function G xyz (u, v, w) solves the above differential equations for with an unfixed function c k+l+n that depends on the sum of the three summation indices. This function is not fixed by momentum space conformal symmetry only. A natural course of action to generate further PDEs is to consider the full set of levelone generators in dual region momentum space instead of only the level-one momentum generator, cf. Table 2. In fact, since the level-zero algebra in dual space is partially broken, it is clear that considering just the level-one momentum generator is no longer sufficient. Considering the full set of dual Yangian level-one generators can of course also be done in momentum space. However, since the level-one generators are scattered among different levels in momentum space, for example, the generator K µ in x-space corresponds to the trilocal level-two momentum generator in p-space, we prefer to work in region momentum space which puts all generators on an equal footing. We will come back to the above integral I m 1 00 3 in Section 7.6 where the function c k+l+n is constrained using the remaining Yangian level-one generators, cf. (7.69):

3 Points, 2 Loops, m 1 00
Consider now the two-loop integral Acting on the above ansatz withK µ m and ∆ 1 = a 1 + a 2 + b 0 − D/2 , (6.26) as follows from the general rule∆ i = ∆ i +∆ i+1 +2s i −2s i+1 , we find exactly the same system of partial differential equations as in the one-loop case (6.19). In fact, as in the one-loop case, those equations do only depend on a 2 , a 3 and D and not on a 1 . TheK µ m equations are solved by the fundamental solution (6.21), with the yet to be determined function c k+l+n . For the one-loop integral this function is fixed by the K equations. Hence, the one-and two-loop integrals (6.16) and (6.23) only differ in the function c k+l+n .

Changing Back: Non-Dual-Conformal Integrals
While we will study dual conformal integrals in the following Section 8, here we consider integrals without imposing any constraints on the propagator powers a j . In particular, this means that these integrals are not invariant under the dual conformal generator K µ at level zero. Despite the absence of dual conformal symmetry, we will see that in comparison with the momentum space conformal symmetry discussed in the previous Section 6, the Yangian level-one generators yield additional constraints for one-loop integrals as discussed in Section 3.1. In particular, we focus on the interplay between the Yangian differential equations and their solutions via hypergeometric functions. We also discuss relations between different cases. Even if integrals with more masses and parameters can in principle be reduced to simpler examples, it is instructive to explicitly discuss different cases with increasing complexity. A useful variable in a more complex example with two masses may be given by u = x 2 12 /m 1 m 2 which in the limit m 2 → 0 diverges, thus obscuring the reduction of the integral to a simpler case.

1 Point, m 1 : Rational
As the simplest example consider the tadpole integral Using a single spacetime point one cannot form a translationally and scaling invariant variable. Hence, the integral is pure weight, i.e.
To fix the propagator weight dependence of the constant c a 1 , we act with P D+1 as described in 3.4, implying Finally, we evaluate the integral numerically at a single point to fix the overall constant and find

2 Points, 00: Rational
Also for two points and two vanishing masses there is no scaling-invariant variable and the one-loop integral collapses into a trivial propagator. This is also known as the group relation, see e.g. [34]: with the constant In Section 8.1 we discuss a similar situation in the dual conformal case with a 1 + a 2 = D, where a massless and a massive propagator are fused into a massive propagator.

2 Points, m 1 0: Gauß 2 F 1
Note that the two-point integral with one mass was explicitly discussed in Section 4.3 for the variable −m 2 1 /x 2 12 as an example to illustrate the bootstrap algorithm. When comparing to limiting cases of other integrals, it can be useful to consider different choices of variables, e.g. we can invert the variable used in Section 4.3 and write the Yangian PDE obtained from invariance under the level-one special conformal generator K reads and is solved by This result can be compared to the below three-point result (7.75) in the coincidence limit of points 2 and 3 in Section 4.3 which yields

2 Points, m 1 m 2 : Kampé de Fériet
Also for two non-vanishing masses different choices of variables lead to different types of functions in which the considered two-point integral can be expressed. For a nice solution in terms of a single Appell F 1 series see [35]. 7 As a first example we choose our variabels to be such that we can write For convenience we set (7.14) Making the series ansatz we can solve the level-one P extra and K y=0 equations to find the fundamental solution We can alternatively express the above series via where for a certain constant C = C(α, β, γ, x, y) and for integer k and n we have f k+x,n+y = C g k+x,n+y .
We find 13 doublets (x, y) that correspond to zeros of the fundamental solution f kn in k and n, which make the above series terminate at an upper or lower bound. Only the following two of these doublets give rise to effective variables u and v: In terms of the two basis functions we can thus make the ansatz φ = c 1 G 00 + c 2 G 0,−α−β . (7.21) Using the covariance under the action of P D+1 , see Section 3.4, the prefactors are determined to be where e 1 and e 2 are fixed by two random numerical configurations to e 1 = π D/2 e 2 = 0. (7.23) This reproduces the generalized Kampé de Fériet hypergeometric function in [29] (modulo conventions).
Equal-Mass Limit. Consider the limit m 2 → m 1 where v → 0. Hence, for lim v→0 v n = δ 0n and assuming α + β < 0, we find The coefficient c 1 is fixed by the below limit The given result agrees with the expression of [32] for the equal-mass two-point integral.
One-Point, One-Mass Limit. Consider now the combined limit m 2 → m 1 , Comparing to the tadpole result of Section 7.1 for a 1 → a 1 + a 2 which reads we can read off (7.25).

2 Points
Let us discuss a second choice of kinematic variables which is more symmetric than (7.12), and which will lead to a solution expressed in terms of Appell hypergeometric functions F 4 . This result will be useful to compare to various limiting cases that were expressed in terms of Gauß' 2 F 1 . In fact, this example will show that more symmetry in the choice of variables does not necessarily lead to a simpler solution in the sense that the resulting expression will be a linear combination of hypergeometric functions rather than a single hypergeometric series as in the previous Section 7. 4. We now write where we choose the symmetric variables Moreover, we define the following abbreviations Linear combinations of the Yangian PDEs obtained from P and K invariance yield the system of two differential equations defining the Appell hypergeometric function F 4 , cf. (6.7): Hence, for small u, v the solution is a linear combination, cf. (6.11), ; u, v , (7.35) Permutation Symmetry. Let us now use further input to constrain the coefficients c αβγγ j . We note that the considered two-point integral is invariant under the permutation (x 1 , m 1 , a 1 ) ↔ (x 2 , m 2 , a 2 ), which translates into (u, γ) ↔ (v, γ ). Under this map we have such that we conclude that permutation symmetry implies the following constraints Acting with P D+1 on the general solution (7.33) and using the linear independence of the g i leads to a set of recurrence relations for the c i in the propagator weights which are solved by Here the e i are constant complex numbers independent of a 1 , a 2 and D and the above contraints (7.39) from permutation invariance imply e 2 = e 3 . Using random points of numerical data we finally fix e 1 = π D/2 , e 2 = e 3 = (−1) −D/2 π D/2 , e 4 = 0. (7.42) Hence, in total we obtain the full expression for the two-mass integral Indeed, the result given in [32] agrees with the above expression obtained from bootstrap (modulo phases due to a different sign convention in the propagator).
One-Mass Limit. For m 2 → 0 we have v → 0, such that due to the reduction formula Taking in addition the conformal limit D → a 1 + a 2 of this expression we have c 1 → 0 and This agrees with the below expression (8.2).
Conformal Limit. We would like to take the limit D → a 1 + a 2 of the above full twomass expression in order to arrive at the dual conformal case presented in terms of associated Legendre functions in (8.12) or in terms of hypergeometric functions in (8.31). In this limit we have where we define a ± = 1 2 (a 1 ± a 2 ). The four basis solutions become whereas the coefficients turn into Numerically, we find the interesting relation 8 which implies that the conformal two-point function becomes ;ũ + (a 1 ↔ a 2 ) , (7.52) To convert this result into the expression (8.13) in terms of associated Legendre functions P and Q that we obtain from bootstrap in the subsequent Section 8, we use the identities [36] as well as which are valid for u < 0 only. Remarkably, these identities imply that the expression (7.52) for the dual conformal two-point integral finally collapses to (see (8.13)): Next we study the one-mass triangle integral We write the three-point integral in terms of a scale invariant function of three arguments: This integral has the full level-one symmetry but no special conformal level-zero symmetry since we do not impose the dual conformal constraint on the propagator powers. Imposing level-one momentum and special conformal symmetry on the above ansatz, we can read off the coefficients of the vectors x µ j . Note that in the case of the special conformal levelone generator K, in order to turn these coeffients into functions of u, v, w modulo overall coefficients, we need to impose the P equations which makes the coefficients of x 2 j with j = 1, 2, 3 vanish.
The resulting PDEs can be turned into recurrence equations for the coefficients g kln in the series ansatz G α 1 α 2 α 3 γ 1 γ 2 xyz = k∈x+Z l∈y+Z n∈z+Z g kln u k v l w n = u x v y w z k∈Z l∈Z n∈Z g k+x,l+y,n+z u k v l w n . (7.59) For convenience we introduce the parameters as well as the depend parameter Modulo an unconstrained overall constant, the recurrence equations are solved by the unique fundamental solution For a fixed constant C = C(α 1 , α 2 , α 3 , γ 1 , γ 2 , x, y, z) we can alternatively express the series (7.59) in terms of f k+x,l+y,n+z = C g k+x,l+y,n+z , Hence, for any triplet (x, y, z) the series (7.59) furnishes a formal solution of the Yangian PDEs. We find 36 zeros of the fundamental solutions, i.e. combinations of x, y, z for which the series terminates. However, only for 2 of these 36 possibilities, u, v and w are the effective variables of the series: (x, y, z) = (0, 0, 0), (x, y, z) = (0, 0, 1 − γ 2 ). (7.65) We note the shift identity which alternatively can be expressed as Hence, we may relate the second series solution specified by (7.65) to a shifted version of the first: Making the ansatz we can fix the coefficients c 1 and c 2 by the two limits discussed in the following: or alternativelyc This result agrees with the expression given in [32].
Two-Point One-Mass Limit. Consider the coincindence limit and thus with lim w→0 w n = δ n0 yields We may thus conclude This can be compared with the two-point integral (7.10) for a 2 + a 3 → a 2 : which fixes Two-Point Zero-Mass Limit. Note that in the limit of a vanishing propagator power a 1 → 0, we have c 1 = 0 and the c 2 term yields the expected propagator type contribution with an overall factor of Gamma functions given in (7.5): This fixes the coefficient c 2 for a 1 = 0 to Note that using the recursions from acting with P D+1 as discussed in Section 3.4 fixes the a 1 dependence of the two coefficients to be , c 2 = f 2 (a 2 , a 3 ). (7.80) This shows that in fact even for a 1 = 0 we have

Dual Conformal Integrals
In this section we systematically apply the Yangian symmetry discussed above to constrain one-loop Feynman integrals with massless and massive propagators. In order to have the full Yangian symmetry, we consider the case of dual conformal integrals, i.e. the Yangian constraints on integrals for which the condition is satisfied by the propagator powers a j entering an n-point vertex in the (region momentum) Feynman graph. Again we start from the simplest examples and increase the complexity step by step. For the non-dual-conformal examples we considered in the previous section we could in principle simply take the dual conformal limit to obtain the solution. However, since the dual conformal integrals are invariant under the whole Yangian and depend on less variables, the resulting constraints allow us to bootstrap more examples than above. Moreover, in the dual conformal case it is natural to employ a different set of (constrained) variables.

2 Points, m 1 0: Rational
For a single massive propagator, the two-point integral has no independent variable. Conformal symmetry fixes it to take the form with an undetermined constant c 1 . This combination is also invariant under the level-one generators. To fix the normalization we can e.g. straightforwardly compute the integral in the incident point limit (or compare with the result given in [37]). From the conformal case we can then read off the coefficient in (8.2) to be We can understand this solution as a fusion rule for a massless and a massive propagator in the dual conformal case a 1 + a 2 = D: Here we have defined This rule is similar to the so-called group relation for two massless propagators in the nonconformal situation, cf. Section 7.2. These relations imply that chains of propagators connected by dual conformal vertices can be reduced according to For the two-point integral with m 1 = 0, m 2 = 0 there is a single conformal variable and one can consider different choices for this variable that lead to different types of well known differential equations. In this subsection we consider a choice that results in two-parameter Legendre functions. This choice is natural since the number of parameters of the class of functions matches the number of free propagator weights of the integral. We choose the conformal variable to be (cf. [3]) Direct Solution of PDEs. For compactness we set and make the conformal ansatz Note that while it may seem unatural at first sight, pulling out the factor (1 − v 2 ) β/2 leads to the canonical form of the below PDE. The systematics behind this prefactor is understood by writing for the matrix v with elements v jk = (x 2 jk + m 2 j + m 2 k )/2m j m k , see also Section 8.5. Acting on this function with the level-one momentum generator P extra leads to the following associated Legendre differential equation: with P β α and Q β α being the associated Legendre function of the first and second kind, respectively. The coefficients can be fixed using numerical data points to find Solution of Recurrence Relations. As an alternative to the above solution, we can apply the steps of the formal boostrap outlined in Section 4.2. We define a functionφ by (8.14) and we make the series ansatzφ such that the level-one momentum constraints turn into the recurrence relation This relation is straightforwardly solved and yields the two fundamental solutions where we abbreviateâ = a/2. Indeed, numerically we find the interesting identity From Legendre to Gauß. We note the relation for v > 1, which implies the alternative representation in terms the Gauß hypergeometric function 2 F 1 : This suggests to introduce the variable Numerically we find for v > 1 that which implies the representation In Section 9 we conjecture a generalization of this expression in u-type variables to higher point integrals. In the subsequent Section 8.3 we will bootstrap the same integral in the variable w = 1/u.

Unit Propagator Powers.
In order to compare with the discussion in [25] we can take the limit a j → 1 or α, β → −1/2 in D = 2 which implies Alternatively, we can solve the above PDE (8.11) directly for a j = 1 which yields Comparing to (8.24) for v > 1, the constants are fixed to As another alternative variable to the above we set and we define the conformal function φ via Here it is convenient to introduce three shorthands α = 1 2 (a 1 − a 2 + 1), β = a 1 , γ = 2α. (8.29) Then acting on the above function with the level-one momentum generator produces the Gauß hypergeometric differential equation 30) which is solved by The coefficients for this dual conformal integral can be fixed from a limit of the non-dualconformal two-point integral (see (7.52)) to find Unit Propagator Powers. It is interesting to evaluate the limit a 1 , a 2 → 1 for D = 2, which corresponds to α → 1/2 and β → 1 and yields One-mass Limit. In the limit where m 2 → 0 we have w → 0 such that 2 F 1 → 1. We are therefore left with where we assumed a 1 > a 2 . This matches the result from Section 8.1 for In the case of the three-point integral with two vanishing masses and evaluated at the dual conformal point D = a 1 + a 2 + a 3 , the integral depends on a single variable which we choose as .

(8.38)
Then level-one momentum invariance of the integral directly implies Gauß' hypergeometric differential equation 39) which is solved by The coefficients are fixed by the below limits to read Coefficients from Star-Triangle Relation. To determine the coefficients c 1 and c 2 , we take the limit m 1 → 0 which implies u → 0 and thus 2 F 1 → 1. This should be compared with the star-triangle relation for the massless three-point integral given in (6.2): We conclude that Two-point limit Another way to relate this integral to a simpler object is the two-point limit in which x 3 → x 2 . In this case, the integral should be given by the conformal one-mass two-point integral from Section 8.1, i.e. we expect lim On the other hand, performing this limit using the general solution (8.40), we find for γ < 1 that This allows us to immediately fix  In Section 4.1 this integral was discussed as an example for how to extract the explicit PDEs in the conformal variables. The P-equations split into two contributions coming from P y=0 and P extra that annihilate the integral separately. We will show that we can find the fundamental solution for the integral by working only with the constraints arising from P y=0 . Reading off the coefficients of x µ jk /m j m k with j, k = 1, 2, 3 yields the following three differential operators that annihilate φ(u, v, w): We make the series ansatz φ(u, v, w) = k,l,n g kln u k v l w n , (8.53) such that (8.50,8.51,8.52) translate into three recurrence equations for the coefficients g kln , e.g. from (8.50) we obtain 0 =2a 3 (k + 1)g k+1,l,n + 2(k + 1)lg k+1,ln + 2(k + 1)ng k+1,ln − (l + 1)(n + 1)g k,l+1,n+1 − (k + 1)(n + 1)g k+1,l,n+1 − (k + 1)(l + 1)g k+1,l+1,n . (8.54) These reccurence equations can be solved in Mathematica, which yields the following fundamental solution that is unique up to overall constants: Remember that here we have D = a 1 + a 2 + a 3 . For a fixed constant C and for integer k, l, n the fundamental solution g kln can be written as f k+x,l+y,n+z = C(a 1 , a 2 , a 3 , x, y, z) g k+x,l+y,n+z , The transpositions of the three external legs of the integral translate into (a 1 ↔ a 2 , l ↔ n), (a 2 ↔ a 3 , k ↔ l), and (a 1 ↔ a 3 , k ↔ n), which are manifest symmetries of the above fundamental solution. The fundamental solution f kln has 29 zeros which corresponds to 29 possible choices for the set (x, y, z) such that the following series terminates at an upper or lower bound Numerical analysis shows that only the choice (x, y, z) = (0, 0, 0) of the 29 possible zeros (x, y, z) of the fundamental solution leads to a series G xyz with effective variables u, v, w. This leads us to an ansatz φ = c 1 G 000 (u, v, w), (8.59) in terms of Srivastava's triple hypergeometric function H C , cf. e.g. [38]: This series is known to converge for Comparing the above ansatz to numerical data from the Feynman parametrization of the integral I m 1 m 2 m 3

3•
we can fix the overall coefficient to find (for D = a 1 + a 2 + a 3 ): ; u, v, w .
We refer to the kinematic region covered by a series representation in these variabels as region B. We note that for Euclidean x 2 jk and m j these variables are never small (as opposed to (8.48)), while at the same time we expect the corresponding series solution to converge only for small u, v, w.
On the support of these equations, the recurrences following from the invariance under P extra take the form 0 = (k + l + a 1 )(k + n + a 2 )g kln − (k + 1)(k + 2)g k+2,l,n , (8.69) 0 = (k + l + a 1 )(l + n + a 3 )g kln − (l + 1)(l + 2)g k,l+2,n , (8.70) 0 = (k + n + a 2 )(l + b + a 3 )g kln − (n + 1)(n + 2)g kl,n+2 . (8.71) While we have not determined the general solution to these equations, we note that they are solved by with the shorthandk = k/2. To motivate this expression we note that in (8.18) we have found that the two-point two-mass integral can be expressed as which shows that (8.72) is the natural three-point generalization of the two-point summand in (8.73).
Based on (8.72) we can now investigate the possible basis series for our ansatz. For a constant C = C(a 1 , a 2 , a 3 , x, y, z) we can write f k+x,l+y,n+z = C g k+x,l+y,n+z .
This function has 17 zeros (x, y, z) in k, l, n but only the triplet (x, y, z) = (0, 0, 0) leads to a series in the effective variables u, v and w. We thus conclude that for small u, v, w the correct Yangian invariant is proportional to this series: The coefficient can be fixed using numerical data points such that we find In the following section we will conjecture an n-point generalization of this expression.

All-Mass Yangian Invariant n-Gon Integrals
Based on the evidence from the previous Section 8, we propose the below conjectural series representations for the dual conformal, i.e. Yangian invariant n-point one-loop integrals with all propagators massive: Note that these variables can become small for real Euclidean x ij and m j , e.g. for large masses. This is relevant since we believe that the below series merely converges for small u ij . We refer to this series representation as the series in region A as opposed to the B-series presented in the following subsection.
When expressed in the above variables we conjecture the dual conformal n-point all-mass integral to be given by the expression where B n = {12, 13, 23, ..., (n − 1, n)} is the set of all ordered pairs of distinct numbers between 1 and n, whereas B n|j is the subset of B n which is comprised of pairs containing j. We note that the Feynman parametrization of the corresponding dual conformal all-mass n-point integrals in terms of the variables (9.2) is given by Let us present our evidence for the correctness of the above series representation at lower points.
Evidence for the Conjecture n = 1. For n = 1, the entire sum collapses, leaving only in agreement with (7.4). n = 2. In the two-point case, the expression reads (a 1 ) k 12 +k 13 (a 2 ) k 12 +k 23 (a 3 ) k 13 +k 23 ( D+1 2 ) k 12 +k 13 +k 23 u k 12 12 u k 13 13 u k 23 9.2 n Points, m 1 . . . m n : Conjecture in Region B The above results suggest a second series representation, which is closely related to a representation given in [24]. This series does not converge in the Euclidean region of the Feynman integrals that we have been focussing on so far, but its analytical continuation is conjectured to agree with the integral. The B-series is expressed in terms of the variables For Euclidean choices of x ij and m j , these variables do not become small, which would correspond to the kinematic region where we expect the below series to converge. The conjectured series reprensentation reads Here again B n = {12, 13, 23, ..., (n − 1, n)} represents the set of all ordered pairs of distinct numbers between 1 and n and B n|j is the subset of B n which contains all pairs containing j. Again we note that the Feynman parameter representation of the conformal massive n-gon integrals in terms of the variables (9.8) takes the form 9 I m 1 ...mn For the series (9.10) this Feynman representation is particularly useful to numerically verify examples of the conjectured representation, since for real values of the m j and x jk the variables v ij in the Euclidean signature do not become small and the series does not converge.
Evidence for the Conjecture n = 1. As in the case of the series representation in region A, for n = 1 the sum collapses and agrees with (7.4). n = 2. The result agrees with(8.18) that we found from Yangian symmetry and with the numerical Feynman parameter integration. n = 3. The result agrees with (8.76) that we found from Yangian symmetry and with the numerical Feynman parameter integration. n = 4. The conjectured series agrees with the numerical evaluation of the above Feynman parametrization. Moreover, it provides a solution to the P Yangian PDEs, cf. Appendix B. n = 5. The conjectured series agrees with the numerical evaluation of the above Feynman parametrization.

Unit Propagator Powers for 2,3 and 4 Points
In this section we compare the above A-and B-series to some lower point expressions for the Yangian invariant all-mass integrals with unit propagator powers.
A-Series for 2 Points. The unit propagator power limit of the result in Section 8.2 in On the other hand, we can evaluate the A-series for a 1 = a 2 = 1 with D = 2 and using v 12 = −2u 12 + 1. This yields the relation where 10 Here the matrix G is defined to have matrix elements G jk = v jk . We can compare this result to the above 3-point result in u-type variables (9.7). Here we note that v jk = −2u jk + 1. and for the prefactor of the series expression in D = 3 have Γ 3/2 /Γ 3 = √ π/4. We thus conclude that in the region of convergence (8.61) for H C we have Expanding the left hand side in Mathematica assuming 0 < u jk < 1 indeed shows that at least up to and including order 8 in the variables u 12 , u 13 , u 23 , the expansions of both sides of this equation coincide.

B-Series for 3 Points: Comparison with Nickel.
Similarly, we can compare the above result (9.15) by Nickel with the B-series (9.10), which is actually formulated in the same vtype variables. Also here we find agreement at leading orders when expanding the two expressions, i.e.
B-Series at 4 Points: Comparison with Murakami-Yano. In [25] the Murakami-Yano formula [41], which gives a compact expression for the volume of a hyperbolic/spherical tetrahedron, has been leveraged to give a concise dilogarithmic expression for the all-mass box integral in four dimensions with unit propagator powers. This result provides a valuable cross check of our series representation (9.9), which we deem worth detailing. The volume of a spherical tetrahedron is most elegantly phrased in terms of its dihedral angles. We therefore begin by making explicit the relation between these angles and our variables v ij . Let G be the matrix whose elements are given by the variables v ij , i.e.
In terms of the function L(z + ), the volume of the tetrahedron described by the matrix G is given by Utilizing this expression, the all-mass box integral can be expressed as where G 0 ij = G ij m i m j . We find this expression to be in perfect numerical agreement with the series representation (9.9) in the Euclidean domain.

Outlook
The results presented in this paper suggest plenty of different directions for further investigation. Let us detail a few.
In the case of one-loop Feynman integrals we have demonstrated that Yangian symmetry is in fact highly constraining. For the simple examples studied in Section 7 and Section 8, this results in a small basis of hypergeometric series whose linear combination yields the integral under study. Here the number of basis elements depends on the chosen variables as can for instance be seen in Section 7.4 and Section 7.5, where the same integral is studied for two different choices of variables. In particular, the solution basis is generically expected to grow with the number of variables, as becomes apparent in the case of the 9-variable massless double box and hexagon integrals considered in [14]. Here the close connection to the Mellin-Barnes approach deserves further study. In particular, it would be desirable to find a symmetry principle that selects the specific subsets of formal Yangian invariants that span the solution, see [15]. Eventually it seems plausible that to fix these linear combinations using only conformal Yangian symmetry, kinematical configurations with on-shell external particles have to be taken into account. These may result in deformations of the symmetry generators similar to the case of scattering amplitudes in N = 4 SYM and ABJM theory and thus in additional relations, cf. [42].
In certain cases and for particular choices of the conformal variables, the Yangian bootstrap selects a single series solution to the symmetry constraints. In these cases merely an overall coefficient remains to be fixed in order to determine a representation of the integral. In particular, this is the case for the all-mass n-gon integrals subject of Section 9 and allowed us to conjecture two different single series representations (9.3) and (9.9) for generic oneloop integrals in D spacetime dimensions. While the outstanding properties of these integrals have been studied in the past for unit propagator powers (cf. e.g. [25]), the novel Yangian symmetry sheds new light on their distinguished role. It would be very interesting to further explore the space of Yangian invariant integrals in order to identify families with similarly beautiful properties. Here the next step is to proceed to two loops. While the analysis of the massless double box in [14] shows the increase of complexity for a higher number of external points, the present paper suggests that it may be beneficial to first consider situations with massive propagators. With regard to the one-loop integrals, it would be interesting to better understand the mathematical properties of the two series representations (9.3) and (9.9). While we have not found a name for the B-series that closely resembles that given in [24], at least for n = 2 and n = 3 points the A-series coincides with Gauß' hypergeometric function 2 F 1 and Srivastava's triple hypergeometric function H C , respectively, and thus can be assumed to represent a useful generalization.
When proceeding to more complicated examples, we note that at loop orders beyond two, the statements about level-one Yangian symmetry are still conjectural, cf. Table 2 and [3]. It would be important to make progress on understanding these cases in detail. Ideally one could find an analytic proof similar to the one in the massless case using the Lax operator formalism [2]. Alternatively, it would be interesting to systematically map out the space of higher loop integrals using advanced numerical integration techniques. Into this direction it would also be interesting to study massive Feynman integrals with particles different from scalars which appear in non-scalar fishnet theories [43]. In fact, in the massless case certain brick wall Feynman graphs including fermionic lines were found to be Yangian invariant [7] and thus represent a natural starting point.
Physical application of our results asks for an extension into two further directions. Firstly, while Feynman integrals with generic propagator powers are clearly of interest, it would be desirable to better understand the considered Yangian approach for integer (in particular for unit) propagator powers, which at one loop order results in polylogarithmic expressions. Here the situation of the massless box integral was explicitly discussed in [14], which underlines the importance of the choice of kinematic variables. Secondly, while we have focussed on the case of Euclidean spacetime signature, there is an obvious demand to explicitly discuss the results in Minkowski signature. Though this step should not modify the Yangian constraints (and thus the solution basis), identifying the correct linear combination needs more care. For the case of the massless box integral, the results of [44] indeed show that also in Minkowski spacetime the integral is spanned by the Yangian invariant building blocks given in [14].
The four-point Basso-Dixon diagrams of [45] represent another example of Feynman integrals with an intriguing connection to integrability, see also [46]. Using sophisticated integrability techniques from AdS/CFT and the Steinmann relations, a conjecture for the polylogarithmic result for these integrals was given at generic loop order. While a Yangian symmetry has not been formulated in this case, these integrals can be understood as coincident point limits of the more generic Yangian invariant fishnet integrals discussed in [2,7]. If one assumes a connection between this Yangian symmetry and the simplicity of the expressions given by Basso and Dixon, one may wonder whether massive propagators can be introduced into their formula.
As mentioned in the introduction, the conformal Yangian and its massive generalization studied here can be considered as the closure of two distinct conformal algebras. As such, it would be very interesting to identify its place within the large landscape of results on the conformal and momentum space conformal bootstrap. Clearly, the Yangian constraints can be studied independently of Feynman integrals and it would be interesting to search for further applications. These might correspond to lifting a conformal setup to an integrable one. Even if one does not consider the full Yangian, it seems natural to generalize the various applications of the momentum space conformal bootstrap (see e.g. [22] and references therein) to the massive extension of the algebra given in (5.10).

A Yangian Level-One Generators
We provide the expressions for all Yangian level-one generators over the (dual) conformal algebra:

B All-Mass Yangian PDEs for 4 Points
Here we give some details on the all-mass Yangian constraints for four points. with u jk as defined in (9.2). From the P invariance equation we can read off the coefficients of the vectors x µ jk /m j m k to find the annihilators of the function φ, e.g.