Inflation Correlators with Multiple Massive Exchanges

The most general tree-level boundary correlation functions of quantum fields in inflationary spacetime involve multiple exchanges of massive states in the bulk, which are technically difficult to compute due to the multi-layer nested time integrals in the Schwinger-Keldysh formalism. On the other hand, correlators with multiple massive exchanges are well motivated in cosmological collider physics, with the original quasi-single-field inflation model as a notable example. In this work, with the partial Mellin-Barnes representation, we derive a simple rule, called family-tree decomposition, for directly writing down analytical answers for arbitrary nested time integrals in terms of multi-variable hypergeometric series. We present the derivation of this rule together with many explicit examples. This result allows us to obtain analytical expressions for general tree-level inflation correlators with multiple massive exchanges. As an example, we present the full analytical results for a range of tree correlators with two massive exchanges.


Introduction
Recent years have witnessed increasing interests in the theoretical study of cosmological correlation functions of large-scale fluctuations, which are believed to be sourced by quantum fluctuations of spacetime and matter fields during cosmic inflation [1]. By observing the correlation functions of the large-scale structure, we can access quantum field theory in inflationary spacetime. This connection has far-reaching consequences to both early-universe cosmology and fundamental particle physics. It has been emphasized that heavy particles produced during cosmic inflation could leave characteristic and oscillatory signals in certain soft limits of correlation functions. Many recent studies have exploited this Cosmological Collider (CC) signal to study particle physics at the inflation scale . At the same time, there are considerable works devoting to the analytical or numerical study of correlation functions of quantum field theory in inflationary spacetime, or inflation correlators for short . These studies have revealed many interesting structures of inflation correlators or wavefunctions, which deepen our understanding of quantum field theory in de Sitter spacetime. On the other hand, explicit analytical results are indispensable for a precise understanding of CC signals and for comparing theoretical predictions of CC models with observational data.
Many explicit analytical results have been obtained in recent years for inflation correlators relating to CC physics [63,64,70,71,[88][89][90][91][92][93]. Most of these results are for the exchange of a single massive particle in the bulk of dS, with a few exceptions at loop orders. However, previous works on CC model building have shown that correlators with multiple exchanges of massive particles could be phenomenologically important. Already in the early studies of quasi-single-field inflation, it was noticed that the correlator with cubic self interaction of a bulk massive scalar can greatly enhance the size of the correlation function. In such models, tree-level graphs exchanging more than one massive scalar make dominant contributions to the 3-point correlator [3,19,21]. However, due to the technical complications, explicit analytical results for inflation correlators with more than one bulk massive field are still beyond our reach at the tree level.
It may come as a surprise to flat-space field theorists that scalar tree graphs are hard to compute. Indeed, setting aside the issues of tensor and flavor structures, the complexity of a scalar Feynman graph in flat spacetime largely increases with the number of loops L: Each loop gives rise to a loop momentum integral, and carrying out these loop integrals are not trivial. However, so long as we stay at the tree level (L = 0), Feynman graphs are simply products of propagators and vertices and are typically rational functions of external momenta. So, increasing the number of vertices and propagators does not generate any difficulty per se.
Things are a little different in inflationary spacetime: Here we normally have full spatial translation and rotation symmetries, but the time translation is usually broken. Accordingly, we Fourier transform only spatial dependence of a function to momentum space, and leave the time dependence untransformed. In this hybrid "time-momentum" representation, we get additional time integrals at all interaction vertices in the Schwinger-Keldysh (SK) formalism [19,[109][110][111][112]. As a result, the complexity of graphs in inflation increases in two directions: either with the number of loops, or with the number of vertices.
Partly for this reason, full analytical computation of tree correlators with multiple massive exchanges remains challenging: In a tree graph, the number of bulk vertices is always equal to the number of bulk propagators plus 1. Thus, a tree graph with I internal legs requires time integrals of (I + 1) layers. Worse still, each bulk propagator D ab (k; τ 1 , τ 2 ) in SK formalism comes in four types, depending on the four choices of SK indices a, b = ± at the two endpoints. The two propagators with same-sign indices D ±± (k; τ 1 , τ 2 ) involve expressions that depend on the time ordering, which make the (I + 1)-layer time integral nested. That is, the integration limit of one layer could depend on the integration variable of the next layer. So, the integration hlquickly becomes intractable with increasing number of bulk lines or vertices.
One may wish to bypass the difficulty of bulk time integrals by taking a boundary approach. For instance, one can try to derive differential equations satisfied by the correlators starting from simple bootstrapping inputs [63-65, 70, 71, 89, 91]. As explored in many previous works, this approach turns out quite successful for single massive exchanges, where the "bootstrap equations" are usually a simple set of second-order ordinary differential equations and usually have wellknown analytical solutions. However, when one goes to the two massive exchanges, the resulting differential equations become much more complicated, and it seems rather nontrivial to directly solve such equations [70].
One can also try other methods such as a full Mellin-space approach, where one still works in the bulk, but rewrites correlators in Mellin space [82][83][84][85]. Then, the time ordering of the same-sign propagators D ab becomes an overall cosecant factor that nests two Mellin variables. While this is enormously simpler than the time-momentum representation, eventually we need to transform the Mellin-space correlators back to a normal time-momentum representation and push the time variables to the boundary: The future boundary is where the observables are naturally defined, and the momentum space is where the cosmological data are presented and analyzed. However, the nested Mellin variables make the inverse Mellin transform nontrivial. Thus, in a sense, in the full Mellin-space approach, we are moving the difficulty of nested time integral to the difficulty of nested Mellin integral.
There are other studies considering inflation correlators or wavefunction coefficients with multiple massive bulk lines. Rather than full analytical computations, most works focused on general properties of such amplitudes, such as the analyticity, unitarity, causality, cutting rules, etc. There is a special case where one does achieve full results for tree graphs with arbitrary number of bulk lines, namely when the bulk field's mass is tuned to the conformal value m 2 = 2H 2 and all couplings are dimensionless. In such cases, the amplitudes reduce to the flat-space results, and one can find nice recursion relations to directly build arbitrary tree amplitudes or even loop integrand [102,103]. However, this result only applies to very special class of theories which are not of direct interest to CC physics. One might want to restore general mass and couplings by integrating the conformal-scalar amplitudes with appropriate weighting functions. However, the complication here is that we encounter fractions of nested energy variables which are hard to integrate.
As we see, no matter what representation we take, there is always a nested part of the amplitude that makes the computation difficult. There is a physical reason behind it: The nested time integrals are from the time ordering of the bulk propagator, and the time-ordered bulk propagator is a solution to the field equation with a local δ-source. Thus, the nested part of the amplitude is closely related to the EFT limit where several or all nested vertices are pinched to a single bulk vertex. Very schematically, we can express this fact with the position-space Feynman propagator D(x, y), which is a solution to the sourced equation of motion (□ x −m 2 )D(x, y) = iδ(x−y). Then, we can make an EFT expansion of D(x, y) ∼ i □x−m 2 δ(x − y). The leading order term is simple, which is just the contact graph with D(x, y) ∼ −iδ(x − y)/m 2 . However, there are higher order terms coming from acting on powers of □ x /m 2 on δ(x − y), which produce a series of momentum ratios when transformed to the momentum-space representation. Technically, as we shall see, such series are typically multi-variable hypergeometric series which in general do not reduce to any well known functions. So, one just has no way to get around with this result; the complication has to show up somewhere. The best we can do is to find a way to write down the analytical result as a convergent hypergeometric series for some kinematic configurations, and then try to find ways to do analytical continuation for other configurations. This is the goal we are going to pursue in this work. Below, we introduce the main results of this work before detailed expositions in subsequent sections.
Summary of main results. In this work, we tackle the problem of analytically computing tree-level inflation correlators with arbitrary number of massive exchanges, via a standard bulk calculation in the SK formalism. The main technical tool is the partial Mellin-Barnes (PMB) representation proposed in [88,89]. The basic idea is very simple: One takes the Mellin-Barnes representation for all factorized bulk propagators, but leaves all the bulk-to-boundary propagators in the original time-momentum representation. Also, one leaves all the time-ordering Heaviside θ-functions untransformed. In this way, one takes the advantage of Mellin-Barnes (MB) representation that it resolves complicated bulk mode functions into simple powers, but still retains the explicit time-domain representation for external modes. As has been shown in several previous works, the PMB representation is suitable for analyzing a range of problems related to inflation correlators, including explicit results at tree and loop levels [88,89], and the analytical properties and on-shell factorizations for arbitrary loop correlators [92,93]. The general procedure of using PMB representation to compute an arbitrary tree-level inflation correlator is detailed in Sec. 2.
As mentioned above, the time orderings are not removed in the PMB representation. So, we still need to deal with them. We solve this problem in Sec. 3. As we will see, the PMB representation greatly simplifies the integrand of nested time integrals. As a result, the most general nested time integral we have to compute takes the following form: where we have time integrals at V vertices, nested arbitrarily by the Heaviside θ functions from the I internal lines. While this integral is still somewhat complicated, it is already in a form that allows us to directly write down the analytical answer. The way to make progress is to recognize that every bulk propagator has a time ordering in a fully nested integral, and we are free to flip the direction of time orderings using a simple relation of the Heaviside θ-function, so that any nested integral can be recast into a partially ordered form. To explain the partial ordering, we adopt this convenient terminology: Whenever we have a factor θ(τ i − τ j ), we call τ j to be τ i 's mother and τ i to be τ j 's daughter. Then, a partially ordered graph simply means that every time variable in the graph can have any number of daughters but must have only one mother, except the earliest member, who is motherless. In plain words, a partially ordered graph can be thought of as a maternal family tree. After rewriting a given nested integral into a partially ordered form, we get new terms with less layers of nested integrals, which can be further rewritten into partially ordered form with additional terms generated. This procedure can be carried out recursively, until all nested integrals are partially ordered. This procedure has a very similar structure with the conventional cluster decomposition in statistical mechanics or quantum field theory. We will call it family-tree decomposition. Then, each of the partially ordered nested integrals is a family tree, which we also call "family" or "family integral" for short. A family is denote by C q 1 ···q N (E 1 , · · · , E N ). The details of this family-tree decomposition will be presented in Sec. 3.1. In practice, this family-tree decomposition takes a very simple form. An example of family-decomposing a graph with 5-layer nested integral is shown in Fig. 1.
As we shall see, the partial order structure allows us to find a simple one-line formula for general family integral C q 1 ···q N (E 1 , · · · , E N ). Working in the configurations where E 1 ≫ E i with i = 2, · · · , N , we find: Here the hatted energy E 1 denotes the maximal energy, which sits at the vertex with the earliest time. On the right hand side, we have N − 1 layers of summations corresponding to the N − 1 descendents of E 1 -site. We use shorthands such as q 1···N ≡ q 1 + · · · + q N . The quantity q j means to take the sum of all q i where either i = j or i is a descendent of j. n j is similarly defined. Explicit application of this formula to the 5-layer graph in Fig. 1 is given in (25)- (28). We give many examples and also a general proof of the formula (2) in Sec. 3.2 and Sec. 3.3. One important point is that the maximal energy variable can be chosen at will: To take the analytical continuation of (2) to kinematic regions where E 1 is no longer maximal, all we need to do is to rearrange the original integral into a different partial order such that the new maximal energy sits at the earliest time. Thus, our method provides a practical way to do analytical continuation of multi-variable hypergeometric series (2) beyond its convergence region. An example of this analytical continuation is given for the example of 5-layer integral in Fig. 2. As will be shown in Sec. 3.4, one can exploit the flexibility of MB representation to rewrite (2) as Taylor series of the sum of several or all energy variables, which further extends the domain of validity of our expressions.
With the formula for general time integrals at hand, the computation of the tree-level inflation correlators becomes a matter of collecting appropriate Mellin poles in the PMB representation. As a demonstration of this procedure, we compute the general tree-level graphs with two bulk massive exchanges in Sec. 4 and present the full analytical result of this type of correlators for the first time. In Sec. 5, we show how to take folded limits of these results, by computing a treelevel 4-point graph with two massive exchanges. We conclude the paper with further discussions in Sec. 6. Useful mathematical formulae on Mellin-Barnes representations and hypergeometric functions are collected in App. A, and some intermediate steps of computing graphs with two massive exchanges are collected in App. B.
Notation and convention. We work in the Poincaré patch of the dS spacetime with inflation coordinates (τ, x) where τ ∈ (−∞, 0) is the conformal time, and x ∈ R 3 is the comoving coordinate. In this coordinate system, the spacetime metric is ds 2 = a 2 (τ )(−dτ 2 + dx 2 ), where a(τ ) = −1/(Hτ ) is the scale factor, and H is the inflation Hubble parameter. We set H = 1 throughout this work for simplicity. We use bold letters such as k to denote 3-momenta and the corresponding italic letter k ≡ |k| to denote its magnitude, which is also called an energy. The energies are often denoted by E i , and the energy ratios such as ϱ ij ≡ E i /E j are often used. We follow the diagrammatic methods reviewed in [19] to compute inflation correlators in SK formalism. We often use shorthand for sums of several indexed quantities. Examples include: Finally, the Mellin integral measures are very often abbreviated in the following way:

Tree Graphs with Partial Mellin-Barnes
In this section, we review the method of PMB representation for a general tree-level inflation correlator with arbitrary massive exchanges. Our starting point is a general B-point connected equal-time correlation function of a bulk field φ in the late time limit: As shown above, the correlation function is defined as an equal-time expectation value of the product of B operators φ k i in 3-momentum space, over a state |Ω⟩ which is taken asymptotic to the Bunch-Davies vacuum state in the early time limit τ → −∞. We assume that the bulk theory of φ is a weakly coupled local quantum field theory. Therefore, after stripping off the momentum-conserving δ-function, the amplitude on the right hand side T (k 1 , · · · , k B ) can be represented as an expansion of connected graphs G(k 1 , · · · , k B ) with increasing number of loops. Thus, the leading contribution is from the tree graphs, which are the focus of this work. We do not specify the type of the field φ, but we do assume that it has a simple mode function φ(k, τ ). More explicitly, if we expand the mode φ k in terms of canonically normalized creation and annihilation operators a k and a † −k , we get mode function φ(k, τ ) as the coefficient: We suppress helicity indices if there are any. We assume that all the time dependence of the mode function φ(k, τ ) can be expressed as an exponential factor e −ikτ times a polynomial of −kτ . This covers essentially all cases relevant to cosmological collider phenomenology where the mode function survives the late-time limit, including the massless spin-0 inflaton field and the massless spin-2 graviton. For instance, the mode function for the inflaton is given by: Our assumption also covers the case where the mode does not survive the late-time limit but is of theoretical interest, such as a conformal scalar ϕ c with mass m c = √ 2 in 3 + 1 dimensions, whose mode function is: The bulk fields appearing in the tree graphs of T (k 1 , · · · , k B ) can be rather arbitrary. In general, they can have arbitrary mass and spin. They can also have dS-boost-breaking dispersion relations, and thus can have nonzero (helical) chemical potential or non-unit sound speed. They can also have rather arbitrary couplings among themselves and to the boundary field φ. In particular, these couplings can break dS boosts and even the dilatation symmetry. However, we do assume that these couplings are well behaved in the infrared so that the diagrammatic expansion remains perturbative in the late-time limit. However, for definiteness, we shall take a fixed type of bulk field, namely a scalar field in the principal series (i.e., with mass m > 3/2), in all the following discussions. Generalization to other cases should be straightforward. For a massive scalar with m > 3/2, it is convenient to introduce a mass parameter ν ≡ m 2 − 9/4. Then, according to the SK formalism [19], we can construct four bulk propagator D ( ν) ab (k; τ 1 , τ 2 ) with a, b = ± for such a field. More explicitly: Then, a general tree graph consisting of massive scalar bulk propagators and massless/conformal scalar bulk-to-boundary propagators can be computed by an integral of the following form: Here we assume that there are V vertices and I bulk propagators in the graph. For each vertex, we have an integral over the conformal time variable τ ℓ (ℓ = 1, · · · , V ). Also, we introduce a factor of ia ℓ as required by the diagrammatic rule [19], and a factor of (−τ ℓ ) p ℓ to account for various types of couplings as well as power factors in the external mode function (such as the ikτ term in the massless mode function (7)). The exponential factor e ia ℓ E ℓ τ ℓ comes from the external mode function, and E ℓ represents the sum of magnitudes of 3-momenta of all external modes. Following the terminology in the literature, we call it the energy at the Vertex ℓ. However, we note that E ℓ is not the total energy at Vertex ℓ since we do not include energies of bulk lines. For each bulk line, we have a bulk propagator D a i1 a i2 (K i , τ i1 , τ i2 ) with momentum K i , which is completely determined by external momenta via the 3-momentum conservation at each vertex. The two time variables τ i1 , τ i2 as well as the two SK variables a i1 , a i2 should be identified with the corresponding time and SK variables at the two vertices to which the bulk propagator attach. The computation of the integral (12) is complicated by the products of Hankel functions, as well as the time-ordering θ-functions in the bulk propagators. To tackle these problems, we use the MB representation for all the bulk propagators, but leave all the bulk-to-boundary propagators untransformed. This is the so-called PMB representation [88,89]. The MB representations for the two opposite-sign bulk propagators (9) and (10) are given by [88,89]: This follows directly from the MB representation of the Hankel function (129), which we collect in App. A. In particular, the Mellin variable s is associated with time τ 1 ands is associated with τ 2 . The same-sign propagators D ±± are obtained by substituting in the above expression into (11). We note that the time-ordering θ-functions are left untransformed. After taking the above PMB representation, the original SK integral (12) becomes: Here we have switched the order of the time integral and the Mellin integral, assuming all integrals are well convergent. With this representation, we see that all SK-index-dependent part goes into the time integral, namely the second line of the above expression. In this time integral, we have used a shorthand (−τ ℓ ) p ℓ −2 ℓ s where ℓ s denotes the sum of all Mellin variables associated to τ ℓ , and the Mellin variables in this summation can be either barred or unbarred. An important fact we shall use below is that the Mellin variables always appear with negative signs in this exponent. Also, we have introduce a function N a 1 ···a V (τ 1 , · · · , τ V ; {s,s}) to represent all combinations of time-ordering θ-functions, as well as the SK-index-dependent phase factor e ∓iπ(s−s) in (13). The reason we introduce the PMB representation is that the time integral now only involves exponentials and powers in its integrand, as shown in the second line of (14). This is significantly simpler than the original time integral, which involves Hankel functions. While this simplification is powerful enough for a single layer time integral, the computation of time-ordered integrals remain nontrivial. In previous works using PMB representation, only the two-layer nested integral was explicitly computed [89]: where 2 F 1 is the dressed hypergeometric function, defined in App. A. For computing inflation correlators with single massive exchange, this result is enough. However, if we wish to go beyond the single massive exchange and consider the most general tree graphs, it is necessary to tackle the problem of computing time integrals of exponentials and powers with arbitrary layers and arbitrary time orderings. We will systematically solve this problem in the next section. From (14), we see that, if the time integral in the second line can be done, then it only remains to finish the Mellin integrals. This is typically done by closing the Mellin contour and collecting the residues of all enclosed poles. So, we need a knowledge about the pole structure of the Mellin integrand. Although the answer to the time integral was not explicitly known in the previous studies, it was proved in [92] that such time integrals, however nested, only contribute right poles of the Mellin integrand. That is, their poles only appear on the right side of the integral contour that goes from −i∞ to +i∞. As a result, all left poles are contributed by the Γ-factors from the bulk propagators, shown in the first line of (14). These are all the poles of the Mellin integrand for a tree graph.
Another important observation is that, if we sum the arguments of all Γ-factors in the first line of (14), we get: That is, all Mellin variables are summed together, with an overall coefficient +2. Here "· · · " denotes s-independent terms, which are irrelevant to our current argument, and happen to be 0 in this particular case. On the other hand, as we shall see from the explicit results in the next section, the right poles contributed by the time integrals are also from Γ-factors of the form Γ[· · · − 2 s]. If we sum over the arguments of all right poles, we will get: which is exactly the s terms from the left-pole Γ-arguments with an additional sign. In this sense, we say that the Mellin variables in the integrand are balanced. In such a balanced situation, the convergence of the Mellin integral is determined by the power factors such as (K i /2) −2s iī in (14). Typically, one can first work in the kinematic region where the internal momenta K i are small (compared to relevant external energies), so that the Mellin integrals will be convergent if we pick up all the left poles, which are all from the bulk propagator Γ-factors in (13). Their poles and residues are well understood. So, if we can finish the time integral, then we only need to collect all left poles from the first line of (14). The result will be a series expansion in K i . So, this result will be valid at least when the bulk momenta K i are not too large. 1 In the opposite limit, when K i becomes large compared to the relevant energy variables, we can instead close the Mellin contour from the right side and pick up all the right poles. In this way, we get analytical continuation of the result from small K i region to large K i region. This will cover most of the parameter space of interest. The narrow intermediate region will be difficult to be expressed by a series solution. Analytically, one needs take the analytical continuation of the series solutions for those intermediate regions, which is a separate mathematical problem. Practically, however, we can use numerical interpolation to bridge the gap between different regions. This strategy has been shown to be workable in previous studies [90]. So, barring possible issues of analytical continuation for special configurations, we can say that, the problem of analytical computation of arbitrary tree-level inflation correlators is solved, if we can compute the arbitrary nested time integral. We will solve the latter problem in the next section.

Time Integrals with Partial Mellin-Barnes
In this section we provide a systematic investigation of arbitrary nested time integrals in the PMB representation. It is clear from the previous section that the most general nested time integral has the following form: Here we are again considering a V -fold time integral with arbitrary nesting. We require that all τ i (1 ≤ i ≤ V ) appear in the θ-factors so that the integral is fully nested. Also, we have used a factor (−τ ℓ ) q ℓ −1 to account for a variety of external modes and couplings, as well as powers of time from the partial MB representation. In the notation of the previous section, we have: The difficulty with time ordering in (18) is easy to understand: A single time integral of exponential with power factors from τ = −∞ to τ = 0 gives rise to a Γ function. However, if there is a time ordering, the integration limit for one time variable would be dependent on another integration variable. As a result, we get incomplete Γ functions after finishing one layer of integral. Then we need to perform time integrals over incomplete Γ functions with integration limits dependent on yet another time variable. This quickly becomes intractable with an increasing number of nested layers.
Our strategy of solving this problem is again the Mellin-Barnes representation: whenever we perform a layer of nested time integration, we take the MB representation of the result so that the integrand for the next layer is still a simple exponential times a power. In this way, the nested time integrals can be done recursively layer by layer, until the last layer, which yields a simple Γ factor. Along the way, we generate many layers of Mellin integrals, which can again be done by closing the contours properly.
As we shall see below, this recursive integration is easiest if the time integral is nested with a partial order, which is not the case in the most general nested integrals. Thus, we should first use a simple relation θ(τ j − τ k ) + θ(τ k − τ j ) = 1 to reorganize the original time integral such that the result is either partially ordered or factorized. This will be called a "family-tree decomposition." Then, we apply the above procedure to the partially ordered integrals to get the explicit results for them. These steps will be carried out in details below. A side remark on notation and terminology: It will be helpful to use a diagrammatic representation for the nested time integral (18). We will use a directional line to denote a θ-function where the direction of the arrow coincides with the direction of the time flow. Two factorized time variables (which is simply associated with a factor of 1) may be connected by a dashed line. So, for instance, we can write the relation θ(τ 1 − τ 2 ) + θ(τ 2 − τ 1 ) = 1 as: Also, to highlight the fact that these diagrams are not the original SK graphs for the inflation correlators, we will use "site" in place of "vertex," and use "line" in place of "propagators." Then, each site τ i is associated with an energy variable E i and an exponent q i , as is clear from (18).

Family-tree decomposition of nested integrals
Now, we describe our family-tree decomposition algorithm in detail. We begin with the most generfamilyal nested time integral (18). After finishing all the time integrals, the result In the following, we shall show that this integral can always be written as a sum over a finite number of terms. Each term is a product of several families. Each family is a multi-variable hypergeometric function of several energy variables E i . Of course, multi-variable hypergeometric functions are not well studied. It is most useful if we can find a fast converging series expansion of this hypergeometric function in terms of any given small energy ratios. Below, we will show that this can be done.
The reduction procedure. Our reduction procedure consists of the following simple steps: Step 1: We start with a particular kinematic region of the integral We want to find an analytical expression for T q 1 ···q V (E 1 , · · · , E V ) as a series in 1/E i , which should be convergent in most of the region where E i remain the largest. We add a hat to the largest energy variable E i to highlight the fact that we are considering a particular kinematic region.
So, if we choose E 1 to be the largest energy, we will write T q 1 ···q V ( E 1 , E 2 , · · · , E V ) to highlight this choice. If, instead, we want to consider the case where E 2 gets larger than E 1 , then we should add the hat on E 2 . The degenerate case where there are multiple maximal energies will be considered in following subsections.
Step 2: We use the relation θ(τ j − τ k ) + θ(τ k − τ j ) = 1 to flip the direction of time flows in some bulk lines, such that the original graph is broken into a sum of several terms. Each term can be represented as a graph, in which all sites are either partially ordered or factorized. As a result, each graph becomes a product of several integrals, each of which has a partial order structure, and is called a family. As a part of the rule, we require that the maximal energy site has the earliest time in a family.
Let us define what is a (partially ordered) family. Clearly, a time-ordered line connects a site with an earlier time to another site with a later time. We call the earlier-time site the mother of the later-time site, and call the later-time site the daughter of the earlier-time site. Then, a partially ordered graph means that every site has a unique mother, except the maximal-energy site, which is the earliest-time site and motherless. On the other hand, a mother can have many daughters.
In this way, all sites within a family integral genuinely belong to a family. Also, for a given site, we call all the sites flowing out of it the descendant sites. Thus, the descendant sites of a given site consist of its daughters, granddaughters, great-granddaughters, etc.
Let us rephrase the above heuristic language into a more rigorous definition. We will use C q 1 ···q N ( E 1 , E 2 , · · · , E N ) to denote a family integral with N sites, where we have highlighted the maximal energy E 1 with a hat. Then, a family integral C q 1 ···q N ( E 1 , E 2 , · · · , E N ) has the following form: with the following restrictions on the θ-function factors: 1. Every time variable τ i (1 ≤ i ≤ N ) appears in time-ordering θ-function. (All sites belong to a family.) 2. In a factor such as θ(τ j − τ k ), let us call τ j to be in the late position and τ k in the early position. Then, it is required that every variable τ i except the maximal energy site appears in the late position once and only once. (Every site has a unique mother except the maximal energy site.) On the other hand, early positions can be taken more than once by a given τ i . (A mother can give birth to more than one daughter.) 3. The maximal energy site τ 1 appears in θ factors only in the early position. (The maximal energy site is motherless, but can have any (including zero) number of daughters.) Step 3: After taking Step 2, each resulting graph is a product of several fully factorized families. The maximal energy site sits in a particular family, which we call the maximal-energy family.
As a consequence, families other than the maximal-energy family are independent of the maximal energy variable E i , and it becomes meaningless to ask for a series expansion in E i for those families. We call them non-maximal energy families. Thus, for each of the nonmaximal energy families, we should further assign a "locally" maximal energy, such that this energy is largest among all energies within the family. Then, we further perform the reduction of Step 2 for all non-maximal energy families and we do this procedure recursively, until, within each family, the locally maximal energy site sits at the earliest time.
Step 4: After taking the above steps, we fully reduce the original integral T q 1 ···q V ( E 1 , E 2 , · · · , E V ) into a sum of products of partially ordered families, and in each family, the locally maximal energy acquires the earliest time.
It then remains to state the rule for directly writing down the answer for arbitrary families. The rule is the following: 1) Within each partially ordered family, we assign a summation variable n i for all sites except the (locally) maximal-energy site. Without loss of generality, we can always relabel the sites within a family such that the (locally) maximal energy is E 1 . Then, for the N -site family defined in (21), the result is: Here, the hatted energy E 1 represents the maximal energy. In the first line, we stripped away a dimensionful factor (iE 1 ) q 1···N so that the resulting integral C q 1 ···q N ( E 1 , E 2 , · · · , E N ) is dimensionless. In the second line, we have defined ϱ jk ≡ E j /E k . Also, n i is defined to be the sum of n i -variables over the site i and all its descendants. q is similarly defined. This completes our reduction of the original time integral into a sum of products of hypergeometric series.
Example. As often happens, it is better to demonstrate an algorithm with examples than mere abstract description. So, now, let us demonstrate the above reduction procedure with a concrete example. Suppose we want to compute a 5-layer time integral: Furthermore, suppose that we want to consider the kinematic region where E 1 is the largest energy among all five energies. Thus, we want to express the final result as a series expansion of 1/E 1 . This is shown on the left hand side of Fig. 1, where the magenta-circled site represents the maximal-energy site. Then, according to the above procedure, we should use the relation θ(τ i − τ j ) + θ(τ j − τ i ) = 1 to change the direction of several lines, such that 1) all sites become either partially ordered or factorized and 2) the maximal-energy site E 1 has the earliest time variable. This is done on the right hand side of Fig. 1. In each diagram on the right hand side, we get a product of one or several partially ordered families.
In all but the last term on the right hand side of Fig. 1, we have families which do not contain the maximal energy site. Thus we should specify locally maximal site for each of them. The one-site family is trivial. The nontrivial non-maximal families appear in the first and third terms on the right hand side of Fig. 1, which can be expressed as C q 3 q 4 (E 3 , E 4 ) and C q 3 q 4 q 5 (E 3 , E 4 , E 5 ), respectively. Thus, we should further assign a maximal energy for these two sites. So, let us further work within the region where E 3 > E 4 , E 5 , so that E 3 is the locally maximal energy, marked with a blue circle in Fig. 1. (On the other hand, the relation between E 3 and E 2 is irrelevant.) Then, we see that E 3 is already in the earliest time site in both families. So, we are done, and the result of our reduction procedure can be expressed as: Here we have also added hats to the locally maximal energy E 3 . In the last term, we added a superscript (iso) to show that this family has a cubic vertex. See the next subsection for details. Next, we assign n 2 , · · · , n 4 for the sites with energy variables E 2 , · · · , E 4 , respectively. It is clear from (24) that there are four independent nontrivial families (i.e., with more than one site)  (24), showing the reduction of a 5-layer time integral into partially ordered families. In this example, we choose The maximal energy site (Site 1) is marked with a magenta circle and the locally maximal energy site (Site 3) is marked with a blue circle.
involved in this example. Applying the formula (22) for each of them, we get On the other hand, the result for the one-site family is trivial: C q (E) = Γ(q). In fact, some of above series can be summed to well known hypergeometric functions, which we shall introduce below. In any case, we have found the series expression for the original 5-layer time integral T q 1 ···q 5 (E 1 , · · · , E 5 ) without actually doing any integrals. The above series solution has a validity range beyond which the summations no longer converge. This happens in particular when any energy E i (i = 2, 3, 4, 5) becomes larger than E 1 . In principle, if we need the result when E 1 is no longer maximal, we need to take analytical continuation of the above series. This analytical continuation can be very conveniently implemented in our procedure. To see this, let us have a second look at the 5-site integral T q 1 ···q 5 (E 1 , · · · , E 5 ) in (23), but now choose E 3 as the maximal energy. Then, according to our procedure, we should do a new family-tree decomposition, as shown in Fig. 2.
Clearly, we do not need to choose any locally maximal energy in this example. The explicit expressions for the above families can be written down directly according to the general formula (22): Thus we have find an expression for the original 5-site integral T q 1 ···q 5 (E 1 , · · · , E 5 ) expanded as powers of 1/E 3 . Let us emphasize that (24) and (29) are just different expansions for the same function T q 1 ···q 5 (E 1 , · · · , E 5 ), with different validity regions.

Partially ordered families: simple examples
Clearly, the only nontrivial step in our family-tree decomposition procedure is the last step, where we directly write down the answer for the family integral (22). The derivation of this result is best illustrated with examples. So in this subsection we will walk the readers through a few simple examples, before presenting a general proof in the next subsection.
One-site family. We begin with the simplest integral, the one-site family, shown in Fig. 3(a): The application of the rule is trivial, and we have the following answer: The answer is obtained by a direct integration of (34). Since there is only one dimensionful variable E involved in the problem, the final answer for the dimensionless family C q (E) must be independent of E.
Two-site family. Next let us look at the simplest nontrivial example, namely the two-site family, shown in Fig. 3(b). The integral is: By design, we take E 1 > E 2 . Now let us try to find the answer for the above integral. It turns out useful to start from the integral of reversed time ordering: Then, the integral over τ 2 can be performed, with the result expressed in terms of an exponential integral E p (z) whose definition is given in (130): At this point, we make use of the following MB representation of E p (z): The details of this MB representation is given in App. A. As explained there, the pole in s from the denominator 1/(s + p − 1) should be interpreted as a left pole, in the sense that the integration contour should go around this pole from the right side. Now, using (39) in (38), we get: Then, the τ 1 integral is trivial, which is simply given by C q 12 −s 2 (E 1 ) = (iE 1 ) −q 12 +s 2 Γ(q 12 − s 2 ). So, finishing the τ 1 integral, we get: Now it remains to finish the Mellin integral over s 2 . Given that ϱ 21 = E 2 /E 1 < 1, we should close the Mellin contour from the left side and collect the residues of all left poles. There are two sets of left poles, one at s 2 = −n 2 with n 2 = 0, 1, 2, · · · , which is from the Γ-factor Γ(s 2 ), and the other at s 2 = q 2 coming from the denominator. Collecting the residues at these poles, we get: Now, we recognize that the last term without any summation is the product of two one-site family: Then, given the relation: we see that the original family integral (36) is: This is exactly what we would get using the rule (22). Incidentally, the above summation can be directly done and the result is the well-known Gauss's hypergeometric function: Here we use the dressed version 2 F 1 instead of the original hypergeometric function 2 F 1 for notational simplicity. The dressed hypergeometric functions are defined in App. A. Figure 4: Two independent family integrals at 3-site level. Now, had we chosen to expand the integral (36) in terms of 1/E 2 , we will get: Clearly, the series expression in the first line in (47) has a different region of convergence from the series expression in (45). However, the two expressions are just two power-series expansions of the same function in two different limits, one at E 2 /E 1 → 0 and the other at E 1 /E 2 → 0. This becomes more transparent after finishing the summations of both series into hypergeometric functions. Indeed, equating (46) with the second line of (47), we just get a transformation-ofvariable formula for the hypergeometric function. Thus, our procedure provides a convenient way to derive many transformation-of-variable formulae for hypergeometric functions, which is particularly convenient for more complicated hypergeometric series, as we shall see below.
Three-site family. Next we consider a slightly nontrivial case with three sites. There is only one topology at the tree level with 3 sites. However, after including the time ordering, there are two independent possibilities, depending on whether the latest site is on the side or in the middle. These two possibilities are shown in Fig. 4. Again, by construction, the latest site is chosen to be the maximal-energy site. So, for the case in Fig. 4(a), we have: We again start from the completely reversed integral: Now, we can repeat the above strategy and finish the three layers of time integrals in the order of τ 3 , τ 2 , τ 1 . The first two layers produce exponential integrals which can then be represented as Mellin integrals. The last layer is again a single-site integral which can be finished directly. Here we show the results after finishing each layer of time integral and taking the MB representation for exponential integrals: We start to observe a pattern here: Let the maximal-energy site be τ i . When carrying out any but the last layer of time integrals, say τ j with j ̸ = i, we are effectively generating a new layer of Mellin integral with Mellin variable s j , a pole-generating factor Γ(s j )/( s j − q j ), and a power of energy ratio ϱ −s j ji . Here s j is the sum of all Mellin variables assigned to the site j and its descendant, and q j is likewise defined.
Then it remains to carry out the Mellin integrals. Unlike the previous case, now we encounter pole-carrying factors involving more than one Mellin variable. In the current case, it is the denominator 1/(s 23 − q 23 ). To avoid any potential complication of such poles, our strategy is that we perform the Mellin integrals in the "anti-chronological" order. In the current case, we integrate out s 3 first, by collecting poles only from Γ(s 3 )/(s 3 −q 3 ). Only after this is done, we then perform the s 2 -integral, by collecting poles from Γ(s 2 )/(s 23 − q 23 ). By this time, the s 3 variables in these factor have already been set to the poles. Thus, we never need to directly deal with poles involving a sum of several Mellin variables. Finishing the Mellin integral in this way, we get: Here we have restored all the dimemsionful energy factors to make clear the following point: The result of Mellin integral is effectively executing the identical transformation θ(τ i −τ j ) = 1−θ(τ j −τ i ) in a line-by-line fashion. Thus, with N lines in a family, we will get 2 N terms. All but one of them are factorized. There is a unique unfactorized term with all line reversed. In the current example it is the first term on the right hand side of (51). This is nothing but the original family integral. Thus: Once again, this is exactly what we would get by applying the simple formula (22). It seems to us that this series does not sum to any widely known special function in general, but it can be represented as a (dressed) Kampé de Fériet function, whose definition is collected in App. A: The lesson we learned from the above example is the following: To find the answer to a given family integral C, all we need to do is to compute another integral R[C] with all time orderings completely reversed. We compute R[C] layer by layer. Each step generates an exponential integral of which we take the MB representation. The last layer of time integral is directly done, and we are left with an (N − 1)-fold Mellin integral. We finish the Mellin integral by retaining poles only from Γ-factors. The result of this term is automatically a sign factor (−1) N −1 times the original family C.
With this lesson learned, we can bypass all steps detailed above, and write down the answers for arbitrary families. Now, let us go on to consider the three-site family in Fig. 4(b), which corresponds to the following integral: The result after finishing all three layers of time integrals for the reversed diagram R[ C] is: Then, we finish the Mellin integral by picking up poles in Γ[s 1 , s 3 ] only. Multiplying the result by a trivial sign factor of (−1) 3−1 = 1, we get the original family: Again, it agrees with what we would get by applying (22). Incidentally, the above two-fold hypergeometric series belongs to the well-known Appell series, which can be summed into the (dressed) Appell F 2 -function: The definition of F 2 is collected in App. A.
Four-site family with a cubic vertex. Finally let us look at four-site families. There are two possible tree topologies with 4 sites. One is the chain graph C (n) q 1 q 2 q 3 q 4 , which is a direct generalization of the case considered above. We are not going to consider this chain graph any more. On the other hand, there is a new topology with a cubic vertex C (iso) q 1 q 2 q 3 q 4 , as shown in Fig.  E As always, we compute the corresponding reversed time-ordering integral: Then, the original integral is obtained by finishing the three-fold Mellin integral in which we only collect poles from Γ[s 1 , s 2 , s 3 ], and multiplying the result by (−1) 4−1 = −1. The result is: (−1) n 123 Γ[n 123 + q 1234 ] (n 1 + q 1 )(n 2 + q 2 )(n 3 + q 3 ) This three-variable series is not covered in commonly known special functions. It is however covered by the so-called (dressed) Lauricella's F A function: The definition of this function is collected in App. A. Next let us look at Fig. 5(b) where the maximal-energy site is on the side. We take it to be E 1 . Then, the corresponding family integral is given by: Without mentioning any details of intermediate steps, we directly provide the final answer to this integral:

General family integrals
Above we have examined enough number of examples. By now, it is clear that why we need to do family-tree decomposition: Our performance of nested integrals requires a partial order of the nested time variables, and this partial order can always be achieved by the family-tree decomposition. Once we have a partially ordered integral, we can always carry out the completely reversed integral, from the originally latest sites (earliest sites in the reversed integral), and to their mothers, and to grandmothers, etc., until the last layer which is the maximal energy site. In this way, a full derivation of the general equation (22) becomes a matter of mathematical induction. Below we complete this proof.
We begin with a general partially ordered family with N sites, C q 1 ···q N ( E 1 , E 2 , · · · , E N ), where we assume that the maximal energy site is τ 1 . Its integral representation is given in (21). As in the previous section, we work with the completely reversed integral R[C q 1 ···q N ( E 1 , E 2 , · · · , E N )], and we integrate every time variable from −∞ to the time variable of her mother.
Suppose that we have finished all the time integrals for the descendants of the site τ j , and now we want to finish the time integral at τ j . Our induction assumption is that, after all the descendants of τ j integrated out, the integral over the variable τ j has the following form: where τ M is the time variable of τ j 's mother, and D(j) denotes the set of labels for all τ j 's descendants. Now, finishing the τ j integral, we get: Here we have used the fact that s j + i∈D(j) s i = s j . Also, we have abbreviated i∈D(j) s i as s i when it appears as an upper or lower index. Now, to go one step further, we should finish the time integral over τ M , which we denote as I M (τ G ) and τ G is the time variable of τ j 's grandmother. To this end, we take products of above integral I j (τ M ) over all daughters τ j of τ M . Then, together with τ M 's own factor (−τ M ) q M −1 e iE M τ M , we get: This is identical to (64) upon a "generation shift" j → M and M → G. So we have shown that the original induction assumption (64) persists to all generations as long as it holds at one generation. On the other hand, it is trivial to check that the induction assumption holds for the initial step, i.e., at any site who has no descendant. Thus, we have proved that the induction assumption (64) holds for all sites. In particular, (64) holds for the maximal energy site τ 1 if we take τ j = τ 1 and τ M = 0. Then, completing this final layer of time integral over τ 1 , we get, for the whole reversed family, As shown many times in the previous subsection, the original family C q 1 ···q N ( E 1 , E 2 , · · · , E N ) is recovered by picking up all poles from Γ[s i ], and including an overall factor (−1) N −1 which comes from reversing the directions of N − 1 bulk lines. Thus, we get: This is exactly the original family formula (22). Thus we have completed the proof.

Alternative representation
The MB representation of a function is not unique. In the previous computations, we have chosen a relatively simple representation (39) for the exponential function E p (z). This representation allows us to find simple expression for the dimensionless family integrals as Taylor series of energy ratios. On the other hand, there exist other MB representations which may be useful in certain cases. One example is the following partially resolved MB representation, which is particularly useful to improve the convergence of the hypergeometric series when there are several energies comparable to the maximal energy: This result can be derived from the MB representation for a confluent hypergeometric function, as discussed in App. A. We note that (69) is not a complete MB representation of E p (z), as there is an exponential factor e −z left. 3 As we shall see, this remaining exponential factor will help us to circumvent the problem of convergence of hypergeometric series with several maximal energies. Let us take monotonic three-site family integral (48) as an example, namely Fig. 4(a), but we do not divide out the dimensionful factor (iE 1 ) q 123 , nor do we assign a maximal energy variable: As before, we compute the integral with all time orderings reversed, but with the new representation (69). The result is: Similar to the Mellin integrals in the previous representation, each Mellin variable s i gets two sets of left poles from the Γ factors, one at s i = −1 − n i (n i = 0, 1, · · · ) from Γ(1 + s i ), and the other more complicated, involving both q i and other Mellin variables from the descendants of s i . We are not going to present a detailed analysis here, but only mention that, similar to the previous case, the original family integral C q 1 q 2 q 3 (E 1 , E 2 , E 3 ) is recovered by picking up poles from all Γ(1 + s i ) factors only and multiplying the result with an appropriate sign factor. Thus: where we have defined ϱ 23T ≡ E 23 /E 123 and ϱ 3T ≡ E 3 /E 123 . Thus we see that, instead of using the inverse of the maximal energy as the expansion parameter, in this representation, we are using the inverse of the total energy E 123 as the expansion parameter. Although it has a somewhat more complicated looking than the previous representation, this representation is a safer choice in certain cases, in particular in the kinematic region of several equal or comparable maximal energies. In any case, it is easy to check numerically that (72) and (52) agree with each other perfectly whenever both series converge. The lesson here is that we can make use of the flexibility of MB representations to get different series solutions for the nested time integrals, expanded either in the inverse power of some energy variable of a given site, in the inverse power of the sum of several energy variables, or even in the inverse power of the total energy. Although the final results may look quite different, these results are just different expansions of the same function. We can thus obtain a large number of transformation-of-variable relations for these multi-variable hypergeometric functions. We leave a more systematic investigation of this topic to future works.

Discussions
We end this section with a discussion of the nested time integral.
Pole structure. In Sec. 2, we mentioned that the time integral in the PMB representation only contains right poles for all Mellin variables, which was proved in [92]. Now that we have explicit results for arbitrary nested time integrals, it is straightforward to check this statement. Indeed, we can rewrite our general result for the family integral (22) in the following way: Then it is clear that all the exponents q i (i = 1, · · · , N ) have the positive coefficients when appearing in the arguments of Γ factors. Now, if we use (19) to rewrite all q's in terms of Mellin variables s, we will see that all Mellin variables s have negative coefficients when appearing in the arguments of all Γ factors in (73). So, we have confirmed with our explicit results that nested time integrals only have right poles in all Mellin variables.
With the explicit results for the nested time integrals, it is also easy to confirm that the Mellin integrand for any tree graph in the PMB representation is well balanced for all Mellin variables, an important fact for the computation of Mellin integrals, as mentioned in Sec. 2. To see this, we only need to derive (17) from our result.
From our result for the time integral in (22), it is trivial to see that the exponents q i (i = 1, · · · , N ) appear in the Γ factor Γ(q 1···N + n 2···N ) as a total sum. Then, let us look at (19), which says that the value of q ℓ at Vertex ℓ receives contributions from all Mellin variables ending at this vertex. Note that, by construction, every Mellin variable is associated to one and only one vertex. Thus, (19) tells us that summing over all q ℓ is equivalent to summing over all Mellin variables. As a result, the argument of the Γ factor Γ(q 1···N + n 2···N ) becomes: which agrees nicely with the general structure in (17). So, we see that the Mellin variables are indeed balanced.
Hard limits. It is interesting to look at different kinematic limits of our result (22). First, it is simple to take a hard limit where one energy E 1 is much greater than any other energies. Obviously, in this limit, we should work with the expression where E 1 is chosen as the maximal energy site. Then, in the series expansion (22), only the leading term with n 2 = · · · = n N = 0 survives the limit. So we get: Apart from the simple numerical factor j 1/ q j , this is very similar to the result of one-site family C q (E) = Γ(q)/(iE) q with E = E 1···N and q = q 1···N . Here we have used the fact that E 1 ≃ E 1···N in the E 1 → ∞ limit. So, in this hard limit, the time integral behaves as if we have pinched all N nested vertices together, with all exponents q i (i = 1, · · · , N ) summed. Thus, this hard limit can in a sense be thought of as an EFT limit, where all internal lines in the family integral shrink into local vertices. From the viewpoint of the cosmological bootstrap [63,89], we know that the EFT part is related to the particular solution of the bootstrap equation with a local source term. This local source term originates exactly from the time-ordering part of the internal propagators. So, there is a close relation between the local EFT limit and the nested integrals, and it is not surprising to get (75). However, we can make a new interesting observation from (75): In the original SK integral (12) for a correlator, we need to sum over all SK indices, which involve all kinds of propagators with arbitrary nesting. This means that the site of E 1 can be nested arbitrarily with other sites. Then, coupled with (75), we see that the E 1 → ∞ limit can generate a power 1/E q 1 ···q N 1 involving exponents q i at any other site. Note that the variable q i contains the Mellin variables of internal lines ending at Site i, and we see that the power 1/E q 1 ···q N 1 could be dependent on Mellin variables of any internal lines not ending at Site 1. Thus, if we finish the Mellin integrals by picking up left poles for those Mellin variables, they can introduce noninteger powers of 1/E 1 . This is exactly the source of local signals. The local signals have been considered mainly for single exchange graphs in previous works [79,89]. Here we see that, in the hard limit E 1 → ∞, the local signal from E 1 can in principle be generated by any internal massive propagators not ending on Site 1. So, the local signal is more subtle and more complicated than the nonlocal signal. This topic will be further explored in a separate work [113].
Soft limits: internal vertices. Now let us look at an opposite limit where one or several energies approach zero, E i → 0. This is a soft limit. Note that, in our general expression for the nested time integral (18), we have assigned an exponential factor e iE i τ i for each site at time τ i . This factor is generally from the bulk-to-boundary propagator of a massless or conformal scalar (or from a massless graviton if nonzero spins are considered). In realistic tree graphs, there are certainly vertices on which only bulk massive propagators end, and no bulk-to-boundary lines are attached. We call such a vertex an internal vertex following [93]. Clearly, we have E i = 0 for such a vertex. So, it is necessary to know how to take soft limits if we want to consider graphs with internal vertices.
Fortunately, our series expression for the time integral makes it very convenient to take a soft limit. For instance, suppose that we want to set E 4 = 0 in the 4-site family (63) which corresponds to Fig. 5(b). Then, the form of the hypergeometric series in (63) allows us to set ϱ 41 = E 4 /E 1 = 0 directly without encountering any singularities. Then, in the summation over n 4 , only the term with n 4 = 0 survives the ϱ 41 → 0 limit, and we get: We take this opportunity to make a general comment on the computation of graphs with internal vertices, as briefly mentioned in Footnote 1. We illustrate the point with a concrete example. Suppose we want to compute the following integral with E 4 = 0: Note that we set E 4 = 0 on the right hand side. Suppose we want the result for this integral with E 1 chosen as the maximal energy. Then, according to our reduction procedure, we should proceed with the following family-tree decomposition: This is shown diagrammatically in Fig. 6. (Note that the last graph in Fig. 6 is exactly the previous example in (76).) Then, applying the general formula (22) to all families here, and setting the summation variable n 4 = 0, we get: (−1) n 23 Γ[n 23 + q 1234 ] (n 2 + q 2 )(n 3 + q 3 )(n 23 + q 234 ) On the other hand, we can as well compute the integral (77) directly. First, we integrate out τ 1 , τ 2 , and τ 3 , and the result is: Now, the final layer integral contains no energy variable since it is from an internal vertex. Finishing this integral, we get a δ function: Now we need to choose a maximal energy. Suppose we choose E 1 without loss of generality. Then, we use the above δ function to integrate out s 1 . Then, (80) becomes: We can finish this integral by collecting the residues of all left poles of the integrand, as before.
The result exactly agrees with (79). There are two lessons to be learnt here. First, when computing a specific nested integral, if we decide to do the time integral directly, we do not have to be as rigid as when we derive the family-tree decomposition. Instead, we can always do the nested integral so long as the integral has a partial order, and the latest (or earliest) site does not have to be the maximal-energy site. The choice of maximal energy can be delayed until we perform the Mellin integral, where we do need a maximal energy to decide how to make the series expansion. On the other hand, the advantage of family-tree decomposition is that we do not have to compute the integral at all; So long as we follow this reduction procedure, we can write down the answer directly.
The second lesson is about the δ function generated from an energy-less time integral, such as the one in (81). When we use (81) to integrate out a Mellin variable, say s 1 , we are effectively setting s 1 = q 1234 − s 23 everywhere in the integrand. Then, all previous left (right) poles of s 1 now become right (left) poles of s 23 . However, this left-right flip is harmless at least for tree graphs. The general rule is that, whenever we have a δ function from an internal vertex, we choose a maximal energy among all energies connected at this vertex, and we use the δ function to integrate out the Mellin variable associated with the maximal energy. Then, we still pick up left poles of other Mellin variables to finish the Mellin integral. In this way, we will end up with a series expansion in terms of small energy ratios, as shown in the above example.
Multiple maximal energies. Finally, there is a more difficult parameter region where the energies at more than one sites become equal or comparable. This case is tractable if the equal energies are not maximal. The only tricky situation is when the equal energies are maximal, so that, in the series solution (22), there is at least one energy ratio ϱ j1 approaching 1. At this point, the series representation is likely divergent. There are several things one can try in this case. First, it is always possible to finish any one layer of summation in the general formula (22) in terms of a (generalized) hypergeometric function p F q . Then, one can study the behavior of this hypergeometric function with argument equal to 1. Such a pure analytical strategy can sometimes be extended to two-variable summation as well. Second, one can switch to the partially resolved representation (69) discussed in Sec. 3.4, so that the result is expanded in powers of 1/E 1···N instead of the inverse of any single energy variable. This helps to improve the convergence of the series in many cases. As mentioned above, this is a practical way to discover many transformationof-variable formula for multi-variable hypergeometric functions, and thus could be particularly useful. Third, when all previous methods fails (such as when all energies become nearly equal), we can do numerical interpolation to sew together disconnected parameter region with convergent series expressions. We leave this somewhat mathematically oriented problem to a future work.

General Two Massive Exchanges
With the nested time integral done in the last section, in principle, we are able to compute arbitrary tree-level inflation correlators with any number of massive exchanges. In this section, we illustrate this procedure with a concrete example, namely a general tree graph with two massive exchanges, as shown in Fig. 7. We follow the diagrammatic representation of [19]. In Figure 7: A general tree graph with two massive exchanges. particular, the external (bulk-to-boundary) propagators can be either conformal scalars with m 2 = 2, massless scalar fields such as the inflaton, or the massless spin-2 graviton. The conformal scalar is technically easiest and is often used as a starting point in a theoretical analysis of inflation correlators. The cases of massless scalar and tensor modes are more relevant to CC phenomenology. On the other hand, the two internal (bulk) propagators represent two massive scalar fields which can be either identical or distinct. There is no difficulty to generalize the bulk lines to massive fields with spins or with helical chemical potentials, but we choose to work with scalars of principle series (m 1 , m 2 > 3/2) for definiteness. Thus, we assign two mass parameters ν 1,2 ≡ m 2 1,2 − 9/4 for the two lines, respectively.

Three-vertex seed integral
Following the diagrammatic rule of SK formalism [19], one can show that the correlators in the form of Fig. 7 can in general be reduced to the following three-vertex seed integral : Here, as before, E i represents the total energies of bulk-to-boundary lines at the Vertex i, while ℓ 1 and ℓ 2 represent the 3-momenta of the two internal lines, respectively. As before, we have included power factors of the form (−τ i ) p i to allow for different choices of external states and coupling types. Finally, the two bulk massive propagators D bc are given in (9), (10), and (11). To minimize unnecessary complications, we will take p 1 , p 2 , p 3 ∈ R. Generalization to complex values of p i is straightforward, although the expressions will be lengthier.
The E and ℓ factors in front of the integral in (83) are included to make the integral dimensionless. The reason we introduce this special combination of energy variables is the following: We can define the dimensionless integration variables z i = E i τ i (i = 1, 2, 3), and use the momentum ratios r 1 = ℓ 1 /E 1 , r 2 = ℓ 1 /E 2 , r 3 = ℓ 2 /E 2 , r 4 = q ℓ /E 3 . Then, one can easily verify: and a similar relation for the ℓ 2 -propagator. Then, the seed integral is manifestly dimensionless and depends only on dimensionless energy ratios: There are simple kinematic constraints for the range of r i variables from the momentum conservation at each vertex. For instance, let there be N external lines ending at Site 1 with 3-momenta k 1 , · · · , k N . Then, by definition, E 1 = k 1 + · · · + k N ≥ |k 1 + · · · + k N | = ℓ 1 . Thus we have 0 < r 1 < 1. Similarly, we have 0 < r 4 < 1. On the other hand, the constraints on r 2 and r 3 are much weaker. In general, these two values can take any nonnegative real values. Many correlators with two massive exchanges can be expressed in terms of I p 1 p 2 p 3 a 1 a 2 a 3 . For example, when the external leges are conformal scalars ϕ c with cubic and quartic direct couplings with two massive scalars σ 1,2 , we can form a 6-point correlator. The Lagrangian is: Here m c = √ 2 is the mass of the conformal scalar ϕ c , and m 1,2 > 3/2 are the mass of two scalars σ 1,2 , respectively. We also include two cubic couplings with dimension-1 coupling constants µ i (i = 1, 2) and a quartic coupling with dimensionless coupling λ. The powers of scale factors a = −1/τ are introduced to make the Lagrangian scale invariant, and the spacetime indices in (86) are contracted by Minkowski metric η µν . Then, the 6-point correlator is shown in Fig. 7 with all black dots removed. With the diagrammatic rule, it is easy to see that the corresponding SK integral reduces to the seed integral in the following way: (−τ f ) 6 k 12 k 34 k 56 64k 1 · · · k 6 ℓ 3 1 ℓ 3 2 a 1 ,a 2 ,a 3 =± where we have introduced a final-time cutoff τ f . We note that this expression is for a single graph G(k 1 , · · · , k 6 ) rather than the whole correlator T (k 1 , · · · , k 6 ) at the same perturbative order. The correlator T can be obtained from the graph G by including suitable permutations, which we do not spell out here.
As another example, we can consider the 4-point correlators of massless inflaton φ with two massive exchanges, as shown in Fig. 9. We assume that the inflaton φ is coupled derivatively, to respect the approximate shift symmetry of the inflaton field, and also to produce a nontrivial result. 4 The relevant Lagrangian is: Then, the 4-point graph in Fig. 9 can be expressed as: Thus the computation of the 4-point correlator (89) requires us to take simultaneous folded limit r 1 → 1 and r 4 → 1 which is a bit nontrivial. We shall take this limit in the next section.

Computing the seed integral
Now we are going to compute the seed integral (83). The computation is rather lengthy and tedious. Here we only outline the main steps, and collect more details in App. B The seed integral I p 1 p 2 p 3 a 1 a 2 a 3 (r 1 , r 2 , r 3 , r 4 ) in (83) has 8 SK branches, depending on the values of the 3 SK indices a 1 , a 2 , a 3 = ±. We only need to compute 4 integrals: I p 1 p 2 p 3 +++ , I p 1 p 2 p 3 ++− , I p 1 p 2 p 3 +−+ , and I p 1 p 2 p 3 −++ . Since all exponents p i (i = 1, 2, 3) are real, the other four can be obtained by taking complex conjugation. As a general rule, for graphs with real couplings, flipping the sign of all SK indices simultaneously brings an integral to its complex conjugate.
As usual, the main difficulty of the computation comes from the time orderings. Thus, our first step is to rewrite the time-ordered propagator D ++ in a more suitable form. For definiteness, we work in the region where E 2 > E 1 and E 2 > E 3 . Then, according to the discussion of the previous section, whenever we have a time ordering between τ 1 and τ 2 , or between τ 3 and τ 2 , we should let τ 2 take the earlier position. Thus, we use the following expression for the two D ++ propagators: After taking this representation, the expressions for the four integrals I p 1 p 2 p 3 +++ , I p 1 p 2 p 3 ++− , I p 1 p 2 p 3 +−+ , and I p 1 p 2 p 3 −++ can be obtained directly. Here we show one example of I p 1 p 2 p 3 +++ . The complete list is given in (147) to (150).
Next, we expand all the integrals and classify the terms according to whether the adjacent two time variables are time-ordered (T) or factorized (F).
The explicit expressions of these integrals are given in (151)-(158). Then, we define the following 4 integrals: After this regrouping of terms, each of the four integrals in {I (FF) , I (FT) , I (TF) , I (TT) } has a definite nesting structure in its time integral, which can then be readily computed using the PMB representation. The procedure is by now standard: We first use the MB representations for the two massive propagators: The assignment of the four Mellin variables s 1 , · · · , s 4 are shown in Fig. 8. Then, the time integrals can be directly done for all branches using results of the previous section. It then remains to finish the integrals over the four Mellin variables s 1 , · · · , s 4 . We work in the region where the two bulk momenta ℓ 1 and ℓ 2 are both softer than all energies E 1 , E 2 , and E 3 . In this region, we have 0 < r i < 1 (i = 1, 2, 3, 4), which means that we should pick up all left poles to finish the Mellin integrals. 5 As already explained in Sec. 2, the left poles of the Mellin integrand are all from the Γ factors in (100) and (101), and are given by: Here n i = 0, 1, 2, · · · (i = 1, · · · , 4), and a i = ± are not SK indices. Thus, by collecting the residues of the integrand at all these poles, we get the final answer for the three-vertex seed integral in (83). Similar to the case of single-exchange graph studied in previous works, it turns out convenient to express the final answer of I p 1 p 2 p 3 abc with all indices a, b, c summed. Then, the result can be written as a sum of four distinct terms, plus trivial permutations: The function F(a, b; z) in (108) is defined in terms of the Gaussian hypergeometric function, given in (140), and the function F 4 denotes the dressed Appell F 4 function, which is defined in (142). Finally, the function B p 1 p 2 p 3 ν 1 ν 2 |a is given by: Figure 9: The tree-level 4-point inflaton correlator with two massive exchanges.
derivatively coupled to two principal massive scalars in the bulk. The Lagrangian is given in (88). We use this example to show how to take folded limits in the three-vertex seed integral. The correlator corresponding to Fig. 9 has been given in (89), which shows that all we need to do is to take the following folded limit of the three-vertex seed integral (103): In the folded limit r 1 → 1 and r 4 → 1, we expect that various individual terms in (103) diverge, but the divergence must cancel out in the full result, as a consequence of choosing the Bunch-Davies initial condition. 7 Specifically, all the hypergeometric functions in the F factors in both (108) and (110) could develop divergent terms when we take their arguments to 1. With the knowledge that these divergent terms must cancel among themselves, we can directly throw them away when evaluating the function F(a, b; z) at z = 1. Using the expansion of hypergeometric function at argument unity, we get the following finite result for F(a, b; z) as z → 1 [91]: where Fin{} means the finite part of the expression within. The case of 1/2 − a − b ∈ Z can be computed by taking the limit. For example, the limit of a term in I FF when p → −2 can be computed as: With all limits of all F factors properly taken as above, we find a finite result for (89), which can again be separated into several pieces according their analytical properties at k 1 → 0 or k 2 → 0: Here r 2 ≡ k 1 /k 34 and r 3 ≡ k 2 /k 34 . The four pieces {T SS , T SB , T BS , T BB } are defined in a similar way as before, according to whether the expression is analytic in r 2 → 0 or r 3 → 0 limit. In particular, the signal-signal piece is nonanalytic when r 2 → 0 and when r 3 → 0, whose explicit expression is: where F 4 is again the dressed Appell F 4 function defined in (142). Next, the signal-background piece T SB is nonanalytic in r 2 → 0 but analytic in r 3 → 0: The background-signal piece T BS is obtained from T SB by switching ν 1 ↔ ν 2 as well as r 2 ↔ r 3 . Finally, the background-background piece is analytic in both r 2 → 0 and r 3 → 0 limit, and its expression is: Again, we have checked that our analytical result (125) for the 4-point graph in Fig. 9 agrees well with a direct numerical integration.

Conclusion and Outlooks
Inflation correlators are important theoretical data for QFTs in dS spacetime, and are promising targets for current and future cosmological observations. Inflation correlators mediated by massive fields are central objects in Cosmological Collider physics. Thus, the analytical computation of inflation correlators deserves a systematic investigation.
Very often in weakly coupled theories, inflation correlators are dominated by tree-level exchanges. However, analytical computation of general tree graphs remains challenging in dS, due to the multi-layer integrals with time orderings in Schwinger-Keldysh formalism. In previous works, it has been shown that the partial Mellin-Barnes representation is useful in the analytical computation of inflation correlators. Even using this method, the complete analytical evaluation of tree graphs is still hampered by the complication of nested time integrals.
In this work, we computed arbitrarily nested time integrals in PMB representation. The result is in general a multi-variable hypergeometric series. With our family-tree decomposition procedure, we can find series representation in terms of the inverse of any desired energy variable, or even the sum of several energy variables. This result largely solves the problem of analytical continuation of the nested time integrals in most physical regions.
With our results, the analytical computation of inflation correlators with arbitrary massive exchanges at the tree level is reduced to a work of pole collecting in the Mellin integrand, which is largely trivial. Thus, barring possible issues with analytical continuation in special kinematics to be commented below, we can say that the problem of analytical computation of tree-level inflation correlators is solved.
At this point, we want to comment on the meaning of analytical computation. As we have seen, most tree-level inflation correlators with massive exchanges have to be expressed in terms of gigantic hypergeometric series which are not yet named. One may say that we can systematically classify hypergeometric functions with increasing number of variables and parameters and give each of them a name. However, given that we know so little about these series in general, this is not super meaningful, and also is not really different from directly giving names to inflation correlators. The meaning of analytical calculation is thus somewhat obscure. Normally, when we say that we obtain an analytical answer, what we really mean is that we have good understanding of this answer in at least two ways: First, we know its analytical properties. This includes how does the answer change with parameters, and where does it blow up or show other singular behavior. Second, we know how to find numerical values of this answer for any choices of parameters with reasonable precision and computation time. Therefore, we can claim that we get an analytical answer only when we have gained sufficient knowledge about this answer. Having an answer is not enough; we need to understand it.
From this viewpoint, it seems to us that our result for the nested time integrals can be called an analytical answer: We know how to write down this answer as Taylor series for most kinematics. As long as there is a largest energy variable and we can use it to form small energy ratios, we can always express our answer as power series of these given small numbers. This means, on the one hand, we know the analytical properties of the answer at any soft energy limit, and also know how to take analytical continuations to different parameter space. On the other hand, having a convergent series often means that we can do fast numerical evaluation of the answer. This is proved true in our examples for two massive exchanges: In many cases, the numerical evaluation of our series solution is way faster than direct numerical integration of the original graph.
Our results have opened new possibilities in the analytical study of inflation correlators. Many interesting problems can be pursued along this direction. We mention some of them as below.
First, a main result of this work is a simple procedure to write down the analytical answer for arbitrary nested time integrals in PMB representation. Since the same time integral also appears in the loop computation, the result here could be useful for the computation of loop correlators as well. Thus it would be interesting to apply our method to make more complete loop computation in PMB representation.
Second, a nice feature of PMB representation is that the energy dependence and the momentum dependence of a graph are separated: The energy dependence is fully in the time integral, while the momentum dependence is fully in the loop momentum integral (or trivially factored out in a tree graph). Thus, our result on the nested time integral will be useful for studying energy dependence of a graph. This is particularly relevant to digging out the local CC signal in a graph, since, by definition, the local signal is a nonanalytic power in the energy ratios. We leave a more systematic study of local signals to a future work.
Third, it is important that the PMB representation does not assume full dS isometries of the problem. Therefore, it is straightforward to apply our results here to fields with dS-breaking dispersions such as non-unit sound speed, helical chemical potential, or even more exotic dispersion relations. Our method is also applicable to correlation functions in more general FRW background. We leave these generalizations to future studies.
Finally, it remains challenging to take analytical continuation of our series expressions to the parameter region where no small energy ratios exist. A pragmatic solution is using numerical interpolation to bridge different parameter regions in which various series expression converge. While this can indeed be implemented in some cases, it is not clear to us if this method works for all possible kinematics. Analytically, we may need more sophisticated methods to take analytical continuation for multi-layer series, which sounds like a nontrivial mathematical problem. We leave these more mathematically oriented problems for future studies as well.

A.1 Mellin-Barnes representation
We use MB representation for quite a few special functions in the main text, which we collect here. All expressions here can be found in standard mathematical handbooks such as [114].
First, the Hankel functions H (j) ν (az) of j'th kind (j = 1, 2) frequently appear. Their MB representations are given by: Next, we use an exponential integral E p (z) defined in the following way: This exponential integral is related to a confluent hypergeometric function U(a, b; z) via E p (z) = z p−1 e −z U(p, p; z), The confluent hypergeometric function U (a, b; z) has the following MB representations: The validities of these expressions put constraints on the range of z, which are always satisfied in the cases studied in this work, and thus we do not spell them out. With these expressions, we can get two different MB representations for the exponential integral E p (z). First, we have a partially resolved representation: Second, we have the following completely resolved representation: We note that the denominator 1/(s + p − 1) in the integrand of the last expression comes from the Γ factors Γ(b − 1 + s)/Γ(a + s) in (133) with a = b = p as required by (131). After taking a = b = p, most left poles of Γ(b − 1 + s) are canceled by the zeros of 1/Γ(a + s), with only one pole left, which is exactly the denominator 1/(s + p − 1) in (135). Thus we see that the pole from 1/(s + p − 1) should be treated as a left pole.
A number of hypergeometric series have been well studied and designated with special names. Several of these hypergeometric functions are used in the main text, and we collect their definitions here. More details about these functions can be found in [115]. First, the (generalized) hypergeometric function p F q is defined by the following way when the series converges: p F q a 1 , · · · , a p b 1 , · · · , b q z = ∞ n=0 (a 1 ) n · · · (a p ) n (b 1 ) n · · · (b q ) n z n n! , where (a) n ≡ Γ(a + n)/Γ(a) is the Pochhammer symbol. In most cases, it turns out simpler to use the following dressed version of hypergeometric function: p F q a 1 , · · · , a p b 1 , · · · , b q z = Γ a 1 , · · · , a p b 1 , · · · , b q p F q a 1 , · · · , a p b 1 , · · · , b q z = ∞ n=0 Γ a 1 + n, · · · , a p + n b 1 + n, · · · , b q + n z n n! .
A special case of Gauss hypergeometric function is frequently used in the main text, and thus we give a particular symbol to it: Next we come to hypergeometric functions of two variables. First, there are four Appell functions F 1 , · · · , F 4 . Two of them are used in the main text. We only present the definition of their dressed versions: Second, a more general class of two-variable hypergeometric functions are called Kampé de Fériet function in the literature, whose definition is: p+q F r+s a 1 , · · · , a p c 1 , · · · , c r b 1 , b ′ 1 ; · · · ; b q , b ′ q d 1 , d ′ 1 ; · · · ; d s , d ′ s x, y = ∞ m,n=0 (a 1 ) m+n · · · (a p ) m+n (c 1 ) m+n · · · (c r ) m+n x m y n m!n! .
Finally, there is a particular n-variable hypergeometric function that appears in the main text, called Lauricella's F A function: Again, we use this function in its dressed form:

B Details of Computing the Three-Vertex Seed Integral
In this appendix we collect some intermediate steps in the computation of the three-vertex seed integral (83) in Sec. 4. First, there are four independent branches in the seed integrals (83), shown below. The other four can be found by taking complex conjugation.
Then, we classify the terms in the above four integrals according to whether the adjacent two time variables are time-ordered (T) or factorized (F). Thus, we get the following eight different terms: I (TF) × D Similarly,