# Inflation Correlators with Multiple Massive Exchanges

Zhong-Zhi Xianyu${}^{\,a}$111Email<EMAIL_ADDRESS>and Jiaju
Zang${}^{\,b}$222Email<EMAIL_ADDRESS>
_a Department of Physics, Tsinghua University, Beijing 100084, China_
_b Zhili College, Tsinghua University, Beijing 100084, China_

###### Abstract

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.

###### Contents

  1. 1 Introduction
  2. 2 Tree Graphs with Partial Mellin-Barnes
  3. 3 Time Integrals with Partial Mellin-Barnes
    1. 3.1 Family-tree decomposition of nested integrals
    2. 3.2 Partially ordered families: simple examples
    3. 3.3 General family integrals
    4. 3.4 Alternative representation
    5. 3.5 Discussions
  4. 4 General Two Massive Exchanges
    1. 4.1 Three-vertex seed integral
    2. 4.2 Computing the seed integral
  5. 5 Four-Point Correlator with Two Massive Exchanges
  6. 6 Conclusion and Outlooks
  7. A Mathematical Appendix
    1. A.1 Mellin-Barnes representation
    2. A.2 Special functions
  8. B Details of Computing the Three-Vertex Seed Integral

## 1 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 [2, 3, 4, 5, 6, 7, 8,
9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28,
29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47,
48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62]. 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 [63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73,
74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92,
93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108].
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 [109, 110,
111, 112, 19]. 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_{\mathsf{a}\mathsf{b}}(k;\tau_{1},\tau_{2})$ in SK formalism comes in four
types, depending on the four choices of SK indices $\mathsf{a},\mathsf{b}=\pm$
at the two endpoints. The two propagators with same-sign indices
$D_{\pm\pm}(k;\tau_{1},\tau_{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,
64, 65, 71, 70, 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 well-known 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_{\mathsf{a}\mathsf{b}}$ 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 $\delta$-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
$(\square_{x}-m^{2})D(x,y)=\mathrm{i}\delta(x-y)$. Then, we can make an EFT
expansion of
$D(x,y)\sim\mbox{$\frac{\,{\mathrm{i}}\,}{\square_{x}-m^{2}}$}\delta(x-y)$.
The leading order term is simple, which is just the contact graph with
$D(x,y)\sim-\mathrm{i}\delta(x-y)/m^{2}$. However, there are higher order
terms coming from acting on powers of $\square_{x}/m^{2}$ on $\delta(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 $\theta$-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:

$\displaystyle\mathbb{T}_{q_{1}\cdots
q_{V}}(E_{1},\cdots,E_{V})=\int\prod_{\ell=1}^{V}\Big{[}{\mathrm{d}}\tau_{\ell}(-\tau_{\ell})^{q_{\ell}-1}e^{\mathrm{i}E_{\ell}\tau_{\ell}}\Big{]}\prod_{i,j}\theta(\tau_{i}-\tau_{j}),$
(1)

where we have time integrals at $V$ vertices, nested arbitrarily by the
Heaviside $\theta$ 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 $\theta$-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
$\theta(\tau_{i}-\tau_{j})$, we call $\tau_{j}$ to be $\tau_{i}$’s mother and
$\tau_{i}$ to be $\tau_{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 $\mathcal{C}_{q_{1}\cdots
q_{N}}(E_{1},\cdots,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 $\mathcal{C}_{q_{1}\cdots
q_{N}}(E_{1},\cdots,E_{N})$. Working in the configurations where $E_{1}\gg
E_{i}$ with $i=2,\cdots,N$, we find:

$\displaystyle{\mathcal{C}}_{q_{1}\cdots q_{N}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},\cdots,E_{N})=\displaystyle\frac{\,{1}\,}{(\mathrm{i}E_{1})^{q_{1\cdots
N}}}\sum_{n_{2},\cdots,n_{N}=0}^{\infty}\Gamma(q_{1\cdots N}+n_{2\cdots
N})\prod_{j=2}^{N}\displaystyle\frac{\,{(-E_{j}/E_{1})^{n_{j}}}\,}{(\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern 2.0mu_{j}+\mkern
2.0mu\widetilde{\mkern-2.0mun\mkern-2.0mu}\mkern 2.0mu_{j})n_{j}!}.$ (2)

Here the hatted energy $\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{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\cdots N}\equiv q_{1}+\cdots+q_{N}$. The quantity $\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern 2.0mu_{j}$ means to take the
sum of all $q_{i}$ where either $i=j$ or $i$ is a descendent of $j$. $\mkern
2.0mu\widetilde{\mkern-2.0mun\mkern-2.0mu}\mkern 2.0mu_{j}$ is similarly
defined. Explicit application of this formula to the 5-layer graph in Fig. 1
is given in (25)-(3.1). 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 tree-level 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
$(\tau,\bm{x})$ where $\tau\in(-\infty,0)$ is the conformal time, and
$\bm{x}\in\mathbb{R}^{3}$ is the comoving coordinate. In this coordinate
system, the spacetime metric is
${\mathrm{d}}s^{2}=a^{2}(\tau)(-{\mathrm{d}}\tau^{2}+{\mathrm{d}}\bm{x}^{2})$,
where $a(\tau)=-1/(H\tau)$ 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 $\bm{k}$ to denote 3-momenta and the corresponding italic
letter $k\equiv|\bm{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 $\varrho_{ij}\equiv 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:

$\displaystyle k_{12}\equiv k_{1}+k_{2},$ $\displaystyle E_{ij}\equiv
E_{i}+E_{j},$ $\displaystyle s_{123}\equiv s_{1}+s_{2}+s_{3},$ $\displaystyle
q_{1\cdots N}\equiv q_{1}+\cdots+q_{N}.$ (3)

Finally, the Mellin integral measures are very often abbreviated in the
following way:

$\int_{s_{1},s_{2}\cdots}\equiv\int_{-\mathrm{i}\infty}^{+\mathrm{i}\infty}\displaystyle\frac{\,{{\mathrm{d}}s_{1}}\,}{2\pi\mathrm{i}}\displaystyle\frac{\,{{\mathrm{d}}s_{2}}\,}{2\pi\mathrm{i}}\cdots.$
(4)

## 2 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 $\varphi$ in the late time limit:

$\lim_{\tau\rightarrow
0}\big{\langle}\Omega\big{|}\varphi_{\bm{k}_{1}}(\tau)\cdots\varphi_{\bm{k}_{B}}(\tau)\big{|}\Omega\big{\rangle}_{\text{C}}=(2\pi)^{3}\delta^{(3)}(\bm{k}_{1}+\cdots+\bm{k}_{B})\mathcal{T}(\bm{k}_{1},\cdots,\bm{k}_{B}).$
(5)

As shown above, the correlation function is defined as an equal-time
expectation value of the product of $B$ operators $\varphi_{\bm{k}_{i}}$ in
3-momentum space, over a state $|\Omega\rangle$ which is taken asymptotic to
the Bunch-Davies vacuum state in the early time limit
$\tau\rightarrow-\infty$. We assume that the bulk theory of $\varphi$ is a
weakly coupled local quantum field theory. Therefore, after stripping off the
momentum-conserving $\delta$-function, the amplitude on the right hand side
$\mathcal{T}(\bm{k}_{1},\cdots,\bm{k}_{B})$ can be represented as an expansion
of connected graphs $\mathcal{G}(\bm{k}_{1},\cdots,\bm{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 $\varphi$, but we do assume that it
has a simple mode function $\varphi(k,\tau)$. More explicitly, if we expand
the mode $\varphi_{\bm{k}}$ in terms of canonically normalized creation and
annihilation operators $a_{\bm{k}}$ and $a_{-\bm{k}}^{\dagger}$, we get mode
function $\varphi(k,\tau)$ as the coefficient:

$\displaystyle\varphi_{\bm{k}}(\tau)=\varphi(k,\tau)a_{\bm{k}}+\varphi^{*}(k,\tau)a_{-\bm{k}}^{\dagger}.$
(6)

We suppress helicity indices if there are any. We assume that all the time
dependence of the mode function $\varphi(k,\tau)$ can be expressed as an
exponential factor $e^{-\mathrm{i}k\tau}$ times a polynomial of $-k\tau$. 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:

$\varphi(k,\tau)=\displaystyle\frac{\,{1}\,}{\sqrt{2k^{3}}}(1+\mathrm{i}k\tau)e^{-\mathrm{i}k\tau}.$
(7)

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
$\phi_{c}$ with mass $m_{c}=\sqrt{2}$ in $3+1$ dimensions, whose mode function
is:

$\phi_{c}(k,\tau)=\displaystyle\frac{\,{\tau}\,}{\sqrt{2k}}e^{-\mathrm{i}k\tau}.$
(8)

The bulk fields appearing in the tree graphs of
$\mathcal{T}(\bm{k}_{1},\cdots,\bm{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 $\varphi$. 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_ $\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu\equiv\sqrt{m^{2}-9/4}$. Then, according to the SK formalism [19], we can
construct four bulk propagator $D_{\mathsf{a}\mathsf{b}}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu)}(k;\tau_{1},\tau_{2})$ with $\mathsf{a},\mathsf{b}=\pm$ for such a
field. More explicitly:

$\displaystyle D_{-+}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu)}(k;\tau_{1},\tau_{2})=$
$\displaystyle~{}\displaystyle\frac{\,{\pi}\,}{4}e^{-\pi\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu}(\tau_{1}\tau_{2})^{3/2}\mathrm{H}_{\mathrm{i}\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu}^{(1)}(-k\tau_{1})\mathrm{H}_{-\mathrm{i}\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern 2.0mu}^{(2)}(-k\tau_{2}),$
(9) $\displaystyle D_{+-}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu)}(k;\tau_{1},\tau_{2})=$ $\displaystyle~{}\Big{[}D_{-+}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu)}(k;\tau_{1},\tau_{2})\Big{]}^{*},$ (10) $\displaystyle
D_{\pm\pm}^{(\mkern 2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu)}(k;\tau_{1},\tau_{2})=$ $\displaystyle~{}D_{\mp\pm}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu)}(k;\tau_{1},\tau_{2})\theta(\tau_{1}-\tau_{2})+D_{\pm\mp}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu)}(k;\tau_{1},\tau_{2})\theta(\tau_{2}-\tau_{1}).$ (11)

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:

$\displaystyle\mathcal{I}=\sum_{\mathsf{a}_{1},\cdots,\mathsf{a}_{V}=\pm}\int\prod_{\ell=1}^{V}\Big{[}{\mathrm{d}}\tau_{\ell}\,\mathrm{i}\mathsf{a}_{\ell}(-\tau_{\ell})^{p_{\ell}}e^{\mathrm{i}\mathsf{a}_{\ell}E_{\ell}\tau_{\ell}}\Big{]}\prod_{i=1}^{I}D_{\mathsf{a}_{i1}\mathsf{a}_{i2}}(K_{i},\tau_{i1},\tau_{i2}).$
(12)

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
$\tau_{\ell}$ ($\ell=1,\cdots,V$). Also, we introduce a factor of
$\mathrm{i}\mathsf{a}_{\ell}$ as required by the diagrammatic rule [19], and a
factor of $(-\tau_{\ell})^{p_{\ell}}$ to account for various types of
couplings as well as power factors in the external mode function (such as the
$\mathrm{i}k\tau$ term in the massless mode function (7)). The exponential
factor $e^{\mathrm{i}a_{\ell}E_{\ell}\tau_{\ell}}$ comes from the external
mode function, and $E_{\ell}$ 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 $\ell$. However, we note that $E_{\ell}$ is not the
total energy at Vertex $\ell$ since we do not include energies of bulk lines.
For each bulk line, we have a bulk propagator
$D_{\mathsf{a}_{i1}\mathsf{a}_{i2}}(K_{i},\tau_{i1},\tau_{i2})$ with momentum
$\bm{K}_{i}$, which is completely determined by external momenta via the
3-momentum conservation at each vertex. The two time variables
$\tau_{i1},\tau_{i2}$ as well as the two SK variables
$\mathsf{a}_{i1},\mathsf{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 $\theta$-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]:

$\displaystyle D_{\pm\mp}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu)}(k;\tau_{1},\tau_{2})=$
$\displaystyle~{}\displaystyle\frac{\,{1}\,}{4\pi}\int_{-\mathrm{i}\infty}^{+\mathrm{i}\infty}\displaystyle\frac{\,{{\mathrm{d}}s}\,}{2\pi\mathrm{i}}\displaystyle\frac{\,{{\mathrm{d}}\bar{s}}\,}{2\pi\mathrm{i}}\,e^{\mp\mathrm{i}\pi(s-\bar{s})}\Big{(}\displaystyle\frac{\,{k}\,}{2}\Big{)}^{-2(s+\bar{s})}(-\tau_{1})^{-2s+3/2}(-\tau_{2})^{-2\bar{s}+3/2}$
$\displaystyle\times\Gamma\Big{[}s-\displaystyle\frac{\,{\mathrm{i}\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu}\,}{2},s+\displaystyle\frac{\,{\mathrm{i}\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu}\,}{2},\bar{s}-\displaystyle\frac{\,{\mathrm{i}\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu}\,}{2},\bar{s}+\displaystyle\frac{\,{\mathrm{i}\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern 2.0mu}\,}{2}\Big{]},$ (13)

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 $\tau_{1}$ and $\bar{s}$ is associated with $\tau_{2}$.
The same-sign propagators $D_{\pm\pm}$ are obtained by substituting in the
above expression into (11). We note that the time-ordering $\theta$-functions
are left untransformed.

After taking the above PMB representation, the original SK integral (12)
becomes:

$\displaystyle\mathcal{I}=$
$\displaystyle\int_{-\mathrm{i}\infty}^{+\mathrm{i}\infty}\prod_{i=1}^{I}\bigg{\\{}\displaystyle\frac{\,{1}\,}{4\pi}\displaystyle\frac{\,{{\mathrm{d}}s_{i}}\,}{2\pi\mathrm{i}}\displaystyle\frac{\,{{\mathrm{d}}\bar{s}_{i}}\,}{2\pi\mathrm{i}}\Big{(}\displaystyle\frac{\,{K}\,}{2}\Big{)}^{-2s_{i\bar{i}}}\Gamma\Big{[}s_{i}-\displaystyle\frac{\,{\mathrm{i}\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu}\,}{2},s_{i}+\displaystyle\frac{\,{\mathrm{i}\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu}\,}{2},\bar{s}_{i}-\displaystyle\frac{\,{\mathrm{i}\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu}\,}{2},\bar{s}_{i}+\displaystyle\frac{\,{\mathrm{i}\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu}\,}{2}\Big{]}\bigg{\\}}$
$\displaystyle\times\bigg{\\{}\sum_{\mathsf{a}_{1},\cdots,\mathsf{a}_{V}=\pm}\int_{-\infty}^{0}\prod_{\ell=1}^{V}\Big{[}{\mathrm{d}}\tau_{\ell}\,\mathrm{i}\mathsf{a}_{\ell}(-\tau_{\ell})^{p_{\ell}-2\sum_{\ell}s}e^{\mathrm{i}\mathsf{a}_{\ell}E_{\ell}\tau_{\ell}}\Big{]}\mathcal{N}_{\mathsf{a}_{1}\cdots\mathsf{a}_{V}}\Big{(}\tau_{1},\cdots,\tau_{V};\\{s,\bar{s}\\}\Big{)}\bigg{\\}}.$
(14)

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 $(-\tau_{\ell})^{p_{\ell}-2\sum_{\ell}s}$ where $\sum_{\ell}s$
denotes the sum of all Mellin variables associated to $\tau_{\ell}$, 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
$\mathcal{N}_{\mathsf{a}_{1}\cdots\mathsf{a}_{V}}(\tau_{1},\cdots,\tau_{V};\\{s,\bar{s}\\})$
to represent all combinations of time-ordering $\theta$-functions, as well as
the SK-index-dependent phase factor $e^{\mp\mathrm{i}\pi(s-\bar{s})}$ in (2).

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 (2). 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]:

$\displaystyle\int_{-\infty}^{0}(-\tau_{1})^{q_{1}-1}(-\tau_{2})^{q_{2}-1}e^{\mathrm{i}E_{1}\tau_{1}+\mathrm{i}E_{2}\tau_{2}}\theta(\tau_{2}-\tau_{1})=\displaystyle\frac{\,{1}\,}{(\mathrm{i}E_{1})^{\mathrm{i}q_{12}}}\,{}_{2}\mathcal{F}_{1}\left[\begin{matrix}q_{2},q_{12}\\\
q_{2}+1\end{matrix}\middle|-\displaystyle\frac{\,{E_{2}}\,}{E_{1}}\right],$
(15)

where ${}_{2}\mathcal{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 (2), 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
$-\mathrm{i}\infty$ to $+\mathrm{i}\infty$. As a result, all _left poles_ are
contributed by the $\Gamma$-factors from the bulk propagators, shown in the
first line of (2). 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
$\Gamma$-factors in the first line of (2), we get:

$+2\sum_{i=1}^{I}(s_{i}+\bar{s}_{i})+\cdots.$ (16)

That is, all Mellin variables are summed together, with an overall coefficient
$+2$. Here “$\cdots$” 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 $\Gamma$-factors of the
form $\Gamma[\cdots-2\sum s]$. If we sum over the arguments of all right
poles, we will get:

$-2\sum_{i=1}^{I}(s_{i}+\bar{s}_{i})+\cdots,$ (17)

which is exactly the $s$ terms from the left-pole $\Gamma$-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\bar{i}}}$ in (2). 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
$\Gamma$-factors in (2). 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 (2). 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.333There is a degenerate situation where a vertex is _internal_ , in
the sense that it is not attached to any bulk-to-boundary propagator. In this
case, the energy variable $E=0$ at this vertex. Then, if we compute a
factorized time integral for this vertex alone, we get a $\delta$ function for
Mellin variables. Consequently, one should integrate out one Mellin variable
using this $\delta$ function instead of picking up poles. We will come back to
this point at the end of Sec. 3. 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.

## 3 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:

$\displaystyle\mathbb{T}_{q_{1}\cdots
q_{V}}(E_{1},\cdots,E_{V})=\int\prod_{\ell=1}^{V}\Big{[}{\mathrm{d}}\tau_{\ell}(-\tau_{\ell})^{q_{\ell}-1}e^{\mathrm{i}E_{\ell}\tau_{\ell}}\Big{]}\prod_{i,j}\theta(\tau_{i}-\tau_{j}).$
(18)

Here we are again considering a $V$-fold time integral with arbitrary nesting.
We require that all $\tau_{i}$ ($1\leq i\leq V$) appear in the
$\theta$-factors so that the integral is fully nested. Also, we have used a
factor $(-\tau_{\ell})^{q_{\ell}-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:

$q_{\ell}-1=p_{\ell}-2\sum\nolimits_{\ell}s.$ (19)

The difficulty with time ordering in (18) is easy to understand: A single time
integral of exponential with power factors from $\tau=-\infty$ to $\tau=0$
gives rise to a $\Gamma$ 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 $\Gamma$ functions after
finishing one layer of integral. Then we need to perform time integrals over
incomplete $\Gamma$ 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 $\Gamma$ 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
$\theta(\tau_{j}-\tau_{k})+\theta(\tau_{k}-\tau_{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 $\theta$-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
$\theta(\tau_{1}-\tau_{2})+\theta(\tau_{2}-\tau_{1})=1$ as:

(20)

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 $\tau_{i}$ is associated
with an energy variable $E_{i}$ and an exponent $q_{i}$, as is clear from
(18).

###  3.1 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 $\mathbb{T}_{q_{1}\cdots
q_{V}}(E_{1},\cdots,E_{V})$ is a function of $V$ energies $E_{1},\cdots,E_{V}$
and $V$ exponents $q_{1},\cdots,q_{V}$. 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
$\mathbb{T}_{q_{1}\cdots q_{V}}(E_{1},\cdots,E_{V})$, where there is a
_largest_ energy, say $E_{i}$, such that $E_{i}>E_{j}$ for all $j\neq i$. We
want to find an analytical expression for $\mathbb{T}_{q_{1}\cdots
q_{V}}(E_{1},\cdots,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 $\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{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
$\mathbb{T}_{q_{1}\cdots q_{V}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{1},E_{2},\cdots,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 $\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{2}$. The degenerate case
where there are multiple maximal energies will be considered in following
subsections.

Step 2:

We use the relation $\theta(\tau_{j}-\tau_{k})+\theta(\tau_{k}-\tau_{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 $\mathcal{C}_{q_{1}\cdots q_{N}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{1},E_{2},\cdots,E_{N})$
to denote a family integral with $N$ sites, where we have highlighted the
maximal energy $\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1}$ with a hat. Then, a family integral $\mathcal{C}_{q_{1}\cdots
q_{N}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},\cdots,E_{N})$ has the following form:

$\mathcal{C}_{q_{1}\cdots q_{N}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},\cdots,E_{N})=\int\prod_{\ell=1}^{N}\Big{[}{\mathrm{d}}\tau_{\ell}(-\tau_{\ell})^{q_{\ell}-1}e^{\mathrm{i}E_{\ell}\tau_{\ell}}\Big{]}\prod_{i,j}\theta(\tau_{i}-\tau_{j}),$
(21)

with the following restrictions on the $\theta$-function factors:

  1. 1.

Every time variable $\tau_{i}$ $(1\leq i\leq N)$ appears in time-ordering
$\theta$-function. (All sites belong to a family.)

  2. 2.

In a factor such as $\theta(\tau_{j}-\tau_{k})$, let us call $\tau_{j}$ to be
in the late position and $\tau_{k}$ in the early position. Then, it is
required that every variable $\tau_{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 $\tau_{i}$. (A mother can give birth to more
than one daughter.)

  3. 3.

The maximal energy site $\tau_{1}$ appears in $\theta$ 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
$\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{i}$, and it
becomes meaningless to ask for a series expansion in $\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{i}$ for those families.
We call them non-maximal energy families. Thus, for each of the non-maximal
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
$\mathbb{T}_{q_{1}\cdots q_{V}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{1},E_{2},\cdots,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:

$\displaystyle\mathcal{C}_{q_{1}\cdots q_{N}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},\cdots,E_{N})=\displaystyle\frac{\,{1}\,}{(\mathrm{i}E_{1})^{q_{1\cdots
N}}}\mkern 2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}\cdots q_{N}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{1},E_{2},\cdots,E_{N}),$
$\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern 2.0mu_{q_{1}\cdots
q_{N}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2}\cdots,E_{N})=\sum_{n_{2},\cdots,n_{N}=0}^{\infty}\Gamma(q_{1\cdots
N}+n_{2\cdots
N})\prod_{j=2}^{N}\displaystyle\frac{\,{(-\varrho_{j1})^{n_{j}}}\,}{(\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern 2.0mu_{j}+\mkern
2.0mu\widetilde{\mkern-2.0mun\mkern-2.0mu}\mkern 2.0mu_{j})n_{j}!}.$ (22)

Here, the hatted energy $\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1}$ represents the maximal energy. In the first line, we stripped away
a dimensionful factor $(\mathrm{i}E_{1})^{q_{1\cdots N}}$ so that the
resulting integral $\mkern 2.0mu\widetilde{\mkern-2.0muC\mkern-2.0mu}\mkern
2.0mu_{q_{1}\cdots q_{N}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{1},E_{2},\cdots,E_{N})$
is dimensionless. In the second line, we have defined $\varrho_{jk}\equiv
E_{j}/E_{k}$. Also, $\mkern 2.0mu\widetilde{\mkern-2.0mun\mkern-2.0mu}\mkern
2.0mu_{i}$ is defined to be the sum of $n_{i}$-variables over the site $i$ and
all its descendants. $\mkern 2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu$ 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:

$\displaystyle\mathbb{T}_{q_{1}\cdots
q_{5}}(E_{1},\cdots,E_{5})\equiv\int_{-\infty}^{0}\prod_{\ell=1}^{5}\Big{[}{\mathrm{d}}\tau_{\ell}(-\tau_{\ell})^{q_{\ell}-1}e^{\mathrm{i}E_{\ell}\tau_{\ell}}\Big{]}\theta(\tau_{2}-\tau_{1})\theta(\tau_{2}-\tau_{3})\theta(\tau_{4}-\tau_{3})\theta(\tau_{3}-\tau_{5}).$
(23)

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 $\theta(\tau_{i}-\tau_{j})+\theta(\tau_{j}-\tau_{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
$\mathcal{C}_{q_{3}q_{4}}(E_{3},E_{4})$ and
$\mathcal{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:

$\displaystyle\mathbb{T}_{q_{1}\cdots q_{5}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{1},\cdots,E_{5})=$
$\displaystyle~{}\mathcal{C}_{q_{1}q_{2}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2})\mathcal{C}_{q_{3}q_{4}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{3},E_{4})\mathcal{C}_{q_{5}}(E_{5})-\mathcal{C}_{q_{1}q_{2}q_{3}q_{4}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},E_{3},E_{4})\mathcal{C}_{q_{5}}(E_{5})$
$\displaystyle-\mathcal{C}_{q_{1}q_{2}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2})\mathcal{C}_{q_{4}q_{3}q_{5}}(E_{4},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{3},E_{5})+\mathcal{C}_{q_{1}q_{2}q_{3}q_{4}q_{5}}^{\text{(iso)}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},E_{3},E_{4},E_{5}).$ (24)

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.

Figure 1: The diagrammatic representation of (3.1), showing the reduction of a
5-layer time integral into partially ordered families. In this example, we
choose $E_{1}>E_{2},E_{3},E_{4},E_{5}$ and $E_{3}>E_{4},E_{5}$. 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.

Next, we assign $n_{2},\cdots,n_{4}$ for the sites with energy variables
$E_{2},\cdots,E_{4}$, respectively. It is clear from (3.1) that there are four
independent nontrivial families (i.e., with more than one site) involved in
this example. Applying the formula (3.1) for each of them, we get

$\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2})=\sum_{n_{2}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{2}}\Gamma(q_{12}+n_{2})}\,}{(q_{2}+n_{2})}\displaystyle\frac{\,{\varrho_{21}^{n_{2}}}\,}{n_{2}!},$
(25) $\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}q_{4}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},E_{3},E_{4})=\sum_{n_{2},n_{3},n_{4}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{234}}\Gamma(q_{1234}+n_{234})}\,}{(q_{234}+n_{234})(q_{34}+n_{34})(q_{4}+n_{4})}\displaystyle\frac{\,{\varrho_{21}^{n_{2}}}\,}{n_{2}!}\displaystyle\frac{\,{\varrho_{31}^{n_{3}}}\,}{n_{3}!}\displaystyle\frac{\,{\varrho_{41}^{n_{4}}}\,}{n_{4}!},$
(26) $\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{4}q_{3}q_{5}}(E_{4},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{3},E_{5})=\sum_{n_{4},n_{5}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{45}}\Gamma(q_{345}+n_{45})}\,}{(q_{4}+n_{4})(q_{5}+n_{5})}\displaystyle\frac{\,{\varrho_{43}^{n_{4}}}\,}{n_{4}!}\displaystyle\frac{\,{\varrho_{53}^{n_{5}}}\,}{n_{5}!},$
(27) $\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}q_{4}q_{5}}^{\text{(iso)}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},E_{3},E_{4},E_{5})$
$\displaystyle=\sum_{n_{2},n_{3},n_{4},n_{5}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{2345}}\Gamma(q_{12345}+n_{2345})}\,}{(q_{2345}+n_{2345})(q_{345}+n_{345})(q_{4}+n_{4})(q_{5}+n_{5})}\displaystyle\frac{\,{\varrho_{21}^{n_{2}}}\,}{n_{2}!}\displaystyle\frac{\,{\varrho_{31}^{n_{3}}}\,}{n_{3}!}\displaystyle\frac{\,{\varrho_{41}^{n_{4}}}\,}{n_{4}!}\displaystyle\frac{\,{\varrho_{51}^{n_{5}}}\,}{n_{5}!}.$
(28)

On the other hand, the result for the one-site family is trivial: $\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q}(E)=\Gamma(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
$\mathbb{T}_{q_{1}\cdots q_{5}}(E_{1},\cdots,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 $\mathbb{T}_{q_{1}\cdots q_{5}}(E_{1},\cdots,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.

$\displaystyle\mathbb{T}_{q_{1}\cdots q_{5}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{1},\cdots,E_{5})=$
$\displaystyle~{}\mathcal{C}_{q_{1}}(E_{1})\mathcal{C}_{q_{2}q_{3}q_{4}}(E_{2},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{3},E_{4})\mathcal{C}_{q_{5}}(E_{5})-\mathcal{C}_{q_{1}}(E_{1})\mathcal{C}_{q_{2}q_{4}q_{5}q_{3}}^{\text{(iso)}}(E_{2},E_{4},E_{5},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{3})$
$\displaystyle-\mathcal{C}_{q_{1}q_{2}q_{3}q_{4}}(E_{1},E_{2},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{3},E_{4})\mathcal{C}_{q_{5}}(E_{5})+\mathcal{C}_{q_{1}q_{2}q_{3}q_{4}q_{5}}^{\text{(iso)}}(E_{1},E_{2},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{3},E_{4},E_{5}).$ (29)

Figure 2: The diagrammatic representation of (3.1), showing the reduction of a
5-layer time integral into partially ordered families. In this example, we
choose $E_{3}>E_{i}$ $(i=1,2,4,5)$. The maximal energy site (Site 3) is marked
with a magenta circle.

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 (3.1):

$\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{2}q_{3}q_{4}}(E_{2},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{3},E_{4})=\sum_{n_{2},n_{4}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{24}}\Gamma(q_{234}+n_{24})}\,}{(q_{2}+n_{2})(q_{4}+n_{4})}\displaystyle\frac{\,{\varrho_{23}^{n_{2}}}\,}{n_{2}!}\displaystyle\frac{\,{\varrho_{43}^{n_{4}}}\,}{n_{4}!},$
(30) $\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{2}q_{4}q_{5}q_{3}}^{\text{(iso)}}(E_{2},E_{4},E_{5},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{3})=\sum_{n_{2},n_{4},n_{5}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{245}}\Gamma(q_{2345}+n_{245})}\,}{(q_{2}+n_{2})(q_{4}+n_{4})(q_{5}+n_{5})}\displaystyle\frac{\,{\varrho_{23}^{n_{2}}}\,}{n_{2}!}\displaystyle\frac{\,{\varrho_{43}^{n_{4}}}\,}{n_{4}!}\displaystyle\frac{\,{\varrho_{53}^{n_{5}}}\,}{n_{5}!},$
(31) $\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}q_{4}}(E_{1},E_{2},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{3},E_{4})=\sum_{n_{1},n_{2},n_{4}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{124}}\Gamma(q_{1234}+n_{124})}\,}{(q_{1}+n_{1})(q_{12}+n_{12})(q_{4}+n_{4})}\displaystyle\frac{\,{\varrho_{13}^{n_{1}}}\,}{n_{1}!}\displaystyle\frac{\,{\varrho_{23}^{n_{2}}}\,}{n_{2}!}\displaystyle\frac{\,{\varrho_{43}^{n_{4}}}\,}{n_{4}!},$
(32) $\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}q_{4}q_{4}}^{\text{(iso)}}(E_{1},E_{2},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{3},E_{4},E_{5})$
$\displaystyle=\sum_{n_{1},n_{2},n_{4},n_{5}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{1245}}\Gamma(q_{12345}+n_{1245})}\,}{(q_{1}+n_{1})(q_{12}+n_{12})(q_{4}+n_{4})(q_{5}+n_{5})}\displaystyle\frac{\,{\varrho_{13}^{n_{1}}}\,}{n_{1}!}\displaystyle\frac{\,{\varrho_{23}^{n_{2}}}\,}{n_{2}!}\displaystyle\frac{\,{\varrho_{43}^{n_{4}}}\,}{n_{4}!}\displaystyle\frac{\,{\varrho_{53}^{n_{5}}}\,}{n_{5}!}.$
(33)

Thus we have find an expression for the original 5-site integral
$\mathbb{T}_{q_{1}\cdots q_{5}}(E_{1},\cdots,E_{5})$ expanded as powers of
$1/E_{3}$. Let us emphasize that (3.1) and (3.1) are just different expansions
for the same function $\mathbb{T}_{q_{1}\cdots q_{5}}(E_{1},\cdots,E_{5})$,
with different validity regions.

###  3.2 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 (3.1). 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.

(a) $\mathcal{C}_{q}(E)$

(b) $\mathcal{C}_{q_{1}q_{2}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{1},E_{2})$

Figure 3: The one-site family and two-site family.

#### One-site family.

We begin with the simplest integral, the one-site family, shown in Fig. 3(a):

$\mkern 2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q}(E)={(\mathrm{i}E)^{q}}\int_{-\infty}^{0}{\mathrm{d}}\tau\,(-\tau)^{q-1}e^{\mathrm{i}E\tau}.$
(34)

The application of the rule is trivial, and we have the following answer:

$\mkern 2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q}(E)=\Gamma(q).$ (35)

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 $\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern 2.0mu_{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:

$\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2})=(\mathrm{i}E_{1})^{q_{12}}\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}{\mathrm{d}}\tau_{2}\,(-\tau_{1})^{q_{1}-1}(-\tau_{2})^{q_{2}-1}e^{\mathrm{i}(E_{1}\tau_{1}+E_{2}\tau_{2})}\theta(\tau_{2}-\tau_{1}).$
(36)

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:

$\displaystyle\mathcal{R}\Big{[}\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2})\Big{]}\equiv(\mathrm{i}E_{1})^{q_{12}}\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}{\mathrm{d}}\tau_{2}\,(-\tau_{1})^{q_{1}-1}(-\tau_{2})^{q_{2}-1}e^{\mathrm{i}(E_{1}\tau_{1}+E_{2}\tau_{2})}\theta(\tau_{1}-\tau_{2}).$
(37)

Then, the integral over $\tau_{2}$ can be performed, with the result expressed
in terms of an exponential integral $\mathrm{E}_{p}(z)$ whose definition is
given in (130):

$\displaystyle\mathcal{R}\Big{[}\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2})\Big{]}=(\mathrm{i}E_{1})^{q_{12}}\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}\,(-\tau_{1})^{q_{12}-1}e^{\mathrm{i}E_{1}\tau_{1}}\mathrm{E}_{1-q_{2}}(-\mathrm{i}E_{2}\tau_{1}).$
(38)

At this point, we make use of the following MB representation of
$\mathrm{E}_{p}(z)$:

$\mathrm{E}_{p}(z)=\int_{-\mathrm{i}\infty}^{+\mathrm{i}\infty}\displaystyle\frac{\,{{\mathrm{d}}s}\,}{2\pi\mathrm{i}}\displaystyle\frac{\,{\Gamma(s)z^{-s}}\,}{s+p-1}.$
(39)

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:

$\displaystyle\mathcal{R}\Big{[}\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2})\Big{]}=(\mathrm{i}E_{1})^{q_{12}}\int_{s_{2}}\displaystyle\frac{\,{\Gamma(s_{2})(\mathrm{i}E)^{-s_{2}}}\,}{s_{2}-q_{2}}\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}\,(-\tau_{1})^{q_{12}-1-s_{2}}e^{\mathrm{i}E_{1}\tau_{1}}.$
(40)

Then, the $\tau_{1}$ integral is trivial, which is simply given by
$\mathcal{C}_{q_{12}-s_{2}}(E_{1})=(\mathrm{i}E_{1})^{-q_{12}+s_{2}}\Gamma(q_{12}-s_{2})$.
So, finishing the $\tau_{1}$ integral, we get:

$\displaystyle\mathcal{R}\Big{[}\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2})\Big{]}=\int_{s_{2}}\displaystyle\frac{\,{\Gamma[s_{2},q_{12}-s_{2}]\varrho_{21}^{-s_{2}}}\,}{s_{2}-q_{2}}.$
(41)

Now it remains to finish the Mellin integral over $s_{2}$. Given that
$\varrho_{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,\cdots$, which is from the
$\Gamma$-factor $\Gamma(s_{2})$, and the other at $s_{2}=q_{2}$ coming from
the denominator. Collecting the residues at these poles, we get:

$\displaystyle\mathcal{R}\Big{[}\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2})\Big{]}=\sum_{n_{2}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{2}+1}\Gamma[n_{2}+q_{12}]}\,}{(n_{2}+q_{2})}\displaystyle\frac{\,{\varrho_{21}^{n_{2}}}\,}{n_{2}!}+\displaystyle\frac{\,{\Gamma[q_{1},q_{2}]}\,}{\varrho_{21}^{q_{2}}}.$
(42)

Now, we recognize that the last term without any summation is the product of
two one-site family:

$\displaystyle\frac{\,{1}\,}{(\mathrm{i}E_{1})^{q_{12}}}\displaystyle\frac{\,{\Gamma[q_{1},q_{2}]}\,}{\varrho_{21}^{q_{2}}}=\displaystyle\frac{\,{\Gamma(q_{1})}\,}{(\mathrm{i}E_{1})^{q_{1}}}\displaystyle\frac{\,{\Gamma(q_{2})}\,}{(\mathrm{i}E_{2})^{q_{2}}}=\mathcal{C}_{q_{1}}(E_{1})\mathcal{C}_{q_{2}}(E_{2}).$
(43)

Then, given the relation:

$\mathcal{R}\Big{[}{\mathcal{C}}_{q_{1}q_{2}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2})\Big{]}+{\mathcal{C}}_{q_{1}q_{2}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2})=\mathcal{C}_{q_{1}}(E_{1})\mathcal{C}_{q_{2}}(E_{2}),$ (44)

we see that the original family integral (36) is:

$\mkern 2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2})=\sum_{n_{2}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{2}}\Gamma[n_{2}+q_{12}]}\,}{(n_{2}+q_{2})}\displaystyle\frac{\,{\varrho_{21}^{n_{2}}}\,}{n_{2}!}.$
(45)

This is exactly what we would get using the rule (3.1). Incidentally, the
above summation can be directly done and the result is the well-known Gauss’s
hypergeometric function:

$\mkern 2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2})=\,{}_{2}\mathcal{F}_{1}\left[\begin{matrix}q_{2},q_{12}\\\
q_{2}+1\end{matrix}\middle|-\varrho_{21}\right].$ (46)

Here we use the dressed version ${}_{2}\mathcal{F}_{1}$ instead of the
original hypergeometric function ${}_{2}F_{1}$ for notational simplicity. The
dressed hypergeometric functions are defined in App. A.

Now, had we chosen to expand the integral (36) in terms of $1/E_{2}$, we will
get:

$\displaystyle\mathcal{C}_{2}(E_{1},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{2})=$
$\displaystyle~{}\displaystyle\frac{\,{\Gamma[q_{1},q_{2}]}\,}{(\mathrm{i}E_{1})^{q_{1}}(\mathrm{i}E_{2})^{q_{2}}}-\displaystyle\frac{\,{1}\,}{(\mathrm{i}E_{2})^{q_{12}}}\sum_{n_{1}}^{\infty}\displaystyle\frac{\,{\Gamma[n_{1}+q_{12}]}\,}{(n_{1}+q_{1})}\displaystyle\frac{\,{(-\varrho_{12})^{n_{1}}}\,}{n_{1}!}$
$\displaystyle=$
$\displaystyle~{}\displaystyle\frac{\,{\Gamma[q_{1},q_{2}]}\,}{(\mathrm{i}E_{1})^{q_{1}}(\mathrm{i}E_{2})^{q_{2}}}-\displaystyle\frac{\,{1}\,}{(\mathrm{i}E_{2})^{q_{12}}}{}_{2}\mathcal{F}_{1}\left[\begin{matrix}q_{1},q_{12}\\\
q_{1}+1\end{matrix}\middle|-\varrho_{12}\right].$ (47)

Clearly, the series expression in the first line in (3.2) 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}\rightarrow 0$ and the other at
$E_{1}/E_{2}\rightarrow 0$. This becomes more transparent after finishing the
summations of both series into hypergeometric functions. Indeed, equating (46)
with the second line of (3.2), we just get a transformation-of-variable
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.

(a) $\mathcal{C}_{q_{1}q_{2}q_{3}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{1},E_{2},E_{3})$

(b) $\mathcal{C}_{q_{1}q_{2}q_{3}}(E_{1},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{2},E_{3})$

Figure 4: Two independent family integrals at 3-site level.

#### 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:

$\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},E_{3})=(\mathrm{i}E_{1})^{q_{123}}\int_{-\infty}^{0}\prod_{i=1}^{3}\Big{[}{\mathrm{d}}\tau_{i}\,(-\tau_{i})^{q_{i}-1}e^{\mathrm{i}E_{i}\tau_{i}}\Big{]}\theta(\tau_{3}-\tau_{2})\theta(\tau_{2}-\tau_{1}).$
(48)

We again start from the completely reversed integral:

$\displaystyle\mathcal{R}\Big{[}\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},E_{3})\Big{]}=(\mathrm{i}E_{1})^{q_{123}}\int_{-\infty}^{0}\prod_{i=1}^{3}\Big{[}{\mathrm{d}}\tau_{i}\,(-\tau_{i})^{q_{i}-1}e^{\mathrm{i}E_{i}\tau_{i}}\Big{]}\theta(\tau_{1}-\tau_{2})\theta(\tau_{2}-\tau_{3}).$
(49)

Now, we can repeat the above strategy and finish the three layers of time
integrals in the order of $\tau_{3},\tau_{2},\tau_{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:

$\displaystyle~{}\mathcal{R}\Big{[}\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},E_{3})\Big{]}$ $\displaystyle=$
$\displaystyle~{}(\mathrm{i}E_{1})^{q_{123}}\int_{s_{3}}\displaystyle\frac{\,{\Gamma(s_{3})(\mathrm{i}E_{3})^{-s_{3}}}\,}{s_{3}-q_{3}}\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}{\mathrm{d}}\tau_{2}\,(-\tau_{1})^{q_{1}-1}(-\tau_{2})^{q_{23}-1-s_{3}}e^{\mathrm{i}(E_{1}\tau_{1}+E_{2}\tau_{2})}\theta(\tau_{1}-\tau_{2})$
$\displaystyle=$
$\displaystyle~{}(\mathrm{i}E_{1})^{q_{123}}\int_{s_{2},s_{3}}\displaystyle\frac{\,{\Gamma[s_{2},s_{3}](\mathrm{i}E_{2})^{-s_{2}}(\mathrm{i}E_{3})^{-s_{3}}}\,}{(s_{23}-q_{23})(s_{3}-q_{3})}\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}\,(-\tau_{1})^{q_{123}-1-s_{23}}e^{\mathrm{i}E_{1}\tau_{1}}$
$\displaystyle=$
$\displaystyle~{}\int_{s_{2},s_{3}}\displaystyle\frac{\,{\Gamma[s_{2},s_{3},q_{123}-s_{23}]\varrho_{21}^{-s_{2}}\varrho_{31}^{-s_{3}}}\,}{(s_{23}-q_{23})(s_{3}-q_{3})}.$
(50)

We start to observe a pattern here: Let the maximal-energy site be $\tau_{i}$.
When carrying out any but the last layer of time integrals, say $\tau_{j}$
with $j\neq i$, we are effectively generating a new layer of Mellin integral
with Mellin variable $s_{j}$, a pole-generating factor $\Gamma(s_{j})/(\mkern
2.0mu\widetilde{\mkern-2.0mus\mkern-2.0mu}\mkern 2.0mu_{j}-\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern 2.0mu_{j})$, and a power of
energy ratio $\varrho_{ji}^{-s_{j}}$. Here $\mkern
2.0mu\widetilde{\mkern-2.0mus\mkern-2.0mu}\mkern 2.0mu_{j}$ is the sum of all
Mellin variables assigned to the site $j$ and its descendant, and $\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern 2.0mu_{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
$\Gamma(s_{3})/(s_{3}-q_{3})$. Only after this is done, we then perform the
$s_{2}$-integral, by collecting poles from $\Gamma(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:

$\displaystyle\displaystyle\frac{\,{1}\,}{(\mathrm{i}E_{1})^{q_{123}}}\mathcal{R}\Big{[}\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},E_{3})\Big{]}=$
$\displaystyle~{}\displaystyle\frac{\,{1}\,}{(\mathrm{i}E_{1})^{q_{123}}}\sum_{n_{2},n_{3}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{23}}\Gamma[n_{23}+q_{123}]}\,}{(n_{23}+q_{23})(n_{3}+q_{3})}\displaystyle\frac{\,{\varrho_{21}^{n_{2}}}\,}{n_{2}!}\displaystyle\frac{\,{\varrho_{31}^{n_{3}}}\,}{n_{3}!}$
$\displaystyle-\displaystyle\frac{\,{\Gamma(q_{1})}\,}{(\mathrm{i}E_{1})^{q_{1}}}\displaystyle\frac{\,{1}\,}{(\mathrm{i}E_{2})^{q_{23}}}\sum_{n_{3}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{3}}\Gamma[n_{3}+q_{23}]}\,}{(n_{3}+q_{3})}\displaystyle\frac{\,{\varrho_{32}^{n_{3}}}\,}{n_{3}!}$
$\displaystyle-\displaystyle\frac{\,{\Gamma(q_{3})}\,}{(\mathrm{i}E_{3})^{q_{3}}}\displaystyle\frac{\,{1}\,}{(\mathrm{i}E_{1})^{q_{12}}}\sum_{n_{2}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{2}}\Gamma[n_{2}+q_{12}]}\,}{(n_{2}+q_{2})}\displaystyle\frac{\,{\varrho_{21}^{n_{2}}}\,}{n_{2}!}$
$\displaystyle+\displaystyle\frac{\,{\Gamma(q_{1})}\,}{(\mathrm{i}E_{1})^{q_{1}}}\displaystyle\frac{\,{\Gamma(q_{2})}\,}{(\mathrm{i}E_{2})^{q_{2}}}\displaystyle\frac{\,{\Gamma(q_{3})}\,}{(\mathrm{i}E_{3})^{q_{3}}}.$
(51)

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
$\theta(\tau_{i}-\tau_{j})=1-\theta(\tau_{j}-\tau_{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 (3.2). This is nothing but the original family integral. Thus:

$\mkern 2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},E_{3})=\sum_{n_{2},n_{3}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{23}}\Gamma[n_{23}+q_{123}]}\,}{(n_{23}+q_{23})(n_{3}+q_{3})}\displaystyle\frac{\,{\varrho_{21}^{n_{2}}}\,}{n_{2}!}\displaystyle\frac{\,{\varrho_{31}^{n_{3}}}\,}{n_{3}!}.$
(52)

Once again, this is exactly what we would get by applying the simple formula
(3.1). 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:

$\mkern 2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}}(\mkern 2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},E_{3})={\;}^{2+1}\mathcal{F}_{1+1}\left[\begin{matrix}q_{123},q_{23}\\\
q_{23}+1\end{matrix}\middle|\begin{matrix}\text{\;-\;},q_{3}\\\
\text{\;-\;},q_{3}+1\end{matrix}\middle|-\varrho_{21},-\varrho_{31}\right].$
(53)

The lesson we learned from the above example is the following: To find the
answer to a given family integral $\mathcal{C}$, all we need to do is to
compute another integral $\mathcal{R}[\mathcal{C}]$ with all time orderings
completely reversed. We compute $\mathcal{R}[\mathcal{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 $\Gamma$-factors. The result of this term is
automatically a sign factor $(-1)^{N-1}$ times the original family
$\mathcal{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:

$\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}}(E_{1},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{2},E_{3})=(\mathrm{i}E_{2})^{q_{123}}\int_{-\infty}^{0}\prod_{i=1}^{3}\Big{[}{\mathrm{d}}\tau_{i}\,(-\tau_{i})^{q_{i}-1}e^{\mathrm{i}E_{i}\tau_{i}}\Big{]}\theta(\tau_{1}-\tau_{2})\theta(\tau_{3}-\tau_{2}).$
(54)

The result after finishing all three layers of time integrals for the
_reversed diagram_ $\mathcal{R}[\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern 2.0mu]$ is:

$\displaystyle\mathcal{R}\Big{[}\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}}(E_{1},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{2},E_{3})\Big{]}=\int_{s_{1},s_{3}}\displaystyle\frac{\,{\Gamma[s_{1},s_{3},q_{123}-s_{13}]\varrho_{12}^{-s_{1}}\varrho_{32}^{-s_{3}}}\,}{(s_{1}-q_{1})(s_{3}-q_{3})}.$
(55)

Then, we finish the Mellin integral by picking up poles in
$\Gamma[s_{1},s_{3}]$ _only_. Multiplying the result by a trivial sign factor
of $(-1)^{3-1}=1$, we get the original family:

$\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}}(E_{1},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{2},E_{3})=\sum_{n_{1},n_{3}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{13}}\Gamma[n_{13}+q_{123}]}\,}{(n_{1}+q_{1})(n_{3}+q_{3})}\displaystyle\frac{\,{\varrho_{12}^{n_{1}}}\,}{n_{1}!}\displaystyle\frac{\,{\varrho_{32}^{n_{3}}}\,}{n_{3}!}.$
(56)

Again, it agrees with what we would get by applying (3.1). 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:

$\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}}(E_{1},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{2},E_{3})=\mathcal{F}_{2}\left[q_{123}\middle|\begin{matrix}q_{1},q_{3}\\\
1+q_{1},1+q_{3}\end{matrix}\middle|-\varrho_{12},-\varrho_{32}\right].$ (57)

The definition of $\mathcal{F}_{2}$ is collected in App. A.

(a)
$\mathcal{C}_{q_{1}q_{2}q_{3}q_{4}}^{\text{(cubic)}}(E_{1},E_{2},E_{3},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{4})$

(b) $\mathcal{C}_{q_{1}q_{2}q_{3}q_{4}}^{\text{(cubic)}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{1},E_{2},E_{3},E_{4})$

Figure 5: Two independent family integrals of 4-site graphs with a cubic
vertex.

#### 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
${\mathcal{C}}_{q_{1}q_{2}q_{3}q_{4}}^{\text{(n)}}$, 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 ${\mathcal{C}}_{q_{1}q_{2}q_{3}q_{4}}^{\text{(iso)}}$, as shown in Fig.
5.444We are borrowing the nomenclature of organic chemistry, where an
unbranched linear carbon chain is dubbed normal (n), while a chain with a
“cubic vertex” is dubbed isometric (iso). Again, there are two independent
ways to assign the maximal-energy site, either at the middle site or at a
boundary site, corresponding to Fig. 5(a) and Fig. 5(b), respectively.

Consider Fig. 5(a) first, where the maximal-energy site is chosen to be the
middle site $\tau_{4}$. The corresponding (dimensionless) family integral is:

$\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}q_{4}}^{\text{(iso)}}(E_{1},E_{2},E_{3},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{4})=(\mathrm{i}E_{4})^{q_{1234}}\int_{-\infty}^{0}\prod_{i=1}^{4}\Big{[}{\mathrm{d}}\tau_{i}\,(-\tau_{i})^{q_{i}-1}e^{\mathrm{i}E_{i}\tau_{i}}\Big{]}\theta(\tau_{1}-\tau_{4})\theta(\tau_{2}-\tau_{4})\theta(\tau_{3}-\tau_{4}).$
(58)

As always, we compute the corresponding reversed time-ordering integral:

$\displaystyle\mathcal{R}\Big{[}\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}q_{4}}^{\text{(iso)}}(E_{1},E_{2},E_{3},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{4})\Big{]}=\int_{s_{1},s_{2},s_{3}}\displaystyle\frac{\,{\Gamma[s_{1},s_{2},s_{3},q_{1234}-s_{123}]\varrho_{14}^{-s_{1}}\varrho_{24}^{-s_{2}}\varrho_{34}^{-s_{3}}}\,}{(s_{1}-q_{1})(s_{2}-q_{2})(s_{3}-q_{3})}.$
(59)

Then, the original integral is obtained by finishing the three-fold Mellin
integral in which we only collect poles from $\Gamma[s_{1},s_{2},s_{3}]$, and
multiplying the result by $(-1)^{4-1}=-1$. The result is:

$\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}q_{4}}^{\text{(iso)}}(E_{1},E_{2},E_{3},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{4})=\sum_{n_{1},n_{2},n_{3}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{123}}\Gamma[n_{123}+q_{1234}]}\,}{(n_{1}+q_{1})(n_{2}+q_{2})(n_{3}+q_{3})}\displaystyle\frac{\,{\varrho_{14}^{n_{1}}}\,}{n_{1}!}\displaystyle\frac{\,{\varrho_{24}^{n_{2}}}\,}{n_{2}!}\displaystyle\frac{\,{\varrho_{34}^{n_{3}}}\,}{n_{3}!}.$
(60)

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:

$\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}q_{4}}^{\text{(iso)}}(E_{1},E_{2},E_{3},\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{4})=\mathcal{F}_{A}\left[q_{1234}\middle|\begin{matrix}q_{1},q_{2},q_{3}\\\
q_{1}+1,q_{2}+1,q_{3}+1\end{matrix}\middle|-\varrho_{14},-\varrho_{24},-\varrho_{34}\right].$
(61)

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:

$\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}q_{4}}^{\text{(iso)}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},E_{3},E_{4})=(\mathrm{i}E_{1})^{q_{1234}}\int_{-\infty}^{0}\prod_{i=1}^{4}\Big{[}{\mathrm{d}}\tau_{i}\,(-\tau_{i})^{q_{i}-1}e^{\mathrm{i}E_{i}\tau_{i}}\Big{]}\theta(\tau_{2}-\tau_{4})\theta(\tau_{3}-\tau_{4})\theta(\tau_{4}-\tau_{1}).$
(62)

Without mentioning any details of intermediate steps, we directly provide the
final answer to this integral:

$\displaystyle\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}q_{4}}^{\text{(iso)}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},E_{3},E_{4})=\sum_{n_{2},n_{3},n_{4}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{234}}\Gamma[n_{234}+q_{1234}]}\,}{(n_{2}+q_{2})(n_{3}+q_{3})(n_{234}+q_{234})}\displaystyle\frac{\,{\varrho_{21}^{n_{2}}}\,}{n_{2}!}\displaystyle\frac{\,{\varrho_{31}^{n_{3}}}\,}{n_{3}!}\displaystyle\frac{\,{\varrho_{41}^{n_{4}}}\,}{n_{4}!}.$
(63)

###  3.3 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 (3.1) becomes a matter of mathematical induction. Below we
complete this proof.

We begin with a general partially ordered family with $N$ sites,
$\mathcal{C}_{q_{1}\cdots q_{N}}(\mkern
2.0mu\widehat{\mkern-2.0muE_{1}\mkern-2.0mu}\mkern 2.0mu,E_{2},\cdots,E_{N})$,
where we assume that the maximal energy site is $\tau_{1}$. Its integral
representation is given in (21). As in the previous section, we work with the
completely reversed integral $\mathcal{R}[\mathcal{C}_{q_{1}\cdots
q_{N}}(\mkern 2.0mu\widehat{\mkern-2.0muE_{1}\mkern-2.0mu}\mkern
2.0mu,E_{2},\cdots,E_{N})]$, and we integrate every time variable from
$-\infty$ to the time variable of her mother.

Suppose that we have finished all the time integrals for the descendants of
the site $\tau_{j}$, and now we want to finish the time integral at
$\tau_{j}$. Our induction assumption is that, after all the descendants of
$\tau_{j}$ integrated out, the integral over the variable $\tau_{j}$ has the
following form:

$\displaystyle\mathcal{I}_{j}(\tau_{M})=\int_{-\infty}^{\tau_{M}}{\mathrm{d}}\tau_{j}\,(-\tau_{j})^{\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu_{j}-1}e^{\mathrm{i}E_{j}\tau_{j}}\prod_{i\in\mathfrak{D}(j)}\int_{s_{i}}\displaystyle\frac{\,{\Gamma[s_{i}](\mathrm{i}E_{i})^{-s_{i}}}\,}{\mkern
2.0mu\widetilde{\mkern-2.0mus\mkern-2.0mu}\mkern 2.0mu_{i}-\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu_{i}}(-\tau_{j})^{-{s}_{i}},$ (64)

where $\tau_{M}$ is the time variable of $\tau_{j}$’s mother, and
$\mathfrak{D}(j)$ denotes the set of labels for all $\tau_{j}$’s descendants.
Now, finishing the $\tau_{j}$ integral, we get:

$\displaystyle\mathcal{I}_{j}(\tau_{M})=$
$\displaystyle~{}\prod_{i\in\mathfrak{D}(j)}\bigg{[}\int_{s_{i}}\displaystyle\frac{\,{\Gamma[s_{i}](\mathrm{i}E_{i})^{-s_{i}}}\,}{\mkern
2.0mu\widetilde{\mkern-2.0mus\mkern-2.0mu}\mkern 2.0mu_{i}-\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu_{i}}\bigg{]}(-\tau_{M})^{\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern 2.0mu_{j}-\sum
s_{i}}\mathrm{E}_{1-\mkern 2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu_{j}+\sum s_{i}}(-\mathrm{i}E_{j}\tau_{M})$ $\displaystyle=$
$\displaystyle~{}\prod_{i\in\mathfrak{D}(j)}\bigg{[}\int_{s_{i}}\displaystyle\frac{\,{\Gamma[s_{i}](\mathrm{i}E_{i})^{-s_{i}}}\,}{\mkern
2.0mu\widetilde{\mkern-2.0mus\mkern-2.0mu}\mkern 2.0mu_{i}-\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu_{i}}\bigg{]}\int_{s_{j}}\displaystyle\frac{\,{\Gamma(s_{j})(\mathrm{i}E_{j})^{-s_{j}}}\,}{s_{j}-\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu_{j}+\sum\limits_{i\in\mathfrak{D}(j)}s_{i}}(-\tau_{M})^{\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern 2.0mu_{j}-s_{j}-\sum s_{i}}$
$\displaystyle=$ $\displaystyle~{}(-\tau_{M})^{\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu_{j}}\prod_{i\in\\{j\\}\cap\mathfrak{D}(j)}\bigg{[}\int_{s_{i}}\displaystyle\frac{\,{\Gamma[s_{i}](\mathrm{i}E_{i})^{-s_{i}}}\,}{\mkern
2.0mu\widetilde{\mkern-2.0mus\mkern-2.0mu}\mkern 2.0mu_{i}-\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu_{i}}(-\tau_{M})^{-s_{i}}\bigg{]}.$ (65)

Here we have used the fact that
$s_{j}+\sum\limits_{i\in\mathfrak{D}(j)}s_{i}=\mkern
2.0mu\widetilde{\mkern-2.0mus\mkern-2.0mu}\mkern 2.0mu_{j}$. Also, we have
abbreviated $\sum\limits_{i\in\mathfrak{D}(j)}s_{i}$ as $\sum s_{i}$ when it
appears as an upper or lower index.

Now, to go one step further, we should finish the time integral over
$\tau_{M}$, which we denote as $\mathcal{I}_{M}(\tau_{G})$ and $\tau_{G}$ is
the time variable of $\tau_{j}$’s grandmother. To this end, we take products
of above integral $\mathcal{I}_{j}(\tau_{M})$ over all daughters $\tau_{j}$ of
$\tau_{M}$. Then, together with $\tau_{M}$’s own factor
$(-\tau_{M})^{q_{M}-1}e^{\mathrm{i}E_{M}\tau_{M}}$, we get:

$\displaystyle\mathcal{I}_{M}(\tau_{G})=$
$\displaystyle\int_{-\infty}^{\tau_{G}}{\mathrm{d}}\tau_{M}\,(-\tau_{M})^{q_{M}-1}e^{\mathrm{i}E_{M}\tau_{M}}\prod_{j\in\text{daughters
of $\tau_{M}$}}(-\tau_{M})^{\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu_{j}}\prod_{i\in\\{j\\}\cap\mathfrak{D}(j)}\bigg{[}\int_{s_{i}}\displaystyle\frac{\,{\Gamma[s_{i}](\mathrm{i}E_{i})^{-s_{i}}}\,}{\mkern
2.0mu\widetilde{\mkern-2.0mus\mkern-2.0mu}\mkern 2.0mu_{i}-\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu_{i}}(-\tau_{M})^{-s_{i}}\bigg{]}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\tau_{G}}{\mathrm{d}}\tau_{M}\,(-\tau_{M})^{\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu_{M}-1}e^{\mathrm{i}E_{M}\tau_{M}}\prod_{j\in\mathfrak{D}(M)}\bigg{[}\int_{s_{j}}\displaystyle\frac{\,{\Gamma[s_{j}](\mathrm{i}E_{j})^{-s_{j}}}\,}{\mkern
2.0mu\widetilde{\mkern-2.0mus\mkern-2.0mu}\mkern 2.0mu_{j}-\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu_{j}}(-\tau_{M})^{-s_{j}}\bigg{]}.$ (66)

This is identical to (64) upon a “generation shift” $j\rightarrow M$ and
$M\rightarrow 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 $\tau_{1}$ if we take $\tau_{j}=\tau_{1}$
and $\tau_{M}=0$. Then, completing this final layer of time integral over
$\tau_{1}$, we get, for the whole reversed family,

$\displaystyle\mathcal{R}\Big{[}\mathcal{C}_{q_{1}\cdots q_{N}}(\mkern
2.0mu\widehat{\mkern-2.0muE_{1}\mkern-2.0mu}\mkern
2.0mu,E_{2},\cdots,E_{N})\Big{]}=$
$\displaystyle\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}\,(-\tau_{1})^{\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu_{1}-1}e^{\mathrm{i}E_{1}\tau_{1}}\prod_{i=2}^{N}\int_{s_{i}}\displaystyle\frac{\,{\Gamma[s_{i}](\mathrm{i}E_{i})^{-s_{i}}}\,}{\mkern
2.0mu\widetilde{\mkern-2.0mus\mkern-2.0mu}\mkern 2.0mu_{i}-\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu_{i}}(-\tau_{1})^{-{s}_{i}}$ $\displaystyle=$
$\displaystyle~{}\prod_{i=2}^{N}\int_{s_{i}}\displaystyle\frac{\,{\Gamma[s_{i}](\mathrm{i}E_{i})^{-s_{i}}}\,}{\mkern
2.0mu\widetilde{\mkern-2.0mus\mkern-2.0mu}\mkern 2.0mu_{i}-\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu_{i}}\displaystyle\frac{\,{\Gamma[q_{1\cdots N}-s_{2\cdots
N}]}\,}{(\mathrm{i}E_{1})^{q_{1\cdots N}-s_{2\cdots N}}}$ $\displaystyle=$
$\displaystyle~{}\displaystyle\frac{\,{1}\,}{(\mathrm{i}E_{1})^{q_{1\cdots
N}}}\prod_{i=2}^{N}\bigg{[}\int_{s_{i}}\displaystyle\frac{\,{\Gamma[s_{i}]\varrho_{i1}^{-s_{i}}}\,}{\mkern
2.0mu\widetilde{\mkern-2.0mus\mkern-2.0mu}\mkern 2.0mu_{i}-\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern
2.0mu_{i}}\bigg{]}\Gamma[q_{1\cdots N}-s_{2\cdots N}].$ (67)

As shown many times in the previous subsection, the original family
$\mathcal{C}_{q_{1}\cdots q_{N}}(\mkern
2.0mu\widehat{\mkern-2.0muE_{1}\mkern-2.0mu}\mkern 2.0mu,E_{2},\cdots,E_{N})$
is recovered by picking up all poles from $\Gamma[s_{i}]$, and including an
overall factor $(-1)^{N-1}$ which comes from reversing the directions of $N-1$
bulk lines. Thus, we get:

$\displaystyle\mathcal{C}_{q_{1}\cdots q_{N}}(\mkern
2.0mu\widehat{\mkern-2.0muE_{1}\mkern-2.0mu}\mkern 2.0mu,E_{2},\cdots,E_{N})=$
$\displaystyle~{}\displaystyle\frac{\,{(-1)^{N-1}}\,}{(\mathrm{i}E_{1})^{q_{1\cdots
N}}}\sum_{n_{2},\cdots,n_{N}=0}^{\infty}\Gamma(q_{1\cdots N}+n_{2\cdots
N})\prod_{i=2}^{N}\displaystyle\frac{\,{(-\varrho_{i1})^{n_{i}}}\,}{(-\mkern
2.0mu\widetilde{\mkern-2.0mun\mkern-2.0mu}\mkern 2.0mu_{i}-\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern 2.0mu_{i})n_{i}!}.$ (68)

This is exactly the original family formula (3.1). Thus we have completed the
proof.

###  3.4 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 $\mathrm{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:

$\displaystyle\mathrm{E}_{p}(z)=e^{-z}\int_{-\mathrm{i}\infty}^{+\mathrm{i}\infty}\displaystyle\frac{\,{{\mathrm{d}}s}\,}{2\pi\mathrm{i}}\Gamma\begin{bmatrix}p+s,1+s,-s\\\
p\end{bmatrix}z^{-s-1}.$ (69)

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 $\mathrm{E}_{p}(z)$, as there is an exponential
factor $e^{-z}$ left.555Incidentally, this is very similar to the partially
resolved MB representation for the Whittaker function used in [89]. 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
$(\mathrm{i}E_{1})^{q_{123}}$, nor do we assign a maximal energy variable:

$\displaystyle{\mathcal{C}}_{q_{1}q_{2}q_{3}}({E}_{1},E_{2},E_{3})=\int_{-\infty}^{0}\prod_{i=1}^{3}\Big{[}{\mathrm{d}}\tau_{i}\,(-\tau_{i})^{q_{i}-1}e^{\mathrm{i}E_{i}\tau_{i}}\Big{]}\theta(\tau_{3}-\tau_{2})\theta(\tau_{2}-\tau_{1}).$
(70)

As before, we compute the integral with all time orderings reversed, but with
the new representation (69). The result is:

$\displaystyle\mathcal{R}\Big{[}{\mathcal{C}}_{q_{1}q_{2}q_{3}}(E_{1},E_{2},E_{3})\Big{]}=$
$\displaystyle\int_{s_{2},s_{3}}\,\Gamma\begin{bmatrix}1+s_{3},1-q_{3}+s_{3},-s_{3}\\\
1-q_{3}\end{bmatrix}\Gamma\begin{bmatrix}1+s_{2},2-q_{23}+s_{23},-s_{2}\\\
2-q_{23}+s_{3}\end{bmatrix}$
$\displaystyle\times(\mathrm{i}E_{3})^{-1-s_{3}}(\mathrm{i}E_{23})^{-1-s_{2}}\displaystyle\frac{\,{\Gamma[q_{123}-2-s_{23}]}\,}{(\mathrm{i}E_{123})^{q_{123}-2-s_{23}}}.$
(71)

Similar to the Mellin integrals in the previous representation, each Mellin
variable $s_{i}$ gets two sets of left poles from the $\Gamma$ factors, one at
$s_{i}=-1-n_{i}$ ($n_{i}=0,1,\cdots$) from $\Gamma(1+s_{i})$, and the other
more complicated, involving both $\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern 2.0mu_{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 $\mkern
2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}q_{2}q_{3}}(E_{1},E_{2},E_{3})$ is recovered by picking up poles
from all $\Gamma(1+s_{i})$ factors only and multiplying the result with an
appropriate sign factor. Thus:

$\displaystyle{\mathcal{C}}_{q_{1}q_{2}q_{3}}(E_{1},E_{2},E_{3})=\displaystyle\frac{\,{1}\,}{(\mathrm{i}E_{123})^{q_{123}}}\sum_{n_{2},n_{3}=0}^{\infty}(-1)^{n_{23}}\Gamma\begin{bmatrix}n_{23}-q_{123},-n_{23}-q_{23},-n_{3}-q_{3}\\\
1-n_{3}-q_{23},1-q_{3}\end{bmatrix}\varrho_{23T}^{n_{2}}\varrho_{3T}^{n_{3}},$
(72)

where we have defined $\varrho_{23T}\equiv E_{23}/E_{123}$ and
$\varrho_{3T}\equiv 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.

###  3.5 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 (3.1) in the following way:

$\mkern 2.0mu\widetilde{\mkern-2.0mu\mathcal{C}\mkern-2.0mu}\mkern
2.0mu_{q_{1}\cdots q_{N}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2}\cdots,E_{N})=\sum_{n_{2},\cdots,n_{N}=0}^{\infty}\Gamma(q_{1\cdots
N}+n_{2\cdots N})\prod_{j=2}^{N}\Gamma\begin{bmatrix}\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern 2.0mu_{j}+\mkern
2.0mu\widetilde{\mkern-2.0mun\mkern-2.0mu}\mkern 2.0mu_{j}\\\ \mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern 2.0mu_{j}+\mkern
2.0mu\widetilde{\mkern-2.0mun\mkern-2.0mu}\mkern
2.0mu_{j}+1\end{bmatrix}\displaystyle\frac{\,{(-\varrho_{j1})^{n_{j}}}\,}{n_{j}!}.$
(73)

Then it is clear that all the exponents $q_{i}$ $(i=1,\cdots,N)$ have the
positive coefficients when appearing in the arguments of $\Gamma$ 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 $\Gamma$ 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 (3.1), it is trivial to see that the
exponents $q_{i}$ ($i=1,\cdots,N$) appear in the $\Gamma$ factor
$\Gamma(q_{1\cdots N}+n_{2\cdots N})$ as a total sum. Then, let us look at
(19), which says that the value of $q_{\ell}$ at Vertex $\ell$ 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_{\ell}$ is equivalent to summing
over all Mellin variables. As a result, the argument of the $\Gamma$ factor
$\Gamma(q_{1\cdots N}+n_{2\cdots N})$ becomes:

$q_{1\cdots N}+n_{2\cdots N}=-2\sum_{i=1}^{I}(s_{i}+\bar{s}_{i})+p_{1\cdots
N}+n_{2\cdots N}+N,$ (74)

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 (3.1).
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 (3.1), only the leading term with $n_{2}=\cdots=n_{N}=0$
survives the limit. So we get:

$\displaystyle\lim_{E_{1}\rightarrow\infty}{\mathcal{C}}_{q_{1}\cdots
q_{N}}(E_{1},E_{2},\cdots,E_{N})=\displaystyle\frac{\,{\Gamma(q_{1\cdots
N})}\,}{(\mathrm{i}E_{1})^{q_{1\cdots
N}}}\prod_{j=2}^{N}\displaystyle\frac{\,{1}\,}{\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern 2.0mu_{j}}.$ (75)

Apart from the simple numerical factor $\prod\limits_{j}1/\mkern
2.0mu\widetilde{\mkern-2.0muq\mkern-2.0mu}\mkern 2.0mu_{j}$, this is very
similar to the result of one-site family
$\mathcal{C}_{q}(E)=\Gamma(q)/(\mathrm{i}E)^{q}$ with $E=E_{1\cdots N}$ and
$q=q_{1\cdots N}$. Here we have used the fact that $E_{1}\simeq E_{1\cdots N}$
in the $E_{1}\rightarrow\infty$ 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,\cdots,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}\rightarrow\infty$ limit can generate a power
$1/E_{1}^{q_{1}\cdots q_{N}}$ 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_{1}^{q_{1}\cdots q_{N}}$
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}\rightarrow\infty$, 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}\rightarrow 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^{\mathrm{i}E_{i}\tau_{i}}$ for each site at time $\tau_{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
$\varrho_{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
$\varrho_{41}\rightarrow 0$ limit, and we get:

$\displaystyle\mathcal{C}_{q_{1}q_{2}q_{3}q_{4}}^{\text{(iso)}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},E_{2},E_{3},E_{4}=0)=\sum_{n_{2},n_{3}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{23}}\Gamma[n_{23}+q_{1234}]}\,}{(n_{2}+q_{2})(n_{3}+q_{3})(n_{23}+q_{234})}\displaystyle\frac{\,{\varrho_{21}^{n_{2}}}\,}{n_{2}!}\displaystyle\frac{\,{\varrho_{31}^{n_{3}}}\,}{n_{3}!}.$
(76)

Figure 6: The family-tree decomposition of a 4-site nested time integral with
an internal vertex. The internal vertex (Vertex 4) with $E_{4}=0$ is marked
with an orange dot.

We take this opportunity to make a general comment on the computation of
graphs with internal vertices, as briefly mentioned in Footnote 3. We
illustrate the point with a concrete example. Suppose we want to compute the
following integral with $E_{4}=0$:

$\displaystyle\mathbb{T}_{q_{1}q_{2}q_{3}q_{4}}(E_{1},E_{2},E_{3})=\int_{-\infty}^{0}\prod_{i=1}^{4}\Big{[}{\mathrm{d}}\tau_{i}\,(-\tau_{i})^{q_{i}-1}e^{\mathrm{i}E_{i}\tau_{i}}\Big{]}\theta(\tau_{4}-\tau_{1})\theta(\tau_{4}-\tau_{2})\theta(\tau_{4}-\tau_{3}).$
(77)

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:

$\displaystyle\mathbb{T}_{q_{1}q_{2}q_{3}q_{4}}(E_{1},E_{2},E_{3})=$
$\displaystyle~{}\mathcal{C}_{q_{1}q_{4}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},0)\mathcal{C}_{q_{2}}(E_{2})\mathcal{C}_{q_{3}}(E_{3})-\mathcal{C}_{q_{1}q_{4}q_{2}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},0,E_{2})\mathcal{C}_{q_{3}}(E_{3})$
$\displaystyle-\mathcal{C}_{q_{1}q_{4}q_{3}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern
2.0mu_{1},0,E_{3})\mathcal{C}_{q_{2}}(E_{2})+\mathcal{C}_{q_{1}q_{4}q_{2}q_{3}}^{\text{(iso)}}(\mkern
2.0mu\widehat{\mkern-2.0muE\mkern-2.0mu}\mkern 2.0mu_{1},0,E_{2},E_{3}).$ (78)

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
(3.1) to all families here, and setting the summation variable $n_{4}=0$, we
get:

$\displaystyle\mathbb{T}_{q_{1}q_{2}q_{3}q_{4}}(E_{1},E_{2},E_{3})$
$\displaystyle=$
$\displaystyle~{}\displaystyle\frac{\,{\Gamma[q_{14},q_{2},q_{3}]/q_{4}}\,}{(\mathrm{i}E_{1})^{q_{14}}(\mathrm{i}E_{2})^{q_{2}}(\mathrm{i}E_{3})^{q_{3}}}-\bigg{\\{}\displaystyle\frac{\,{\Gamma[q_{3}]}\,}{(\mathrm{i}E_{1})^{q_{142}}(\mathrm{i}E_{3})^{q_{3}}}\sum_{n_{2}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{2}}\Gamma[q_{142}+n_{2}]}\,}{(n_{2}+q_{2})(n_{2}+q_{24})}\displaystyle\frac{\,{\varrho_{21}^{n_{2}}}\,}{n_{2}!}$
$\displaystyle+(2\leftrightarrow
3)\bigg{\\}}+\displaystyle\frac{\,{1}\,}{(\mathrm{i}E_{1})^{q_{1234}}}\sum_{n_{2},n_{3}=0}^{\infty}\displaystyle\frac{\,{(-1)^{n_{23}}\Gamma[n_{23}+q_{1234}]}\,}{(n_{2}+q_{2})(n_{3}+q_{3})(n_{23}+q_{234})}\displaystyle\frac{\,{\varrho_{21}^{n_{2}}}\,}{n_{2}!}\displaystyle\frac{\,{\varrho_{31}^{n_{3}}}\,}{n_{3}!}.$
(79)

On the other hand, we can as well compute the integral (77) directly. First,
we integrate out $\tau_{1}$, $\tau_{2}$, and $\tau_{3}$, and the result is:

$\displaystyle\mathbb{T}_{q_{1}q_{2}q_{3}q_{4}}(E_{1},E_{2},E_{3})=\prod_{i=1}^{3}\bigg{[}\int_{s_{i}}\displaystyle\frac{\,{\Gamma[s_{i}](\mathrm{i}E_{i})^{-s_{i}}}\,}{(s_{i}-q_{i})}\bigg{]}\int_{-\infty}^{0}{\mathrm{d}}\tau_{4}(-\tau_{4})^{q_{1234}-1-s_{123}}.$
(80)

Now, the final layer integral contains no energy variable since it is from an
internal vertex. Finishing this integral, we get a $\delta$ function:

$\int_{-\infty}^{0}{\mathrm{d}}\tau_{4}(-\tau_{4})^{q_{1234}-1-s_{123}}=(2\pi)\delta\big{[}\mathrm{i}(q_{1234}-s_{123})\big{]}.$
(81)

Now we need to choose a maximal energy. Suppose we choose $E_{1}$ without loss
of generality. Then, we use the above $\delta$ function to integrate out
$s_{1}$. Then, (80) becomes:

$\displaystyle\mathbb{T}_{q_{1}q_{2}q_{3}q_{4}}(E_{1},E_{2},E_{3})=\int_{s_{2},s_{3}}\displaystyle\frac{\,{\Gamma[q_{1234}-s_{23},s_{2},s_{3}](\mathrm{i}E_{1})^{s_{23}-q_{1234}}(\mathrm{i}E_{2})^{-s_{2}}(\mathrm{i}E_{3})^{-s_{3}}}\,}{(q_{234}-s_{23})(s_{2}-q_{2})(s_{3}-q_{3})}.$
(82)

We can finish this integral by collecting the residues of all left poles of
the integrand, as before. The result exactly agrees with (3.5).

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 $\delta$ 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 $\delta$ function from an internal vertex, we choose
a maximal energy among all energies connected at this vertex, and we use the
$\delta$ 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 (3.1), there is at least one
energy ratio $\varrho_{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 (3.1) in terms of a (generalized) hypergeometric
function ${}_{p}\mathcal{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\cdots 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 transformation-of-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.

## 4 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
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 $\mkern 2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{1,2}\equiv\sqrt{m_{1,2}^{2}-9/4}$ for the two lines, respectively.

Figure 7: A general tree graph with two massive exchanges.

###  4.1 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_ :

$\displaystyle\mathcal{I}_{\mathsf{a}_{1}\mathsf{a}_{2}\mathsf{a}_{3}}^{p_{1}p_{2}p_{3}}\equiv$
$\displaystyle-\mathrm{i}\mathsf{a}_{1}\mathsf{a}_{2}\mathsf{a}_{3}E_{1}^{p_{1}+1}E_{2}^{p_{2}+1}E_{3}^{p_{3}+1}\ell_{1}^{3}\ell_{2}^{3}\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}{\mathrm{d}}\tau_{2}{\mathrm{d}}\tau_{3}(-\tau_{1})^{p_{1}}(-\tau_{2})^{p_{2}}(-\tau_{3})^{p_{3}}$
$\displaystyle\times
e^{\mathrm{i}\mathsf{a}_{1}E_{1}\tau_{1}+\mathrm{i}\mathsf{a}_{2}E_{2}\tau_{2}+\mathrm{i}\mathsf{a}_{3}E_{3}\tau_{3}}D_{\mathsf{a}_{1}\mathsf{a}_{2}}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{1})}(\ell_{1};\tau_{1},\tau_{2})D_{\mathsf{a}_{2}\mathsf{a}_{3}}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{2})}(\ell_{2};\tau_{2},\tau_{3}).$ (83)

Here, as before, $E_{i}$ represents the total energies of _bulk-to-boundary_
lines at the Vertex $i$, while $\bm{\ell}_{1}$ and $\bm{\ell}_{2}$ represent
the 3-momenta of the two internal lines, respectively. As before, we have
included power factors of the form $(-\tau_{i})^{p_{i}}$ to allow for
different choices of external states and coupling types. Finally, the two bulk
massive propagators $D_{\mathsf{a}\mathsf{b}}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern 2.0mu_{1})}$ and
$D_{\mathsf{b}\mathsf{c}}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern 2.0mu_{2})}$ are given in
(9), (10), and (11). To minimize unnecessary complications, we will take
$p_{1},p_{2},p_{3}\in\mathbb{R}$. Generalization to complex values of $p_{i}$
is straightforward, although the expressions will be lengthier.

The $E$ and $\ell$ factors in front of the integral in (4.1) 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}\tau_{i}$ $(i=1,2,3)$, and use
the momentum ratios $r_{1}=\ell_{1}/E_{1}$, $r_{2}=\ell_{1}/E_{2}$,
$r_{3}=\ell_{2}/E_{2}$, $r_{4}=q_{\ell}/E_{3}$. Then, one can easily verify:

$\ell_{1}^{3}D_{\mathsf{a}\mathsf{b}}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{1})}(\ell_{1};\tau_{1},\tau_{2})=\mkern
2.0mu\widehat{\mkern-2.0muD\mkern-2.0mu}\mkern
2.0mu_{\mathsf{a}\mathsf{b}}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{1})}(r_{1}z_{1},r_{2}z_{2}),$ (84)

and a similar relation for the $\ell_{2}$-propagator. Then, the seed integral
is manifestly dimensionless and depends only on dimensionless energy ratios:

$\displaystyle\mathcal{I}_{\mathsf{a}_{1}\mathsf{a}_{2}\mathsf{a}_{3}}^{p_{1}p_{2}p_{3}}(r_{1},r_{2},r_{3},r_{4})\equiv$
$\displaystyle\int_{-\infty}^{0}\prod_{i=\ell}^{3}\Big{[}\mathrm{i}\mathsf{a}_{\ell}{\mathrm{d}}z_{\ell}(-z_{\ell})^{p_{\ell}}e^{\mathrm{i}\mathsf{a}_{\ell}z_{\ell}}\Big{]}\mkern
2.0mu\widehat{\mkern-2.0muD\mkern-2.0mu}\mkern
2.0mu_{\mathsf{a}_{1}\mathsf{a}_{2}}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{1})}(r_{1}z_{1},r_{2}z_{2})\mkern
2.0mu\widehat{\mkern-2.0muD\mkern-2.0mu}\mkern
2.0mu_{\mathsf{a}_{2}\mathsf{a}_{3}}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{2})}(r_{3}z_{2},r_{4}z_{3}).$ (85)

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 $\bm{k}_{1},\cdots,\bm{k}_{N}$.
Then, by definition,
$E_{1}=k_{1}+\cdots+k_{N}\geq|\bm{k}_{1}+\cdots+\bm{k}_{N}|=\ell_{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
$\mathcal{I}_{\mathsf{a}_{1}\mathsf{a}_{2}\mathsf{a}_{3}}^{p_{1}p_{2}p_{3}}$.
For example, when the external leges are conformal scalars $\phi_{c}$ with
cubic and quartic direct couplings with two massive scalars $\sigma_{1,2}$, we
can form a 6-point correlator. The Lagrangian is:

$\displaystyle{\mathscr{L}}\supset-\displaystyle\frac{\,{1}\,}{2}a^{2}(\partial_{\mu}\phi_{c})^{2}-\displaystyle\frac{\,{1}\,}{2}a^{4}m_{c}^{2}\phi_{c}^{2}-\displaystyle\frac{\,{1}\,}{2}\sum_{i=1}^{2}\Big{[}a^{2}(\partial_{\mu}\sigma_{i})^{2}+a^{4}m_{i}^{2}\sigma_{i}+a^{4}\mu_{i}\phi_{c}^{2}\sigma_{i}\Big{]}-\displaystyle\frac{\,{1}\,}{2}a^{4}\lambda\phi_{c}^{2}\sigma_{1}\sigma_{2}.$
(86)

Here $m_{c}=\sqrt{2}$ is the mass of the conformal scalar $\phi_{c}$, and
$m_{1,2}>3/2$ are the mass of two scalars $\sigma_{1,2}$, respectively. We
also include two cubic couplings with dimension-1 coupling constants $\mu_{i}$
$(i=1,2)$ and a quartic coupling with dimensionless coupling $\lambda$. The
powers of scale factors $a=-1/\tau$ are introduced to make the Lagrangian
scale invariant, and the spacetime indices in (86) are contracted by Minkowski
metric $\eta_{\mu\nu}$. 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:

$\displaystyle\mathcal{G}(\bm{k}_{1},\cdots,\bm{k}_{6})=-\mu_{1}\mu_{2}\lambda\displaystyle\frac{\,{(-\tau_{f})^{6}k_{12}k_{34}k_{56}}\,}{64k_{1}\cdots
k_{6}\ell_{1}^{3}\ell_{2}^{3}}\sum_{\mathsf{a}_{1},\mathsf{a}_{2},\mathsf{a}_{3}=\pm}\mathcal{I}_{\mathsf{a}_{1}\mathsf{a}_{2}\mathsf{a}_{3}}^{-2,-2,-2}\bigg{(}\displaystyle\frac{\,{\ell_{1}}\,}{k_{12}},\displaystyle\frac{\,{\ell_{1}}\,}{k_{56}},\displaystyle\frac{\,{\ell_{2}}\,}{k_{56}},\displaystyle\frac{\,{\ell_{2}}\,}{k_{34}}\bigg{)},$
(87)

where we have introduced a final-time cutoff $\tau_{f}$. We note that this
expression is for a single graph $\mathcal{G}(\bm{k}_{1},\cdots,\bm{k}_{6})$
rather than the whole correlator $\mathcal{T}(\bm{k}_{1},\cdots,\bm{k}_{6})$
at the same perturbative order. The correlator $\mathcal{T}$ can be obtained
from the graph $\mathcal{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 $\varphi$ with two massive exchanges, as shown in Fig. 9. We assume
that the inflaton $\varphi$ is coupled derivatively, to respect the
approximate shift symmetry of the inflaton field, and also to produce a
nontrivial result.666Directly coupled 2-point vertex (the mass mixing) is
trivial in the sense that it can be rotated away by diagonalizing the mass
matrix. The relevant Lagrangian is:

${\mathscr{L}}\supset-\displaystyle\frac{\,{1}\,}{2}a^{2}(\partial_{\mu}\varphi)^{2}-\sum_{i=1}^{2}\bigg{[}\displaystyle\frac{\,{1}\,}{2}a^{2}(\partial_{\mu}\sigma_{i})^{2}+\displaystyle\frac{\,{1}\,}{2}a^{4}m_{i}^{2}\sigma_{i}+a^{3}\mu_{i}\varphi^{\prime}\sigma_{i}\bigg{]}-\displaystyle\frac{\,{\lambda}\,}{2}a^{2}(\varphi^{\prime})^{2}\sigma_{1}\sigma_{2}.$
(88)

Then, the 4-point graph in Fig. 9 can be expressed as:

$\displaystyle\mathcal{G}(\bm{k}_{1},\cdots,\bm{k}_{4})=-\displaystyle\frac{\,{\mu_{1}\mu_{2}\lambda}\,}{16k_{1}^{3}k_{2}^{3}k_{3}k_{4}k_{34}}\sum_{\mathsf{a}_{1},\mathsf{a}_{2},\mathsf{a}_{3}=\pm}\mathcal{I}_{\mathsf{a}_{1}\mathsf{a}_{2}\mathsf{a}_{3}}^{-2,0,-2}\Big{(}1,\displaystyle\frac{\,{k_{1}}\,}{k_{34}},\displaystyle\frac{\,{k_{2}}\,}{k_{34}},1\Big{)}.$
(89)

Thus the computation of the 4-point correlator (89) requires us to take
simultaneous folded limit $r_{1}\rightarrow 1$ and $r_{4}\rightarrow 1$ which
is a bit nontrivial. We shall take this limit in the next section.

###  4.2 Computing the seed integral

Now we are going to compute the seed integral (4.1). The computation is rather
lengthy and tedious. Here we only outline the main steps, and collect more
details in App. B

The seed integral
$\mathcal{I}_{\mathsf{a}_{1}\mathsf{a}_{2}\mathsf{a}_{3}}^{p_{1}p_{2}p_{3}}(r_{1},r_{2},r_{3},r_{4})$
in (4.1) has 8 SK branches, depending on the values of the 3 SK indices
$\mathsf{a}_{1},\mathsf{a}_{2},\mathsf{a}_{3}=\pm$. We only need to compute 4
integrals: $\mathcal{I}_{+++}^{p_{1}p_{2}p_{3}}$,
$\mathcal{I}_{++-}^{p_{1}p_{2}p_{3}}$, $\mathcal{I}_{+-+}^{p_{1}p_{2}p_{3}}$,
and $\mathcal{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 $\tau_{1}$ and
$\tau_{2}$, or between $\tau_{3}$ and $\tau_{2}$, we should let $\tau_{2}$
take the earlier position. Thus, we use the following expression for the two
$D_{++}$ propagators:

$\displaystyle D_{++}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{1})}(\ell_{1};\tau_{1},\tau_{2})=$ $\displaystyle~{}D_{+-}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{1})}(\ell_{1};\tau_{1},\tau_{2})+\Big{(}D_{-+}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{1})}(\ell_{1};\tau_{1},\tau_{2})-D_{+-}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{1})}(\ell_{1};\tau_{1},\tau_{2})\Big{)}\theta(\tau_{1}-\tau_{2}),$ (90)
$\displaystyle D_{++}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{2})}(\ell_{2};\tau_{2},\tau_{3})=$ $\displaystyle~{}D_{-+}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{2})}(\ell_{2};\tau_{2},\tau_{3})+\Big{(}D_{+-}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{2})}(\ell_{2};\tau_{2},\tau_{3})-D_{-+}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{2})}(\ell_{2};\tau_{2},\tau_{3})\Big{)}\theta(\tau_{3}-\tau_{2}).$ (91)

After taking this representation, the expressions for the four integrals
$\mathcal{I}_{+++}^{p_{1}p_{2}p_{3}}$, $\mathcal{I}_{++-}^{p_{1}p_{2}p_{3}}$,
$\mathcal{I}_{+-+}^{p_{1}p_{2}p_{3}}$, and
$\mathcal{I}_{-++}^{p_{1}p_{2}p_{3}}$ can be obtained directly. Here we show
one example of $\mathcal{I}_{+++}^{p_{1}p_{2}p_{3}}$. The complete list is
given in (B) to (B).

$\displaystyle\mathcal{I}_{+++}^{p_{1}p_{2}p_{3}}=$
$\displaystyle-\mathrm{i}E_{1}^{p_{1}+1}E_{2}^{p_{2}+1}E_{3}^{p_{3}+1}\ell_{1}^{3}\ell_{2}^{3}\int_{-\infty}^{0}{\mathrm{d}}\tau_{1}{\mathrm{d}}\tau_{2}{\mathrm{d}}\tau_{3}(-\tau_{1})^{p_{1}}(-\tau_{2})^{p_{2}}(-\tau_{3})^{p_{3}}e^{\mathrm{i}(E_{1}\tau_{1}+E_{2}\tau_{2}+E_{3}\tau_{3})}$
$\displaystyle\times\Big{[}D_{+-}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{1})}(\ell_{1};\tau_{1},\tau_{2})+\Big{(}D_{-+}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{1})}(\ell_{1};\tau_{1},\tau_{2})-D_{+-}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{1})}(\ell_{1};\tau_{1},\tau_{2})\Big{)}\theta(\tau_{1}-\tau_{2})\Big{]}$
$\displaystyle\times\Big{[}D_{-+}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{2})}(\ell_{2};\tau_{2},\tau_{3})+\Big{(}D_{+-}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{2})}(\ell_{2};\tau_{2},\tau_{3})-D_{-+}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{2})}(\ell_{2};\tau_{2},\tau_{3})\Big{)}\theta(\tau_{3}-\tau_{2})\Big{]}.$
(92)

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).

$\displaystyle\mathcal{I}_{+++}^{p_{1}p_{2}p_{3}}=$
$\displaystyle~{}\mathcal{I}_{+++}^{\text{(FF)}}+\mathcal{I}_{+++}^{\text{(FT)}}+\mathcal{I}_{+++}^{\text{(TF)}}+\mathcal{I}_{+++}^{\text{(TT)}},$
(93) $\displaystyle\mathcal{I}_{++-}^{p_{1}p_{2}p_{3}}=$
$\displaystyle~{}\mathcal{I}_{++-}^{\text{(F)}}+\mathcal{I}_{++-}^{\text{(T)}},$
(94) $\displaystyle\mathcal{I}_{++-}^{p_{1}p_{2}p_{3}}=$
$\displaystyle~{}\mathcal{I}_{-++}^{\text{(F)}}+\mathcal{I}_{-++}^{\text{(T)}}.$
(95)

The explicit expressions of these integrals are given in (B)-(B). Then, we
define the following 4 integrals:

$\displaystyle\mathcal{I}^{\text{(FF)}}=$
$\displaystyle~{}\Big{[}\mathcal{I}_{+++}^{\text{(FF)}}+\mathcal{I}_{++-}^{\text{(F)}}+\mathcal{I}_{-++}^{\text{(F)}}+\mathcal{I}_{+-+}\Big{]}+\text{c.c.};$
(96) $\displaystyle\mathcal{I}^{\text{(FT)}}=$
$\displaystyle~{}\Big{[}\mathcal{I}_{+++}^{\text{(FT)}}+\mathcal{I}_{-++}^{\text{(T)}}\Big{]}+\text{c.c.};$
(97) $\displaystyle\mathcal{I}^{\text{(TF)}}=$
$\displaystyle~{}\Big{[}\mathcal{I}_{+++}^{\text{(TF)}}+\mathcal{I}_{++-}^{\text{(T)}}\Big{]}+\text{c.c.};$
(98) $\displaystyle\mathcal{I}^{\text{(TT)}}=$
$\displaystyle~{}\mathcal{I}_{+++}^{\text{(TT)}}+\text{c.c.}.$ (99)

After this regrouping of terms, each of the four integrals in
$\\{\mathcal{I}^{\text{(FF)}},\mathcal{I}^{\text{(FT)}},\mathcal{I}^{\text{(TF)}},\mathcal{I}^{\text{(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:

$\displaystyle D_{\pm\mp}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{1})}(\ell_{1};\tau_{1},\tau_{2})=$
$\displaystyle~{}\displaystyle\frac{\,{1}\,}{4\pi}\int_{s_{1},s_{2}}e^{\mp\mathrm{i}\pi(s_{1}-s_{2})}\Big{(}\displaystyle\frac{\,{\ell_{1}}\,}{2}\Big{)}^{-2s_{12}}(-\tau_{1})^{-2s_{1}+3/2}(-\tau_{2})^{-2s_{2}+3/2}$
$\displaystyle\times\Gamma\Big{[}s_{1}-\displaystyle\frac{\,{\mathrm{i}\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{1}}\,}{2},s_{1}+\displaystyle\frac{\,{\mathrm{i}\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{1}}\,}{2},s_{2}-\displaystyle\frac{\,{\mathrm{i}\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{1}}\,}{2},s_{2}+\displaystyle\frac{\,{\mathrm{i}\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern 2.0mu_{1}}\,}{2}\Big{]},$
(100) $\displaystyle D_{\pm\mp}^{(\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{2})}(\ell_{2};\tau_{2},\tau_{3})=$
$\displaystyle~{}\displaystyle\frac{\,{1}\,}{4\pi}\int_{s_{3},s_{4}}e^{\mp\mathrm{i}\pi(s_{3}-s_{4})}\Big{(}\displaystyle\frac{\,{\ell_{2}}\,}{2}\Big{)}^{-2s_{34}}(-\tau_{2})^{-2s_{3}+3/2}(-\tau_{3})^{-2s_{4}+3/2}$
$\displaystyle\times\Gamma\Big{[}s_{3}-\displaystyle\frac{\,{\mathrm{i}\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{2}}\,}{2},s_{3}+\displaystyle\frac{\,{\mathrm{i}\mkern
2.0mu\widetilde{\mkern-2.0mu\nu\mkern-2.0mu}\mkern
2.0mu_{2}}\,}{2},s_{4}-\displaystyle\frac{\,{\mathrm{i}\mkern