Near mass-shell double boxes

Two-loop multi-leg form factors in off-shell kinematics require knowledge of planar and nonplanar double box Feynman diagrams with massless internal propagators. These are complicated functions of Mandelstam variables and external particle virtualities. The latter serve as regulators of infrared divergences, thus making these observables finite in four space-time dimensions. In this paper, we use the method of canonical differential equations for calculation of (non)planar double box integrals in the near mass-shell kinematical regime, i.e., where virtualities of external particles are much smaller than the Mandelstam variables involved. We deduce a basis of master integrals with uniform transcendental weight based on the analysis of leading singularities by means of the Baikov representation as well as an array of complementary techniques. We dub the former asymptotically canonical since it is valid in the near mass-shell limit of interest. We iteratively solve resulting differential equations up to weight four in terms of multiple polylogarithms.


Introduction
Infrared structure of off-shell observables in massless gauge theories attracted attention in the past couple of years.Within the context of the maximally supersymmetric Yang-Mills theory (aka N = 4 sYM) this kinematical regime can be addressed in a fully gauge-invariant fashion by studying the theory on its Coulomb branch [1].With a proper choice of vacuum expectation values for the model's scalar fields, one can mimic the off-shellness of the unbroken gauge symmetry phase with nonvanishing masses for external particles only, while keeping all states propagating in internal lines of Feynman graphs massless.This regime is of phenomenological interest in physical theories like QCD.As opposed to the fully massless case where infrared singularities arise as poles in the dimensional regularization parameter ε = (4 − D)/2, in the nearly mass-shell regime of virtual amplitudes, they are replaced by the logarithms of external states' virtualities.An orthodox universality of infrared physics would suggest that critical exponents in both cases will be given by the very same function of the Yang-Mills coupling constant.However, recently this was demonstrated to be far from the truth.
Four- [2] and five-leg [3] scattering amplitudes as well as two-particle form factors [4,5] were the first few examples to exhibit a novel feature of the near mass-shell kinematics as opposed to the fully massless regime.While the infrared physics in the latter was known, since the inception of QCD [6], to be governed by the cusp anomalous dimension [7,8] the former involved a different function of the coupling, the so-called octagon anomalous dimension [9][10][11].To further elucidate its role, one needs to address more complicated observables containing more scales.They are of interest for several reasons.First, it is desirable to test < l a t e x i t s h a 1 _ b a s e 6 4 = " p 3 V z 8 r Q 1 g e K B 8 n 2 Q 0 V 6 W e 0 C J 2 g o = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l K / T g W v H i s a G u h D W W z 3 b R L N 5 u w O x F K 6 E / w 4 k E R r / 4 i b / 4 b t 2 0 O 2 v p g 4 P H e D D P z g k Q K g 6 7 7 7 R T W 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R 2 8 S p Z r z F Y h n r T k A N l 0 L x F g q U v J N o T q N A 8 s d g f D P z H 5 + 4 N i J W D z h J u B / R o R K h Y B S t d J / 0 a / 1 y x a 2 6 c 5 B V 4 u W k A j m a / f J X b x C z N O I K m a T G d D 0 3 Q T + j G g W T f F r q p Y Y n l I 3 p k H c t V T T i x s / m p 0 7 J m V U G J I y 1 L Y V k r v 6 e y G h k z C Q K b G d E c W S W v Z n 4 n 9 d N M b z 2 M 6 G S F L l i i 0 V h K g n G Z P Y 3 G Q j N G c q J J Z R p Y W 8 l b E Q 1 Z W j T K d k Q v O W X V 0 m 7 V v U u q x d 3 9 U q j n s d R h B M 4 h X P w 4 A o a c A t N a A G D I T z D K 7 w 5 0 n l x 3 p 2 P R W v B y W e O 4 Q + c z x 8 A D o 2 V < / l a t e x i t > p2 < l a t e x i t s h a 1 _ b a s e 6 4 = " F 9 e K K h R Z g X 8 D 6 E T m H L q 4 R / q y W S E = " > x 6 c F + f d + V i 0 5 p x s 5 h j + w P n 8 A d q n j P E = < / l a t e x i t > q < l a t e x i t s h a 1 _ b a s e 6 4 = " a a W Z j  < l a t e x i t s h a 1 _ b a s e 6 4 = " a a W Z j 1 K l z W N u V K z Y r s 2 W T h z + a L 4 = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e x q f B w D X j x G N A 9 I l j A 7 6 U 2 G z M 4 u M 7 N C W P I J X j w o 4 t U v 8 u b f O E n 2 o N G C h q K q m + 6 u I B F c G 9 f 9 c g o r q 2 v r G 8 X N 0 t b 2 z u 5 e e f + g p e N U M W y y W M S q E 1 C N g k t s G m 4 E d h K F N A o E t o P x z c x v P 6 L S P J Y P Z p K g H 9 G h 5 C F n 1 F j p P u m f 9 8 s V t + r O Q f 4 S L y c V y N H o l z 9 7 g 5 i l E U r D B N W 6 6 7 m J 8 T O q D G c C p 6 V e q j G h b E y H 2 L V U 0 g i 1 n 8 1 P n Z I T q w x I G C t b 0 p C 5 + n M i o 5 H W k y i w n R E 1 I 7 3 s z c T / v G 5 q w m s / 4 z J J D U q 2 W B S m g p i Y z P 4 m A 6 6 Q G T G x h D L F 7 a 2 E j a i i z N h 0 S j Y E b / n l v 6 R 1 V v U u q x d 3 t U q 9 l s d R h C M 4 h l P w 4 A r q c A s x 6 c F + f d + V i 0 5 p x s 5 h j + w P n 8 A d q n j P E = < / l a t e x i t > q < l a t e x i t s h a 1 _ b a s e 6 4 = " r a 3 8 U a F q g c g 0 the infrared factorization and off-shell universality in circumstances which involve multiple scales at the same time.Second, for near mass-shell scattering amplitudes with more than five legs and form factors with more than two, there is a residual finite contribution free from infrared logarithms but depending on Mandelstam-like variables.These are known as remainder functions.A natural question arises whether they are the same both in on-and off-shell regimes.Given that the critical exponents are different, one would expect them to differ as well.But one needs explicit verifications.
Remainder functions in the massless case of planar N = 4 sYM are endowed with a stringy description [12,13] in terms of an effective two-dimensional world-sheet [14,15].The string in question is the so-called GKP string [16] with its energy density determined by the cusp anomalous dimension.The string supports a set of elementary excitations of a T-dual theory, with their dispersion relations and scattering matrices known exactly in 't Hooft coupling constant [17].One may wonder then, given the vacuum of off-shell observables is determined by the octagon anomalous dimension, whether the spectrum of excitations that live on it and their interactions remain the same as on the GKP background.This is a long-term goal.On the way toward it, one needs 'experimental' data from explicit field theoretical calculations to confirm or deny this expectation.A first step will be undertaken in this paper.
In this work, we calculate Feynman integrals of double box families, see Fig. 1, relevant for the problem of three-leg form factors ⟨p 1 , p 2 , p 3 |O(q)|0⟩ at two-loop perturbative order.As we advocated above, we are particularly interested in the kinematical situation where the off-shellnesses of external particles with momenta p ℓ (ℓ = 1, 2, 3) are equal and small compared to the virtuality q 2 associated with the operator O, -the lowest super-component of the stress-tensor multiplet.To this end, we will rely on the method of differential equations [18,19] in its modern day reincarnation that employs canonical bases [20].To date, this is the most powerful and efficient technique to tackle multi-scale Feynman integrals.Recently, it was successfully applied to the planar double box integrals (left panel in Fig. 1) with three [21] and four [22] external squared momenta being off the light cone and all internal lines being massless.A basis of uniformly transcendental (UT) integrals was established and its symbol alphabet was analyzed.The latter was shown to be populated by letters expressed in terms of multiple square roots but no solution to the differential equations was offered.We will fill in this gap below for our kinematical situation and supplement it with an analysis of nonplanar graphs as well (right panel in Fig. 1) which are far more complex.We will demonstrate that the near mass-shell limit provides sufficient simplification of the structure of canonical differential equations to offer a solution in terms of multiple polylogarithms [23] with all square roots gone from their arguments.
Our subsequent presentation is organized as follows.In the next section, we set up the kinematics and classes of Feynman integrals to be studied both for planar and nonplanar families.Next, we move on to the construction of canonical bases.We will discuss the two cases in parallel providing necessary technical details as needed so that a curious reader could reproduce all results if desired.We start in Section 3 by building an initial basis of master integrals (MIs) making use of the integration-by-parts (IBP) to deduce a primary basis.Then in Section 4, we use leading singularities and a variety of other available techniques to cast them in the canonical form for the planar case.The nonplanar family if far more complex and while we manage to build canonical basis for lower sectors, we encounter elliptic cases and thus turn to the asymptotic limit in question where these degenerate into poles.The asymptotically canonical basis for non-planar graphs is presented in Section 5. We provide solutions to the resulting differential equations and determine corresponding integration constants in Section 6 making use of a variety of criteria which bypass the necessity for explicit evaluation of parametric integrals.In both situations, explicit results are given up to weight-four in terms of multiple polylogarithms and they are further recast in terms of classical Euler polylogarithms and an additional two-argument polylogarithm Li 2,2 at weight-four.Mathematica notebooks and ancillary files attached with this submission spill sufficient 'secrets of the trade' for a newbie to familiarize oneself with the subject and can be regarded as 'blueprints' for construction and solution of canonical differential equations for any (non-elliptic) family of Feynman graphs.

Setting up conventions
To begin with, let us establish our notations.The kinematics of interest corresponds to the momentum flow from the operator source O(q) to off-shell external particles, obeying the conservation condition q = p 1 + p 2 + p 3 .We introduce three Mandelstam variables according to which are related by the equation Nevertheless, throughout our subsequent analysis, we will treat u, v and w as independent since this provides stringent checks on the correctness of our derivations.Above, we introduced equal Euclidean virtualities for all particle legs p 2 ℓ ≡ −µ.Since the overall scale of a Feynman integral can be always unambiguously restored on dimensional grounds, we will set q 2 = −1 in what follows.
The families of graphs which take centre stage in this paper are shown in Fig. 1.In spite of the fact that Feynman integrals contributing to physical observables are finite for nonvanishing off-shellness, nevertheless we will work with a dimensionally regularized theory in order to be able to apply IBP reduction which requires a D-dimensional setup to render bases of sought after MIs complete.The two-loop non-and planar integrals are determined by a set of massless propagators 1/D j (for j = 2, . . ., 7) and two irreducible scalar products D 8,9 defined according to the momentum flow exhibited in Fig. 1 as with the remaining denominator corresponding to the planar and non-planar cases, respectively.All indices a i are integers with a 8,9 ≤ 0. Let us turn to these two families one by one.

Primary basis of MIs and differential equations
Let us begin with the planar graph as a case study.It was previously addressed in Refs.[21,22], however, we will use this more familiar family to set up our formalism so that we could be more concise in our following presentation of the non-planar graph, which is computationally more demanding but does not bring anything new to the table to a certain degree.
Preliminary counting of MIs can be done with a variety of available tools, say with Mint [24] or the modular component of FIRE [25], which was the go-to tool in our analysis.Constructing a list of integrals in the top, i.e., level-seven sector, obtained by inequivalent permutations of two indices equal to two G 221111100 , we prepare start files and generate symmetry relations with LiteRed [26,27].A modular arithmetic IBP then yields an initial set of 74 MIs.We give preference to Laporta-reducible values of indices being equal to 2 since experience with canonical bases have taught us that these more likely than not be endowed with single leading singularities and thus serve as UT candidates.Next, we use FindRules command of FIRE to deduce 10 symmetry equations between MIs in our preferential basis, thus reducing their number to 64.However, this is not the end of the story.We can further determine 'hidden' relations as well.To accomplish this, we create lists of integrals sufficiently close to the preliminary set and containing these as a subset: it includes integrals with none, one and two indices set to 2. Then an IBP reduction reveals additional two relations among them reducing their total number to 62.At this step, it is always advisable to verify that thus obtained basis does not yield the so-called 'bad' denominators according to the nomenclature of Ref. [28].In fact, we find none.But at level five, IBP yields quite lengthy denominator polynomials in the Mandelstam variables and virtuality and these can be traded however for a significantly more compact ones.This provides us with a solid starting set of 62 preliminary MIs I for our subsequent analytical analysis which we will use from now on as the option #masters for IBP reduction with FIRE.All of the above steps are presented in Sections 1 through 10 of the attached Mathematica notebook A2Zdbox.nb and output is saved in the subdirectory dbox.An analysis identical to the one just discussed is performed then for the non-planar family to give us a set of 97 MIs I and stored in the subdirectory nbox.
After these preparatory studies, we move on to the construction of the derivatives in the Mandelstam variables and virtuality by performing differentiations with LiteRed.The differential dI = du ∂ u I + dv ∂ v I + dw ∂ w I + dµ ∂ µ I then needs to be IBP-reduced back to the MIs I thus generating the sought-after differential equations with i = u, v, w, µ.While the analytic IBP reduction for the planar family takes a matter of hours on a typical machine, the non-planar case is far more computationally demanding.For instance, the reduction of level-seven integrals from the left-hand side of the differential equations (3.1) given in the accompanying file intsde-nbox2.mdown to MIs #3 G 010110000 and #4 G 101001000 from pr-nbox2.mtakes a staggering 7 and 10 days, respectively.To crosscheck that the resulting tables are indeed correct, we relied on modular arithmetic runs with FIRE with subsequent balanced rational reconstruction developed in Ref. [29].With an MPI parallelization of 1024 cores of ASU's Sol supercluster [30], the 7 → #3 IBP took 58 hours with Flint [31] and 46 hours with Symbolica [32] but indeed confirmed our earlier analytical findings.We provide detailed account of the derivation in Section 11 of the accompanying notebook A2Zdbox.nb for the planar graph.

Canonical basis
Now the main task at hand is to transform the basis of MIs I = T •J such that the differential equations (3.1) admit their canonical form with each element of the A-matrices being Fuchsian, i.e., possessing simple poles only, and εindependent [20].To this end, we need to determine viable UT candidates from our primary list of MIs.Provided this procedure is successful, a differentiation of these pure UT integrals will reduce their transcendental weight by one and, thus, the right-hand side of (3.1) will have to be proportional to ε, which carries the transcendentality weight −1.To practically implement this strategy, we rely on the well-known conjecture that connects uniform weight integrals with the properties of their integrands, namely, that the singularities of an integrand are locally of logarithmic type [33,34], see, e.g., [35] for a comprehensive review.
As the calculation of unitarity cuts is in general downright easier than solution of integrals per se, the idea is to use the former for the identification of Feynman integrals that correspond to pure functions.To perform multidimensional unitarity cuts efficiently, one has to rely on an appropriate parametrization.Since Feynman integrals possess integrands which are rational functions of propagators and ISPs, it is only natural to choose these as integration parameters z i ≡ D i .To date this is the most concise framework which is known as the Baikov representation [36,37], see Refs.[38,39] for comprehensive reviews.This form of integrals trivializes computation of unitarity cuts.The so-called leading singularities correspond to taking the maximal cut, i.e., successive residues in all z i = 0, followed by residues in composite singularities emerging along the way from any Jacobian factors [40].This completely localizes all integrations and provides a function of external kinematical variables, which once being divided out from the Feynman integral in question yields a pure UT candidate with constant leading singularity.Of course, for a given integral, there could be multiple ways to localize all integrations depending on the order of taking residues and this can yield different leading singularities.Only integrals with a single leading singularity can be autonomously recast as UT, while in cases where there are more that one, linear combinations of these have to be studied as well.It is important to realize that UT candidates found this way may not correspond to MIs of the traditional Laporta algorithm.This is the reason why we chose from the very beginning to favor MIs having indices equal to 2 in our IBP studies.Leading singularities analysis completely fixes the diagonal blocks of the A-matrices, which do not mix MIs at different levels.Then we move on to study sub-maximal cuts to find corrections from lower-level subsectors.
In our analysis, we relied on the Mathematica implementation of the Baikov parametrization via Baikov.mpackage developed in Ref. [41].Starting from the lowest sector, we iteratively constructed the canonical basis of MIs for the planar graph.Since a more general case was already studied in the literature [22], we provide only sporadic details in Sections 12-14 of the ancillary file A2Zdbox.nb with final results for the canonical basis and all A-matrices given in dboxCan62.mand AuPC.m,AvPC.m,AwPC.m,AmPC.m, respectively.

Asymptotically canonical basis
The analysis alluded to above immediately convinces us that the non-planar graph (right panel in Fig. 1) possesses elliptic sectors [42], see Ref. [43] for a thorough review, implying that some leading singularities reside on elliptic curves [44,45] rather being merely algebraic.However, they smoothly degenerate into the latter as we send the virtuality µ down to zero.Since at the end of the day, all we are after is the asymptotic behavior of our MIs as µ → 0 up to terms which vanish in µ, we can implement this limit on differential equations for the primary set of 97 MIs and then construct what we call as the asymptotically canonical basis.
Thus, we change the strategy for basis construction in the non-planar case.Namely, we tend to assume generic values for all variables (u, v, w as well as µ) as long as we encounter only algebraic leading singularities and swiftly pass to the asymptotic consideration when it is no longer the case.In particular, this occurs in the two level-six sectors G 111101100 and G 111111000 , and the top level-seven sector G 111111100 .For these, we require the following properties to be fulfilled by the differential equations: (i) the 'virtuality' matrix M µ can be cast in the form with A 0 µ being a matrix of rational numbers, i.e., its elements are strictly independent of u, v and w; (ii) M i matrices for the Mandelstam variables i = (u, v, w) have well-defined finite limit as µ → 0 and do not possess elements with square roots; (iii) last but not least, resulting differential equations in u, v and w are canonical or can be made canonical with an appropriate similarity transformation.
As usual, we focus on diagonal blocks first but use now the leading order form of the elements of the matrices M i (with i = u, v, w, µ) as µ goes to zero for a proper choice of asymptotically canonical elements.It is easier to demonstrate it with a example of the G 111101100 sector.The primary set of MIs defining this sector is G 112101100 , G 121101100 } .
To start with we utilize the form of M µ = M 0 µ /µ + O(µ), with M 0 µ being a function of u, v, w and ε, to conclude that if we are to multiply the elements 2,3,5,6 with µ, after a similarity transformation the virtuality matrix will take the required form (5.1).We know a priori, however, that the ε-dependence of this seed basis will have to be adjusted later since, as rule of thumb, one associates ε 4 with MIs without any twos and ε 2−n for MIs with n twos in the first 7 positions.For now, it will do the job however.Next, we change the basis by multiplying each element in it with an unknown function of the Mandelstam variables Enforcing the ε-form on the resulting differential equations for this new basis, we solve the arising differential equations on the functions f (u, v, w) and get, after a gentle mixture of elements with each other and re-arrangement, The last two elements of this naive basis contain square roots of Mandelstam variables in functions f 5,6 .This is hardly a surprise since we completely ignored up to now the correct ε-dependence of the basis elements which resulted in erroneous differential equations for f 5,6 .If we do it in a proper manner, we observe a violation of the canonical form of the differential equations in u, v and w.Details on this calculation can be found in Section 1 of the accompanying Mathematica notebook AsyClevel6.nb.
The above finding instructs us to look further for a better choice of the elements 5 and 6 of this sector.The method of trial and error is quite tedious and exhausting, so we turn to the massless non-planar double box for inspiration.In the strict limit µ → 0, the first four elements if (5.4) vanish and one is left with just two elements.The massless non-planar box analysis demonstrates that indeed it possesses two level-six sectors with one of them containing two primary elements G 111101100 and G 111101200 .Construction of the canonical basis in this sector is performed in Section 2 of the notebook AsyClevel6.nb and offers two options for UT elements, namely, or We then build upon (5.5) to lift the analysis to the off-shell case in the asymptotic limit and construct the final form of the diagonal block of this sector in Section 3 of AsyClevel6.nb such that the last two elements in Eq. (5.4) are replaced with + εµv(u and respectively, and we simultaneously restored a proper relative ε-normalization of these elements.All other sectors can be analyzed in the same fashion.Though all diagonal sectors were brought by us to the ε-form, not all elements of the Amatrices are Fuchsian.Moreover the off-diagonal blocks are not even close to the required ε-form.However, their transformation to the canonical form is now purely algorithmic.The Fuchsian form is easily obtained making use of the code Canonica.m[46][47][48].The latter cannot handle, however, transformation of all off-diagonal elements in the differential equations to the ε-form using the Lee's trick [49].The latter is implemented in a powerful package Libra.m [50], which calls for an external Fermat [51] computer algebra system, though it works just as fine1 even with built-in Mathematica commands.This systematic procedure is demonstrated step-by-step in the ancillary notebook AsyCnbox.nbproving the output in the file nboxAsyCan97.m in the Canonica format.This culminates our quest to reach the asymptotically canonical basis in the non-planar case.

Integration and determination of integration constants
To summarize, in the previous section we determined the asymptotically canonical form of the differential equation with the one in the off-shellness µ being up to terms vanishing as µ → 0, with A 0 µ being a purely numerical matrix.Solving this leading order equation is simple and it provides a transformation to the µ-independent basis J 0 via a matrix exponent whose entries are expressed as linear combinations of µ −nε -terms accompanied by rational number coefficients.The basis J 0 solves in turn the asymptotically canonical equations where the elements of A 0 i = A i | µ=0 are expressed in terms of rational functions of the Mandelstam invariants u, v and w only.
The main advantage of the canonical form of differential equations (6.3) is that their solution can be written in terms of a path-ordered exponential and explicitly calculated via Chen's iterative integrals theory [52], where A 0 = duA 0 u + dvA 0 v + dwA 0 w and J 00 is a vector of integration constants.At each order in the ε-expansion it receives an independent set of unknowns (6.5) The solution (6.4) is independent of the choice of the path γ since the integrability of the differential equations is a zero-curvature condition, dA 0 − εA 0 ∧ A 0 = 0.In our analysis, we chose a piece-wise path In more practical terms, we take the first equation with i = u in (6.3) and solve it with respect to u.Then we turn to the v-variable.To eliminate the u-dependence from the differential equation, we form a difference between the right-hand side of (6.3) for i = v and the derivative in v of the solution found at the previous step.We then find its primitive in v. Finally, we repeat the procedure for w.This procedure is cast in a Mathematica module Integrator in the attached notebook AsySolCnbox.nb.The result is then given in terms of multiple polylogarithms (MPLs) [23].The latter are defined recursively via an integral iteration, e.g., However, in order to have a better handle on the analytical structure of our results, we recast them in terms of classical Euler polylogarithms2 whenever possible.It is well known that for the weight up to (and including) three, all MPLs can be traded to Li n 's, see, e.g., [53][54][55] (Chapter 2 on MPLs of the latter reference is available at [56]).At weight four, one needs to include an additional two variable MPL Li 2,2 to the minimal basis of classical polylogarithms as was observed in Ref. [57].This basis transformation was worked out and is conveniently implemented into the routine gtolrules.mdevised in Ref. [58].We just need to make sure that a proper 'branch' of MPLs is taken into account at a given point in the (u, v, w)-space since a single expression in terms of Li n 's and Li 2,2 is not sufficient to cover the entire space of (complexified) Mandelstam variables.At each p-order of the ε-expansion, arrays of integration constants c (p) have to be determined from a set of boundary conditions.However, we would like to avoid an explicit calculation of any Feynman integrals since, even in some corners of the phase space, they are very complex and, which is worse, quite numerous.Instead, we relied on several criteria to fix them such as (i) numerology, (ii) cancellation of unphysical poles, (iii) absence of imaginary parts and (iv) finite integrals.
(i) The first condition is self-explanatory.Using the fact that our asymptotically canonical MIs obey the property of being UT, we cast the integration constants into a product of rational numbers times values of the Riemann zeta function ζ p = ζ(p) of the same transcendental weight, with the employed convention ζ 0 = 1 and ζ 1 = 0 for the first two values of p.Then, computing the MIs at a random point for the Mandelstam variables with FIESTA [59], we confronted its results against numerical evaluation of our solutions.In this manner, we managed to confidently determine 84 r (p) j for p = 0, 2, 81 r j 's and 79 r j 's.The monotonically decreasing number of rationally reconstructed constants with increasing p is related to the loss of FIESTA's numerical precision and emergence of large rationals at higher orders in ε.
(ii) To alleviate the aforementioned problem and cross check correctness of previous numerological findings, we employed conditions for unphysical pole cancellation in the righthand side of the canonical differential equations, namely, at u + v = 0, v + w = 0 and w + u = 0.Then, decomposing the A-matrices in explicitly Fuchsian form we imposed the following equations on our basis a i,α J 0 | α=0 = 0 .(6.9) These provided a further set of 10,7,9 and 10 identifications/relations between integration constants at levels p = 0, 2, 3 and 4, respectively.
(iii) To further constrain the integration constants at order p, we used solutions at order p + 1 and required vanishing of imaginary parts as one approaches unphysical poles in Eqs.(6.9).This provided the value on the last r (3) 90 element from the solution at order ε 4 .The solution at fifth order in ε was used in conjunction with high-precision numerical computations of MPLs with the C++ package GiNaC [60] making use of a Mathematica interface from Ref. [61] and subsequent reconstruction of analytical expressions with the PSLQ algorithm [62].This allowed us to deduce 3 equations for r (4) j (j = 86, 92, 97).(iv) The implementation of the above three conditions fixed all but 3,6,7 and 8 integration constants for p = 0, 2, 3 and 4, correspondingly.Then, by a judicious choice, we found a set of 26 which were reduced with IBP identities to our set of 97 MIs.The resulting relations are divergent and pole cancellation in the Laurent ε-expansion provided an ultimate set of equations to completely fix the solutions at orders one through three.At fourth order, we obtained the last five integration constants whose numerical value to O(10 −3 ) were determined to be r To rationalize these, one has to either perform an analytic calculation of a very large set of Feynman integrals or increase the accuracy of their numerical evaluation to twelve decimal places with FIESTA or any other program.Currently, alas, this is beyond our reach.However, a particular combination of these constants shows up in the three-leg form factor [63], which allows us to eliminate one of them from the above list.It reads All steps of the above analysis for the non-planar family are thoroughly presented in the accompanying Mathematica notebook AsySolCnbox.nb.For the planar graph, it suffices to use just the first two conditions (i) and (ii).All solutions up to the same order in ε are quoted in the ancillary file AsySolCdbox.nb.

Conclusions
With this paper, we initiate a series of studies of multi-scale two-loop observables in N = 4 sYM.Currently, we constructed bases of UT MIs for double box planar and non-planar graphs in the kinematical limit of small virtualities of three external particles and an arbitrary invariant mass for the last leg.This is a preparatory study for a full-fledged analysis of the three-leg form factor of the stress tensor multiplet to be published separately [63].The bases in question were used to determine canonical form of differential equations in the Mandelstam variables u, v, w as well as the off-shellness µ.Solutions to these equations were constructed up to terms vanishing as µ → 0. The results for all master integrals were obtained as a double Laurent/Taylor expansion in ε/log µ up to (and including) weight four contributions.All integration constants were successfully fixed analytically for the exception of 5 coefficients at level four where they were found numerically with the accuracy 10 −3 .Future more precise numerical studies with FIESTA could potentially fix them unambiguously as rational coefficients accompanying the value of ζ 4 .
Consideration performed in this work will be generalized in a number of avenues.From the point of view of identifying two-dimensional integrable physics of the octagon flux-tube, form factors of super-descendants of the stress-tensor multiplet could provide simpler circumstances for its elucidation since they are sensitive to contribution from single charged flux-tube excitations [64] as opposed to singlet pairs determining the lowest half-BPS operator [15].So far the limitation to form factor observables was solely driven by a lowermultiplicity requirement on the number of external legs in a graph to attain a non-trivial remainder function.It is well-known that in the case of scattering amplitudes, nontrivial remainder functions emerge starting from six legs.Thus, it is important to analyze these in the near mass-shell kinematics introduced in this paper and compare them both functionally as well as from the microscopic stringy point of view.
Regarding development of computational techniques of multi-loop Feynman integrals per se, we are currently capable to break free from the simplifying assumption of the near mass-shell limit and uplift our asymptotically canonical basis for the non-planar graph to the situation of arbitrary virtualities.Solution to the resulting equations is a very different issue though.
A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e x K f B w D X j x G N A 9 I l j A 7 6 U 2 G z M 4 u M 7 N C W P I J X j w o 4 t U v 8 u b f O E n 2 o I k F D U V V N 9 1 d Q S K 4 N q 7 7 7 R T W 1 j c 2 t 4 r b p Z 3 d v f 2 D 8 u F R S 8 e p Y t h k s Y h V J 6 A a B Z f Y N N w I 7 C Q K a R Q I b A f j 2 5 n f f k K l e S w f z S R B P 6 J D y U P O q L H S Q 9 L 3 + u W K W 3 X n I K v E y 0 k F c j T 6 5 a / e I G Z p h N I w Q b X u e m 5 i / I w q w 5 n A a a m X a k w o G 9 M h d i 2 V N E L t Z / N T p + T M K g M S x s q W N G S u / p 7 I a K T 1 J A p s Z 0 T N S C 9 7 M / E / r 5 u a 8 M b P u E x S g 5 I t F o W p I C Y m s 7 / J g C t k R k w s o U x x e y t h I 6 o o M z a d k g 3 B W 3 5 5 l b Q u q t 5 V 9 f K + V q n X 8 j i K c A K n c A 4 e X E M d 7 q A B T W A w h G d 4 h T d H O C / O u / O x a C 0 4 + c w x / I H z + Q P + e 4 2 U < / l a t e x i t > p1 < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 n F x 5 w z L S S 5 m B b p P c 6 j y o 0 r S o A Y = " > A A A B 6 H i c b V D L T g J B E O z F F + I L 9 e h l I j H x R H Y N P o 4 k X j x C I o 8 E N m R 2 6 I W R 2 d l 1 Z t a E E L 7 A i w e N 8 e o n e f N v H G A P C l b S S a W q O 9 1 d Q S K 4 N q 7 7 7 e T W 1 j c 2 t / L b h Z 3 d v f 2 D 4 u F R U 8 e p Y t h g s Y h V O 6 A a B Z f Y M N w I b C c K a R Q I b A W j 2 5 n f e k K l e S z v z T h B P 6 I D y U P O q L F S / b F X L L l l d w 6 y S r y T 6 5 a / e I G Z p h N I w Q b X u e m 5 i / I w q w 5 n A a a m X a k w o G 9 M h d i 2 V N E L t Z / N T p + T M K g M S x s q W N G S u / p 7 I a K T 1 J A p s Z 0 T N S C 9 7 M / E / r 5 u a 8 M b P u E x S g 5 I t F o W p I C Y m s 7 / J g C t k R k w s o U x x e y t h I 6 o o M z a d k g 3 B W 3 5 5 l b Q u q t 5 V 9 f K + V q n X 8 j i K c A K n c A 4 e X E M d 7 q A B T W A w h G d 4 h T d H O C / O u / O x a C 0 4 + c w x / I H z + Q P + e 4 2 U < / l a t e x i t > p1 < l a t e x i t s h a 1 _ b a s e 6 4 = " p 3 V z 8 r Q 1 g e K B 8 n 2 Q 0 V 6 W e 0 C J 2 g o = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l K / T g W v H i s a G u h D W W z 3 b R L N 5 u w O x F K 6 E / w 4 k E R r / 4 i b / 4 b t 2 0 O 2 v p g 4 P H e D D P z g k Q K g 6 7 7 7 R T W 1 j c 2 t 4 r b p 4 h X P w 4 A o a c A t N a A G D I T z D K 7 w 5 0 n l x 3 p 2 P R W v B y W e O 4 Q + c z x 8 A D o 2 V < / l a t e x i t > p2 N a A K D I T z B C 7 w 6 w n l 2 3 p z 3 R W v B y W c O 4 R e c j 2 8 B k o 2 W < / l a t e x i t > p3 < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 n F x 5 w z L S S 5 m B b p P c 6 j y o 0 r S o A Y = " > A A A B 6 H i c b V D L T g J B E O z F F + I L 9 e h l I j H x R H Y N P o 4 k X j x C I o 8 E N m R 2 6 I W R 2 d l 1 Z t a E E L 7 A i w e N 8 e o n e f N v H G A P C l b S S a W q O 9 1 d Q S K 4 N q 7 7 7 e T W 1 j c 2 t / L b h Z 3 d v f 2 D 4 u F R U 8 e p Y t h g s Y h V O 6 A a B Z f Y M N w I b C c K a R Q I b A W j 2 5 n f e k K l e S z v z T h B P 6 I D y U P O q L F S / b F X L L l l d w 6 y S r y M l C B D r V f 8 6 v Z j l k Y o D R N U 6 4 7 n J s a f U G U 4 E z g t d F O N C W U j O s C O p Z J G q P 3 J / N A p O b N K n 4 S x s i U N m a u / J y Y 0 0 n o c B b Y z o m a o l 7 2 Z + J / X S U 1 4 4 0

Figure 1 .
Figure 1.Planar and non-planar double box graphs in the left and right panels, respectively.