Cosmological Perturbation Theory Using the FFTLog: Formalism and Connection to QFT Loop Integrals
Marko Simonović Tobias Baldauf Matias Zaldarriaga John Joseph Carrasco Juna A. Kollmeier

















 
 



.1mps*


 

C[1]>p#1
=8pt














[1]















   1.1

@width









equationsection


    ḍ d Tr
            Ø OD̅ḍ d Tr WD̅Φ̅ϕ̅R̅M̅F̅α̇β̇δ̇γ̇M_ plḍłℓ aḇH̋ GΓ^† C 0k̨qpsv̌xyzPŭn∇k_ NLf_ NL M A I J K
                [4]
                =15.5pt empty
                
                2028 Cosmological Perturbation Theory Using the FFTLog: Formalism and Connection to QFT Loop Integrals
                
                
                
                
                1330 Marko Simonović,^ a Tobias Baldauf,^ b Matias Zaldarriaga,^ a
                
                John Joseph Carrasco^ c and Juna A. Kollmeier^ d
                
                 8pt
                ^ aSchool of Natural Sciences, Institute for Advanced Study, 1 Einstein Drive,
                 Princeton, NJ 08540, United States 
                 2pt
                ^ bDAMTP, Center for Mathematical Sciences,
                Wilberforce Road,
                
                Cambridge, CB3 0WA, United Kingdom
                 2pt
                ^ cInstitute of Theoretical Physics (IPhT), CEA-Saclay and University of Paris-Saclay
                
                F-91191 Gif-sur-Yvette cedex, France
                 2pt
                ^ dThe Observatories of the Carnegie Institution for Science, 813 Santa Barbara St,
                
                Pasadena, CA 91101, United States 7pt
                
                
                
                 Abstract
                
                We present a new method for calculating loops in cosmological perturbation theory. This method is based on approximating a ΛCDM-like cosmology as a finite sum of complex power-law universes. The  decomposition is naturally achieved using an FFTLog algorithm. For power-law cosmologies, all loop integrals are formally equivalent to  loop integrals of massless quantum field theory. These integrals have analytic  solutions in terms of generalized hypergeometric functions. We provide explicit formulae for the one-loop and the two-loop power spectrum and the one-loop bispectrum. A chief advantage of our approach is that the difficult part of the calculation is cosmology independent, need be done only once, and can be recycled for any relevant predictions. Evaluation of standard loop diagrams then boils down to a simple matrix multiplication. We demonstrate the promise of this method for applications to higher multiplicity/loop correlation functions.
                 10pt
                
                
                
                
                § INTRODUCTION
                
                Cosmological perturbation theory (PT) <cit.> offers a complementary approach to predicting and interrogating large-scale structure (LSS) observables in the weakly non-linear regime.  While many investigations in the literature are aimed at understanding and exploiting the highly non-linear regime (e.g. galaxy and cluster formation), it is clear that upcoming experiments (e.g. CMB-S4, DESI, CHIME) will also provide voluminous datasets probing  matter distribution on very large scales. As we continue in this “era of precision cosmology", it is critical to exploit this investment of resources, to the fullest extent possible, in order to deliver on the promise of these upcoming surveys to measure cosmological parameters and potentially probe new physics.  The PT  approach, when recognized as a classical effective field theory <cit.>, has a number of important advantages.  It converges to the correct answer for clustering statistics on large scales as more orders are included.  Its errors are parametrically controlled, so the typical size of the deviations from the correct answer can be easily estimated.  For the scales of relevance, it is rapidly computable compared to the necessarily large and high-resolution cosmological simulations required to attempt adequate comparison of theory and  large-scale data.  Simulation boxes are computed with a single cosmology and thus many computations are required to investigate small changes in cosmological parameters, the effects of cosmic variance, and to perform consistency checks of scheme independence. PT works around these issues entirely.  Indeed it can serve as a spectacular large-scale (IR) complement to simulations, allowing them to focus their power in the incredibly non-linear smaller-scale (UV) regimes where they excel.
                
                
                In the PT approach one treats dark matter and baryons as non-ideal self-gravitating fluids. At early times or on large scales these fluids are nearly homogeneous with small density fluctuations. This allows for the equations of motion to be solved perturbatively. The rigorous foundation of PT as an effective field theory of  large scale structure (EFTofLSS) was recently made <cit.>, although many important results were known for a long time (for a review see <cit.>). One feature of these perturbative solutions is that they convolve initial density fields. Therefore observables, such as correlation functions of density-contrast, or  overdensity, are efficiently written as momentum integrals over a certain number of initial power spectra. These integrals are refereed to as loop integrals. They admit a graphical organization which is why the atomic units of loop contributions to correlation functions are often referred to as loop diagrams.  Higher multiplicity/loop correlation calculations are critical not only to extending the scale of relevance of the analytic approach, but to break degeneracies and optimally extract cosmological parameters from the real data. 
                
                
                While calculating loop integrals is a straightforward task in principle, the computational cost of exact solution becomes prohibitive for higher points (multiplicity) as well as higher loop order corrections. Higher multiplicity kernels quickly become  complicated and every loop brings additional three-dimensional integral. The linear power spectrum that appears in the integrand for real universe applications is known only as a numerical function which makes analytic solution of non-trivial integrals  impossible. Naive numerical integration, to desired precision, on the other hand, quickly becomes slow even with advanced Monte Carlo methods, due to the growth in dimensionality.  This poses a direct challenge to our ability to interrogate large datasets and one that merely more and faster computers will not address.
                
                
                In order to simplify and speed up loop calculations we require new ideas, new strategies, to approach the problem. One inspiring idea, developed in <cit.> and <cit.>, is to use Fast Fourier Transform (FFT) for efficient evaluation of the one-loop power spectrum. After first “deconvolving” the lowest order PT solutions, and performing all angular integrals, the one-loop expressions reduce to a set of simple one-dimensional integrals that can be efficiently evaluated using FFT. Unfortunately, deconvolving higher order perturbative solutions and extending this approach to the one-loop bispectrum or the two-loop power spectrum proves to be challenging <cit.>. 
                
                
                In this paper we build on ideas of <cit.> but choose a slightly different strategy which allows us to go beyond the one-loop power spectrum. Let us briefly sketch the main idea behind our proposal. Prior to doing any integrals, the linear power spectrum is expanded as a superposition of ideal self-similar power-law cosmologies. This is naturally accomplished using FFT in log k. Given some range of wavenumbers of interest, from k_ min to k_ max, the approximation for the linear power spectrum with N sampling points is <cit.>P̅_ lin(k_n) = ∑_m=-N/2^m=N/2 c_m k_n^ν+iη_m , 
                
                where the coefficients c_m and the frequencies η_m are given by
                
                c_m = 1/N∑_l=0^N-1 P_ lin(k_l)  k_l^-νk_ min^-iη_m  e^-2π i m l /N ,   η_m = 2π m/log (k_ max/k_ min) .
                
                Notice that the we denote the approximation for the linear power spectrum with P̅_ lin(k), while eq. (<ref>) uses the exact linear power spectrum P_ lin(k) to calculate the coefficients c_m. We will keep using the same notation throughout the paper. The parameter ν is an arbitrary real number. As we will see, the simplest choice ν=0 is insufficient in some applications, so we will use the more general form of the Fourier transform. In the terminology of <cit.> we call this ν parameter  bias. Note that the powers in the power-law expansion are complex numbers. In practice, even a small number of power-laws, 𝒪(100), is enough to capture all features of the linear power spectrum including the BAO wiggles. One important thing to keep in mind is that the Fourier transform produces the power spectrum that is periodic in log k. Therefore, we will take care to choose k_ min and k_ max such that we cover the range of scales where we actually care about the value of the power spectrum. In other words we are choosing the momentum range where the loop integrals have the most of the support. However, one always has to be careful about possible contributions particularly from high k modes or short scales.  
                
                
                Is this a limitation?  Absolutely not.  At the heart of the EFT understanding is the simple recognition that the PT idealized description of satisfying fluid-like equations of motion can only be valid at certain scales.  This is much the same as the hydrodynamic description of liquid water is only valid at certain scales. Attempting to integrate this approximation over scales outside of its validity introduces non-parametrically controlled errors.  Instead the information in the linear approximation must be supplemented by small-scale UV physics. This data is encoded in physical parameters like speed of sound or viscosity – potentially any dimensionally consistent operators. Such EFT parameters serve two roles.  They must eat up any cutoff-dependence, by definition non-physical, and they must accurately represent the integrating out of any small-scale degrees of freedom.  So baked into the framework that places PT on a rigorous footing is the realization that any integrals of the linear approximation should only be performed over a range of scales consistent with its validity.
                
                
                Notice that the decomposition (<ref>) reduces the evaluation of a loop diagram for an arbitrary cosmology to evaluation of the same diagram for a set of different power-law universes with numerical coefficients. Power-law momentum integrals  can be done analytically. The final answer is a sum of familiar special functions which are straightforward to evaluate. In the simplest case of the one-loop power spectrum, the momentum integral for a power-law universe can be expressed entirely in terms of gamma functions <cit.>. Looking at higher order correlators an interesting pattern emerges. The perturbation theory kernels can always be written such that the general form of loop integrals in a power-law cosmology is formally identical to the one of a massless Quantum Field Theory (QFT) with cubic interactions.[More precisely, it is a QFT in three dimensions with the Euclidean signature.] This is just a formal relationship, but it should prove rather useful. Many results recently developed in the theory of scattering amplitudes can be applied to LSS correlation functions as well. Some steps in this direction have already been made for the one-loop bispectrum in <cit.>. 
                
                
                In this work we derive formulas for the one-loop bispectrum and the two-loop power spectrum in power-law cosmologies which are suitable for effecient numerical evaluation. Generically, the higher multiplicity/loop correlation functions are expressed in terms of the generalized hypergeometric functions. One thing to keep in mind is that the powers ν+iη_m are complex and one has to be careful about the analytic continuation of all results to the entire complex plane.
                
                
                Before diving into the details, let us comment on one important virtue of our method. The decomposition (<ref>) is useful because it separates the cosmology dependent portion, encoded entirely in the coefficients c_m, from the loop calculations which have been reduced to that of much more tractable, ideal, cosmologies. This means that for the fixed value of bias ν, the momentum range (k_ min,k_ max) and the number of sampling points N the difficult part of the calculation which involves momentum integrals can be done only once, saved as a table of numbers and then used for any cosmology. As we will see, the evaluation of the contribution of loop diagrams reduces to a simple (small) matrix multiplication and it is very fast. This opens up the possibility of using our method in Markov chain Monte Carlo parameter estimation.
                
                
                In the rest of the paper we focus on three examples: the one-loop power spectrum, the one-loop bispectrum and the two-loop power spectrum. We will present our calculations in detail and for one-loop diagrams compare them with the standard numerical results. We leave a detailed comparison with the numerical two-loop power spectrum for future work. A  Mathematica notebook used to produce plots is available as an auxiliary file associated with the preprint of the paper on  arXiv. 
                
                
                
                § ONE-LOOP POWER SPECTRUM
                
                Let us first consider the simplest case—the one-loop power spectrum. In perturbation theory there are two different one-loop contributions. Using the usual approximation in which the time dependence is separated from k dependence (for a review see <cit.>), the one-loop power spectrum reads
                
                P_ 1-loop(k,τ) = D^4(τ)[P_22(k)+P_13(k)]  ,
                
                where τ  is conformal time, D(τ) is the growth factor for matter fluctuations and the two terms in the square brackets are given by
                
                P_22(k) = 2 ∫_ F_2^2(, -̨) P_ lin(q) P_ lin( |-̨| )  ,
                
                P_13(k) = 6P_ lin(k) ∫_ F_3(, -, )̨ P_ lin(q)  ,
                
                where ∫_≡∫d^3q/(2π)^3. Diagrammatic representation of these two contributions is shown in Fig. <ref>. The explicit form of kernels F_n can be calculated using well-known recursion relations <cit.>. One important point is that it is always possible to expand kernels in (<ref>) and (<ref>) in integer powers of k^2, q^2 and |-̨|^2. For example,
                
                    F_2(,-̨)     = 5/14 +3 k^2/28 q^2 +3 k^2/28 |-̨|^2-5 q^2/28 |-̨|^2-5 |-̨|^2/28 q^2 +k^4/14 |-̨|^2 q^2.
                
                A similar expression can be found for F_3(,-,)̨.[In the expansion of F_3(,-,)̨ some terms contain |+̨|^2. Given that the kernels are always integrated over , one is allowed to do the following change of coordinates →- and bring these terms to the same form as in (<ref>)] If we further decompose P_ lin(k) in power laws using (<ref>), the one-loop power spectrum becomes a sum of simple momentum integrals of the following form
                ∫_1/q^2ν_1|-̨|^2ν_2≡ k^3-2ν_12(ν_1,ν_2)  ,
                
                where ν_1 and ν_2 are in general complex numbers. 
                
                
                
                
                As we already mentioned, the form of the integral is identical to the one-loop massless two point function in QFT. The only difference is that in this case the powers of the “propagators" are complex numbers rather than integers. Still, the unknown dimensionless function (ν_1,ν_2) can be easily calculated using the standard technique with Feynman parameters. The result is a well known expression <cit.>(ν_1,ν_2) = 1/8π^3/2Γ(32 - ν_1) Γ(32 - ν_2) Γ(ν_12 - 32)/Γ(ν_1) Γ(ν_2) Γ(3 - ν_12) ,
                
                were ν_12=ν_1+ν_2 (throughout the paper we adopt the following notation ν_1...n≡ν_1+⋯ +ν_n). Notice that, thanks to the analytic continuation, (ν_1,ν_2) gives a finite answer even for the values of parameters for which the integral is formally divergent. In practice, breaking the loop calculation into many pieces can lead to some divergent terms. However, as long as the  total sum is well defined and finite, for at least some power-law cosmology P_ lin(k)∼ k^ν, by analytic continuation it is guaranteed that eq. (<ref>) gives the correct answer. 
                
                
                Sometimes the condition that the integral at hand is convergent for at least some power-law power spectrum cannot be met, and one has to use eq. (<ref>) with some care. For example, the function (ν_1,ν_2) vanishes if one of the arguments is zero (or a negative integer). Applying (<ref>) blindly would lead in these cases to paradoxical results. For instance, after power-law decomposition of the linear power spectrum, eq. (<ref>) would imply 
                ∫_0^∞ dq P_ lin(q) = 0  ,
                
                which is obviously the wrong answer. This is a consequence of the well known statement that in dimensional regularization all power-law divergences vanish: ∫_ q^ν = 0.[More precisely, this integral is related to a delta function <cit.>. A  change of coordinates relates
                ∫_1/q^3+2ν_1 = i/2π^2δ(ν_1)  .
                
                To get the consistent results one can use this equation. In practice, there is a much simpler way, as described in the main text. 
                ]
                
                Similar issues can appear in calculating loop diagrams. Luckily, for a ΛCDM-like cosmology, they can be always easily fixed. Let us imagine that the integral we are interested in is divergent for a given bias ν. Then, if the integral diverges in the UV(IR), one has to find the UV(IR) limit of the integrand. This can be easily done fixing all external momenta and sending the loop momentum to infinity(zero). This limit always has the form of eq. (<ref>) and it would be set to zero by dimensional regularization. Therefore, to get the correct answer, one simply has to add the UV(IR) contribution by hand. In the following sections we will give more details for each specific case we consider.  
                
                
                Let us also point out that all UV divergences have a well defined momentum dependence. This momentum dependence is the same as for the counterterms in the EFTofLSS. Therefore, one can proceed without explicitly adding the UV-dependent terms to the loop calculation. The only effect of this choice is to change the usual values of the counterterms. In this sense we can say that eq. (<ref>) calculates only the “finite” part of the loop integral. As expected, the counterterms absorb all UV-dependent pieces. 
                 10pt
                
                
                
                
                 §.§ Symmetries of (ν_1,ν_2) and Recursion Relations
                
                Before moving on to applications, it is instructive to take a closer look at symmetries of the integral in (<ref>). For the one-loop power spectrum this is a straightforward exercise. We use it to introduce some notation and derive a couple of results that will be very useful in more complicated cases, such as the one-loop bispectrum or the two-loop power spectrum. 
                
                
                The most obvious symmetry of the integral (<ref>) is invariance under the shift →-̨. Consequently, the function (ν_1,ν_2) is symmetric in ν_1 and ν_2(ν_1,ν_2) = (ν_2,ν_1)  .
                
                As we will see later, there are similar transformations for higher multiplicity/order diagrams and they always lead to some permutation of parameters ν_i. We are going to call these kind of identities  translation formulas, because they are derived using translations in momentum space. 
                
                
                The integral (<ref>) preserves its form under rescaling of momenta, but this does not lead to any non-trivial condition on (ν_1,ν_2). However, a more complicated rescaling, such as inversion of momenta, does lead to interesting results. For simplicity, let us choose $̨ to be a unit vector:=̨̨̂,̨̂^2=1. Under an inversion→/q^2different factors in the integrand transform in the following way
                
                d^3q → d^3q/q^6 ,    q^2 → q^-2 ,    |̨̂-|^2 → |̨̂-|^2/q^2 .
                  
                Using  these transformations we can write
                (ν_1,ν_2) = ∫_1/q^2ν_1|̨̂-|^2ν_2 = ∫_1/q^2(3-ν_12)|̨̂-|^2ν_2 ,
                
                which immediately implies the  inversion formula
                (ν_1,ν_2) = (3-ν_12,ν_2) .
                
                
                There is one more method to find nontrivial identities for massless loop integrals which is based on the following relation between real and momentum space
                1/q^2ν = Γ(32-ν)/Γ(ν)π^-3/2 2^-2ν∫_1/x^3-2ν e^-i· ,
                
                where∫_≡∫ d^3x. Let us illustrate the main idea behind this method. The starting point is to close all external lines in order to form additional loops. This is equivalent to integrating over all external momenta. For example, for the one-loop power spectrum we can start from
                ∫_1/k^3 (ν_1,ν_2) = ∫_,̨1/q^2ν_1|-̨|^2ν_2k^2(3-ν_12) .
                
                Notice that this expression has a form of a two-loop vacuum diagram. We have chosen to multiply(ν_1,ν_2)with the1/k^3factor such that the whole integral is dimensionless. On the l.h.s. the function(ν_1,ν_2)does not depend onkand the integral trivially reduces to
                ∫_1/k^3 (ν_1,ν_2) = (ν_1,ν_2) 1/(2π)^3∫d^3k/k^3 .
                
                On the r.h.s. we can use (<ref>) and integrate overand$̨. The momentum integrals lead to two delta functions which can be used to do two integrals in real space. At the end of the day, we are left with the following expression
                (ν_1,ν_2) 1/(2π)^3∫d^3k/k^3 = 1/64π^9/2Γ(32 - ν_1) Γ(32 - ν_2) Γ(ν_12 - 32)/Γ(ν_1) Γ(ν_2) Γ(3 - ν_12)∫d^3x/x^3 ,
                
                from which result (<ref>) immediately follows. As we can see, the one-loop two-point function is simple enough that relation (<ref>) is sufficient to fix its form. For higher multiplicity/loop correlation functions this is not the case. The reason is that in those cases the real space integrals are not trivial anymore. However, it turns out that they alway have the same structure as the original momentum integrals. The only difference is that the parameters are shifted: ν_i→ν̃_i≡ 32-ν_i. It follows that there is always an identity which relates two functions with parameters ν_i and ν̃_i. We will refer to these identities as  reflection formulas.
                
                
                Finally, let us present recursion relations which connect functions whose parameters differ by an integer. These relations can be always derived using the fact that the integral of a total derivative vanishes. For instance, 
                ∫_∂/∂ q_i( q_i/q^2ν_1|-̨|^2ν_2) = 0  .
                
                Expanding the derivative we find
                
                (3-2ν_1-ν_2)(ν_1,ν_2)+ν_2[(ν_1,ν_2+1)-(ν_1-1,ν_2+1)]=0  ,
                
                and a similar relation in which ν_1 and ν_2 are exchanged. The importance of these identities is that they relate different terms in the expansion of kernels, such as (<ref>). As we will see in the explicit calculations of P_22 and P_13 diagrams, thanks to the recursion relations all terms in the expansion of the kernels (for fixed ν_1 and ν_2) can be evaluated using a single function (ν_1,ν_2).  
                 10pt
                
                
                
                
                 §.§ Numerical Evaluation of the One-loop Power Spectrum
                
                In this section we will apply eq. (<ref>) to the calculation of the one-loop power spectrum. We will first separately discuss P_22 and P_13 diagrams (see Fig. <ref>).
                P_22 diagram.—Let us begin by reviewing some properties of the P_22 diagram in a cosmology with P_ lin(k)∼ k^ν. In particular, we are interested in finding the powers ν for which the integral is convergent. In order to do that we have to find the asymptotic form of the integrand in the UV and the IR regime. The behavior of the F_2 kernel in these two limits is 
                
                F_2(,-̨) →k/q ,      q→ 0  ,
                
                F_2(,-̨) →k^2/q^2 ,      q→∞ .
                
                It follows that P_22 diagram is convergent if -1<ν<1/2. If we choose bias in FFT to be in this range, then the integral in P_22 is finite for each term in the sum (<ref>) and using (<ref>) we are guaranteed to get the same answer as with the usual numerical evaluation. 
                
                
                Before turning to results, let us write the explicit formula for P_22 diagram. Using (<ref>) and (<ref>) we can write the approximation to the P_22 diagram in the following way 
                P̅_22(k) = 2 ∑_m_1,m_2 c_m_1 c_m_2∑_n_1,n_2 f_22(n_1,n_2) k^-2(n_1+n_2)∫_1/q^2ν_1-2n_1|-̨|^2ν_2-2n_2 .
                
                In this expression c_m_1 and c_m_2 are the coefficients in (<ref>) and n_1 and n_2 are integer powers of q^2 and |-̨|^2 in the expansion of F_2^2(,-̨). Corresponding rational coefficients in this expansion are labeled by f_22(n_1,n_2) and they can be read off from (<ref>). The complex numbers ν_1 and ν_2 are given by
                ν_1 = -12(ν+iη_m_1)    and   ν_2 = -12(ν+iη_m_2)  .
                
                Using the solution for the momentum integral, expression (<ref>) can be further simplified and written in the following way
                P̅_22(k) = k^3∑_m_1,m_2 c_m_1k^-2ν_1· M_22(ν_1,ν_2) · c_m_2k^-2ν_2 ,
                
                where the matrix M_22(ν_1,ν_2) is given by
                
                M_22(ν_1,ν_2)= (32 - ν_12) (12 - ν_12) [ν_1ν_2 (98ν_12^2 - 14ν_12+ 36 ) - 91 ν_12^2 + 3ν_12 + 58]/196 ν_1 (1+ν_1) (12 - ν_1) ν_2 (1+ν_2) (12-ν_2) (ν_1,ν_2).
                
                As we already pointed out, only a single function (ν_1,ν_2) is sufficient to calculate the full diagram. Thanks to the recursion relations (<ref>), all terms from the expansion of F_2 kernels are encoded in the ν-dependent prefactor in matrix M_22(ν_1,ν_2).
                
                
                One can use eq. (<ref>) to calculate the P_22 diagram. The result is shown in Fig. <ref>. As expected, the agreement with the usual numerical integration is excellent. An important thing to notice is that the only cosmology dependence in (<ref>) is in the coefficients c_m. For a given number of sampling points N, bias ν and k_ min and k_ max, the matrix M_22(ν_1,ν_2) is fixed. This means that it can be calculated only once and saved as a table of numbers. The evaluation of the P_22 diagram for an arbitrary cosmology then boils down to doing one FFT to determine coefficients c_m, calculating a vector c_mk^-2ν for each k and a simple matrix multiplication (<ref>).
                
                
                So far we restricted ourselves to biases in the range -1<ν< 12. Outside this range one has to be more careful because the integrals are not convergent anymore and eq. (<ref>) is not guaranteed to give the correct answer. For example, for biases ν<-1 the integrals are divergent in the IR. The leading piece of the P_22 diagram in this limit can be calculated by fixing k and sending → 0 in the integrand. The result is
                
                P_22^ IR (k) = P_ lin(k) k^2 σ_v^2  ,
                
                where σ_v^2≡16π^2∫_0^∞ dq P_ lin(q). The integral in σ_v^2 would be set to zero by eq. (<ref>) and missed in the final answer. Therefore, to get a correct result, one simply has to add P_22^ IR (k) to eq. (<ref>) at the end of the calculation. Notice that we have kept only the leading IR divergence, which is enough for biases in the range -3<ν<-1. If the bias was even smaller, one would have to keep track of sufficient number of subleading IR divergences. Similarly, when 12<ν<32, the leading UV divergence that has to be added on the r.h.s. of eq. (<ref>) to get the correct result is 
                
                P_22^ UV (k) = 9/196π^2k^4 ∫_0^∞ dq P_ lin(q)/q^2 . 
                
                However, for very high or low values of biases, the momentum range (k_ min,k_ max) has to be very wide for the integrals to converge to correct values. This implies a large number of frequencies and it is not practical.
                P_13 diagram.—The asymptotic behavior of the F_3 kernel in the UV and the IR limit is the same 
                
                F_3(, -, )̨→k^2/q^2 ,      q→ 0    or   q→∞ .
                
                Consequently, P_13 diagram is divergent in the UV for ν>-1 and divergent in the IR for ν<-1. In other words, P_13 diagram is never finite in a power-law cosmology. For ν>-1, the only possible mismatch between the true answer and eq. (<ref>) comes from the UV part of the integral. For fixed k and taking the limit →∞
                P_13^ UV (k) = -61/105 P_ lin(q) k^2 σ_v^2 .
                
                On the other hand, for ν<-1, the possible error comes from the IR limit[Notice that P_13 diagram has two IR divergences →0 and →$̨, which are combined in a single expression.
                P_13^ IR (k) = -P_ lin(k) k^2 σ_v^2  .
                
                As we already mentioned, in dimensional regularization this type of the integrals would be set to zero by (<ref>). This implies that to get the correct values forP_13, depending on the choice of bias, we have to add either the UV or IR term by hand. Notice that here we are writing down only the leading IR and UV parts of theP_13diagram. For biasesν<-3orν>1, one would have to include the corresponding subleading terms as well.
                
                
                    < g r a p h i c s >
                
                    < g r a p h i c s >
                Two contributions to the one-loop power spectrum calculated using direct numerical integration and eq. (<ref>) and eq. (<ref>). Both plots are produced using ν=-1.6, N=150, k_ min=3· 10^-4 h Mpc^-1 and k_ max=180 h Mpc^-1. For this value of bias both P_22 and P_13 are very different from their standard values.
                
                Let us see how the formulas above work in practice. With the same notation as for theP_22diagram, we can write 
                P̅_13(k) = 6P_ lin(k) ∑_m_1 c_m_1∑_n_1,n_2 f_13(n_1,n_2) k^-2(n_1+n_2)∫_1/q^2ν_1-2n_1|-̨|^-2n_2 .
                
                Solving the momentum integral, this expression can be further simplified 
                P̅_13(k) = k^3 P_ lin(k) ∑_m_1 c_m_1k^-2ν_1· M_13(ν_1) ,
                
                where the vectorM_13(ν_1)is given by
                
                M_13(ν_1) = 1 + 9 ν_1/4tan(ν_1 π)/28 π (ν_1 + 1) ν_1 (ν_1 - 1) (ν_1 - 2) (ν_1 - 3) .
                
                Notice that to eq. (<ref>) one has to add the UV or the IR part of the integral. For example, forν>-1, we plot the result in Fig. <ref>. As expected, onceP_13^ UV(k)is added to eq. (<ref>), the agreement with the usual numerical result is excellent.
                 The full one-loop power spectrum.—So far we were trying to reproduce the usual numerical results for separate pieces of the one-loop power spectrum. However, only their sum is a well defined observable. Thanks to the Equivalence Principle the IR divergences cancel and the  total one-loop power spectrum is convergent for the range of power laws-3<ν<-1<cit.>. This means that with the choice of bias in this range, the formulas above should lead to the correct answer forP_ 1-loop(k), without having to deal with the IR divergences explicitly. In Fig. <ref> we plot the one-loop power spectrum calculated in this way and show that our method indeed agrees with the usual numerical result. As expected, the separate termsP_13andP_22are rather different from their usual values (see Fig. <ref>). However, the “mistake” that eq. (<ref>) makes in assigning some finite values to divergent integrals has to cancel between the two contributions in the same way the IR divergences cancel. Indeed, the IR limit of theP_22diagram is exactly the same as the IR limit ofP_13diagram, but with the opposite sign (see eq. (<ref>) and eq. (<ref>)).
                
                
                
                    < g r a p h i c s >
                
                The full one-loop power spectrum calculated summing up contributions from Fig. <ref>. 
                
                In conclusion, the one-loop power spectrum can be easily calculated using decomposition (<ref>). For practical applications, the most optimal choice of bias is in the range close to zero-0.5<ν<0because it requires the least number of frequencies to reproduce the linear power spectrum on relevant scales. For this range theP_22diagram can be evaluated directly applying our method. To get the correctP_13diagram, one has to addP_13^ UVterm to the r.h.s of eq. (<ref>). 
                 10pt
                
                
                
                
                 §.§ One-loop Power Spectrum of Biased Tracers
                
                The method described above can be also applied to the one-loop power spectrum of biased tracers <cit.> (for a review see <cit.>). In this section we give explicit formulas for all relevant one-loop contributions. We will follow the notation of <cit.>.
                
                
                The density contrast of biased tracers, such as the dark matter halosδ_h, is a local function of the underlying dark matter field. The functional dependence is expressed through all possible operators built from the tidal tensor∂_i∂_jΦ(and its derivatives), whereΦcan be either gravitational potentialΦ_gor velocity potentialΦ_v. These two potentials are the same at leading order in perturbation theory but starting from second order they differ. To calculate the one-loop power spectrum of biased tracers one has to keep in bias expansion all operators up to third order in perturbation theory[Notice that <cit.> is using b_δ_2=b_22 and b_δ_3=b_36.]δ_h = b_1 δ +b_2/2δ^2 +b__2_2 + b_3/6δ^3 + b__3_3 + b_(_2δ)_2δ + b_Γ_3Γ_3 .
                
                The operators_2,_3andΓ_3are defined as
                
                    _2(Φ)    ≡ (∂_i∂_jΦ)^2-(∂^2Φ)^2  , 
                    _3(Φ)    ≡ -∂_i∂_j Φ ∂_j∂_kΦ ∂_k ∂_i Φ -1/2 (∂^2Φ)^2 +3/2 (∂_i∂_jΦ)^2∂^2Φ , 
                    Γ_3(Φ_g,Φ_v)    ≡_2(Φ_g) - _2(Φ_v)   .
                
                However, only four renormalized operators contribute to the one-loop power spectrum. These areδ,[δ^2],[_2]and[Γ_3]. The final answer is given in terms of four corresponding renormalized bias parameters and six different momentum integrals <cit.>
                    < g r a p h i c s >
                
                    < g r a p h i c s >
                
                    < g r a p h i c s >
                
                    < g r a p h i c s >
                Four different contributions to the one-loop power spectrum of biased tracers. All plots are produced using ν=-1.6, N=150, k_ min=10^-5 h Mpc^-1 and k_ max=5 h Mpc^-1. For these values of parameters the difference with respect to the usual numerical calculation is less than 0.1% at all scales.
                    P_h(k,τ) =    b_1^2(P_ lin(k,τ)+ P_ 1-loop(k,τ))  
                        + b_1 b_2  ℐ_δ^2(k,τ) +  2b_1 b__2 ℐ__2(k,τ) + ( 2b_1 b__2 + 4/5 b_1 b_Γ_3) ℱ__2(k,τ) 
                        + 1/4 b_2^2  ℐ_δ^2δ^2(k,τ) + b__2^2  ℐ__2_2 (k,τ) + 1/2 b_2b__2 ℐ_δ_2_2(k,τ)  .
                
                In principle, at this order in perturbation theory one has to add higher derivative operators such as∂^2δ. However the contribution from this operator is trivial and it does not lead to a loop integral.
                
                The time dependence of all momentum integrals inP_h(k,τ)isD(τ)^4. The explicitk-dependences are 
                
                    ℐ_δ^2(k)     = 2∫_ F_2(,-̨)P_ lin(q)P_ lin(|-̨|)  ,    (-3<ν<- 12) 
                    ℐ__2(k)     = 2 ∫_σ^2(,-̨)F_2(,-̨) P_ lin(q)P_ lin(|-̨|)  ,    (-3<ν< 12) 
                    ℱ__2(k)     = 4 P_ lin(k) ∫_σ^2(,-̨)F_2(,̨-) P_ lin(q)  ,    (-3<ν<-1) 
                    ℐ_δ^2δ^2(k)    = 2 ∫_P_ lin(q)P_ lin(|-̨|)  ,    (-3<ν<- 32) 
                    ℐ__2_2 (k)     = 2 ∫_ (σ^2(,-̨))^2 P_ lin(q)P_ lin(|-̨|)  ,    (-3<ν< 12) 
                    ℐ_δ_2_2(k)    =  2 ∫_σ^2(,-̨) P_ lin(q)P_ lin(|-̨|)  ,    (-3<ν<- 12 )
                
                whereσ^2(_̨1,_̨2)=(_̨1·_̨2/k_1k_2)^2-1. For each term we give a range of power laws for which the integral is convergent. Following the same steps as in the case of the one-loop power spectrum of matter fluctuations, we find that matrices analogous toM_22andM_13are given by 
                
                    M_ℐ_δ^2(ν_1,ν_2)    = (3-2ν_12)(4-7ν_12)/17ν_1ν_2 (ν_1,ν_2) , 
                    
                    M_ℐ__2(ν_1,ν_2)    = -(3-2ν_12)(1-2ν_12)(6+7ν_12)/28ν_1(1+ν_1)ν_2(1+ν_2) (ν_1,ν_2) , 
                    
                    M_ℱ__2(ν_1)    = -15tan(ν_1π)/28π(ν_1+1)ν_1(ν_1-1)(ν_1-2)(ν_1-3) , 
                    
                    M_ℐ_δ^2δ^2(ν_1,ν_2)    = 2 (ν_1,ν_2) , 
                    
                    M_ℐ__2_2(ν_1,ν_2)    = (3-2ν_12)(1-2ν_12)/ν_1(1+ν_1)ν_2(1+ν_2) (ν_1,ν_2) , 
                    
                    M_ℐ_δ^2_2(ν_1,ν_2)    = 3-2ν_12/ν_1ν_2 (ν_1,ν_2) .
                
                In Fig. <ref> we plot some of the shapes and compare our method with the standard numerical evaluation.  Notice that forℐ_δ^2δ^2(k)shape we subtract the constant shot-noise part and plot just the differenceℐ_δ^2δ^2(k)-ℐ_δ^2δ^2(0). This difference is convergent even forν>- 32. 
                
                
                One important point to make is that the full one-loop power spectrum of biased tracers requires only a single function(ν_1,ν_2)with a single bias in the range-3<ν<- 32. This range can be extended to higher biases by adding the corresponding UV parts of the integrals in the same way as for the matter power spectrum. 
                 10pt
                
                
                
                § ONE-LOOP BISPECTRUM
                
                In perturbation theory there are four different diagrams that contribute to the one-loop bispectrum and their sum can be schematically written like <cit.>
                B_ 1-loop(_̨1,_̨2,_̨3,τ) = D^4(τ)[B_222+B_321^I+B_321^II+B_411] .
                
                From translational invariance it follows that_̨1+_̨2+_̨3 = 0. The individual terms in square brackets are given by the following integrals
                
                    B_222 = 8 ∫_ F_2 (,_̨1-)     F_2(_̨1-,_̨2+) F_2(_̨2+,-) 
                       × P_ lin(q) P_ lin(|_̨1-|) P_ lin(|_̨2+|)  ,
                
                B_321^I = 6 P_ lin(k_1) ∫_ F_3(,_̨2-,-_̨1) F_2(,_̨2-) P_ lin(q) P_ lin(|_̨2-|) + 5  perms ,
                
                B_321^II = F_2(_̨1,_̨2) P_ lin(k_1) P_13(k_2) + 5  perms ,
                
                B_411 = 12 P_ lin(k_1) P_ lin(k_2) ∫_ F_4(,-,-_̨1,-_̨2) P_ lin(q) + 2  cyclic  perms .  
                
                The diagrammatic representation of all these contributions is shown in Fig. <ref>.
                
                
                    < g r a p h i c s >
                
                Diagrammatic representation of four contributions to the one-loop bispectrum.
                
                To evaluate the one-loop bispectrum we can follow the same steps as for the one-loop power spectrum. After expanding the kernels and decomposing the linear power spectrum in power laws, all terms in the sums are proportional to the integral of the following form <cit.>∫_1/q^2ν_1|_̨1-|^2ν_2 |_̨2+|^2ν_3≡ k_1^3-2ν_123 (ν_1,ν_2,ν_3;x,y)  ,
                
                wherex≡ k_3^2/k_1^2andy≡ k_2^2/k_1^2. Parametersν_1,ν_2andν_3have the same form as before (see (<ref>)). The overall scaling of the integral with momentum is fixed and here we choose to express that scaling in terms ofk_1. The rest defines a function(ν_1,ν_2,ν_3;x,y)which depends only on the ratiosxandy. Triangle inequality|k_2-k_3|≤ k_1≤ k_2+k_3implies that the physically  allowed region in(x,y)plane is given by inequalities|√(x)-√(y)|≤ 1and√(x)+√(y)≥ 1and we will focus on evaluating the function in this region (see Fig. <ref>). Before giving the explicit expression suitable for numerical evaluation we present some important symmetry properties of(ν_1,ν_2,ν_3;x,y)which can be derived from its integral representation.
                
                
                
                
                    < g r a p h i c s >
                
                 Left panel: Domain of (ν_1,ν_2,ν_3;x,y) allowed by the triangle inequality. Six different regions correspond to six permutations of external momenta. All bispectrum configurations can be evaluated focusing on one of these regions. For example, in this paper we choose x≤ y≤ 1 which corresponds to k_3≤ k_2≤ k_1. Solid lines split the domain in three different parts that will be relevant for the evaluation of the two-loop power spectrum (see Section <ref>).  Right panel: The same as left panel in more conventional variables k_3/k_1 and k_2/k_1. 10pt
                
                
                
                 §.§ Symmetries of (ν_1,ν_2,ν_3;x,y) and Recursion Relations
                
                As in the case of the one-loop power spectrum, the simplest identities follow from shifts of the momentum. There are two basic translation formulas for function(ν_1,ν_2,ν_3;x,y). The first one follows from→_̨1-and reads 
                (ν_1,ν_2,ν_3;x,y) = (ν_2,ν_1,ν_3;y,x)  .
                
                If we do a different shift,→ - _̨2, we get 
                (ν_1,ν_2,ν_3;x,y)=x^3/2-ν_123(ν_3,ν_2,ν_1;1x,yx)  .
                
                These two formulas are sufficient to generate identities involving all six permutations of parametersν_1,ν_2andν_3. These are
                (ν_1,ν_2,ν_3; x,y)=   (ν_2,ν_1,ν_3; y,x)
                
                =   x^3/2-ν_123(ν_3,ν_2,ν_1; 1x,yx)=x^3/2-ν_123(ν_2,ν_3,ν_1; yx,1x)
                
                =   y^3/2-ν_123(ν_3,ν_1,ν_2; 1y,xy)=y^3/2-ν_123(ν_1,ν_3,ν_2; xy,1y)  .
                 
                An intuitive way to understand these symmetries is to realize that they map_̨1,_̨2and_̨3into each other, preserving the shape of the triangle. The six equations then correspond to nothing other but six possible permutations of three external momenta. From another point of view, for evaluation of the bispectrum we can always choose a small “corner” in the(x,y)plane (see Fig. <ref>). The choice that we make in this paper isx≤ y≤ 1which corresponds to the following ordering of momenta:k_3≤ k_2≤ k_1. 
                
                
                Let us now derive the inversion formula for(ν_1,ν_2,ν_3;x,y). We can start from 
                (ν_1,ν_2,ν_3;x,y) = ∫_1/q^2ν_1|̨̂_1-|^2ν_2 |√(y)̨̂_2+|^2ν_3 ,
                
                wherę̂_1and̨̂_2are unit vectors. Under inversion→/q^2, apart from transformations described in (<ref>), we also get 
                
                |√(y)̨̂_2+|^2  →  y/q^2 |̨̂_2/√(y) +  |^2  .
                
                The whole integral then changes to 
                (ν_1,ν_2,ν_3;x,y) = ∫_y^-ν_3/q^2(3-ν_123)|̨̂_1-|^2ν_2 |̨̂_2/√(y)+|^2ν_3 .
                
                It is easy to read off the inversion formula from this expression. One only has to keep in mind that, due to√(y)appearing in the denominator in the last term, the arguments of the function change to xyand 1y. We finally get
                (ν_1,ν_2,ν_3;x,y) =  y^-ν_3(3-ν_123,ν_2,ν_3; xy, 1y )  . 
                
                Combining this inversion formula with translation formulas (<ref>), we get a set of three identities for functions evaluated at the  same point(x,y)(ν_1,ν_2,ν_3;x,y)=     x^3/2-ν_23 (ν_3,3-ν_123,ν_1;x,y)
                
                =     x^3/2-ν_23  y^3/2-ν_13 (ν_2,ν_1,3-ν_123;x,y) 
                
                =     y^3/2-ν_13 (3-ν_123,ν_3,ν_2;x,y)  .
                
                
                Whenν_123=3the previous equation implies 
                (ν_1,ν_2,3-ν_12) = x^3/2-ν_23  y^3/2-ν_13 (ν_1,ν_2)  .
                
                This expression is sometimes referred to as the  star-triangle duality. The conditionν_123=3does not correspond to a generic situation but it can be used to derive another generic formula for(ν_1,ν_2,ν_3;x,y). The idea is to “split” one of the parameters using the one-loop integral, such thatν_123=3is satisfied. For example, we can choose to write 
                1/|_̨2+|^2ν_3 = 1/(ν_123- 32, 3-ν_12)∫_1/s^2ν_123-3|_̨2+-|^6-2ν_12 .
                
                Now one can start with the integral representation of(ν_1,ν_2,ν_3;x,y)and use the previous formula. Notice that we have chosen parameters such that the integration incan be then easily done using (<ref>). The remaining integral inhas again the form of the one-loop bispectrum. Following these steps one derives the star-triangle formula 
                (ν_1,ν_2,ν_3;x,y) = Γ(ν̃_1)Γ(ν_1)Γ(ν̃_2)Γ(ν_2)Γ(ν̃_3)Γ(ν_3)Γ(3-ν̃_123)Γ(3-ν_123) (ν̃_2,ν̃_1,3-ν̃_123 ;x,y)  ,
                
                whereν̃_i = 32-ν_i. Finally, using the method described after eq. (<ref>), it is possible to derive the following reflection formula 
                (ν_1,ν_2,ν_3;x,y) =  Γ(ν̃_1)Γ(ν_1)Γ(ν̃_2)Γ(ν_2)Γ(ν̃_3)Γ(ν_3)Γ(3-ν̃_123)Γ(3-ν_123) x^3/2-ν_23 y^3/2-ν_13(ν̃_1,ν̃_2, ν̃_3 ; x,y)  .
                 
                This is not an independent relation because it follows from eq. (<ref>) and eq. (<ref>). 
                
                
                For arbitrary choice ofνeq. (<ref>) and eq. (<ref>) are not very useful, because they relate two functions with two  different biases. However, there are some special choices ofνfor which this is not the case. Let us remember that, up to an integer, the structure of parameters is 
                ν_i = -ν/2 - i η_m_i/2 ,
                
                which implies 
                ν̃_i = 3/2 + ν/2 + i η_m_i/2 .
                
                Up to an integer, the transformationν→ν̃does not change the bias ifνis a odd integer multiple of- 12. For example, let us imagine thatν=-32. In this case
                ν_i = 3/4 - i η_m_i/2    ⇒    ν̃_i = 3/4 + i η_m_i/2 .
                
                In other words, forν=-32,ν̃_iis just a complex conjugate ofν_i. For example, assumingν=-32, the reflection formula becomes
                
                    (ν_1+n_1,    ν_2+n_2,ν_3+n_3 ;x,y) =  Γ(ν̃_1-n_1)Γ(ν_1+n_1)Γ(ν̃_2-n_2)Γ(ν_2+n_2)Γ(ν̃_3-n_3)Γ(ν_3+n_3)Γ(3-ν̃_123+n_123)Γ(3-ν_123-n_123)
                       × x^3/2-ν_23-n_23 y^3/2-ν_13-n_13^*(ν_1-n_1,ν_2-n_2, ν_3-n_3 ; x,y)  ,
                
                wheren_iare integers coming from the expansion of the kernels. This equation provides a simple relation between two functions in the same point and with the same parameters, but with the opposite sign of integer part ofν_i. For example, such pairs of functions do exist in the expansion of kernels inB_222and each of them can be calculated for the price of a single evaluation. Similar identity can be derived forν=- 12and the same conclusions apply to a set of formulas (<ref>). 
                
                
                Let us conclude showing that the function(ν_1,ν_2,ν_3;x,y)satisfies a set of recursion relations <cit.>. As in the case of the one-loop power spectrum we can start from the following identity
                ∫_∂/∂ q_i( q_i/q^2ν_1|_̨1-|^2ν_2 |_̨2+|^2ν_3) = 0  .
                
                After expanding the derivative in the integrand, the previous equation can be rewritten as 
                
                    ν_2  (ν_1,ν_2+1,ν_3) + ν_3y (ν_1,ν_2,ν_3+1) = 
                          (ν_1+ν_123 -3) (ν_1,ν_2,ν_3) + ν_2  (ν_1-1,ν_2+1,ν_3) +ν_3  (ν_1-1,ν_2,ν_3+1) ,
                
                where we suppressed the(x,y)argument in all functions to avoid clutter. There are two other similar expressions that can be derived replacingq_iin the numerator of the integral in (<ref>) with(_̨1-)_ior(_̨2+)_i. These are
                
                    ν_1  (ν_1+1,ν_2,ν_3) + ν_3x (ν_1,ν_2,ν_3+1) = 
                          (ν_2+ν_123 -3) (ν_1,ν_2,ν_3) + ν_1  (ν_1+1,ν_2-1,ν_3) +ν_3  (ν_1,ν_2-1,ν_3+1) , 
                       ν_1y  (ν_1+1,ν_2,ν_3) + ν_2 x (ν_1,ν_2+1,ν_3) = 
                          (ν_3+ν_123 -3) (ν_1,ν_2,ν_3) + ν_1  (ν_1+1,ν_2,ν_3-1) +ν_2  (ν_1,ν_2+1,ν_3-1) .
                
                Notice that we wrote these equations such that the sum of the arguments in each function on the l.h.s. isν_123+1and  the sum of the arguments in each function on the r.h.s. isν_123. This splitting suggests the following interpretation of the recursion relations.
                They can be thought of as a system of three linear equations where three unknown functions are those in which one of the parameters is increased by 1. If we denote the r.h.s. of previous equations withA_1,A_2andA_3A_1≡     (ν_1+ν_123 -3) (ν_1,ν_2,ν_3) + ν_2  (ν_1-1,ν_2+1,ν_3) +ν_3  (ν_1-1,ν_2,ν_3+1)  , 
                    
                    A_2≡     (ν_2+ν_123 -3) (ν_1,ν_2,ν_3) + ν_1  (ν_1+1,ν_2-1,ν_3) +ν_3  (ν_1,ν_2-1,ν_3+1) ,
                    
                    A_3≡     (ν_3+ν_123 -3) (ν_1,ν_2,ν_3) + ν_1  (ν_1+1,ν_2,ν_3-1) +ν_2  (ν_1,ν_2+1,ν_3-1) ,
                
                then the solution of the system is given by <cit.>
                    (ν_1+1,ν_2,ν_3)    = 1/2ν_1y(-A_1 x + A_2 y + A_3 )  , 
                    (ν_1,ν_2+1,ν_3)    = 1/2ν_2 x(A_1 x - A_2 y + A_3 )  , 
                    (ν_1,ν_2,ν_3+1)    = 1/2ν_3 x y(A_1 x + A_2 y - A_3 )  .
                
                In other words, seven different functions whose parameters live on the planeν_123= const.determine three extra integrals on the planeν_123+1= const. These identities are very useful. For example, in the expansion of kernels inB_222they reduce the number of independent terms by roughly a factor of 2 (from 72 to 38). They are also very important for simplifying the two-loop calculation as we are going to see in the following sections. 
                 10pt
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                
                 §.§ Evaluation of (ν_1,ν_2,ν_3;x,y)
                
                After making these general remarks based on the integral representation, let us turn to the explicit expression for(ν_1,ν_2,ν_3;x,y). Unlike(ν_1,ν_2), this function cannot be simply expressed in a closed form in terms of gamma functions. Starting from (<ref>) and using Feynman parameters we get (see Appendix <ref>)
                
                    (ν_1,ν_2, ν_3 ;x,y)     = 1/8π^3/2Γ(ν_123-3/2)/Γ(ν_1)Γ(ν_2)Γ(ν_3)
                         ×∫_0^1 du ∫_0^1 dv u^ν_1-1(1-u)^ν_2-1v^1/2-ν_3(1-v)^ν_3-1/( u v (1-u) + u(1-v) y + (1-u) (1-v) x )^ν_123-3/2 .
                
                The expression on the r.h.s. belongs to the class of hypergeometric functions of two variables. In particular,(ν_1,ν_2, ν_3 ;x,y)can be written as a linear combination of AppellF_4functions <cit.>. These special functions can be evaluated using their series representations. The region of convergence is given by√(x)+√(y)< 1, which unfortunately covers only the unphysical part of the(x,y)plane. As usual, this kind of problems can be circumvented by performing the analytic continuation. This can be done in several ways, depending on the region of parameter space that one wants to cover <cit.>. Although all results are formally equivalent and  can be related to each other, distinct expressions can be very different from the point of view of practical calculation. A series representation of(ν_1,ν_2, ν_3 ;x,y), optimized for numerical evaluation of the bispectrum, is given by the following formula 
                
                    (ν_1,ν_2,    ν_3 ;x,y) = (πν_23)/8√(π)Γ(ν_1) Γ(ν_2) Γ(ν_3) Γ(3-ν_123)
                       [ x^3/2-ν_23∑_n=0^∞  a_n(ν_1,ν_2,ν_3)·  x^n _2F_1(ν_1+n,  32-ν_2+n, 3-ν_23+2n, 1-y ) . 
                       . -y^3/2-ν_13∑_n=0^∞ a_n(ν̃_1,ν̃_2,ν̃_3)·  x^n _2F_1(ν̃_1+n,  32-ν̃_2+n, 3-ν̃_23+2n, 1-y ) ]  ,
                
                where
                
                a_n(ν_1,ν_2,ν_3)=Γ(ν_1+ n)Γ(3 -ν_123 + n)/Γ(5/2-ν_23+ n)n!Γ(3/2-ν_3+n)Γ(3/2-ν_2+n)/Γ(3-ν_23+2n) .
                
                The functions_2F_1(…,1-y)that appear in the result are standard Gauss hypergeometric functions. We review their definition and some important properties in Appendix <ref>. The derivation of eq. (<ref>) is given in Appendix <ref>.  
                
                
                Let us make some comments about expression (<ref>). The first thing to notice is that the series is always convergent if we restrict ourselves to the regionx≤ y≤ 1. The minimal allowed value of1-yin the given region is34, which corresponds to folded triangles. The smaller1-yis, the easier it gets to calculate the hypergeometric functions using their power series representation. For higher values ofyand smaller values ofxthe convergence is very fast. In the limits ofx→ 0andy→ 1,  corresponding  to squeezed triangles, only a few terms need be kept in the sum. The slowest convergence is for high values ofx. The limiting case isx=1andy=1, which corresponds to equilateral triangles. Even in this case only a relatively modest number of terms,∼𝒪(50), need be kept in the sum to reach satisfactory precision.
                
                
                Another important point to keep in mind is that dependences onxandyare explicitly separated in our formula. Furthermore, thex-dependence is trivial. This means that in practice, for a giveny, calculation for anyk_1andxcan be done evaluating the hypergeometric functions only once. This can speed up any full bispectrum calculation significantly. 
                
                
                There are additional optimizations which can exploit many well-known properties of hypergeometric functions. One such property is that_2F_1(a+n_1,b+n_2,c+n_3,z)for any set of integers(n_1,n_2,n_3)can be always written as a linear combination of just two hypergeometric functions such as, for example,_2F_1(a,b,c,z)and_2F_1(a+1,b,c,z). Using this property one can prove the following recursion relation[Notice that in equation (10) of reference <cit.> there is a typo. The sign between the two terms in square brackets should be + instead of -.]<cit.>
                f_n-1 = f_n +c(1-a-b-2n)+2ab-2n(n-1)(c+2n-2)(c+2n) z   f_n - (a+n)(b+n)(c-a+n)(c-b+n)(c+2n)^2((c+2n)^2-1)z^2 f_n+1 ,
                  
                wheref_n=_2F_1(a+n,b+n,c+2n,z). This equations gives a way to recursively calculate all hypergeometric functions in the power series of (<ref>). In practice one should exercise some caution, as certain values of arguments in this recursion relation can be numerically unstable. 
                
                
                In some special cases eq. (<ref>) further simplifies. For example, for the case of isosceles triangles,y=1and all hypergeometric functions are equal to one. As expected from the symmetry properties in (<ref>), the result in this case becomes symmetric inν_2andν_3.
                
                
                Finally, let us point out that if one of the parameters is a negative integer or zero, the sum in (<ref>) truncates. To see this explicitly, considerν_1=-N, withN≥0a non-negative integer. In the limitν_1→-N,1/Γ(ν_1)in the normalization goes to zero. If there were no terms in the sum which diverge in the same limit, the result would be zero. The hypergeometric functions are always regular. Therefore, we have to look at the coefficients. By inspection we see that alla_n(ν̃_1,ν̃_2,ν̃_3)coefficients are regular as well. The only divergence comes fromΓ(ν_1+n)in the coefficienta_n(ν_1,ν_2,ν_3), forn≤ N. Given that 
                Γ(-N+n)/Γ(-N) = (-1)^n Γ(N+1)/Γ(N-n+1) , 
                
                we can rewrite the final answer in the following way
                
                    (-N,ν_2,    ν_3 ;x,y) = (-1)^N+1√(π)(πν_23)(πν_3)/8 Γ(ν_2) Γ(ν_3) Γ(3+N-ν_23)
                       ∑_n=0^N ∑_m=0^N-nN! (-1)^m+n/(N-m-n)!Γ( 32-ν_2+n+m)/Γ(5/2-ν_23+ n)Γ(ν_3 - N - 12+m)x^3/2-ν_23+n/n!y^m/m! .
                
                In conclusion, when one of the arguments of(ν_1,ν_2,ν_3 ;x,y)is a negative integer or zero, the function becomes a simple polynomial inxandyof degreeN. If two of the arguments are negative integers or zero, then the function vanishes.
                 10pt
                
                
                
                 §.§ Numerical Evaluation of the One-loop Bispectrum
                
                Let us now turn to the numerical evaluation of the bispectrum. We will consider each term in (<ref>) separately. 
                
                
                    < g r a p h i c s >
                
                Equilateral B_222 diagram as a function of k. The calculation is done with bias ν=-0.25 and N=50 sampling points for the power spectrum. This is enough to reach the sub-percent precision.B_222 term.—We begin with theB_222contribution. We will fist find the range of biases for which the integral is convergent. Two of threeF_2kernels, have the same structure and asymptotic behavior as in theP_22diagram. The third kernel tends to𝒪(1)constant in the→ 0limit. In the UV limit
                
                F_2(_̨1-,_̨2+) →k^2/q^2 ,      q→∞ .
                
                Combining all these limits it follows that the integral inB_222diagram is convergent for power laws in the range-1<ν<1. Therefore, choosing a bias close to zero, we expect our method to reproduce the results of the usual numerical integration. After expanding the kernels and linear power spectra in power laws, we can write the result as a matrix multiplication
                B̅_222(k_1,k_2,k_3) = k_1^3∑_m_i M_222(ν_1,ν_2,ν_3;x,y) · c_m_1k_1^-2ν_1· c_m_2k_1^-2ν_2· c_m_3k_1^-2ν_3 ,
                
                where the matrixM_222(ν_1,ν_2,ν_3;x,y)is given by
                
                M_222(ν_1,ν_2,ν_3;x,y) = 8∑_n_i f_222(n_1,n_2,n_3;x,y) (ν_1-n_1,ν_2-n_2,ν_3-n_3;x,y)  .
                
                As before,c_m_iare the coefficients of the expansion (<ref>) andn_1,n_2andn_3are integer powers ofq^2,|_̨1-|^2and|_̨2+|^2in the expansion of the three kernels in the integrand. The rational coefficients in this expansion are labeled byf_222(n_1,n_2,n_3;x,y). We give the explicit expression forM_222in Appendix <ref>. 
                
                
                Let us make a couple of comments about the formulas above. As in the case of the one-loop power spectrum, the evaluation of the bispectrum boils down to a simple matrix multiplication. The matrixM_222is cosmology independent, so it has to be calculated only once. Notice thatM_222depends only on the shape of the triangle formed by the three external momenta, and not on its absolute size. It is clear from eq. (<ref>) that for fixedxandyone can calculate all triangles with differentk_1using the sameM_222. The size of this matrix isN^3, whereNis the number of sampling points of the linear power spectrum. In practice, in order to reach the sub-percent precision, it is enough to use𝒪(50)points. We plot theB_222diagram in Fig. <ref> and as expected it is in a very good agreement with the usual numerical evaluation. 
                
                
                    < g r a p h i c s >
                
                Equilateral B_321^I diagram as a function of k. The calculation is done with bias ν=-0.25 and N=50 sampling points for the power spectrum. This is enough to reach the sub-percent precision.B_321^I term.—Let us now turn toB_321^Idiagram. The asymptotic behavior of theF_3kernel in the integrand is the same in the UV and the IR 
                
                F_3(,_̨2-,-_̨1) →k/q     q→∞   and   q→ 0  .
                
                Therefore, the integral is convergent for power laws in the range-1<ν<0. For small negative bias we can use our method without dealing with the possible UV divergences. As in the previous case we can write
                B̅_321^I(k_1,k_2,k_3) = k_1^3 P_ lin(k_1) ∑_m_i M_321(ν_1,ν_2;x,y) · c_m_1k_1^-2ν_1· c_m_2k_1^-2ν_2 + 5  perms ,
                
                where the matrixM_321(ν_1,ν_2;x,y)is given by
                
                M_321(ν_1,ν_2;x,y) = 6∑_n_i f_321(n_1,n_2,n_3;x,y) (⋯;x,y)  .
                
                We do not explicitly specify the argument becauseν_1andν_2can be at different positions in different terms. The explicit form of the matrix can be found in the  Mathematica notebook file associated with the preprint of the paper on  arXiv. 
                
                The same conclusions as in the previous case apply here as well. The matrixM_321is cosmology independent and hasN^2elements. This makes it numerically much less challenging  thanM_222. In Fig. <ref> we compare our method forN=50sampling points in the power spectrum with the usual numerical evaluation and find an excellent agreement between the two.
                
                
                    < g r a p h i c s >
                
                Equilateral B_411 diagram as a function of k. The calculation is done with bias ν=-0.25 and N=50 sampling points for the power spectrum. This is enough to each the sub-percent precision. B_411 term.—Unlike for previous diagrams, theB_411integral does not converge for any bias. The situation is similar toP_13diagram and the problem can be solved in a similar way. For small negative biases the integral is UV divergent. Therefore, to get the usual numerical result, one has to add the UV part of the loop integral which is given by <cit.>
                    B_411^ UV =    -P_ lin (k_2) P_ lin (k_3) σ_v^2/226380 k_2^2k_3^2( 12409 k_1^6 + 20085 k_1^4(k_2^2 + k_3^2) 
                       - k_1^2 (44518 k_2^4 - 76684 k_2^2k_3^2 +44518 k_3^4) + 12024 (k_2^2 - k_3^2)^2 (k_2^2 + k_3^2) )  + 2   perm.
                
                As in the case of the power spectrum, dimensional regularization would set these terms to zero. As expected, the structure of the UV part of theB_411diagram is such that it can be reabsorbed by the bispectrum counterterms in the EFT approach to LSS <cit.>.
                
                
                The regular terms can be organized in a vectorM_411which hasNelements and is cosmology independent. The approximation to theB_411diagram can be then written as
                B̅_411(k_1,k_2,k_3) = k_1^3 P_ lin(k_1) P_ lin(k_2) ∑_m M_411(ν;x,y) · c_mk_1^-2ν + 2  perms ,
                
                where
                
                M_411(ν;x,y) = 12∑_n_i f_411(n_1,n_2,n_3;x,y) (⋯;x,y)  .
                
                The explicit form of this vector can be found in the  Mathematica notebook file associated with the preprint of the paper on  arXiv. In Fig. <ref> we find an excellent agreement of our method (includingB_411^ UVterms) with the usual numerical evaluation.
                
                
                    < g r a p h i c s >
                
                Equilateral one-loop bispectrum as a function of k. The calculation is done with bias ν=-0.25 and N=50 sampling points for the power spectrum. With these parameters it is possible to reach sub-percent precision on all scales.  The full one-loop bispectrum.—In Fig. <ref> we plot the equilateral bispectrum as a function ofkand compare our method with the usual numerical result. Given that typically there are no large cancellations between different diagrams, we achieve a similar precision for the full result as for each individual term in the sum. 
                
                
                In summary, for each bispectrum shape given by ratiosxandy, one has to calculate three matricesM_222,M_321andM_411. These matrices haveN^3,N^2andNelements respectively whereN∼𝒪(50)is sufficient to achieve sub-percent precision on relevant scales. Computation of these matrices is relatively fast and it can be further optimized using properties of hypergeometric functions. The most practical approach for data analysis is to precompute all matrices and evaluate the one-loop bispectrum with different cosmological parameters as a simple matrix multiplication. 
                
                
                
                
                
                § TWO-LOOP POWER SPECTRUM
                
                Now we can turn to the most complicated case of the two-loop power spectrum. There are four different contributions at this order in perturbation theory <cit.> (for the EFTofLSS treatment of the two-loop power spectrum see <cit.>)
                
                P_ 2-loop(k,τ) = D^4(τ)[P_33^I(k) + P_33^II(k) + P_24(k) + P_15(k) ]  .
                
                The explicit form of the four terms in the square brackets is
                
                P_33^I (k) = 9P_ lin(k) ∫_ F_3(,̨,-) P_ lin(q) ∫_ F_3(-,̨, -) P_ lin(p)  ,
                
                P_33^II (k) = 6 ∫_∫_ F_3(,,-̨-) F_3(-,-,+-)̨ P_ lin(q) P_ lin(p) P_ lin(|-̨-|)  ,
                
                P_24 (k) = 24 ∫_∫_ F_2(, -̨) F_4(,-,-,-)̨ P_ lin(q) P_ lin(p) P_ lin(|-̨|)  ,
                
                P_15 (k) = 30 P_ lin(k) ∫_∫_ F_5(,̨,-,,-) P_ lin(q) P_ lin(p)  .
                
                The corresponding diagrams are shown in Fig. <ref>. In the first contribution,P_33^I (k), two integrals have the same structure as theP_13(k)part of the one-loop calculation. In other cases the integrals are not separable. 
                
                After expanding kernels and linear power spectra in power laws, all terms in the sum can be written in the following form
                
                    ∫_1/q^2ν_4 |-̨|^2ν_5∫_1/p^2ν_1  |-̨|^2ν_2 | - |^2ν_3≡ k^6-2ν_12345 (ν_1,…,ν_5)  .
                
                One important point to make is that at most three of five parametersν_1, …,ν_5are generic complex numbers. The reason for this is that there are at most three linear power spectra in the two-loop integrals. The other two parameters must be integers coming from the expansion of kernels. As we are going to see, this simplifies evaluation of some diagrams significantly. 
                
                
                
                    < g r a p h i c s >
                
                Diagrammatic representation of four contributions to the two-loop power spectrum.
                
                
                Before we move on, let us notice that the second integral in (<ref>) has identical structure as the one-loop bispectrum. Therefore, choosing the following change of coordinatesx=|-̨|^2/k^2andy=q^2/k^2, the function(ν_1,…,ν_5)that we are ultimately interested in can be written as follows 
                (ν_1,…,ν_5) = 1/16π^2∫_x,y x^-ν_5 y^-ν_4(ν_1,ν_2,ν_3;x,y)  .
                
                We are going to use this equation and the series representation of(ν_1,ν_2,ν_3;x,y)to find the explicit formula for(ν_1,…,ν_5). One difficulty is that the region of integration is rather complicated:√(x)+√(y)≥ 1and|√(x)-√(y)|≤ 1. This reflects the constraints that physical momenta in the two-loop diagram have to satisfy. Before we discuss the explicit expression, we derive some symmetry properties of the function(ν_1,…,ν_5). These properties will prove useful in practical applications to the full two-loop integral. 
                
                
                
                
                 §.§ Symmetries of (ν_1,…,ν_5) and Recursion Relations
                
                The two-loop diagram of massless theory is known to have a lot of symmetries which translates to many symmetries of the function(ν_1,…,ν_5). All symmetry transformations were first derived in <cit.>. The full symmetry group isZ_2× S_6and it has2× 6! =1440elements <cit.>. We review here some of the relevant symmetry transformations. To derive them it is enough to use the integral representation
                
                    (ν_1,…,ν_5) = ∫_1/q^2ν_4 |̨̂- |^2ν_5∫_1/p^2ν_1  |̨̂- |^2ν_2 | - |^2ν_3 ,
                
                wherę̂is the unit vector. Two obvious symmetries are the following. First, exchanging integration variablesandleads to(ν_1,ν_2) ↔ (ν_4,ν_5), leavingν_3in the same position
                (ν_1, ν_2, ν_3, ν_4, ν_5) = (ν_4, ν_5, ν_3, ν_1, ν_2)  .
                
                Second, we can simultaneously shift both momenta→-̨and→-̨. Effectively, this produces the simultaneous exchange(ν_1 ↔ν_2)and(ν_4 ↔ν_5)(ν_1, ν_2, ν_3, ν_4, ν_5) = (ν_2, ν_1, ν_3, ν_5, ν_4)  .
                
                Additional translation formulas can be derived using translation formulas of(ν_1,ν_2,ν_3;x,y). Plugging the transformations (<ref>) into the integral representation (<ref>), it is straightforward to get the following extra identities
                (ν_1, ν_2, ν_3, ν_4, ν_5) = (ν_3,ν_1,ν_2, ν_5, ν_6)  ,
                (ν_1, ν_2, ν_3, ν_4, ν_5) = (ν_3,ν_2,ν_1, ν_4, ν_6)  ,
                
                whereν_6≡ 92-ν_12345. Transformations (<ref>) to (<ref>) form the symmetric group of degree fourS_4. This group has4!=24different elements. We recover  these elements by starting from(ν_1, ν_2, ν_3, ν_4, ν_5)and successively applying identities (<ref>) to (<ref>). Invoking notation(1,2,3,4,5)≡(ν_1,…,ν_5), then the 24 equivalent functions are
                
                    (1,2,3,4,5)      (2,1,3,5,4)      (4,5,3,1,2)      (5,4,3,2,1) 
                        (3,1,2,5,6)      (1,3,2,6,5)      (5,6,2,3,1)      (6,5,2,1,3) 
                        (3,2,1,4,6)      (2,3,1,6,4)      (4,6,1,3,2)      (6,4,1,2,3) 
                        (3,5,4,1,6)      (5,3,4,6,1)      (1,6,4,3,5)      (6,1,4,5,3) 
                        (3,4,5,2,6)      (4,3,5,6,2)      (2,6,5,3,4)      (6,2,5,4,3) 
                        (5,2,6,4,1)      (2,5,6,1,4)      (4,1,6,5,2)      (1,4,6,2,5)  .
                
                
                Applying inversion (<ref>) and star-triangle (<ref>) formulae for(ν_1,ν_2,ν_3;x,y)in the integral representation (<ref>), provides analogous results for(ν_1,…,ν_5):
                (ν_1,ν_2,ν_3,ν_4,ν_5) = (ν_2,ν_1,3-ν_123,ν_134- 32 ,ν_235 -32 ) ,
                (ν_1,ν_2,ν_3,ν_4,ν_5) = Γ(ν̃_1)Γ(ν_1)Γ(ν̃_2)Γ(ν_2)Γ(ν̃_3)Γ(ν_3)Γ(3-ν̃_123)Γ(3-ν_123) (ν̃_2,ν̃_1,3-ν̃_123,ν_4,ν_5)  .
                
                Finally, combining these with translation leads to the following reflection formula
                
                    (ν_1,   ν_2,ν_3,ν_4,ν_5) = g (ν̃_1,ν̃_2,ν̃_3,ν̃_4,ν̃_5)  ,
                
                where the prefactorgis given by
                
                g = Γ(ν̃_1)Γ(ν_1)Γ(ν̃_2)Γ(ν_2)Γ(ν̃_3)Γ(ν_3)Γ(ν̃_4)Γ(ν_4)Γ(ν̃_5)Γ(ν_5)Γ(ν̃_6)Γ(ν_6)Γ(3-ν̃_123)Γ(3-ν_123)Γ(3-ν̃_345)Γ(3-ν_345)Γ(ν̃_235-3/2)Γ(ν_235-3/2)Γ(ν̃_134-3/2)Γ(ν_134-3/2) .
                 
                Applying successively all these transformations one can generate the entireZ_2× S_6group. However, almost all of these transformations fail to preserve the bias. As  in the case of the one-loop bipsectrum, only  special choices ofνtrivially offer compact relations between functions with different integer parts of the parameters. 
                
                
                Finally, let us write down an example of a recursion relation that(ν_1,…,ν_5)satisfies. Again, the simplest way to derive these identities is to use the analogous results for(ν_1,ν_2,ν_3;x,y). For instance, using recursion relation (<ref>) one can immediately write
                
                    (3-ν_1-ν_123)    (ν_1,ν_2, ν_3,ν_4,ν_5) 
                    
                    + ν_2    (ν_1,ν_2+1,ν_3,ν_4,ν_5) - ν_2 (ν_1-1,ν_2+1,ν_3,ν_4,ν_5)  
                    
                    +  ν_3    (ν_1,ν_2,ν_3+1,ν_4-1,ν_5) - ν_3 (ν_1-1,ν_2,ν_3+1,ν_4,ν_5) = 0   .
                
                Other similar expressions can be found using the symmetry transformations discussed above. We will give some more details about the practical application of these equations in the following section. 
                 10 pt
                
                
                
                 §.§ Evaluation of (ν_1,…,ν_5)
                
                As we already mentioned, finding an explicit formula for(ν_1,…,ν_5)is not straightforward. In the integral
                (ν_1,…,ν_5) = 1/16π^2∫_x,y x^-ν_5 y^-ν_4(ν_1,ν_2,ν_3,x,y)  ,
                
                the region of integration is complicated and we lack a simple representation of(ν_1,ν_2,ν_3;x,y)convergent in the entire domain. Our strategy to find the full solution is to first concentrate on a part of the region of integration and then use symmetries (<ref>) to find solutions elsewhere. One possibility is to start with the integral
                
                K(ν_1,…,ν_5)  ≡∫_𝒟_1 x^-ν_5 y^-ν_4 J(ν_1,ν_2,ν_3,x,y)  ,
                
                over the following domain𝒟_1 ={ (x,y)  | √(x)+√(y)≥ 1,  x≤1,  y≤ 1 }. Given that bothxandyare smaller than 1, we can use result (<ref>) to evaluate the integral. The details of this derivation are given in Appendix <ref>. Here we report only the final formula. We write the result in a form which resembles the series representation of(ν_1,ν_2,ν_3;x,y)K(ν_1,   ν_2,ν_3,ν_4,ν_5)  = (πν_23)/8√(π)Γ(ν_1) Γ(ν_2) Γ(ν_3) Γ(3-ν_123)
                       ×∑_n=0^∞[ a_n(ν_1,ν_2,ν_3)  κ(  32-ν_235+n ,-ν_4, ν_1+n ,  32-ν_2+n, 3-ν_23+2n ) . 
                           . -  a_n(ν̃_1,ν̃_2,ν̃_3)  κ(-ν_5+n, 32-ν_134, ν̃_1+n,  32 -ν̃_2+n, 3-ν̃_23+2n ) ]  .
                 
                The complicated part of the answer is the functionκ(⋯)which is given in terms of generalized hypergeometric functions
                
                    κ (α,β,a,b,c)     = 11+α[ 11+β _3F_2([  1,a,b;; 2+β,c ; ]1) . - 2·_4F_3([   a,b,1+β,3/2+β;; 1-d ,5/2+γ,3+γ ; ]1)
                    
                         . -2·_4F_3([ c-a,c-b,1+β+d,3/2+β+d;;    1+d,d+γ+5/2,3+γ +d; ]1) ]  ,
                
                whered≡ c-a-bandγ≡α+β. However, this is not the end of the story. We have found only one piece of the final answer, which corresponds to the integral over the region𝒟_1. If one splits the remaining part of domain of integration in (<ref>) in the following way (see Fig. <ref>)
                
                    𝒟_2     ={ (x,y)  | √(x)-√(y)≤ 1,  x≥1,  y≤ x } , 
                    𝒟_3     ={ (x,y)  | √(y)-√(x)≥ 1,  y≥1,  x≤ y } ,
                
                then we find that the integrals over𝒟_2and𝒟_3can be mapped to an integral of the form (<ref>). The proof is based on symmetries of(ν_1,ν_2,ν_3;x,y)given in eq. (<ref>). The full solution is then just a sum of three terms given by (<ref>), with slightly different parameters
                
                    (ν_1,…,ν_5) =    1/16π^2( K(ν_1,ν_2,ν_3,ν_4,ν_5) +  K(ν_1,ν_3,ν_2, ν_6, ν_5) + K(ν_2,ν_3,ν_1, ν_6 ,ν_4) )  .
                
                Obviously, the final result for(ν_1,…,ν_5)is very complicated and not very illuminating. Symmetries of the two-loop diagram that we discussed in the previous section are not manifest at all. It is actually quite remarkable that this messy formula satisfies all functional identities of(ν_1,…,ν_5). This is also an indication that there may exist a much simpler and elegant representation. However, for the time being, it remains elusive. 
                
                
                The difficulties with the result (<ref>) are not only aesthetic but also practical. The biggest problem is that the sum on the r.h.s. of (<ref>) is not convergent for all values of parameters. Even when it is, sometimes there are big cancellations between different terms in the sum. Summing up many large numbers which eventually leads to a small answer can be numerically quite challenging. On the other hand, there is a large region of parameter space where the sum converges very rapidly. Using many of the symmetries of(ν_1,…,ν_5)such as (<ref>), it is always possible to evaluate the sum efficiently for any choice of parameters. Understanding the radius of convergence more quantitatively is very important for knowing ahead of time which symmetry transformation to use. Due to complexity of the final answer making some analytic progress is hard and we leave it for the future work.
                
                
                Another problem with eq. (<ref>) is the appearance of generalized hypergeometric functions whose argument is equal to 1. At this point the hypergeometric series which is usually used to calculate the function is either divergent or converges very slowly for a generic set of complex parameters. In order to calculate the hypergeometric functions efficiently one can use some of the functional identities such as recursion relations. In all examples when the integral in definition of(ν_1,…,ν_5)is convergent, the eq. (<ref>) agrees with the result of numerical integration. Our power series representation is typically several orders of magnitude faster. 
                
                
                So far we considered the general case where the function(ν_1,…,ν_5)depends on five arbitrary complex numbers. However, as we already pointed out, at least two of(ν_1,…,ν_5)are integers which come from the expansion of perturbation theory kernels. In this more specialized case some of the formulas above simplify. There are two different situations that we meet in practice.
                 One of integer parameters is zero or negative.—The simplest case is when one of integer parameters is zero. The integral (<ref>) becomes a product of two one-loop expressions and the result can be written in terms of gamma functions. For example, let us imagine thatν_1=0. It follows
                
                    (0,ν_2,…,ν_5) = k^-6+2ν_2345∫_1/q^2ν_4 |-̨|^2ν_5∫_1/|-̨|^2ν_2 | - |^2ν_3
                          =  k^-6+2ν_2345(ν_2,ν_3) ∫_ |-̨|^3 - 2 ν_23/q^2ν_4 |-̨|^2ν_5 =  (ν_2,ν_3)  (ν_4,ν_235-32)  .
                
                The next simplest case is when none of integer parameters is zero, but rather one of them is negative. Let us imagine thatν_1=-N, whereN>0. In this case the infinite sum in(ν_1,ν_2,ν_3;x,y)truncates (see eq. (<ref>)). The integral (<ref>) can again be expressed in terms of gamma functions only. It is straightforward to get
                
                    (-N,ν_2, ν_3,ν_4,ν_5) = (-1)^N+1√(π)(πν_23)(πν_3)/8 Γ(ν_2) Γ(ν_3) Γ(3+N-ν_23)∑_n=0^N ∑_m=0^N-n
                           N! (-1)^m+n/(N-m-n)! n! m!Γ( 32-ν_2+n+m)/Γ(5/2-ν_23+ n)Γ(ν_3 - N - 12+m) (ν_4-m, ν_235-32-n)  .
                
                Notice that forN=0this expression reduces to (<ref>). This is in agreement with results of <cit.> (see Appendix F of <cit.>). In practice, the sums always have at most a few terms. For all diagrams in the two-loop power spectrumN≤ 5. The cases in which other parameters are non-positive integers can be easily evaluated using translation formulas (<ref>). 
                
                
                Let us point out that in the expansion of the perturbation theory kernels in theP_33,P_24andP_15diagrams, most of the terms do have a negative(zero) integer parameter. For example, the expansion ofF_5kernel in theP_15contribution has several thousand terms. Only ten of them have two positive integer parameters. In other words, the largest part of the two-loop result can be written in terms of gamma functions. Given that in dimensional regularization we do not expect different terms to have very different magnitudes, even neglecting contributions with two positive integer parameters may not affect the result considerably. 
                 Both integer parameters are positive.—Finally, let us discuss the option in which both integer parameters are positive. For the two-loop integral the only possibility is that both of these parameters are equal to one. This comes from the fact that we can have multiple inverse Laplacians in the perturbation theory kernels, but we never have a square of the inverse Laplacian. Under translation formulas (<ref>),  all the apparently different cases reduce to the following two cases
                (1,ν_2,ν_3,ν_4,1)      and    (1,ν_2,ν_3,1,ν_5)  .
                
                Both of these terms can be calculated using eq. (<ref>). It is worth noting that improved numerical stability may be found by considering the related functions provided by the reflection formula (<ref>)(12,ν̃_2,ν̃_3,ν̃_4,12)      and    (12,ν̃_2,ν̃_3,12,ν̃_5)  .
                
                
                As in the case of the one-loop power spectrum and the one-loop bispectrum, from the expansion ofF_nkernels we get a lot of terms whereν_iparameters differ just by an integer. Many of those are related by recursion relations. Let us see how the recursion relations look like in the special case when two parameters are equal to one. For example, ifν_1=ν_4=1, then eq. (<ref>) becomes
                
                    (1-ν_23)    (1,ν_2, ν_3,1,ν_5)
                    + ν_2 (1,ν_2+1,ν_3,1,ν_5) - ν_2 (0,ν_2+1,ν_3,1,ν_5) 
                    
                    +  ν_3    (1,ν_2,ν_3+1,0,ν_5) - ν_3 (0,ν_2,ν_3+1,1,ν_5) = 0   .
                
                Notice that in three of the five terms one of the arguments is equal to zero. Therefore, they can be written in terms of gamma functions. In this way we get a simple functional identity which relates(1,ν_2,ν_3,1,ν_5)and(1,ν_2+1,ν_3,1,ν_5). Whenν_1=ν_5=1we get
                
                    (1-ν_23)    (1,ν_2, ν_3,ν_4,1)
                    + ν_2 (1,ν_2+1,ν_3,ν_4,1) - ν_2 (0,ν_2+1,ν_3,ν_4,1) 
                    
                    +  ν_3    (1,ν_2,ν_3+1,ν_4-1,1) - ν_3 (0,ν_2,ν_3+1,ν_4,1) = 0   .
                
                This equation is slightly more complicated because only two terms have one zero parameter, but it is still very useful. Similar recursion relations can be found exploiting symmetry properties of(ν_1,…,ν_5). In practice, these relations can reduce the number of terms that one has to evaluate by a factor of a few. 
                
                
                The bottom line is that by using explicit expressions for(ν_1,…,ν_5), and its symmetry properties, it is possible to calculate all contributions to the two-loop power spectrum. As before, all information can be compressed in three cosmology independent matricesM_33,M_24andM_15which correspond toP^II_33,P_24andP_15diagrams. Evaluation of these matrices is not trivial because of the convergence properties of the series (<ref>). However, these matrices have to be calculated only once and once they are known, the evaluation of the two-loop power spectrum for any cosmology is just a simple matrix multiplication. These matrices have at mostN^3elements whereN∼𝒪(100). Therefore, evaluation of the two-loop power spectrum in onekbin is significantly faster than the usual numerical techniques. We leave the implementation and testing of our algorithm for the two-loop power spectrum for future work. 
                
                
                
                
                § CONCLUSIONS
                
                
                In this paper we demonstrate the path forward for the efficient computation of higher multiplicity/loop correlation functions in cosmological perturbation theory. Our starting point is similar to recent proposals for fast evaluation of the one-loop power spectrum <cit.> and it is based on representing the linear power spectrum as a sum of complex power laws. However, our implementation and generalization to higher order correlators is different. When comparison is possible, all methods agree. 
                
                
                We mainly focus on deriving relevant analytic expressions for the one-loop and the two-loop power spectrum and the one-loop bispectrum. All one-loop diagrams evaluated using our method are in excellent agreement with the usual numerical results. We leave writing a dedicated code for the two-loop power spectrum (and possibly higher order correlation functions) for future work. 
                
                
                Our method splits the computation of loop diagrams in two parts. The first, more “difficult” part is related to solving momentum integrals for power-law power spectra and it is cosmology independent. The second part is a simple matrix multiplication which evaluates the loops for aΛCDM-like cosmology. The matrices can be precomputed, they are cosmology independent and they are relatively small. For example, for the two-loop power spectrum, the largest matrix hasN^3elements, whereN∼𝒪(100). The number of operations needed for evaluation of the power spectrum or the bispectrum is significantly smaller than using direct numerical integration. Furthermore, the same building blocks used to calculate dark matter correlation functions can be also used for correlators of biased tracers. There are no fundamental obstacles in applying our method in redshift space as well. 
                
                
                One interesting aspect of the method described in this paper is that it relies on evaluation of loop integrals that are formally identical to those of a massless QFT. This is a new bridge between cosmology and particle physics and the full potential of this connection is still to be explored. This remains the major direction for future work. One hope is that many developments in the theory of scattering amplitudes will prove useful for going beyond the lowest order statistics discussed in this paper. The first step in this direction is a practical analytic formula for the one-loop trispectrum. In principle, that would allow the calculation of the three-loop power spectrum or the two-loop bispectrum. In practice, following procedure described in this paper may turn out to be too difficult or impractical. After all, the trispectrum is a function of six variables, which makes it much more complicated than examples we considered so far. 
                
                
                However, there are many alternative representations of loop integrals that may be more useful for higher loop diagrams. In this paper we have insisted on finding well-behaved power series representations for functions such as(ν_1,ν_2)or(ν_1,ν_2,ν_3;x,y). Alternative ways to evaluate these integrals include numerical integration using Mellin-Barnes representation of loop integrals (see for instance <cit.>), projecting onto a basis of other known higher-loop integrals as in <cit.>, or solving numerically partial differential equations that the loop integrals satisfy <cit.>.
                
                
                In some important situations things simplify. One such example is the one-loop covariance of the power spectrum. Given that there are only two independent vectors_̨1and_̨2, this special case of the one-loop four-point function depends only on three variables (before integrating over the angle between_̨1and_̨2). Furthermore, a set of diagrams in the one-loop covariance matrix which give the largest contribution to the final answer (see <cit.>) have the same structure as the one-loop bispectrum. These diagrams can be easily calculated using our function(ν_1,ν_2,ν_3;x,y). We leave application of our method to the covariance matrix and more generally one-loop four point function for future work.
                
                
                At the end, let us stress that the idea of representing theΛCDM-like cosmology as a set of power-law universes can be also very useful outside the context of PT. One can benefit form this decomposition whenever some numerically heavy integral has a simple solution for a power-law universe. One example of this kind is projection of the power spectrum or the bispectrum on the sky, which is difficult due to many integrals over highly oscillatory spherical Bessel functions. It was shown in <cit.> that decomposition (<ref>) can be used to find the solutions of these integrals very accurately and efficiently. It would be interesting to think of other similar applications in the future. 
                 10 pt
                
                
                
                § ACKNOWLEDGMENTS
                
                We would like to thank Nima Arkani-Hamed, Valentin Assassi, Diego Blas, Jonathan Blazek, Paolo Creminelli, Guido D'Amico, Chris Hirata, Lam Hui, Mikhail Ivanov, David Kosower, Marcel Schmittfull, Roman Scoccimarro, Leonardo Senatore, Sergey Sibiryakov, Kris Sigurdson, Zachary Slepian, Gabriele Trevisan and Zvonimir Vlah for many useful discussions.  J. J. M. C. is supported by the
                European Research Council under ERC-STG-639729,  preQFT: Strategic
                  Predictions for Quantum Field Theories. M.S. gratefully acknowledges support from the Institute for Advanced Study and the Raymond and Beverly Sackler Foundation. M.Z. is supported by NSF grants AST-1409709 and PHY-1521097 and by the Canadian Institute for Advanced Research (CIFAR) program on Gravity and the Extreme Universe.
                
                
                
                § HYPERGEOMETRIC FUNCTIONS
                
                The hypergeometric function _2F_1(a,b,c,z)is usually defined as a solution of Euler's hypergeometric equation:
                
                z(1-z) f”(z)+(c-(a+b+1)z) f'(z) -ab f(z) = 0 ,
                
                wherea,bandcare arbitrary complex numbers. The hypergeometric function has the power series representation:
                 _2F_1(a,b,c,z)  =  Γ(c)/Γ(a)Γ(b)∑_n=0^∞Γ(a+n)Γ(b+n)/Γ(c+n)n! z^n  ,
                
                which is convergent inside the unit circle in the complex plane|z|<1. This power series can be used for numerical evaluation. The series is convergent at the pointz=1only when the parameters satisfy Re(c-a-b)>0. It should be stressed that the convergence sometimes may be slow or the series has large cancellations, particularly for parameters with large imaginary parts. In order to avoid such issues or evaluate the hypergeometric function outside the unit disc, one can use many functional identities. For example, one such identity is 
                
                    _2F_1(a,b,c,1-z)     =   Γ(c) Γ(c - a - b)/Γ(c - a) Γ(c - b) _2F_1(a, b, a + b - c + 1, z) 
                            + Γ(c) Γ(a + b - c)/Γ(a) Γ(b) z^c - a - b _2F_1(c - a, c - b, 1 - a - b + c, z)  ,
                
                which maps points close to|z|=1to a region aroundz=0where the series converges rapidly. Outside the unit disc the hypergeometric function can be calculated using 
                _2F_1(a,b,c,1/z) = Γ(b - a) Γ(c)/Γ(b) Γ(c - a) (-z)^a _2F_1(a, a - c + 1, a - b + 1, z)   +   (a↔ b) .
                
                These two identities are sufficient to evaluate the hypergeometric function in the entire complex plane.
                
                
                It is possible to generalize the basic hypergeometric series (<ref>) and use it to define generalized hypergeometric functions
                _pF_q([ a_1 , a_2 ,…, a_p;  b_1 , b_2 ,…,b_q ];  z )  ≡ Γ(b_1)⋯Γ(b_q)/Γ(a_1)⋯Γ(a_p)∑_n=0^∞Γ(a_1+n)⋯Γ(a_p+n)/Γ(b_1+n)⋯Γ(b_q+n)z^n/n! ,
                
                wherepandqare positive integers. In this paper we use two generalized hypergeometric functions_4F_3and_3F_2. In these cases whenp=q+1the generalized hypergeometric series (<ref>) converges for|z|<1. A the point isz=1which is of special interest in evaluation of the two-loop power spectrum the series converges only when Re(b_1+⋯+b_q-a_1-⋯ a_q+1)>0. This condition is not always satisfied in practice. One simple way out is to use recursion relations which increase the real part of one ofb_iuntil the series becomes convergent.
                
                
                § EXPLICIT FORM OF THE M_222 MATRIX
                
                In this appendix we give the explicit form of theM_222matrix. The starting point are threeF_2kernels in theB_222diagram
                
                8 F_2 (,_̨1-) F_2(_̨1-,_̨2+) F_2(_̨2+,-)  .
                
                Expanding this expression in powers ofq^2,|_̨1-|^2and|_̨2+|^2for a single set of parameters(ν_1,ν_2,ν_3)and factoring outk_1dependence we get a sum which can be rearranged in the following way
                
                    M_222 = 25x^21372_-2,2,2 -15 (y+1) x^21372_-1,2,2 -(10 y^2-9 y+10)x^21372_0,2,2 +3x^2y (y+1)686_1,2,2 +x^2y^2343_2,2,2 +75x2744_-2,1,2
                        +75x2744_-2,2,1-5 x(20 x+9 y+9)2744_-1,1,2 -5x(20x+9y+9)2744_-1,2,1 +3x(-10 y^2+9 y+10 x (2 y-1)-10)2744_0,1,2
                         -3x(10 y^2-9 y+10 x (y-2)+10)2744_0,2,1 +x y (20 y x+9 x+9 y+9)1372_1,1,2 +x(9 y (y+1)+x (9 y+20))1372_1,2,1 + 3x(x+1) y^2686_2,1,2
                        +3xy(x+y)686_2,2,1 -1252744_-2,0,2 + 1251372_-2,1,1 - 1252744_-2,2,0 + 1251372_-1,-1,2 - 1251372_-1,0,1 -75 (2 x-y-1)2744_-1,0,2
                        -1251372_-1,1,0 -75 (2 x+y+1)1372_-1,1,1 + 1251372_-1,2,-1 - 75 (2 x-y-1)2744_-1,2,0 - 1252744_0,-2,2 - 1251372_0,-1,1
                        + 75(x-2y+1)2744_0,-1,2 + 3751372_0,-1,2 + 75(5x+5y-4)2744_0,0,1 + 5 (10 x^2+9 (2 y-1) x+10 y^2-9 y+10)2744_0,0,2 - 1251372_0,1,-1
                        + 75 (5 x-4 y+5)2744_0,1,0 + 5 (40 x^2+9 (y+1) x-20 y^2+18y-20)2744_0,1,1 -1252744_0,2,-2 + 75 (x+y-2)2744_0,2,-1
                        + 5(10 x^2-9 (y-2) x+10 y^2-9 y+10)2744_0,2,0 + 1251372_1,-2,1 + 75 y2744_1,-2,2 - 1251372_1,-1,0 - 75 (x+2 y+1)1372_1,-1,1
                        - 5y(9x+20y+9)2744_1,-1,2 -1251372_1,0,-1 - 75(4 x-5 (y+1))2744_1,0,0 - 5 (20x^2-9 (y+2) x-40 y^2-9 y+20)2744_1,0,1
                        + 3 y (-10 x^2+(20 y+9) x-10 (y+1))2744_1,0,2 + 1251372_1,1,-2 - 75(x+y+2)1372_1,1,-1 - 5(20 x^2-9 (2
                       y+1) x+20 y^2-9 y-40)2744_1,1,0
                        + 3(10 (y+1) x^2+(10 y^2+9 y+10)x+10y (y+1))1372_1,1,1 + 752744_1,2,-2 - 5 (9 x+9 y+20)2744_1,2,-1
                        - 3(10 x^2-(9 y+20)x+10y(y+1))2744_1,2,0 - 1252744_2,-2,0 + 75y2744_2,-2,1 + 25 y^21372_2,-2,2 + 1251372_2,-1,-1
                        + 75 (x-2 y+1)2744_2,-1,0 - 5 y (9 x+20 y+9)2744_2,-1,1 - 15(x+1)y^21372_2,-1,2 - 1252744_2,0,-2 + 75 (x+y-2)2744_2,0,-1
                        + 5 (10 x^2-9 (y+1) x+2 (5 y^2+9y +5))2744_2,0,0 - 3 y (10 x^2+(10 y-9) x-20 y+10)2744_2,0,1 -(10 x^2-9 x+10) y^21372_2,0,2
                        + 752744_2,1,-2 - 5(9x+9y+20)2744_2,1,-1 - 3(10 x^2+(10-9 y) x+10 (y-2) y)2744_2,1,0 + y(9 x^2+9 (y+1) x+20 y)1372_2,1,1
                        + 251372_2,2,-2 - 15 (x+y)1372_2,2,-1 - (10 x^2-9 y x+10 y^2)1372_2,2,0 .
                
                We are using shorten notation in which_n_1,n_2,n_3≡(ν_1+n_1,ν_2+n_2,ν_3+n_3;x,y). Coefficientsf_222can be easily read off form this expression. Notice that there are 72 terms in the sum, but not all of them are independent. Using recursion relations (<ref>) one can further reduce this expression to a sum of 38 different_n_1,n_2,n_3functions. We do not write this sum explicitly because the coefficients multiplying_n_1,n_2,n_3functions become too cumbersome. Nevertheless, these new coefficients are still only rational functions that depend onx,y,ν_1,ν_2andν_3and application of recursion relations effectively reduces the cost of evaluating theB_222diagram by roughly a factor of 2. 
                
                
                
                § DERIVATION OF (Ν_1,Ν_2, Ν_3 ;X,Y)
                 
                Let us begin with the usual Feynman parametrization
                
                    ∫_1/q^2ν_1|_̨1-|^2ν_2 |_̨2+|^2ν_3 =Γ(ν_123)/Γ(ν_1) Γ(ν_2) Γ(ν_3)
                         ×∫_0^1 du_1 ∫_0^1 du_2 ∫_0^1 du_3 ∫_ u_1^ν_1-1u_2^ν_2-1u_3^ν_3-1δ^(D)(1-u_1-u_2-u_3)/(u_1q^2 +u_2|_̨1-|^2 + u_3 |_̨2+|^2 )^ν_123
                        = Γ(ν_123)/Γ(ν_1) Γ(ν_2) Γ(ν_3)∫_0^1 du_1 ∫_0^1-u_1 du_2 ∫_ u_1^ν_1-1u_2^ν_2-1(1-u_1-u_2)^ν_3-1/(u_1q^2 +u_2|_̨1-|^2 + (1-u_1-u_2) |_̨2+|^2 )^ν_123 .
                
                Next, we do the following change of variables:u_1=uvandu_2=(1-u)v. This transforms(1-u_1-u_2)intovand now both integrals inuandvhave the same boundaries[0,1]∫_1/q^2ν_1|_̨1-|^2ν_2 |_̨2+|^2ν_3 = Γ(ν_123)/Γ(ν_1)Γ(ν_2)Γ(ν_3)
                               ×∫_0^1 du ∫_0^1 dv ∫_ u^ν_1-1(1-u)^ν_2-1v^ν_12-1(1-v)^ν_3-1/( uvq^2 +(1-u)v|_̨1-|^2 +(1-v)|_̨2+|^2 )^ν_123 .
                
                At this point the momentum integral can be done easily. In the denominator we first complete the square
                
                    uvq^2 +(1-u)v|_̨1-|^2 +(1-v)|_̨2+|^2 
                          = (  - (1-u) v_̨1 + (1-v)_̨2 )^2 + v ( u v (1-u) k_1^2 + u(1-v)k_2^2 + (1-u) (1-v) k_3^2  )  ,
                
                and use the following identity to do the integral in∫_1/(q^2 +m^2)^ν_123 = 1/8π^3/2Γ(ν_123-3/2)/Γ(ν_123)1/(m^2)^ν_123-3/2 .
                
                The expression for one-loop bispectrum simplifies and we are left with two integrals inuandv∫_1/q^2ν_1|_̨1-|^2ν_2 |_̨2+|^2ν_3 = k_1^3-2ν_123/8π^3/2Γ( ν_123-3/2)/Γ(ν_1)Γ(ν_2)Γ(ν_3)
                               ×∫_0^1 du ∫_0^1 dv u^ν_1-1(1-u)^ν_2-1v^1/2-ν_3(1-v)^ν_3-1/( u v (1-u) + u(1-v) y + (1-u) (1-v) x )^ν_123-3/2 ,
                
                from which we can read off(ν_1,ν_2, ν_3 ;x,y)(ν_1,ν_2, ν_3 ;x,y)     = 1/8π^3/2Γ(ν_123-3/2)/Γ(ν_1)Γ(ν_2)Γ(ν_3)
                         ×∫_0^1 du ∫_0^1 dv u^ν_1-1(1-u)^ν_2-1v^1/2-ν_3(1-v)^ν_3-1/( u v (1-u) + u(1-v) y + (1-u) (1-v) x )^ν_123-3/2 .
                
                Notice that the denominator is linear invand that the integral invis nothing but the hypergeometric function
                
                    (ν_1,ν_2,ν_3 ;x,y) = Γ(3/2-ν_3) Γ( ν_123-3/2)/4π^2Γ(ν_1) Γ(ν_2)∫_0^1 du u^ν_1-1(1-u)^ν_2-1
                         × (x(1-u)+y u)^3/2-ν_123 _2F_1(3/2 - ν_3, ν_123-3/2, 3/2, 1- u(1-u)/x(1-u)+y u)  .
                
                At this point it is useful to transform this expression using (<ref>)
                    (ν_1,ν_2,ν_3;x,y) = Γ(3/2-ν_3) Γ( ν_123-3/2)/4π^2Γ(ν_1) Γ(ν_2)∫_0^1 du  u^ν_1-1(1-u)^ν_2-1(x(1-u)+y u)^3/2-ν_123
                         [ Γ( 3/2) Γ(ν_12-3/2)/Γ( 3/2-ν_3 ) Γ( ν_123-3/2)u^3/2-ν_12(1-u)^3/2-ν_12/(x(1-u)+y u)^3/2-ν_12 _2F_1(ν_3, 3-ν_123, 5/2-ν_12, u(1-u)/x(1-u)+y u) . 
                         . + Γ( 3/2) Γ( 3/2- ν_12)/Γ(ν_3)Γ(3-ν_123) _2F_1(3/2 - ν_3, ν_123-3/2, ν_12- 1/2, u(1-u)/x(1-u)+y u)]  .
                
                The reason is that0≤ u(1-u)/(x(1-u)+y u)≤1for anyxandy, and one can use the power series representation of hypergeometric functions in order to solve the integral inu.
                Notice that this power series keeps the integral as simple as  possible, because only powers oru,(1-u)and(x(1-u)+y u)appear in the expression. Simplifying the gamma functions we get
                
                    (ν_1,ν_2,   ν_3;x,y) = (πν_12)/8√(π)Γ(ν_1) Γ(ν_2) Γ(ν_3)Γ(3-ν_123) 
                         [ ∑_n=0^∞Γ(ν_3+n) Γ(3-ν_123+n)/Γ(5/2-ν_12+n )n!∫_0^1 du u^1/2-ν_2+n(1-u)^1/2-ν_1+n/(x(1-u)+y u)^ν_3+n. 
                         . - ∑_n=0^∞Γ(3/2 - ν_3+n ) Γ( ν_123-3/2+n)/Γ( ν_12- 1/2+n )n!∫_0^1 du u^ν_1-1+n(1-u)^ν_2-1+n/(x(1-u)+y u)^ν_123-3/2+n]  .
                
                The integration inuleads to another hypergeometric function. The result can be written in the following way
                
                    (ν_1,    ν_2,ν_3;x,y) = (πν_12)/8√(π)Γ(ν_1) Γ(ν_2) Γ(ν_3)Γ(3-ν_123) 
                        [ ∑_n=0^∞ b_n(ν_1,ν_2,ν_3)· x^-ν_3-n _2F_1(3/2 - ν_2+n, ν_3+n, 3-ν_12+2n, 1-y/x) . 
                         . - ∑_n=0^∞ b_n(ν̃_1,ν̃_2,ν̃_3)·  x^-ν_2-n y^3/2-ν_13 _2F_1(3/2-ν_3+n, ν_2+n, ν_12+2n, 1-y/x)  ]  ,
                
                where the coefficientsb_nare given by
                
                b_n(ν_1,ν_2,ν_3) = Γ(ν_3+n) Γ(3-ν_123+n)/Γ(5/2-ν_12+n )n!Γ(3/2-ν_1+n ) Γ(3/2-ν_2+n )/Γ(3-ν_12+2n ) .
                
                One last step is to use the identity 
                (ν_1,ν_2,ν_3;x,y)=x^3/2-ν_123(ν_3,ν_2,ν_1;1/x, y/x) ,
                
                in order to bring the result to its final form
                
                    (ν_1,ν_2,    ν_3;x,y) = (πν_23)/8√(π)Γ(ν_1) Γ(ν_2) Γ(ν_3) Γ(3-ν_123)
                       ∑_n=0^∞[ x^3/2-ν_23· b_n(ν_3,ν_2,ν_1)   x^n _2F_1(ν_1+n, 3/2-ν_2+n, 3-ν_23+2n, 1-y ) . 
                         . -y^3/2-ν_13· b_n(ν̃_3,ν̃_2,ν̃_1)   x^n _2F_1(  ν_2+n, 3/2 -ν_1+n, ν_23+2n, 1-y  ) ]  .
                
                This precisely matches eq. (<ref>) wherea_n(ν_1,ν_2,ν_3) = b_n(ν_3,ν_2,ν_1). 
                
                
                
                
                § DERIVATION OF (Ν_1,…,Ν_5)
                
                We are interested in calculating
                
                K(ν_1,…, ν_5) ≡∫_𝒟_1 x^-ν_5 y^-ν_4 J(ν_1,ν_2,ν_3,x,y)  , 
                
                where the region of integration is given by
                𝒟_1 ={ (x,y)  | √(x)+√(y)≥ 1,  x≤1,  y≤ 1 } .
                
                In this domain the power series representation of(ν_1,ν_2,ν_3;x,y)is uniformly convergent, and we can use it to rewrite the integral in the following way
                
                    K(ν_1,   …, ν_5) = (πν_23)/8√(π)Γ(ν_1) Γ(ν_2) Γ(ν_3) Γ(3-ν_123)∫_0^1dy∫_(1-√(y))^2^1 d x  
                       [ x^3/2-ν_23∑_n=0^∞  a_n(ν_1,ν_2,ν_3)·  x^n _2F_1(ν_1+n,  32-ν_2+n, 3-ν_23+2n, 1-y ) . 
                       . -y^3/2-ν_13∑_n=0^∞ a_n(ν̃_1,ν̃_2,ν̃_3)·  x^n _2F_1(ν̃_1+n,  32-ν̃_2+n, 3-ν̃_23+2n, 1-y ) ]  .
                
                Therefore, the basic integral that we want to solve is the integral over the hypergeometric function. For simplicity, let us define 
                
                    κ(α,β,a,b,c) ≡∫_0^1 dy ∫_(1-√(y))^2^1 d x  x^α y^β _2F_1(a,b,c,1-y)  .
                
                The integral inxis straightforward, leading to
                
                    κ(α,β,a,b,c)     = 1/1+α[∫_0^1 dy   y^β _2F_1(a,b,c,1-y) . 
                         . - 2 ∫_0^1 dt   t^2β+1 (1-t)^2+2α _2F_1(a,b,c,1-t^2) ]  ,
                
                where in the second integral we did a change of variablesy=t^2. Both integrals can be expressed in terms of higher order hypergeometric functions. It is not difficult to find
                
                    κ (α,β,a,b,c)     = 1/1+α[ 1/1+β _3F_2([  1,a,b;; 2+β,c ; ]1) . - 2·_4F_3([   a,b,1+β,3/2+β;; 1-d ,5/2+γ,3+γ ; ]1)
                    
                       . -2·_4F_3([ c-a,c-b,1+β+d,3/2+β+d;;    1+d,d+γ+5/2,3+γ +d; ]1) ]  ,
                
                whered=c-a-bandγ=α+β. The integration of_2F_1(⋯,1-y)function is straightforward. To integrate_2F_1(⋯,1-t^2)we first have to use (<ref>), expand the hypergeometric functions in power series, integrate in series, then resum the result. 
                
                
                99
                Lifshitz:1946a
                  E. M. Lifshitz,
                  “On the gravitational stability of the expanding universe,"
                  Zhurnal Eksperimentalnoi i Teoreticheskoi Fiziki   587, 16 (1946)
                
                
                Fry:1983cj 
                  J. N. Fry,
                  “The Galaxy correlation hierarchy in perturbation theory,”
                  Astrophys. J.  279, 499 (1984).
                  
                  
                  
                Goroff:1986ep 
                  M. H. Goroff, B. Grinstein, S. J. Rey and M. B. Wise,
                  “Coupling of Modes of Cosmological Mass Density Fluctuations,”
                  Astrophys. J.  311, 6 (1986).
                  
                  
                  
                Scoccimarro:1995if 
                  R. Scoccimarro and J. Frieman,
                  “Loop corrections in nonlinear cosmological perturbation theory,”
                  Astrophys. J. Suppl.  105, 37 (1996)
                  
                  [astro-ph/9509047].
                  
                  
                Baumann:2010tm 
                  D. Baumann, A. Nicolis, L. Senatore and M. Zaldarriaga,
                  “Cosmological Non-Linearities as an Effective Fluid,”
                  JCAP  1207, 051 (2012)
                  
                  [arXiv:1004.2488 [astro-ph.CO]].
                  
                  
                Carrasco:2012cv 
                  J. J. M. Carrasco, M. P. Hertzberg and L. Senatore,
                  “The Effective Field Theory of Cosmological Large Scale Structures,”
                  JHEP  1209, 082 (2012)
                  
                  [arXiv:1206.2926 [astro-ph.CO]].
                  
                  
                Lewandowski:2014rca 
                  M. Lewandowski, A. Perko and L. Senatore,
                  “Analytic Prediction of Baryonic Effects from the EFT of Large Scale Structures,”
                  JCAP  1505, no. 05, 019 (2015)
                  
                  [arXiv:1412.5049 [astro-ph.CO]].
                  
                  
                Bernardeau:2001qr 
                  F. Bernardeau, S. Colombi, E. Gaztanaga and R. Scoccimarro,
                  “Large scale structure of the universe and cosmological perturbation theory,”
                  Phys. Rept.  367, 1 (2002)
                  
                  [astro-ph/0112551].
                  
                  
                McEwen:2016fjn 
                  J. E. McEwen, X. Fang, C. M. Hirata and J. A. Blazek,
                  “FAST-PT: a novel algorithm to calculate convolution integrals in cosmological perturbation theory,”
                  JCAP  1609, no. 09, 015 (2016)
                  
                  [arXiv:1603.04826 [astro-ph.CO]].
                  
                  
                  
                
                Schmittfull:2016jsw 
                  M. Schmittfull, Z. Vlah and P. McDonald,
                  “Fast large scale structure perturbation theory using one-dimensional fast Fourier transforms,”
                  Phys. Rev. D  93, no. 10, 103528 (2016)
                  
                  [arXiv:1603.04405 [astro-ph.CO]].
                  
                  
                  
                
                Schmittfull:2016yqx 
                  M. Schmittfull and Z. Vlah,
                  “FFT-PT: Reducing the two-loop large-scale structure power spectrum to low-dimensional radial integrals,”
                  Phys. Rev. D  94, no. 10, 103530 (2016)
                  
                  [arXiv:1609.00349 [astro-ph.CO]].
                  
                  
                  
                
                Hamilton:1999uv 
                  A. J. S. Hamilton,
                  “Uncorrelated modes of the nonlinear power spectrum,”
                  Mon. Not. Roy. Astron. Soc.  312, 257 (2000)
                  
                  [astro-ph/9905191].
                  
                  
                Scoccimarro:1996jy 
                  R. Scoccimarro,
                  “Cosmological perturbations: Entering the nonlinear regime,”
                  Astrophys. J.  487, 1 (1997)
                  
                  [astro-ph/9612207].
                  
                  
                  
                
                Pajer:2013jj 
                  E. Pajer and M. Zaldarriaga,
                  “On the Renormalization of the Effective Field Theory of Large Scale Structures,”
                  JCAP  1308, 037 (2013)
                  
                  [arXiv:1301.7182 [astro-ph.CO]].
                  
                  
                Gorishnii:1984te 
                  S. G. Gorishnii and A. P. Isaev,
                  “On an Approach to the Calculation of Multiloop Massless Feynman Integrals,”
                  Theor. Math. Phys.  62, 232 (1985)
                  [Teor. Mat. Fiz.  62, 345 (1985)].
                  
                  
                  
                McDonald:2009dh 
                  P. McDonald and A. Roy,
                  “Clustering of dark matter tracers: generalizing bias for the coming era of precision LSS,”
                  JCAP  0908, 020 (2009)
                  
                  [arXiv:0902.0991 [astro-ph.CO]].
                  
                  
                Senatore:2014eva 
                  L. Senatore,
                  “Bias in the Effective Field Theory of Large Scale Structures,”
                  JCAP  1511, no. 11, 007 (2015)
                  
                  [arXiv:1406.7843 [astro-ph.CO]].
                  
                  
                Assassi:2014fva 
                  V. Assassi, D. Baumann, D. Green and M. Zaldarriaga,
                  “Renormalized Halo Bias,”
                  JCAP  1408, 056 (2014)
                  
                  [arXiv:1402.5916 [astro-ph.CO]].
                  
                  
                  
                
                Mirbabayi:2014zca 
                  M. Mirbabayi, F. Schmidt and M. Zaldarriaga,
                  “Biased Tracers and Time Evolution,”
                  JCAP  1507, no. 07, 030 (2015)
                  
                  [arXiv:1412.5169 [astro-ph.CO]].
                  
                  
                  
                
                Desjacques:2016bnm 
                  V. Desjacques, D. Jeong and F. Schmidt,
                  “Large-Scale Galaxy Bias,”
                  arXiv:1611.09787 [astro-ph.CO].
                  
                  
                Baldauf:2014qfa 
                  T. Baldauf, L. Mercolli, M. Mirbabayi and E. Pajer,
                  “The Bispectrum in the Effective Field Theory of Large Scale Structure,”
                  JCAP  1505, no. 05, 007 (2015)
                  
                  [arXiv:1406.4135 [astro-ph.CO]].
                  
                  
                  
                
                Angulo:2014tfa 
                  R. E. Angulo, S. Foreman, M. Schmittfull and L. Senatore,
                  “The One-Loop Matter Bispectrum in the Effective Field Theory of Large Scale Structures,”
                  JCAP  1510, no. 10, 039 (2015)
                  
                  [arXiv:1406.4143 [astro-ph.CO]].
                  
                  
                Davydychev:1992xr 
                  A. I. Davydychev,
                  “Recursive algorithm of evaluating vertex type Feynman integrals,”
                  J. Phys. A  25, 5587 (1992).
                  
                  
                Exton:1994ab 
                  H. Exton,
                  “On the system of partial differential equations associated with Appell's function F_4,”
                  Journal of Physics A Mathematical General
                 28, 631-641 (1995).
                
                
                Morita:1996ab 
                  T. Morita,
                  “Use of the Gauss Contiguous Relations in Computing the Hypergeometric Functions F(n+1/2,n+1/2,m;z),”
                  Interdisciplinary Information Sciences
                 2, 63-74 (1996).  
                
                
                Carrasco:2013sva 
                  J. J. M. Carrasco, S. Foreman, D. Green and L. Senatore,
                  “The 2-loop matter power spectrum and the IR-safe integrand,”
                  JCAP  1407, 056 (2014)
                  
                  [arXiv:1304.4946 [astro-ph.CO]].
                  
                  
                Carrasco:2013mua 
                  J. J. M. Carrasco, S. Foreman, D. Green and L. Senatore,
                  “The Effective Field Theory of Large Scale Structures at Two Loops,”
                  JCAP  1407, 057 (2014)
                  
                  [arXiv:1310.0464 [astro-ph.CO]].
                  
                  
                Baldauf:2015aha 
                  T. Baldauf, L. Mercolli and M. Zaldarriaga,
                  “Effective field theory of large scale structure at two loops: The apparent scale dependence of the speed of sound,”
                  Phys. Rev. D  92, no. 12, 123007 (2015)
                  
                  [arXiv:1507.02256 [astro-ph.CO]].
                  
                  
                Cataneo:2016suz 
                  M. Cataneo, S. Foreman and L. Senatore,
                  “Efficient exploration of cosmology dependence in the EFT of LSS,”
                  JCAP  1704, no. 04, 026 (2017)
                  
                  [arXiv:1606.03633 [astro-ph.CO]].
                  
                  
                Broadhurst:1986bx 
                  D. J. Broadhurst,
                  “Exploiting the 1.440 Fold Symmetry of the Master Two Loop Diagram,”
                  Z. Phys. C  32, 249 (1986).
                  
                  
                  
                  
                
                Barfoot:1987kg 
                  D. T. Barfoot and D. J. Broadhurst,
                  “Z(2) X S(6) Symmetry of the Two Loop Diagram,”
                  Z. Phys. C  41, 81 (1988).
                  
                  
                  
                Czakon:2005rk 
                  M. Czakon,
                  “Automatized analytic continuation of Mellin-Barnes integrals,”
                  Comput. Phys. Commun.  175, 559 (2006)
                  
                  [hep-ph/0511200].
                  
                  
                Johansson:2012sf 
                  H. Johansson, D. A. Kosower and K. J. Larsen,
                  “An Overview of Maximal Unitarity at Two Loops,”
                  PoS LL  2012, 066 (2012)
                  [arXiv:1212.2132 [hep-th]].
                  
                  
                Henn:2014qga 
                  J. M. Henn,
                  “Lectures on differential equations for Feynman integrals,”
                  J. Phys. A  48, 153001 (2015)
                  
                  [arXiv:1412.2296 [hep-ph]].
                  
                  
                Mohammed:2016sre 
                  I. Mohammed, U. Seljak and Z. Vlah,
                  “Perturbative approach to covariance matrix of the matter power spectrum,”
                  Mon. Not. Roy. Astron. Soc.  466, no. 1, 780 (2017)
                  
                  [arXiv:1607.00043 [astro-ph.CO]].
                  
                  
                Barreira:2017kxd 
                  A. Barreira and F. Schmidt,
                  “Response Approach to the Matter Power Spectrum Covariance,”
                  arXiv:1705.01092 [astro-ph.CO].
                  
                  
                
                Assassi:2017lea 
                  V. Assassi, M. Simonović and M. Zaldarriaga,
                  “Efficient Evaluation of Cosmological Angular Statistics,”
                  arXiv:1705.05022 [astro-ph.CO].]
