Phase-Space Integration and Monte Carlo Estimators

Phase-space integration turns a finite weighted matrix element into a number. For many final-state particles the domain is high-dimensional and the integrand is sharply structured by propagators, thresholds, soft and collinear remnants, cuts, and subtraction terms. Monte Carlo methods remain useful because their canonical $N^{-1/2}$ root-mean-square convergence rate does not worsen explicitly with dimension, although its variance prefactor certainly can. That advantage is meaningful only when the sampling density, variance, correlations, and convergence diagnostics are stated.

Required background. Lorentz-Invariant Phase Space supplies the measure and recursive factorizations. Measurement Functions and Inclusive Observables supplies the weight being integrated.

Helpful background. Probabilistic Convergence, Laws of Large Numbers, and Central Limit Theorems explains estimator convergence and its hypotheses. Product Measures, Fubini–Tonelli, and Change of Variables supplies the measure-theoretic basis for mappings and iterated phase space.

Importance sampling

Write the desired finite integral as

$$I=\int_\Omega \mathrm dx\,f(x).$$

Choose a normalized density $g(x)>0$ wherever $f(x)\ne0$, draw independent samples $X_k\sim g$, and define weights

$$w(X_k)=\frac{f(X_k)}{g(X_k)}.$$

The estimator

$$\widehat I_N=\frac1N\sum_{k=1}^N w(X_k)$$

is unbiased when $f$ is absolutely integrable. When the second moment is finite, its variance is

$$\operatorname{Var}(\widehat I_N) =\frac1N\left[ \int_\Omega \mathrm dx\,\frac{f(x)^2}{g(x)}-I^2 \right].$$

Thus $g$ should approximate $|f|$, not merely reproduce the phase-space volume. The unattainable zero-variance choice for a nonnegative integrand is $g=f/I$. Real calculations use analytic knowledge and adaptive estimates to approach this shape.

If samples are independent and the weight has finite variance, the central limit theorem motivates

$$\begin{aligned} \widehat I_N&\pm\frac{s_w}{\sqrt N},\\ s_w^2&=\frac{1}{N-1}\sum_k(w_k-\widehat I_N)^2. \end{aligned}$$

This one-standard-error estimate is a sampling diagnostic, not a theory uncertainty or an automatically calibrated confidence interval. Heavy tails, adaptive reuse of samples, quasi-random sequences, or correlated Markov samples require a different error analysis.

Mapping physical structures

A phase-space generator maps uniform random variables $u\in[0,1]^d$ to momenta obeying on-shell conditions and momentum conservation. The transformed integrand includes the Jacobian,

$$I=\int_{[0,1]^d}\mathrm du\, J(u)\,f(x(u)).$$

Good coordinates flatten known structures. For a resonance with invariant mass $s$, for example, a mapping based on

$$s=M^2+M\Gamma\tan y$$

absorbs the Breit–Wigner denominator into $ds/dy$ over the allowed range. A collinear channel benefits from variables adapted to an angle or transverse momentum; a threshold benefits from the appropriate velocity variable. The map must still cover the complete integration region with a known, nonzero density.

Recursive phase-space factorization gives

$$\begin{aligned} \mathrm d\Phi_n(P;p_1,\ldots,p_n) ={}&\frac{\mathrm ds_Q}{2\pi}\, \mathrm d\Phi_2(P;Q,p_n)\\ &\times\mathrm d\Phi_{n-1}(Q;p_1,\ldots,p_{n-1}), \end{aligned}$$

which naturally parameterizes sequential decays. It is a change of variables, not a narrow-width approximation: $s_Q$ remains integrated over its physical range unless an approximation is explicitly invoked.

Multi-channel sampling

One map rarely resolves all propagator and singular structures. Let normalized channel densities $g_c(x)$ cover complementary regions and choose nonnegative coefficients $\alpha_c$ with $\sum_c\alpha_c=1$. The mixture

$$g(x)=\sum_c\alpha_c g_c(x)$$

gives the weight

$$w(x)=\frac{f(x)}{\sum_c\alpha_cg_c(x)}.$$

Every sampled point must be divided by the full mixture density, even when it was generated through one channel. Otherwise overlapping channels are double counted. Channel weights can be adapted to reduce variance, but their update rule and the independence of the final error sample must be documented. Multi-channel optimization for scattering integrals is developed in Kleiss and Pittau 1994, pp. 141–146.

Adaptive stratification, such as VEGAS, learns a separable approximation to the marginal structure of the integrand Lepage 1978, pp. 192–203. It is effective when peaks align with coordinates, but a rotated or strongly correlated ridge still needs a physical mapping or multiple channels.

A two-body benchmark

For a decay of invariant mass $M$ to particles of masses $m_1,m_2$, integration over solid angle gives

$$\int \mathrm d\Phi_2 =\frac{\lambda^{1/2}(M^2,m_1^2,m_2^2)}{8\pi M^2},$$

with

$$\begin{aligned} \lambda(x,y,z)={}&x^2+y^2+z^2\\ &-2xy-2xz-2yz. \end{aligned}$$

For two massless daughters, the result is $1/(8\pi)$. A phase-space generator should reproduce this normalization for a constant integrand before it is trusted on a matrix element. Further checks include rotational invariance, the endpoint at $M=m_1+m_2$, and exact four-momentum conservation at every generated point.

A nonconstant moment also checks the estimator variance. Generate a uniform direction with $c=\cos\theta=2u-1$ and $\phi=2\pi v$, and weight the massless two-body phase space by

$$W(c)=1+c^2.$$

The weight is even under exchange of two identical final particles. Relative to the constant-integrand result $\Phi_2=1/(8\pi)$,

$$\frac{\Phi_2[W]}{\Phi_2} =\frac12\int_{-1}^{1}\mathrm dc\,(1+c^2) =\frac43.$$

For $N$ independent uniform samples, the ratio estimator is the mean of $1+c_k^2$. Its variance is fixed analytically:

$$\begin{aligned} \operatorname{Var}(1+c^2) &=\frac12\int_{-1}^{1}\mathrm dc\,(1+c^2)^2-\frac{16}{9}\\ &=\frac4{45}, \end{aligned}$$

so

$$\operatorname{SE}\!\left(\frac{\widehat\Phi_{2,N}[W]}{\Phi_2}\right) =\frac{2}{\sqrt{45N}}.$$

At $N=10^4$ the predicted standard error is $2.98\times10^{-3}$. This bounded fixture simultaneously checks the $4\pi$ angular Jacobian, the mean, the variance calculation, and the observed $N^{-1/2}$ scaling. It does not test a peaked matrix element or an adaptive proposal.

Convergence and failure diagnostics

Record more than the final standard error:

• estimates and variances for statistically independent batches;
• the largest absolute weights and their phase-space locations;
• effective sample sizes or channel occupancies;
• stability under alternative mappings and channel coefficients;
• scaling of the error with $N$ over a resolved range;
• cancellation loss when positive and negative subtraction contributions are combined;
• results for analytic constant or low-degree benchmark integrands;
• reproducible random seeds or a documented seed policy.

For signed integrands, $\int |f|$ can be much larger than $|I|$. Then a small absolute error on two large cancelling contributions may still be a large relative error on $I$. Integrating correlated pieces together can reduce variance, but the covariance must be included.

Unweighting converts positive weighted samples into equal-weight events by accepting with probability $w/w_{\max}$. It requires a valid upper bound and is inefficient for broad weight distributions; negative fixed-order weights do not admit this elementary interpretation. Event-generator construction is beyond this page.

No runnable phase-space laboratory is currently available. The analytic two-body normalization and batch diagnostics above provide the static reproducibility path available here.

Common pitfalls

Allowing $g=0$ where $f\ne0$. Importance sampling is then biased because part of the integral is never sampled. Every channel mixture must have complete support.

Quoting $1/\sqrt N$ without checking finite variance. Rare extreme weights can invalidate the observed error estimate. Inspect tails and repeat independent batches.

Reusing adaptive samples as if the density were fixed. Adaptation creates dependence between the learned map and sampled values. Separate training from the final integration or use a justified adaptive estimator.

Hiding normalization in the generator. State whether flux, symmetry factors, Jacobians, spin/color averages, and unit conversions are in $f$, $g$, or an outer factor.

Exercises

Use the ideal one-dimensional density

$$g(c)=\frac38(1+c^2), \qquad -1\le c\le1,$$

for the two-body moment above. Show that it is normalized and that its importance weight is constant.

Solution

Since $\int_{-1}^{1}(1+c^2)\,\mathrm dc=8/3$, the prefactor $3/8$ normalizes $g$. Relative to the uniform angular density $1/2$, the estimator for the moment ratio has weight

$$\frac{(1/2)(1+c^2)}{(3/8)(1+c^2)}=\frac43.$$

Every sample has the exact mean weight, so this idealized one-dimensional fixture has zero variance. Constructing or sampling the ideal density is not generally possible for a realistic unknown integrand.

Where to continue

• Fixed-Order Organization and Scale Dependence distinguishes integration error from perturbative truncation diagnostics.
• Validation and Theory Uncertainties records benchmark, convergence, and reproducibility checks together.
• Unstable-Particle Observables and Controlled Resonance Approximations explains when a resonance mapping is combined with a pole approximation and when it is not.

References

• Kleiss, Ronald, and Roberto Pittau. “Weight Optimization in Multichannel Monte Carlo.” Computer Physics Communications 83 (1994): 141–146. doi:10.1016/0010-4655(94)90043-4. Open PDF.
• Lepage, G. Peter. “A New Algorithm for Adaptive Multidimensional Integration.” Journal of Computational Physics 27 (1978): 192–203. doi:10.1016/0021-9991(78)90004-9.