Massive Euclidean multiloop graphs#

integrate_graph accepts a connected scalar graph instead of a named topology. It constructs the Symanzik polynomials, checks UV subgraphs, applies Ramanujan’s Master Theorem to the common Schwinger scale, and numerically integrates the remaining parameters in Hepp sectors. This is an experimental hybrid evaluator: the remaining integration uses scrambled Sobol quadrature, not a multivariate Ramanujan/Method-of-Brackets residue sum.

For a multivariate residue series, the separate Multivariate Ramanujan series for the D=2 sunset backend now supports the unit-power massive \(D=2\) sunset in two explicitly justified convergence charts. This hybrid evaluator remains an independent reference.

Installation and first calculation#

pip install ".[multiloop]"    # from the repository; adds SciPy

from ramanujan_qft import ScalarEdge, ScalarGraph, integrate_graph

# Three parallel propagators between two vertices: a two-loop sunset.
# mass_squared and power may be integers, decimal strings, or rationals.
graph = ScalarGraph([
    ScalarEdge("a", "b", mass_squared=1),
    ScalarEdge("a", "b", mass_squared=1),
    ScalarEdge("a", "b", mass_squared=1),
])
result = integrate_graph(graph, D=2, epsrel=1e-4, seed=0)
print(graph.loops)
print(round(result.value, 4))
print(result.converged)

2
2.3439
True

result also contains error_estimate, standard_error, evaluations, sectors, samples_per_replicate, replicates, seed, and method.

For nonzero external momenta, specify vectors entering each vertex. Missing vertices have zero momentum, and the supplied components must sum to zero exactly. For example:

moving = ScalarGraph([
    ScalarEdge(0, 1, 1), ScalarEdge(0, 1, 4), ScalarEdge(0, 1, 9),
], external_momenta={0: [2, 1], 1: [-2, -1]})
parameters, U, F = moving.symanzik()
assert moving.loops == 2
x0, x1, x2 = parameters
assert U == x0*x1 + x0*x2 + x1*x2
# The external Euclidean invariant is p_E^2 = 2^2 + 1^2 = 5.

External vectors specify dot products independently of the real integration dimension D. Their component count does not constrain D; this permits analytic continuation in dimension inside the convergence domain. Decimal strings and rational numbers are stored exactly for graph algebra and power counting. Other numeric inputs are interpreted through their decimal string representation. Use explicit strings for reproducible decimal inputs.

Conventions and supported domain#

The returned integral is

\[ I_G=\int\prod_{\ell=1}^{L}\frac{d^Dk_\ell}{\pi^{D/2}} \prod_{e=1}^{E}(q_e^2+m_e^2)^{-\nu_e}. \]

All masses squared, propagator powers, and D must be finite and strictly positive real numbers. The graph must be connected; parallel edges, self-loops, bridges, and tree graphs are supported. No tensor numerator is included.

This Euclidean positive-propagator convention differs from the existing LoopTools-style functions. For integer powers and spacelike invariants, the corresponding Minkowski scalar integral in the package’s convention has an additional factor \((-1)^{\sum_e\nu_e}\). There is no automatic conversion or physical-region analytic continuation.

Every cyclic edge subset \(H\) is checked for

\[ \sum_{e\in H}\nu_e-\frac D2 L_H>0. \]

The check includes overlapping and disconnected subsets, and uses exact rational arithmetic. Testing only the whole graph would miss subdivergences. graph.uv_margins(D) exposes these margins with zero-based edge indices. Positive masses eliminate infrared singularities in this Euclidean domain. An integral failing a UV check raises ValueError; its analytic continuation or renormalized finite part is not substituted for the divergent integral.

Where Ramanujan’s theorem enters#

After Gaussian loop integration, write Schwinger parameters as \(\alpha_e=\rho x_e\). Define \(A=\sum_e\nu_e\) and \(s=A-LD/2\). The scale dependence is \(\rho^{s-1}\exp(-\rho F/U)\). Its Taylor coefficient continuation is \(\varphi(n)=(F/U)^n\), so RMT gives

\[ \int_0^\infty \rho^{s-1}e^{-\rho F/U}\,d\rho =\Gamma(s)(U/F)^s. \]

The code evaluates the equivalent unit-scale RMT factor with \(\varphi(n)=1\) and applies the scale power separately. This leaves the standard projective parameter integral

\[ I_G=\frac{\Gamma(s)}{\prod_e\Gamma(\nu_e)} \int_{x_e\geq0}\delta\left(1-\sum_e x_e\right) \prod_e x_e^{\nu_e-1}\frac{U^{s-D/2}}{F^s}\,d^Ex. \]

Graph construction follows spanning-tree and spanning-two-forest formulas: \(U=\sum_T\prod_{e\notin T}x_e\) and \(F=U\sum_e m_e^2x_e+\sum_{T_1\cup T_2}p_{T_1}^2\prod_{e\notin T_1\cup T_2}x_e\). See Bogner and Weinzierl, Feynman graph polynomials. The generalized RMT treatment of Feynman diagrams is described by González, Moll and Schmidt. The separate Multivariate Ramanujan series for the D=2 sunset implementation generates exact bracket families and justifies their summation for a restricted sunset domain.

Sector integration and errors#

For each ordering \(\pi\), projective invariance allows the largest parameter to be set to one. Set \(x_{\pi(0)}=1\) and \(x_{\pi(j)}=t_1\cdots t_j\). Common monomials are factored out of \(U\) and \(F\). Their exponents agree because every internal mass is positive. The remaining polynomials are bounded away from zero on the closed cube.

For the suffix \(H_j=\{\pi(j),\ldots,\pi(E-1)\}\), the endpoint weight is \(t_j^{\beta_j-1}\), with \(\beta_j=\sum_{e\in H_j}\nu_e-DL_{H_j}/2>0\). The map \(t_j=z_j^{\lceil\beta_j\rceil/\beta_j}\) replaces this with a nonnegative integer power of \(z_j\), including its Jacobian. This avoids introducing endpoint cusps when \(\beta_j>1\).

Each Sobol replicate estimates the sum over all sectors. Standard errors are computed across independently scrambled replicates, so sector correlations are included. Sample counts double until the diagnostic

max(3 * standard_error, change_since_previous_level, roundoff_floor)

is at most max(epsabs, epsrel * abs(value)). At least two levels are required. Single-edge integrals need no sampling and are evaluated directly.

These are statistical diagnostics, not rigorous error bounds or certified digits. Adaptive stopping does not provide a fixed confidence level. Extremely localized features or large hierarchies can fool sampling estimates; use independent methods and additional seeds for important results. The backend uses float64; dps is not an option, and epsrel < 1e-12 is refused. Common mass scales are factored out, but hierarchies exceeding the floating-point coefficient range are rejected. This does not limit the arbitrary-precision APIs elsewhere in the package.

If tolerance is not reached, MultiloopConvergenceError.result contains the partial estimate. For exploratory calculations, strict=False returns the same result with converged=False:

estimate = integrate_graph(moving, D=2, epsrel=1e-12,
                           min_power=3, max_power=4, replicates=4, strict=False)
assert not estimate.converged
assert estimate.evaluations == 6 * 4 * 16

Work is bounded by max_power, max_evaluations, max_sectors, and max_forests. UV enumeration is also capped at 65,535 edge subsets. There are \(E!\) sectors and \(2^E-1\) subsets: this initial implementation targets small graphs. Defaults permit at most 5,040 sectors (seven edges), but the default two-million-evaluation budget only permits the required initial refinement through six edges. Larger graphs need explicit budget increases; they can become impractical quickly. No scalability or speed superiority is claimed.

Validation and path to broader coverage#

tests/test_multiloop.py checks exact graph polynomials against a weighted Laplacian, numerical bubbles and a general triangle against independent simplex quadrature, the equal-mass two-loop vacuum sunset against a Clausen value, and the three-loop banana against \(7\zeta(3)\). One-dimensional two-, three-, and four-loop bananas and a kite are checked against independent position-space integrals. Other checks cover normalization, fractional powers and dimension, UV subdivergences, mass scaling, permutations, reproducibility, and resource exhaustion. These checks validate those examples, not every graph in the accepted domain.

The next steps toward a general multiloop package are:

  1. Extend the Multivariate Ramanujan series for the D=2 sunset charts across the uncovered Euclidean region, including equal-mass vacuum, then add topologies with independently derived convergence and pole-collision rules. Exact bracket algebra alone does not establish convergence.

  2. More efficient sector strategies and graph factorization to replace factorial enumeration; benchmark nontrivial five- to seven-edge graphs.

  3. Endpoint subtraction and Laurent expansion in \(D=D_0-2\epsilon\), including overlapping UV/IR divergences and a clearly defined renormalization scheme.

  4. Massless domains, causal continuation and contour deformation for physical thresholds, plus independent pySecDec/FIESTA reference data.

  5. Tensor numerators and integration-by-parts reduction, followed by a broader external validation campaign.

The legacy MoBDynamicEngine remains an incomplete bracket-analysis tool. Its unsupported evaluations still raise; this backend does not silently replace them with different conventions or a different method.