format: text

text:
# Two different quasiparticle scattering rates in vortex line liquid phase of layered 𝑑-wave superconductors

## Abstract

We carry out a quantum mechanical analysis of the behavior of nodal quasiparticles in the vortex line liquid phase of planar $`d`$-wave superconductors. Applying a novel path integral technique we calculate a number of experimentally relevant observables and demonstrate that in the low-field regime the quasiparticle scattering rates deduced from photoemission and thermal transport data can be markedly different from that extracted from tunneling, specific heat, superfluid stiffness or spin-lattice relaxation time.

In recent years, the physics of nodal quasiparticles in planar $`d`$-wave superconductors such as the high-$`T_c`$ cuprates has attracted a lot of both, theoretical and experimental, attention. The spectroscopic and transport properties of these Dirac-like quasiparticles with a linear dispersion have been studied in quite some detail, and elaborate analyses of various mechanisms of elastic scattering in the uniform superconducting state have been carried out, including the effects of potential , Kondo-like , and extended impurities  as well as twin boundaries .

In the mixed state, the recent quantum mechanical generalization of the earlier semiclassical approach  proposed in Ref. allows one to account for the non-uniformity of the local $`d`$-wave order parameter $`\mathrm{\Delta }_𝐩(𝐫)=\mathrm{\Delta }\mathrm{cos}(2\theta _𝐩)\mathrm{exp}(i\varphi (𝐫))`$ by means of a singular gauge transformation from the physical electrons $`c_\sigma (𝐫)`$ of spin $`\sigma `$ to the new fermionic quasiparticles. Unlike electrons, the latter are subject to the effective magnetic field with zero mean which, besides the physical field, also includes the supercurrent circulating outside vortex cores.

Applying the gauge transformation of Ref. to the electronic states with energies small compared to the maximum gap one can represent them in terms of the Nambu operators creating the auxiliary fermions with the momenta near the nodes of $`\mathrm{\Delta }_𝐩(𝐫)`$

$`\left(\begin{array}{c}c_\sigma (𝐫)\\ c_\sigma ^{}(𝐫)\end{array}\right)={\displaystyle \underset{n=1,2;\pm }{}}e^{\pm i𝐤_F^n𝐫}\left(\begin{array}{c}e^{i\varphi _A(𝐫)}u_{n\sigma }(𝐫)\\ e^{i\varphi _B(𝐫)}v_{n\sigma }(𝐫)\end{array}\right),`$ (1)
where $`n=1,2`$ labels the pairs of the opposite nodes, while the choice of the phases $`\varphi _{A,B}(𝐫)`$ is only restricted by the condition $`\varphi _A(𝐫)+\varphi _B(𝐫)=\varphi (𝐫)`$.

In the case of a regular vortex lattice, the representation (1) was used to demonstrate that the structure of the quasiparticle energy spectrum is that of the energy bands, rather than the Landau levels .

In the present paper, we extend the analysis based on the representation (1) to the experimentally well-documented vortex line liquid (VLL) phase  where the vortices are distributed totally randomly due to their strong pinning by columnar or other defects.

Any disorder, including that induced by random vortices, is expected to affect the behavior of the $`d`$-wave quasiparticles most strongly at the lowest energies, possibly resulting in the complete quasiparticle localization which, however, is still awaiting for experimental confirmation . In what follows, we focus on the ballistic regime of quasiparticle energies large compared to the localization scale (see below) which is readily accessible by a number of standard probes such as angular-resolved photoemission (ARPES), thermal transport, tunneling, specific heat, muon spin rotation ($`\mu `$SR), and spin-lattice relaxation.

In the quantum-mechanical approach of Ref., the combined effect of external magnetic field $`𝐇=\times 𝐀(𝐫)`$ and swirling supercurrent characterized by the superfluid velocity $`𝐯_s(𝐫)=\frac{\mathrm{}}{2}(\varphi _A+\varphi _B)\frac{e}{c}𝐀`$ does not amount solely to the semiclassical Doppler shift $`ϵϵ𝐤_F𝐯_s(𝐫)`$ of the quasiparticle energies . The latter is to be complemented by the vector potential, $`𝐚(𝐫)=\frac{\mathrm{}}{2}(\varphi _A\varphi _B)`$, that couples to the quasiparticles via their momentum, $`𝐤𝐤𝐚(𝐫)`$, and accounts for the quantum mechanical Berry phase corresponding to their Bohm-Aharonov (BA) scattering by the vortices. Thus, the complete Hamiltonian of the non-interacting nodal quasiparticles contains both the scalar- and the vector-like random terms

$$=\underset{n,\sigma }{}𝑑𝐫\overline{\psi }_\sigma [\widehat{\gamma }_0𝐤_F^n𝐯_s(𝐫)+v_i\widehat{\gamma }_i(p_ia_i(𝐫))]\psi _\sigma .$$
(2)
In Eq.(2), we used the $`4\times 4`$ representation for the matrices $`\widehat{\gamma }_\mu =(\widehat{\sigma }_2,i\widehat{\sigma }_1,i\widehat{\sigma }_3)\widehat{\sigma }_3`$ acting in the space of the Dirac bi-spinors composed of the Nambu spinors: $`\psi _\sigma =[(u_{1\sigma },ϵ_{\sigma \sigma ^{}}v_{1\sigma ^{}}),(u_{2\sigma },ϵ_{\sigma \sigma ^{}}v_{2\sigma ^{}})(\widehat{\sigma }_1+\widehat{\sigma }_3)/\sqrt{2}]`$.

In order to keep our discussion and formulas relatively simple we consider the case of isotropic quasiparticle dispersion and use the units where $`v_i=hc/2e=k_B=1`$. Moreover, because of the predominantly small-angle nature of the quasiparticle scattering by the vortices we choose to neglect the processes of inter-node scattering. While anticipating that neither of these simplifying assumptions will affect our main conclusions, we recognize that a combination of the above factors in the case of the real cuprates may give rise to an additional one-to-two dimensional crossover regime .

As in the previous studies of the VLL phase , we average over different vortex configurations by assuming the Gaussian distributions

$$\mathrm{g}_{(v,a)}^i(𝐪)\mathrm{g}_{(v,a)}^j(𝐪)=w_{(v,a)}(𝐪)\left(\delta _{ij}\frac{q_iq_j}{𝐪^2}\right),$$
(3)
for both $`𝐠_v=𝐯_s`$ and $`𝐠_a=𝐚`$. The striking difference between $`w_v(𝐪)=\alpha /(𝐪^2+\alpha )`$ and $`w_a(𝐪)=\alpha /𝐪^2`$ which are both proportional to the areal density of vortices $`\alpha =2\pi H`$ reflects the presence of screening for the scalar-like Doppler potential and its absence for the vector-like BA scattering. The longitudinal ($`q_iq_j/𝐪^2`$) part of the correlator (3) proves to contribute negligibly to all the quantities of interest (besides, it is suppressed by the Coulomb interactions).

In contrast to the previous analyses which focused solely on the effect of Doppler scattering , we find that vortex disorder has a profound effect on the quasiparticle spectrum which can not be adequately modelled by a constant quasiparticle width. In order to illustrate this point, we apply the standard self-consistent Born equation

$$\mathrm{\Sigma }(ϵ,𝐩)=\frac{d𝐪}{(2\pi )^2}\frac{ϵ+\mathrm{\Sigma }(ϵ,𝐪)}{𝐪^2(ϵ+\mathrm{\Sigma }(ϵ,𝐪))^2}w_{(v,a)}(𝐩+𝐪)$$
(4)
to the separate contributions of the two scattering mechanisms towards the total quasiparticle width $`\mathrm{Im}\mathrm{\Sigma }_v+\mathrm{Im}\mathrm{\Sigma }_a`$.

In the case of scalar disorder Eq.(4) yields $`\mathrm{Im}\mathrm{\Sigma }_v\alpha ^{1/2}`$ for small energies and momenta ($`ϵ,p\alpha ^{1/2}`$) while at $`\mathrm{max}(ϵ,p)\alpha ^{1/2}`$ it behaves as $`\alpha /\mathrm{max}(ϵ,p)`$ which dominates over the scattering by the vortex cores whose rate is estimated as $`\alpha /\mathrm{\Delta }`$.

At first sight, the effect of the BA scattering may seem to be much stronger, since a naive solution $`\mathrm{\Sigma }_a`$ of Eq.(4) with the singular kernel $`w_a(𝐪)`$ is plagued with a logarithmic infrared divergence of the momentum integral.

In order to avoid this spurious divergence which stems from the non-gauge invariant nature of the auxiliary fermion propagator one has to proceed directly with computing the manifestly gauge-invariant (retarded) Green function of the physical electrons

$$G^R(ϵ,𝐫)=_0^{\mathrm{}}𝑑te^{iϵt}c(t,𝐫)c^{}(0,\mathrm{𝟎})=\underset{n=1,2;\pm }{}e^{\pm i𝐤_F^n𝐫}$$

$$\psi _n(𝐫)\mathrm{exp}[i_C(𝐯_s+𝐚\sigma _3\mathrm{𝟏})𝑑𝐫^{}]\overline{\psi }_n(\mathrm{𝟎}).$$
(5)
It turns out that the exponential decay of this function (see below) makes it largely independent of the contour $`C`$ which can then be chosen as the straight path between the end points $`𝐫`$ and $`\mathrm{𝟎}`$.

To compute the amplitude (5) we apply the path-integral representation of Ref. to the propagator of the auxiliary Dirac fermions. First, for a fixed vortex configuration, we cast Eq.(5) in the form of a functional integral over the space-time coordinate $`𝐫(\tau )`$ and the conjugate momentum $`𝐩(\tau )`$ parameterized by the proper time $`\tau `$

$$𝒢^R(ϵ,𝐫|𝐯_s,𝐚)=_0^{\mathrm{}}𝑑\tau _{𝐫(0)=0}^{𝐫(\tau )=𝐫}D𝐫D𝐩e^{i\widehat{S}_0[\tau ]}$$

$$\mathrm{exp}[i_0^\tau d\tau ^{}(𝐤_F𝐯_s+\frac{d𝐫}{d\tau }𝐚)i_C(𝐯_s+𝐚)d𝐫^{})],$$
(6)
where we dropped, for the sake of compactness, the sum over the nodal points and introduced the free fermion action

$$\widehat{S}_0[\tau ]=_0^\tau 𝑑\tau ^{}[ϵ\widehat{\gamma }_0+𝐩(\frac{d𝐫}{d\tau ^{}}\widehat{\gamma })].$$
(7)
Averaging over the disorder variables $`𝐯_s`$ and $`𝐚`$ with the use of Eq.(3) results in the electron Green function $`G^R(ϵ,𝐫)=𝒢^R(ϵ,𝐫|𝐯_s,𝐚)`$ which is given by Eq.(6) where, instead of the exponential phase factor, the integrand contains a product of two attenuation factors

$$W_{(v,a)}[𝐫(\tau )]=\mathrm{exp}[\frac{1}{2}\frac{d𝐪}{(2\pi )^2}_0^\tau d\tau _1_0^\tau d\tau _2$$

$$𝐮_i(\tau _1)𝐮_i(\tau _2)w_{(v,a)}(𝐪)e^{i𝐪(𝐫(\tau _1)𝐫(\tau _2))}],$$
(8)
with $`𝐮_v=𝐯_F𝐫/\tau `$ and $`𝐮_a=d𝐫/d\tau 𝐫/\tau `$. Thus, the presence of the exponent of the line interal taken along the contour $`C`$ in Eq.(5) strongly reduces the effect of both $`𝐚`$ and $`𝐯_s`$, as compared to the case of the gauge-variant propagator of the auxiliary fermions. This observation seems to have been overlooked in the earlier studies of this and related problems where the phase factor in question would either not appear at all  or be averaged separately from the fermion propagator computed in a different approximation .

Proceeding along the lines of the previous analyses of the problem of non-relativistic fermions subject to a random vector potential , one can show that in the ballistic regime (which in the present case is defined by the condition $`ϵ\alpha ^{1/2}`$) the path integral (6) is dominated by the fermion trajectories which only slightly depart from the straight line $`𝐫_0(\tau )=𝐯\tau `$. Evaluating the factor $`W_v`$ for such a trajectory, one obtains

$$W_v[𝐫(\tau )]\mathrm{exp}\left[r\left(\sqrt{\frac{\tau }{r}}\sqrt{\frac{r}{\tau }}\right)^2_0^{\mathrm{}}\frac{dq}{2\pi }w_v(q)\right],$$
(9)
while the integral in $`W_a`$ turns out to be proportional to the so-called Amperian area of the closed contour composed of a fermion trajectory $`𝐫(\tau )`$ and the ”return” path $`𝐫_0(\tau )`$. Although this purely geometrical term vanishes for the saddle-point trajectory $`𝐫(\tau )=𝐫_0(\tau )`$, its expansion to first order in the transverse deviation $`r_{}(\tau )`$ yields

$$W_a[𝐫(\tau )]\mathrm{exp}[\frac{\alpha }{2}_0^\tau 𝑑\tau ^{}|r_{}(\tau ^{})|].$$
(10)
By analogy with the non-relativistic problem studied in Refs. the path integral (6) with the $`W_i`$ factors given by Eqs.(9) and (10) can be related to the resolvent of the Schroedinger equation describing the transverse motion of the Dirac fermion

$$[_x^2+(ϵ^2q^2)+(\alpha ϵ/q)^2(|x|+x_0))^2+i\alpha \mathrm{sign}x]$$

$$\times g(ϵ,q|x,x^{})=\delta (xx^{}),$$
(11)
where $`x_0=(\sqrt{|ϵ/q|}\sqrt{|q/ϵ|})^2/2\alpha ^{1/2}`$.

By analogy with the results of Refs. the averaged physical electron propagator $`G^R(ϵ,𝐩)`$ can be obtained by convoluting the kernel $`1/(𝐩^2q^2)^{3/2}`$ with the solution of (11) taken at $`x=x^{}=0`$ which is given by the formula

$$g(ϵ,q|0,0)=\left[\frac{d}{dx}\mathrm{ln}(U_+U_{})\right]_{x=0}^1$$
(12)
where $`U_\pm =U(a_\pm ,\lambda [x_0\pm x])`$ is the parabolic cylinder function of the parameter $`a_\pm =(ϵ^2q^2\pm i\alpha )/\lambda ^2`$ and $`\lambda =(2i\alpha |ϵ/q|)^{1/2}`$.

Turning now to the applications of Eq.(5), we first discuss the electron spectral function satisfying the dispersion relation $`G^R(ϵ,𝐩)=A(ϵ^{},𝐩)𝑑ϵ^{}/\pi (ϵϵ^{}+i\delta )`$. Near the maximum, $`|ϵ^2𝐩^2|\alpha `$, it can be shown to take the form

$$A(ϵ,𝐩)\frac{(ϵ\widehat{\gamma }_0𝐩\widehat{\gamma })\sqrt{\beta }}{[(ϵ^2p^2)^2+\beta ^2(1+(|ϵ/p|1)^4/4)]^{3/4}}$$
(13)
where $`\beta \pi \alpha /8`$, thereby demonstrating a replacement of the bare pole by a branch cut of the function $`(zz_0)^{3/2}`$ resulting from the above convolution procedure. According to Eq.(13), the decay of the electron propagator in real space ($`G^R(ϵ,𝐫)ϵ(\mathrm{\Gamma }/r)^{1/2}e^{\mathrm{\Gamma }r}`$ for $`r1/\mathrm{\Gamma }`$) is governed by $`\mathrm{\Gamma }(ϵ)\alpha /ϵ`$ which should be thought of as the actual (energy-dependent) quasiparticle width.

This direct experimenta prediction can be tested by performing ARPES measurements in the VLL phase of the cuprates under the weak-field conditions ($`\sqrt{H}T\mathrm{\Delta }`$).

Also, comparing Eq.(13) with the estimate for $`\mathrm{Im}\mathrm{\Sigma }_v`$ obtained from Eq.(4), we conclude that in the ballistic regime both the BA and the Doppler scattering mechanisms appear to be equally important, contrary to the conclusions drawn in Ref..

Next, we compute thermal conductivity given by the averaged product of two electron propagators

$$\kappa _{xx}=_0^{\mathrm{}}\frac{(ϵ/T)^2dϵ}{\mathrm{cosh}^2(ϵ/2T)}\frac{d𝐫}{2\pi }\mathrm{Tr}\widehat{\gamma }_1𝒢^A(ϵ,𝐫)\widehat{\gamma }_1𝒢^R(ϵ,𝐫).$$
(14)
The corresponding path integral reads as

$$\widehat{\gamma }𝒢^A(ϵ,𝐫)\widehat{\gamma }𝒢^R(ϵ,𝐫)=_0^{\mathrm{}}𝑑\tau _1𝑑\tau _2_{𝐫_{1,2}(0)=0}^{𝐫_{1,2}(\tau _{1,2})=\pm 𝐫}$$

$$\underset{\alpha ,\beta =1,2}{}D𝐫_\alpha D𝐩_\beta \widehat{\gamma }e^{i\widehat{S}_0[\tau _1]}\widehat{\gamma }e^{i\widehat{S}_0[\tau _2]}\underset{i=a,v}{}W_i(𝐫_\alpha 𝐫_\beta ),$$
(15)
where the factors $`W_i(𝐫_\alpha 𝐫_\beta )`$ with $`\alpha \beta `$ account for the vertex corrections, alongside the self-energy ones ($`\alpha =\beta `$). Thus, the path-integral method of computing the Dirac fermion conductivity is capable of proceeding beyond the conventional (non-crossing and fan-shaped) ladder series of the vertex corrections to the bare fermion bubble (the latter suffice only if the number $`N`$ of the Dirac species is large, while in the $`d`$-wave problem $`N=_\sigma 1=2`$).

Upon integrating over the ”center of mass” variables $`𝐫_1+𝐫_2`$ and $`𝐩_1+𝐩_2`$ and rescaling the ones describing relative motion, we again arrive at Eq.(11). This time around, it is formulated in terms of the transverse relative coordinate $`(r_1r_2)`$, and its solution yields the average $`𝒢^A𝒢^Rϵ(\mathrm{\Gamma }_{tr}/r)^{1/2}e^{\mathrm{\Gamma }_{tr}r}`$, where $`\mathrm{\Gamma }_{tr}=2\mathrm{\Gamma }`$ is now playing the role of the momentum relaxation rate, in agreement with the above estimate for the quasiparticle spectral width. Plugging this asymptote into Eq.(14) we find that in the low-field regime the thermal conductivity behaves as

$$\kappa _{xx}(T,H)\frac{7\pi ^2}{15}\frac{T^3}{H}.$$
(16)
While being in agreement with the estimate based on the kinetic equation $`\kappa T\nu /\mathrm{\Gamma }_{tr}`$ proportional to the linear density of states (DOS) $`\nu T`$ and $`\mathrm{\Gamma }_{tr}H/T`$, Eq.(16) is strikingly different from the result ($`\kappa T^2/H^{1/2}`$) that one would obtain by naively assuming that the quasiparticle width remains constant ($`\mathrm{\Gamma }^{}H^{1/2}`$) up to the energies $`ϵT`$ (cf. with Refs.).

The analysis of the data of Ref. taken in $`YBa_2Cu_3O_{6.99}`$ shows that Eq.(16) should be expected to hold for $`T<30K`$ and $`0.1\sqrt{H}/T1`$ where the vortex-induced $`\mathrm{\Gamma }_{tr}`$ dominates over the other mechanisms of scattering, including potential impurities.

Our approach also enables one to compute other averages, such as $`𝒢^R𝒢^Rϵ(\mathrm{\Gamma }_{tr}/r)^{1/2}e^{2iϵr\mathrm{\Gamma }_{tr}r}`$, which controls the effect of vortex disorder on superfluid stiffness measured by $`\mu `$SR

$$\rho _s(0)\rho _s(T)=\frac{1}{\pi k_F}\mathrm{Im}𝑑ϵ\mathrm{tanh}\left(\frac{ϵ}{2T}\right)$$
(17)
$$𝑑𝐫\mathrm{Tr}𝒢^R(ϵ,𝐫)𝒢^R(ϵ,𝐫)\frac{2\mathrm{ln}2}{\pi }T+\frac{H}{8\pi T}\mathrm{ln}\frac{T^2}{H},$$

Notably, in contrast to the spectral and transport characteristics whose behavior is determined by the structure of the electron spectral function near its maximum, Eq.(17) is governed by the overall momentum integral of the solution of the two-particle analog of Eq.(11).

This behavior is common amongst the thermodynamic quantities associated with the averages of local bi-linear operators $`\overline{\psi }_\sigma (𝐫)\psi _\sigma ^{}(𝐫)`$ which are invariant under the gauge transformation (1). In fact, such averages as, e.g., $`G^R(ϵ,\mathrm{𝟎})`$ must be computed differently, since now the trajectories contributing to the path integral (6) may deviate very strongly from the semiclassical one, $`𝐫_0(\tau )=\mathrm{𝟎}`$, although the mean square of the distance over which a typical trajectory $`𝐫(\tau )`$ ventures from the origin still scales quadratically with time, $`𝐫^2(\tau )\tau ^2`$. For $`1/\mathrm{\Delta }\tau `$ this allows one to evaluate the damping factors as

$$W_i[𝐫(\tau )]\mathrm{exp}\left[\frac{\tau ^2}{2}\frac{d𝐪}{(2\pi )^2}w_i(𝐪)\right].$$
(18)
Plugging (18) instead of (8) into (6) and computing the resulting (quadratic) path-integral, we arrive at the correction to the linear DOS corresponding to the clean limit

$$\nu (ϵ)=\mathrm{ImTr}\left[\widehat{\gamma }_0G^R(ϵ,\mathrm{𝟎})\right]=\frac{\sigma _H}{\pi ^{3/2}}\left(\frac{ϵ}{\sigma _H}\right),$$
(19)
where $`\sigma _H^2H\mathrm{ln}(\mathrm{\Delta }^2/H)`$ and $`(x)=\pi ^{1/2}x\mathrm{Erf}(x)+\mathrm{exp}(x^2)`$. The effect of vortex disorder is most pronounced at small energies $`ϵH^{1/2}`$, and it appears to be stronger than in the semiclassical (Doppler-only) approximation (cf. with ).

Directly, this DOS correction can be extracted from the tunneling conductance $`G(V)\nu (V)`$. Indirectly, it can also be manifested through the correction to electronic specific heat

$$C(T)=_0^{\mathrm{}}\frac{(ϵ/T)^2\nu (ϵ)dϵ}{\mathrm{cosh}^2(ϵ/2T)}\frac{18\zeta (3)T^2}{\pi }+\frac{\sigma _H^4}{16\pi T^2}$$
(20)
Notably, the correction $`\mathrm{\Delta }CH^2\mathrm{ln}(\mathrm{\Delta }^2/H)/T^2`$ is smaller than the result ($`\mathrm{\Delta }CH`$) obtained under in the situation where the inter-vortex repulsion is stronger than random pinning, and therefore the VLL is partially ordered .

Although the smallness of the disorder-induced term in Eq.(20) might hinder its detection, an alternate possibility is offered by the spin-lattice relaxation time

$$\frac{1}{T_1(T)}_0^{\mathrm{}}\frac{\nu ^2(ϵ)dϵ}{\mathrm{cosh}^2(ϵ/2T)}\frac{2}{3}T^3+\frac{\sigma _H^3}{3\pi ^{5/2}},$$
(21)
where we dropped the overall prefactor proportional to the ion-specific matrix elements.

In summary, we carried out a fully quantum mechnical analysis of the quasiparticle properties of the VLL phase of layered $`d`$-wave superconductors. We demonstrated that both the semiclassical Doppler shift and the intrinsically quantum mechanical BA scattering have comparable effects on all the observables. Our path-integral approach enabled us to identify the energy-dependent effective quasiparticle width $`\mathrm{\Gamma }(ϵ)H/ϵ`$ describing the near-maximum ($`ϵ^2𝐩^2`$) behavior of the distinctly non-Lorentzian electron spectral function which can be directly measured by ARPES.

Also, we exposed the striking difference between the two rates: $`\mathrm{\Gamma }_{tr}\mathrm{\Gamma }(T)H/T`$ displayed by the transport characteristics and $`\mathrm{\Gamma }^{}\mathrm{\Gamma }(H^{1/2})H^{1/2}`$ manifested by tunneling and thermodynamic quantities.

In contrast, the real-space averaging procedure applied in the previous studies of the disordered vortex states  is bound to deliver the latter rate, because, focusing solely on the energy distribution, it does not faithfully represent the momentum dispersion of the averaged electron spectral function.

We emphasize that the origin (that is, a non-trivial energy dependence) of the different apparent quasiparticle rates, as revealed by the different measurements, must be distinguished from the conventional juxtaposition of quasiparticle lifetime versus transport time in, e.g., photoemission and transport experiments.

We expect that, albeit derived under a number of simplifying assumptions, our main conclusions are robust against including such factors as spatial anisotropy and non-linearity of the fermion dispersion as well as inter-node scattering and, most importantly, modifying the distributions (3) with the purpose of describing partially ordered ”vortex glass” states .

Lastly, the results of this paper can also be used to describe the effects of thermal phase fluctuations controlled by $`\alpha T`$ in the pseudogap phase of the cuprates . To this end, we predict that different measurements may return different values of the effective quasiparticle width, which might explain the inconsistency between the widths deduced from the tunneling and ARPES data .

This research was supported by the NSF under Grant No. DMR-0071362.


text:
# 1 Introduction

## 1 Introduction

The reliable description of multi-jet production at the LHC  will be an important issue. This is not only related to the study of QCD in multi-parton final states but it is also very important in order to estimate several backgrounds for new physics effects. For instance, new unstable massive particles that decay to many partons may be discovered at the LHC only when a reliable description of these final states is established.

In this respect, apart from the problem of computing scattering matrix elements with many particles, also the efficient phase space generation is of great importance, because the scattering amplitudes in QCD exhibit strong peaking structures in phase space, which have to be taken into account by the generation algorithm. Flat phase space generators, like RAMBO , will not be adequate for this task. In the last years several methods to efficiently integrate the peaking structures of the scattering amplitudes have emerged, and have been used in several contexts . For instance, PHEGAS  is an example where an efficient, automated, mapping of all possible peaking structures of a given scattering process has been established. The algorithm is based on the “natural” mappings dictated by the Feynman graphs contributing to the given process, so that the number of kinematical channels used to generate the phase space is equal to the number of Feynman graphs. Using adaptive methods, like multi-channel optimization  and by throwing away channels that are negligible, we may end up with a few channel generator exhibiting high efficiency, as is indeed the case in $`n(+\gamma )`$-fermion production in $`e^+e^{}`$ collisions. In contrast, the QCD scattering amplitudes point towards the opposite direction: large number of Feynman graphs which means large number of kinematical channels which, moreover, contribute equally to the result.

A way out off this problem may be based on the long-standing remark that $`n+2`$-gluon amplitude may be described by a very compact expression when special helicities are assigned to the gluons, which, combined with the leading color approximation, results to

$$\underset{c}{}||^2=8\left(\frac{N_\mathrm{c}}{2}\right)^n(N_\mathrm{c}^21)\underset{1i<j}{\overset{n+2}{}}(p_ip_j)^4\underset{P(2,\mathrm{},n+2)}{}A_{n+2}(p_1,\mathrm{},p_{n+2}),$$
(1)
where $`N_\mathrm{c}`$ refers to the number of colors,

$$A_{n+2}(p_1,\mathrm{},p_{n+2}):=[(p_1p_2)(p_2p_3)\mathrm{}(p_{n+1}p_{n+2})(p_{n+2}p_1)]^1,$$
(2)
and the sum over all permutations of the $`2^{\text{nd}}`$ to the $`(n+2)^{\text{nd}}`$ argument of this function is taken, with the exception of those that are equivalent under reflection $`in+4i`$ .

SARGE  is the first known example of a phase space generator that deals with the momentum structures entering the above expression, namely with (2), known as antenna structures. The algorithm is based on the “democratic” strategy to generate the $`n`$ body phase space, as is the case for RAMBO, and it makes use of the scale symmetry of the antenna to achieve the required goal.

In this paper, we study the “hierarchical” strategy for phase space generation in order to efficiently map the momentum antenna structures. The idea is as follows. Using the standard two-body phase space (neglecting factors of $`2\pi `$)

$$d\mathrm{\Phi }_2(P;s_1,s_2;p_1,p_2):=d^4p_1\delta _+(p_1^2s_1)d^4p_2\delta _+(p_2^2s_2)\delta ^4(Pp_1p_2),$$
(3)
we decompose the phase space as

$`d\mathrm{\Phi }_n(P;p_1\mathrm{},p_n):=`$ $`\left({\displaystyle \underset{i=1}{\overset{n}{}}}d^4p_i\delta _+(p_i^2\sigma _i)\right)\delta ^4\left({\displaystyle \underset{i=1}{\overset{n}{}}}p_iP\right)`$
$`=`$ $`ds_{n1}d\mathrm{\Phi }_2(Q_n;\sigma _n,s_{n1};p_n,Q_{n1})`$
$`\times `$ $`ds_{n2}d\mathrm{\Phi }_2(Q_{n1};\sigma _{n1},s_{n2};p_{n1},Q_{n2})`$
$`\mathrm{}`$
$`\times `$ $`ds_2d\mathrm{\Phi }_2(Q_3;\sigma _3,s_2;p_3,Q_2)d\mathrm{\Phi }_2(Q_2;\sigma _2,\sigma _1;p_2,p_1).`$ (4)
The task is to express the phase space in terms of the invariants $`p_ip_j`$ appearing in the antenna structure (2), so that, using a suitable mapping, we can construct a density that, apart from constant and soft terms, will be identical to this antenna structure.

In the first section, we describe the basic building block of the algorithm, which is the expression of the two-body phase space in terms of the scaled invariants. In the second section, we demonstrate how this basic building block can be used in a sequential way to produce the full antenna. Finally in the third section, some details concerning the numerical implementation of the algorithm as well as comparisons with known generators is given. The Appendices present all relevant generation algorithms from which the exact functioning of the generator can be reconstructed.

## 2 The hierarchical antenna

### 2.1 The basic building block

To illustrate the idea we consider the generation of the 2-body phase space (3) when two massless antenna momenta $`q_1,q_2`$ are given. The momentum $`P`$ can be decomposed as

$$P^\mu =rq_1^\mu +A_1^\mu =:A_1^\mu +B_1^\mu \mathrm{with}r=P^2/(2Pq_1),$$

so that the Sudakov parameterization of $`p_1`$ is given by

$$p_1^\mu =a_1A_1^\mu +b_1B_1^\mu +k_1^\mu .$$

where the variables $`a_1`$ and $`b_1`$ are given by

$$a_1=\frac{p_1B_1}{A_1B_1},b_1=\frac{p_1A_1}{A_1B_1}.$$

The same can be done in terms of $`p_2`$ and $`q_2`$, and in the center-of-mass frame (CMF) of $`P`$, where $`P=(\sqrt{s},0,0,0)`$, $`\mathrm{cos}(\mathrm{}(\stackrel{}{q}_1,\stackrel{}{q}_2))=\mathrm{c}`$ and $`\mathrm{s}=\sqrt{1\mathrm{c}^2}`$, one can choose

$`A_1`$ $`=\frac{1}{2}\sqrt{s}(1,0,0,1)`$ $`,A_2`$ $`=\frac{1}{2}\sqrt{s}(1,0,\mathrm{s},\mathrm{c}),`$
$`B_1`$ $`=\frac{1}{2}\sqrt{s}(1,0,0,1)`$ $`,B_2`$ $`=\frac{1}{2}\sqrt{s}(1,0,\mathrm{s},\mathrm{c}),`$
$`k_1`$ $`=(0,x_1,y_1,0)`$ $`,k_2`$ $`=(0,x_2,y_2\mathrm{c},y_2\mathrm{s}).`$
The phase space can now be completely expressed in terms of $`a_1`$ and $`a_2`$, leading to

$$d\mathrm{\Phi }_2(P;p_1,p_2)=da_1da_2\mathrm{\Pi }^{(1/2)}\mathrm{\Theta }(\mathrm{\Pi }),$$
(5)
with

$`\mathrm{\Pi }(a_1,a_2)`$ $`=4\mathrm{s}^2[(1a_2+\overline{s}_2\overline{s}_1)a_2\overline{s}_2]`$
$`[(12a_1\overline{s}_1+\overline{s}_2)+(12a_2\overline{s}_1+\overline{s}_2)\mathrm{c}]^2,`$ (6)
where $`\overline{s}_{1,2}=s_{1,2}/s`$, and where $`\mathrm{\Theta }`$ is the step function. In terms of Lorentz invariants, the parameters $`a_1`$ and $`a_2`$ are given by

$$a_1=\frac{q_1p_1}{q_1P},a_2=\frac{p_2q_2}{Pq_2}.$$
(7)
So in order to obtain a two-body phase space with a density which depends on the invariants $`a_1,a_2`$ following some given function $`f(a_1,a_2)`$, one has to generate $`a_1,a_2`$ following a density proportional to $`f(a_1,a_2)\times \mathrm{\Pi }^{(1/2)}(a_1,a_2)`$ in the region where $`\mathrm{\Pi }(a_1,a_2)>0`$, and construct the momenta following the Sudakov parameterization. Explicitly, the direct construction is given by

$`p_1^0`$ $`(s+s_1s_2)/(2\sqrt{s}),`$
$`p_1^3`$ $`p_1^0\sqrt{s}a_1,`$
$`p_1^2`$ $`((\sqrt{s}p_1^0\sqrt{s}a_2)+\mathrm{c}p_1^3)/\mathrm{s},`$
$`p_1^1`$ $`ϵ((p_1^0)^2s_1(p_1^2)^2(p_1^3)^2)^{1/2},`$
where $`ϵ`$ should be a fair random variable which can take values $`+1`$ and $`1`$. For more details about this procedure, we refer the reader to the Appendix A.

### 2.2 Antenna generation

In the hierarchical/sequential approach, the generation strategy proceeds through a sequence of two-body phase space generations following the decomposition (4). At each two-body generation, one final-state momentum $`p_k`$ is generated, together with the sum $`Q_{k1}`$ of the remaining final-state momenta to be generated. This suggest to label the momenta in a way opposite to the order of generation, so first $`p_n,Q_{n1}`$ are generated, then $`p_{n1},Q_{n2}`$ and so on. The starting point is the CMF of the initial momenta $`q_1`$ and $`q_2`$ with $`Q_n=q_1+q_2`$ being the overall momentum. The CMF of momentum $`Q_k`$ we denote by CMFk. The pair $`p_k,Q_{k1}`$ is generated by generating variables $`a_1^{(k)},a_2^{(k)}`$ and constructing the momenta as described before. These variables are now equal to

$$a_1^{(k)}=\frac{p_{k+1}p_k}{p_{k+1}Q_k}\mathrm{and}a_2^{(k)}=\frac{q_2Q_{k1}}{q_2Q_k}.$$

This happens in CMFk, so in order to obtain $`p_k,Q_{k1}`$, the constructed momenta have to be boosted such that $`(\sqrt{Q_k^2},0,0,0)`$ is transformed to $`Q_k`$.

We would like to generate the momenta following a density that is proportional to

$$A_{n+2}(q_1,p_n\mathrm{},p_1,q_2)=[(q_1p_n)(p_np_{n1})\mathrm{}(p_1q_2)]^1.$$
(8)
Since the integrand is infra-red singular, a cutoff on the invariants is necessary. Therefore, we define a symmetric matrix $`\sigma _{ij}`$ which encodes the restrictions on the momenta through

$$\sigma _{ii}=\sigma _i=p_i^2,\sqrt{\sigma _i\sigma _j}\sigma _{ij}p_ip_j\text{and}\mathrm{\Sigma }_k:=\underset{i=1}{\overset{k}{}}\sigma _i.$$
(9)
Before we proceed, we do three observations. Firstly, we have

$$p_{k+1}Q_k=(Q_{k+1}^2Q_k^2p_{k+1}^2)/2.$$

Secondly, we have

$$\frac{s_{k+1}\mathrm{\Sigma }_{k+1}}{\left(s_{k+1}\sigma _{k+1}s_k\right)\left(s_k\mathrm{\Sigma }_k\right)}=\frac{d}{ds_k}\mathrm{log}\left(\frac{s_k\mathrm{\Sigma }_k}{s_{k+1}\sigma _{k+1}s_k}\right)=:g_{k+1}(s_k),$$
(10)
and thirdly, we can write

$$A_{n+2}(q_1,p_n\mathrm{},p_1,q_2)=\frac{1}{2^{n1}}(s_n\mathrm{\Sigma }_n)(q_1Q_n)(q_2Q_n)\left(\underset{k=n}{\overset{3}{}}g_k(s_{k1})\frac{1}{a_1^{(k)}a_2^{(k)}}\right)\frac{1}{a_1^{(2)}a_2^{(2)}},$$

with $`s_n=Q_n^2`$, $`p_{n+1}=q_1`$ and $`Q_1=p_1`$. These observations suggest that the phase space generation

$`ds_{n1}g_n(s_{n1})`$ $`da_1^{(n)}{\displaystyle \frac{1}{a_1^{(n)}}}`$ $`da_2^{(n)}{\displaystyle \frac{1}{a_2^{(n)}}}\mathrm{\Pi }_{(n)}^{(1/2)}\mathrm{\Theta }(\mathrm{\Pi }_{(n)})`$
$`ds_{n2}g_{n1}(s_{n2})`$ $`da_1^{(n1)}{\displaystyle \frac{1}{a_1^{(n1)}}}`$ $`da_2^{(n1)}{\displaystyle \frac{1}{a_2^{(n1)}}}\mathrm{\Pi }_{(n1)}^{(1/2)}\mathrm{\Theta }(\mathrm{\Pi }_{(n1)})`$
$`\mathrm{}`$
$`ds_2g_3(s_2)`$ $`da_1^{(3)}{\displaystyle \frac{1}{a_1^{(3)}}}`$ $`da_2^{(3)}{\displaystyle \frac{1}{a_2^{(3)}}}\mathrm{\Pi }_{(3)}^{(1/2)}\mathrm{\Theta }(\mathrm{\Pi }_{(3)})`$
$`da_1^{(2)}{\displaystyle \frac{1}{a_1^{(2)}}}`$ $`da_2^{(2)}{\displaystyle \frac{1}{a_2^{(2)}}}\mathrm{\Pi }_{(2)}^{(1/2)}\mathrm{\Theta }(\mathrm{\Pi }_{(2)}),`$ (11)
will lead to a density for the momenta that is proportional to $`A_{n+2}`$. Three variables are generated in each CMFk, namely $`s_{k1}`$, $`a_1^{(k)}`$ and $`a_2^{(k)}`$. Just as the integration of $`s_{k1}`$ (10), also the integration of $`a_1^{(k)},a_2^{(k)}`$ results in a volume factor that depends on the corresponding variables generated in CMFk+1. As we will show in Appendix A, however, these factors are logarithmic functions of their arguments and exhibit a non-singular behavior, and we call them soft factors. The total actual density will therefore be the product of $`n1`$ soft factors times the antenna structure under consideration.

In the end, we want to generate all permutations in the momenta of (8). Those for which $`q_1`$ and $`q_2`$ each appear in two factors (none of which is $`q_1q_2`$) cannot be obtained by simple re-labeling. In order to obtain these, we observe that they can be decomposed into two antennas, namely

$$A_{m+2}(q_1,p_m,p_{m1},\mathrm{},p_2,p_1,q_2)\times A_{nm+2}(q_2,p_n,p_{n1},\mathrm{},p_{m+2},p_{m+1},q_1)$$
(12)
and each of these can be generated after the decomposition,

$`d\mathrm{\Phi }_n(P;p_1\mathrm{},p_n)`$ $`=ds_mds_{nm}d\mathrm{\Phi }_2(Q_n;s_m,s_{nm};Q_m,Q_{nm})`$
$`\times d\mathrm{\Phi }_m(Q_m;p_1,\mathrm{},p_m)d\mathrm{\Phi }_{nm}(Q_{nm};p_{m+1},\mathrm{},p_n).`$ (13)
In order to combine the two sub-antennas to the required antenna structure, we have to take into account in the first decomposition a density that is proportional to

$$\frac{\mathrm{\Theta }(\sqrt{s_n}\sqrt{s_m}\sqrt{s_{nm}})}{(q_1Q_m)(q_1Q_{nm})(q_2Q_m)(q_2Q_{nm})s_ms_{nm}}.$$
(14)
The case of $`m=1`$ is still special. Then, the first step in (11) should be replaced by

$$ds_{n1}g^{}(s_{n1})da_1^{(n)}da_2^{(n)}\mathrm{\Pi }_{(n)}^{(1/2)}\mathrm{\Theta }(\mathrm{\Pi }_{(n)})\frac{1}{a_1^{(n)}a_2^{(n)}(1a_1^{(n)})(1a_2^{(n)})},$$
(15)
and the rest of the sequence should go on with the replacement of $`p_n`$ by $`p_{n+1}`$ at the second step. For the density $`g^{}`$ we refer to Appendix A. There, we have collected all integrals and generation algorithms related to the antenna generation described so far.

### 2.3 Open antennas

As it will be clear from the numerical analysis presented in the next section, the soft factors appearing in the description of the QCD antenna contribute to a certain extent to the variance of the Monte Carlo integration. There is an alternative approach, that still follow the hierarchical/sequential generation strategy, and give better results. It based on the observation that the production of an ‘open’ antenna structure, namely one where the last product $`q_2p_1`$ is missing, is simpler, since it can be constructed without using the variables $`a_2^{(k)}`$. They can be replaced by flatly generated azimuthal angles, to that $`da_2^{(k)}\mathrm{\Pi }_{(k)}^{(1/2)}\mathrm{\Theta }(\mathrm{\Pi }_{(k)})d\phi ^{(k)}`$. The basic decomposition therefore becomes

$$\begin{array}{c}ds_{n1}g_n(s_{n1})d\phi ^{(n)}da_1^{(n)}\frac{1}{a_1^{(n)}}\times ds_{n2}g_{n1}(s_{n2})d\phi ^{(n1)}da_1^{(n1)}\frac{1}{a_1^{(n1)}}\times \mathrm{}\hfill \\ \hfill \mathrm{}\times ds_2g_3(s_2)d\phi ^{(3)}da_1^{(3)}\frac{1}{a_1^{(3)}}\times d\phi ^{(2)}da_1^{(2)}\frac{1}{a_1^{(2)}}.\end{array}$$
(16)
This way, we will get the antenna density (8) without the factor $`p_1q_2`$ in the denominator, which is the reason why we call this an ‘open’ antenna. A ‘closed’ antenna can be obtained using the fact that, by combining two open antennas, one can choose for another factor from the antenna string to be missing. Then, a multi-channeling procedure can be performed with these different choices, leading to a density that is, roughly speaking, proportional to

$$\frac{(q_1p_n)+(p_np_{n1})+\mathrm{}+(p_3p_2)+(p_2p_1)+(p_1q_2)}{(q_1p_n)(p_np_{n1})\mathrm{}(p_3p_2)(p_2p_1)(p_1q_2)}.$$

To get the different choices, a first splitting of $`Q_n`$ into $`Q_m`$ and $`Q_{nm}`$ has to be performed, after which open antennas are generated from each of these, one with $`q_1`$ and the other with $`q_2`$ as initial momentum. For details, we refer to Appendix B.

## 3 Results

In this section, we present results obtained by SARGE and HAAG333HAAG stands for: Hierarchical AntennA Generation., the program that implements the hierarchical algorithm described before. In order to be as general as possible, the only cut we apply is

$$(p_i+p_j)^2s_0,$$

where $`i,j(ij)`$ runs from $`1`$ to $`n+2`$ where $`n`$ is the number of final-state particles. Unless explicitly mentioned differently, we use $`s_0=900\mathrm{GeV}^2`$ and the total energy $`\sqrt{s}=1000\mathrm{GeV}`$. Moreover, all particles are assumed to be massless in order to compare with SARGE, with which only massless particles can be treated.

As it was mentioned in the introduction, we are interested in integrating sums of QCD antenna structures (2). We start by considering the simplest case, namely integrating the function

$$s^2[(p_1p_3)(p_3p_4)(p_4p_2)(p_2p_5)\mathrm{}(p_{n+2}p_1)]^1$$
(17)
that corresponds to a given permutation of the momenta, namely $`(1,3,4,2,5,\mathrm{},n+2)`$. In Table 1 we give the results for SARGE, HAAG with open antenna generation, and HAAG(C) with closed antenna generation. In all three codes the same single channel, corresponding to (17), has been used in the generation. $`N_{\mathrm{gen}}`$ and $`N_{\mathrm{acc}}`$ are the number of generated and accepted events, and by $`f`$ we define

$$f:=\frac{V_2}{I^2},$$

where $`V_2`$ is the quadratic variance and $`I`$ is the estimated integral. $`f`$ is clearly a measure of the efficiency of the generator. Moreover $`\epsilon `$, defined as

$$\epsilon :=\frac{<w>}{w_{\mathrm{max}}},$$

is the usual generation efficiency related for instance to ‘unweighted’ events in a realistic simulation.

The results agree well, and exhibit the fact that the generated densities of the generators the hierarchical type are much closer to the integrand. Moreover, the closed antenna algorithm HAAG(C) becomes less efficient compared to the open one as the number of particles increases. The same picture is reproduced for an arbitrary permutation.

For a realistic QCD calculation, the integrated function may be approximated by a sum over permutations. Therefore, an efficient generator has to include all possible channels, where each channel corresponds to a given permutation of the momenta. In that case, a multi-channeling optimization procedure can be applied, which is incorporated in HAAG. In order to study the efficiency of this optimization we consider the same integration as before, but with all channels contributing to the generation and allowing for optimization. In this optimization procedure, we discard channels that contribute less than a certain pre-determined fraction to the set of available channels. It is expected, of course, that in end the right permutation will be ‘chosen’ by the optimization. This is indeed the case and the results are presented in Table 2. We see that the optimization results to a picture close to the one obtained with the single channel generation, with some noticeable improvement in the case of SARGE. We also include results with SARGE.n, a slightly different version, described in Appendix C.

As is the case for any multi-channel generator, a computational complexity problem arises when the number of channels increases. For instance, in our case we are facing a number of $`\frac{1}{2}(n+1)!`$ channels! On the other hand, it is also clear that the channels we are considering have a large overlap in most of the available phase space. It is therefore worth to investigate the dependence of the integration efficiency on the number of channels used. This is presented in Table 3, where the full antenna

$$s^2\underset{P(2,\mathrm{},n+2)}{}[(p_1p_3)(p_3p_4)(p_4p_2)(p_2p_5)\mathrm{}(p_{n+2}p_1)]^1$$
(18)
is integrated, using a number of channels that has been selected on a random basis. We see the rather interesting phenomenon that a decent description can be achieved with a much smaller number of channels.

Variations of this technique of using only subsets of channels, for example choosing another subset after each step of multi-channel optimization, lead to the same picture.

The complete results of the integration of the full antenna are presented in Table 4.

We see that HAAG has a much better $`f`$ factor than SARGE. On the other hand the $`\epsilon `$ exhibits a less dramatic effect. This is related to the fact that SARGE generates a phase space that is much larger than the one defined by the cut on $`s_0`$. In that sense, if the main time consumption in a given computation is spent over the evaluation of the integrand (matrix element squared), it is more fair to compare the square of the estimated expected error, normalized by the number of accepted events $`N_{\mathrm{acc}}`$. In that case we see that HAAG is still 2-3 times more efficient, and if we consider a smaller cut, namely $`\sqrt{s_0}=10`$ GeV, this gain goes up to an order of magnitude (Table 5).

For a realistic calculation of the cross section of a QCD process, one may assume that the time it takes to perform one evaluation of the integrand is much larger than the time it takes to generate one accepted event and to calculate the weight. This means that the computing time is completely determined by the number of accepted events $`N_{\mathrm{acc}}`$. We introduce

$$\frac{N_{\mathrm{acc}}f}{N_{\mathrm{gen}}}$$

as a measure of the computing time. For a realistic calculation, one has to multiply this number by the evaluation time of the integrand, and devide by the square of the relative error one wants to reach. Figure 1 shows this quantity as function of the number of produced partons using the data of Table 5. According to this graph, a calculation with SARGE would take $`10`$ times longer than the calculation with HAAG.

Figure 1: $`N_{\mathrm{acc}}f/N_{\mathrm{gen}}`$ (a measure of computing time) as function of the number of produced partons.

Finally, Figure 2 shows the dependence of the result on the value of the infrared cut-off $`\sqrt{s_0}`$ and the number of produced partons. The function clearly exhibits a negative power behavior for $`s_0`$. Moreover, the curve becomes steeper with increasing number of jets, suggesting that at least the leading power is related to the final-state multiplicity.

Figure 2: The integral of the full antenna as function of $`\sqrt{s_0}`$ for different values of the number of produced partons.

## 4 Conclusions

HAAG exhibits an improved efficiency compared to SARGE for multi-parton calculations. It is also more powerful in describing densities where a partial symmetrization over the permutation space is considered. Finally, HAAG makes no fundamental distinction among massless and massive particles, so it can be used for an arbitrary multi-partonic process.

## Acknowledgments

The authors would like to thank Ronald Kleiss for proof reading this paper. The research has been financially supported by the European Union under contract number HPRN-CT-2000-00149.

## Appendices

The following appendices contain details about the generation of the various random variables necessary to build phase space with the desired density. The techniques used to achieve this are inversion, rejection and multi-channeling. For details about these techniques, we refer to . We only want to mention that inversion is applicable if one has an analytic expression (with reasonable complexity) of the inverse of the indefinite integral of the density. Rejection can be applied if one knows a function which is strictly larger than the density and is proportional to a density one is able to generate. The efficiency is given by one divided by the integral of that function. Multi-channeling can be used if the density can be written as the weighted sum of densities, each of which one is able to generate, and if the weights are positive.

### A Closed antenna

In this appendix we present all relevant algorithms for the generation of the closed antenna. More specifically, in the following we describe the generation of the three variables $`s_{k1}`$, $`a_1^{(k)}`$ and $`a_2^{(k)}`$ that are needed to describe the antenna at the $`k`$th CMF. Most of the time the superscript $`(k)`$ is omitted for convenience.

#### A1 Generation of $`𝒔_{𝒌\mathbf{}\mathrm{𝟏}}`$

Since the indefinite integral is given in (10), we see that $`s_{k1}`$ can be generated with inversion. The limits are given by

$$s_k\mathrm{\Delta }_{k1}s_{k1}\mathrm{\Lambda }_{k1}\text{with}\mathrm{\Lambda }_k=\mathrm{\Sigma }_k+\underset{ij}{\overset{k}{}}\sigma _{ij},\mathrm{\Delta }_k=\sigma _{k+1}+2\underset{i=1}{\overset{k}{}}\sigma _{k+1,i},$$
(19)
and the weight factor is

$$\mathrm{log}\left(\frac{s_k\mathrm{\Sigma }_{k1}\mathrm{\Delta }_{k1}}{\mathrm{\Delta }_{k1}\sigma _k}\right)+\mathrm{log}\left(\frac{s_k\sigma _k\mathrm{\Lambda }_{k1}}{\mathrm{\Lambda }_{k1}\mathrm{\Sigma }_{k1}}\right).$$

In the case of $`m=1`$ (15) we simply have

$$𝑑s_{k1}\frac{1}{s_{k1}\mathrm{\Sigma }_{k1}}=\mathrm{log}(s_{k1}\mathrm{\Sigma }_{k1}),$$
(20)
so that the weight factor is

$$\mathrm{log}\left(s_k\mathrm{\Delta }_{k1}\mathrm{\Sigma }_{k1}\right)\mathrm{log}\left(\mathrm{\Lambda }_{k1}\mathrm{\Sigma }_{k1}\right).$$

#### A2 Generation of $`𝒂_\mathrm{𝟏}\mathbf{,}𝒂_\mathrm{𝟐}`$

The general integral corresponding to the generation of $`a_1,a_2`$ is given by

$$𝑑a_1𝑑a_2\frac{\mathrm{\Theta }(\mathrm{\Pi }(a_1,a_2))\mathrm{\Theta }(a_1a_1^{(0)})\mathrm{\Theta }(a_2a_2^{(0)})}{a_1a_2\sqrt{\mathrm{\Pi }(a_1,a_2)}},$$
(21)
where $`\mathrm{\Pi }`$ is defined in (6), and $`a_1^{(0)},a_2^{(0)}`$ are possibly necessary infra-red cut-offs. We shall analyze this integral by first integrating over the $`a_2`$-variable, and then over the $`a_1`$-variable. The generation has then to be performed in the opposite order: first $`a_1`$ and then $`a_2`$. As we shall see, the inclusion of the cut-off on $`a_1`$ does not lead to complications, but the inclusion of both cut-offs does. For that case, we see two solutions. Firstly, we can replace the integral by

$$𝑑a_1𝑑a_2\frac{\mathrm{\Theta }(\mathrm{\Pi }(a_1,a_2))\mathrm{\Theta }(a_1a_1^{(0)})}{a_1(a_2+h)\sqrt{\mathrm{\Pi }(a_1,a_2)}},$$

where $`h`$ is related to $`a_2^{(0)}`$. This, of course, changes the actual density with which $`a_1`$ and $`a_2`$ are generated, but may, for small $`h`$, still be considered suitable for the desired antenna structure. Also, the final antenna generator will cover phase space less efficiently (will generate ‘to much’ phase space), simply because less cuts are included analytically. In the second solution, we write the integrand as the sum

$$\frac{\mathrm{\Theta }(\mathrm{\Pi }(a_1,a_2))\mathrm{\Theta }(a_1a_1^{(0)})}{a_1(a_2+a_1)\sqrt{\mathrm{\Pi }(a_1,a_2)}}+\frac{\mathrm{\Theta }(\mathrm{\Pi }(a_1,a_2))\mathrm{\Theta }(a_2a_2^{(0)})}{(a_1+a_2)a_2\sqrt{\mathrm{\Pi }(a_1,a_2)}},$$

which is exactly equal to the original integrand on the phase space for which both cuts are included. Both integrands can be integrated analytically, so that the multi-channeling procedure can be applied to generate their sum. For this solution, the only problem is that, again, ‘to much’ phase space is generated.

##### A2.1 The $`𝒂_\mathrm{𝟐}`$-variable

In order to integrate over the $`a_2`$-variable, $`\mathrm{\Pi }(a_1,a_2)`$ is more conveniently written as

$$\mathrm{\Pi }(a_1,a_2)=4(a_2^+(a_1)a_2)(a_2a_2^{}(a_1)),$$
(22)
with

$$\begin{array}{c}a_2^\pm (a_1):=\frac{1}{2}\left(1+\frac{s_{k1}\sigma _k}{s_k}+\left(12a_1\frac{s_{k1}\sigma _k}{s_k}\right)\mathrm{c}\right)\hfill \\ \hfill \pm \mathrm{s}\left(a_1\left(1a_1\frac{s_{k1}\sigma _k}{s_k}\right)\frac{\sigma _k}{s_k}\right)^{1/2}.\end{array}$$
(23)
Then, the general $`a_2`$-integral is given by

$$_{a_2^{}}^{a_2^+}𝑑a_2\frac{1}{(a_2+h)\sqrt{(a_2^+a_2)(a_2a_2^{})}},$$

with $`h=0`$ if no cut-off on $`a_2`$ is desirable, and $`h`$ related to $`a_2^{(0)}`$ or $`h=a_1`$, depending on the solution mentioned above if a cut-off is desirable. By substituting $`a_2a_2h`$, this integral can be written as

$$_{\overline{a}_2^{}}^{\overline{a}_2^+}𝑑a_2\frac{1}{a_2\sqrt{(\overline{a}_2^+a_2)(a_2\overline{a}_2^{})}}=\frac{2}{\sqrt{\overline{a}_2^+\overline{a}_2^{}}}\left[\mathrm{arctan}\left(\frac{\overline{a}_2^{}(\overline{a}_2^+a_2)}{\overline{a}_2^+(a_2\overline{a}_2^{})}\right)\right]_{\overline{a}_2^{}}^{\overline{a}_2^+}=\frac{\pi }{\sqrt{\overline{a}_2^+\overline{a}_2^{}}},$$

with $`\overline{a}_2^\pm =a_2^\pm +h`$. The explicit indefinite integral shows that the variable $`a_2`$ can be generated by inversion.

##### A2.2 The $`𝒂_\mathrm{𝟏}`$-variable

We start this section by mentioning that, in the case $`p_{k+1}^2=\sigma _{k+1}0`$, the variables

$$a_1^{(k)}=\frac{p_kp_{k+1}}{Q_kp_{k+1}}$$

can be expressed in terms of the ‘massless’ $`\stackrel{~}{a}_1^{(k)}`$ defined in terms of the ‘long’ component $`p_{k+1}^{(L)}`$ of $`p_{k+1}`$ in the CMFk: if $`p_{k+1}^{(L)}=(p_{k+1}^0,\beta _k^1\stackrel{}{p}_{k+1})`$ and $`p_{k+1}^{(S)}=(p_{k+1}^0,\beta _k^1\stackrel{}{p}_{k+1})`$ with

$$\beta _k:=\frac{|\stackrel{}{p}_{k+1}|}{p_{k+1}^0}=\frac{\sqrt{\lambda (s_{k+1},\sigma _{k+1},s_k)}}{s_{k+1}\sigma _{k+1}s_k},$$
(24)
then $`p_{k+1}=\frac{1+\beta _k}{2}p_{k+1}^{(L)}+\frac{1\beta _k}{2}p_{k+1}^{(S)}`$, and

$$\stackrel{~}{a}_1^{(k)}:=\frac{p_kp_{k+1}^{(L)}}{Q_kp_{k+1}^{(L)}}=\beta _k^1a_1^{(k)}h_1\text{with}h_1:=\frac{1\beta _k}{2\beta _k}\left(1+\frac{\sigma _ks_{k1}}{s_k}\right).$$

Since the relation between $`a_1^{(k)}`$ and $`\stackrel{~}{a}_1^{(k)}`$ is linear, the two-body phase space is still expressible in terms of $`a_1^{(k)}`$ scaled by $`\beta _k`$. The generic $`a_1`$-integral is given by

$$_{a_1^{\mathrm{min}}}^{a_1^{\mathrm{max}}}𝑑a_1\frac{1}{(a_1+h_1)\sqrt{(a_2^+(a_1)+h(a_1))(a_2^{}(a_1)+h(a_1))}},$$
(25)
with $`a_2^\pm `$ as defined in (23), $`h`$ is a constant (for this integral) or $`h(a_1)=a_1`$. The kinematical integration limits, coming from the requirement that $`a_2^\pm (a_1)`$ are real, are given by

$$a_1^\pm =\frac{1}{2}\left(1+\frac{\sigma _ks_{k1}}{s_k}\pm \sqrt{\lambda (1,\frac{\sigma _k}{s_k},\frac{s_{k1}}{s_k})}\right).$$

In the massless case we get $`a_1^{}=0,a_1^+=1s_{k1}/s_k`$ and we have to impose a lower bound on $`a_1`$ given by $`a_1^{(0)}=\sigma _{k+1,k}/p_{k+1}Q_k`$. In all cases for the form of $`h(a_1)`$, the $`a_1`$-integral is of the type

$$_{a_1^{\mathrm{min}}}^{a_1^{\mathrm{max}}}da_1\frac{1}{(a_1+h_1)f(a_1)}\text{with}f(a_1):=\sqrt{a_1^2+2va_1+w^2},v^2<w^2,$$
(26)
and the indefinite integral is given by

$$𝑑a_1\frac{1}{(a_1+h_1)f(a_1)}=\frac{1}{f(h_1)}\mathrm{log}\left(\frac{a_1+h_1+f(a_1)f(h_1)}{a_1+h_1+f(a_1)+f(h_1)}\right),$$

which is analytically invertible. The definite integral (times $`\pi /4`$ from the $`a_2`$-integral and a factor $`1/\beta _k`$ if $`p_{k+1}`$ is massive) gives the weight factors in the generation of $`a_1,a_2`$ for the closed antenna.

##### A2.3 The $`𝒂_\mathrm{𝟏}`$-variable in the case $`𝒎\mathbf{=}\mathrm{𝟏}`$

This refers to (15), where we have $`\mathrm{c}=1`$ and $`a_2=a_1+\mu `$ with

$$\mu :=\frac{s_{k1}\sigma _k}{s_k},\text{so that}_{a_1^{\mathrm{min}}}^{a_1^{\mathrm{max}}}da_1\frac{1}{a_1(1a_1)(a_1+\mu )(1\mu a_1)}$$
(27)
is the integral to be performed. Since the integrand is equal to

$$\frac{1}{a_1(a_1+\mu )}+\frac{1}{a_1(1\mu a_1)}+\frac{1}{(a_1+\mu )(1a_1)}+\frac{1}{(1a_1)(1\mu a_1)},$$

we see that the generation of $`a_1`$ can be done easily using the multi-channel technique with four channels with weight

$$w_i=\frac{1}{g_i^{(2)}e_ig_i^{(1)}d_i}\left(\mathrm{log}\left(\frac{g_i^{(1)}a_1^{\mathrm{max}}+d_i}{g_i^{(2)}a_1^{\mathrm{max}}+e_i}\right)\mathrm{log}\left(\frac{g_i^{(1)}a_1^{\mathrm{min}}+d_i}{g_i^{(2)}a_1^{\mathrm{min}}+e_i}\right)\right),$$

where

| $`i`$ | $`d_i`$ | $`e_i`$ | $`g_i^{(1)}`$ | $`g_i^{(2)}`$ |
| --- | --- | --- | --- | --- |
| 1 | 0 | $`\mu `$ | $`+`$ | $`+`$ |
| 2 | 0 | $`1\mu `$ | $`+`$ | $``$ |
| 3 | $`\mu `$ | $`1`$ | $`+`$ | $``$ |
| 4 | $`1`$ | $`1\mu `$ | $``$ | $``$ |

The soft factor will simply be equal to $`_{i=1}^4w_i`$. The integration over the azimuthal angle, replacing the $`a_2`$-integration, gives an extra factor of $`2\pi `$.

#### A3 Antenna split

With $`q_1(1,0,01)`$ and $`q_2(1,0,0,1)`$, the two-body phase space integral (14) assumes the form

$$𝑑s_1𝑑s_2𝑑Q_z\frac{\mathrm{\Theta }(\sqrt{s}\sqrt{s_1}\sqrt{s_2})\mathrm{\Theta }(s_1s_1^{(0)})\mathrm{\Theta }(s_2s_2^{(0)})}{s_1s_2(E_1(s_1,s_2)^2Q_z^2)(E_2(s_1,s_2)^2Q_z^2)},$$
(28)
with the energies

$$E_1(s_1,s_2):=\frac{s+s_1s_2}{2\sqrt{s}},E_2(s_1,s_2):=\sqrt{s}E_1(s_1,s_2),$$

and where $`s_{1,2}^{(0)}`$ are the sums of the matrix elements $`\sigma _{ij}`$ (9) corresponding with the momenta in two antennas to be generated. $`Q_z`$ is integrated between the kinematical limits $`\pm \sqrt{E_1^2s_1^2}`$, and can be treated in a way similar to the one described in the previous paragraph, by multi-channeling over four channels with

| $`i`$ | $`d_i`$ | $`e_i`$ | $`g_i^{(1)}`$ | $`g_i^{(2)}`$ |
| --- | --- | --- | --- | --- |
| 1 | $`E_1`$ | $`E_2`$ | $`+`$ | $`+`$ |
| 2 | $`E_1`$ | $`E_2`$ | $`+`$ | $``$ |
| 3 | $`E_1`$ | $`E_2`$ | $`+`$ | $``$ |
| 4 | $`E_1`$ | $`E_2`$ | $``$ | $``$ |

The final weight is the sum of the channel-weights, divided by the (dimensionful) factor $`4E_1E_2`$. The generation of $`s_1,s_2`$ is a specific case of a more general problem described in Appendix B2.

### B Open antenna

The three variables needed to describe the open antenna are $`s_{k1}`$, $`a_1^{(k)}`$ and $`\phi ^{(k)}`$. In each CMFk, the $`\phi ^{(k)}`$-variable should be generated with uniform distribution between $`0`$ and $`2\pi `$, and the $`a_1^{(k)}`$-variable should be distributed following $`1/a_1^{(k)}`$ between

$`a_{1,+}^{(k)}`$ $`={\displaystyle \frac{s_k+\sigma _ks_{k1}+\beta _k\sqrt{\lambda (s,\sigma _k,s_{k1})}}{2s_k}}`$
$`a_{1,}^{(k)}`$ $`=\mathrm{max}[{\displaystyle \frac{s_k+\sigma _ks_{k1}\beta _k\sqrt{\lambda (s,\sigma _k,s_{k1})}}{2s_k}},{\displaystyle \frac{\sigma _{k+1,k}}{p_{k+1}Q_k}}],`$
with $`\beta _k`$ as defined in (24). The normalized density for these generations is equal to

$$\frac{2\beta _k}{\pi L_ka_1^{(k)}}=\frac{\sqrt{\lambda (s_{k+1},\sigma _{k+1},s_k)}}{\pi L_k(p_{k+1}p_k)}\text{with}L_k=\mathrm{log}(a_{1,+}^{(k)}/a_{1,}^{(k)}).$$

This suggests to use

$`g_{k+1}(s_k):=`$ $`{\displaystyle \frac{\sqrt{\lambda (s_{k+1},\sigma _{k+1},\mathrm{\Sigma }_k)}}{(s_k\mathrm{\Sigma }_k)\sqrt{\lambda (s_{k+1},\sigma _{k+1},s_k)}}}`$
$`=`$ $`{\displaystyle \frac{d}{ds_k}}\mathrm{log}\left({\displaystyle \frac{s_k\mathrm{\Sigma }_k+\sqrt{\lambda (s_{k+1},\sigma _{k+1},s_k)}\sqrt{\lambda (s_{k+1},\sigma _{k+1},\mathrm{\Sigma }_k)}}{s_k\mathrm{\Sigma }_k+\sqrt{\lambda (s_{k+1},\sigma _{k+1},s_k)}+\sqrt{\lambda (s_{k+1},\sigma _{k+1},\mathrm{\Sigma }_k)}}}\right)`$
instead of (10) for the case that $`\beta _k1`$. The logarithm contributes to the weight as a soft factor again. For small values of the squared masses $`\sigma _k`$, the factor in the numerator will be canceled by $`s_{k+1}\mathrm{\Sigma }_{k+1}`$ in the denominator of $`g_{k+2}(s_{k+1})`$. Just as in the case of the closed antenna, we end up with one remaining, and desirable, factor $`s_2\mathrm{\Sigma }_2=2(p_2p_1)`$ in the denominator of the open antenna density, which cannot be achieved by the generation of the $`a_1^{(k)}`$-variables.

Let us denote the soft factor coming from the $`s_k`$-generation by

$$G_{k+1}:=_{\mathrm{\Lambda }_k}^{s_{k+1}\mathrm{\Delta }_k}𝑑s_kg_{k+1}(s_k),$$

then the complete density resulting from the decomposition (16) is

$$\frac{dD_0(Q_n;p_n,\mathrm{},p_1)}{d\mathrm{\Phi }_n(Q_n;p_n,\mathrm{},p_1)}=\frac{(q_1Q_n)\sqrt{\lambda (s_n,\sigma _n,\mathrm{\Sigma }_{n1})}B_n}{(q_1p_n)(p_np_{n1})\mathrm{}(p_3p_2)(p_2p_1)},$$

with soft factor

$`B_n:={\displaystyle \frac{1}{(4\pi )^{n1}G_3L_3L_2}}{\displaystyle \underset{k=n}{\overset{4}{}}}{\displaystyle \frac{1}{G_kL_k}}\sqrt{1{\displaystyle \frac{4\sigma _{k1}\mathrm{\Sigma }_{k2}}{(s_{k1}\mathrm{\Sigma }_{k1})^2}}}.`$
The factor $`p_1q_2`$ is missing in the denominator of the density, which is the reason why we call this an open antenna. Other open antennas can be obtained by starting with a decomposition of $`Q_n`$ into two momenta, from each of which open antennas of the above type are generated, with initial momentum $`q_1`$ for the one, and $`q_2`$ for the other.

In order to digress about this procedure, let us extend the labeling a bit. With a set $`\{I_1,I_2,\mathrm{},I_n\}`$ of $`n`$ non-equal labels, we write

$$Q_{I_k}:=\underset{m=1}{\overset{k}{}}p_{I_m},s_{I_k}:=Q_{I_k}^2$$

and so on. Now take $`n=M+N`$ and let

$$\{I_1,\mathrm{},I_N\}:=\{M+1,M+2,\mathrm{},M+N\},\{J_1,\mathrm{},J_M\}:=\{M,M1,\mathrm{},2,1\}.$$

We introduce $`p_{I_{N+1}}=q_1`$, $`p_{J_{M+1}}=q_2`$, and the decomposition

$`ds_{J_M}ds_{I_N}d\mathrm{\Phi }_2(Q_n;s_{I_N},s_{J_M};Q_{I_N},Q_{J_M})`$
$`\times dD_0(Q_{I_N};p_{I_N},p_{I_{N1}},\mathrm{},p_{I_1})dD_0(Q_{J_M};p_{J_M},p_{J_{M1}},\mathrm{},p_{J_1}),`$
which produces the density

$`{\displaystyle \frac{dD_M(Q_n;p_n,\mathrm{},p_1)}{d\mathrm{\Phi }_n(Q_n;p_n,\mathrm{},p_1)}}`$ $`={\displaystyle \frac{2Q_n^2B_{J_M}B_{I_N}}{\pi \sqrt{\lambda (Q_n^2,s_{J_M},s_{I_N})}}}\times {\displaystyle \frac{(p_{M+1}p_M)}{(q_1p_n)(p_np_{n1})\mathrm{}(p_2p_1)(p_1q_2)}}`$
$`\times (q_1Q_{I_N})(q_2Q_{J_M})\sqrt{\lambda (s_{I_N},\sigma _{I_N},\mathrm{\Sigma }_{I_{N1}})\lambda (s_{J_M},\sigma _{J_M},\mathrm{\Sigma }_{J_{M1}})}.`$ (29)
In order to cancel the ‘undesirable’ factors on the second line, we need to take care of the generation of $`s_{I_N},s_{J_M},Q_{I_N},Q_{J_M}`$.

#### B1 Generation of $`𝑸_{𝑰_𝑵}\mathbf{,}𝑸_{𝑱_𝑴}`$

Since $`q_1(1,0,0,1)`$ and $`q_2(1,0,0,1)`$, we can write

$$a_1=\frac{q_1Q_{I_N}}{q_1Q_n}\text{and}\frac{Q_{J_M}q_2}{Q_nq_2}=a_1+\mu \text{with}\mu =\frac{s_{J_M}s_{I_N}}{s_n},$$

so that the generation of an azimuthal angle $`\phi `$ between $`0`$ and $`2\pi `$ with the uniform distribution, and the generation of $`a_1`$ with a density proportional to

$$\frac{1}{a_1(a_1+\mu )}\text{between}a_1^\pm =\frac{s_n+s_{I_N}s_{J_M}\pm \sqrt{\lambda (s_n,s_{I_N},s_{J_M})}}{2s_n}$$

leads to the total density

$$d\mathrm{\Phi }_2(Q_n;s_{I_N},s_{J_M};Q_{I_N},Q_{J_M})\frac{2\mu }{\pi \mathrm{log}\frac{a_1^+(a_1^{}+\mu )}{a_1^{}(a_1^++\mu )}}\frac{(q_1Q_n)(Q_nq_2)}{(q_1Q_{I_N})(Q_{J_M}q_2)}.$$

#### B2 Generation of $`𝒔_{𝑰_𝑵}\mathbf{,}𝒔_{𝑱_𝑴}`$

##### B2.1 If $`𝑴\mathbf{=}\mathrm{𝟏}`$,

then $`s_{J_M}=\sigma _{J_M}`$, and we only need to generate $`s_{I_N}`$ with a density proportional to

$$\frac{1}{\sqrt{\lambda (s_{I_N},\sigma _{I_N},\mathrm{\Sigma }_{I_{N1}})}}=\frac{d}{ds_{I_N}}\mathrm{log}\left(s_{I_N}\sigma _{I_N}\mathrm{\Sigma }_{I_{N1}}+\sqrt{\lambda (s_{I_N},\sigma _{I_N},\mathrm{\Sigma }_{I_{N1}})}\right).$$

This density cancels the corresponding factor in the total antenna density. Something similar can be done if $`N=1`$.

##### B2.2 If $`𝑵\mathbf{>}\mathrm{𝟏}`$ and $`𝑴\mathbf{>}\mathrm{𝟏}`$,

then both $`s_{I_N}`$ and $`s_{J_M}`$ have to be generated, in the region where $`\sqrt{s_n}\sqrt{s_{I_N}}\sqrt{s_{J_M}}>0`$. This is far more complicated than the previous case, and we restrict ourselves to a density with a denominator proportional to $`(s_{I_N}\mathrm{\Sigma }_{I_N})(s_{J_M}\mathrm{\Sigma }_{J_M})`$. Because $`\sigma _{I_N}\mathrm{\Sigma }_{I_{N1}}`$ and $`\sigma _{J_M}\mathrm{\Sigma }_{J_{M1}}`$ may be considered small, it still cancels the factor

$$\sqrt{\lambda (s_{I_N},\sigma _{I_N},\mathrm{\Sigma }_{I_{N1}})\lambda (s_{J_M},\sigma _{J_M},\mathrm{\Sigma }_{J_{M1}})}$$

in (29). We shall write $`s_1,s_2`$ instead of $`s_{I_N},s_{J_M}`$ from now on, and denote $`m:=\sqrt{s_n}`$,

$$m_1^2:=\mathrm{\Sigma }_{I_N},c_1^2:=\underset{K,L=1}{\overset{N}{}}\sigma _{I_K,I_L}\text{and}m_2^2:=\mathrm{\Sigma }_{J_M},c_2^2:=\underset{K,L=1}{\overset{M}{}}\sigma _{J_K,J_L}.$$

We choose first to generate $`s_1`$, and then $`s_2`$, so that the integral, corresponding with the generation, becomes

$$_{c_1^2}^{(mm_2)^2}\frac{ds_1}{s_1m_1^2}_{c_2^2}^{(m\sqrt{s_1})^2}\frac{ds_2}{s_2m_2^2}.$$

The $`s_2`$-integral is simple, and shows that $`s_2`$ can easily be obtained by inversion. After integration over $`s_2`$, the $`s_1`$-integral becomes

$$_{c_1^2}^{(mm_2)^2}\frac{ds_1}{s_1m_1^2}\left[\mathrm{log}\left((m\sqrt{s_1})^2m_2^2\right)\mathrm{log}\left(c_2^2m_2^2\right)\right].$$

The $`s_1`$-variable distributed following this integrand can be obtained with high efficiency by rejection from the density proportional to $`1/(s_1m_1^2)`$. The total integral can be calculated, and is given by

$$\begin{array}{c}\underset{\rho ,ϵ=\pm 1}{}\left[\mathrm{Li}\left(\frac{x_1+m_1}{m+\rho m_2+ϵm_1}\right)+\mathrm{log}\left(\frac{m+\rho m_2x_1}{m}\right)\mathrm{log}\left(\frac{x_1+m_1}{m+\rho m_2+ϵm_1}\right)\right]_{x_1=c_1}^{mm_2}\hfill \\ \hfill \mathrm{log}\left(\frac{c_2^2m_2^2}{m^2}\right)\mathrm{log}\left(\frac{(mm_2)^2m_1^2}{c_1^2m_1^2}\right),\end{array}$$
(30)
where

$$\mathrm{Li}(x):=\underset{n=1}{\overset{\mathrm{}}{}}\frac{(1x)^n}{n^2}.$$

### C SARGE

In this Appendix, we give a short review about SARGE and present the adaptations applied in SARGE.n. We need to start with the establishment of some notation.

$`_p`$ is a Lorentz transformation that boosts momentum $`p`$ to $`(\sqrt{p^2},0,0,0)`$. $`_p`$ is a Lorentz transformation that rotates momentum $`p`$ to $`(p^0,0,0,|\stackrel{}{p}|)`$. The standard representation of a unit vector in terms of parameters $`z[1,1]`$ and $`\phi [0,2\pi ]`$ is denoted by

$$\widehat{n}(z,\phi ):=(\sqrt{1z^2}\mathrm{sin}\phi ,\sqrt{1z^2}\mathrm{cos}\phi ,z).$$

$`𝐏_n[1,1]^n`$ is the $`n`$-dimensional polytope, which consists of the support of the indicator function

$$\vartheta _{𝐏_n}(x_1,x_2,\mathrm{},x_n):=\underset{i=1}{\overset{n}{}}\mathrm{\Theta }(1|x_i|)\underset{j,k=1}{\overset{n}{}}\mathrm{\Theta }(1|x_jx_k|).$$

The algorithm for the generation of $`n`$ final-state momenta with a single antenna structure without initial-state momenta as given in  is

###### Algorithm C.1

1. generate two massless momenta $`q_1,q_n`$, back-to-back;
2. generate $`(x_1,\mathrm{},x_{2n4})𝐏_{2n4}`$ and $`(\phi _2,\mathrm{},\phi _{n1})[0,2\pi ]^{n2}`$, all uniformly distributed;
3. for $`i=2,\mathrm{},n2`$ construct $`q_i`$ following (with $`x_0:=0`$)

$$b_1q_{i1}+q_n,b_2_{b_1}q_{i1},$$

$$\xi _1e^{(x_{2i3}x_{2i4})\mathrm{log}\xi _\mathrm{m}},\xi _2e^{(x_{2i2}x_{2i4})\mathrm{log}\xi _\mathrm{m}},$$

$$v^0\sqrt{\frac{(q_{i1}q_n)}{2}}(\xi _2+\xi _1),z\frac{\xi _2\xi _1}{\xi _2+\xi _1},\stackrel{}{v}q^0_{b_2}^1\widehat{n}(z,\phi _i),$$

$$q_i_{b_1}^1v;$$
4. put $`p_iu_Qq_i`$ for $`i=1,\mathrm{},n`$, where $`Q=_{i=1}^nq_i`$ and $`u=\sqrt{s/Q^2}`$ .

The density with which the generated momenta are distributed is then given by

$$d\mathrm{\Phi }_n(P;p_1,\mathrm{},p_n)\frac{s^2}{2\pi ^{n2}}g_n(p_1,p_2,\mathrm{},p_n)A_n(p_1,p_2,\mathrm{},p_n),$$
(31)
with $`P=(\sqrt{s},0,0,0)`$, where $`A_n(p_1,p_2,\mathrm{},p_n)`$ is defined in (2), and where

$$g_n(p_1,p_2,\mathrm{},p_n):=\frac{\vartheta _{𝐏_{2n4}}(x_1,x_2,\mathrm{},x_{2n4})}{(2n3)(\mathrm{log}\xi _\mathrm{m})^{2n4}},$$

with

$$x_{2i3}=\frac{1}{\mathrm{log}\xi _\mathrm{m}}\mathrm{log}\frac{(p_{i1}p_i)}{(p_1p_n)}\text{and}x_{2i4}=\frac{1}{\mathrm{log}\xi _\mathrm{m}}\mathrm{log}\frac{(p_ip_n)}{(p_1p_n)},i=2,\mathrm{},n2.$$

The density obtained before step 4 of Algorithm C.1 is invariant under simultaneous scaling of all momenta, but the sum of the momenta is not at rest. In order to achieve this, the scaling symmetry has to be broken. In order to include the initial-state momenta in the density, more symmetries have to be broken: the symmetry under cyclic permutations of the momenta, a part of the simultaneous rotation symmetry of the momenta, and finally the symmetry under the reflection permutation $`(1,2,\mathrm{},n)(n,n1,\mathrm{},1)`$.

Let us denote $`k=kmodn`$. The cyclic symmetry is broken through

###### Algorithm C.2

1. choose a $`k\{0,2,\mathrm{},n1\}`$ with relative probability $`(p_kp_{k+1})`$ ;
2. put $`\{p_1,p_2,\mathrm{},p_n\}\{p_{1+k},p_{2+k},\mathrm{},p_{n+k}\}`$ .

As a result, the function $`A_n`$ in the density (31) is replaced by

$$\frac{n}{(p_np_1)+(p_1p_2)+\mathrm{}+(p_{n1}p_n)}\times (p_np_1)A_n(p_1,\mathrm{},p_n).$$

The factor $`(p_np_1)`$ in the denominator of $`A_n`$ is replaced by the average of all factors. Next, part of the rotation symmetry and the symmetry under the reflection permutation are broken through

###### Algorithm C.3

1. generate $`\phi [0,2\pi ]`$ uniformly distributed, and $`z[1,1c]`$ following $`1/(1z)`$ ;
2. rotate all momenta such that $`\stackrel{}{p}_1`$ lies along $`\widehat{n}(z,\phi )`$ .

If we denote $`q_1:=\frac{1}{2}\sqrt{s}(1,0,0,1)`$, this leads to an extra factor

$$\frac{1}{2\pi \mathrm{log}(2/c)}\times \frac{\frac{1}{2}\sqrt{s}p_1^0}{(q_1p_1)}$$

in the density. The cut-off $`c`$ can be taken equal to $`s_0/(\sqrt{s}p_1^0)`$, where $`s_0`$ should be a cut-off on the invariant mass $`(q_1+p_1)^2`$. Notice that $`p_1^0=(Pp_1)/\sqrt{s}`$, so that the new factor in the density becomes

$$\frac{1}{4\pi \mathrm{log}(\mathrm{\hspace{0.17em}2}(Pp_1)/s_0)}\times \frac{(Pp_1)}{(q_1p_1)}.$$

In order to include a factor proportional to $`(p_nq_2)`$ in the denominator of the density, where $`q_2:=\frac{1}{2}\sqrt{s}(1,0,0,1)`$, Algorithm C.3 can be preceded by

###### Algorithm C.4

1. choose with equal probabilities whether to rotate all momenta such that $`\stackrel{}{p}_1`$ lies along $`\widehat{n}(z,\phi )`$, or that $`\stackrel{}{p}_n`$ lies along $`\widehat{n}(z,\phi )`$.

The total density (31) becomes then such that $`A_n`$ is replaced by

$$\frac{n}{4\pi }\times \frac{\frac{(Pp_1)(p_nq_2)}{\mathrm{log}(\mathrm{\hspace{0.17em}2}(Pp_1)/s_0)}+\frac{(q_1p_1)(Pp_n)}{\mathrm{log}(\mathrm{\hspace{0.17em}2}(Pp_n)/s_0)}}{(p_np_1)+(p_1p_2)+\mathrm{}+(p_{n1}p_n)}\times A_{n+2}(q_1,p_1,p_2,\mathrm{},p_{n1},p_n,q_2).$$

In the end, we want to obtain all permutations in the momenta $`\{p_1,p_2,\mathrm{},p_n,q_2\}`$ of the above density. Some of them can be obtained by relabeling, and for the others, which change the position of $`q_2`$, we can do the following. If the cyclic permutation in Algorithm C.2 is chosen with relative probability $`(p_kp_{k+1})(p_{k+i}p_{k+1+i})`$, then the factor $`(p_np_1)(p_ip_{i+1})`$ in the denominator of $`A_n`$ is replaced by the average over its cyclic permutations. Then we can choose with equal relative probabilities whether to rotate all final-state momenta such that $`\stackrel{}{p}_1`$ lies along $`\widehat{n}`$, $`\stackrel{}{p}_i`$ lies along $`\widehat{n}`$, $`\stackrel{}{p}_{i+1}`$ lies along $`\widehat{n}`$, or $`\stackrel{}{p}_n`$ lies along $`\widehat{n}`$. This will lead to density that is proportional to

$$A_{n+2}(q_1,p_1,p_2,\mathrm{},p_{i1},p_i,q_2,p_{i+1},p_{i+2},\mathrm{},p_{n1},p_n).$$

### D Permutations

For the analyses of the the multi-channel procedure over the antenna densities for the different permutations of the momenta, we need an enumeration of the permutations. In other words, we need a mapping

$$\{1,2,\mathrm{},n!\}\mathrm{𝑆𝑦𝑚}_n.$$

During a computation, one could, of course, just store an enumeration in the memory of the computer, but this costs an amount of memory of $`𝒪(n!)`$. Every algorithm that delivers all permutations supplies a mapping a priori, but this algorithm has an a priori computational complexity of $`𝒪(n!)`$.

We propose the following algorithm as a solution. Any number $`l\{1,2,\mathrm{},n!\}`$ can uniquely be written in the basis of the lower factorials: $`l=1+l_2+l_32!+l_43!+\mathrm{}+l_n(n1)!`$, where $`l_k\{0,1,\mathrm{},k1\}`$. Let $`\gamma _i^{(k)}`$ denote the $`i`$-fold cyclic permutation of the first $`k`$ elements of $`(1,2,\mathrm{},n)`$, for example $`\gamma _2^{(3)}(1,2,3,4)=(2,3,1,4)`$. Then any permutation of $`(1,2,\mathrm{},n)`$ can be written as $`\gamma _{l_2}^{(2)}\gamma _{l_3}^{(3)}\mathrm{}\gamma _{l_n}^{(n)}(1,2,\mathrm{},n)`$. This leads to a mapping of the required kind, which has a complexity of $`𝒪(n^2)`$. In a loop delivering all permutations, it can easily be reduced to complexity $`𝒪(n)`$ at the cost of an amount of memory of $`𝒪(n)`$.