* Hook et al. (2004) Hook I. M., Jørgensen I., Allington-Smith J. R., Davies R. L., Metcalfe N., Murowinski R. G., Crampton D., 2004, PASP, 116, 425
  * Howard et al. (2010) Howard A. W., et al., 2010, ApJ, 721, 1467
  * Humphreys & Davidson (1994) Humphreys R. M., Davidson K., 1994, PASP, 106, 1025
  * Hurley et al. (2002) Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897
  * Iglesias-Marzoa et al. (2015) Iglesias-Marzoa R., López-Morales M., Jesús Arévalo Morales M., 2015, PASP, 127, 567
  * Israelian et al. (1999) Israelian G., Rebolo R., Basri G., Casares J., Martín E. L., 1999, Nature, 401, 142
  * Ivanova (2018) Ivanova N., 2018, ApJ, 858, L24
  * Ivanova et al. (2015) Ivanova N., Justham S., Podsiadlowski P., 2015, MNRAS, 447, 2181
  * Janssens et al. (2022) Janssens S., et al., 2022, A&A, 658, A129
  * Jayasinghe et al. (2021) Jayasinghe T., et al., 2021, MNRAS, 504, 2577
  * Jayasinghe et al. (2022) Jayasinghe T., Rowan D. M., Thompson T. A., Kochanek C. S., Stanek K. Z., 2022, arXiv e-prints, p. arXiv:2207.05086
  * Jeans (1924) Jeans J. H., 1924, MNRAS, 85, 2
  * Justham et al. (2006) Justham S., Rappaport S., Podsiadlowski P., 2006, MNRAS, 366, 1415
  * Justham et al. (2014) Justham S., Podsiadlowski P., Vink J. S., 2014, ApJ, 796, 121
  * Kaufer et al. (1999) Kaufer A., Stahl O., Tubbesing S., Nørregaard P., Avila G., Francois P., Pasquini L., Pizzella A., 1999, The Messenger, 95, 8
  * Kawahara et al. (2018) Kawahara H., Masuda K., MacLeod M., Latham D. W., Bieryla A., Benomar O., 2018, AJ, 155, 144
  * Khokhlov et al. (2018) Khokhlov S. A., et al., 2018, ApJ, 856, 158
  * Klencki et al. (2021) Klencki J., Nelemans G., Istrate A. G., Chruslinska M., 2021, A&A, 645, A54
  * Kochanek et al. (2017) Kochanek C. S., et al., 2017, PASP, 129, 104502
  * Kolbl et al. (2015) Kolbl R., Marcy G. W., Isaacson H., Howard A. W., 2015, AJ, 149, 18
  * Kremer et al. (2018) Kremer K., Ye C. S., Chatterjee S., Rodriguez C. L., Rasio F. A., 2018, ApJ, 855, L15
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kruse & Agol (2014) Kruse E., Agol E., 2014, Science, 344, 275
  * Kurucz (1970) Kurucz R. L., 1970, SAO Special Report, 309
  * Kurucz (1979) Kurucz R. L., 1979, ApJS, 40, 1
  * Kurucz (1993) Kurucz R. L., 1993, SYNTHE spectrum synthesis programs and line data
  * Lacy et al. (2020) Lacy M., et al., 2020, PASP, 132, 035001
  * Lam et al. (2022) Lam C. Y., et al., 2022, ApJ, 933, L23
  * Langer et al. (2020) Langer N., et al., 2020, A&A, 638, A39
  * Laycock et al. (2015) Laycock S. G. T., Maccarone T. J., Christodoulou D. M., 2015, MNRAS, 452, L31
  * Lejeune et al. (1997) Lejeune T., Cuisinier F., Buser R., 1997, A&AS, 125, 229
  * Lejeune et al. (1998) Lejeune T., Cuisinier F., Buser R., 1998, A&AS, 130, 65
  * Lennon et al. (2021) Lennon D. J., Dufton P. L., Villaseñor J. I., Evans C. J., Langer N., Saxton R., Monageng I. M., Toonen S., 2021, arXiv e-prints, p. arXiv:2111.12173
  * Maeder (1987) Maeder A., 1987, A&A, 178, 159
  * Mahy et al. (2022) Mahy L., et al., 2022, arXiv e-prints, p. arXiv:2207.07752
  * Mamajek & Hillenbrand (2008) Mamajek E. E., Hillenbrand L. A., 2008, ApJ, 687, 1264
  * Marigo et al. (2017) Marigo P., et al., 2017, ApJ, 835, 77
  * Marshall et al. (2008) Marshall J. L., et al., 2008, in McLean I. S., Casali M. M., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II. p. 701454 (arXiv:0807.3774), doi:10.1117/12.789972
  * Martin et al. (2005) Martin D. C., et al., 2005, ApJ, 619, L1
  * Mashian & Loeb (2017) Mashian N., Loeb A., 2017, MNRAS, 470, 2611
  * Masseron et al. (2014) Masseron T., et al., 2014, A&A, 571, A47
  * Masseron et al. (2016) Masseron T., Merle T., Hawkins K., 2016, BACCHUS: Brussels Automatic Code for Characterizing High accUracy Spectra, Astrophysics Source Code Library, record ascl:1605.004 (ascl:1605.004)
  * Masuda & Hotokezaka (2019) Masuda K., Hotokezaka K., 2019, ApJ, 883, 169
  * McClintock & Remillard (1986) McClintock J. E., Remillard R. A., 1986, ApJ, 308, 110
  * McClintock & Remillard (2006) McClintock J. E., Remillard R. A., 2006, in , Vol. 39, Compact stellar X-ray sources. pp 157–213
  * McClure & Woodsworth (1990) McClure R. D., Woodsworth A. W., 1990, ApJ, 352, 709
  * McMillan (2017) McMillan P. J., 2017, MNRAS, 465, 76
  * Merloni et al. (2003) Merloni A., Heinz S., di Matteo T., 2003, MNRAS, 345, 1057
  * Michaely & Perets (2016) Michaely E., Perets H. B., 2016, MNRAS, 458, 4188
  * Miller-Jones et al. (2021) Miller-Jones J. C. A., et al., 2021, Science, 371, 1046
  * Morris & Naftilan (1993) Morris S. L., Naftilan S. A., 1993, ApJ, 419, 344
  * Mroz et al. (2022) Mroz P., Udalski A., Gould A., 2022, arXiv e-prints, p. arXiv:2207.10729
  * Naoz et al. (2016) Naoz S., Fragos T., Geller A., Stephan A. P., Rasio F. A., 2016, ApJ, 822, L24
  * Narayan & Yi (1995) Narayan R., Yi I., 1995, ApJ, 452, 710
  * Ness et al. (2015) Ness M., Hogg D. W., Rix H. W., Ho A. Y. Q., Zasowski G., 2015, ApJ, 808, 16
  * Nidever et al. (2002) Nidever D. L., Marcy G. W., Butler R. P., Fischer D. A., Vogt S. S., 2002, ApJS, 141, 503
  * Padmanabhan et al. (2008) Padmanabhan N., et al., 2008, ApJ, 674, 1217
  * Parent et al. (2022) Parent E., et al., 2022, ApJ, 924, 135
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Penoyre et al. (2022) Penoyre Z., Belokurov V., Evans N. W., 2022, MNRAS, 513, 5270
  * Pfahl et al. (2003) Pfahl E., Rappaport S., Podsiadlowski P., 2003, ApJ, 597, 1036
  * Phinney & Verbunt (1991) Phinney E. S., Verbunt F., 1991, MNRAS, 248, 21P
  * Plez (2012) Plez B., 2012, Turbospectrum: Code for spectral synthesis, Astrophysics Source Code Library (ascl:1205.004)
  * Podsiadlowski et al. (2002) Podsiadlowski P., Rappaport S., Pfahl E. D., 2002, ApJ, 565, 1107
  * Podsiadlowski et al. (2003) Podsiadlowski P., Rappaport S., Han Z., 2003, MNRAS, 341, 385
  * Portegies Zwart et al. (1997) Portegies Zwart S. F., Verbunt F., Ergma E., 1997, A&A, 321, 207
  * Prochaska et al. (2020) Prochaska J. X., Hennawi J. F., Westfall K. B., Cooke R. J., Wang F., Hsyu T., Davies F. B., Farina E. P., 2020, arXiv e-prints, p. arXiv:2005.06505
  * Qian et al. (2008) Qian S. B., Liao W. P., Fernández Lajús E., 2008, ApJ, 687, 466
  * Qin et al. (2019) Qin Y., Marchant P., Fragos T., Meynet G., Kalogera V., 2019, ApJ, 870, L18
  * Quataert & Narayan (1999) Quataert E., Narayan R., 1999, ApJ, 520, 298
  * Raghavan et al. (2010) Raghavan D., et al., 2010, ApJS, 190, 1
  * Ramírez-Tannus et al. (2021) Ramírez-Tannus M. C., et al., 2021, A&A, 645, L10
  * Ramírez et al. (2012) Ramírez I., Fish J. R., Lambert D. L., Allende Prieto C., 2012, ApJ, 756, 46
  * Rasio et al. (2000) Rasio F. A., Pfahl E. D., Rappaport S., 2000, ApJ, 532, L47
  * Remillard & McClintock (2006) Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49
  * Rice & Brewer (2020) Rice M., Brewer J. M., 2020, ApJ, 898, 119
  * Rizzuto et al. (2013) Rizzuto A. C., et al., 2013, MNRAS, 436, 1694
  * Sahlmann (2019) Sahlmann J., 2019, Johannes-Sahlmann/pystrometry: Release for Zenodo, Zenodo, doi:10.5281/zenodo.3515526
  * Sahu et al. (2022) Sahu K. C., et al., 2022, ApJ, 933, 83
  * Sana et al. (2014) Sana H., et al., 2014, ApJS, 215, 15
  * Sana et al. (2017) Sana H., Ramírez-Tannus M. C., de Koter A., Kaper L., Tramper F., Bik A., 2017, A&A, 599, L9
  * Schoettler et al. (2019) Schoettler C., Parker R. J., Arnold B., Grimmett L. P., de Bruijne J., Wright N. J., 2019, MNRAS, 487, 4615
  * Shahaf et al. (2019) Shahaf S., Mazeh T., Faigler S., Holl B., 2019, MNRAS, 487, 5610
  * Shahaf et al. (2022) Shahaf S., Bashi D., Mazeh T., Faigler S., Arenou F., El-Badry K., Rix H.-W., 2022, arXiv e-prints, p. arXiv:2209.00828
  * Shao & Li (2019) Shao Y., Li X.-D., 2019, ApJ, 885, 151
  * Sharma et al. (2007) Sharma P., Quataert E., Hammett G. W., Stone J. M., 2007, ApJ, 667, 714
  * Sheinis et al. (2002) Sheinis A. I., Bolte M., Epps H. W., Kibrick R. I., Miller J. S., Radovan M. V., Bigelow B. C., Sutin B. M., 2002, PASP, 114, 851
  * Shenar (2022) Shenar T., 2022, arXiv e-prints, p. arXiv:2208.02614
  * Shenar et al. (2020) Shenar T., et al., 2020, A&A, 639, L6
  * Shenar et al. (2022a) Shenar T., et al., 2022a, Nature Astronomy,
  * Shenar et al. (2022b) Shenar T., et al., 2022b, arXiv e-prints, p. arXiv:2207.07674
  * Shikauchi et al. (2020) Shikauchi M., Kumamoto J., Tanikawa A., Fujii M. S., 2020, PASJ, 72, 45
  * Shikauchi et al. (2022) Shikauchi M., Tanikawa A., Kawanaka N., 2022, ApJ, 928, 13
  * Skrutskie et al. (2006) Skrutskie M. F., et al., 2006, AJ, 131, 1163
  * Skumanich (1972) Skumanich A., 1972, ApJ, 171, 565
  * Smith (2014) Smith N., 2014, ARA&A, 52, 487
  * Sukhbold et al. (2016) Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H. T., 2016, ApJ, 821, 38
  * Takeda et al. (2013) Takeda Y., Honda S., Ohnishi T., Ohkubo M., Hirata R., Sadakane K., 2013, PASJ, 65, 53
  * Tanikawa et al. (2022) Tanikawa A., Hattori K., Kawanaka N., Kinugawa T., Shikauchi M., Tsuna D., 2022, arXiv e-prints, p. arXiv:2209.05632
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration et al., 2021, arXiv e-prints, p. arXiv:2111.03606
  * Thompson et al. (2019) Thompson T. A., et al., 2019, Science, 366, 637
  * Trimble & Thorne (1969) Trimble V. L., Thorne K. S., 1969, ApJ, 156, 1013
  * Valenti & Fischer (2005) Valenti J. A., Fischer D. A., 2005, ApJS, 159, 141
  * Vernet et al. (2011) Vernet J., et al., 2011, A&A, 536, A105
  * Vogt et al. (1994) Vogt S. S., et al., 1994, in Crawford D. L., Craine E. R., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 2198, Instrumentation in Astronomy VIII. p. 362, doi:10.1117/12.176725
  * Wang et al. (2021) Wang S., et al., 2021, Research in Astronomy and Astrophysics, 21, 292
  * Webbink (1984) Webbink R. F., 1984, ApJ, 277, 355
  * Wiktorowicz et al. (2020) Wiktorowicz G., Lu Y., Wyrzykowski Ł., Zhang H., Liu J., Justham S., Belczynski K., 2020, ApJ, 905, 134
  * Wong et al. (2022) Wong K. W. K., Breivik K., Farr W. M., Luger R., 2022, arXiv e-prints, p. arXiv:2206.04062
  * Woosley & Heger (2006) Woosley S. E., Heger A., 2006, ApJ, 637, 914
  * Woosley et al. (1995) Woosley S. E., Langer N., Weaver T. A., 1995, ApJ, 448, 315
  * Wright et al. (2010) Wright E. L., et al., 2010, AJ, 140, 1868
  * Yee et al. (2017) Yee S. W., Petigura E. A., von Braun K., 2017, ApJ, 836, 77
  * Yoon et al. (2006) Yoon S. C., Langer N., Norman C., 2006, A&A, 460, 199
  * van de Kamp (1975) van de Kamp P., 1975, ARA&A, 13, 295

##  Appendix A Observations and data reduction

###  A.1 MagE

We observed Gaia BH1 4 times using the Magellan Echellette spectrograph (MagE;
Marshall et al., 2008) on the 6.5m Magellan Clay telescope at Las Campanas
Observatory. All observations were carried out with the 0.7 arcsec slit and
exposure times ranging from 300 to 600 seconds. This yielded a typical SNR of
100 per pixel at 5,500 Å, spectral resolution $R\approx 5,400$, and wavelength
coverage of $3,500-11,000$ Å. To minimize drift in the wavelength solution due
to flexure, we obtained a ThAr arc on-sky immediately after each science
exposure.

We reduced the spectra using PypeIt (Prochaska et al., 2020), which performs
bias and flat field correction, cosmic ray removal, wavelength calibration,
sky subtraction, extraction of 1d spectra, and heliocentric RV corrections. To
assess the stability of the wavelength solution, we observed several standard
stars, measured their RVs using the same procedure we used for the science
targets, and compared the results to their measured Gaia DR3 RVs. We found
that the MagE-measured RVs were $1.5\pm 2\,\rm km\,s^{-1}$ larger than those
measured by Gaia, so we adopted a conservative uncertainty of $3\,\rm
km\,s^{-1}$ for all the MagE RVs.

###  A.2 GMOS

We observed Gaia BH1 4 times using the Gemini Multi-Object Spectrograph (GMOS;
Hook et al., 2004) on the 8.1m Gemini-North telescope on Mauna Kea (program
GN-2022B-DD-202). We used the R831_G5302 grating with a 0.5 arcsecond slit and
central wavelength of 5200 Å, leading to wavelength coverage from 4000 to 6400
Å and resolution $R\approx 4,400$. We obtained a CuAr arc on-sky immediately
after each science exposure. We used 342 second exposures, yielding a typical
SNR of 90 per pixel at 6000 Å.

We reduced the data using PypeIt, which required construction of a new
template for the R831 grating at blue wavelengths. Cross-correlation with a
template yielded formal RV uncertainties of order $1\,\rm km\,s^{-1}$, but we
found variations of up to $10\,\rm km\,s^{-1}$ between RVs measured at the
blue and red ends of the spectrum. We adopted a RV uncertainty of $4\,\rm
km\,s^{-1}$ for the GMOS RVs, which was also validated by observations of Gaia
BH1 with higher-resolution instruments that occurred nearly concurrently with
some of the GMOS observations.

###  A.3 X-shooter

We observed Gaia BH1 4 times using the X-shooter spectrograph (Vernet et al.,
2011) on the 8.2m UT3 telescope at the VLT on Cerro Paranal (program
2109.D-5047). We used the 0.5, 0.4, and 0.4 arcsec slits on the UVB, VIS, and
NIR arms. This setup yielded a resolution of $R\approx 9,700$ in the UVB arm,
$R\approx 18,400$ in the VIS arm, and $R\approx 11,600$ in the NIR arm, with
near-continuous wavelength coverage from 3100-24000 Å. We reduced the spectra
using pypeit.

We generally measured RVs using the NIR spectra, whose wavelength solution is
calculated from telluric OH lines in the science spectra, yielding a maximally
stable wavelength solution. For the first epoch, we inadvertently used a 0.9
arcsec slit in the NIR arm, yielding a lower spectral resolution. We measured
the RV from the VIS arm for this epoch. We adopt an uncertainty of $1\,\rm
km\,s^{-1}$ on all the X-shooter RVs.

###  A.4 HIRES

We observed Gaia BH1 9 times using the High Resolution Echelle Spectrometer
(HIRES; Vogt et al., 1994) on the 10m Keck I telescope on Maunakea. The data
were obtained and reduced using the standard California Planet Survey setup
(CPS; Howard et al., 2010), including use of the C2 decker (0.86 arcseconds
$\times$ 14 arcseconds), which yields spectra with $R\approx 55,000$ and
wavelength coverage over most of 3700–8000 Å, with 2 gaps 100–200 Å wide. We
used 600 second exposures, yielding SNR of order 40 per pixel at 6000 Å. The
CPS reduction includes sky-subtraction using the long C2 decker.

As part of the CPS pipeline, absolute RVs are measured using the method
described by Chubak et al. (2012): to correct for instrumental shifts in the
wavelength zero-point, the telluric A and B molecular oxygen absorption bands
of each star are cross-correlated with those of RV standard stars, whose
adopted RVs are matched to the scale of Nidever et al. (2002) and other RV
standard catalogs. The resulting RVs are robust to $0.1\,\rm km\,s^{-1}$. We
use these RVs throughout our analysis. Most of our observations were performed
through an iodine cell, making it possible to measure higher-precision RVs
from these data in the future, after a higher-SNR template spectrum has been
obtained and the method of Butler et al. (1996) can be applied.

###  A.5 FEROS

We observed Gaia BH1 17 times with the Fiberfed Extended Range Optical
Spectrograph (FEROS; Kaufer et al., 1999) on the 2.2m ESO/MPG telescope at La
Silla Observatory. We used $2\times 2$ binning in order to reduce readout
noise (allowing us to measure RVs down to $G\approx 16$); this resulted in an
effective resolution $R\approx 37000$ measured from sky lines. We took 2400
second exposures and achieved a typical SNR of 20 per pixel. We reduced the
data with the standard ESO MIDAS pipeline, which performs bias-subtraction,
flat fielding, wavelength calibration, and sky subtraction, and merges the
spectra from different orders.

On some nights, the wavelength solution drifted from the one established by
the afternoon arcs by up to a few $\rm km\,s^{-1}$. We corrected this by
cross-correlating an arc spectrum (which is taken during each exposure with a
second fiber) with a reference arc spectrum taken on the first night of our
observations. Based on observations of an RV standard with the same setup, we
adopted and uncertainty of $0.3\,\rm km\,s^{-1}$ for all the FEROS
observations, except in one case where poor conditions yielded a formal
uncertainty of $0.5\,\rm km\,s^{-1}$. Near-concurrent FEROS and HIRES
observations of the source allowed us to verify that the two RV zeropoints are
consistent within this tolerance.

###  A.6 ESI

We observed Gaia BH1 one time with the Echellette Spectrograph and Imager
(ESI; Sheinis et al., 2002) on the 10m Keck II telescope on Maunakea. We used
a 300 second exposure with the 0.3 arcsec slit, yielding a resolution
$R\approx 12000$ and SNR = 60, with useful wavelength coverage of 3900-10000
Å. We reduced the data using the MAuna Kea Echelle Extraction (MAKEE)
pipeline, which performs bias-subtraction, flat fielding, wavelength
calibration, and sky subtraction, and we refined the wavelength solution using
sky lines. Based on RV measurements of RV standards on the same night, we
adopted an uncertainty of $1.5\,\rm km\,s^{-1}$.

##  Appendix B Gaia observations

As discussed in Section 4.2.1, the astrometric uncertainties for the
parameters describing the orbit of Gaia BH1 are unusually large for the star’s
apparent magnitude. Here we show that this is likely a consequence of the Gaia
scanning law.

Epoch-level astrometric data is not published in DR3, but one can determine
when a source was observed using the Gaia observation scheduling tool
(GOST).777https://gaia.esac.esa.int/gost/ Given a source position, GOST
returns a list of observation times when the scanning law predicts that a
source will transit across the Gaia focal plane. It is not guaranteed that
Gaia actually obtains data on each source at the predicted times, as the tool
does not account for gaps between CCDs and issues that cause temporary gaps in
the datastream, such as micrometeor impacts. However, for Gaia BH1, GOST
predicts 21 visibility periods (i.e., groups of observations separated from
other groups by at least 4 days), and the gaia_source table reports that 21
visibility periods were indeed used in the astrometric solution. This suggests
that the scan times predict by GOST are likely a good approximation of
reality.

We visualize these times in Figure 10. It is clear that although Gaia DR3
covered more than 5 orbital period of Gaia BH1, the phase coverage is far from
uniform, and only about half of the orbit in phase was covered by astrometric
data. This may be a consequence of the fact that the orbital period is a near-
integer multiple of the Gaia precession period, and close to half a year. Our
RV follow-up shows that despite this limitation, Gaia was able to measure an
astrometric orbit that was basically correct. However, the incomplete phase
coverage make it unsurprising that the uncertainties in the astrometric
parameters are larger than usual for a source of this brightness.

GOST predicts that the north side of the orbit will remain poorly sampled even
in Gaia DR4. We note that a few strategically timed observations with HST,
JWST, GRAVITY+, or MICADO could fill in this part of the orbit and yield
improved astrometric precision if analyzed together with the epoch-level Gaia
astrometric data that will become available in DR4.

Figure 10: Predicted observation times of Gaia BH1 from the Gaia observation
scheduling tool (GOST). Black line shows the best-fit photocenter orbit. Red
points show the predicted photocenter positions at the times when GOST
predicts Gaia would have observed the source. Note that we do not have access
to the actual measured $\Delta\rm RA$ and $\Delta\rm Dec$ values; only to the
predicted scan times. GOST predicts that Gaia observed the source in 21
visibility periods which all fall on one side of the orbit. This may be a
result of the fact that the orbital period is a near-integer multiple of
Gaia’s 63-day precession period, and almost exactly half a year.

##  Appendix C Comparison of the astrometric and RV solutions

As discussed in Section 4, the tightest constraints on the Gaia BH1 orbit can
be obtained through joint fitting of the RVs and astrometric constraints.
However, it is also useful to compare the RVs predicted by the Gaia solution
alone to the measured RVs, in order to assess whether or not the astrometric
solution is consistent with the RVs. Figure 11 compares the measured RVs to
predictions based on the Gaia astrometric solution, which are generated from
posterior draws of the fit without RVs. Because the astrometric solution does
not constrain the center-of-mass velocity, we set it to $\rm 45\,km\,s^{-1}$
(consistent with the measured value) in all cases.

The measured RVs are basically consistent with the Gaia-only solution, in the
sense that some posterior draws go through all the RVs. Figure 12 and Table 3
also show that all inferred parameters are consistent within $\lesssim
2\sigma$. This suggests that the Gaia solution and its uncertainties are
reliable, despite the incomplete orbit coverage discussed in Appendix B.

Parameter | units | Astrometry-only constraint | Joint RVs+astrometry constraint | How many sigma discrepant?
---|---|---|---|---
a_thiele_innes | mas | $-0.26\pm 0.17$ | $-0.01\pm 0.05$ | 1.1
b_thiele_innes | mas | $2.93\pm 0.18$ | $2.62\pm 0.02$ | 1.6
f_thiele_innes | mas | $1.52\pm 0.16$ | $1.61\pm 0.02$ | 0.5
g_thiele_innes | mas | $0.53\pm 0.55$ | $-0.36\pm 0.05$ | 1.5
period | days | $185.77\pm 0.31$ | $185.59\pm 0.05$ | 0.5
eccentricity |  | $0.49\pm 0.07$ | $0.45\pm 0.01$ | 0.5
t_periastron | days | $-12.0\pm 6.3$ | $-1.1\pm 0.7$ | 1.6
Table 3: Comparison of the Gaia astrometric constraints on Gaia BH1’s
photocenter ellipse to constraints from the joint RVs+astrometry fit. All
parameters are consistent at the 1.6$\sigma$ level.  Figure 11: Comparison of
measured RVs (same as Figure 2) to predictions from the Gaia astrometric
solution alone. The center-of-mass RV (which is not constrained by astrometry)
is set to the value inferred from RVs, but the cyan lines otherwise do not
know about the measured RVs. The uncertainty in the predictions is larger than
when RVs are included in the fit, but the measured RVs are generally
consistent with the predictions, validating the Gaia solution.  Figure 12:
Comparison between constraints obtained from joint fitting of RVs and
astrometry (black) and constraints from astrometry alone (cyan). The two sets
of constraints are generally consistent, but those with RVs are tighter.
Numbers on the diagonals reflect the joint astrometry + RVs constraints. Both
sets of constraints are listed in Table 1.

##  Appendix D Radial velocities

Radial velocities are listed in Table 4. All times are reported as mid-
exposure times in the heliocentric frame (HJD UTC). RVs are also in the
heliocentric frame and their uncertainties include uncertainties in wavelength
calibration. The reported SNR is at the wavelength range we used for RV
measurement, which in most cases was near $5400$ Å.

HJD UTC | Radial velocity ($\rm km\,s^{-1}$) | Instrument | Resolution $R$ | SNR
---|---|---|---|---
2457881.2899 | $20.0\pm 4.0$ | LAMOST | 1800 | 6
2458574.3663 | $8.9\pm 5.6$ | LAMOST | 1800 | 55
2459767.6226 | $63.8\pm 3$ | MagE | 5400 | 110
2459791.9186 | $131.90\pm 0.1$ | HIRES | 55000 | 40
2459795.6461 | $141.4\pm 3$ | MagE | 5400 | 100
2459796.4995 | $142.7\pm 3$ | MagE | 5400 | 100
2459798.8399 | $140.6\pm 4$ | GMOS | 4300 | 80
2459805.5101 | $127.7\pm 1.0$ | X-shooter | 18400 | 55
2459808.7388 | $118.0\pm 4$ | GMOS | 4300 | 75
2459813.6045 | $90.5\pm 0.3$ | FEROS | 37000 | 20
2459814.5874 | $86.1\pm 0.3$ | FEROS | 37000 | 15
2459815.5927 | $81.9\pm 0.3$ | FEROS | 37000 | 19
2459817.5278 | $74.5\pm 0.3$ | FEROS | 37000 | 20
2459818.5266 | $71.0\pm 0.3$ | FEROS | 37000 | 20
2459818.7870 | $67.8\pm 4$ | GMOS | 4300 | 85
2459819.5543 | $67.0\pm 0.3$ | FEROS | 37000 | 20
2459820.5465 | $64.0\pm 0.3$ | FEROS | 37000 | 21
2459821.5669 | $60.7\pm 0.3$ | FEROS | 37000 | 22
2459822.5745 | $57.8\pm 0.3$ | FEROS | 37000 | 22
2459823.5422 | $54.8\pm 0.3$ | FEROS | 37000 | 21
2459823.8525 | $53.76\pm 0.1$ | HIRES | 55000 | 40
2459824.5305 | $52.1\pm 0.3$ | FEROS | 37000 | 18
2459824.8516 | $51.18\pm 0.1$ | HIRES | 55000 | 40
2459825.5361 | $49.8\pm 0.3$ | FEROS | 37000 | 20
2459826.7920 | $46.59\pm 0.1$ | HIRES | 55000 | 40
2459828.5677 | $42.2\pm 0.3$ | FEROS | 37000 | 20
2459829.5373 | $42.1\pm 3$ | MagE | 5400 | 80
2459829.5768 | $40.4\pm 0.3$ | FEROS | 37000 | 17
2459830.6452 | $38.5\pm 0.3$ | FEROS | 37000 | 10
2459831.6223 | $36.5\pm 0.5$ | FEROS | 37000 | 5
2459833.7523 | $33.23\pm 0.1$ | HIRES | 55000 | 40
2459834.5509 | $32.3\pm 0.3$ | FEROS | 37000 | 11
2459834.7691 | $31.74\pm 0.1$ | HIRES | 55000 | 40
2459835.7678 | $30.14\pm 0.1$ | HIRES | 55000 | 40
2459838.8082 | $26.35\pm 0.1$ | HIRES | 55000 | 40
2459838.7208 | $27.5\pm 4$ | GMOS | 4300 | 78
2459840.7729 | $24.20\pm 0.1$ | HIRES | 55000 | 40
2459845.5069 | $19.4\pm 1.0$ | X-shooter | 11600 | 45
2459855.5012 | $14.2\pm 1.0$ | X-shooter | 11600 | 45
2459868.5128 | $9.3\pm 1.0$ | X-shooter | 11600 | 40
2459877.6978 | $10.5\pm 1.5$ | ESI | 12000 | 60
Table 4: Radial velocities.

##  Appendix E Other candidates

Figure 13: Parameter space within which we selected BH candidates. Left panel
shows a simulation, in which we model a population of luminous (main-sequence)
binaries, triples, and compact object + normal star binaries subject to the
cuts imposed on binary solutions published in DR3. The relative number of
objects in the three populations is chosen arbitrarily. Binaries containing
dark companions fall at larger separations at fixed period than binaries
containing two (or three) luminous stars. Right panel shows the same
quantities for binaries published in DR3. Red star shows Gaia BH1; black
points show other initially promising candidates, which we discuss in the
text.

Figure 13 illustrates our strategy for selecting BH candidates among the Gaia
astrometric solutions. We show that binaries containing a BH and a low-mass
main-sequence star are expected to have much larger photocenter wobble than
binaries containing two luminous stars.

The left panel shows a simulation. Here we begin with all sources in the
gaia_source catalog with parallax_over_error > 50\. We then attached a
simulated binary population to these sources, drawing periods from the
canonical lognormal distribution for solar-type stars (Raghavan et al., 2010)
and mass ratios from a uniform distribution. For each simulated binary, we
then calculate the expected photocenter semi-major axis and its uncertainty,
and we retain only the binaries that satisfy the cuts imposed on the binary
solutions published in DR3. We also generate a population of hierarchical
triples, drawing both the inner and outer mass ratios from uniform
distributions. We calculate the flux ratios using the $G-$band mass-magnitude
relation tabulated by Janssens et al. (2022). Finally, we generate a
population of BH + normal star and neutron star + normal binaries, by
assigning $10\,M_{\odot}$ and $1.4\,M_{\odot}$ dark companions to a random
subset of the luminous sources, again drawing periods from the canonical
lognormal. This is not meant to represent a realistic distribution of orbital
periods, masses, or mass ratios, but only to show the expected photocenter
orbit size at a given period.

In the right panel, we show the same observables for all objects in DR3 with
astrometric solutions (Orbital and AstroSpectroSB1 solutions). As expected,
most solutions fall in the region of parameter space that we expect to be
populated by luminous binaries and triples. For each binary, we calculate
$\tilde{M}$, a quantity which is related to, but not equal to, the total
dynamical mass:

$\tilde{M}=\frac{4\pi^{2}}{P_{{\rm
orb}}^{2}G}\left(\frac{a_{0}}{\varpi}\right)^{3}.$ (13)

For a binary with a luminous star and a dark companion,
$\frac{a_{0}}{\varpi}=a\frac{q}{1+q}$, where $q=M_{2}/M_{\star}$ and $a$ is
the semi-major axis. Thus, $\tilde{M}=M_{{\rm
tot}}\left(\frac{q}{1+q}\right)^{3}$. For BH companions to low-mass stars,
$q/(1+q)\approx 1$, and thus $\tilde{M}$ is only slightly less than $M_{\rm
tot}$. On the other hand, for a binary with two luminous stars,
$\frac{a_{0}}{\varpi}=a\delta_{q\ell}$, where

$\displaystyle\delta_{q\ell}=\left|\frac{q/\left(1+\ell\right)-\ell/\left(1+\ell\right)}{\left(1+q\right)}\right|,$
(14)

and $\ell=F_{2}/F_{1}$ is the ratio of the two stars’ $G-$band fluxes (e.g.
Penoyre et al., 2022). For binaries with two main-sequence stars,
$\delta_{q\ell}\lesssim 0.35$ (e.g. Shahaf et al., 2019), and thus in general
$\tilde{M}\lesssim 0.05M_{\rm tot}$. $\tilde{M}$ is thus a good quantity with
which to select massive dark companions to low-mass stars, though more
sophisticated methods exist that are sensitive to lower-mass BHs and BH
companions to massive stars (e.g. Shahaf et al., 2019, 2022).

With this in mind, we selected objects from Gaia DR3 with
$\tilde{M}>3\,M_{\odot}$. This yielded 6 objects, which are highlighted in the
right panel of Figure 13 and enumerated in Table 5. Sources brighter than
$G_{\rm RVs}=12$ have multi-epoch RVs measured (but not yet published) by
Gaia. The min-to-max range of these RVs is published as rv_amplitude_robust.
Comparing this to the expected full RV amplitude implied by the Gaia solution
provides a useful diagnostic of whether the solution is likely to be correct.

Each solution is accompanied by a goodness_of_fit value. As discussed by
Halbwachs et al. (2022), this quantity is formally expected to follow a normal
distribution, $\mathcal{N}(0,1)$, with values significantly larger than 1
indicative of a problematic solution. However, this may not hold in practice,
since the uncertainties in epoch-level measurements can be underestimated. We
find that the distribution of goodness_of_fit values displays a sharp
discontinuity at $G=13$. Fainter values approximately follow a
$\mathcal{N}(0,1.5)$ distribution, such that a value much larger than 1 is
indeed unusual. However, for brighter solutions, the distribution peaks at
$\approx 3.5$, and 13% of solutions of have goodness_of_fit > 10\. We thus
expect significantly larger goodness_of_fit values to be (potentially)
consistent with a good solution for bright sources.

Gaia DR3 source ID | $G$ | $\tilde{M}$ | $P_{\rm orb}$ | $(a_{0}\times d)$ | GoF | rv_amplitude_robust | expected RV amplitude | Verdict
---|---|---|---|---|---|---|---|---
| [mag] | $[M_{\odot}]$ | $[\rm days]$ | $[\rm au]$ |  | [$\rm km\,s^{-1}$] | [$\rm km\,s^{-1}$] |
3640889032890567040 | 9.2 | $122\pm 47$ | $1076\pm 12$ | $10.2\pm 1.3$ | 10.3 | 12.9 | $674\pm 154$ | ✗
4467000291193143808 | 15.5 | $119\pm 71$ | $1647\pm 520$ | $13.4\pm 2.2$ | 4.9 |  | $358\pm 182$ | ✗
4373465352415301632 | 13.8 | $11.5\pm 2.5$ | $185.8\pm 0.3$ | $1.44\pm 0.10$ | 0.3 |  | $165\pm 20$ | ✓
6281177228434199296 | 11.3 | $10.8\pm 1.6$ | $153.9\pm 0.4$ | $1.24\pm 0.06$ | 8.0 | 20.1 | $171\pm 10$ | ✗
5870569352746779008 | 12.3 | $6.7\pm 0.5$ | $1352\pm 45$ | $4.52\pm 0.13$ | 3.1 | 37.0 | $50\pm 2$ |
3664684869697065984 | 11.6 | $3.5\pm 2.6$ | $1220\pm 233$ | $3.35\pm 1.13$ | 3.6 | 18.3 | $54\pm 14$ | ✗
Table 5: BH candidates highlighted in Figure 13. GoF is the goodness_of_fit
parameter, for which large values indicate a potentially problematic solution
(see text for details). rv_amplitude_robust is the min-to-max range in the
individual-epoch Gaia RV measurements (only available for bright sources);
“expected RV amplitude” is the min-to-max RV amplitude predicted by the
astrometric solution. Sources are discussed individually in the text.

We comment briefly on each source:

3640889032890567040: Poor GoF and implausibly high companion mass. Ruled out
by low rv_amplitude_robust.

4467000291193143808: Period of 1647 days is significantly longer than Gaia DR3
baseline. Implausibly high companion mass; poor GoF for the source magnitude.

4373465352415301632: The subject of this paper.

6281177228434199296: Ruled out by our RV follow-up. Low rv_amplitude_robust
also makes the astrometric solution implausible.

5870569352746779008: Most plausible candidate beside Gaia BH1, but still
uncertain. $P_{\rm orb}$ is longer than the Gaia DR3 observing baseline. The
GoF is significantly larger than 1, but is not unusual for sources in the
relevant range of $G$ magnitude. The CMD position is consistent with a lower
giant, or possibly a red clump star. The inclination is relatively low ($35\pm
1$ deg), so the nature of the companion hinges on the correctness of the
astrometric solution. That is, if the Gaia solution is correct, RVs alone will
only yield a spectroscopic mass function near $1\,M_{\odot}$, which is not
high enough to unambiguously rule out a luminous companion to a red giant. The
reliability of the astrometric solution will be easier to assess in DR4, when
epoch-level astrometric data are published and a full orbit has been covered.
The measured rv_amplitude_robust is not too far from the expected peak-to-peak
amplitude. We have initiated a spectroscopic follow-up campaign, on which we
will report in future work. Thus far, the RVs are consistent with the
AstroSpectroSB1 solution. Further RV monitoring is warranted and necessary.
This source was recently investigated in detail by Tanikawa et al. (2022), who
found it to be the only credible BH candidate in Gaia DR3 with both
astrometric and spectroscopic Gaia data, but did not carry out follow-up. If
the solution turns out to be credible, the source will present a similar
puzzle to Gaia BH1: the BH progenitor would likely not have fit inside the
current orbit, but the orbit is much too wide to be a result of common
envelope evolution. The current orbital separation is such that if the
luminous star is a core helium burning red clump star, mass transfer via winds
and/or stable Roche lobe overflow would have occurred prior to the helium
flash.

3664684869697065984: Ruled out by APOGEE RVs.

For a $10\,M_{\odot}$ BH companion, the cut of $\tilde{M}>3\,M_{\odot}$ would
be satisfied for luminous stars with $M_{\star}\lesssim 8\,M_{\odot}$. For
higher-mass luminous stars, the AMRF, $\mathcal{A}$ statistic designed by
Shahaf et al. (2019) is more sensitive. We verified that no objects in the
gaiadr3.nss_two_body_orbit catalog with photometric masses above
$5\,M_{\odot}$ have $\mathcal{A}>0.6$, as would be required to identify a dark
companion with high confidence.

##  Appendix F HIRES Cannon labels

The labels inferred by the Cannon fit to the HIRES spectrum of Gaia BH1 are
listed in Table 6. The method is not expected to yield reliable formal
uncertainties on labels, because at high SNR the errors are dominated by
systematics in the spectral model rather than photon noise. We thus report the
cross-validation error achieved by Rice & Brewer (2020) for each label when
testing the model, which provides a rough lower limit on the uncertainties in
each label.

The atmospheric parameters and abundances reported by the Cannon are quite
similar to those we find with BACCHUS (Section 5.1), which we adopt as our
fiducial values. Note that the Cannon reports [X/H], while for the BACCHUS fit
we reported [X/Fe] = [X/H] - [Fe/H].

Parameter | Best-fit | Cross-validation error
---|---|---
$T_{\rm eff}\,[\rm K]$ | 5863 | 56
$\log\left(g/{\rm cm\,s^{-2}}\right)$ | 4.43 | 0.09
$v\sin i\,\left[{\rm km\,s^{-1}}\right]$ | 1.3 | 0.87
$\rm[C/H]$ | -0.24 | 0.05
$\rm[N/H]$ | -0.37 | 0.08
$\rm[O/H]$ | -0.10 | 0.07
$\rm[Na/H]$ | -0.33 | 0.05
$\rm[Al/H]$ | -0.28 | 0.04
$\rm[Mg/H]$ | -0.23 | 0.04
$\rm[Si/H]$ | -0.26 | 0.03
$\rm[Ca/H]$ | -0.26 | 0.03
$\rm[Ti/H]$ | -0.21 | 0.04
$\rm[V/H]$ | -0.22 | 0.06
$\rm[Cr/H]$ | -0.30 | 0.04
$\rm[Mn/H]$ | -0.50 | 0.05
$\rm[Fe/H]$ | -0.29 | 0.03
$\rm[Ni/H]$ | -0.34 | 0.04
$\rm[Y/H]$ | -0.31 | 0.08
Table 6: Parameters of the G star, from the HIRES cannon. The cross-validation
error is taken from Rice & Brewer (2020) and represents an estimate of the
uncertainty in each label.

##  Appendix G Dependence of Gaia astrometric uncertainties on observables

In Figure 14, we show the mean uncertainty in $a_{0}$ and $\varpi$ for DR3
astrometric binary solutions as a function of apparent magnitude and $a_{0}$.
We only include binaries with best-fit $P_{\rm orb}<1000$ days, which is
roughly the time baseline of observations used to calculate DR3 solutions. The
figure shows that typical astrometric uncertainties – especially the
uncertainty in $a_{0}$ – depend both on apparent magnitude and on $a_{0}$. At
fixed apparent magnitude, binaries with $a_{0}\approx 3\,\rm mas$ have
uncertainties in $a_{0}$ and $\varpi$ that are respectively about
$2.5\,\times$ and $1.7\,\times$ larger than sources with $a_{0}\approx
0.3\,\rm mas$.

It is not immediately obvious whether this trend reflects an actual dependence
of the astrometric uncertainties on $a_{0}$ among all solutions that were
initially fit, or is mainly a result of the astrometric SNR cuts imposed on
the published solutions (Equations 10-12). Figure 15 shows that removal of
solutions with larger uncertainties and short periods likely has a significant
effect on trends with $a_{0}$: there is a “wedge” of missing solutions with
small $a_{0}$ and larger uncertainties.

To determine whether an inherent dependence of the astrometric uncertainties
on $a_{0}$ (or another parameter that scales with it) is required to explain
the trends in Figure 14, we preformed the following experiment. Using the
simulation setup described in Appendix E, we begin with all sources in the
gaia_source catalog with parallax_over_error > 10 and create a simulated
binary population that we then subject to the DR3 SNR cuts. To predict the
uncertainties in $a_{0}$ and $\varpi$ for each simulated binary, we first
tabulate the mean and standard deviation of $\sigma_{a_{0}}$ and
$\sigma_{\varpi}$ as a function of apparent $G-$band magnitude, multiply this
mean by an $a_{0}$-dependent scaling factor, and then draw uncertainties for
each source from a Gaussian with the appropriate mean and standard deviation.
Finally, we compare the simulated distribution of $\sigma_{a_{0}}$ and
$\sigma_{\varpi}$ for binaries surviving the quality cuts to the observed
distribution shown in Figure 14.

We parameterize the scaling factor as $f=a\log(a_{0}/{\rm mas})+b$. We find
that when $a=0$ and $b=1$ (i.e., no dependence of uncertainties on $a_{0}$),
the predicted correlation between $a_{0}$ and the astrometric uncertainties is
weaker than observed. We can best reproduce the observed distribution for
$\sigma_{a_{0}}$ when $a=1.3$ and $b=1.4$; for $\sigma_{\varpi}$, we find
$a=0.64$ and $b=1.2$. We adopt these values when predicting the expected
astrometric uncertainties of a simulated binary population in Section 9.

Figure 14: Mean uncertainty in photocenter semi-major axis ($a_{0}$; left) and
parallax (right) for all Gaia DR3 astrometric solutions with $P_{\rm
orb}<1000\,\rm days$. Fainter sources and sources with large $a_{0}$ have
larger astrometric uncertainties. The correlation between $a_{0}$ and the
astrometric uncertainties appears to be only partially a result of period-
dependent quality cuts imposed on sources in DR3 (see text). We incorporate an
explicit $a_{0}$-dependence when predicting the astrometric uncertainties of a
hypothetical binary observed by Gaia.  Figure 15: Uncertainties in parallax
(top) and photocenter semi-major axis ($a_{0}$; bottom) for all sources with
astrometric solutions published in DR3. Left column shows uncertainties as a
function of $a_{0}$; right column shows them as a function of $P_{\rm orb}$.
Red star shows Gaia BH1. A correlation is apparent between the astrometric
uncertainties and $a_{0}$. This is in part because sources with small $a_{0}$
and large uncertainties were removed from the catalog through period-dependent
quality cuts (Equations 10-12). More structure is apparent as a function of
$P_{\rm orb}$: binaries with periods close to 0.5 years and 1 year, as well as
harmonics of the 63 day precession period, have larger uncertainties.
Uncertainties in $a_{0}$ also increase at $P_{\rm orb}>700$ days, likely due
to incomplete phase coverage.


CERN-TH-2021-230

UUITP-66/21

#  Evaluation of multi-loop multi-scale Feynman integrals for precision
physics

Ievgen Dubovyk  Institute of Physics, University of Silesia, Katowice, Poland
Ayres Freitas  Pittsburgh Particle physics, Astrophysics & Cosmology Center
(PITT PACC),
Department of Physics & Astronomy, University of Pittsburgh, Pittsburgh, PA
15260, USA  Janusz Gluza  Institute of Physics, University of Silesia,
Katowice, Poland  Krzysztof Grzanka  Institute of Physics, University of
Silesia, Katowice, Poland  Martijn Hidding  Department of Physics and
Astronomy, Uppsala University, SE-75120 Uppsala, Sweden  Johann Usovitsch
Theoretical Physics Department, CERN, 1211 Geneva, Switzerland

###### Abstract

Modern particle physics is increasingly becoming a precision science that
relies on advanced theoretical predictions for the analysis and interpretation
of experimental results. The planned physics program at the LHC and future
colliders will require three-loop electroweak and mixed electroweak-QCD
corrections to single-particle production and decay processes and two-loop
electroweak corrections to pair production processes. This article presents a
new semi-numerical approach to multi-loop multi-scale Feynman integrals
calculations which will be able to fill the gap between rigid experimental
demands and theory. The approach is based on differential equations with
boundary terms specified at Euclidean kinematic points. These Euclidean
boundary terms can be computed numerically with high accuracy using sector
decomposition or other numerical methods. They are then mapped to the physical
kinematic configuration with a series solution of the differential equation
system. The method is able to deliver 8 or more digits precision, and it has a
built-in mechanism for checking the accuracy of the obtained results. Its
efficacy is illustrated with state-of-the-art examples for three-loop self-
energy and vertex integrals and two-loop box integrals.

(Dedicated to the memory of Tord Riemann)

###### pacs:

14.40.Be, 13.40.Gp, 13.66.Bc, 13.40.Em

#### I. Introduction.

With the discovery of the Higgs boson at the Large Hadron Collider (LHC), all
building blocks of the Standard Model (SM) have been experimentally confirmed,
with the only exception of the Higgs self-coupling, which still awaits direct
measurement. However, the SM does not account for important phenomena such as
dark matter and the matter-antimatter asymmetry, so that physics beyond the SM
is needed. It is reasonable to expect that this new physics couples to the
electroweak and/or Higgs sector of the SM, since there are important model-
building constraints for couplings to the strong force 111A new physics
particle with renormalizable QCD couplings but no weak or Higgs interactions
would be stable and thus cosmologically excluded..

Therefore, possible evidence for such new physics can be explored in precision
studies of electroweak and Higgs physics at the high-luminosity run of the LHC
(HL-LHC) or one of several proposed future $e^{+}e^{-}$ colliders: FCC-ee
Abada:2019zxq , CEPC CEPCStudyGroup:2018ghi , ILC Baer:2013cma ;
Bambade:2019fyw , CLIC Linssen:2012hp ; Charles:2018vfv . Through their high
integrated luminosities of several ab-1, these machines will be sensitive to
very small deviations between the measured value and the SM expectation for a
given observable. Thus they can probe extremely feebly coupled new particles
or very large new physics scales of tens of TeV.

The SM predictions for these precision analyses are obtained by computing
higher-order quantum corrections. At the HL-LHC, some of the most interesting
precision studies are Higgs boson production and lepton pair (Drell-Yan)
production. For the former, one of the largest sources of theoretical
uncertainty stems from mixed QCD-electroweak corrections Anastasiou:2008tj ;
Anastasiou:2016cez . While some partial results at this order have been
computed Bonetti:2017ovy ; Anastasiou:2018adr ; Bonetti:2020hqh ;
Becchetti:2020wof ; Becchetti:2021axs , contributions from electroweak
diagrams with internal top quarks, both for 3-loop Higgs production and 2-loop
Higgs+jet production, are still needed to complete this missing piece. For
Drell-Yan production, 2-loop electroweak corrections for the full process
$pp\to\ell^{+}\ell^{-}$, not just on the Z-boson resonance, are important
since LHC measurements cover a broad range of invariant mass ATLAS:2018gqq ;
CMS:2018ktx .

Similarly, electroweak 2-loop corrections for several different pair
production processes will be essential for the physics goals of future
$e^{+}e^{-}$ colliders Freitas:2019bre : $e^{+}e^{-}\to W^{+}W^{-}$,
$e^{+}e^{-}\to ZH$, and $e^{+}e^{-}\to f\bar{f}$. Measurements of these cross-
sections will allow us to determine the W-boson mass with high precision,
constrain anomalous couplings between gauge bosons and/or the Higgs boson, and
probe heavy neutral vector bosons ($Z^{\prime}$ bosons). Currently, some
results for mixed QCD-electroweak 2-loop corrections are available
Gong:2016jys ; Sun:2016bel ; Heller:2020owb ; Bonciani:2021zzf , but so far no
complete electroweak 2-loop calculation for any pair production process has
been carried out. Even higher-order corrections will be needed for studies of
Z-boson production and decay at these future $e^{+}e^{-}$ colliders, as well
as the indirect prediction of the $W$ boson mass from the Fermi constant. To
match the expected experimental precision, 3-loop and partial 4-loop self-
energy and vertex corrections will be required Blondel:2018mad ;
Freitas:2019bre , which is one order of perturbation theory beyond the current
state of the art Dubovyk:2019szj .

It should be emphasized that these are loop corrections in the full SM,
involving many massive particles inside the loops. The currently most advanced
techniques for analytically computing such multi-loop Feynman integrals first
reduce them to a small set of master integrals, which then are solved by
constructing suitable differential equations (DEs), see Ref. Kotikov:2021tai
for a recent review. Both of these steps require integration-by-parts (IBP)
equation systems Chetyrkin:1981qh ; Laporta:2000dsw  that become intractable
for multi-loop integrals with many masses. Instead, one must resort to
numerical integration techniques.

The recent calculation of full 2-loop corrections to Z-boson production and
decay Dubovyk:2016aqv ; Dubovyk:2018rlg ; Dubovyk:2019szj  made use of
numerical evaluations based on sector decomposition (SD) Binoth:2000ps ;
Borowka:2017idc ; Borowka:2018goh ; Smirnov:2008py ; Smirnov:2015mct  and
Mellin-Barnes (MB) representations Gluza:2007rt ; Czakon:2005rk ; Usovitsch-
phd:2018qmt ; Czakon:2005rk ; Dubovyk:2016ocz ; Usovitsch:2018shx ;
Dubovyk:2016aqv . However, these methods require large amounts of computing
resources and do not always converge to the required level of accuracy, so
that a straightforward extension to more loops and/or legs is not possible.
Due to numerical cancellations between individual loop integrals, at least 8
digits precision are required in many cases for practical applications
Dubovyk:2018rlg . Another interesting approach, which allows to tackle a wider
class of problems, evaluates a set of master integrals by numerically solving
a DE system, either in terms of kinematic parameters Mandal:2018cdj ;
Czakon:2020vql  or in terms of an auxiliary flow variable Liu:2017jxz ;
Liu:2018dmc ; Bronnum-Hansen:2020mzk .

This article introduces an efficient but still very general approach that can
be applied to many challenging 2- and 3-loop problems with multiple mass and
momentum scales 222In principle, it is also applicable at higher loop orders,
but we have not studied any examples beyond three loops.. The key elements are
a system of DEs, with boundary terms evaluated at one or more Euclidean
(space-like) kinematic points (which can be reliably determined to high
precision with numerical methods). The DEs are then solved, using series
expansions, to obtain the final result at the physical Minkowski (time-like)
kinematic point. This approach,which is already fully automated in its main
parts, will be described in more detail in the next section. In section III we
will apply this technique to examples of SM self-energy and vertex Feynman
integrals that occur in three-loop Z-decay corrections. The chosen examples
are very difficult to evaluate with other analytical or numerical methods. A
summary and outlook are given in the final section. The supplemental material
included with this submission contains input parameters for the example
integrals, additional examples for the application of our calculation
technique, e.g. a 2-loop four-scale box diagram, and miscellaneous
implementational details and remarks for a more detailed discussion of the
method.

#### II. Description of the method.

Solving Feynman integrals from DEs is an approach initiated in the last decade
of the last century Kotikov:1990kg ; Kotikov:1991hm ; Kotikov:1991pm ;
Remiddi:1997ny . Many families of Feynman integrals admit a choice of master
integrals for which the system of DEs has a particularly simple “canonical”
form Henn:2013pwa , which in many cases can be straightforwardly solved in
terms of multiple polylogarithms.

More generally, not all Feynman integrals are of polylogarithmic type, and it
can become increasingly difficult to find a closed set of analytic functions
in terms of which the DEs can be solved. In such cases, it is often still
possible to find precise numerical solutions. For example, the DE system can
be integrated numerically Mandal:2018cdj ; Czakon:2020vql ; Liu:2017jxz ;
Liu:2018dmc ; Bronnum-Hansen:2020mzk . In this work, we use the approach of
iterated series expansions, which is both fully automatable and numerically
efficient.

Let us give a brief overview of the method. Consider a basis of master
integrals (MIs), $\vec{F}(x,\epsilon)$, depending on a single scale $x$. We
work in dimensional regularization, with $D=4-2\epsilon$ space-time
dimensions. We may then derive DEs of the form:

$\frac{d}{dx}\vec{F}(x,\epsilon)=\hat{M}(x,\epsilon)\vec{F}(x,\epsilon)\,,$
(1)

where $\hat{M}(x,\epsilon)$ is a block-triangular matrix. Each block is
associated with a sector of integrals. If we denote such a sector by
$\vec{f}_{i}(x,\epsilon)$, we can decompose the DEs in the form:

$\frac{d}{dx}\vec{f}_{i}(x,\epsilon)=M_{i}(x,\epsilon)\vec{f}_{i}(x,\epsilon)+B_{i}(x,\epsilon)\vec{g}_{i}(x,\epsilon)\,,$
(2)

where $M_{i}(x,\epsilon)$ denotes the diagonal block of $\hat{M}(x,\epsilon)$
corresponding to the sector $i$, and
$B_{i}(x,\epsilon)\vec{g}_{i}(x,\epsilon)$ captures the off-diagonal terms.
One can then expand the integrals and matrices in $\epsilon$:

$\displaystyle\vec{f}_{i}(x,\epsilon)$
$\displaystyle=\sum_{j=-k}^{\infty}\vec{f}_{i}^{(j)}(x,\epsilon)\,\epsilon^{j}\,,$
$\displaystyle M_{i}(x,\epsilon)$
$\displaystyle=\sum_{j=0}^{\infty}M_{i}^{(j)}(x,\epsilon)\,\epsilon^{j}\,,$
(3)

and solve the system order by order in $\epsilon$. For a given basis, the
condition that $M_{i}(x,\epsilon)$ is finite in $\epsilon$ is not always
manifest. We obtain such a form by rescaling individual master integrals with
powers of $\epsilon$, until the matrix is finite. For further ideas and
software to help facilitate the choice of MIs, see Refs. Prausa:2017ltv ;
Gituliar:2017vzm ; vonManteuffel:2017hms ; Adams:2017tga ; Henn:2020lye ;
Dlapa:2020cwj ; Usovitsch:2020jrk ; Smirnov:2020quc .

The DE system in eq. (2) fixes the master integrals up to some boundary
conditions. It turns out that, in the case of our automated DE approach, a
convenient choice for the boundary terms are MIs which are finite in the
dimensional regulator $\epsilon$. We use the package Reduze Panzer:2015ida ;
von_Manteuffel_2015 ; von_Manteuffel_2016 ; vonManteuffel:2012np  to identify
these MIs. They can be evaluated efficiently for Euclidean kinematics using
the method of SD, since only a small number of sectors is needed for finite
integrals and no contour deformation is required to avoid Minkowskian
thresholds. We employ the package pySecDec Borowka:2017idc ; Borowka:2018goh
for this purpose. The derivation of a DE system is done with the help of the
IBP reduction program Kira Maierhoefer:2017hyi ; Klappert:2020nbg ;
Klappert:2019emp ; Klappert:2020aqs . With the boundary terms fixed
numerically and the DE system derived analytically, we transport the Euclidean
point to the Minkowski point with the aid of the method of series expansions
of the DE system Pozzorini:2005ff ; Moriello:2019yhu ; Bonciani:2019jyb 333
See also Ref. Fael:2021kyg  for a similar approach for single-scale problems.
as implemented in DiffExp Hidding:2020ytt .

As demonstrated in Fig. 1 we may choose different Euclidean points to fix the
boundary terms numerically. This allows us to obtain a numerical error
estimate of our automated method by taking the difference of two generated
results for the same final Minkowski point. A more detailed discussion of the
error estimate is provided in the supplemental material.

Typically, the transport from the Euclidean boundary point to the physical
Minkowski kinematics requires several steps since the convergence radius of
the series expansion at the boundary point is not large enough to reach the
target point. The program DiffExp automatically determines the convergence
radius and the number of required transport steps.

In general, the complexity of the multi-loop computation increases with the
number of loops and independent scales and the number of MIs involved. In our
automated approach, the largest investment of computing resources is required
for the IBP reduction with Kira and the numerical evaluation of the boundary
terms with pySecDec. However, the former needs to be done only once for a
given Feynman integral family, and the latter only once for a given choice of
mass parameter values. The transport to the Minkowski region with DiffExp is
very fast, so that one can easily evaluate results for multiple different
kinematic points, as needed e.g. for phase-space integrations. Quantitative
information on the run time for our approach is given in the supplemental
material. There the reader can also find a description for how our method can
be extended to problems with multiple time-like momentum scales.

Figure 1: Illustration of the DE transport method. The boundary conditions for
the integral $f_{i}$ are evaluated at one or several Euclidean points $\rm
A_{k}$, where the integral is purely real and one can obtain robustly
converging numerical results with the SD method (using the package pySecDec in
our case). The boundary value(s) are transported to the physical kinematic
point, using solutions of the DE system eq. (2) derived with DiffExp, yielding
the final result indicated by the red dot in the figure. The numerical
uncertainty of a boundary value translates to an error estimate of the final
result, as illustrated by the error bars in a zoomed-in area in the dotted
circle, which permits a non-trivial cross-check if several boundary values
$\rm A_{k}$ are employed.

#### III. Results and discussion.

To demonstrate the power and broad applicability of our method, in the
following and in the supplemental material, we present examples for 3-loop
self-energy and vertex integrals and 2-loop box integrals. As discussed in the
introduction, these are all examples of key theory ingredients for the physics
program of future $e^{+}e^{-}$ colliders and/or the HL-LHC. The 3-loop
integrals are needed for currently unknown third-order corrections to
electroweak precision observables connected with Z-boson producton and decay,
whereas two-loop box integrals are important to improve the precision of
several $2\to 2$ processes, such as $W^{+}W^{-}$, $ZH$ or $f\bar{f}$
production 444See Refs. Song:2021vru ; Liu:2021wks  for recent independent
efforts on electroweak two-loop box integrals..

The technique described in this article allows one to compute the desired
integrals to, in principle, arbitrary order in the dimension regularization
parameter $\epsilon=(4-D)/2$ with multi-digit precision. To achieve a certain
order $\epsilon^{k}$, some boundary terms need to be evaluated to higher
orders $k^{\prime}>k$ in $\epsilon$. The required order $k^{\prime}$ is
determined automatically from the IBP relations. For the examples shown below,
some simple boundary-term integrals have to be computed to ${\cal
O}(\epsilon^{7})$, whereas no more than ${\cal O}(\epsilon^{3})$ is needed for
more complicated boundary terms. When evaluating the boundary terms with SD as
implemented in pySecDec, the computing time grows approximately linear with
the order in $\epsilon$.

All of the following numerical examples are based on the input parameters
given in the supplemental material.

|
---|---
lhNp1 |  vtwPl
Figure 2: Three loop self-energy non-planar and planar vertex diagrams which
correspond to integrals in (4) and (8), respectively. W, Z and t stand for the
W-boson, Z-boson and top quark, respectively.

Example 1. As part of the 3-loop ${\cal O}(\alpha^{2}\alpha_{\rm s})$
corrections to electroweak precision observables, one encounters the following
scalar non-planar self-energy integral with eight propagators and only one
massive W- or Z-boson internal line 555In this and the following examples, the
large number of massless propagators occurs because all SM fermions except the
top quark are taken to be massless. (see Fig. 2 left):

$\displaystyle
I_{\text{lhNp1}}[D,\\{a_{i}\\},p^{2},M_{a}^{2}]=\int\frac{\mathfrak{D}q_{1}\mathfrak{D}q_{2}\mathfrak{D}q_{3}}{[(q_{1}-q_{2})^{2}]^{a_{1}}[q_{2}^{2}]^{a_{2}}}$
$\displaystyle\times\frac{1}{[(q_{1}-q_{3})^{2}]^{a_{3}}[(q_{2}-q_{3})^{2}-M_{a}^{2}]^{a_{4}}[q_{3}^{2}]^{a_{5}}}$
$\displaystyle\times\frac{[q_{1}^{2}]^{-a_{9}}}{[(q_{1}+p)^{2}]^{a_{6}}[(q_{1}-q_{2}+p)^{2}]^{a_{7}}[(q_{3}+p)^{2}]^{a_{8}}}\,.$
(4)

where $\mathfrak{D}q_{n}\equiv\frac{d^{D}q_{n}}{i\pi^{D/2}}$ and $a=\rm W,Z$.
This example, for the parameter point $p^{2}=M_{\rm Z}^{2}$ and $M_{a}=M_{\rm
Z}$ belongs to a group of integrals which are difficult to evaluate with SD
due to threshold effects. Using pySecDec with $10^{7}$ integration points we
obtain a result with less than two digits precision:

$\displaystyle I_{\text{lhNp1}}^{{\tt
pySecDec}}[4-2\epsilon,1,1,1,1,1,1,1,1,0,M_{\rm Z}^{2},M_{\rm Z}^{2}]$
$\displaystyle=0.460-19.164\;i\pm(0.298+0.281\;i).$ (5)

Increasing the number of integration points does not improve the accuracy
substantially. On the other hand, pySecDec can deliver accurate results for
Euclidean parameter points, $p^{2}<0$, which are used as boundary terms for
our automated DE transport. We thus obtain stable and precise results at the
physical point:

$\displaystyle I_{\text{lhNp1}}[4-2\epsilon,1,1,1,1,1,1,1,1,0,M_{\rm
Z}^{2},M_{\rm Z}^{2}]$ $\displaystyle=-0.000000000-19.1262302\;i$
$\displaystyle\phantom{=}+(151.51529-150.40641\;i)\;\epsilon+{\cal
O}(\epsilon^{2}),$ (6) $\displaystyle
I_{\text{lhNp1}}[4-2\epsilon,1,1,1,1,1,1,1,1,0,M_{\rm Z}^{2},M_{\rm W}^{2}]$
$\displaystyle=(5.1112260-18.5692007\;i)$
$\displaystyle\phantom{=}+(194.660753-78.842016\;i)\;\epsilon+{\cal
O}(\epsilon^{2}).$ (7)

Here and in all the following results, we show all significant digits, i.e.
the numerical error only affects digits beyond the ones shown in the
equations. The error estimation will be described in more detail in the
supplemental material. The integral family $I_{\text{lhNp1}}$ (4) involves 30
master integrals and is considered simple in the context of our method.

Example 2. The next example is a family of 3-loop vertex integrals with one
massive top quark and two massive W-boson propagators, see Fig. 2 (right),
defined as

$\displaystyle
I_{\text{vtwPl}}[D,\\{a_{i}\\},p^{2},M_{\text{W}}^{2},m_{\text{t}}^{2}]=\int\frac{\mathfrak{D}q_{1}\mathfrak{D}q_{2}\mathfrak{D}q_{3}}{[q_{3}^{2}-M_{\rm
W}^{2}]^{a_{1}}[q_{2}^{2}]^{a_{2}}}$
$\displaystyle\times\frac{1}{[q_{1}^{2}]^{a_{3}}[(q_{1}-p)^{2}]^{a_{4}}[(q_{2}-p)^{2}]^{a_{5}}[(q_{3}-p)^{2}-M_{\text{W}}^{2}]^{a_{6}}}$
$\displaystyle\times\frac{[(q_{1}-q_{3})^{2}]^{-a_{10}}[(q_{1}-p_{2})^{2}]^{-a_{11}}[(q_{2}-p_{2})^{2}]^{-a_{12}}}{[(q_{3}-p_{1})^{2}]^{a_{7}}[(q_{2}-q_{3})^{2}-m_{\text{t}}^{2}]^{a_{8}}[(q_{1}-q_{2})^{2}]^{a_{9}}},$
(8)

where $p=p_{1}+p_{2}$ and $p_{1}^{2}=p_{2}^{2}=0$. These integrals also appear
in so far unknown ${\cal O}(\alpha^{2}\alpha_{\rm s})$ corrections to Z-pole
electroweak precision observables, constituting their most difficult parts.

With pySecDec we are unable to obtain a numerical result for the Minkowski
point $p^{2}=M_{\rm Z}^{2}$. The problem already starts with the contour
deformation which is necessary for SD with Minkowski kinematics and which
fails to complete in a reasonable time. Similar to the SD method, the MB
technique fails to deliver high-accuracy results for the considered integrals
for $p^{2}=M_{\rm Z}^{2}$.

Using our automated DE transport method, the calculation requires the
numerical evaluation of 77 master integrals with Euclidean kinematics,
$p^{2}<0$, for the boundary terms. For the purpose of the present example,
they have been evaluated with pySecDec to 10-digit accuracy. After the
transport to the physical point $p^{2}=M_{\rm Z}^{2}$, we get at least eight
significant digits for integrals of the family (8) up to tensor rank-3 (i.e.
$-3\leq a_{10}+a_{11}+a_{12}\leq 0$). We here give numerical result for one
rank-3 case:

$\displaystyle I_{\text{vtwPl}}[1,1,1,1,1,1,1,1,1,-1,-1,-1,M_{\rm
Z}^{2},M_{\rm W}^{2},M_{\rm t}^{2}]=$ $\displaystyle
0.0833333333/\epsilon^{3}+0.636273147/\epsilon^{2}$
$\displaystyle+(0.63462699+0.77044487\;i)/\epsilon$
$\displaystyle+(5.5847828+6.1606031\;i)+{\cal O}(\epsilon),$ (9)

Additional examples, a 3-loop self-energy diagram with many massive
propagators and a two-loop box diagram with four scales, are discussed in the
supplemental material.

#### IV. Summary and Outlook.

In this work, we have proposed an efficient and versatile approach for the
evaluation of a wide class of massive multi-loop, multi-scale Feynman
integrals numerically, with typically 8 or more digits precision. It is based
on the method of DEs with boundary terms specified for Euclidean kinematics,
which are transported to the physical Minkowski kinematics using series
solutions of the DE. The Euclidean boundary term integrals avoid all threshold
singularities and thus can be straightforwardly evaluated numerically. Our
implementation combines the public programs Kira, Reduze, pySecDec and DiffExp
in a way that allows us to automatically construct the required integral
families and the transport from the Euclidean boundary point to the physical
kinematic point.

In principle, the technique can be extended to higher numerical accuracy and
to wider classes of integrals with more loops and more external legs. A major
bottleneck are the integration-by-parts (IBP) reductions that are needed to
construct the DE system. A significant speed-up of this step is achieved when
using numerical values for the relevant mass and kinematic parameters. In
addition, the evaluation of the boundary terms for Euclidean kinematics can be
time-consuming if a high level of precision is required. Fortunately, there
are ongoing improvements to the sector decomposition (SD) and Mellin-Barnes
(MB) methods, see e.g. Refs. Borowka:2018goh ; Borinsky:2020rqs  and
Ananthanarayan:2020fhl .

It is worth mentioning that the 3-loop examples presented in this article are
very difficult to solve with existing analytical techniques (e.g. using IBP
and DE) and general numerical methods (such as SD or MB methods). The proposed
new technique is sufficiently general to provide the foundation for the
computation of the required 3-loop corrections needed for electroweak and
Higgs precision studies at the HL-LHC and future $e^{+}e^{-}$ colliders, which
are key elements of the physics program of these machines Blondel:2018mad ;
Anastasiou:2016cez . Other applications include flavor physics at Belle-II and
low-energy precision tests of the Standard Model.

#### Acknowledgements

This work has been supported in part by the Polish National Science Center
(NCN) under grant 2017/25/B/ST2/01987, the Research Excellence Initiative of
the University of Silesia in Katowice, the U.S. National Science Foundation
under grant nos. PHY-1820760 and PHY-2112829, MH is supported by the European
Research Council under ERC-STG-804286 UNISCAMP.

##  I Supplemental Material

###  I.1 Numerical results for relevant Master Integrals

A set of auxiliary files are provided together with this article, which
contain numerical results for the master integrals evaluated at the physical
point (i.e. for Minkowski kinematics). Here are the building blocks for
computing the results in eqs. (6),(7),(9) in the main text and eqs. (11)-(15).
The following files correspond to the different examples described in the
text:

Example 1 (main text): | master_lhNp1_minkowski.m
---|---
Example 2 (main text): | master_vtwPl_minkowski.m
Example 3 (main text): | master_box2l_minkowski.m
Example 4 (appendix): | master_taNp1_minkowski.m

The objects lhNp1D4[...], lhNp1D6[...], lhNp1D8[...], lhNp1D10[...],
taNp1D4[...], taNp1D6[...], vtwPlD4[...], vtwPlD6[...], vtwPlD8[...] and
box2lD6[...] define the master integrals for the four topologies lhNp1, taNp1,
vtwPl and box2l, where D4, D6, D8, D10 denote that the master integral is
evaluated around $\epsilon=0$ for $d=\\{4,6,8,10\\}-2\epsilon$, respectively.

The notation which we use comes from chosen classification of the 3-loop self-
energy (SE) and vertex SM diagrams, in general we have

    ta = SE diagrams with top quark and photon
    lh = SE diagrams with light quark(s) W/Z/H
    vtw = 3-loop vertex topologies with top and W
    Pl = planar ladder topologies
    Np = non-planar topologies

###  I.2 Input parameters

The input parameters for the lhNp1, taNp1 and vtwPl topologies are summarized
in Tab. 1, where each parameter $M_{\rm Z}=m_{\rm Z}/m_{\rm Z}$, $M_{\rm
W}=m_{\rm W}/m_{\rm Z}$ and $M_{\rm t}=m_{\rm t}/m_{\rm Z}$ represents a
physical scale rescaled by $m_{\rm Z}$. They are written as rational numbers
suitable for IBP reductions with Kira. The dimensionful values for $m_{\rm
Z}$, $m_{\rm W}$, $m_{\rm t}$ are the same as in Ref. Dubovyk:2019szj , i.e.
$m_{\rm Z}=91.1876$ GeV, $m_{\rm W}=80.385$ GeV, $m_{\rm t}=173.2$ GeV.

Table 1: Input values, which are used for the 3-loop self-energy and vertex examples in this paper, are rescaled such that the Z-boson mass is 1. Mass | Input value
---|---
$M_{\rm Z}$ | $1$
$M_{\rm W}$ | $401925/455938=80.385/91.1876$
$M_{\rm t}$ | $433000/227969=173.2/91.1876$

###  I.3 3-loop self-energy integrals with top quarks

taNp1

Figure 3: Three loop self-energy non-planar integral defined in (10). Z and t
stands for the massive SM Z gauge boson and the top quark, respectively.

$\displaystyle
I_{\text{taNp1}}[D,\\{a_{i}\\},p^{2},m_{\text{t}}^{2}]=\int\mathfrak{D}q_{1}\mathfrak{D}q_{2}\mathfrak{D}q_{3}$
$\displaystyle\frac{1}{[(q_{1}+p)^{2}]^{a_{1}}[(q_{1}-q_{2})^{2}-M_{\rm
t}^{2}]^{a_{2}}[(q_{1}-q_{2}+p)^{2}-M_{\rm t}^{2}]^{a_{3}}}$
$\displaystyle\times\frac{1}{[q_{2}^{2}-M_{\rm
t}^{2}]^{a_{4}}[(q_{1}-q_{3})^{2}-M_{\rm
t}^{2}]^{a_{5}}[(q_{2}-q_{3})^{2}]^{a_{6}}}$
$\displaystyle\times\frac{[q_{1}^{2}]^{-a_{9}}[(q_{2}+p)^{2}]^{-a_{10}}}{[q_{3}^{2}-M_{\rm
t}^{2}]^{a_{7}}[(q_{3}+p)^{2}-M_{\rm t}^{2}]^{a_{8}}},$ (10)

see Fig. 3, include the top quark which circulates through six propagators.
These integrals are computed with at least ten digits accuracy in the
Minkowski point around $\epsilon=0$ for $D=4-2\epsilon$

$\displaystyle I_{\text{taNp1}}[1,1,1,1,1,1,1,1,-2,-1,M_{\rm Z}^{2},M_{\rm
t}^{2}]=$ $\displaystyle
8.27490485938/\epsilon^{3}-34.9869281045/\epsilon^{2}$
$\displaystyle+102.43077689/\epsilon-253.5072352,$ (11) $\displaystyle
I_{\text{taNp1}}[4-2\epsilon,1,1,1,1,1,1,1,1,-1,-2,M_{\rm Z}^{2},M_{\rm
t}^{2}]=$ $\displaystyle
9.47745432492/\epsilon^{3}-40.4955852564/\epsilon^{2}$ (12)
$\displaystyle+116.63419570/\epsilon-273.3763275,$ (13) $\displaystyle
I_{\text{taNp1}}[1,1,1,1,1,1,1,1,0,-3,M_{\rm Z}^{2},M_{\rm t}^{2}]=$
$\displaystyle 19.8715753165/\epsilon^{3}-74.436608700/\epsilon^{2}$ (14)
$\displaystyle+239.02713087/\epsilon-540.2221570.$ (15)

Note the tenth propagator in (10) is linearly dependent and can be written in
terms of the first nine propgators. We included it as an auxiliary propagator
to the definition to improve the readability of the final results in
(11)-(15).

###  I.4 Two-loop box diagram

Figure 4: Two-loop box diagram with four scales: $s,t,m_{1},m_{2}$.

We consider a family of box integrals that are part of the ${\cal
O}(\alpha^{2})$ corrections to massive $e\mu$ scattering Banerjee:2020tdt ,
see Fig. 4:

$\displaystyle
I_{\text{box2l}}[D,\\{a_{i}\\},s,t,m_{1}^{2},m_{2}^{2}]=\int\frac{\mathfrak{D}q_{1}\mathfrak{D}q_{2}}{[(q_{1}-p_{1})^{2}]^{a_{1}}[q_{1}^{2}-m_{1}^{2}]^{a_{2}}}$
$\displaystyle\times\frac{1}{[(q_{1}-p_{1}+p_{3})^{2}-m_{1}^{2}]^{a_{3}}[(q_{1}-q_{2})^{2}-m_{1}^{2}]^{a_{4}}}$
$\displaystyle\times\frac{[(q_{2}-p_{1})^{2}]^{-a_{8}}[(q_{1}-p_{2})^{2}]^{-a_{9}}}{[(q_{2}-p_{1}+p_{3})^{2}]^{a_{5}}[(q_{2}+p_{2})^{2}-m_{2}^{2}]^{a_{6}}[q_{2}^{2}]^{a_{7}}}$
(16)

where $p_{1}^{2}=p_{3}^{2}=m_{1}^{2}$, $p_{2}^{2}=p_{4}^{2}=m_{2}^{2}$,
$p_{1}\,p_{2}=(s-m_{1}^{2}-m_{2}^{2})/2$, $p_{1}\,p_{3}=(2m_{1}^{2}-t)/2$,
$p_{2}\,p_{3}=(s+t-m_{1}^{2}-m_{2}^{2})/2$. With input parameters $s=2$,
$t=5$, $m_{1}^{2}=4$ and $m_{2}^{2}=16$ we obtain

$\displaystyle I_{\text{box2l}}[2,1,1,1,1,1,1,0,0,s,t,m_{1}^{2},m_{2}^{2}]=$
$\displaystyle+0.000328707579/\epsilon^{2}$
$\displaystyle-(0.0014129475-0.0020653306\;i)/\epsilon$
$\displaystyle-(0.005702737-0.000485980\;i)+{\cal O}(\epsilon).$ (17)

The $I_{\text{box2l}}$ integral family involves 55 master integrals. The input
parameters for the box2l topology are $s=2$, $t=5$, $m_{1}^{2}=4$ and
$m_{2}^{2}=16$.

###  I.5 Tracing the run time of the calculations

For the calculations demonstrated in this paper it makes sense to categorize
the run time in three parts. The first part includes the preparation of the DE
system suitable for later use with DiffExp. For this step, the major
bottleneck is the run time of Kira. For the most complicated example vtwPl we
need 4 hours to construct the DE system on a 12 core, 2.7 GHz Intel Xenon
processor cluster with 128 GB of RAM. The second part is the computation of
Euclidean boundary terms, which we evaluate with the sector decomposition
program pySecDec. Since the boundary terms are finite integrals by
construction and are free of thresholds in the Euclidean region, the Monte-
Carlo (MC) or quasi-Monte-Carlo (QMC) scaling rules apply straightforwardly.
For QMC the theoretical bound on the error is proportional to
$\mathcal{O}(1/N)$, where $N$ is the number of points. For traditional MC, the
bound is $\mathcal{O}(1/\sqrt{N})$. For the most complicated example vtwPl, we
ran the QMC integrator configured with a maximum of $10^{7}$ points. We ran
the computation on a machine equipped with a 16 core Threadripper Pro 3955WX.
Two different points in the Euclidean region were chosen, and both took
approximately 3 days to complete. The maximum relative error reported by
pySecDec was approximately $2\cdot 10^{-9}$. We expect that the precision of
our result can be further improved by running on a more powerful machine, or
by allocating a longer running time. The final part of the computation is the
transport of the boundary terms to the Minkowski regions. The run time of this
part of the computation depends on the size of the differential equations, the
order of the homogeneous parts of the differential equations, and the number
of line segments. In our case we are interested in the evaluation of the
Feynman integral in one particular point. With DiffExp this is accomplished in
about 3.5 hours on a single CPU core for the most complicated example vtwPl.
Note that once the computation with DiffExp is finished the integral is
available for a fast evaluation (in terms of milliseconds) at an arbitrary
point along the line of the transport from Euclidean to Minkowski point.

###  I.6 Tracing the error of the calculations

The numerical error estimation has two main ingredients, the numerical series
expansion of the differential equation (DE) system with DiffExp and the
numerical evaluation of the boundary terms with pySecDec. The numerical error
from the series expansion can always be rendered negligible compared to the
error from the boundary terms by evaluating the expansion to sufficiently high
order. On the other hand, the numerical errors of the initial boundary terms
can usually be directly mapped to the final result. However, instead of
relying on the error estimate from pySecDec, we verify the accuracy by
carrying out separate transports from two different initial boundary points to
the same final Minkowski point and taking the difference as an error estimate.
This cross-check was sufficient for the examples presented here. Also, if
necessary, it is relatively straightforward to increase the accuracy of the
boundary terms by using more Monte-Carlo integration points, since the SD
integrand is well-behaved in the Euclidean region.

###  I.7 Iterative approach for multi-scale kinematics

In the case of problems with multiple time-like momentum scales, one may use
our method in an iterative way. For instance, we may derive the DE system just
in one scale $m_{1}$, while setting all other scales $\\{m_{i}\\}\setminus
m_{1}$ to numerical values, to simplify the IBP reduction. Then one can carry
out the transport with respect to $m_{1}$ from the numerical point $a_{1}$,
where the boundary conditions are supplied with all scales taken to be
Euclidean (space-like), to a new point of interest denoted as $a_{2}$, where
$m_{1}$ is time-like. Next we derive a new set of DEs with respect to the next
scale $m_{2}$, and again we set all other scales $\\{m_{i}\\}\setminus m_{2}$
to numerical values, and transport from the point $a_{2}$ to a new point of
interest $a_{3}$. One can continue in this way until all ($n$) relevant
momentum scales have been transported from the Euclidean to the time-like
domain. The downside of the iterative approach is that one has to repeat the
IBP reduction whenever we want to reach a new kinematic point that differs in
any of the $n-1$ first transported variables. Thus new improvements in the IBP
reduction programs are highly anticipated.

## References

  * (1) A new physics particle with renormalizable QCD couplings but no weak or Higgs interactions would be stable and thus cosmologically excluded.
  * (2) A. Abada, et al., FCC-ee: The Lepton Collider: Future Circular Collider Conceptual Design Report Volume 2, Eur. Phys. J. ST 228 (2) (2019) 261–623.  doi:10.1140/epjst/e2019-900045-4.
  * (3) CEPC Study Group, J. Guimarães da Costa, Y. Gao, S. Jin, J. Qian, C. Tully, C. Young, L. Wang, M. Ruan, H. Zhu, Q. Ouyang et al., CEPC Conceptual Design Report: Volume 2 - Physics & Detector. http://cds.cern.ch/record/2678417. arXiv:1811.10545.
  * (4) H. Baer, T. Barklow, K. Fujii, Y. Gao, A. Hoang, S. Kanemura, J. List, H. E. Logan, A. Nomerotski, M. Perelstein, et al., The International Linear Collider Technical Design Report - Volume 2: Physics, arXiv:1306.6352.
  * (5) P. Bambade, et al., The International Linear Collider: A Global Project, arXiv:1903.01629.
  * (6) Physics and Detectors at CLIC: CLIC Conceptual Design ReportarXiv:1202.5940, doi:10.5170/CERN-2012-003.
  * (7) T. Charles, et al., The Compact Linear Collider (CLIC) - 2018 Summary Report 2/2018.  arXiv:1812.06018, doi:10.23731/CYRM-2018-002.
  * (8) C. Anastasiou, R. Boughezal, F. Petriello, Mixed QCD-electroweak corrections to Higgs boson production in gluon fusion, JHEP 04 (2009) 003.  arXiv:0811.3458, doi:10.1088/1126-6708/2009/04/003.
  * (9) C. Anastasiou, C. Duhr, F. Dulat, E. Furlan, T. Gehrmann, F. Herzog, A. Lazopoulos, B. Mistlberger, High precision determination of the gluon fusion Higgs boson cross-section at the LHC, JHEP 05 (2016) 058.  arXiv:1602.00695, doi:10.1007/JHEP05(2016)058.
  * (10) M. Bonetti, K. Melnikov, L. Tancredi, Three-loop mixed QCD-electroweak corrections to Higgs boson gluon fusion, Phys. Rev. D 97 (3) (2018) 034004.  arXiv:1711.11113, doi:10.1103/PhysRevD.97.034004.
  * (11) C. Anastasiou, V. del Duca, E. Furlan, B. Mistlberger, F. Moriello, A. Schweitzer, C. Specchia, Mixed QCD-electroweak corrections to Higgs production via gluon fusion in the small mass approximation, , JHEP 03 (2019) 162.  arXiv:1811.11211, doi:10.1007/JHEP03(2019)162.
  * (12) M. Bonetti, E. Panzer, V. A. Smirnov, L. Tancredi, Two-loop mixed QCD-EW corrections to $gg\to Hg$, JHEP 11 (2020) 045.  arXiv:2007.09813, doi:10.1007/JHEP11(2020)045.
  * (13) M. Becchetti, R. Bonciani, V. Del Duca, V. Hirschi, F. Moriello, A. Schweitzer, Next-to-leading order corrections to light-quark mixed QCD-EW contributions to Higgs boson production, Phys. Rev. D 103 (5) (2021) 054037.  arXiv:2010.09451, doi:10.1103/PhysRevD.103.054037.
  * (14) M. Becchetti, F. Moriello, A. Schweitzer, Two-loop amplitude for mixed QCD-EW corrections to $gg\to Hg$, arXiv:2112.07578.
  * (15) Measurement of the effective leptonic weak mixing angle using electron and muon pairs from $Z$-boson decay in the ATLAS experiment at $\sqrt{s}=8$ TeV, ATLAS-CONF-2018-037.
  * (16) A. M. Sirunyan, et al., Measurement of the weak mixing angle using the forward-backward asymmetry of Drell-Yan events in pp collisions at 8 TeV, Eur. Phys. J. C 78 (9) (2018) 701.  arXiv:1806.00863, doi:10.1140/epjc/s10052-018-6148-7.
  * (17) A. Freitas, et al., Theoretical uncertainties for electroweak and Higgs-boson precision measurements at FCC-ee, arXiv:1906.05379.
  * (18) Y. Gong, Z. Li, X. Xu, L. L. Yang, X. Zhao, Mixed QCD-EW corrections for Higgs boson production at $e^{+}e^{-}$ colliders, Phys. Rev. D 95 (9) (2017) 093003.  arXiv:1609.03955, doi:10.1103/PhysRevD.95.093003.
  * (19) Q.-F. Sun, F. Feng, Y. Jia, W.-L. Sang, Mixed electroweak-QCD corrections to e+e-$\rightarrow$HZ at Higgs factories, Phys. Rev. D 96 (5) (2017) 051301.  arXiv:1609.03995, doi:10.1103/PhysRevD.96.051301.
  * (20) M. Heller, A. von Manteuffel, R. M. Schabinger, H. Spiesberger, Mixed EW-QCD two-loop amplitudes for $q\bar{q}\to\ell^{+}\ell^{-}$ and $\gamma_{5}$ scheme independence of multi-loop corrections, JHEP 05 (2021) 213.  arXiv:2012.05918, doi:10.1007/JHEP05(2021)213.
  * (21) R. Bonciani, L. Buonocore, M. Grazzini, S. Kallweit, N. Rana, F. Tramontano, A. Vicini, Mixed strong$-$electroweak corrections to the Drell$-$Yan process, arXiv:2106.11953.
  * (22) A. Blondel, et al., Standard model theory for the FCC-ee Tera-Z stage, in: Mini Workshop on Precision EW and QCD Calculations for the FCC Studies : Methods and Techniques, Vol. 3/2019 of CERN Yellow Reports: Monographs, CERN, Geneva, 2018.  arXiv:1809.01830, doi:10.23731/CYRM-2019-003.
  * (23) I. Dubovyk, A. Freitas, J. Gluza, T. Riemann, J. Usovitsch, Electroweak pseudo-observables and Z-boson form factors at two-loop accuracy, JHEP 08 (2019) 113.  arXiv:1906.08815, doi:10.1007/JHEP08(2019)113.
  * (24) A. V. Kotikov, Differential equations and Feynman integrals, in: Antidifferentiation and the Calculation of Feynman Amplitudes, 2021.  arXiv:2102.07424.
  * (25) K. Chetyrkin, F. Tkachov, Integration by parts: The algorithm to calculate $\beta$ functions in 4 loops, Nucl. Phys. B192 (1981) 159–204.  doi:10.1016/0550-3213(81)90199-1.
  * (26) S. Laporta, High precision calculation of multiloop Feynman integrals by difference equations, Int. J. Mod. Phys. A 15 (2000) 5087–5159.  arXiv:hep-ph/0102033, doi:10.1016/S0217-751X(00)00215-7.
  * (27) I. Dubovyk, A. Freitas, J. Gluza, T. Riemann, J. Usovitsch, The two-loop electroweak bosonic corrections to $\sin^{2}\theta_{\rm eff}^{\rm b}$, Phys. Lett. B762 (2016) 184–189.  arXiv:1607.08375, doi:10.1016/j.physletb.2016.09.012.
  * (28) I. Dubovyk, A. Freitas, J. Gluza, T. Riemann, J. Usovitsch, Complete electroweak two-loop corrections to $Z$ boson production and decay, Phys. Lett. B783 (2018) 86–94.  arXiv:1804.10236, doi:10.1016/j.physletb.2018.06.037.
  * (29) T. Binoth, G. Heinrich, An automatized algorithm to compute infrared divergent multi-loop integrals, Nucl. Phys. B585 (2000) 741–759.  arXiv:hep-ph/0004013v.2.
  * (30) S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner, J. Schlenk, T. Zirke, pySecDec: a toolbox for the numerical evaluation of multi-scale integrals, Comput. Phys. Commun. 222 (2018) 313–326.  arXiv:1703.09692, doi:10.1016/j.cpc.2017.09.015.
  * (31) S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner, J. Schlenk, A GPU compatible quasi-Monte Carlo integrator interfaced to pySecDec, Comput. Phys. Commun. 240 (2019) 120–137.  arXiv:1811.11720, doi:10.1016/j.cpc.2019.02.015.
  * (32) A. V. Smirnov, M. N. Tentyukov, Feynman Integral Evaluation by a Sector decomposiTion Approach (FIESTA), Comput. Phys. Commun. 180 (2009) 735–746.  arXiv:0807.4129, doi:10.1016/j.cpc.2008.11.006.
  * (33) A. V. Smirnov, FIESTA 4: Optimized Feynman integral calculations with GPU support, Comput. Phys. Commun. 204 (2016) 189–199.  arXiv:1511.03614, doi:10.1016/j.cpc.2016.03.013.
  * (34) J. Gluza, K. Kajda, T. Riemann, AMBRE: A Mathematica package for the construction of Mellin-Barnes representations for Feynman integrals, Comput. Phys. Commun. 177 (2007) 879–893.  arXiv:0704.2423, doi:10.1016/j.cpc.2007.07.001.
  * (35) M. Czakon, Automatized analytic continuation of Mellin-Barnes integrals, Comput. Phys. Commun. 175 (2006) 559–571, mathematica program MB.m version 1.2 (Jan 2, 2009), available at the MB Tools webpage, http://projects.hepforge.org/mbtools/.  arXiv:hep-ph/0511200, doi:10.1016/j.cpc.2006.07.002.
  * (36) J. Usovitsch, Numerical evaluation of Mellin-Barnes integrals in Minkowskian regions and their application to two-loop bosonic electroweak contributions to the weak mixing angle of the $Z{\bar{b}}b$-vertex, Ph.D. thesis, Humboldt-Universität zu Berlin, Inst. Math. (2018), http://dx.doi.org/10.3204/PUBDB-2017-12744, http://dx.doi.org/10.3204/PUBDB-2018-05160.
  * (37) I. Dubovyk, J. Gluza, T. Riemann, J. Usovitsch, Numerical integration of massive two-loop Mellin-Barnes integrals in Minkowskian regions, PoS LL2016 (2016) 034. https://pos.sissa.it/260/034/pdf.  arXiv:1607.07538, doi:10.22323/1.260.0034.
  * (38) J. Usovitsch, I. Dubovyk, T. Riemann, MBnumerics: Numerical integration of Mellin-Barnes integrals in physical regions, PoS LL2018 (2018) 046.  arXiv:1810.04580, doi:10.22323/1.303.0046.
  * (39) M. K. Mandal, X. Zhao, Evaluating multi-loop Feynman integrals numerically through differential equations, JHEP 03 (2019) 190.  arXiv:1812.03060, doi:10.1007/JHEP03(2019)190.
  * (40) M. L. Czakon, M. Niggetiedt, Exact quark-mass dependence of the Higgs-gluon form factor at three loops in QCD, JHEP 05 (2020) 149.  arXiv:2001.03008, doi:10.1007/JHEP05(2020)149.
  * (41) X. Liu, Y.-Q. Ma, C.-Y. Wang, A Systematic and Efficient Method to Compute Multi-loop Master Integrals, , Phys. Lett. B 779 (2018) 353–357.  arXiv:1711.09572, doi:10.1016/j.physletb.2018.02.026.
  * (42) X. Liu, Y.-Q. Ma, Determining arbitrary Feynman integrals by vacuum integrals, Phys. Rev. D 99 (7) (2019) 071501.  arXiv:1801.10523, doi:10.1103/PhysRevD.99.071501.
  * (43) C. Brønnum-Hansen, C.-Y. Wang, Contribution of third generation quarks to two-loop helicity amplitudes for W boson pair production in gluon fusion, JHEP 01 (2021) 170.  arXiv:2009.03742, doi:10.1007/JHEP01(2021)170.
  * (44) In principle, it is also applicable at higher loop orders, but we have not studied any examples beyond three loops.
  * (45) A. Kotikov, Differential equations method: New technique for massive Feynman diagrams calculation, Phys. Lett. B254 (1991) 158–164.  doi:10.1016/0370-2693(91)90413-K.
  * (46) A. Kotikov, Differential equations method: The Calculation of vertex type Feynman diagrams, Phys. Lett. B259 (1991) 314–322.
  * (47) A. V. Kotikov, Differential equation method: The Calculation of N point Feynman diagrams, Phys. Lett. B267 (1991) 123–127.
  * (48) E. Remiddi, Differential equations for Feynman graph amplitudes, Nuovo Cim. A110 (1997) 1435–1452.  arXiv:hep-th/9711188.
  * (49) J. M. Henn, Multiloop integrals in dimensional regularization made simple, Phys. Rev. Lett. 110 (25) (2013) 251601.  arXiv:1304.1806, doi:10.1103/PhysRevLett.110.251601.
  * (50) M. Prausa, epsilon: A tool to find a canonical basis of master integrals, Comput. Phys. Commun. 219 (2017) 361–376.  arXiv:1701.00725, doi:10.1016/j.cpc.2017.05.026.
  * (51) O. Gituliar, V. Magerya, Fuchsia: a tool for reducing differential equations for Feynman master integrals to epsilon form, Comput. Phys. Commun. 219 (2017) 329–338.  arXiv:1701.04269, doi:10.1016/j.cpc.2017.05.004.
  * (52) A. von Manteuffel, L. Tancredi, A non-planar two-loop three-point function beyond multiple polylogarithms, JHEP 06 (2017) 127.  arXiv:1701.05905, doi:10.1007/JHEP06(2017)127.
  * (53) L. Adams, E. Chaubey, S. Weinzierl, Simplifying Differential Equations for Multiscale Feynman Integrals beyond Multiple Polylogarithms, Phys. Rev. Lett. 118 (14) (2017) 141602.  arXiv:1702.04279, doi:10.1103/PhysRevLett.118.141602.
  * (54) J. Henn, B. Mistlberger, V. A. Smirnov, P. Wasser, Constructing d-log integrands and computing master integrals for three-loop four-particle scattering, JHEP 04 (2020) 167.  arXiv:2002.09492, doi:10.1007/JHEP04(2020)167.
  * (55) C. Dlapa, J. Henn, K. Yan, Deriving canonical differential equations for Feynman integrals from a single uniform weight integral, JHEP 05 (2020) 025.  arXiv:2002.02340, doi:10.1007/JHEP05(2020)025.
  * (56) J. Usovitsch, Factorization of denominators in integration-by-parts reductions, arXiv:2002.08173.
  * (57) A. V. Smirnov, V. A. Smirnov, How to choose master integrals, Nucl. Phys. B 960 (2020) 115213.  arXiv:2002.08042, doi:10.1016/j.nuclphysb.2020.115213.
  * (58) E. Panzer, Feynman integrals and hyperlogarithms, Ph.D. thesis, Humboldt U. (2015).  arXiv:1506.07243, doi:10.18452/17157.
  * (59) A. von Manteuffel, E. Panzer, R. M. Schabinger, A quasi-finite basis for multi-loop feynman integrals, Journal of High Energy Physics 2015 (2).  doi:10.1007/jhep02(2015)120.
URL http://dx.doi.org/10.1007/JHEP02(2015)120

  * (60) A. von Manteuffel, E. Panzer, R. M. Schabinger, Computation of form factors in massless qcd with finite master integrals, Physical Review D 93 (12).  doi:10.1103/physrevd.93.125014.
URL http://dx.doi.org/10.1103/PhysRevD.93.125014

  * (61) A. von Manteuffel, C. Studerus, Reduze 2 - Distributed Feynman Integral Reduction. arXiv:1201.4330.
  * (62) P. Maierhöfer, J. Usovitsch, P. Uwer, Kira - a Feynman integral reduction program, Comput. Phys. Commun. 230 (2018) 99–112.  arXiv:1705.05610, doi:10.1016/j.cpc.2018.04.012.
  * (63) J. Klappert, F. Lange, P. Maierhöfer, J. Usovitsch, Integral reduction with Kira 2.0 and finite field methods, Comput. Phys. Commun. 266 (2021) 108024.  arXiv:2008.06494, doi:10.1016/j.cpc.2021.108024.
  * (64) J. Klappert, F. Lange, Reconstructing rational functions with FireFly, Comput. Phys. Commun. 247 (2020) 106951.  arXiv:1904.00009, doi:10.1016/j.cpc.2019.106951.
  * (65) J. Klappert, S. Y. Klein, F. Lange, Interpolation of dense and sparse rational functions and other improvements in FireFly, Comput. Phys. Commun. 264 (2021) 107968.  arXiv:2004.01463, doi:10.1016/j.cpc.2021.107968.
  * (66) S. Pozzorini, E. Remiddi, Precise numerical evaluation of the two loop sunrise graph master integrals in the equal mass case, Comput. Phys. Commun. 175 (2006) 381–387.  arXiv:hep-ph/0505041, doi:10.1016/j.cpc.2006.05.005.
  * (67) F. Moriello, Generalised power series expansions for the elliptic planar families of Higgs + jet production at two loops, JHEP 01 (2020) 150.  arXiv:1907.13234, doi:10.1007/JHEP01(2020)150.
  * (68) R. Bonciani, V. Del Duca, H. Frellesvig, J. M. Henn, M. Hidding, L. Maestri, F. Moriello, G. Salvatori, V. A. Smirnov, Evaluating a family of two-loop non-planar master integrals for Higgs + jet production with full heavy-quark mass dependence, JHEP 01 (2020) 132.  arXiv:1907.13156, doi:10.1007/JHEP01(2020)132.
  * (69) See also Ref. Fael:2021kyg  for a similar approach for single-scale problems.
  * (70) M. Hidding, DiffExp, a Mathematica package for computing Feynman integrals in terms of one-dimensional series expansions, Comput. Phys. Commun. 269 (2021) 108125\.  arXiv:2006.05510, doi:10.1016/j.cpc.2021.108125.
  * (71) See Refs. Song:2021vru ; Liu:2021wks  for recent independent efforts on electroweak two-loop box integrals.
  * (72) In this and the following examples, the large number of massless propagators occurs because all SM fermions except the top quark are taken to be massless.
  * (73) M. Borinsky, Tropical Monte Carlo quadrature for Feynman integrals, arXiv:2008.12310.
  * (74) B. Ananthanarayan, S. Banik, S. Friot, S. Ghosh, Multiple Series Representations of N-fold Mellin-Barnes Integrals, Phys. Rev. Lett. 127 (15) (2021) 151601.  arXiv:2012.15108, doi:10.1103/PhysRevLett.127.151601.
  * (75) P. Banerjee, et al., Theory for muon-electron scattering @ 10 ppm: A report of the MUonE theory initiative, Eur. Phys. J. C 80 (6) (2020) 591.  arXiv:2004.13663, doi:10.1140/epjc/s10052-020-8138-9.
  * (76) M. Fael, F. Lange, K. Schönwald, M. Steinhauser, A semi-analytic method to compute Feynman integrals applied to four-loop corrections to the $\overline{\rm MS}$-pole quark mass relation, arXiv:2106.05296, doi:10.1007/JHEP09(2021)152.
  * (77) Q. Song, A. Freitas, On the evaluation of two-loop electroweak box diagrams for $e^{+}e^{-}\to HZ$ production, JHEP 04 (2021) 179.  arXiv:2101.00308, doi:10.1007/JHEP04(2021)179.
  * (78) X. Liu, Y.-Q. Ma, Multiloop corrections for collider processes using auxiliary mass flow, arXiv:2107.01864.