# The accretion of a solar mass per day by a 17-billion solar mass black hole

Christian Wolf, Samuel Lai, Christopher A. Onken, Neelesh Amrutha, Fuyan

Bian, Wei Jeat Hon, Patrick Tisserand, Rachel L. Webster

► To cite this version:

Christian Wolf, Samuel Lai, Christopher A. Onken, Neelesh Amrutha, Fuyan Bian, et al.. The accretion of a solar mass per day by a 17-billion solar mass black hole. Nature Astronomy, 2024, 8, pp.520-529. 10.1038/s41550-024-02195-x. insu-04851568

HAL Id: insu-04851568

https://insu.hal.science/insu-04851568v1

Submitted on 4 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The accretion of a solar mass per day by a 17-billion solar mass black hole

Received: 19 April 2023

Accepted: 4 January 2024

Published online: 19 February 2024

<figure><img src="image_2.png" /></figure>

Check for updates

Christian Wolf Ⓓ1,2✉, Samuel Lai Ⓓ1, Christopher A. Onken Ⓓ1, Neelesh Amrutha Ⓓ1, Fuyan Bian³, Wei Jeat Hon Ⓓ4, Patrick Tisserand Ⓓ5 & Rachel L. Webster Ⓓ4

Around a million quasars have been catalogued in the Universe by probing deeper and using new methods for discovery. However, the hardest ones to find seem to be the rarest and brightest specimens. Here we study the properties of the most luminous of all quasars found so far. These have been overlooked until recently, which demonstrates that modern all-sky surveys have much to reveal. The black hole in this quasar accretes around one solar mass per day onto an existing mass of ~17 billion solar masses. In this process, the accretion disk alone releases a radiative energy of $2 \times 10^{41}$ W. If the quasar is not strongly gravitationally lensed, then its broad-line region is expected to have the largest physical and angular diameter occurring in the Universe and this will allow the Very Large Telescope Interferometer to image its rotation and measure its black-hole mass directly. This will be an important test for broad-line region size-luminosity relationships, whose extrapolation has underpinned common black-hole mass estimates at high redshift.

In 1963, Maarten Schmidt identified the first quasar<fnref n="1" />, known as 3C 273. It appeared as a remarkably bright star of 12th magnitude, and its redshift suggested that it was among the most distant objects known in the Universe at the time. The two facts together implied an implausibly huge output of light and, ever since then, newly found quasars have impressed with their immense energy release from a small region of space. This could only be explained as gravitational energy being converted into heat and light within a highly viscous accretion disk around a super-massive black hole²,³ (SMBH). Quasars are signposts of fast growth in SMBHs on public display and allow the study of these growth processes.

Finding large samples of quasars then provides population and growth statistics to explain the origin of SMBHs in the Universe⁴⁻⁷. Generally, the most luminous quasars contain the fastest-growing SMBHs, although the relationship between mass accretion rate and luminosity is affected by the mass and spin of the black hole as well as the structure and viewing angle of the accretion disk and disk winds⁸⁻¹⁴.

Today, around a million quasars are known¹⁵, although some specimens stand out from the crowd: in 2015, the ultra-luminous quasar J0100+2802 at redshift $z = 6.3$ was identified¹⁶ with an SMBH of 10 billion

solar masses¹⁷. In 2018, an even more luminous object¹⁸ was found, J2157-3602 at $z = 4.7$, with an SMBH of 24 billion solar masses¹⁴. Although their luminosity implies rapid growth, their existence is hard to explain: when black holes start from the remnant of a stellar collapse and grow episodically within the Eddington limit, they are not expected to reach the evident masses in the time from the Big Bang to the epoch of their observation, which has triggered a search for alternative scenarios¹⁹⁻²².

Although exceptionally rare, the most extremely luminous quasars are interesting for several reasons beyond their intrinsic nature, as discussed later. In this paper, we present the properties of the recently discovered²³ quasar SMSS J052915.80-435152.0, hereafter called J0529-4351, which is a 16th magnitude object at redshift $z = 3.962$ (Fig. 1), and we reveal it to be the most luminous quasar currently known in the Universe (Fig. 2).

## Results

When quasars appear extremely bright, it may be suspected that their observed brightness is magnified by gravitational lensing from a massive galaxy on the line of sight. Strong lensing causes multiple separate

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹Research School of Astronomy and Astrophysics, Australian National University, Canberra, Australian Capital Territory, Australia. ²Centre for Gravitational Astrophysics, Australian National University, Canberra, Australian Capital Territory, Australia. ³European Southern Observatory, Santiago, Chile. ⁴School of Physics, University of Melbourne, Parkville, Victoria, Australia. ⁵Sorbonne Universités, CNRS, UMR 7095, Institut d'Astrophysique de Paris, Paris, France.</footnote>

<footnote>e-mail: christian.wolf@anu.edu.au</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1** | The exceptional quasar J0529–4351. This image is from the Dark Energy Camera Legacy Survey DR10 with size 20 × 20 arcsec²; north is up and east is to the left; a neighbouring M star appears in red.</figcaption></figure>

images of a quasar in the sky<sup>24,25</sup>. Two other quasars with a redshift and apparent brightness similar to J0529–4351 are known to be strongly magnified by lensing: the double-image APM 08279+5255 at z = 3.91, with a separation of 0.5 arcsec (ref. 26), and the quadruply imaged B 1422+231 at z = 3.62, with separations up to 1.5 arcsec (ref. 27). Estimated magnification factors for these two objects range from 40 to 100, which implies that these quasars are not intrinsically extreme, but are members of the bulk population<sup>26,28</sup>.

J0529–4351, in contrast, shows no sign of strong lensing; data from the European Space Agency (ESA) Gaia satellite suggest that it is a point source, in terms of object morphology and astrometric excess noise (AEN), which has been used to find dual quasars or lensed quasars that appear unresolved to Gaia<sup>29</sup> (see Fig. 2 and online methods for more detail). We also find no strong foreground absorber system, which our high signal-to-noise spectrum probes in Mg II at z > 1.15; the strongest system, at z = 2.118, has an equivalent width of EW(2796) = 0.8 Å, which suggests an impact parameter of 20 kpc (or 2.4 arcsec) from the line of sight<sup>30</sup> to the quasar, whereas a plausible image separation in a lensed scenario is ~0.1 arcsec. We can also estimate the probability of lensing a source at z ≈ 4 using common models for the galaxy distribution as isothermal masses<sup>31,32</sup>, and obtain *p* ≈ 1.3 × 10<sup>-3</sup>. An image separation of 0.2 arcsec or less reduces it further to *p* ≈ 2 × 10<sup>-4</sup>. These estimates change by only a factor of a few when changing the input galaxy velocity dispersion function or the bright end of the quasar luminosity function. Even if a very steep intrinsic quasar luminosity function at the bright end enhanced the magnification bias, the probability that this source is strongly lensed will be less than 1%. We thus take the strength of the quasar emission at face value, although final confirmation from a high-resolution space-based or adaptive-optics image would still be desirable.

We also investigate the recent history in the brightness of J0529-4351 to see whether it may have been previously overlooked due to extreme variability. The 0.5 m NASA (National Aeronautics and Space Administration) Asteroid Terrestrial-impact Last Alert System<sup>33</sup> (ATLAS) telescope has been providing a light curve since mid-2017 (Fig. 3). Brightness variations of ~15% have been found over the last six years, which are not unexpected for luminous quasars. Earlier records, from photographic plates observed in 1980 and 1998 (measured by the SuperCOSMOS Sky Surveys), found the R-band brightness to be consistent with recent observations by the SkyMapper Southern Survey<sup>34</sup> (SMSS). The somewhat longer light curve from the Wide-field Infrared Survey Explorer<sup>35</sup> (WISE) shows also only modest variability.

The quasar is undetected in the Rapid ASKAP Continuum Survey<sup>36</sup> (RACS), and hence has a flux of less than 1 µJy (4σ) in the broad 887 MHz band. Using a common definition of radio loudness<sup>37</sup>, we find the object to be safely in the radio-quiet regime (R < 1). Thus, we have no reason to suspect that its luminosity is affected by jet emission, let alone being relativistically boosted.

For further analysis, we use a spectrum of optical and near-infrared light (Fig. 4) from the 8.2 m Very Large Telescope (VLT) at the European Southern Observatory (ESO). We split the spectrum with a publicly available spectrum-fitting code<sup>38</sup> into an accretion disk continuum and emission-line contributions. From the disk continuum we quantify the disk luminosity of the quasar and get a first proxy for the accretion rate. We find monochromatic luminosities of log($\lambda L_{135}/\text{erg s}^{-1}$) = 47.93 and log($\lambda L_{300}/\text{erg s}^{-1}$) = 47.76 at rest-frame wavelengths of $\lambda$ = 135 and 300 nm, respectively. An approximate estimate of the full radiative output from the accretion disk using standard bolometric corrections<sup>39</sup> and a 0.75 anisotropy factor<sup>40</sup> yields log($L_{\text{bol}}/\text{erg s}^{-1}$) = 48.37. Using a standard radiative efficiency value<sup>41</sup> of 0.1, this translates to an accretion rate of ~413 solar masses per year. This result makes J0529-4351 the most luminous quasar and, by inference, the fastest-growing black hole in the Universe known to date (in terms of mass growth per unit time).

A more refined analysis of the bolometric luminosity and accretion rate requires modelling the possible spectra of the disk continuum over a grid of black-hole masses and spins and simulating their observation for different accretion rates and viewing angles<sup>42–45</sup>. Using slim-disk models as appropriate for SMBHs with high accretion rates<sup>46,47</sup> and a publicly available code<sup>14,48</sup> (online methods), we find a best-fitting solution to our spectrum being an Eddington-accreting disk around an SMBH with a moderate spin of $a \approx 0.4$ viewed at an intermediate angle of $i \approx 45^\circ$ (Fig. 5). The intrinsically emitted bolometric luminosity is log($L_{\text{bol}}/\text{erg s}^{-1}$) = 48.16, which is 62% of the value derived above; the lower value mostly results from generic bolometric corrections being overestimates for very massive black holes with their colder disks and suppressed ultraviolet (UV) emission. The best-fit slim-disk model yields an accretion rate of 370 $M_\odot$ per year and a radiative efficiency of ~0.09. Given the lack of higher-frequency data in the spectral energy distribution (SED), the integrated model SED may still underpredict the bolometric flux, whereas the standard bolometric corrections must be overcorrecting for black holes of over a billion solar masses. The best estimate is expected to lie between the two, and we thus choose the average of log($L_{\text{bol}}/\text{erg s}^{-1}$) = 48.27.

When modelling SEDs for a fixed mass accretion rate, different black-hole spins lead to different bolometric corrections and radiative efficiencies, whereas the monochromatic UV-optical luminosities are only modestly affected. The main uncertainty arises from the unknown viewing angle, which affects both the apparent monochromatic luminosity and the derived accretion rate. Given the broad confidence intervals for spin and viewing angle in the model fit, we consider the full plausible angle range from pole-on to $i = 60^\circ$ as a 95% error margin; this corresponds to a 95% range in luminosity of $\pm 0.12$ dex and an accretion rate of 280 to 490 solar masses per year.

## Black-hole mass

We estimate the mass of the black hole powering this quasar using two fundamentally different methods.

1. Assuming that the continuum emission is affected at its blue end by the inner truncation of the accretion disk due to the innermost stable circular orbit around the black hole<sup>8</sup>, we can infer a combined estimate for mass and spin of the black hole. More massive black holes impose larger truncation radii and move the peak of the continuum emission to cooler temperatures and longer wavelengths<sup>49</sup>. With this method, the mass of the black hole in 3C 273 was found to be between 200 and 500 million solar masses<sup>50</sup>, 35 years before a value of ~300 million solar masses
<figure><img src="image_5.png" /><figcaption>**Fig. 2 | J0529-4351 in relation to other quasars.** Left: J0529-4351 versus the population from SDSS and from XQ-100, which has the highest-quality spectra⁵⁹ for luminous quasars at z = 4. The luminosity axis shows monochromatic luminosity with standard bolometric corrections applied. Right: AEN of</figcaption></figure>

J0529-4351 compared with known bright quasars of Gaia magnitude $G_{\text{RP}} < 18$ at $z > 1$; the known Gaia-unresolved lensed or dual quasars are all at AEN > 1 mas. Zero astrometric excess noise was replaced with a value of 0.01.

<figure><img src="image_6.png" /><figcaption>**Fig. 3 | Light curves of J0529-4351.** Top: data from NASA ATLAS are combined into weekly (orange passband) or moon-period (cyan passband) median averages of typically eight data points. Bottom: data from WISE are median averages of typically 25 visits completed within less than two weeks. Error bars show interquartile ranges within the averaging period.</figcaption></figure>

was measured by interferometric observations⁵¹ of the broad emission-line region in 3C 273. Since then, the continuum-fitting method has not only proved useful for estimating masses using thin-disk models⁴³⁻⁴⁵ but has also been applied¹⁴ with slim-disks models⁴⁶ that are expected to be a more realistic description of the near-Eddington accretion disks of fast-growing SMBHs. From the continuum shape of J0529-4351, we find a mass of $log M/M_{\odot} = 10.28^{+0.17}_{-0.10}$ (or ~19 billion solar masses).

(2) Assuming that the broad emission-line region represents virialized gas moving at the velocity of Keplerian orbits around the SMBH, we can use the width of emission lines and the continuum luminosity to infer the mass of the black hole. This method is known as the virial single-epoch method⁵²,⁵³ and has also been used to estimate the mass of the SMBH in J0100+2802 (ref. 17) at $z=6.3$. Its application at the high-luminosity end of the quasar population is mainly limited by a lack of independent calibrations for the SMBH mass and relies on extrapolations from lower-luminosity quasars calibrated with reverberation mapping⁵⁴.

Our spectra offer two emission lines for this method: the triple-ionized carbon line C IVλ1549 and the singly-ionized magnesium line Mg IIλ2799. The C IVλ1549 line appears asymmetric and blueshifted relative to the Mg IIλ2799 line ($\Delta\nu = -3,120 \pm 80 \text{ km s}^{-1}$) and a line full-width at half-maximum (FWHM) of $7,245 \pm 175 \text{ km s}^{-1}$. The Mg

IIλ2799 line turns out to be hard to calibrate at this redshift because of atmospheric absorption, but we do measure a FWHM of 4,395 ± 435 km s⁻¹. Using calibrations commonly applied to high-luminosity quasars, the line properties and luminosity translate into SMBH mass estimates from $log M/M_{\odot} = 10.03 \pm 0.06$ to $10.45 \pm 0.02$ for the C IV line⁵³,⁵⁵ and from $log M/M_{\odot} = 10.03 \pm 0.09$ to $10.36 \pm 0.09$ for the Mg II line⁵⁶,⁵⁷ (errors are standard deviations due to propagated observational uncertainties and the scatter among the values is in line with large systematic calibration uncertainties). Table 1 summarizes all our estimates of masses and luminosities.

The line-based estimates are consistent with each other and with the continuum-based estimate. Given that the systematic uncertainties in these methods are larger than the statistical error propagation and may be as high as 0.4 dex, we combine the mean mass estimates of the two lines and that from the SED without weighting into a final result for the mass of $log M/M_{\odot} = 10.24 \pm 0.02$ (simple mean and standard deviation), or ~17 billion solar masses. This also implies that the SMBH is accreting near the Eddington limit (Eddington ratio of ~0.9).

## Discussion

In terms of luminosity and likely growth rate, J0529-4351 is the most extreme quasar known. The accretion of J0529-4351 is near the Eddington limit, which is common among quasars of the highest luminosity¹⁶,⁵⁸,⁵⁹. The growth rates are mostly uncertain owing to the unknown viewing angle. Assuming persistent accretion at its current Eddington ratio, the mass doubling time is ~30 Myr. However, with 17 billion solar masses, the black hole in J0529-4351 at $z=3.962$ is not the largest SMBH found in the most luminous quasars. It has over 50% more mass than the black hole in J0100+2802 at $z=6.3$ but one third less mass than that in J2157-3602 at $z=4.7$. The age of the Universe at these three redshifts is 859 Myr, 1,244 Myr and 1,530 Myr, assuming a flat Universe with a concordance cosmology (a cosmological constant of $\Omega_{\Lambda} = 0.3$). Given that J0529-4351 is observed at a later epoch in the Universe than J2157-4351 and J0100+2802, it is less of a challenge to models of early SMBH growth.

It may be tempting to speculate on the spin of extreme SMBHs: although it has been argued that growing black holes should spin up over time, this would also increase radiative feedback and slow down accretion, making it harder to grow the most massive SMBHs within the age of the Universe. Instead, one way of growing black holes from stellar seeds to the greatest masses we measure is 'chaotic accretion' with randomly changing orientation that keeps the black-hole spins and radiative feedback low²².

A long-standing question has been what mechanism fuels the high accretion rate, which also must have persisted for some time previously, although not necessarily in the immediate past. The mass
<figure><img src="image_7.png" /><figcaption>**Fig. 4** | Spectrum of the J0529-4351 from the ESO 8.2 m VLT. Top: observed-frame wavelength range of 0.42–2.27 µm; a power law $f_v \propto v^{-0.5}$ is shown as a dashed line for comparison. Centre: decomposition, emission-line fits and noise level (grey line, very low in the C IVλ1549 part) lead to derived black-hole masses</figcaption></figure>

of $\log M/M_{\odot} = 10.24$ for C IVλ1549 (left) and 10.20 for Mg IIλ2799 (right). Bottom: earlier optical-only spectra from AllBRICQS and from Gaia DR3 in arbitrary units.

and kinematics of the gas in the host galaxy of J0529-4351 could be observed with the Atacama Large Millimetre Array (ALMA). ALMA has already revealed the largest spiral galaxy in the early Universe$^{60}$, in the quasar BRI 1335-0417 at z = 4.4, which is approximately ten times less luminous than J0529-4351. In contrast, ALMA observations have revealed a merger signature in the host galaxy of a quasar at z = 7.54 (ref. 61). Further observations of extremely luminous quasars are progressing$^{62}$. If extreme quasars were caused by unusual host galaxy gas flows, ALMA should see this; if nothing unusual was found in the host

gas, then this would sharpen the well-known puzzle of how to sustain high accretion rates for long enough to form such extreme SMBHs.

**Have we overlooked still more extreme quasars?**

Finding rare and exceptionally bright quasars such as J0529-4351 does not require large telescopes, but is instead a needle-in-the-haystack problem that needs precise data with discriminative power across large areas of sky. Such objects are often hiding in plain sight and are mistaken for stars. For example, the quasar J1144-4308 at redshift
<figure><img src="image_8.png" /></figure>

**Fig. 5 |** Results of the continuum fitting for J0529−4351 using Markov chain Monte Carlo and slim-disk models. Top right: spectrum (blue line), continuum (blue squares) and SED fit (red line). Smaller panels: parameter matrix with posteriors. Estimates are given as the peak of the probability distribution and confidence intervals are the central 68-percentiles. The highest likelihood model

is a disk viewed at an inclination angle ($\theta_{\text{inc}}$) of ~45° around an intermediate-spin SMBH with ~19 billion solar masses and an accretion rate of 370 $M_{\odot}$ per year; the constraints around viewing angle and spin are, however, weak. Spin is denoted by $a$, and $L_{\text{disk}}$ is in units of Eddington luminosity.

$z = 0.83$ was recognized as one of the brightest quasars in the night sky (at 14 mag) only in 2022, despite having been imaged on photographic plates since the nineteenth century⁶³.

In the 1960s, finding the first quasars was driven by radio detections of objects that otherwise appeared like regular nearby stars in our Milky Way galaxy. Later quasar surveys employed colour selection, template-fitting and Bayesian methods⁶⁴⁻⁷². But the main barrier for obtaining complete samples has always been contamination of quasar candidate selections with stars from our Milky Way, which appear similar if no discriminating information is available and vastly outnumber true quasars in the extreme regime of the brightest quasars. With a million quasars known by now and the wealth of data from modern all-sky surveys, machine-learning approaches are now most popular⁷³⁻⁷⁵. These tend to get the classification right for the majority of objects, although training samples imply that they perform less well for rare extreme cases. In the case of extremely luminous quasars, the obvious bias of a training sample is simply that they do not seem to exist until they are found.

A low-resolution spectrum of J0529-4351 revealing its quasar nature and redshift is part of the public all-sky Gaia DR3 data set

published on 13 June 2022. The machine-learning classification of this data set by the Gaia Discrete Source Classifier (DSC) has assigned J0529-4351 a 99.98% probability of being a Milky Way star⁷⁶, although a human astronomer looking at the Gaia spectrum would recognize the quasar and redshift at first sight.

These days, combining data from Gaia and WISE makes an all-sky search for bright quasars straightforward. An absence of parallax and proper motions removes most of the bright galactic stars from the search, and the mid-infrared photometry from WISE is then sufficient to discriminate the disk continuum and dust emission of quasars from the Rayleigh-Jeans tail of stellar photospheres⁶⁹. The All-sky Bright Complete Quasar Survey²³ (AllBRICQS), which reported J0529-4351, selects its candidates with a simple heuristic WISE colour cut of $W1 - W2 > 0.3$ and an absence of >4$\sigma$ evidence for stellar motion. Although it aims for completeness and its recall of known quasars is nearly perfect, its efficiency is still high, given that 97% of its candidates turn out to be quasars²³. Given that the AllBRICQS follow-up is complete at $G_{\text{RP}} < 16$ in both hemispheres, we doubt that a quasar of higher UV-optical luminosity will be found in the future, unless it is hiding behind the galactic plane: AllBRICQS did not search at galactic latitude of $|b| < 10°$,
<table><caption>**Table 1 | Luminosity, black-hole mass and accretion rate for the quasar J052915.80-435152.0**</caption><thead><tr><th>Luminosity</th><th>log L/(ergs<sup>-1</sup>)</th></tr></thead><tbody><tr><td>λL<sub>135</sub></td><td>47.93</td></tr><tr><td>λL<sub>300</sub></td><td>47.76</td></tr><tr><td>L<sub>bol</sub> (from 135 nm)</td><td>48.39</td></tr><tr><td>L<sub>bol</sub> (from 300 nm)</td><td>48.35</td></tr><tr><td>L<sub>bol</sub> (slim-disk model)</td><td>48.16</td></tr><tr><td>L<sub>bol</sub> (best estimate)</td><td>48.27±0.06</td></tr><tr><th>Black hole mass</th><th>log M/M<sub>☉</sub></th></tr><tr><td>C ivλ1549 line (average)</td><td>10.24±0.15</td></tr><tr><td>Mg iiλ2799 line (average)</td><td>10.20±0.08</td></tr><tr><td>Continuum fit (slim disk<sup>47</sup>)</td><td>10.28<sup>+0.17</sup><sub>-0.10</sub></td></tr><tr><td>All methods (best estimate)</td><td>10.24±0.02</td></tr><tr><th>Mass accretion rate</th><th>log ˙M/(M<sub>☉</sub> yr<sup>-1</sup>)</th></tr><tr><td>for plausible inclination range</td><td>2.57±0.06</td></tr></tbody></table>

Given are best-fit estimates and 68% confidence intervals. Systematic uncertainties for the mass estimates are proposed to be as high as 0.4 dex.

where source crowding and higher dust extinction makes the discovery of quasars still very challenging<sup>77</sup>.

## Future observations with bright quasars

Extremely luminous quasars enable further observations for specific challenging quests: after decades of instrumental development to increase the spatial resolution of telescopic observations, the VLTI/ GRAVITY Collaboration<sup>51</sup> made headlines with their spatiokinematic mapping of the broad-line region (BLR) in the iconic nearby quasar 3C 273 (redshift $z = 0.157$). This observation revealed the orbital rotation of the disk-like BLR in 3C 273 in a spatially resolved pattern and thus provided a direct measurement of the black-hole mass from the BLR orbits. The measurement was possible as 3C 273 was bright enough and displayed its Paschen-α emission line in GRAVITY's K-band window.

Although the redshift of J0529-4351 imposed challenges for the work presented here, as nearly all strong UV emission lines are in places heavily affected by atmospheric absorption, it does place its Hβ line at 2.4 µm in the K-band window of the GRAVITY instrument. We predict the size of the Hβ-emitting BLR in J0529-4351 by extrapolating the radius-luminosity relationship<sup>78</sup> and find a radius of 2.2 pc. This implies an angular diameter of 0.64 milli-arcseconds, which is an order of magnitude larger than the BLR in 3C 273 and thus the largest-appearing quasar BLR in the Universe. We expect the soon-to-be-upgraded VLTI/GRAVITY+ to obtain a superbly well-resolved picture of the BLR rotation in J0529-4351 and, therefore, a much more reliable measurement of its black-hole mass. Crucially, black-hole masses at the high end are estimated from relationships that have been extrapolated by orders of magnitude such that the whole scale for objects like this is at risk. Getting a direct mass measurement for a black hole that is likely to have 50 times the mass of the black hole in 3C 273 would be extremely valuable for constraining the relationships commonly used to estimate the masses of black holes in the early Universe. The true mass scale of the earliest SMBH would also affect the answer to the question of how hard exactly it is to form them.

Future plans include watching the Universe expand with repeated observations of quasars lasting for a decade: the expansion shifts the redshift of individual gas clouds observed as Lyman-α forest absorption lines in quasar spectra<sup>79-81</sup> (Sandage test). The signal of such observing campaigns depends on the availability of a sufficient number of very bright quasars and will still require advanced facilities such as the forthcoming ESO Extremely Large Telescope (ELT). J0529-4351 will quite obviously be an important part of this long-term endeavour.

## Methods

### Evidence for point-source geometry

Although Gaia provides no images as such, we use the source characterization from its Data Release 3 to provide evidence for a point-source geometry in J0529-4351. First, there is no noticeable flux outside an aperture with 0.175 arcsec radius around the object centroid in the high-resolution imaging<sup>82</sup>, as determined by the excess ratio from comparing the fluxes in the BP and RP channels with the G-band aperture<sup>23</sup>. The closest neighbouring object is 3.2 arcsec away; it is not detected by Gaia, but by the DECam component of the DESI Legacy Imaging Surveys<sup>83</sup>, with $g = 23.581 \pm 0.039$, $r = 22.319 \pm 0.015$, $i = 20.686 \pm 0.007$ and $z = 20.099 \pm 0.007$ in DR10, which is a good match to the colours of an M4 star in the Milky Way.

Second, astrometric excess noise is absent in J0529-4351, although 319 good observations were available to determine astrometric solutions. All Gaia parameters that may hint at unresolved multiple sources, including frac_multipeak (0), astrometric_excess_noise (0), astrometric_sigma5d_max (0.069), ruwe (1.013) and ipd_gof_harmonic_amplitude (0.0167) are in line with the average values (0.132, 0.042, 0.074, 1.015 and 0.045) for bright quasars with $G_{RP} < 16$ at redshift $z > 2$ and far below any selection cuts used to search for Gaia-unresolved binary stars<sup>84</sup> and dual or lensed quasars<sup>29</sup>. Figure 2 shows the astrometric excess noise of quasars at $G_{RP} < 18$ including all dual and lensed quasars<sup>25,29</sup> with image separations below 0.5 arcsec from the list at ref. 85 that are unresolved by Gaia; even at the smallest known image separation of 0.18 arcsec, they show an AEN > 1.0 mas.

For the lensing calculation, we followed a procedure used in the analysis of a recently discovered high-redshift lensed quasar<sup>32</sup>, and adjusted the source redshift to the $z \approx 4$ of J0529-4351.

### Spectroscopic observations and data processing

The optical-to-near-infrared spectrum of J0529-4351 was obtained with the X-Shooter instrument<sup>86</sup> on ESO's Very Large Telescope (UT3) on 14 January 2023. The air mass at the time of observation was 1.06 and the estimated seeing was about 1.5 arcsec. In the UVB, visible and near-infrared arms, the slit widths and total exposure times were: 1.0 arcsec and 500 s; 0.9 arcsec and 400 s; and 0.9 arcsec and 800 s, respectively.

The data were reduced using Pypelt, a Python-based spectroscopic data reduction pipeline<sup>87</sup>. Because it was not possible to a priori know the underlying true shape of the C ivλ1549 emission line beneath the atmospheric A-band absorption, we could not reliably fit the telluric model to the J0529-4351 visible spectrum alone. To obtain an improved correction, we took the telluric model generated from a standard star observed on the same night and modified it by adopting a power-law intrinsic stellar spectrum across a broad span of wavelengths unaffected by strong stellar atmospheric lines. Applying this modified model to J0529-4351 yielded a final spectrum with a smooth C ivλ1549 profile containing only narrow absorption lines and no spurious emission from overcorrected telluric absorption.

The observed spectrum was corrected for slit losses by calibration to external photometric data, in particular, the quasar's J-band magnitude of 14.812 ± 0.003 mag (Vega) in the VISTA Hemisphere Survey (VHS) DR6. With the VHS calibration anchoring the near-infrared arm, the visible and UVB data were sequentially matched to the spectroscopic flux levels in the wavelength regions of overlap between the arms. This J-band calibration is consistent (within the photometric errors) with that of the i<sub>SDSS</sub> X-Shooter acquisition image, as calibrated by synthetic photometry for a neighbouring, non-variable star with Gaia low-resolution spectroscopy and the flux scale measured in the Gaia spectrum of J0529-4351 itself.

### Spectral decomposition and emission-line fitting

The reduced and telluric-corrected X-Shooter spectrum was then transformed into the rest frame using a redshift of $z = 3.962$. We modelled the
broad emission-line profile of both the C IVλ1549 line and Mg IIλ2799 line using the code PyQSpecFit³⁸, a Python-based spectral modelling package that is designed for rest-frame UV and optical quasar spectra. We modelled each line individually, selecting independent windows in wavelength to constrain the quasar continuum and the emission-line flux. We masked the narrow absorption features present in our spectrum by applying a boxcar sigma-clipping routine with a width of 50 pixels and a 3σ threshold.

For C IVλ1549, we constrained the continuum on either side of the feature with a power-law model over line-free wavelength windows of 1,445 Å–1,455 Å and 1,973 Å–1,983 Å. The emission-line model was composed of a maximum of three broad Gaussian components with a minimum FWHM of 1,000 km s⁻¹. We fitted the C IVλ1549 profile between 1,480 Å and 1,528 Å and between 1,537 Å and 1,565 Å, avoiding the narrow absorption feature between 1,528 Å and 1,537 Å. We also measured the monochromatic luminosity at 1,450 Å using the power-law model of the continuum. We estimated the errors in the spectral fitting by creating 50 realizations of the spectrum after randomly redistributing the flux data according to each point's Gaussian uncertainties. The standard deviation in each of the derived properties of the C IVλ1549 emission line was adopted as the statistical error.

For Mg IIλ2799, the continuum model was composed of a power law and a template of the UV flux contribution from broad Fe II emission, which was constrained across the wavelength ranges 1,973 Å–1,983 Å, 2,060 Å–2,340 Å, 2,600 Å–2,740 Å and 2,840 Å–3,100 Å. The emission-line model for Mg IIλ2799 was composed of one narrow and up to three broad Gaussian components, where the narrow and broad distinction was set at a FWHM of 1,000 km s⁻¹. We fitted the Mg IIλ2799 line between 2,700 Å and 2,870 Å and the monochromatic luminosity at 3,000 Å was also measured from the continuum power-law model. We derived the uncertainty in each of the line properties by adopting four different templates of the UV Fe II emission⁸⁸⁻⁹¹. As with C IVλ1549, we resampled the spectrum 50 times and measured the standard deviation in the Mg IIλ2799 line properties, summing in quadrature with the uncertainty derived from Fe II templates. More information about the PyQSpecFit modelling of C IVλ1549 and Mg IIλ2799 lines is detailed in a study of the now second-most luminous quasar, SMSS J2157–3602 (ref. 14), and also in a study of line properties in the luminous quasar sample of the European Southern Observatory Large Programme XQ-100 (ref. 59), which is at a comparable redshift range to J0529–4351.

## Luminosity calculations

For calculating monochromatic continuum luminosities, we corrected the spectra for extinction by dust in the Milky Way. We used the estimate of $E(B-V) = 0.041 \pm 0.0027$ from the Schlegel maps⁹² with the correction factor⁹³ of 0.86 and the extinction law by Fitzpatrick⁹⁴, which yields rest-frame absorption values of $A_{1,450} = 0.074$ mag and $A_{3,000} = 0.025$ mag. We then applied bolometric corrections for mean quasar SEDs³⁹ of $k_{1,350} = 3.81$ and $k_{3,000} = 5.15$ and applied an anisotropy correction factor of 0.75 assuming mean orientation⁴⁰. The accretion disk continuum fitting described below determines the bolometric luminosity directly from the integrated disk-model SEDs and implies a spin-and inclination-dependent anisotropy correction. We then took the average of the luminosity measurements.

## Spectral decomposition and continuum fitting

We used a publicly available code⁴⁸ to model the shape of the accretion disk continuum to SED models predicted for slim disks with slimbh⁴⁷, which is a grid of synthetic spectra from ray-traced numerical solutions of slim accretion disk equations. The free parameters are mass and spin of the black hole as well as accretion rate and disk inclination. A Markov chain Monte Carlo method was used to map out likelihood contours across the four-dimensional parameter space. We collected the Milky Way extinction-corrected spectrum and extended the infrared coverage with WISE photometry, specifically the CatWISE2020 (ref. 95) W1

and W2 photometric passbands. We then created synthetic data points to represent the accretion-disk flux within selected line-free windows in the observed spectrum. For the WISE wavelengths, we adopted the spectrum of the Selsing high-luminosity quasar template⁹⁶, scaled independently to the flux of W1 and W2 bandpasses, and created synthetic line-free data points in the same fashion. Furthermore, we used the X-Shooter UVB arm to estimate the continuum flux through the Lyman-α forest and create one more synthetic data point to constrain the continuum. We then used a Markov chain Monte Carlo method to infer the Bayesian posterior probability distributions of the intrinsic black-hole properties using the set of synthetic data points to represent the flux of the accretion disk continuum.

Further detailed information on the continuum-fitting method is presented in an in-depth study of the quasar SMSS J2157–3602 (ref. 14). However, unlike this previous study, the static and smooth thermal accretion disk models are unable to fully reproduce the hardening of the continuum in the visible arm or the flux passing through the Lyman-α forest. It has been shown previously that quasars with high C IVλ1549 blueshift, which is indicative of strong outflows and winds, tend towards a bluer UV continuum to the blue side of $\lambda_{rest} \approx 200$ nm and a slightly redder continuum on the red side⁹⁷, which thus produces a spectral break that is missing in the thermal-disk models. Therefore, our slimbh synthetic spectra utilize the disk atmosphere model, BHSPEC⁹⁸,⁹⁹, to help to reproduce the emerging Compton-hardened radiation.

We also attempted to fit the data with alternative models. For example, we adopted kerubb thin-disk models¹⁰⁰ despite the fact that the slim-disk model reproduces the thin-disk SED at low Eddington ratios, and indeed found similar results for the mass of the black hole and the luminosity of the disk, although the colour of the UV continuum was not properly reproduced in the absence of Comptonization considerations. We also examined the option of J0529–4351 being lensed by demagnifying the spectrum and found that it did not improve the quality of fit.

## Host galaxy dust extinction

We assumed no dust extinction in the quasar host galaxy given that typical host reddening levels are found to be consistent with $E(B-V) \approx 0.0$ among luminous quasars¹⁰¹, whereas less than 1% of quasars seem to have $E(B-V) > 0.1$. Also, an object appearing as the brightest object in the Universe has a low probability of being extinguished by notable levels of dust. If dust were present, it would make the continuum redder and the measured luminosity lower. This would lead to the black-hole mass being overestimated by the continuum-fitting method and underestimated by the virial single-epoch method. Such a discrepancy was not observed here.

## Light-curve construction

The light curve with photometry from the NASA ATLAS³³ 0.5 m telescope was obtained from the ATLAS website. For any object in the ATLAS observing area, up to four observations per night are available, depending on weather. The orange passband is observed in all clear nights, while the cyan passband is only used during the half-period around a new moon. For slowly varying objects, noise can be suppressed by combining observations from longer periods; we determined median measurements to reduce the influence of outlier measurements, using per-week intervals in the orange passband and per-moon-period intervals in the cyan passband. Error bars on the medians express the level of variability within the interval by showing interquartile ranges of the values.

## Data availability

Data from Data Release 3 (DR3) of the European Space Agency's Gaia mission are publicly available (https://gea.esac.esa.int/archive/). NASA ATLAS data are available from https://fallingstar-data.com/forced-phot/. The SkyMapper Southern Survey data are available from https://
skymapper.anu.edu.au/(https://doi.org/10.25914/5f14eded2d116). The raw spectrum and calibration files from the ESO/VLT are available in the ESO archive at http://archive.eso.org/. A reduced spectrum is available from the authors on reasonable request.

## Code availability

The spectral fitting code and the quasar continuum-fitting code were written by S.L. in Python and are publicly available on GitHub at https://github.com/samlaihei/PyQSpecFit (ref. 38) and https://github.com/samlaihei/BADFit (ref. 48).

## References

1. Schmidt, M. 3C 273: a star-like object with large red-shift. *Nature* **197**, 1040 (1963).

2. Hoyle, F. & Fowler, W. A. Nature of strong radio sources. *Nature* **197**, 533–535 (1963).

3. Shakura, N. I. & Sunyaev, R. A. Black holes in binary systems. Observational appearance. *Astron. Astrophys.* **24**, 337–355 (1973).

4. Soltan, A. Masses of quasars. *Mon. Not. R. Astron. Soc.* **200**, 115–122 (1982).

5. Boyle, B. J. et al. The 2dF QSO Redshift Survey – I. The optical luminosity function of quasi-stellar objects. *Mon. Not. R. Astron. Soc.* **317**, 1014–1022 (2000).

6. Richards, G. T. et al. The Sloan Digital Sky Survey quasar survey: quasar luminosity function from Data Release 3. *Astron. J.* **131**, 2766–2787 (2006).

7. Tucci, M. & Volonteri, M. Constraining supermassive black hole evolution through the continuity equation. *Astron. Astrophys.* **600**, A64 (2017).

8. Bardeen, J. M. Kerr metric black holes. *Nature* **226**, 64–65 (1970).

9. Davis, S. W. & Laor, A. The radiative efficiency of accretion flows in individual active galactic nuclei. *Astrophys. J.* **728**, 98 (2011).

10. Dexter, J. & Agol, E. Quasar accretion disks are strongly inhomogeneous. *Astrophys. J. Lett.* **727**, L24 (2011).

11. Slone, O. & Netzer, H. The effects of disc winds on the spectrum and black hole growth rate of active galactic nuclei. *Mon. Not. R. Astron. Soc.* **426**, 656–664 (2012).

12. Campitiello, S., Ghisellini, G., Sbarrato, T. & Calderone, G. How to constrain mass and spin of supermassive black holes through their disk emission. *Astron. Astrophys.* **612**, A59 (2018).

13. Starkey, D. A., Huang, J., Horne, K. & Lin, D. N. C. Rimmed and rippled accretion disc models to explain AGN continuum lags. *Mon. Not. R. Astron. Soc.* **519**, 2754–2768 (2023).

14. Lai, S., Wolf, C., Onken, C. A. & Bian, F. Characterising SMSS J2157–3602, the most luminous known quasar, with accretion disc models. *Mon. Not. R. Astron. Soc.* **521**, 3682–3698 (2023).

15. Flesch, E. W. Identification confusion and blending concealment in the SDSS-DR16 quasar catalogues – 40 new quasars and 82 false quasars identified. *Mon. Not. R. Astron. Soc.* **504**, 621–635 (2021).

16. Wu, X.-B. et al. An ultraluminous quasar with a twelve-billion-solar-mass black hole at redshift 6.30. *Nature* **518**, 512–515 (2015).

17. Eilers, A.-C. et al. EIGER. III. JWST/NIRCam observations of the ultraluminous high-redshift quasar J0100+2802. *Astrophys. J.* **950**, 68 (2023).

18. Wolf, C. et al. Discovery of the most ultra-luminous QSO using Gaia, SkyMapper, and WISE. *Publ. Astron. Soc. Aust.* **35**, e024 (2018).

19. Volonteri, M. Formation of supermassive black holes. *Astron. Astrophys. Rev.* **18**, 279–315 (2010).

20. Volonteri, M., Silk, J. & Dubus, G. The case for supercritical accretion onto massive black holes at high redshift. *Astrophys. J.* **804**, 148 (2015).

21. Amarantidis, S. et al. The first supermassive black holes: indications from models for future observations. *Mon. Not. R. Astron. Soc.* **485**, 2694–2709 (2019).

22. Zubovas, K. & King, A. High-redshift SMBHs can grow from stellar-mass seeds via chaotic accretion. *Mon. Not. R. Astron. Soc.* **501**, 4289–4297 (2021).

23. Onken, C. A. et al. AllBRICQS: the All-sky BRIght, Complete Quasar Survey. *Publ. Astron. Soc. Aust.* **40**, e010 (2023).

24. Walsh, D., Carswell, R. F. & Weymann, R. J. 0957+561 A, B: twin quasistellar objects or gravitational lens? *Nature* **279**, 381–384 (1979).

25. Lemon, C. et al. Gravitationally lensed quasars in Gaia – IV. 150 new lenses, quasar pairs, and projected quasars. *Mon. Not. R. Astron. Soc.* **520**, 3305–3328 (2023).

26. Irwin, M. J., Ibata, R. A., Lewis, G. F. & Totten, E. J. APM 08279+5255: an ultraluminous broad absorption line quasar at a redshift z=3.87. *Astrophys. J.* **505**, 529–535 (1998).

27. Patnaik, A. R., Browne, I. W. A., Walsh, D., Chaffee, F. H. & Foltz, C. B. B 1422+231: a new gravitationally lensed system at z=3.62. *Mon. Not. R. Astron. Soc.* **259**, 1P–4 (1992).

28. Egami, E. et al. APM 08279+5255: Keck near- and mid-infrared high-resolution imaging. *Astrophys. J.* **535**, 561–574 (2000).

29. Chen, Y.-C. et al. Varstrometry for off-nucleus and dual subkiloparsec AGN (VODKA): Hubble Space Telescope discovers double quasars. *Astrophys. J.* **925**, 162 (2022).

30. Chen, H.-W. et al. An empirical characterization of extended cool gas around galaxies using Mg II absorption features. *Astrophys. J.* **714**, 1521–1541 (2010).

31. Mason, C. A. et al. Correcting the z~8 galaxy luminosity function for gravitational lensing magnification bias. *Astrophys. J.* **805**, 79 (2015).

32. Yue, M., Fan, X., Yang, J. & Wang, F. Revisiting the lensed fraction of high-redshift quasars. *Astrophys. J.* **925**, 169 (2022).

33. Tonry, J. L. et al. Atlas: a high-cadence all-sky survey system. *Publ. Astron. Soc. Pac.* **130**, 064505 (2018).

34. Onken, C. A. et al. SkyMapper Southern Survey: second data release (DR2). *Publ. Astron. Soc. Aust.* **36**, e033 (2019).

35. Wright, E. L. et al. The Wide-field Infrared Survey Explorer (WISE): mission description and initial on-orbit performance. *Astron. J.* **140**, 1868–1881 (2010).

36. McConnell, D. et al. The Rapid ASKAP Continuum Survey I: design and first results. *Publ. Astron. Soc. Aust.* **37**, e048 (2020).

37. Kellermann, K. I., Sramek, R., Schmidt, M., Shaffer, D. B. & Green, R. VLA observations of objects in the Palomar Bright Quasar Survey. *Astron. J.* **98**, 1195–1207 (1989).

38. Lai, S. samlaihei/pyqspecfit: Pyqspecfit v1.0.0. Zenodo https://doi.org/10.5281/zenodo.7772752 (2023).

39. Richards, G. T. et al. Spectral energy distributions and multiwavelength selection of type 1 quasars. *Astrophys. J. Suppl. Ser.* **166**, 470–497 (2006).

40. Runnoe, J. C., Brotherton, M. S. & Shang, Z. Updating quasar bolometric luminosity corrections. *Mon. Not. R. Astron. Soc.* **422**, 478–493 (2012).

41. Yu, Q. & Tremaine, S. Observational constraints on growth of massive black holes. *Mon. Not. R. Astron. Soc.* **335**, 965–976 (2002).

42. Calderone, G., Ghisellini, G., Colpi, M. & Dotti, M. Black hole mass estimate for a sample of radio-loud narrow-line Seyfert 1 galaxies. *Mon. Not. R. Astron. Soc.* **431**, 210–239 (2013).

43. Capellupo, D. M., Netzer, H., Lira, P., Trakhtenbrot, B. & Mejía-Restrepo, J. Active galactic nuclei at z~1.5–1. Spectral energy distribution and accretion discs. *Mon. Not. R. Astron. Soc.* **446**, 3427–3446 (2015).

44. Mejía-Restrepo, J. E., Lira, P., Netzer, H., Trakhtenbrot, B. & Capellupo, D. M. The effect of nuclear gas distribution on the mass determination of supermassive black holes. *Nat. Astron.* **2**, 63–68 (2018).
45. Campitiello, S., Celotti, A., Ghisellini, G. & Sbarrato, T. Estimating black hole masses: accretion disk fitting versus reverberation mapping and single epoch. *Astron. Astrophys.* **640**, A39 (2020).

46. Abramowicz, M. A., Czerny, B., Lasota, J. P. & Szuszkiewicz, E. Slim accretion disks. *Astrophys. J.* **332**, 646–658 (1988).

47. Sadowski, A. et al. Relativistic slim disks with vertical structure. *Astron. Astrophys.* **527**, A17 (2011).

48. Lai, S. samlaihei/badfit: Badfit v1.0.0. Zenodo https://doi.org/10.5281/zenodo.7772748 (2023).

49. Laor, A. & Davis, S. W. Cold accretion discs and lineless quasars. *Mon. Not. R. Astron. Soc.* **417**, 681–688 (2011).

50. Malkan, M. A. The ultraviolet excess of luminous quasars. II. Evidence for massive accretion disks. *Astrophys. J.* **268**, 582–590 (1983).

51. GRAVITY Collaboration. Spatially resolved rotation of the broad-line region of a quasar at sub-parsec scale. *Nature* **563**, 657–660 (2018).

52. McLure, R. J. & Jarvis, M. J. Measuring the black hole masses of high-redshift quasars. *Mon. Not. R. Astron. Soc.* **337**, 109–116 (2002).

53. Vestergaard, M. & Peterson, B. M. Determining central black hole masses in distant active galaxies and quasars. II. Improved optical and UV scaling relationships. *Astrophys. J.* **641**, 689–709 (2006).

54. Peterson, B. M. et al. Central masses and broad-line region sizes of active galactic nuclei. II. A homogeneous analysis of a large reverberation-mapping database. *Astrophys. J.* **613**, 682–699 (2004).

55. Coatman, L. et al. Correcting C IV-based virial black hole masses. *Mon. Not. R. Astron. Soc.* **465**, 2120–2142 (2017).

56. Vestergaard, M. & Osmer, P. S. Mass functions of the active black holes in distant quasars from the Large Bright Quasar Survey, the Bright Quasar Survey, and the color-selected sample of the SDSS Fall Equatorial Stripe. *Astrophys. J.* **699**, 800–816 (2009).

57. Shen, Y. et al. A catalog of quasar properties from Sloan Digital Sky Survey Data Release 7. *Astrophys. J. Suppl. Ser.* **194**, 45 (2011).

58. Shen, Y., Greene, J. E., Strauss, M. A., Richards, G. T. & Schneider, D. P. Biases in virial black hole masses: an SDSS perspective. *Astrophys. J.* **680**, 169–190 (2008).

59. Lai, S. et al. Virial black hole mass estimates of quasars in the XQ-100 legacy survey. *Mon. Not. R. Astron. Soc.* **526**, 3230–3247 (2023).

60. Tsukui, T. & Iguchi, S. Spiral morphology in an intensely star-forming disk galaxy more than 12 billion years ago. *Science* **372**, 1201–1205 (2021).

61. Bañados, E. et al. The z=7.54 quasar ULAS J1342+0928 is hosted by a galaxy merger. *Astrophys. J. Lett.* **881**, L23 (2019).

62. Bischetti, M. et al. The WISSH quasars project. IX. Cold gas content and environment of luminous QSOs at z~2.4-4.7. *Astron. Astrophys.* **645**, A33 (2021).

63. Onken, C. A. et al. Discovery of the most luminous quasar of the last 9 Gyr. *Publ. Astron. Soc. Aust.* **39**, e037 (2022).

64. Webster, R. L., Francis, P. J., Petersont, B. A., Drinkwater, M. J. & Masci, F. J. Evidence for a large undetected population of dust-reddened quasars. *Nature* **375**, 469–471 (1995).

65. Warren, S. J., Hewett, P. C. & Foltz, C. B. The KX method for producing K-band flux-limited samples of quasars. *Mon. Not. R. Astron. Soc.* **312**, 827–832 (2000).

66. Wisotzki, L. et al. The Hamburg/ESO survey for bright QSOs. III. A large flux-limited sample of QSOs. *Astron. Astrophys.* **358**, 77–87 (2000).

67. Fan, X. et al. High-redshift quasars found in Sloan Digital Sky Survey commissioning data. III. A color-selected sample at *I*\*<20 in the Fall Equatorial Stripe. *Astron. J.* **121**, 31–53 (2001).

68. Wolf, C. et al. The evolution of faint AGN between z=1 and z=5 from the COMBO-17 survey. *Astron. Astrophys.* **408**, 499–514 (2003).

69. Stern, D. et al. Mid-infrared selection of active galaxies. *Astrophys. J.* **631**, 163–168 (2005).

70. Richards, G. T. et al. Efficient photometric selection of quasars from the Sloan Digital Sky Survey. II. ~1,000,000 quasars from Data Release 6. *Astrophys. J. Suppl. Ser.* **180**, 67–83 (2009).

71. Ivezić, Ž. et al. in *Multiwavelength AGN Surveys and Studies* Vol. 304 (eds Mickaelian, A. M. & Sanders, D. B.) 11–17 (IAUS, 2014).

72. Reed, S. L. et al. Eight new luminous z≥6 quasars discovered via SED model fitting of VISTA, WISE and Dark Energy Survey Year 1 observations. *Mon. Not. R. Astron. Soc.* **468**, 4702–4718 (2017).

73. Calderone, G. et al. Finding the brightest cosmic beacons in the southern hemisphere. *Astrophys. J.* **887**, 268 (2019).

74. Yang, Q. & Shen, Y. A southern photometric quasar catalog from the Dark Energy Survey Data Release 2. *Astrophys. J. Suppl. Ser.* **264**, 9 (2023).

75. Storey-Fisher, K. et al. Quaia, the Gaia-unWISE quasar catalog: an all-sky spectroscopic quasar sample. Preprint at https://doi.org/10.48550/arXiv.2306.17749 (2023).

76. Delchambre, L. et al. Gaia DR3: Apsis III – non-stellar content and source classification. *Astron. Astrophys.* **674**, A31 (2023).

77. Fu, Y. et al. Finding quasars behind the galactic plane. II. Spectroscopic identifications of 204 quasars at |b|<20°. *Astrophys. J. Suppl. Ser.* **261**, 32 (2022).

78. Bentz, M. C. et al. The low-luminosity end of the radius-luminosity relationship for active galactic nuclei. *Astrophys. J.* **767**, 149 (2013).

79. Sandage, A. The change of redshift and apparent luminosity of galaxies due to the deceleration of selected expanding universes. *Astrophys. J.* **136**, 319–333 (1962).

80. Liske, J. et al. Cosmic dynamics in the era of extremely large telescopes. *Mon. Not. R. Astron. Soc.* **386**, 1192–1218 (2008).

81. Cristiani, S. et al. Spectroscopy of QUBRICS quasar candidates: 1672 new redshifts and a golden sample for the Sandage test of the redshift drift. *Mon. Not. R. Astron. Soc.* **522**, 2019–2028 (2023).

82. Gaia Collaboration et al. Gaia Early Data Release 3. Summary of the contents and survey properties. *Astron. Astrophys.* **649**, A1 (2021).

83. Dey, A. et al. Overview of the DESI Legacy Imaging Surveys. *Astron. J.* **157**, 168 (2019).

84. Fabricius, C. et al. Gaia Early Data Release 3. Catalogue validation. *Astron. Astrophys.* **649**, A5 (2021).

85. Lemon, C. Gravitationally Lensed Quasar Database (University of Cambridge); https://research.ast.cam.ac.uk/lensedquasars/

86. Vernet, J. et al. X-shooter, the new wide band intermediate resolution spectrograph at the ESO Very Large Telescope. *Astron. Astrophys.* **536**, A105 (2011).

87. Prochaska, J. X. et al. Pypeit: the Python spectroscopic data reduction pipeline. *J. Open Source Softw.* **5**, 2308 (2020).

88. Vestergaard, M. & Wilkes, B. J. An empirical ultraviolet template for iron emission in quasars as derived from I Zwicky 1. *Astrophys. J. Suppl. Ser.* **134**, 1–33 (2001).

89. Tsuzuki, Y. et al. Fe II emission in 14 low-redshift quasars. I. Observations. *Astrophys. J.* **650**, 57–79 (2006).

90. Mejía-Restrepo, J. E., Trakhtenbrot, B., Lira, P., Netzer, H. & Capellupo, D. M. Active galactic nuclei at z~1.5 – II. Black hole mass estimation by means of broad emission lines. *Mon. Not. R. Astron. Soc.* **460**, 187–211 (2016).

91. Bruhweiler, F. & Verner, E. Modeling Fe II emission and revised Fe II (UV) empirical templates for the Seyfert 1 galaxy I Zw 1. *Astrophys. J.* **675**, 83–95 (2008).

92. Schlegel, D. J., Finkbeiner, D. P. & Davis, M. Maps of dust infrared emission for use in estimation of reddening and cosmic microwave background radiation foregrounds. *Astrophys. J.* **500**, 525–553 (1998).
93. Schlafly, E. F. & Finkbeiner, D. P. Measuring reddening with Sloan Digital Sky Survey stellar spectra and recalibrating SFD. *Astrophys. J.* **737**, 103 (2011).

94. Fitzpatrick, E. L. Correcting for the effects of interstellar extinction. *Publ. Astron. Soc. Pac.* **11**, 63–75 (1999).

95. Marocco, F. et al. The CatWISE2020 catalog. *Astrophys. J. Suppl.* Ser. **253**, 8 (2021).

96. Selsing, J., Fynbo, J. P. U., Christensen, L. & Krogager, J. K. An X-Shooter composite of bright $1<z<2$ quasars from UV to infrared. *Astron. Astrophys.* **585**, A87 (2016).

97. Temple, M. J., Banerji, M., Hewett, P. C., Rankine, A. L. & Richards, G. T. Exploring the link between C IV outflow kinematics and sublimation-temperature dust in quasars. *Mon. Not. R. Astron. Soc.* **501**, 3061–3073 (2021).

98. Davis, S. W., Blaes, O. M., Hubeny, I. & Turner, N. J. Relativistic accretion disk models of high-state black hole X-ray binary spectra. *Astrophys. J.* **621**, 372–387 (2005).

99. Davis, S. W. & Hubeny, I. A grid of relativistic, non-LTE accretion disk models for spectral fitting of black hole binaries. *Astrophys. J. Suppl.* Ser. **164**, 530–535 (2006).

100. Li, L.-X., Zimmerman, E. R., Narayan, R. & McClintock, J. E. Multitemperature blackbody spectrum of a thin accretion disk around a Kerr black hole: model computations and comparison with observations. *Astrophys. J. Suppl.* Ser. **157**, 335–370 (2005).

101. Krawczyk, C. M. et al. Mining for dust in type 1 quasars. *Astron. J.* **149**, 203 (2015).

102. Astropy Collaboration et al. The Astropy project: sustaining and growing a community-oriented open-source project and the latest major release (v5.0) of the core package. *Astrophys. J.* **935**, 167 (2022).

## Acknowledgements

This work was supported by the Australian Research Council (ARC) through Discovery Project DP190100252 (C.W., F.B., C.A.O. and S.L.). S.L. is grateful to the Research School of Astronomy and Astrophysics at the Australian National University (ANU) for funding his Ph.D. studentship. We thank G. Ferrami from the University of Melbourne for discussing solutions for strong gravitational lensing. Data for this project were obtained at the European Southern Observatory through DDT proposal 2110.B-5032. This work has made use of data from the European Space Agency mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular, the institutions participating in the Gaia Multilateral Agreement. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, and NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology. WISE and NEOWISE are funded by NASA. SuperCOSMOS Sky Survey material is based on photographic data originating from the UK, Palomar and ESO Schmidt telescopes and is provided by the Wide-Field Astronomy Unit, Institute for Astronomy, University of Edinburgh. This work made use of Astropy (http://www.astropy.org), a community-developed core Python package and an ecosystem

of tools and resources for astronomy¹⁰². This research has made use of the SVO Filter Profile Service (http://svo2.cab.inta-csic.es/theory/fps/) supported from the Spanish MINECO through grant AYA2017-84089. The national facility capability for SkyMapper has been funded through ARC LIEF grant LE130100104 from the Australian Research Council, awarded to the University of Sydney, the ANU, Swinburne University of Technology, the University of Queensland, the University of Western Australia, the University of Melbourne, Curtin University of Technology, Monash University and the Australian Astronomical Observatory. SkyMapper is owned and operated by The ANU's Research School of Astronomy and Astrophysics. The survey data were processed and provided by the SkyMapper Team at ANU. The SkyMapper node of the All-Sky Virtual Observatory is hosted at the National Computational Infrastructure. Development and support of the SkyMapper node of the All-Sky Virtual Observatory has been funded in part by Astronomy Australia Limited and the Australian Government through the Commonwealth's Education Investment Fund and the National Collaborative Research Infrastructure Strategy, in particular, the National eResearch Collaboration Tools and Resources and the Australian National Data Service Projects. This work uses data from the University of Hawaii's ATLAS project, funded through NASA grants NN12AR55G, 80NSSC18K0284 and 80NSSC18K1575, with contributions from the Queen's University Belfast, the Space Telescope Science Institute, the South African Astronomical Observatory and the Millennium Institute of Astrophysics, Chile.

## Author contributions

All authors contributed to data collection. S.L. led the data analysis with contributions from C.A.O. and C.W. C.W. selected the quasar candidates and led the drafting and editing of the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

**Correspondence and requests for materials** should be addressed to Christian Wolf.

**Peer review information** *Nature Astronomy* thanks the anonymous reviewers for their contribution to the peer review of this work.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

© The Author(s), under exclusive licence to Springer Nature Limited 2024