Logarithmic corrections for jet production at the LHC

E. P. Byrne1$\star$

1 Higgs Centre for Theoretical Physics, University of Edinburgh,

Peter Guthrie Tait Road, Edinburgh, EH9 3FD, UK

*<EMAIL_ADDRESS>
50th International Symposium on Multiparticle Dynamics
(ISMD2021)
12-16 July 2021 10.21468/SciPostPhysProc.?

## Abstract

Several important processes and analyses at the LHC are sensitive to higher-
order perturbative corrections beyond what can currently be calculated at
fixed order. One important class of logarithmic corrections are those which
appear when the centre-of-mass energy of a QCD collision is much larger than
the transverse momenta of the observed jets. We describe the High Energy Jets
(HEJ) framework, which includes the dominant high-energy logarithms to provide
all-order predictions for several LHC processes including Higgs, $W$, or $Z$
boson production in association with jets. We will summarise some recent
developments, in particular the first matching of HEJ to a NLO calculation.

##  1 Introduction

Standard Model predictions for the LHC are typically produced using general
purpose Monte Carlo event generators which describe many aspects of a proton-
proton collision, from PDFs to hadronisation. Here we restrict our attention
to the description of the hard matrix element which describes the scattering
between one parton from each proton. The hard matrix element is typically
computed via fixed-order perturbation theory. For QCD amplitudes, the coupling
$\alpha_{s}$ decreases with increasing centre of mass energy. However, there
are regions of phase space where an expansion in $\alpha_{s}$ is unstable. For
example, consider the cross section for 2$\to$2 scattering of partons in the
region where the Mandelstam variable $s$ is much larger than $|t|$. At order
$\alpha_{s}^{n+2}$ there are terms which are of the form $\log^{n}(s/|t|)$. At
energies accessible by the LHC the quantity
$\left(\alpha_{s}\log(s/|t|)\right)^{n}$ can be of order one. Logarithms of
this type are called high-energy logarithms and must be summed to all orders
to ensure stability of our perturbative predictions. This motivates an
alternative classification of accuracy, where NmLL refers to all terms in the
perturbative series of the form $\alpha_{s}^{n+2}\log(s/|t|)^{n-m}$.

The LL contribution to QCD amplitudes, to all orders in $\alpha_{s}$, takes a
simple factorised form [1, 2, 3]. Rather than describing the hard matrix
element of a collision at the LHC via fixed-order perturbation theory we can
instead take as our starting point matrix elements which are LL accurate. High
Energy Jets (HEJ) is a framework that builds on this LL accuracy. However, HEJ
utilises Monte Carlo integration for phase space integration, which means
minimal approximations to the amplitude need to be made, and no formal
approximations to phase space. This in turn allows HEJ to make all-order
predictions for not only asymptotically large energies, but also at the scales
of energy accessible by the LHC and for arbitrary cuts and analyses.

In Section 2 of we review the simple factorised form that amplitudes take in
the high-energy limit. In Section 3 we show how the HEJ framework builds upon
this factorised picture to obtain an improved LL accurate description, and we
discuss NLL improvements to this framework. In Section 4 we discuss how this
framework can be matched to existing fixed-order calculations. We review a
recent application of this framework to the production of a $W$ boson in
association with at least two jets before stating our conclusions.

##  2 Factorisation of amplitudes in the high-energy limit

To LL accuracy, the amplitudes for $2\to n$ processes in QCD (with $n\geq 2$)
take a simple factorised form [1, 2, 3]. These amplitudes become exact in the
so-called Multi-Regge-Kinematic (MRK) limit, where for momenta $p_{a},p_{b}\to
p_{1},\cdots,p_{n}$ the outgoing momenta are strictly ordered in rapidity
$y_{i}$, while the transverse momenta $\mathbf{p}_{i}$ are of the same finite
magnitude:

$\displaystyle\begin{split}&y_{i}\gg y_{i+1}\quad\forall
i\in\\{1,\cdots,n-1\\},\quad|\mathbf{p}_{i}|\approx|\mathbf{p}_{j}|\quad\forall
i,j\in\\{1,\cdots,n\\}.\end{split}$ (1)

As an illustration of the simplicity of these amplitudes, we give the squared
amplitude for the scattering of distinguishable quarks $q$ and $Q$ to $n$
partons:

$\displaystyle\begin{split}\left|\mathcal{M}^{\mathrm{LL\ MRK}}_{qQ\to
n}\right|^{2}&=\frac{4s^{2}}{(N_{C}^{2}-1)}\left(\frac{g_{s}^{2}C_{F}}{|\mathbf{p}_{1}|^{2}}\right)\times\left(\prod_{i=2}^{n-1}\frac{4g_{s}^{2}C_{A}}{|\mathbf{p}_{i}|^{2}}\right)\times\left(\prod_{i=1}^{n-1}e^{2\alpha(\mathbf{q}_{i})\left(y_{j+1}-y_{j}\right)}\right)\times\left(\frac{g_{s}^{2}C_{F}}{|\mathbf{p}_{n}|^{2}}\right).\end{split}$
(2)

This expression makes full use of the simplifications allowed by the MRK
region eq. 1. The only real-emission corrections which are relevant to LL
order are emissions of gluons within the rapidity span of the outgoing quarks.
Virtual corrections exponentiate in this limit and are described by the Regge
trajectory of the associated $t$-channel momentum $\mathbf{q}_{i}$,

$\displaystyle\begin{split}\alpha(\mathbf{q}_{i})=-g_{s}^{2}C_{A}\frac{\Gamma(1-\epsilon)}{(4\pi)^{2+\epsilon}}\frac{1}{\epsilon}\left(\frac{\mathbf{q}_{i}^{2}}{\mu^{2}}\right)^{\epsilon}.\end{split}$
(3)

The IR singularity here can be regularised by the procedure described in [4].

While expressions such as eq. 2 are useful for probing the asymptotic limit of
QCD, the kinematic approximations are too severe to be directly applied to the
LHC. This motivates constructing a framework which maintains the logarithmic
accuracy of eq. 2 while relaxing many of the kinematic approximations used to
obtain it.

##  3 The High Energy Jets framework

###  3.1 HEJ at LL

We begin by noting that at LO, the amplitude for $qQ\to qQ$ scattering already
has the form of high-energy factorisation without taking any kinematic
approximations:

$i\mathcal{M}^{(0)}_{qQ\rightarrow qQ}=\left(ig_{s}T_{1a}^{C}j_{q\to
q}^{\mu}(p_{a},p_{1})\right)\left(\frac{-i}{t}\right)\left(ig_{s}T_{2b}^{C}j_{Q\to
Q\ \mu}(p_{b},p_{2})\right),$ (4)

where we have defined the factorised spinor current

$\displaystyle\begin{split}j_{q\to
q}^{\mu}(p_{1},p_{a})&=\bar{u}^{\lambda_{1}}(p_{1})\gamma^{\mu}u^{\lambda_{a}}(p_{a}).\end{split}$
(5)

If we use these quark currents as our factorised building blocks we can
construct a framework which is LO accurate for this process without loosing LL
accuracy [5]. The leading-order amplitude for $qg\to qg$ contains $s$\- and
$u$-channel diagrams as well as a $t$-channel diagram However in ref. [6] it
was shown that for the dominant helicity configurations, the $qg\to qg$
amplitude can be written exactly in the same form as $qQ\to qQ$ with $C_{F}$
replaced with a momentum-dependent factor $K_{g}$.

The HEJ form of LL amplitudes for pure-QCD processes is

$\displaystyle\begin{split}\left|\mathcal{M}^{\mathrm{LL\
\texttt{HEJ}}}_{f_{a}f_{b}\to f_{a}\cdot(n-2)g\cdot
f_{b}}\right|^{2}&=\frac{1}{4(N_{C}^{2}-1)}\left|j_{f_{a}\to f_{a}}\cdot
j_{f_{b}\to
f_{b}}\right|^{2}\left(g_{s}^{2}K_{f_{a}}(p_{a},p_{1})\frac{1}{t_{1}}\right)\left(g_{s}^{2}K_{f_{b}}(p_{b},p_{n})\frac{1}{t_{n-1}}\right)\\\
&\times\left(\prod_{i=1}^{n-1}e^{\omega(t_{i})\left(y_{j+1}-y_{j}\right)}\right)\times\left(\prod_{i=2}^{n-1}\frac{-g_{s}^{2}C_{A}}{t_{i}t_{i+1}}\left|V_{g}(q_{i},q_{i+q})\right|^{2}\right),\end{split}$
(6)

which should be compared to eq. (2). While all-orders virtual corrections are
described in the same manner, significantly more information about the
emission of real radiation is contained via the HEJ form of the so-called
Lipatov vertex [2, 5]

$\displaystyle\begin{split}V_{g}^{\rho}(q_{i},q_{i+1})=-(q_{i}+q_{i+1})^{\rho}&+\frac{p_{a}^{\rho}}{2}\left(\frac{q_{i}^{2}}{p_{i+1}\cdot
p_{a}}+\frac{p_{i+1}\cdot p_{b}}{p_{a}\cdot p_{b}}+\frac{p_{i+1}\cdot
p_{n}}{p_{a}\cdot p_{n}}\right)+p_{a}\leftrightarrow p_{1}\\\
&-\frac{p_{b}^{\rho}}{2}\left(\frac{q_{i+1}^{2}}{p_{i+1}\cdot
p_{b}}+\frac{p_{i+1}\cdot p_{a}}{p_{b}\cdot p_{a}}+\frac{p_{i+1}\cdot
p_{1}}{p_{b}\cdot p_{1}}\right)-p_{b}\leftrightarrow p_{n}.\end{split}$ (7)

The HEJ framework can be extended to describe the production of a $W$, $Z$, or
Higgs boson or same-sign $W$-pair production in association with jets [5, 7,
8, 9, 10].

###  3.2 NLL improvements to HEJ

The HEJ framework may be systematically improved by either matching to fixed-
order calculations (see section 4), or by increasing the logarithmic accuracy
of the all-orders amplitudes. Leading powers in $s/|t|$ in matrix elements
lead to leading logarithms after integration. Regge theory gives the scaling
of matrix elements in terms of the spin of the particles exchanged in the
$t$-channels of the planar Feynman diagrams that can be drawn for a given
process when the legs are ordered in rapidity. A flavour/momentum
configuration contributes at LL if it permits an exchange of a gluon in all
$t$-channels. Beyond that, a flavour/momentum configuration contributes at
NmLL if it permits an exchange of a gluon in all but $m$ of the $t$-channels
(see ref. [11] for further details).

The first NLL component to be included in HEJ was the description of a gluon
emission where that gluon was more extreme in rapidity than one of the
outgoing quarks [12]. In ref. [11], the necessary pieces have been calculated
to include all-order corrections to all configurations at 3$j$ and above whose
leading contribution is at NLL order. This also includes the potential
emission of a $W$ boson from these pieces. This is a gauge-invariant subset of
the full NLL correction to inclusive 2$j$ or $W+2j$ production.

##  4 NLO matching

We now describe the method used in ref. [11] to increase the fixed-order
accuracy of HEJ to NLO. We write the NLO calculation of the 2$j$ cross section
as

$\sigma^{\text{NLO}}_{2j}=f_{2j}^{(2)}\alpha_{s}^{2}+\left(f_{2j}^{(3)}+f_{3j}^{(3)}\right)\alpha_{s}^{(3)}.$
(8)

Eq. (8) is finite, after cancellation of IR singularities between the virtual
corrections and unresolved real emissions. We may write the predictions of HEJ
as an all-orders expression in $\alpha_{s}$

$\sigma^{\text{HEJ}}_{2j}=h_{2j}^{(2)}\alpha_{s}^{2}+(h_{2j}^{(3)}+h_{3j}^{(3)})\alpha_{s}^{3}+(h_{2j}^{(4)}+h_{3j}^{(4)}+h_{4j}^{(4)})\alpha_{s}^{4}+(h_{2j}^{(5)}+h_{3j}^{(5)}+h_{4j}^{(5)}+h_{5j}^{(5)})\alpha_{s}^{5}+\mathcal{O}(\alpha_{s}^{6}).$
(9)

Within HEJ [13], $h^{(n)}_{nj}$ is matched to $f^{(n)}_{nj}$ where this is
available. We may truncate the above series to $\alpha_{s}^{3}$, which we will
refer to as HEJ@NLO:

$\displaystyle\begin{split}\sigma^{\text{HEJ@NLO}}_{2j}&=h_{2j}^{(2)}\alpha_{s}^{2}+(h_{2j}^{(3)}+h_{3j}^{(3)})\alpha_{s}^{3}.\\\
\end{split}$ (10)

In order to ensure that eq. (9) agrees with the fixed-order approach up to
$\alpha_{s}^{3}$ without damaging the logarithmic accuracy of HEJ, we multiply
$\sigma^{\text{HEJ}}_{2j}$ by the ratio

$\displaystyle\frac{\sigma^{\text{NLO}}_{2j}}{\sigma^{\text{HEJ@NLO}}_{2j}}$
$\displaystyle=1+(f_{2j}^{(3)}-h_{2j}^{(3)})\prod_{n=0}^{\infty}(-1)^{n}\alpha_{s}^{(n+1)}\frac{(f_{3j}^{(3)}+h_{2j}^{(3)})^{n}}{(f_{2j}^{(2)})^{(n+1)}}.$
(11)

In the limit of $s\gg|t|$ this ratio tends to one, where in particular,
$f_{2j}^{(3)}\to h_{2j}^{(3)}$. Applying this reweighting factor to eq. (9)
gives

$\displaystyle\sigma^{\text{HEJ}}_{2j}\left(\frac{\sigma^{\text{NLO}}_{2j}}{\sigma^{\text{HEJ@NLO}}_{2j}}\right)$
$\displaystyle=f_{2j}^{(2)}\alpha_{s}^{2}+\left(f_{2j}^{(3)}+f_{3j}^{(3)}\right)\alpha_{s}^{3}+\left(h_{2j}^{(4)}+h_{3j}^{(4)}+f_{4j}^{(4)}\right)\alpha_{s}^{4}+\mathcal{O}(\alpha_{s}^{5}),$
(12)

and we see that the resulting quantity is indeed NLO accurate in all of phase
space.

So far we have discussed cross-sections, but the same argument holds for any
particular bin in a distribution of any observable in which we might be
interested. For the practical implementation of this method, we must obtain
separate histograms for NLO, HEJ and HEJ@NLO. In the study ref. [11] we used
Sherpa [14] with OpenLoops [15] to obtain pure NLO results. For HEJ@NLO we
truncate the all-orders prediction of the relevant HEJ amplitudes. We further
separate the full HEJ result into two pieces:

  * •

The LL resummed HEJ predictions for 2$j$ inclusive (in practice we observe
convergence after 6 jets) including LL resummation applied to the NLL
configurations discussed in section 3.2, plus $3j$ configurations which are
included at LO only via fixed-order matching.

  * •

Configurations at $\geq 4j$ which are included at LO only via fixed-order
matching.

This separation allows us to avoid applying NLO-matching to the fixed-order
component of the HEJ result which begins at 4$j$. Finally, we can obtain the
NLO-matched bin weight

$w_{\mathrm{\texttt{HEJ}}2\mathrm{NLO}}=w_{\mathrm{\texttt{HEJ}}(\geq 2j\
\mathrm{Resummed}+3j\
\mathrm{FO})}\frac{w_{\mathrm{\texttt{HEJ}\text{@NLO}}}}{w_{\mathrm{NLO}}}+w_{\mathrm{\texttt{HEJ}}(\geq
4j\ \mathrm{FO})}.$ (13)

Figure 1: The NLO-matching corrections for distributions presented in [16, 17]
for inclusive $W^{+}+2j$ production (blue, solid) and $W^{-}+2j$ production
(red, dashed).

Fig. 1 shows the matching corrections,
($w_{\mathrm{\texttt{HEJ}\text{@NLO}}}/w_{\mathrm{NLO}}$), for two
distributions corresponding to measurements in refs. [17, 16]. Such
$p_{\perp}$-based observables are of interest as they deviate from the
transverse MRK condition (eq. (1)). As such, they are regions where we expect
HEJ to particularly benefit from matching to NLO. We apply the matching
procedure separately for $W^{+}$ and $W^{-}$ production but observe no
significant difference. All calculations presented use a central scale of
$\mu_{r}=\mu_{f}=H_{T}/2$. The scale variation bands are obtained by varying
$\mu_{r}$ and $\mu_{f}$ by a factor of two around this central scale, while
keeping their ratio between 0.5 and 2.

The application of these matching corrections are show in fig. 2. The blue
line shows the results of the original LL HEJ framework. The green line
benefits from both the inclusion of the class of NLL improvements discussed in
sec. 3.2 and the application of the matching corrections shown in fig. 1. The
pure NLO predictions are shown in red. We see that the NLO-matched HEJ
predictions give a favourable description of the data, in particular at large
$p_{\perp}$ where NLO fixed-order and LL corrections are both significant.
Furthermore, we see that the NLO-matched HEJ predictions benefit from a
significant reduction in the scale variation band.

Figure 2: Predictions for $pp\to(W\to l\nu)+\geq 2j$ for the LHC, compared to
ATLAS data at 7 TeV [16] and 8 TeV [17]. The content of the lines is discussed
in the text.

##  5 Conclusion

In this proceeding we have introduced the High Energy Jets (HEJ) framework,
which provides all-orders predictions for jet production at the LHC that
include high-energy effects and restore stability to the perturbative
expansion. As an example, we have applied this to the process
$pp\to(W\to\ell\nu)+\geq 2j$. We have seen that the combined effect of two
improvements to the HEJ framework, namely the inclusion of a gauge invariant
class of NLL improvements and the matching to a NLO calculation, leads to
close agreement with data. This highlights the importance of utilising the
best of fixed-order and all-order calculations in order to provide precise and
stable predictions for the LHC.

## Acknowledgements

We thank the organisers for an interesting and engaging workshop. We would
further like to thank the other members of the HEJ collaboration for
discussions during this work. We gratefully acknowledge funding received from
the ERC Starting Grant 715049 "QCDforfuture".

## References

  * [1] V. S. Fadin, E. A. Kuraev and L. N. Lipatov,  _On the Pomeranchuk Singularity in Asymptotically Free Theories_ ,  Phys. Lett. B 60, 50 (1975),  10.1016/0370-2693(75)90524-9.
  * [2] L. N. Lipatov,  _Reggeization of the Vector Meson and the Vacuum Singularity in Nonabelian Gauge Theories_ ,  Sov. J. Nucl. Phys. 23, 338 (1976).
  * [3] E. A. Kuraev, L. N. Lipatov and V. S. Fadin,  _The Pomeranchuk Singularity in Nonabelian Gauge Theories_ ,  Sov. Phys. JETP 45, 199 (1977).
  * [4] J. R. Andersen and J. M. Smillie,  _Multiple Jets at the LHC with High Energy Jets_ ,  JHEP 06, 010 (2011),  10.1007/JHEP06(2011)010,  1101.5394.
  * [5] J. R. Andersen and J. M. Smillie,  _Constructing All-Order Corrections to Multi-Jet Rates_ ,  JHEP 01, 039 (2010),  10.1007/JHEP01(2010)039,  0908.2786.
  * [6] J. R. Andersen and J. M. Smillie,  _The Factorisation of the t-channel Pole in Quark-Gluon Scattering_ ,  Phys. Rev. D 81, 114021 (2010),  10.1103/PhysRevD.81.114021,  0910.5113.
  * [7] J. R. Andersen, T. Hapola and J. M. Smillie,  _W Plus Multiple Jets at the LHC with High Energy Jets_ ,  JHEP 09, 047 (2012),  10.1007/JHEP09(2012)047,  1206.6763.
  * [8] J. R. Andersen, J. J. Medley and J. M. Smillie,  _$Z/\gamma$ * plus multiple hard jets in high energy collisions_,  JHEP 05, 136 (2016),  10.1007/JHEP05(2016)136,  1603.05460.
  * [9] J. R. Andersen, J. D. Cockburn, M. Heil, A. Maier and J. M. Smillie,  _Finite Quark-Mass Effects in Higgs Boson Production with Dijets at Large Energies_ ,  JHEP 04, 127 (2019),  10.1007/JHEP04(2019)127,  1812.08072.
  * [10] J. R. Andersen, B. Ducloué, C. Elrick, A. Maier, G. Nail and J. M. Smillie,  _Logarithmic corrections to the QCD component of same-sign W-pair production for VBS studies_ (2021),  2107.06818.
  * [11] J. R. Andersen, J. A. Black, H. M. Brooks, E. P. Byrne, A. Maier and J. M. Smillie,  _Combined subleading high-energy logarithms and NLO accuracy for W production in association with multiple jets_ ,  JHEP 04, 105 (2021),  10.1007/JHEP04(2021)105,  2012.10310.
  * [12] J. R. Andersen, T. Hapola, A. Maier and J. M. Smillie,  _Higgs Boson Plus Dijets: Higher Order Corrections_ ,  JHEP 09, 065 (2017),  10.1007/JHEP09(2017)065,  1706.01002.
  * [13] J. R. Andersen, T. Hapola, M. Heil, A. Maier and J. Smillie,  _HEJ 2: High Energy Resummation for Hadron Colliders_ ,  Comput.Phys.Commun. (2019),  10.1016/j.cpc.2019.06.022,  1902.08430.
  * [14] E. Bothmann _et al._ ,  _Event Generation with Sherpa 2.2_ ,  SciPost Phys. 7(3), 034 (2019),  10.21468/SciPostPhys.7.3.034,  1905.09127.
  * [15] F. Buccioni, S. Pozzorini and M. Zoller,  _On-the-fly reduction of open loops_ ,  Eur. Phys. J. C 78(1), 70 (2018),  10.1140/epjc/s10052-018-5562-1,  1710.11452.
  * [16] G. Aad _et al._ ,  _Measurements of the W production cross sections in association with jets with the ATLAS detector_ ,  Eur. Phys. J. C75(2), 82 (2015),  10.1140/epjc/s10052-015-3262-7,  1409.8639.
  * [17] M. Aaboud _et al._ ,  _Measurement of differential cross sections and $W^{+}/W^{-}$ cross-section ratios for $W$ boson production in association with jets at $\sqrt{s}=8$ TeV with the ATLAS detector_,  JHEP 05, 077 (2018),  10.1007/JHEP05(2018)077,  1711.03296.


# Driven-dissipative Bose-Einstein condensation of photons in a vertical-
cavity surface-emitting laser

M. Pieczarka<EMAIL_ADDRESS>Department of Experimental Physics,
Faculty of Fundamental Problems of Technology, Wrocław University of Science
and Technology, Wyb. Wyspiańskiego 27, 50-370 Wrocław, Poland  M. Gębski
Institute of Physics, Lodz University of Technology, ul. Wólczańska 217/221,
93-005 Łódź, Poland  A. N. Piasecka  Department of Experimental Physics,
Faculty of Fundamental Problems of Technology, Wrocław University of Science
and Technology, Wyb. Wyspiańskiego 27, 50-370 Wrocław, Poland  J. A. Lott
Institute of Solid State Physics and Center of Nanophotonics, Technical
University Berlin, Hardenbergstraße 36, 10623 Berlin, Germany  A. Pelster
Department of Physics and Research Center OPTIMAS, Rheinland-Pfälzische
Technische Universität Kaiserslautern-Landau, Erwin Schrödinger Straße 46,
67663 Kaiserslautern, Germany  M. Wasiak  Institute of Physics, Lodz
University of Technology, ul. Wólczańska 217/221, 93-005 Łódź, Poland  T.
Czyszanowski  Institute of Physics, Lodz University of Technology, ul.
Wólczańska 217/221, 93-005 Łódź, Poland

††preprint: VCSEL-BEC

Many bosons can occupy a single quantum state without a limit, which is
described by quantum-mechanical Bose-Einstein statistics and allows the
formation of a Bose-Einstein condensate at low temperatures and high particle
densities. Photons, historically the first considered bosonic gas, were among
the last to show this phenomenon, which was observed in rhodamine-filled
microlaser cavities or doped fiber cavities. These more recent findings have
raised the natural question as to whether condensation is common in laser
systems, potentially implying its technological application. Here, we show the
Bose-Einstein condensation of photons in a semiconductor vertical-cavity
surface-emitting laser with positive cavity-gain energy detuning. We observed
a Bose-Einstein condensate in the ground mode at the critical phase-space
density. Experimental results follow the equation of state for a two-
dimensional gas of bosons in thermal equilibrium, although the extracted
spectral temperatures were lower than those of the device. This is interpreted
as originating from the driven-dissipative nature of the condensate being not
in full thermal equilibrium with the device. In contrast, nonequilibrium
lasing action is observed in higher-order modes in a negatively detuned
device. Our work enables the potential exploration of superfluid physics of
interacting photons mediated by semiconductor optical non-linearities. It also
shows great promise for single-mode high-power emission from a large aperture
device.

At the beginning of the 20th century, Albert Einstein extended the statistical
theory of Satyendra Nath Bose to describe massive particles and made a
pioneering prediction of the Bose-Einstein Condensate (BEC) below a critical
temperature [1]. BEC is characterized by both saturation of occupation in the
excited states and condensation in the ground state of the system [2]. Seventy
years after its theoretical prediction, this macroscopic quantum phenomenon
was first directly observed in dilute clouds of atomic gases at temperatures
close to absolute zero [3, 4]. The reason for such a small critical
temperature is that it turns out to be inversely proportional to the mass of
the boson, thereby a heavy particle gas must be extremely cold to reach the
transition point. However, considering the mass as the parameter of the energy
dispersion and not as the gravitational one, one can find a bosonic
quasiparticle described with a dispersion of large curvature, and hence a
quite small effective mass, which enables condensation at elevated
temperatures. This concept has been realized in a variety of bosonic
quasiparticle systems, such as magnons [5], excitons [6, 7, 8], and plasmons
[9], as well as hybrid excitations of strongly coupled systems of exciton and
photons, namely cavity polaritons [10, 11].

Figure 1: Basic properties of a VCSEL. (a) Scheme of the investigated VCSEL
devices with all main components indicated by arrows. (b) Simplified picture
of the conduction and valence subbands confined in the QWs expressed in the
in-plane wavevector (left). The occupations of the conduction band $f_{c}$ and
valence band states $f_{v}$ expressed with Fermi-Dirac distributions of
different quasi-Fermi levels $\mu_{c}$ and $\mu_{v}$, respectively. (c) Output
power-current-voltage (LIV) characteristics of the BEC device.

Photons, on the other hand, have been out of the picture for many years
because they represent a massless gas with linear energy dispersion and a
trivial, null ground state. In principle, the number of particles is not
conserved, i.e. in a blackbody radiation model in thermal equilibrium the
chemical potential vanishes, and therefore condensation cannot occur.
Nevertheless, after many years of research in laser and atomic BEC physics,
many analogies between them have been drawn, yielding ultimately a more
detailed understanding of these two worlds. Eventually, a system that meets
all the requirements of an equilibrium photon BEC was obtained in a laboratory
tabletop system of a cavity filled with rhodamine solution [12]. Remarkably,
this system has undoubtedly demonstrated many textbook effects of a non-
interacting condensate of bosons, reaching from thermodynamic and caloric
properties [13, 14] to quantum-statistical ones [15, 16]. Moreover, the driven
dissipative nature of this system beyond equilibrium has also been
demonstrated [17], and the phase boundaries between photon BEC and non-
equilibrium lasing have been extensively investigated in this system [18, 19].
However, rhodamine-based photon BECs are limited by their weak and slow
thermo-optical nonlinearity, having so far prevented the observation of static
or dynamic superfluid effects. Pioneering observations have stimulated the
search for BEC conditions in other laser systems, such as fiber cavities [20]
and semiconductor lasers [21, 22, 23], to enable its true technological
application outside of the laboratory environment and to find a material
system with non-negligible and fast non-linearities.

Here, we demonstrate a photon BEC in a well-established semiconductor device,
a large aperture electrically driven vertical-cavity surface-emitting laser
(VCSEL) at room temperature. By testing devices with different energy
detunings between the cavity fundamental mode $\varepsilon_{0}$ and the
quantum well (QW) fundamental transition $\varepsilon_{\rm QW}$, defined as
$\Delta=\varepsilon_{0}-\varepsilon_{\rm QW}$, we found a homogeneous BEC of
photons with a thermalized spectrum. This occurred for $\Delta>0$ and standard
non-equilibrium laser operation at higher-order modes in another device of the
same geometry but with $\Delta<0$. In the BEC regime, we found that the
photonic gas thermalizes to temperatures below the temperature of the VCSEL,
suggesting that it is not fully equilibrated with the optical gain medium.
Nevertheless, the extracted temperatures and chemical potentials allowed us to
experimentally determine the equation of state (EOS), which follows the
behavior of a 2D Bose gas in thermal equilibrium.

The device under study is an epitaxially-grown, oxide-confined 980 nm VCSEL
with a large 23-µm diameter aperture. The VCSEL is designed for simultaneous
high bandwidth, high optical output power, and moderate to high wall plug
efficiency [24, 25] (see Methods for details). We drive our semiconductor
device at room temperature with direct current by applying a constant voltage
across the laser diode; see Fig. 1(a). This sets the non-equilibrium
distribution of carriers in the QW region, as the separation of the quasi-
Fermi levels for electrons in the conduction band states $\mu_{c}$ and holes
in the valence band states $\mu_{v}$ is proportional to the applied voltage.
Due to the sub-picosecond relaxation of carriers within the bands [26], the
electrons and holes are in equilibrium with the device; hence both gases can
be described with separate Fermi distributions with different quasi-Fermi
levels setting the occupation in both bands, Fig. 1(b) [27].

Figure 2: Momentum-space and real-space spectra of the BEC and a laser device.
(a) Scheme of the experimental setup used for momentum-space imaging. The back
focal plane of the microscope objective is imaged onto the entrance slit of
the monochromator and then it is dispersed to the CCD camera enabling a
probing of the spectral information at the center cut of the momentum space.
(b) Momentum-space spectrum of the BEC device below and (c) above the
condensation threshold showing the narrowing from thermal distribution to the
ground state $k_{\parallel}\approx 0$. (d) Momentum-space spectrum of the
lasing device in the higher-order mode above the threshold. (e),(f) Real-space
spectra of the BEC device below and above the threshold showing homogeneity of
the gas. (g) Real-space spectrum of the lasing device that presents the
domination of the higher-order mode. All color scales are logarithmic. Insets
in (b-g) represent the normalized energy-integrated spectra.

Let us assess the essential conditions to obtain the photon BEC in a VCSEL. In
electrically driven semiconductors, excited electrons and holes can recombine,
emitting photons. Thus, the condition of chemical equilibrium can be
established if the chemical potential of photons is equal to
$\mu=\mu_{c}-\mu_{v}$, in close analogy to a photochemical reaction [28]. This
well-defined chemical potential of a photon gas is essential to obtain a BEC
at equilibrium. Another key ingredient is the detailed balance condition
between emission and absorption, which is common in organic laser dyes [12].
This condition is also met in semiconductors, where the ratio between emission
and absorption rates is expressed as the van Roosbroeck-Schockley relation
$R_{\rm abs}(\varepsilon)/R_{\rm
em}(\varepsilon)=\exp(\frac{\varepsilon-\mu}{k_{B}T})$ (see Methods) [27, 29].
Hence, the thermalization of light occurs after a few cycles of spontaneous
emission and absorption events of the photon before escaping the cavity
through the top mirror. Such an energy exchange with the active medium enables
the photon gas to establish both a chemical potential and a temperature.
Eventually, it leads to a modified Bose-Einstein (BE) distribution of photons,
which can be derived from the laser rate equations (see Methods):

$N(\varepsilon)=\frac{1}{\exp(\frac{\varepsilon-\mu}{k_{B}T})-1+\Gamma(\varepsilon)}\,.$
(1)

Here the correction parameter $\Gamma(\varepsilon)=\gamma(\varepsilon)/R_{\rm
em}(\varepsilon)$ represents the ratio of photon decay rate from the passive
cavity $\gamma$ and the spontaneous emission rate $R_{\rm em}$ to the photon
mode at given energy $\varepsilon$. Consequently, this correction parameter
can be treated as a measure of the thermalization degree. It is expected to be
small if many photon emission-absorption cycles occur before escaping the
cavity. In this limit, Equation (1) approaches the Bose-Einstein distribution.
Based on our numerical modeling and experimental measurements, we estimated
this ratio for the fundamental mode and found it to be
$\Gamma(\varepsilon_{0})\approx 0.008$ (see Methods and Supplementary
Information for details), ensuring that we obtained, indeed, a thermalized
photon gas in our system.

According to standard semiconductor laser theory, the Bernard-Durraffourg
condition [30], essential for non-equilibrium lasing, is met when the value of
the chemical potential exceeds the energy of an optical mode
$\mu>\varepsilon$. This creates a positive optical gain at this energy [27],
so thermalization is expected to dominate below reaching this limit.
Therefore, we probed devices with different cavity-QW gain detunings $\Delta$.
To this end, we used the side effect of epitaxial growth that the resulting
layers are not homogeneous throughout the entire wafer and have a tendency to
become thinner towards the edge [31, 32, 33]. Phenomenologically, this affects
the cavity energy shifts more than the spectral shifts of the gain. Thus,
close to the center of the two-inch wafer we probe the device with $\Delta<0$,
which is the standard designed detuning for high-power and high-temperature
lasing operation, which we denote as the lasing device. In contrast, close to
the wafer edge, the detuning becomes positive, and the device is expected to
operate in the thermalized BEC regime, which we denote as the BEC device.
Although we cannot directly measure the precise value of the detuning, we
observed a stark contrast in the performance of devices from these distant
positions on the sample, supporting our assumptions. The electrical and total
output power characteristics on the driving current of the BEC device are
shown in Fig. 1(c) (the results of the lasing device are provided in the
Supplementary Information). It shows all the standard features of a laser, the
electrical characteristics of a diode, and the emission threshold current
$I_{\rm th}$, however, the device is characterized by significant spontaneous
emission below $I_{\rm th}$. Therefore, the information contained in the
spectral characteristics of the device must be examined to distinguish between
BEC and lasing.

Figure 3: Experimental energy distributions and thermodynamic quantities (a)
Solid lines represent energy distributions extracted from the momentum spectra
for different driving currents. (b) The same data is represented in
logarithmic form (see text). In (a),(b) the energy scale is expressed with
respect to the energy in the ground mode. The dashed lines are the fits of the
BE distribution to the experimental data. The error bars (95% confidence
interval) are represented as shaded regions. (c) Population of the ground
state ($N_{0}$) and excited states ($N_{T}$) extracted from the experimental
spectra. The dashed line is the linear fit above the condensation threshold to
calculate the critical density ($N_{C}$). The inset shows the zoom-in into the
low-number region of the main plot. (d) Thermodynamic quantities, effective
chemical potential $\mu_{\text{eff}}$ and temperature $T_{\textbf{eff}}$
extracted from fitting the experimental distributions, as a function of
driving current. $T=293~{}\text{K}$ is the temperature of the heat sink. Error
bars in (c),(d) represent the 95% confindence intervals.

We performed an analysis of the VCSEL spectral features, especially the
distribution of occupations in energy states. The investigated devices have
large electrical apertures of 23 µm in diameter, resulting in a quasi-
continuum of transversal modes (or in-plane energy states). Thus, photons in
the resonator can be described by a parabolic dispersion in the in-plane
direction $\varepsilon_{k}=\varepsilon_{0}+\frac{\hbar^{2}k^{2}}{2m_{\rm ph}}$
with an effective mass $m_{\rm ph}$. We employed the back focal plane (Fourier
space or far-field) imaging technique to directly access the momentum
dispersion, as is schematically shown in Fig. 2(a). The image is directed onto
the monochromator slit allowing us to analyze the momentum dispersions
spectrally. The momentum dispersion below $I<I_{\rm th}$ is presented in Fig.
2(b). One observes a thermalized distribution in momentum space, following the
expected parabolic dispersion. The most distinguishing feature is observed
above the threshold $I>I_{\rm th}$, see Fig. 2(b), where the ground mode at
$k_{\parallel}=0$ dominates the spectrum. This is unusual behavior for such a
large aperture resonator, as lasing in higher-order modes is commonplace. We
obtain this standard behavior in our lasing device with negative detuning,
where right above the threshold current lasing in an excited mode is detected,
with a distinctive splitting in the momentum space; see Fig. 2(d). This
crucial difference between the BEC and the lasing devices is confirmed in the
spatially resolved spectra (near field), as for the BEC behavior, we deal with
a spatially homogeneous gas of photons, presented in Figs. 2(e),(f), where the
condensation occurs in the fundamental mode of the system. In contrast, the
lasing device emits stronger near the aperture edges, where current density
and optical gain are higher, see Fig. 2(g) [34, 35, 36].

We further explore the thermodynamic properties of the photon gas in the BEC
device by extracting the occupancies in energy states. We integrate the
momentum space data considering the density of states, estimated photon
lifetimes, and efficiency of the optical setup (see Methods for details). The
experimental energy distributions at different driving currents are presented
in Fig. 3(a). All data have been successfully fitted with the BE distributions
(eq. (1)), assuming negligible $\Gamma\approx 0$. Additional verification of
the BEC distribution is also carried out, representing the data in logarithmic
form, by transforming the eq. (1) as
$\ln(1+1/N(\varepsilon))=\varepsilon/(k_{B}T)-\mu/(k_{B}T)$, which results as
a linear function of energy; see Fig. 3(b). Furthermore, the data resemble the
textbook behavior of a Bose condensed gas, such as the massive occupation and
threshold-like dependency of the ground state occupancy as a function of the
total number of particles along with the saturation of the excited states seen
in the distributions in Fig. 3(a) and $N_{0}$, $N_{T}$ values summarized in
Fig. 3(c). However, the thermal tails do not have the same slopes, which is
more evident in Fig. 3(b), which implies that, although the photons seem
equilibrated, the temperature of the gas is not equal at different driving
currents. Therefore, we denote the fitting parameters of the BE distribution
as an effective chemical potential $\mu_{\text{eff}}$ and temperature
$T_{\text{eff}}$, because these may not be equal to those set by the device
conditions. Importantly, the geometry of the device imposes the inhomogeneous
current density across the aperture, thus the chemical potential set by the
quasi-Fermi levels and the temperature slightly vary spatially (see
Supplementary Information). Thus, the thermodynamic properties of the photon
gas are a result of the spatially averaged overlap of the optical modes with
the inhomogeneous QWs active medium [37]. The results of the fitting to the
experimental data are presented in Fig. 3(d). The chemical potential is always
negative with respect to the ground mode and approaches zero when condensation
to the ground state occurs, supporting the BEC behavior for an ideal gas. On
the other hand, the temperature is a monotonic function of the driving current
and saturates above the condensation transition to $T_{\text{eff}}\approx
250~{}\text{K}$, which is approximately $T_{\text{eff}}/T\approx 0.85$
compared to the temperature of the heat sink $T=293~{}\text{K}$ (the actual
temperature of the active region is expected to be slightly higher, see
Supplementary Information). From the data in Fig. 3(c), we experimentally
extracted the critical condensation value of the particles
$N_{C}^{\text{exp}}=2604\pm 91$, which is the value of the condensation
temperature $T\approx 220~{}\text{K}$, in line with the experimental value
extracted at the condensation threshold $T_{\rm eff}\approx 223~{}\text{K}$,
Fig. 3(d) (see Supplementary Information).

Figure 4: Measurement of the EOS Points are extracted from the experimental
data. The dashed line is the theoretical EOS for a 2D Bose gas in the
thermodynamic limit. The solid line is calculated by taking into account the
finite collection angle of the optical setup. Error bars represent the 95%
confidence intervals.

These results suggest that we are dealing with a photonic gas that is not in
full thermal and chemical equilibrium with the reservoir, which is the active
region of the device. Equilibration to temperatures lower than the reservoir
by stimulated cooling has recently been predicted for driven-dissipative
bosonic condensates in the fast thermalization limit considering correlations
between all quantum states [38]. An experimental indication for the stimulated
cooling effect can be seen in our data, as there is a saturation of
$T_{\text{eff}}$ above the condensation threshold in Fig. 3(d). Therefore, it
is interesting to examine what is the EOS of the probed photon condensate and
if it follows the EOS for a 2D Bose gas, which is written in the thermodynamic
limit as:

$\mathcal{D}=-\ln\left(1+\exp\left(\frac{\mu}{k_{B}T}\right)\right),$ (2)

where $\mathcal{D}=n\lambda_{T}^{2}$ is the dimensionless phase space density
defined by the photon density $n=N/(\pi R^{2})$ ($\pi R^{2}$ is the surface
area of the aperture and $R$ its radius) and the thermal de Broglie wavelength
of photons $\lambda_{T}=\sqrt{(2\pi\hbar^{2})/(m_{ph}k_{B}T)}$. The EOS is
expressed in normalized quantities, like $\mathcal{D}$ and
$\tilde{\mu}=\mu/(k_{B}T)$, hence the properties of the 2D bosonic gas are
expected to be universal [39, 40]. The measured EOS, expressed by the
experimental values of $\mu_{\text{eff}}$, $T_{\text{eff}}$, is presented in
Fig. 4. The data follows the equilibrium EOS, but with a larger slope in
comparison to the thermodynamic limit. This can be explained by the finite
collection angle (NA of the microscope objective) of the collection optics in
our setup. We cannot detect energies emitted beyond the maximal angle.
Nonetheless, numerical calculations taking this fact into account confirm the
observations (see Supplementary Information), as shown in Fig. 4 as a solid
line, which perfectly matches the experimental data. The observed discrepancy
is analogous to previous reports [41], where the finite trap depth was the
reason for the lower than expected phase space density. In our case, the mode
energies in the mesa are expected to go beyond the value dictated by the
objective NA.

We have demonstrated that emission from a positively detuned VCSEL has the
properties of a homogeneous 2D Bose-Einstein condensed gas of photons in a
finite system. The measured nonequilibrium nature of the gas can be a
signature of reaching the fast thermalization limit, because the cavity is
characterized by a relatively short photon lifetime. Photon condensation in
semiconductor resonators offers the possibility of observing the superfluidity
of a weakly interacting Bose gas. Photon interactions are expected to be
mediated by semiconductor non-linearity, which is significantly enhanced by
the cavity and has a subpicosecond-order response time [42, 43]. In our data
there are no clear indications of such interactions because the cavity energy
shifts are dominated by the current- and temperature-induced changes of the
refractive index, so further studies must be performed, like probing directly
the hydrodynamics of the condensed photons by perturbing them from the steady
state [44, 45]. Nevertheless, in addition to material non-linearities, the
dissipative nature of the photon gas encourages further studies of phase
ordering [46] and universal scaling in a 2D geometry [47, 48] and signs of
non-Hermitian effects [49]. Another direction is to test the fluctuations of
the non-equilibrium BEC and compare it to the BEC in thermal equilibrium [50,
51] and to standard VCSEL operation [52, 53]. Lastly, the mature technology of
semiconductor VCSELs offers the possibility of utilizing the BEC regime to
achieve single-mode emission from large aperture devices characterized by
excellent beam quality, without the need for sophisticated additional
fabrication and processing of the laser mesa. Furthermore, BEC VCSELs can also
be applied in more complex lattice geometries to study topological effects in
well-controlled current-operated devices at room temperature [54].

## Methods

### Thermalization of photons in a semiconductor laser

The principles of light absorption and recombination in an excited
semiconductor QW, depicted in Fig. 1(b), can be described by the following
transition rates [37, 27] for emission:

$R_{\rm
em}(\varepsilon)=R(\varepsilon)f_{c}(\varepsilon,T,\mu_{c})\big{(}1-f_{v}(\varepsilon,T,\mu_{v})\big{)}$
(3)

and absoprtion:

$R_{\rm
abs}(\varepsilon)=R(\varepsilon)f_{v}(\varepsilon,T,\mu_{v})\big{(}1-f_{c}(\varepsilon,T,\mu_{c})\big{)},$
(4)

where
$f_{c,v}=\left[\exp\left((\varepsilon-\mu_{c,v})/(k_{B}T)\right)+1\right]^{-1}$
are the thermalized Fermi-Dirac distributions of electrons in the conduction
and holes in the valence bands. $R(\varepsilon)$ is the transition rate at
energy $\varepsilon$, taking into account the photonic and electronic density
of states, the overlap of the optical modes with the active medium, and the
intrinsic properties of the active medium itsel [37]. The natural consequence
in semiconductors is the van Roosbroeck-Schockley relation, which appears,
after some algebra, from the relation:

$R_{\rm abs}(\varepsilon)/R_{\rm
em}(\varepsilon)=\exp\left(\frac{\varepsilon-\mu}{k_{B}T}\right)$ (5)

with $\mu=\mu_{c}-\mu_{v}$ [28, 37, 27].

Now, the rate equations for the occupation of a mode at $\varepsilon$ is
expressed as:

$\frac{d}{dt}N(\varepsilon)=R_{\text{em}}(\varepsilon)(N(\varepsilon)+1)-R_{\text{abs}}(\varepsilon)N(\varepsilon)-\gamma(\varepsilon)N(\varepsilon)$
(6)

and the steady-state solution gives

$N(\varepsilon)=\frac{R_{\text{em}}(\varepsilon)}{\gamma(\varepsilon)-(R_{\text{em}}(\varepsilon)-R_{\text{abs}}(\varepsilon))},$
(7)

where $\gamma(\varepsilon)=1/\tau(\varepsilon)$ is the decay rate of a photon
from an empty cavity at $\varepsilon$. After dividing the nominator and
denominator by $R_{\text{em}}(\varepsilon)$ and using the van Roosbroeck-
Schockley relation we have:

$N(\varepsilon)=\frac{1}{\exp(\frac{\varepsilon-\mu}{k_{B}T})-1+\Gamma(\varepsilon)}$

as in the main text.

We measure $\Gamma(\varepsilon_{0})$ for the ground mode of the device using
the decay time calculated from the realistic model of the empty cavity
$\gamma(\varepsilon_{0})=1/\tau_{0}=1/3.04~{}\text{ps}^{-1}\approx
0.33~{}\text{ps}^{-1}$. We were able to determine the value of
$R_{\text{em}}(\varepsilon_{0})$ close to the threshold, by measuring the
linewidth dependence of the ground mode occupation below the condensation
threshold $R_{\text{em}}(\varepsilon_{0})=42\pm 3~{}\text{ps}^{-1}$ [55].

### Sample

The VCSEL epitaxial structure is designed for high-speed data communication at
980 nm. The epitaxial structure is monolithically grown on an n-doped GaAs
substrate. The multi-quantum well (MQW) active region is composed of 5 $\rm
In_{0.23}Ga_{0.77}As$ QWs and 6 $\rm GaAs_{0.86}P_{0.14}$ barriers centered in
$\mathrm{Al}_{x}\mathrm{Ga}_{1-x}\mathrm{As}$ cavity graded from $x=0.38$ to
$0.80$ with an optical cavity thickness of $\lambda/2$. The cavity is
sandwiched by 15.5-pair GaAs/$\rm Al_{0.9}Ga_{0.1}As$ top and 37.5-pair bottom
distributed Bragg reflector (DBR) mirrors. The top and bottom DBRs are C-doped
for the p-type and Si-doped for the n-type, respectively. In both mirrors,
graded interfaces are incorporated for lower electrical resistance of the
structure. Importantly, two 20 nm thick $\rm Al_{0.98}Ga_{0.02}As$ layers are
placed to form oxide apertures in the first nodes of the standing wave at the
top and bottom of the cavity. These oxide layers are halfway in the optical
cavity and halfway in the first pair of layers in the DBRs.

The VCSELs are processed using standard top-down photolithography. In the
first step, the Ti/Pt/Au p-type contact rings are deposited with the use of
electron beam deposition (E-beam). The mesa structures are then patterned and
etched using inductively coupled plasma reactive-ion etching (ICP-RIE) in a
$\rm Cl_{2}/BCl_{3}$-based plasma. After etching, current confinement
apertures are formed by selective wet thermal oxidation of $\rm
Al_{0.98}Ga_{0.02}As$ layers in an oxidation oven in a nitrogen atmosphere
with overpressure of water vapor and at high temperature (420°C). In the
following step, horseshoe Ni/AuGe/Au n-type contact pads are deposited and
annealed in a rapid thermal processing (RTP) furnace. The structures were then
planarized with the use of a spin-on dielectric polymer of benzocyclobutene
(BCB). The BCB layer is patterned with the use of photolithography and RIE
etching in a $\rm CF_{4}$-based plasma to selectively open surface areas to
the bias to the p- and n-type contacts. In the final step, the ground-signal-
ground (GSG) Cr/Pt/Au contact pads are deposited.

### Experimental setup

The sample used in this study was the fully processed quarter of the whole
epitaxial wafer. The sample was placed on a thermo-electrically cooled plate
(Thorlabs PTC1) with an additional temperature sensor buried inside a custom
heatsink plate mounted on the top to ensure control of the temperature in the
close vicinity of the sample. The temperature of the heatsink was set to 20 °C
throughout all experiments. The temperature-controlled plate was placed on a
manual translation stage. The sample was contacted by a microwave probe (GGB
Picoprobe 40A) located in an additional manual translation stage. The devices
were driven with a direct current by a stabilized precise source/measure unit
(Keysight B2901B).

The device emission was collected using an infinity-corrected objective of NA
= 0.65 (Mitutoyo M Plan Apo NIR HR 50x). As written in the main text, to
measure the momentum spectra (far field), we imaged the back focal plane of
the objective with a set of achromatic lenses onto the 0.3 m-focal length
monochromator entrance slit (Princeton Instruments Acton SP-2300i) and
dispersed by a grating (1200 grooves/mm) onto the electron-multiplied charge-
coupled device (EMCCD – Teledyne Princeton Instruments ProEM-HS:1024BX3). To
record the spatially resolved spectra (near field) one of the lenses was
removed from the optical path, which enabled one to project the real-space
image onto the monochromator slit. This lens was mounted on a flip mount,
allowing quick and convenient switching between the two measurement modes of
the setup.

### Analysis of the momentum space

Taking advantage of homogeneous emission from the BEC device, we determined
the thermodynamic properties of the photon gas from the momentum space. We
extracted the mean photon occupation distribution by integrating the momentum
space emission, using the standard procedure used in cavity-polariton physics
[56].

The mean number of photons collected at a CCD pixel row representing a chosen
$k$-state is represented as follows:

$N_{\rm ph}(k)=\eta\frac{dN_{\rm CCD}(k)}{dt}\tau(k),$ (8)

where $\eta$ is the calibrated collection efficiency of our setup, $dN_{\rm
CCD}(k)/dt$ is the count rate per second on the CCD camera and $\tau(k)$ is
the photon lifetime at $k$. The photon lifetime was estimated from the
experiment by extracting the emission linewidth
$\Delta\varepsilon_{k}=\hbar/\tau(k)$ [57] by fitting a Lorentzian function to
the data from a $k$-state pixel row.

Subsequently, the occupation number at the $k$-state is calculated taking into
account the number of states subtended by a pixel at $k$-position in
cylindrical coordinates $N_{\rm st}(k)=k\Delta k\Delta\phi(4\pi/S)^{-1}$,
where $S$ is the surface area of the device aperture. The number of states in
momentum space was confirmed by numerical simulations of the optical modes
confined in the device (see Supplementary Information). The final expression
is as follows:

$N(\varepsilon(k))=\frac{N_{\rm ph}(k)}{N_{\rm
st}(k)}=\frac{4\pi^{2}\eta}{2k\Delta k\Delta\phi S}\frac{dN_{\rm
CCD}(k)}{dt}\tau(k),$ (9)

which considers also the spin degeneracy 2 of all states, as our experiment
was not polarization resolved. The energy $\varepsilon(k)$ of the measured
$k$-state is extracted from the fitted Lorentzian peak.

### Acknowledgments

We gratefully thank Maciej Dems for his support in improving the numerical
simulation codes used in this work. MP and ANP acknowledge support from the
Polish National Science Center, grant SONATA no. 2020/39/D/ST3/03546. AP
acknowledges financial support by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) via the Collaborative Research Center SFB/TR185
(Project No. 277625399).

## References

  * Einstein [1925] A. Einstein, Quantentheorie des einatomigen idealen Gases, Sitzungsberichte der Preussischen Akademie der Wissenschaften 1, 3 (1925).
  * Pitaevskii and Stringari [2016] L. P. Pitaevskii and S. Stringari,  _Bose-Einstein Condensation and Superfluidity_ , 1st ed. (Oxford University Press, 2016) p. 567.
  * Anderson  _et al._ [1995] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Observation of Bose-Einstein condensation in a dilute atomic vapor, Science 269, 198 (1995).
  * Davis  _et al._ [1995] K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Dufree, D. M. Kurn, and W. Ketterle, Bose-Einstein condensation in a gas of sodium atoms, Physical Review Letters 75, 3969 (1995).
  * Demokritov  _et al._ [2006] S. O. Demokritov, V. E. Demidov, O. Dzyapko, G. A. Melkov, A. A. Serga, B. Hillebrands, and A. N. Slavin, Bose–Einstein condensation of quasi-equilibrium magnons at room temperature under pumping, Nature 443, 430 (2006).
  * Butov  _et al._ [2002] L. V. Butov, A. C. Gossard, and D. S. Chemla, Macroscopically ordered state in an exciton system, Nature 418, 751 (2002).
  * Wang  _et al._ [2019] Z. Wang, D. A. Rhodes, K. Watanabe, T. Taniguchi, J. C. Hone, J. Shan, and K. F. Mak, Evidence of high-temperature exciton condensation in two-dimensional atomic double layers, Nature 574, 76 (2019).
  * Alloing  _et al._ [2014] M. Alloing, M. Beian, M. Lewenstein, D. Fuster, Y. González, L. González, R. Combescot, M. Combescot, and F. Dubin, Evidence for a Bose-Einstein condensate of excitons, EPL (Europhysics Letters) 107, 10012 (2014).
  * Hakala  _et al._ [2018] T. K. Hakala, A. J. Moilanen, A. I. Väkeväinen, R. Guo, J. P. Martikainen, K. S. Daskalakis, H. T. Rekola, A. Julku, and P. Törmä, Bose-Einstein condensation in a plasmonic lattice, Nature Physics 14, 739 (2018).
  * Kasprzak  _et al._ [2006] J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. Keeling, F. M. Marchetti, M. H. Szymánska, R. André, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, L. S. Dang, K. J., _et al._ , Bose-Einstein condensation of exciton polaritons, Nature 443, 409 (2006).
  * Sun  _et al._ [2017] Y. Sun, P. Wen, Y. Yoon, G. Liu, M. Steger, L. N. Pfeiffer, K. West, D. W. Snoke, and K. A. Nelson, Bose-Einstein condensation of long-lifetime polaritons in thermal equilibrium, Physical Review Letters 118, 016602 (2017).
  * Klaers  _et al._ [2010] J. Klaers, J. Schmitt, F. Vewinger, and M. Weitz, Bose-Einstein condensation of photons in an optical microcavity., Nature 468, 545 (2010).
  * Damm  _et al._ [2016] T. Damm, J. Schmitt, Q. Liang, D. Dung, F. Vewinger, M. Weitz, and J. Klaers, Calorimetry of a Bose-Einstein-condensed photon gas, Nature Communications 7, 11340 (2016).
  * Damm  _et al._ [2017] T. Damm, D. Dung, F. Vewinger, M. Weitz, and J. Schmitt, First-order spatial coherence measurements in a thermalized two-dimensional photonic quantum gas, Nature Communications 8, 158 (2017).
  * Schmitt  _et al._ [2014] J. Schmitt, T. Damm, D. Dung, F. Vewinger, J. Klaers, and M. Weitz, Observation of grand-canonical number statistics in a photon Bose-Einstein condensate, Phys. Rev. Lett. 112, 030401 (2014).
  * Walker  _et al._ [2020] B. T. Walker, J. D. Rodrigues, H. S. Dhar, R. F. Oulton, F. Mintert, and R. A. Nyman, Non-stationary statistics and formation jitter in transient photon condensation, Nature Communications 11, 1390 (2020).
  * Walker  _et al._ [2018] B. T. Walker, L. C. Flatten, H. J. Hesten, F. Mintert, D. Hunger, A. A. Trichet, J. M. Smith, and R. A. Nyman, Driven-dissipative non-equilibrium Bose–Einstein condensation of less than ten photons, Nature Physics 14, 1173 (2018).
  * Kirton and Keeling [2013] P. Kirton and J. Keeling, Nonequilibrium model of photon condensation, Phys. Rev. Lett. 111, 100404 (2013).
  * Radonjić  _et al._ [2018] M. Radonjić, W. Kopylov, A. Balaž, and A. Pelster, Interplay of coherent and dissipative dynamics in condensates of light, New Journal of Physics 20, 055014 (2018).
  * Weill  _et al._ [2019] R. Weill, A. Bekker, B. Levit, and B. Fischer, Bose–Einstein condensation of photons in an Erbium–Ytterbium co-doped fiber cavity, Nature Communications 10, 747 (2019).
  * Barland  _et al._ [2021] S. Barland, P. Azam, G. L. Lippi, R. A. Nyman, and R. Kaiser, Photon thermalization and a condensation phase transition in an electrically pumped semiconductor microresonator, Optics Express 29, 8368 (2021).
  * Kammann  _et al._ [2012] E. Kammann, H. Ohadi, M. Maragkou, A. V. Kavokin, and P. G. Lagoudakis, Crossover from photon to exciton-polariton lasing, New Journal of Physics 14, 105003 (2012).
  * Bajoni  _et al._ [2007] D. Bajoni, P. Senellart, A. Lemaître, and J. Bloch, Photon lasing in GaAs microcavity: Similarities with a polariton condensate, Physical Review B 76, 201305 (2007).
  * Haghighi and Lott [2019] N. Haghighi and J. A. Lott, State-of-the-art VCSEL epitaxial designs (Invited), paper C-01, 10th International Conference on Materials for Advanced Technologies (ICMAT), Singapore (23-28 June 2019) (2019).
  * Haghighi  _et al._ [2019] N. Haghighi, P. Moser, and J. A. Lott, Power, bandwidth, and efficiency of single VCSELs and small VCSEL arrays, IEEE Journal of Selected Topics in Quantum Electronics 25, 1 (2019).
  * Knox  _et al._ [1986] W. H. Knox, C. Hirlimann, D. A. B. Miller, J. Shah, D. S. Chemla, and C. V. Shank, Femtosecond excitation of nonthermal carrier populations in GaAs quantum wells, Phys. Rev. Lett. 56, 1191 (1986).
  * Coldren  _et al._ [2012] L. A. Coldren, S. W. Corzine, and M. L. Masanovic,  _Diode Lasers and Photonic Integrated Circuits_ , 2nd ed. (John Wiley & Sons, 2012).
  * Wurfel [1982] P. Wurfel, The chemical potential of radiation, J. Phys. C: Solid State Phys 15, 3967 (1982).
  * Bhattacharya  _et al._ [2012] R. Bhattacharya, B. Pal, and B. Bansal, On conversion of luminescence into absorption and the van Roosbroeck-Shockley relation, Appl. Phys. Lett 100, 222103 (2012).
  * Bernard and Duraffourg [1961] M. G. A. Bernard and G. Duraffourg, Laser conditions in semiconductors, Physica Status Solidi B 1, 669 (1961).
  * Mogg  _et al._ [2004] S. Mogg, N. Chitica, U. Christiansson, R. Schatz, P. Sundgren, C. Asplund, and M. Hammar, Temperature sensitivity of the threshold current of long-wavelength InGaAs-GaAs VCSELs with large gain-cavity detuning, IEEE Journal of Quantum Electronics 40, 453 (2004).
  * Balili  _et al._ [2006] R. B. Balili, D. W. Snoke, L. Pfeiffer, and K. West, Actively tuned and spatially trapped polaritons, Appl. Phys. Lett 88, 31110 (2006).
  * Hegarty  _et al._ [1999] S. P. Hegarty, G. Huyet, J. G. Mcinerney, and K. D. Choquette, Pattern formation in the transverse section of a laser with a large fresnel number, Physical Review Letters 82, 1434 (1999).
  * Michalzik  _et al._ [1998] R. Michalzik, M. Grabherr, and K. J. Ebeling, High-power VCSELs: modeling and experimental characterization, Proc. SPIE 3286, 206 (1998).
  * Degen  _et al._ [1999] C. Degen, I. Fischer, W. Elsaesser, J. l Vey, and P. Schnitzer, Transverse modes in oxide confined VCSELs: Influence of pump profile, spatial hole burning, and thermal effects, Optics Express 5, 38 (1999).
  * Ackemann  _et al._ [2000] T. Ackemann, S. Barland, M. Cara, S. Balle, J. R. Tredicce, R. Jäger, M. Grabherr, M. Miller, and K. J. Ebeling, Spatial mode structure of bottom-emitting broad-area vertical-cavity surface-emitting lasers, Journal of Optics B: Quantum and Semiclassical Optics 2, 406 (2000).
  * Loirette-Pelous and Greffet [2023] A. Loirette-Pelous and J.-J. Greffet, Photon Bose-Einstein condensation and lasing in semiconductor cavities (2023), arXiv:2301.13073 [cond-mat.quant-gas] .
  * Shishkov  _et al._ [2022] V. Y. Shishkov, E. S. Andrianov, A. V. Zasedatelev, P. G. Lagoudakis, and Y. E. Lozovik, Exact analytical solution for the density matrix of a nonequilibrium polariton Bose-Einstein condensate, Physical Review Letters 128, 065301 (2022).
  * Hung  _et al._ [2011] C. L. Hung, X. Zhang, N. Gemelke, and C. Chin, Observation of scale invariance and universality in two-dimensional Bose gases, Nature 470, 236 (2011).
  * Yefsah  _et al._ [2011] T. Yefsah, R. Desbuquois, L. Chomaz, K. J. Günter, and J. Dalibard, Exploring the thermodynamics of a two-dimensional Bose gas, Physical Review Letters 107, 130401 (2011).
  * Busley  _et al._ [2022] E. Busley, L. E. Miranda, A. Redmann, C. Kurtscheid, K. K. Umesh, F. Vewinger, M. Weitz, and J. Schmitt, Compressibility and the equation of state of an optical quantum gas in a box, Science 375, 1403 (2022).
  * Kriso  _et al._ [2021] C. Kriso, T. Bergmeier, N. Giannini, A. R. Albrecht, M. Sheik-Bahae, S. Benis, S. Faryadras, E. W. V. Stryland, D. J. Hagan, M. Koch, G. Mette, and A. Rahimi-Iman, Probing the ultrafast gain and refractive index dynamics of a VECSEL, Applied Physics Letters 119, 191105 (2021).
  * Yüce  _et al._ [2016] E. Yüce, G. Ctistis, J. Claudon, J.-M. Gérard, and W. L. Vos, Optimal all-optical switching of a microcavity resonance in the telecom range using the electronic Kerr effect, Optics Express 24, 239 (2016).
  * Estrecho  _et al._ [2021] E. Estrecho, M. Pieczarka, M. Wurdack, M. Steger, K. West, L. N. Pfeiffer, D. W. Snoke, A. G. Truscott, and E. A. Ostrovskaya, Low-energy collective oscillations and Bogoliubov sound in an exciton-polariton condensate, Physical Review Letters 126, 75301 (2021).
  * Peinke  _et al._ [2021] E. Peinke, T. Sattler, G. M. Torelly, P. L. Souza, S. Perret, J. Bleuse, J. Claudon, W. L. Vos, and J.-M. Gérard, Tailoring the properties of quantum dot-micropillars by ultrafast optical injection of free charge carriers, Light: Science & Applications 10, 215 (2021).
  * Gladilin and Wouters [2020] V. N. Gladilin and M. Wouters, Vortices in nonequilibrium photon condensates, Phys. Rev. Lett. 125, 215301 (2020).
  * Fontaine  _et al._ [2022] Q. Fontaine, D. Squizzato, F. Baboux, I. Amelio, A. Lemaître, M. Morassi, I. Sagnes, L. L. Gratiet, A. Harouri, M. Wouters, I. Carusotto, A. Amo, M. Richard, A. Minguzzi, L. Canet, S. Ravets, and J. Bloch, Kardar–Parisi–Zhang universality in a one-dimensional polariton condensate, Nature 608, 687 (2022).
  * Bloch  _et al._ [2022] J. Bloch, I. Carusotto, and M. Wouters, Non-equilibrium Bose–Einstein condensation in photonic systems, Nature Reviews Physics 4, 470 (2022).
  * Khurgin [2020] J. B. Khurgin, Exceptional points in polaritonic cavities and subthreshold Fabry–Perot lasers, Optica 7, 1015 (2020).
  * Christensen  _et al._ [2021] M. B. Christensen, T. Vibel, A. J. Hilliard, M. B. Kruk, K. Pawłowski, D. Hryniuk, K. Rzążewski, M. A. Kristensen, and J. J. Arlt, Observation of microcanonical atom number fluctuations in a Bose-Einstein condensate, Phys. Rev. Lett. 126, 153601 (2021).
  * Öztürk  _et al._ [2023] F. E. Öztürk, F. Vewinger, M. Weitz, and J. Schmitt, Fluctuation-dissipation relation for a Bose-Einstein condensate of photons, Phys. Rev. Lett. 130, 033602 (2023).
  * Agarwal [1972] G. S. Agarwal, Fluctuation-dissipation theorems for systems in non-thermal equilibrium and application to laser light, Physics Letters A 38, 93 (1972).
  * Aßmann  _et al._ [2009] M. Aßmann, F. Veit, M. Bayer, M. V. D. Poel, and J. M. Hvam, Higher-order photon bunching in a semiconductor microcavity, Science 325, 297 (2009).
  * Dikopoltsev  _et al._ [2021] A. Dikopoltsev, T. H. Harder, E. Lustig, O. A. Egorov, J. Beierlein, A. Wolf, M. Emmerling, C. Schneider, S. Höfling, M. Segev, and S. Klembt, Topological insulator vertically-emitting laser array, Science 373, 1514 (2021).
  * Pérez  _et al._ [2014] P. Pérez, A. Valle, I. Noriega, and L. Pesquera, Measurement of the intrinsic parameters of single-mode VCSELs, Journal of Lightwave Technology 32, 1601 (2014).
  * Pieczarka  _et al._ [2019] M. Pieczarka, M. Boozarjmehr, E. Estrecho, Y. Yoon, M. Steger, K. West, L. N. Pfeiffer, K. A. Nelson, D. W. Snoke, A. G. Truscott, and E. A. Ostrovskaya, Effect of optically induced potential on the energy of trapped exciton-polaritons below the condensation threshold, Physical Review B 100, 85301 (2019).
  * Pollnau and Eichhorn [2020] M. Pollnau and M. Eichhorn, Spectral coherence, part i: Passive-resonator linewidth, fundamental laser linewidth, and Schawlow-Townes approximation, Progress in Quantum Electronics 72, 100255 (2020).


#  Analytic and numerical approaches for depictive 3-loop integrals using
sector decomposition

Elise de Donckera, Tadashi Ishikawab, Kiyoshi Katoc,111e-mail :
<EMAIL_ADDRESS>, Fukuko Yuasab
aWestern Michigan University, Kalamazoo MI 49008, USA,
bHigh Energy Accelerator Research Organization (KEK), 1-1 Oho Tsukuba,
Ibaraki 305-0801, Japan,
cKogakuin University, Nishi-Shinjuku 1-24, Shinjuku, Tokyo 163-8677, Japan

###### Abstract

Four 3-loop two-point functions are studied analytically and numerically using
a simplified sector decomposition method. The coefficients of the ultraviolet
divergent part are determined analytically, and those of the finite part are
computed numerically. The energy dependence of the integrals is shown
explicitly, and a discussion of its behavior is presented.

##  1 Introduction

As high-precision experimental data accumulates, it becomes increasingly
important to estimate the theoretical radiative correction that matches its
accuracy. Additionally, by examining the contribution of virtual particles
included in the radiative correction, it would be possible to predict unknown
particles. The evaluation of such radiative corrections requires the
calculation of loop diagrams with masses as an essential element.

The sector decomposition (SD) method [1, 2, 3, 4, 5, 6] is well-established
for calculating loop integrals. In this paper, we use a simple form of the SD
method. Conventionally, the number of regions increases depending on the
number of particles and the degree of the loops, which has been a problem.
Therefore, various techniques [7, 8, 9, 10, 11] have been used to suppress
this increase. However, with the advances in computer capabilities, it has
been found that the method can achieve precise calculations up to a certain
level by combining it with automatic processing of the integration. From the
basic formula derived by this method, it is possible to obtain a specific
expression of the coefficients that represent the ultraviolet (UV) divergence
analytically, and numerical calculations with sufficient accuracy can be
performed even for the non-divergent terms.

We have developed fully numerical approaches [12, 13, 14] using non-adaptive
or adaptive integration with transformations and/or extrapolation to manage
the singular behavior of loop integrals. The method consists of an accurate
integration and estimation of the coefficients. The former includes the
numerical integration techniques, such as adaptive integration [15, 16, 17],
the double-exponential formula [18], and Quasi-Monte Carlo integration based
on lattice rules [19], and the latter uses nonlinear extrapolation [20, 21] or
a linear solver. We have succeeded to manage the singularity from both the UV-
divergence and in the physical kinematics below and beyond the threshold in
the squared momentum. To achieve this, we use a double extrapolation technique
and have already shown some results for the 2-loop [22] and for the 3-loop
order [23].

In the following sections, as a concrete example of this method, we provide
results for several 3-loop diagrams in Fig. 1 that are of a form we call
“complete”. We deal with the UV singularity by dimensional regularization,
where the space-time dimension is set to $n=4-2\varepsilon$. A 3-loop integral
${\bf I}$ is given by its Laurent series in $\varepsilon$:

${\bf I}=\frac{{\bf C}_{-3}}{\varepsilon^{3}}+\frac{{\bf
C}_{-2}}{\varepsilon^{2}}+\frac{{\bf C}_{-1}}{\varepsilon}+{\bf
C}_{0}+\cdots\,.$ (1)

As referenced in Section 2, the Loop (I), sometimes called banana diagram, is
studied in the literature [24, 25, 26, 27] based on a general framework. The
Loop (I), (II) and (IV) diagrams are studied in [28, 29] using integration-by-
parts (IBP) relations for the generic mass case.

Further in Section 2, we give the formulation of our calculation approach and
analytical results for the divergent part, i.e., for ${\bf C}_{-3},{\bf
C}_{-2}$ and ${\bf C}_{-1}$. In Section 3, we describe the numerical methods
and the results by numerical computation for the coefficients up to and
including ${\bf C}_{0}$. Finally, we summarize and provide future prospects in
Section 4.

##  2 Formalism and analytic results

We define the scalar integral with numerator $1$ as

${\bf
I}=\prod_{k=1}^{L}\int\frac{d^{n}\ell_{k}}{(2\pi)^{n}}\frac{1}{D_{1}\ldots
D_{N}}$ (2)

where $D_{j}=p_{j}^{2}-m_{j}^{2}$ is the inverse propagator of the $j$-th
line, $N$ is the number of internal lines and $L$ is the number of loops.

Subsequently we consider 3-loop two-point functions, so that $L=3$. The
external momentum that flows into the diagram is $P$ and we denote $s=P^{2}$.
Using the Nakanishi-Cvitanović-Kinoshita formalism [30, 31, 32, 33], the
integral is given by

${\bf I}=(-1)^{N}\frac{\Gamma(N-3n/2)}{(4\pi)^{3n/2}}\int\prod
dx_{j}\frac{\delta(1-\sum x_{j})}{U^{n/2}V^{N-3n/2}}\ .$ (3)
$V=M^{2}-s\frac{W}{U},\qquad M^{2}=\sum m_{j}^{2}x_{j}-i\rho$ (4)

where $U$ ($W$) is a sum of monomials and each monomial is a product of three
(four) different $x$’s with coefficient $1$. Here, $i\rho$, the infinitesimal
imaginary part of mass, is clearly expressed to be used in the numerical study
of the next section. Factors are separated from the integral as follows:

${\bf
I}=(-1)^{N}\left(\frac{1}{(4\pi)^{2-\varepsilon}}\right)^{3}\Gamma(1+3\varepsilon)\times\tilde{I},\qquad\tilde{I}=\frac{\Gamma(N-6+3\varepsilon)}{\Gamma(1+3\varepsilon)}\times
I$ (5)

$I=\int d{\bf\Gamma}\frac{V^{6-N-3\varepsilon}}{U^{2-\varepsilon}},\qquad
d{\bf\Gamma}=\prod_{j=1}^{N}dx_{j}~{}\delta(1-\sum x_{j})$ (6)

The UV singularity appears in the Gamma function if $N\leq 6$, and in the
integral when $U=0$, so that the UV singularity of $I$ is
$O(\varepsilon^{-2})$ or less:

$I=\frac{C_{-2}}{\varepsilon^{2}}+\frac{C_{-1}}{\varepsilon^{1}}+C_{0}+C_{1}\varepsilon+\cdots$
(7)

and the ${\bf C}$’s in Eq. (1) can be calculated from the $C_{j}$.

(I)$P$(II)$P$(III)$P$(IV)$P$123412345123456123456 Figure 1: 3-loop diagrams

We consider the (I), (II), (III), (IV) 3-loop diagrams in Fig. 1 where
$N=4,5,6,6.$ The Feynman parameter $x_{j}$ is assigned to the line $j$ in the
figure. The $U,W$ functions are as given below. For the Feynman parameters, we
use the notation $x_{ab\cdots}=x_{a}+x_{b}+\cdots$.

$\begin{array}[]{ll}\mathrm{Loop~{}(I)}&U=x_{1}x_{2}x_{3}+x_{1}x_{2}x_{4}+x_{1}x_{3}x_{4}+x_{2}x_{3}x_{4},\\\
&W=x_{1}x_{2}x_{3}x_{4}.\\\
\mathrm{Loop~{}(II)}&U=x_{1}x_{2}x_{3}+(x_{1}x_{2}+x_{1}x_{3}+x_{2}x_{3})x_{45},\\\
&W=(x_{12}x_{3}x_{4}+x_{1}x_{2}x_{34})x_{5}.\\\
\mathrm{Loop~{}(III)}&U=x_{1}x_{2}x_{3}+(x_{1}x_{2}+x_{1}x_{3}+x_{2}x_{3})x_{456},\\\
&W=(x_{1}x_{2}x_{3}+(x_{1}x_{2}+x_{1}x_{3}+x_{2}x_{3})x_{45})x_{6}.\\\
\mathrm{Loop~{}(IV)}&U=x_{12}x_{34}x_{56}+(x_{12}+x_{34})x_{5}x_{6},\\\
&W=(x_{1}x_{2}x_{3}+x_{1}x_{2}x_{4}+x_{1}x_{3}x_{4}+x_{2}x_{3}x_{4})x_{56}+x_{13}x_{24}x_{5}x_{6}.\\\
\end{array}$ (8)

We introduce the following variable transformations. Here $\bar{a}=1-a$.

$\begin{array}[]{llllll}\mathrm{Loop~{}(I)}&x_{1}=y_{1},&x_{2}=y_{2},&x_{3}=y_{3},&x_{4}=y_{4}.\\\
\mathrm{Loop~{}(II)}&x_{1}=y_{1},&x_{2}=y_{2},&x_{3}=y_{3},&x_{4}=\bar{w}y_{4},&x_{5}=wy_{4}.\\\
\mathrm{Loop~{}(III)}&x_{1}=y_{1},&x_{2}=y_{2},&y_{3}=y_{3},\\\
&x_{4}=y_{4}w\bar{z},&x_{5}=y_{4}wz,&x_{6}=y_{4}\bar{w}.\\\
\mathrm{Loop~{}~{}(IV)}&x_{1}=y_{1}\bar{w}_{1},&x_{2}=y_{1}w_{1},&x_{3}=y_{2}\bar{w}_{2},&x_{4}=y_{2}w_{2},\\\
&x_{5}=y_{3},&x_{6}=y_{4}.\\\ \end{array}$ (9)

After the transformation, all $U$ functions are of the form

$U=y_{1}y_{2}y_{3}+y_{1}y_{2}y_{4}+y_{1}y_{3}y_{4}+y_{2}y_{3}y_{4}.$ (10)

We consider the $U$ function being “complete” in the sense that all possible
monomials in $y_{1},y_{2},y_{3}$, and $y_{4}$ are present. A consequence of
this property is described later. The $W$ function is also expressed in terms
of the $y$’s and additional variables ($z,w$).

The phase space can be written as

$\begin{array}[]{ll}\mathrm{Loop~{}(I)}&d{\bf\Gamma}=d\tilde{\Gamma}\\\
\mathrm{Loop~{}(II)}&d{\bf\Gamma}=dw\,d\tilde{\Gamma}\,y_{4}\\\
\mathrm{Loop~{}(III)}&d{\bf\Gamma}=w\,dw\,dz\,d\tilde{\Gamma}\,y_{4}^{2}\\\
\mathrm{Loop~{}(IV)}&d{\bf\Gamma}=dw_{1}\,dw_{2}\,d\tilde{\Gamma}\,y_{1}y_{2}\\\
\end{array}$ (11)

where

$d\tilde{\Gamma}=dy_{1}\,dy_{2}\,dy_{3}\,dy_{4}\,\delta(1-\sum y_{j})$ (12)

and

$R=\\{(y_{1},y_{2},y_{3},y_{4})~{}|~{}\sum_{j}y_{j}=1\\}$ (13)

is the full phase space in the $y$’s.

Hereafter, we assume that all masses are the same,

$m=m_{j}\quad(j=1,\ldots,N).$ (14)

Extension to the general case is possible. In Loop (III), if $m_{4}=m_{5}$,
the $z$-dependence disappears.

The calculation of Eq. (6) is done by the following procedure. The region $R$
is divided into regions $R(klm)$ in the following three steps where $k,l,m,n$
is a permutation of $1,2,3,4$.

First step: The region $R$ is divided into four regions as $R=\cup_{j}R(j)$
where $R(k)=\\{(y_{1},y_{2},y_{3},y_{4})$ $\in
R~{}|~{}y_{k}=\mathrm{max}_{j}\\{y_{j}\\}\\}$. In $R(k)$ we perform the
transformation $y_{j}=s_{k}t_{j},\,t_{k}=1$ and we integrate over $s_{k}$ to
eliminate the delta-function. Then $R(k)$ is the 3-dimensional unit cube of
$(t_{l},t_{m},t_{n})$. The phase space is

$d\tilde{\Gamma}=s_{k}^{4}dt_{l}dt_{m}dt_{n}$ (15)

and $s_{k}^{-1}=\sum_{j=1}^{4}t_{j}$.

Second step: The region $R(k)$ is divided into three regions as
$R(k)=\cup_{j}R(kj)$ where $R(kl)=\\{(t_{l},t_{m},t_{n})\in
R(k)~{}|~{}t_{l}=\mathrm{max}_{j}\\{t_{j}\\}\\}$. In $R(kl)$ we perform the
transformation $t_{l}=t_{l},\,t_{m}=t_{l}u_{m},\,t_{n}=t_{l}u_{n}$. The phase
space is $dt_{l}dt_{m}dt_{n}=t_{l}^{2}dt_{l}du_{m}du_{n}$.

Third step: The region $R(kl)$ is divided into two regions as
$R(kl)=\cup_{j}R(klj)$ where $R(klm)=\\{(t_{l},t_{l}u_{m},t_{l}u_{n})\in
R(kl)~{}|~{}u_{m}=\mathrm{max}_{j}\\{u_{j}\\}\\}$. In $R(klm)$ we perform the
transformation $u_{m}=u_{m},\,u_{n}=u_{m}v_{n}$. The phase space is
$dt_{l}dt_{m}dt_{n}=t_{l}^{2}u_{m}dt_{l}du_{m}dv_{n}$.

In summary, in the region $R(klm)$, the $y$-variables are transformed as

$y_{k}=s_{k},\quad y_{l}=s_{k}t_{l},\quad y_{m}=s_{k}t_{l}u_{m},\quad
y_{n}=s_{k}t_{l}u_{m}v_{n}$ (16)

and the $U,\,V$ functions are

$U=s_{k}^{3}t_{l}^{2}u_{m}f,\qquad W=s_{k}^{4}t_{l}^{2}u_{m}q,\quad
V/s_{k}=G_{k}=m^{2}s_{k}^{-1}-s\frac{q}{f}-i\rho\,,$ (17)

where the $f$ and $q$ functions are introduced.

The integral in the region $R(klm)$ is denoted as $I(klm)$ and $I$ is
expressed as the sum of the integrals $I(klm)$. Hereafter, we omit the
subscript from $t,u,v,G$ when it is redundant.

$\begin{array}[]{lll}\mathrm{Loop~{}(I)}&I(klm)={\displaystyle\int
d\Gamma\frac{h}{t^{2-2\varepsilon}u^{1-\varepsilon}}H},&h=1,\\\ {}&{}&{}\\\
\mathrm{Loop~{}(II)}&I(klm)={\displaystyle\int_{0}^{1}dw\,\int
d\Gamma\frac{h}{t^{2-2\varepsilon}u^{1-\varepsilon}}H},&h=y_{4}/s_{k},\\\
{}&{}&{}\\\ \mathrm{Loop~{}(III)}&I(klm)={\displaystyle\int_{0}^{1}w\,dw\,\int
d\Gamma\frac{h}{t^{2-2\varepsilon}u^{1-\varepsilon}}H},&h=y_{4}^{2}/s_{k}^{2},\\\
{}&{}&{}\\\
\mathrm{Loop~{}(IV)}&I(klm)={\displaystyle\int_{0}^{1}dw_{1}\int_{0}^{1}dw_{2}\,\int
d\Gamma\frac{h}{t^{2-2\varepsilon}u^{1-\varepsilon}}H},&h=y_{1}y_{2}/s_{k}^{2}\,.\\\
\end{array}$ (18)

$\int d\Gamma=\int_{0}^{1}dt\int_{0}^{1}du\int_{0}^{1}dv\,,\qquad
H=\frac{G^{6-N-3\varepsilon}}{f^{2-\varepsilon}}.$ (19) $f=1+v+uv+tuv,\qquad
G=m^{2}(1+t(1+u(1+v)))-s\frac{q}{f}-i\rho.$ (20)

Eq. (18) is our basic formula to calculate the loop integral analytically and
numerically. It should also be noted that $s_{k}$ does not occur in $h$ or
$q$. Since $U$ is “complete”, the form of the $f$ function after the three
transformation steps is shown in Eq. (20), so that it is a regular function,
and the UV singularity from the $t$\- and $u$-integrals is under control as
the explicit factor in Eq. (18). Also, the $f$ function is common in all
regions by the “completeness”. If $U$ is not “complete”, we need another
transformation step in several regions and the $f$ function is determined in
each region.

Although there are 24 regions $R(klm)$, since we consider the case of equal
masses, the integral $I$ is expressed as follows using the symmetry of the
diagram:

$\begin{array}[]{ll}\mathrm{Loop~{}(I)}&I=24I(123)\\\
\mathrm{Loop~{}(II)}&I=6(I(123)+I(124)+I(142)+I(412))\\\
\mathrm{Loop~{}(III)}&I=6(I(123)+I(124)+I(142)+I(412))\\\
\mathrm{Loop~{}(IV)}&I=4(I(123)+I(132)+I(312)+I(341)+I(314)+I(134))\end{array}$
(21)

In each region, $h$ and $q$ can be determined by Eq. (16).

A procedure for the computation of $I(klm)$ is sketched in Appendices A and B.
We present the analytic results for the first three coefficients of $I$ in the
Laurent series in $\varepsilon$.

Loop (I)

$I=\frac{1}{\varepsilon^{2}}12(m^{2})^{2}-\frac{1}{\varepsilon}m^{2}\left(2s+8m^{2}+36m^{2}\log
m^{2}\right)$ $+\frac{1}{6}s^{2}-(6\pi^{2}+48)(m^{2})^{2}+m^{2}\log
m^{2}\left(6s+24m^{2}+54m^{2}\log m^{2}\right)+O(\varepsilon)$ (22)

Loop (II)

$I=\frac{1}{\varepsilon^{2}}3m^{2}+\frac{1}{\varepsilon}6m^{2}\left(-\frac{1}{3}-\frac{3}{4}(j_{s}^{(1)}(s)+\log
m^{2})\right)+\frac{3}{4}m^{2}(1+12\log m^{2})j_{s}^{(1)}(s)$
$+\frac{9}{4}m^{2}j_{s}^{(2)}(s)+\frac{1}{4}s+\frac{3}{4}m^{2}(7\log
m^{2}+3(\log m^{2})^{2})-6m^{2}(\frac{\pi^{2}}{4}+2)+O(\varepsilon)$ (23)

Loop (III)

$I=-\frac{1}{\varepsilon}\,\frac{1}{4}(9j_{t}^{(0)}(s)+1)+6\left(-\frac{3}{16}+\frac{1}{8}j_{s}^{(1)}(s)+\left(\frac{3}{4}\log
m^{2}-\frac{17}{16}\right)j_{t}^{(0)}(s)+\frac{3}{8}j_{t}^{(1)}(s)\right)+O(\varepsilon)$
(24)

In this loop instance, the $O(\varepsilon^{-2})$ term happens to be vanishing.

Loop (IV)

$I=\frac{1}{\varepsilon^{2}}+\frac{1}{\varepsilon}(1-3j_{s}^{(1)}(s))+\left(\frac{3}{2}j_{s}^{(2)}(s)+3(j_{s}^{(1)}(s))^{2}-3j_{s}^{(1)}(s)-\frac{\pi^{2}}{2}+1\right)+O(\varepsilon)$
(25)

In the above formulae, we use the following functions where
$D(x)=m^{2}-sx\bar{x}$. 222 Here, $j_{t}^{(0)}=(2j_{s}^{(1)}+4-2\log
m^{2})/(4m^{2}-s)$ and $j_{t}^{(1)}=m^{2}(j_{s}^{(2)}+4j_{s}^{(1)}-(\log
m^{2})^{2})/(4m^{2}-s)$.

$j_{s}^{(k)}(s)=\int_{0}^{1}dx\,[\log D(x)]^{k},\quad
j_{t}^{(k)}(s)=\int_{0}^{1}dx\,\frac{m^{2}[\log D(x)]^{k}}{D(x)}.$ (26)

In the present approach, we sometimes encounter a pseudo-threshold. As an
example, in the calculation of the $O(1)$ term of Loop (III), the following
term appears:

$\int_{0}^{1}dx\,\frac{(2m^{2}-sx\bar{x})\log(2m^{2}-sx\bar{x})}{D(x)},$ (27)

which has the threshold at $s=8m^{2}$. While such terms exist in several
regions, in the sum they cancel each other out, so that the threshold at
$s=8m^{2}$ vanishes. In the calculation of the $O(\varepsilon)$ terms, we find
more complicated pseudo-thresholds. The numerical integration is done in each
region, so that the appearance of a pseudo-threshold might affect the high-
precision computation. Further study will be required to clarify this
behavior.

##  3 Numerical results

The numerical calculation of Eq. (18) is done in the following way.

We consider the coefficients in Eq. (7) and there are two approaches. In the
first, $\varepsilon$ no longer stands for an infinitesimal quantity but is a
(small) finite value. We prepare a descending positive-number sequence
$\\{\varepsilon_{j}\\}$ and calculate the $I(klm)$ taking
$\varepsilon=\varepsilon_{j}$ to obtain a value $I_{j}$. From the sequence of
values $\\{I_{j}\\}$, we estimate the $C^{\prime}s$ in Eq. (7) by a linear or
nonlinear extrapolation method [34, 35, 13]. If the singularity at $t=0$ or
$u=0$ is of type $\int_{0}^{1}dx\,x^{-1+\varepsilon}f(x)$ where $f$ is a non-
singular function, the calculation is just a numerical integration since
$\varepsilon$ is a positive finite quantity. If the behavior of the integrand
is of type $\int_{0}^{1}dx\,x^{-2+\varepsilon}f(x)$, we use the treatment
described in Appendix A. Furthermore, we can use the numerical integration
routine DQAGSE from QUADPACK [15, 12], which performs a built-in
extrapolation. Another approach is, after the separation of the UV singularity
of Eq. (18) as shown in Appendix B, by expanding the integrand in a Maclaurin
series in $\varepsilon$ and calculating each term using numerical integration
to obtain the coefficients of Eq. (7).

Above the threshold, the imaginary part appears in $I$. In this case, we no
longer treat the $\rho$ in Eq. (4) as infinitesimal. We prepare a geometric
sequence of numerical values ${\rho_{j}}$ and calculate the integrals $I(klm)$
with $\rho=\rho_{j}$ to obtain a value $I_{j}$. Then we perform the estimation
of $I$ as the limiting value for $\rho\rightarrow 0$. The extrapolations with
respect to $\rho$ and $\varepsilon$ can be combined to constitute a double
extrapolation [22, 23].

Figure 2: Expansion coefficients of Loop (I), $C_{-2}$ (dark-orange solid
line), $C_{-1}$ (green solid line), $C_{0}$ (red solid line) and $C_{1}$ (real
part in blue dots and imaginary part in blue square boxes) as a function of
$s$. The imaginary part of $C_{1}$ appears for $16m^{2}<s$ and it is
multiplied by 100 to fit in the plot.

In Figures 2 \- 5, we present $C_{-2},C_{-1},C_{0},C_{1}$ of Loop (I), (II),
(III), and (IV) as a function of $s$. Here we take the numerical value of
$m=1$, so that the horizontal axis shows $s$ in units of $m^{2}$. In these
figures, the first three coefficients are shown in lines and the values of
$C_{1}$ are shown in dots.

(a)  (b)

Figure 3: Expansion coefficients of Loop (II), $C_{-2}$ (dark-orange solid
line), $C_{-1}$ (green solid line), $C_{0}$ (red solid line) and $C_{1}$ (blue
dots) as a function of $s$. The coefficients around $s=4m^{2}$ are shown in
the small inserted plot. (a) Real part; (b) Imaginary part.

(a)  (b)

Figure 4: Expansion coefficients of Loop (III), $C_{-1}$ (green solid line),
$C_{0}$ (red solid line) and $C_{1}$ (blue dots) as a function of $s$. (a)
Real part, $C_{1}$ is multiplied by 0.5 to fit in the plot; (b) Imaginary
part.

(a)  (b)

Figure 5: Expansion coefficients of Loop (IV), $C_{-2}$ (dark-orange solid
line), $C_{-1}$ (green solid line), $C_{0}$ (red solid line) and $C_{1}$ (blue
dots) as a function of $s$. The coefficients around $s=4m^{2}$ are shown in
the small inserted plot. (a) Real part; (b) Imaginary part.

In Fig. 2, we show Loop (I) and a moderate four-body threshold at $s=16m^{2}$
can be seen. The first three coefficients are real. Except for Loop (I), we
can observe a significant threshold at $s=4m^{2}$, though the behavior depends
on each loop instance. In Loop (II) and (IV), the real part shows a finite
peak at $s=4m^{2}$, while it is divergent in Loop (III). Any other possible
threshold is numerically not very clear. In the figures for Loop (II), (III)
and (IV), the behavior near $s=4m^{2}$ is shown as the insertion.

For the first three terms in Eq. (7), the formulae in the previous section and
the numerical results agree well. By the numerical method we can compute the
fourth coefficient, $C_{1}$. The values obtained also agree with those by
pySecDec [36, 37, 38] in the range of $s$ shown in the figures. For example,
when $s=8m^{2}$, the number of digits in agreement for the real part of Loop
(I), (II), (III) and (IV) is 8, 4, 3 and 3, respectively. Similarly, for the
imaginary part, we obtain 6, 4, 3, and 4 digits, respectively. Particularly
for $s=m^{2}$, the coefficients of Loop (I), (II) and (IV) show good agreement
(10, 10 and 12 digits, respectively, for the $C_{1}$ term) with those by
Laporta [39].

##  4 Conclusion

In this paper, we have evaluated 3-loop two-point diagrams with reasonable
accuracy. The first four coefficients of the Laurent series in $\varepsilon$
are computed as a function of $s$ beyond the threshold. The coefficients of
the UV divergent terms are given as explicit analytic expressions. The
analysis uses a unique property, which we call “complete”, where the $f$
function factored from the $U$ function is common in all regions.

Many multi-loop integrals have no such property. Also, when the masses of the
internal lines are not the same, the integrals for each region should be
computed one by one. However, the method introduced in the current paper can
be applied to many diagrams and to general mass configurations. A preliminary
study of some other 3-loop two-point functions for $N\leq 6$ shows that they
can be handled by the method discussed in the text.

The numerical method developed by the authors is able to calculate the
coefficients of the Laurent series in the kinematical region both below and
above the threshold. The method shows good performance for the calculation of
the 3-loop integrals studied.

## Acknowledgement

We acknowledge the support by JSPS KAKENHI under Grant Numbers JP20K11858,
JP20K03941 and JP21K03541, and by the National Science Foundation, Award
Number 1126438 that funded initial work on numerical integration for the
ParInt package. We would like to express our deep gratitude to Dr. T. Kaneko
for showing us with a clear explanation how to handle the SD method.

## Appendix

##  Appendix A Separation of the UV singularity

We briefly explain how the UV singularity is extracted in Eq. (18). The basic
tool is a subtraction to separate the $\varepsilon$-pole. Suppose one
considers the integral

$J=\int_{0}^{1}dx\,\frac{1}{x^{1-\varepsilon}}f(x,\varepsilon)$ (28)

where $f(x,\varepsilon)$ is regular at $x=0$. In the following, we write
$f(x,\varepsilon)$ as $f(x)$ for brevity. Then the subtraction is done as

$J=\int_{0}^{1}dx\,\frac{1}{x^{1-\varepsilon}}[f(0)+(f(x)-f(0))]=\frac{1}{\varepsilon}f(0)+\int_{0}^{1}dx\,\frac{1}{x^{1-\varepsilon}}[f(x)-f(0)].$
(29)

If the integral is singular in two variables, the above procedure is repeated:

$\displaystyle J$ $\displaystyle=$
$\displaystyle\int_{0}^{1}dx\,\int_{0}^{1}dy\,\frac{1}{x^{1-\varepsilon}y^{1-\varepsilon}}f(x,y)=\frac{1}{\varepsilon^{2}}f(0,0)$
(30)
$\displaystyle+\frac{1}{\varepsilon}\int_{0}^{1}dx\,\frac{1}{x^{1-\varepsilon}}(f(x,0)-f(0,0))+\frac{1}{\varepsilon}\int_{0}^{1}dy\,\frac{1}{y^{1-\varepsilon}}(f(0,y)-f(0,0))$
$\displaystyle+\int_{0}^{1}dx\,\int_{0}^{1}dy\,\frac{1}{x^{1-\varepsilon}y^{1-\varepsilon}}(f(x,y)-f(x,0)-f(0,y)+f(0,0)).$

We also encounter an integral of the following type:

$J=\int_{0}^{1}dx\,\frac{1}{x^{2-\varepsilon}}f(x),$ (31)

which needs analytic continuation for the singularity at $x=0$. We use two
methods to handle the integral.

Method-1: Expansion and subtraction

$J=\frac{1}{-1+\varepsilon}f(0)+\frac{1}{\varepsilon}f^{\prime}(0)+\int_{0}^{1}dx\,\frac{1}{x^{2-\varepsilon}}(f(x)-f(0)-xf^{\prime}(0))$
(32)

Method-2: Integration by parts

$J=\left[\frac{x^{-1+\varepsilon}}{-1+\varepsilon}f(x)\right]_{0}^{1}-\int_{0}^{1}dx\,\frac{x^{-1+\varepsilon}}{-1+\varepsilon}f^{\prime}(x)$
$=\frac{1}{-1+\varepsilon}f(1)+\frac{1}{(1-\varepsilon)\varepsilon}f^{\prime}(0)+\frac{1}{1-\varepsilon}\int_{0}^{1}dx\,\frac{1}{x^{1-\varepsilon}}(f^{\prime}(x)-f^{\prime}(0))$
(33)

Numerical evaluation by the two methods gives the same result. From a
pedagogical point of view, this is an explicit example to demonstrate the
uniqueness of the analytic continuation.

##  Appendix B Expansion in $\boldsymbol{\varepsilon}$

In Eq. (29), Eq. (30), Eq. (32) and Eq. (33), the UV singularity is explicitly
shown as a factor $\varepsilon^{-1}$ or $\varepsilon^{-2}$ and the integrals
are no longer singular. Each term in the right-hand side can be expanded as a
series in $\varepsilon$ up to the required order. Then we obtain an expression
to be computed analytically or numerically.

As an example, explicit computation of Loop (IV) is shown here. We show the
way to compute $I(klm)$ in Eq. (7) up to the $O(\varepsilon)$. In the
following equations, terms beyond $O(\varepsilon)$ are omitted.

We use the following notations.

$H=\frac{G^{-3\varepsilon}}{f^{2-\varepsilon}},\quad
H_{a}=\left.H\right|_{t=0},\quad H_{b}=\left.H\right|_{u=0},\quad
H_{0}=\left.H\right|_{t=0,u=0}.$ (34) $f=1+v+uv+tuv,\quad
f_{a}=\left.f\right|_{t=0}=1+v+uv,\ f_{b}=\left.f\right|_{u=0},\
f_{0}=\left.f\right|_{t=0,u=0}=f_{b}=1+v~{},$ (35)

and similar notation of the suffix-$a,b,0$ is used for $G$ in Eq. (20). While
$f$ is common in all regions, $G$ depends on the region.

The notation for the integral is as follows.

$\int
d\Gamma_{X}=\int_{0}^{1}dw_{1}\int_{0}^{1}dw_{2}\int_{0}^{1}dt\int_{0}^{1}du\int_{0}^{1}dv,\quad\int
d\Gamma_{a}=\int_{0}^{1}dw_{1}\int_{0}^{1}dw_{2}\int_{0}^{1}du\int_{0}^{1}dv,$
$\int
d\Gamma_{b}=\int_{0}^{1}dw_{1}\int_{0}^{1}dw_{2}\int_{0}^{1}dt\int_{0}^{1}dv,\quad\int
d\Gamma_{0}=\int_{0}^{1}dw_{1}\int_{0}^{1}dw_{2}\int_{0}^{1}dv.$ (36)

$I(123)$ In this integral, $t$\- and $u$\- divergences appear. Here $h=t$,
$G_{a}=G_{0}=m^{2}-s\bar{w}_{1}w_{1}$ and
$G_{b}=m^{2}-s\bar{w}_{1}w_{1}+t(m^{2}-s\bar{w}_{2}w_{2})$.

$I(123)=I(123)_{0}+I(123)_{1}+I(123)_{2}+I(123)_{3}$ (37)

where $0,1,2,3$ denote the integrals of
$H_{0},H_{a}-H_{0},H_{b}-H_{0},H-H_{a}-H_{b}+H_{0}$, respectively.

$I(123)_{0}=\frac{1}{2\varepsilon^{2}}\int
d\Gamma_{0}\frac{1}{f_{0}^{2}}\Bigl{[}1+\varepsilon\left(\log f_{0}-3\log
G_{0}\right)+\frac{1}{2}\varepsilon^{2}\left(\log f_{0}-3\log
G_{0}\right)^{2}$ $+\frac{1}{6}\varepsilon^{3}\left(\log f_{0}-3\log
G_{0}\right)^{3}\Bigr{]}\,.$ (38)

$I(123)_{1}=\frac{1}{2\varepsilon}\int
d\Gamma_{a}\frac{1}{u}\Bigl{[}\Bigl{(}\frac{1}{f_{a}^{2}}-\frac{1}{f_{0}^{2}}\Bigr{)}$
$+\varepsilon\Big{(}\frac{1}{f_{a}^{2}}(\log u+\log f_{a}-3\log
G_{0})-\frac{1}{f_{0}^{2}}(\log u+\log f_{0}-3\log G_{0})\Big{)}$
$+\frac{1}{2}\varepsilon^{2}\Big{(}\frac{1}{f_{a}^{2}}(\log u+\log f_{a}-3\log
G_{0})^{2}-\frac{1}{f_{0}^{2}}(\log u+\log f_{0}-3\log
G_{0})^{2}\Big{)}\Bigr{]}\,.$ (39)

$I(123)_{2}=\frac{1}{\varepsilon}\int
d\Gamma_{b}\frac{1}{t}\frac{1}{f_{0}^{2}}\Bigl{[}(-3)\varepsilon(\log
G_{b}-\log G_{0})$ $+\varepsilon^{2}\Bigl{(}\frac{9}{2}((\log G_{b})^{2}-(\log
G_{0})^{2})-3(2\log t+\log f_{0})(\log G_{b}-\log G_{0})\Bigr{)}\Bigr{]}\,.$
(40)

$I(123)_{3}=\int
d\Gamma_{X}\frac{1}{tu}\Bigl{[}\Bigl{(}\frac{1}{f^{2}}-\frac{1}{f_{a}^{2}}\Bigr{)}+\varepsilon\Bigl{(}\frac{1}{f^{2}}(2\log
t+\log u+\log f-3\log G)$ $-\frac{1}{f_{a}^{2}}(2\log t+\log u+\log
f_{a}-3\log G_{0})+3\frac{1}{f_{0}^{2}}(\log G_{b}-\log
G_{0})\Bigr{)}\Bigr{]}\,.$ (41)

$I(132),I(134)$ These two integrals have divergences from $t=0$ but they have
no divergence from $u=0$. Here $h=tu,tuv$, where the factor $u$ cancels out
$1/u$ in Eq. (18), and $G_{a}=m^{2}-s\bar{w}_{1}w_{1}$.

$I(132,134)=I(132,134)_{0}+I(132,134)_{1}$ (42)

where $0,1$ indicate the integrals of $H_{a},H-H_{a}$, respectively.

$I(132,134)_{0}=\frac{1}{2\varepsilon}\int
d\Gamma_{a}\frac{(1,v)}{f_{a}^{2}}\Bigl{[}1+\varepsilon(\log u+\log
f_{a}-3\log G_{a})+\frac{1}{2}\varepsilon^{2}(\log u+\log f_{a}-3\log
G_{a})^{2}\Bigr{]}\,.$ (43)

$I(132,134)_{1}=\int
d\Gamma_{X}\frac{(1,v)}{t}\Bigl{[}\Bigl{(}\frac{1}{f^{2}}-\frac{1}{f_{a}^{2}}\Bigr{)}+\varepsilon\Bigl{(}\frac{1}{f^{2}}(2\log
t+\log u+\log f-3\log G)$ $-\frac{1}{f_{a}^{2}}(2\log t+\log u+\log
f_{a}-3\log G_{a})\Bigr{)}\Bigr{]}\,.$ (44)

$I(312),I(341),I(314)$ These three integrals have no divergent term. Here
$h=t^{2}u,t^{2}u^{2}v,t^{2}uv$.

$I(312,341,314)=\int
d\Gamma_{X}\frac{(1,uv,v)}{f^{2}}\left[1+\varepsilon(2\log t+\log u+\log
f-3\log G)\right]\,.$ (45)

The integrals presented above for the $O(\varepsilon^{-2})$,
$O(\varepsilon^{-1})$ and $O(1)$ terms are computed analytically to obtain Eq.
(25). The $O(\varepsilon)$ terms are evaluated numerically as shown in Section
3.

## References

  * [1] N. Nakanishi.  Graph Theory and Feynman Integrals.  Gordon and Breach, New York, 1971.
  * [2] T. Binoth and G. Heinrich.  An automized algorithm to compute infrared divergent multi-loop integrals.  Nucl. Phys. B, 585:741–759, 2000.  arXiv:hep-ph/0004013.
  * [3] T. Binoth and G. Heinrich.  Numerical evaluation of multi-loop integrals by sector decomposition.  Nucl. Phys. B, 680:375, 2004.  arXiv:hep-ph/0305234.
  * [4] T. Binoth and G. Heinrich.  Numerical evaluation of phase space integrals by sector decomposition.  Nucl. Phys. B, 693:134, 2004.  arXiv:hep-ph/0402265.
  * [5] G. Heinrich.  Sector Decomposition.  Int. J. Mod. Phys. A, 23:1457, 2008.  arXiv:0803.4177.
  * [6] G. Heinrich.  Collider physics at the precision frontier.  Physics Reports, 922:1, 2021.  arXiv:2009.00516.
  * [7] C. Bogner and S. Weinzierl.  Resolution of singularities for multi-loop integrals.  Comput. Phys. Commun., 178:596, 2008.  arXiv:0709.4092.
  * [8] C. Bogner and S. Weinzierl.  Blowing up Feynman integrals.  Nuclear Physics B (Proc. Suppl.), 183:256, 2008.  arXiv:0806.4307.
  * [9] A. V. Smirnov and M. Tentyukov.  Feynman integral evaluation by a sector decomposition approach (fiesta).  Comput. Phys. Commun., 180:735, 2009.  arXiv:0807.4129.
  * [10] A. V. Smirnov and V. A. Smirnov.  Hepp and speer sectors within modern strategies of sector decomposition.  JHEP, 05:004, 2009.  arXiv:0812.4700.
  * [11] T. Kaneko and T. Ueda.  A geometric method of sector decomposition.  Comput. Phys. Commun., 181:1352, 2010.  arXiv:0908.2897.
  * [12] E. de Doncker, J. Fujimoto, N. Hamaguchi, T. Ishikawa, Y. Kurihara, Y. Shimizu, and F. Yuasa.  Quadpack computation of Feynman loop integrals.  Journal of Computational Science (JoCS), 3(3):102–112, 2011.  DOI:10.1016/j.jocs.2011.06.003.
  * [13] E. de Doncker, F. Yuasa, K. Kato, T. Ishikawa, J. Kapenga, and O. Olagbemi.  Regularization with numerical extrapolation for finite and UV-divergent multi-loop integrals.  Comput. Phys. Commun., 224:164–185, 2018.  arXiv:1702.04904.
  * [14] E. de Doncker, A. Almulihi, and F. Yuasa.  High speed evaluation of loop integrals using lattice rules.  The Journal of Physics: Conf. Series (JPCS), IOP Series, 1085(052005), 2018.  http://iopscience.iop.org/article/10.1088/1742-6596/1085/5/052005.
  * [15] R. Piessens, E. de Doncker, C. W. Überhuber, and D. K. Kahaner.  QUADPACK, A Subroutine Package for Automatic Integration, volume 1 of Springer Series in Computational Mathematics.  Springer-Verlag, 1983.
  * [16] ParInt.  http://www.cs.wmich.edu/parint.
  * [17] O. E. Olagbemi and E. de Doncker.  Scalable algorithms for multivariate integration with ParAdapt and CUDA.  In Proc. 2019 Int. Conf. on Comp. Science and Comp. Intelligence. IEEE Computer Society, 2019.  https://american-cse.org/csci2019/pdfs/CSCI2019-14dQVW1stBtXVEInMQPd3t/558400a481/558400a481.pdf.
  * [18] M. Mori.  An IMT double exponential formula for numerical integration.  Publ. RIMS Kyoto Univ., 14:713–729, 1978.  https://doi.org/10.2977/prims/1195188835.
  * [19] I.H. Sloan and S. Joe.  Lattice Methods for Multiple Integration.  Oxford University Press, 1994.
  * [20] D. Shanks.  Non-linear transformations of divergent and slowly convergent sequences.  J. Math. and Phys., 34:1–42, 1955.
  * [21] P. Wynn.  On a device for computing the $e_{m}(s_{n})$ transformation.  Mathematical Tables and Aids to Computing, 10:91–96, 1956.
  * [22] E. de Doncker, F. Yuasa, T. Ishikawa, and K. Kato.  Loop integral computation in the euclidean or physical kinematical region using numerical integration and extrapolation.  In Proceedings of ACAT 2022 (21st International Workshop on Advanced Computing and Analysis Techniques in Physics Research), 2022.  To be published, https://indico.cern.ch/event/1106990/contributions/4997231/.
  * [23] E. de Doncker, F. Yuasa, T. Ishikawa, and K. Kato.  3-loop feynman diagram computations in quantum field theory using sector decomposition and numerical extrapolation.  In Proceedings of CCP2023 (34th IUPAP Conference on Computational Physics), 2023.  Submitted for publication.
  * [24] S. Bloch, M. Kerr, and P. Vanhove.  A feynman integral via higher normal functions.  Compos. Math., 151:2329, 2015.  arXiv:1406.2664.
  * [25] J. Broedel, C. Duhr, and N. Matthes.  Meromorphic modular forms and the three-loop equal-mass banana integral.  JHEP, 02:184, 2022.  arXiv:2109.15251.
  * [26] D. Kreimer.  Bananas: multi-edge graphs and their feynman integrals.  Lett. Math. Phys., 113:38, 2023.  arXiv:2202.05490.
  * [27] S. Pögel, X. Wang, and S. Weinzierl.  The three-loop equal-mass banana integral in $\epsilon$-factorised form with meromorphic modular forms.  JHEP, 2022, 2022.  arXiv:2207.12893.
  * [28] S. P. Maritin.  Evaluation of three-loop self-energy master integrals with four or five propagators.  Phys. Rev. D, 107:053005, 2023.  arXiv:2211.16539v2.
  * [29] S. P. Maritin.  Renormalized $\epsilon$-finite master integrals and their virtues: The three-loop self-energy case.  Phys. Rev. D, 105:056014, 2022.  arXiv:2112.07694v2.
  * [30] N. Nakanishi.  General integral formula of perturbation term in the quantized field theory.  Prog. Theor. Phys., 17:401–418, 1957.
  * [31] P. Cvitanović and T. Kinoshita.  Feynman-dyson rules in parametric space.  Phys. Rev., D10:3978, 1974.
  * [32] P. Cvitanović and T. Kinoshita.  New approach to the separation of ultraviolet and infrared divergences of feynman-parametric integrals.  Phys. Rev., D10:3991, 1974.
  * [33] P. Cvitanović and T. Kinoshita.  Sixth-order magnetic moment of the electron.  Phys. Rev., D10:4007, 1974.
  * [34] E. de Doncker, J. Fujimoto, N. Hamaguchi, T. Ishikawa, Y. Kurihara, M. Ljucovic, Y. Shimizu, and F. Yuasa.  Extrapolation algorithms for infrared divergent integrals, 2010.  arXiv:1110.3587.
  * [35] E. de Doncker, F. Yuasa, and Y. Kurihara.  Regularization of IR-divergent loop integrals.  Journal of Physics: Conf. Ser., 368(012060), 2012.  doi:10.1088/1742-6596/368/1/012060.
  * [36] S. Borowka, G. Heinrich, S. Jahn, P. Jones, S. M. Kerner, J. Schlenk, and T. Zirke.  pySecDec: A toolbox for the numerical evaluation of multi-scale integrals.  Comput. Phys. Commun., 222:313–326, 2018.  arXiv:1703.09692 [hep-ph].
  * [37] S. Borowka, G. Heinrich, S. Jahn, P. Jones, S. M. Kerner, and J. Schlenk.  A GPU compatible quasi-Monte Carlo integrator interfaced to pySecDec.  Comput. Phys. Commun., 240:120–137, 2019.  arXiv:1811.11720.
  * [38] G. Heinrich, P. Jones, S. , M. Kerner, V. Magerya, A. Olsson, and J. Schlenk.  Numerical scattering amplitudes with pySecDec.  Comput. Phys. Commun., 295, 2024.  arXiv:2305.19768.
  * [39] S. Laporta.  High-precision calculation of multi-loop Feynman integrals by difference equations.  Int. J. Mod. Phys. A, 15:5087–5159, 2000.  arXiv:hep-ph/0102033.