# EIGER III. JWST/NIRCam observations of the ultra-luminous high-redshift
quasar J0100+2802

Anna-Christina Eilers MIT Kavli Institute for Astrophysics and Space Research,
77 Massachusetts Avenue, Cambridge, 02139, Massachusetts, USA  Robert A.
Simcoe MIT Kavli Institute for Astrophysics and Space Research, 77
Massachusetts Avenue, Cambridge, 02139, Massachusetts, USA  Minghao Yue MIT
Kavli Institute for Astrophysics and Space Research, 77 Massachusetts Avenue,
Cambridge, 02139, Massachusetts, USA  Ruari Mackenzie Department of Physics,
ETH Zürich, Wolfgang-Pauli-Strasse 27, Zürich, 8093, Switzerland  Jorryt
Matthee Department of Physics, ETH Zürich, Wolfgang-Pauli-Strasse 27, Zürich,
8093, Switzerland  Dominika Ďurovčíková MIT Kavli Institute for Astrophysics
and Space Research, 77 Massachusetts Avenue, Cambridge, 02139, Massachusetts,
USA  Daichi Kashino Institute for Advanced Research, Nagoya University, Nagoya
464-8601, Japan  Department of Physics, Graduate School of Science, Nagoya
University, Nagoya 464-8602, Japan  Rongmon Bordoloi Department of Physics,
North Carolina State University, Raleigh, 27695, North Carolina, USA  Simon J.
Lilly Department of Physics, ETH Zürich, Wolfgang-Pauli-Strasse 27, Zürich,
8093, Switzerland Anna-Christina Eilers<EMAIL_ADDRESS>
###### Abstract

We present the first rest-frame optical spectrum of a high-redshift quasar
observed with JWST/NIRCam in Wide Field Slitless (WFSS) mode. The observed
quasar, J0100+2802, is the most luminous quasar known at $z>6$. We measure the
mass of the central supermassive black hole (SMBH) by means of the rest-frame
optical H$\beta\,$ emission line, and find consistent mass measurements of the
quasar’s SMBH of $M_{\bullet}\approx 10^{10}\,M_{\sun}$ when compared to the
estimates based on the properties of rest-frame UV emission lines C IV and Mg
II, which are accessible from ground-based observatories. To this end, we also
present a newly reduced rest-frame UV spectrum of the quasar observed with
X-Shooter/VLT and FIRE/Magellan for a total of $16.8$ hours. We readdress the
question whether this ultra-luminous quasar could be effected by strong
gravitational lensing making use of the diffraction limited NIRCam images in
three different wide band filters (F115W, F200W, F356W), which improves the
achieved spatial resolution compared to previous images taken with the Hubble
Space Telescope by a factor of two. We do not find any evidence for a
foreground deflecting galaxy, nor for multiple images of the quasar, and
determine the probability for magnification due to strong gravitational
lensing with image separations below the diffraction limit of
$\Delta\theta\lesssim 0.05\arcsec$ to be $\lesssim 2.2\times 10^{-3}$. Our
observations therefore confirm that this quasar hosts a ten billion solar mass
black hole less than $1$ Gyr after the Big Bang, which is challenging to
explain with current black hole formation models.

dark ages, early universe — quasars: emission lines, supermassive black holes
— methods: data analysis — gravitational lensing: strong gravitational lensing

††software: numpy (Harris et al., 2020), scipy (Jones et al., 2001),
matplotlib (Hunter, 2007), astropy (Astropy Collaboration et al., 2018),
PypeIt (Prochaska et al., 2020), photutils (Bradley et al., 2020), psfMC
(Mechtley, 2014)††thanks: Pappalardo Fellow

##  1 Introduction

Since the discovery of the first luminous quasars nearly six decades ago
(Schmidt, 1965) astronomers have detected and observed more than a million
quasars in the universe up to redshifts of $z\gtrsim 7.5$ (e.g. Bañados et
al., 2016; Lyke et al., 2020; Wang et al., 2021). Quasars are the most
luminous, non-transient objects in the universe, powered by accretion onto a
central supermassive black hole (SMBH) of millions to several billions of
solar masses in size (e.g. Mazzucchelli et al., 2017; Yang et al., 2021; Wu &
Shen, 2022). The most luminous, high-redshift quasar ($z>6$) known to date was
discovered by Wu et al. (2015), i.e. J0100+2802 at a redshift of $z=6.3270\pm
0.0005$ (Wang et al., 2019). Its extreme luminosity of $L_{\rm bol}\sim
10^{48}\rm\,erg\,s^{-1}$ implies that the quasar hosts a highly accreting SMBH
in its center with a mass of $M_{\bullet}\sim 10^{10}\,M_{\sun}$ (Wu et al.,
2015).

Thus, J0100+2802 harbors the most massive SMBH known at $z\gtrsim 6$ – which
is among the $0.2\%$ of the most massive SMBHs at all redshifts (Wu & Shen,
2022) – at a time when the universe is only $\sim 800$ Myr old. Its existence
poses significant challenges to current models aiming to explain the rapid
growth of SMBHs in the early universe (e.g. Bañados et al., 2018). Adding to
the challenge, the quasar’s rest-frame ultraviolet (UV) spectrum exhibits a
very small proximity zone in the vicinity of the quasar, which indicates a
short UV luminous quasar lifetime of $t_{\rm Q}\sim 10^{5}$ yr, during which
accretion onto the black hole is expected to occur (e.g. Eilers et al., 2017;
Davies et al., 2020; Morey et al., 2021).

To explain the anomaly of J0100+2802’s SMBH it has been suspected that this
ultra-luminous quasar could possibly be strongly gravitationally lensed.
Arguments in favor of this hypothesis were largely based on the clumpy
morphology of observations with the Atacama Large Millimetre Array (ALMA;
Fujimoto et al., 2020). However, high-resolution images taken with the Hubble
Space Telescope (HST) did not reveal any evidence for multiple images due to
strong gravitational lensing (Fujimoto et al., 2020, Yue et al. in prep.).
Furthermore, detailed analyses of the flux transmission in the quasar’s
proximity zone (Davies et al., 2020), the expected ratios of X-ray luminosity
to rest-frame UV and IR luminosities (Connor et al., 2021), as well as the
implications for the bright-end slope of the quasar luminosity function
(Pacucci & Loeb, 2020) indicate that a significant magnification of the
quasar’s luminosity is unlikely.

The quasar’s extremely massive SMBH has also called into question the accuracy
of all black hole mass estimates of high-redshift quasars, which are obtained
using scaling relations between the width of broad emission lines observed in
single-epoch quasar spectra, the quasars’ luminosities and their black hole
masses (e.g. Vestergaard & Osmer, 2009). These scaling relations are
calibrated based on a small sample of quasars at very low redshifts, i.e.
$z\lesssim 0.2$, for which black hole masses can be precisely determined by
means of reverberation mapping (RM) measurements (e.g. Vestergaard & Peterson,
2006). However, RM measurements require long monitoring campaigns and thus for
most quasars – especially at higher redshifts where time delays are longer due
to time dilation and their generally more massive black holes – this method is
unfeasible. Thus, masses of SMBHs are inferred from single-epoch spectra using
aforementioned scaling relations, most commonly calibrated to the, at low-
redshift easily observable, rest-frame optical H$\beta\,$ emission line (e.g.
Vestergaard & Peterson, 2006; Grier et al., 2017). Since the H$\beta\,$ line
is very challenging to observe for quasars at $z>4$ with ground-based
observatories, additional uncertainty and possibly biases in the black hole
mass estimates for high-redshift quasars could be introduced by re-scaling the
H$\beta\,$ emission line properties to properties of rest-frame UV emission
lines, such as C IV and Mg II (e.g. Vestergaard & Osmer, 2009; Coatman et al.,
2017). These uncertainties have motivated recent efforts to avoid such scaling
relations altogether and derive black hole mass estimates and other quasar
properties from the spectra themselves by means of data-driven modeling of the
quasars’ spectral features (e.g. Eilers et al., 2022).

In this Letter we present the first spectrum of the ultra-luminous quasar
J0100+2802 at rest-frame optical wavelengths observed with the NIRCam grism on
JWST, as well as a newly reduced high signal-to-noise spectrum at rest-frame
UV wavelengths taken with ground-based facilities (§ 2). This spectrum covers
the wavelength region between $0.6\mu{\rm m}\leq\lambda\leq 4\mu$m (with a gap
between $2.3\mu{\rm m}\leq\lambda\leq 3.1\mu$m), which enables us to obtain
measurements of the quasar’s black hole mass based on the rest-frame UV C IV
and Mg II broad emission lines, as well as the rest-frame optical H$\beta\,$
emission (§ 3). Furthermore, we will show NIRCam Imaging observations that set
constraints on the possibility of strong gravitational lensing of the quasar
(§ 4), before summarizing our results (§ 5). Throughout this work, we adopt a
flat $\Lambda$CDM cosmology with $\Omega_{m}=0.31$ and
$H_{0}=67.7~{}\mathrm{km~{}s^{-1}~{}Mpc^{-1}}$ (Planck Collaboration et al.,
2018).

##  2 Data

###  2.1 Ground-based Spectroscopy

The ground-based optical and NIR spectroscopic data for J0100+2802 (RA:
$01^{\rm h}00^{\rm m}13^{\rm s}.020$; DEC:
$+28^{\circ}02{\arcmin}25{\arcsec}.840$) have been obtained using both the
X-Shooter spectrograph (Vernet et al., 2011) on the Very Large Telescope (VLT)
as well as the Folded-port InfraRed Echellette instrument (FIRE; Simcoe et
al., 2013) on the Magellan Telescope in the years of 2015 and 2016. The target
was observed for a total of $16.8$ hours, of which $11$ hours were observed
with VLT/X-Shooter (program ID: 096.A-0095; PI: Pettini) and $5.8$ hours
observed with Magellan/FIRE (PI: Simcoe).

The ground-based optical and NIR data are reduced consistently with the open-
source python-based spectroscopic data reduction pipeline
PypeIt111https://github.com/pypeit/PypeIt version 1.7.1 (Prochaska et al.,
2020). We derive the wavelength solution from the night sky OH lines, in order
to have on-sky wavelength calibrations and to reduce the overheads of our
observations. The sky subtraction is performed on the 2D images by
differencing exposures dithered along the slit in the standard A-B mode and
fitting a $b$-spline to further eliminate sky line residuals following
Bochanski et al. (2009). We then perform an optimal extraction (Horne, 1986)
to obtain the 1D spectra. The spectra are flux calibrated using sensitivity
functions of standard stars observed during the same observing run. For a few
runs no standard stars were observed and hence a sensitivity function is
created from an A0 star with known magnitude.

We co-add the flux-calibrated 1D spectra from each night and correct for
telluric absorption features by jointly fitting an atmospheric model and a
quasar model. The telluric model grids are produced using the Line-By-Line
Radiative Transfer Model (LBLRTM; Clough et al., 2005), while we fit a third
order polynomial as the quasar model for this very featureless spectrum. The
individual telluric-corrected 1D spectra observed in different observing runs
or with different telescopes and instruments are then co-added weighted by the
average square SNR of the exposure, and stacked on a common wavelength grid.
Note that this steps avoids any interpolation of the data in order to prevent
correlated noise properties.

In order to combine the optical and NIR part of the quasar spectrum observed
with X-Shooter, the two parts are stitched together by matching the median
flux at the intersection between $9900$ Å and $10100$ Å. As a last step we
apply an absolute flux calibration to the fully reduced and co-added quasar
spectrum by requiring their integrated flux within the UKIRT-J filter to match
the observed $J$-band magnitude of the quasar, i.e. $J_{\rm AB}=17.64\pm
0.02$.

Finally, the spectrum is re-binned to a wavelength grid sampled linearly in
velocity space with a chosen step size of $\Delta v\approx 50~{}\rm
km\,s^{-1}$ per pixel. We then divide all native pixels into the respective
wavelength bins, and determine the stacked flux in each bin as the mean from
all native pixels. Note that the final wavelength grid is the weighted average
of the individual wavelengths used for each exposure that fall into a given
wavelength bin in the input wavelength grid, and hence not necessarily
linearly sampled anymore.

The resulting ground-based spectrum shown in the top panel of Fig. 1 has a
${\rm SNR}\approx 205$ per $\Delta v\approx 50~{}\rm km\,s^{-1}$ pixel
measured at a rest-frame wavelength of $\lambda_{\rm rest}=1280\pm 20$ Å.

Figure 1: Quasar spectrum of J0100+2802 observed with ground-based facilities
(top) and NIRCam WFSS on JWST (bottom) with the filter F356W shown as the
orange shaded region. Grey shaded areas show masked wavelengths affected by
foreground absorption systems, regions of significant telluric absorption,
regions with strong sky line residuals, as well as regions at the very edge of
the filter transmission. The red curve shows the best fit to the quasar
spectrum, while the colored curves in the insets show the individual
components of the spectral fit (power-law continuum: blue dashed; Balmer
continuum: dark green dotted; iron-template: green; Gaussian line components:
yellow, grey or purple).

###  2.2 Observations with JWST/NIRCam WFSS

The quasar field of J0100+2802 was observed with JWST/NIRCam in Imaging and
WFSS mode on August 22, 2022, as part of the Emission-line galaxies and
Intergalactic Gas in the Epoch of Reionization (EIGER) GTO program (program
ID: 1243, PI: S. Lilly). The quasar field is observed in a mosaic containing
four individual visits such that the final field of view spans approximately
$3\arcmin\times 6\arcmin$ corresponding to $7.6\,{\rm cMpc}\times 15.3\,{\rm
cMpc}$ at $z\approx 6$. The central $40\arcsec\times 40\arcsec$ area around
the quasar is part of every visit and thus this region has the deepest
observations with a total exposure time of $8760\,$s per visit, i.e.
approximately $9.7$ hours in total on source. The NIRCam WFSS observations of
the J0100+2802 field use the grism “R” in the F356W filter. Simultaneously, we
obtain direct imaging observations in two short wavelengths filters, i.e.
F115W and F200W (see § 2.3). The observational setup as well as the data
reduction process are described in detail Kashino et al. (2022) and Matthee et
al. (2022). We will summarize the main steps briefly below.

The data is reduced using the jwst pipeline222https://jwst-
pipeline.readthedocs.io/en/latest/ version 1.8.0 provided by STScI. We first
run the Detector1Pipeline step to obtain the uncalibrated NIRCam exposures
(*_rate.fits). We then apply two additional reduction steps from the
Spec2Pipeline to transfer the pixel coordinates to astronomical coordinates
(AssignWcsStep) and afterwards apply the flat fielding to the exposures
(FlatFieldStep). We then extract the 2D spectrum of the quasar from each of
the flat-fielded exposures using customized scripts, which are making use of
the grismconf333https://github.com/npirzkal/GRISMCONF module using the latest
(V4) trace models. The spectra are normalized by the filter’s sensitivity
curve and rectified to a common observed wavelength grid, ranging from
$3.0\mu{\rm m}\leq\lambda_{\rm obs}\leq 4.2\mu$m in bins of
$\Delta\lambda=9.75$ Å. We correct the trace of the quasar in each exposure
for small curvatures and then perform an optimal extraction of the quasar
spectrum from each 2D exposure. This results in $96$ individual 1D spectra,
which are co-added using the average square SNR of each exposure as weights to
a combined 1D spectrum. We estimate the uncertainty on each pixel flux via
bootstrapping: to this end we re-sample the $96$ individual spectra with
replacement to create $1000$ co-added 1D spectra from which we take the $16$th
and $84$th percentile at each pixel as the uncertainty. The final spectrum is
shown in the bottom panel of Fig. 1 and has a ${\rm SNR}\approx 235$ per
$9.75\rm{\AA}$ pixel ($\Delta v\sim 80~{}\rm km\,s^{-1}$) measured at a rest-
frame wavelength of $\lambda_{\rm rest}=5100\pm 20$ Å .

The spectrum of the quasar from both the ground-based observatories as well as
JWST/NIRCam WFSS will be made publicly available upon publication.

###  2.3 Imaging Data from JWST/NIRCam

The NIRCam imaging data was reduced using the jwst pipeline version 1.8.2. In
order to avoid saturation in the central pixels of bright sources in the
quasar field we run the Detector1Pipeline without suppressing the computations
for saturated ramps with only one good, unsaturated sample in the ramp-fitting
step. We then run the Image2Pipeline to obtain calibrated images (*_cal.fits),
and apply the $1/f$ noise reduction, remove any snowballs and other artifacts
from the data, and subtract the sky and any wisp artifacts. We use
tweakwcs444https://github.com/spacetelescope/tweakwcs to align all exposures
relative to each other, and then stack all exposure with a global alignment
calibrated to Gaia stars. Any offsets are applied to the calibrated images
directly, before all exposures are ultimately combined to a final mosaic by
the Image3Pipeline. The exposure time of the images in the two short-
wavelengths filters F115W and F200W add up to $4380\,$s per visit, i.e.
approximately $4.9$ hours in total, while the direct images with the F356W
filter sum up to $1578\,$s per visit, i.e. $1.8$ hours total. Individual
exposures are stacked weighted by their inverse variance.

The final images (for one of the four visits) in the three broad band NIRCam
filters F115W, F200W and F356W of a $5\arcsec\times 5\arcsec$ region around
the quasar are shown in the left column of Fig. 3. Note that we do not take
the quasar images from the full mosaic for modeling the PSF (see § 4), but
only the images from one visit (i.e. visit 1), since it results in a better
PSF subtraction.

Since the noise properties in the short wavelength filters that are reported
by the jwst pipeline seem to overestimate the pixel-to-pixel variance in the
science image, we re-scale the inverse variance in the science images by a
factor of $2.0$ in F115W and $1.8$ in F200W. We do not apply any re-scaling in
the long wavelength filter F356W.

##  3 Mass estimates of the quasar’s supermassive black hole

Under the assumption that the dynamics in the quasar’s broad line region (BLR)
are dominated by the gravitational pull of the black hole and the system is in
virial equilibrium, we can estimate the mass of the SMBH using the FWHM of
broad emission lines, which gives an estimate for the velocity of the gas
clouds in the BLR orbiting the black hole, and the quasar’s luminosity, which
acts as a proxy for the radius of the BLR, i.e.

$\frac{M_{\bullet}}{M_{\sun}}=10^{\alpha}\left(\frac{\rm FWHM_{\rm
line}}{1000\,\rm km\,s^{-1}}\right)^{\beta}\left(\frac{\lambda
L_{\lambda}}{10^{44}\,\rm erg\,s^{-1}}\right)^{\gamma}.$ (1)

Traditionally, these scaling relations are calibrated based on a relatively
small sample of low-redshift quasars using the properties of the H$\beta\,$
emission line. However, as H$\beta\,$ is not observable with ground-based
observatories for quasars at $z>4$, these scaling relations have been re-
calibrated using properties of the rest-frame UV broad emission lines, such as
Mg II and C IV instead. The best-fit values for the scaling parameters
$\alpha$, $\beta$ and $\gamma$ dependent on the respective emission lines and
monochromatic luminosity $L_{\lambda}$, and are listed in Tab. 1.

Table 1: Parameters of the black hole mass scaling relations.  emission line |  $\alpha$ | $\beta$ |  $\gamma$ | $\lambda$ | reference
---|---|---|---|---|---
C IV (corrected) | 6.71 | 2 | 0.53 | 1350 Å | Coatman et al. (2017)
Mg II | 6.86 | 2 | 0.50 | 3000 Å | Vestergaard & Osmer (2009)
H$\beta\,$ | 6.91 | 2 | 0.50 | 5100 Å | Vestergaard & Peterson (2006)

###  3.1 Spectral Fitting

To estimate the mass of the SMBH in J0100+2802, we fit the ground-based and
NIRCam/WFSS quasar spectrum in the wavelength regions between $0.95~{}\mu$m
and $2.25~{}\mu$m, as well as $3.15~{}\mu$m and $3.99~{}\mu$m, respectively.
We mask all wavelengths that are affected by foreground metal absorption
systems, regions affected by significant telluric absorption (at
$\lambda\approx 14,000$ Å and $\lambda\approx 19,000$ Å), regions with strong
sky line residuals (at $\lambda\approx 20,000$ Å), as well as regions at the
very edge of the filter transmission (at $\lambda<31,500$ Å and
$\lambda>40,000$ Å). We fit the observed ground- and space-based spectrum
separately, with a combination of (A) a power-law continuum, (B) an iron
template (provided by Vestergaard & Wilkes (2001) for the rest-frame UV iron
emission and by Park et al. (2022) in the rest-frame optical) that is
broadened to match the widths of the broad emission lines through a
convolution of the template with a corresponding Gaussian kernel, (C) a Balmer
continuum below the Balmer edge (which only applies to the ground-based
spectrum; see De Rosa et al., 2014; Schindler et al., 2020, for details), and
(D) one or multiple Gaussian components for the broad emission lines.

We use a single Gaussian component to fit the Si IV $\lambda 1304$, the
unresolved Si II $\lambda\lambda 1393,1402$ doublet, Mg II $\lambda 2799$,
H$\gamma\,$ $\lambda 4341$, [O III] $\lambda 4363$, and He II $\lambda 4868$
emission lines, while we use two Gaussian components to capture the more
complex line profile of the C IV $\lambda\lambda 1548,1550$, and H$\beta\,$
$\lambda 4861$ emission lines, as well as the [O III] $\lambda\lambda
4959,5007$ doublet. This results in $19$ and $32$ free parameters of the
spectral fit for the ground-based and NIRCam/WFSS spectrum, respectively,
which we estimate using the Markov Chain Monte Carlo (MCMC) algorithm emcee
(Foreman-Mackey et al., 2013). We apply flat priors for each parameter and
take the median of the posterior probability distribution as the best
parameter estimate. The final spectral fits are shown in red in Fig. 1.

###  3.2 Black Hole Mass Estimates

In order to estimate the mass of the quasar’s SMBH using scaling relations, we
need to estimate the FWHM of the broad emission lines as well as the
monochromatic luminosity $L_{\lambda}$ (see Eqn. 1). In order to derive the
FWHM of the emission lines with multiple Gaussian components, we first combine
both profiles to form a joint line profile and measure the width of this joint
profile at the half maximum. The estimates for the monochromatic luminosities
$L_{\lambda}$ are determined from the continuum flux $f_{\lambda}$, for which
we evaluate the power-law continuum fit at the respective wavelength $\lambda$
(we evaluate $f_{1350}$ and $f_{3000}$ from the power-law fit to the ground-
based spectrum, i.e. $f_{\lambda}\propto\lambda^{-1.6}$, and $f_{5100}$ based
on the power-law fit to the grism spectrum, i.e.
$f_{\lambda}\propto\lambda^{-2.4}$). The uncertainties on both the FWHM and
$L_{\lambda}$ are estimated from the $16$th and $84$th percentile of $10,000$
random draws from the MCMC posterior.

In order to estimate the black hole mass based on the C IV emission line we
follow the procedure suggested in Coatman et al. (2017). It is well known that
the C IV emission line in high-redshift quasars is often blue shifted due to
strong winds and outflows in the BLR (e.g. Meyer et al., 2019; Schindler et
al., 2020), and thus Coatman et al. (2017) report better consistency in their
black hole mass estimates once they “correct” the FWHM of the C IV line by its
velocity shift from the systemic redshift of $6.3270\pm 0.0005$ which was
determined from the sub-mm [C II] emission using observations from the Atacama
Large Millimetre Array (ALMA) (Wang et al., 2019). However, determining the
velocity shift of the C IV emission line is challenging for for J0100+2802
since it is affected by significant foreground and telluric absorption. Thus,
to reduce the impact of narrow absorption systems or telluric effects on the
emission-line profile we define a “pseudo continuum” within the wavelength
interval $1450-1600$ Å by applying a median filter to the quasar spectrum, as
suggested by (Coatman et al., 2016). Pixels within the wavelength range around
the C IV emission line profile that lie more than $2\sigma$ below the pseudo-
continuum are deemed to be affected by absorption and are additionally masked.
All measured spectral properties are listed in Tab. 2.

These spectral properties allow us to obtain three different estimates of the
black hole mass based on the C IV, Mg II and H$\beta\,$ emission using Eqn. 1.
Reassuringly, we find consistent results between the black hole mass
measurements derived from the rest-frame UV as well as the rest-frame optical
emission lines. All estimates point to a SMBH mass of J0100+2802 between
$9.7\leq\log_{10}(M_{\bullet}/M_{\odot})\leq 10.2$ (see Tab. 2). The
observational uncertainties on these measurements are obtained using again the
$16$th and $84$th percentile of black hole mass estimates from $10,000$ random
draws for the spectral fit parameters from their MCMC posteriors. While the
statistical errors on these measurements are small, there are large systematic
uncertainties in the applied scaling relations of approximately $0.4-0.5$ dex
(e.g. Vestergaard & Peterson, 2006), and thus all estimates of the SMBHs are
consistent within the expected uncertainties. The measurements show that this
quasar indeed hosts a ten billion solar mass black hole at $z=6.3270$, which
is the most massive SMBH known in the early universe.

Figure 2: Median rest-frame optical quasar spectra from low-redshift ($0\leq
z\leq 1$) SDSS quasars, binned by the quasars’ bolometric luminosity (left) or
EW(H$\beta\,$) (right). All spectra are normalized to unity at $5100$ Å. The
spectra show a clear decrease in the strength of the narrow [O III] emission
lines with increasing $L_{\rm bol}$ and decreasing EW. The quasar J0100+2802
is with a bolometric luminosity of $\log_{10}(L_{\rm bol}/\rm
erg\,s^{-1})\approx 48$ more luminous than any of the low-redshift SDSS
quasars. However, low-redshift SDSS quasars with a similar EW(H$\beta\,$)
compared to J0100+2802 (i.e. EW(H$\beta\,$)$=215\pm 2$ Å) exhibit clearly
narrow [O III] emission, which is not present in the spectrum of J0100+2802.

###  3.3 Equivalent Width of the Quasar’s [O III] Emission

We will now take a closer look at the spectral shape of the quasar around the
H$\beta\,$+[O III] emission line complex. We estimate an equivalent width (EW)
in the observed wavelength frame for the [O III] doublet emission of EW([O
III]$)\approx 86\pm 1$ Å. Interestingly, nearly all of this emission arises
from the broad [O III] components, since we do not find any narrow [O III]
emission (i.e. EW([O III]${}_{\rm narrow})<1$ Å) in the spectrum of J0100+2802
(see inset in Fig. 1). It is not unusual to observe weaker narrow [O III]
emission lines with increasing bolometric luminosity (e.g. Baldwin, 1977) and
decreasing equivalent width of the H$\beta\,$ emission in low-redshift quasar
spectra from the Sloan Digital Sky Survey (SDSS; Ahumada et al., 2020; Lyke et
al., 2020), as shown in the median low-redshift SDSS quasar spectra in Fig. 2.
However, these low-redshift quasars at $0\leq z\leq 1$ clearly exhibit narrow
[O III] emission lines even if the quasars are bright (i.e. $L_{\rm
bol}\gtrsim 10^{47}\rm erg\,s^{-1}$), and have a H$\beta\,$ equivalent width
comparable to J0100+2802 (i.e. EW(H$\beta\,$)$=215\pm 2$ Å). At higher-
redshifts, however, this situation changes. Vietri et al. (2018) study the
rest-frame UV and optical spectra of 18 hyper-luminous quasars at $2\lesssim
z\lesssim 4$ and show that weak [O III] emission is common in about $70\%$ of
their sample and is strongly correlated with large blueshifts in the C IV
emission line, as also seen here for J0100+2802. In a similar study using
$330$ luminous quasars at $1.5\leq z\leq 4.0$ Coatman et al. (2019) confirm
the observed anti-correlation between the [O III] EW and C IV blueshifts, and
find that a signifcant fraction of approximately $10\%$ of their quasars
exhibit very weak narrow [O III] emission (i.e. EW$<1$ Å) as seen in
J0100+2802, which is $10$ times higher than among the lower-redshift and
lower-luminosity SDSS quasars. Thus, the absence of narrow [O III] emission in
the hyper-luminous quasar J0100+2802 could point towards evolutionary effects
in the quasar spectra and might be a sign that the narrow line region (NLR)
around the quasar has already been cleared by strong quasar-driven winds on a
relatively short timescale (Vietri et al., 2018; Coatman et al., 2019).
Alternatively, it has been suggested recently that under-luminous [O III]
emission in quasar spectra might be linked to a low gas content in the NLR of
the quasar (Agostino et al., 2022). Upcoming additional rest-frame optical
spectra of high-redshift quasars observed with JWST will be required to shed
more light on the evolution of the NLR at high redshifts.

Table 2: Spectral properties of the quasar J0100+2802.  spectral property | measurement
---|---
$\rm FWHM_{\rm C\,IV}~{}[\rm km\,s^{-1}]$ (corrected) | 3190±2070
$\rm FWHM_{\rm Mg\,II}~{}[\rm km\,s^{-1}]$ | 3890±40
$\rm FWHM_{\rm H\beta}~{}[\rm km\,s^{-1}]$ | 6730±170
$\Delta v_{\rm C\,IV}~{}[\rm km\,s^{-1}]$ | -8850±2970
$\Delta v_{\rm Mg\,II}~{}[\rm km\,s^{-1}]$ | -980±20
$\Delta v_{\rm H\beta}~{}[\rm km\,s^{-1}]$ | -30±130
$1350{\rm{\AA}}\,L_{1350\rm{\AA}}~{}[10^{46}\,\rm erg\,s^{-1}]$ | 40.9±0.1
$3000{\rm{\AA}}\,L_{3000\rm{\AA}}~{}[10^{46}\,\rm erg\,s^{-1}]$ | 25.7±0.1
$5100{\rm{\AA}}\,L_{5100\rm{\AA}}~{}[10^{46}\,\rm erg\,s^{-1}]$ | 20.1±0.1
$\log_{10}(M_{\bullet}/M_{\sun})$ (C IV) | 9.6±0.9
$\log_{10}(M_{\bullet}/M_{\sun})$ (Mg II) | 9.7±0.1
$\log_{10}(M_{\bullet}/M_{\sun})$ (H$\beta\,$) | 10.2±0.1

##  4 No Evidence for Strong Gravitational Lensing

Figure 3: PSF modeling of the quasar images in the NIRCam filters F115W (top),
F200W (middle) and F356W (bottom). The left columns show a $5\arcsec\times
5\arcsec$ cutouts of the quasars, the middle columns show the quasar model,
i.e. a single point source convolved with the PSF model, while the right
columns show the residuals of the quasar after subtracting the model,
normalized by the peak pixel value of the quasar. The black dashed circle
indicates the most probable image separation of $\Delta\theta\approx
0.8\arcsec$ for a lensed source at the quasar’s redshift (see Fig. 4).

All estimates of the SMBH mass rely on the assumption that the observed
luminosity of the quasar is intrinsic to the quasar itself. If the quasar’s
luminosity was, however, magnified by a factor of $\mu$ due to strong
gravitational lensing, the black hole masses would be overestimated by a
factor of $\mu^{\gamma}$ (with $\gamma\approx 0.5$; see Tab. 1).

Previous studies searched for multiple images of the quasar or deviations from
a single point spread function (PSF) due to strong gravitational lensing
effects based on high spatial resolution HST imaging, but have been
unsuccessful. The observed quasar images are consistent with a single PSF
indicating no evidence for multiple images or arcs with a spatial separation
of the diffraction limit of $\theta\approx 0.10\arcsec$, using the F850LP
filter at a central wavelength of $\lambda=0.914\,\mu$m on HST’s Advanced
Camera for Surveys (ACS; Fujimoto et al., 2020, Yue et al. in prep.).

The diameter of JWST’s mirror is, at $6.5$ meters, nearly $3$-times larger
compared to the $2.4$-meter mirror of HST, which reduces the diffraction limit
by a factor of two to $\theta\approx 0.05\arcsec$ in the F115W filter at a
central wavelength of $\lambda=1.154\,\mu$m. In the F200W (F356W) filter with
a longer central wavelength $\lambda=1.989\,\mu$m ($\lambda=3.568\,\mu$m) the
diffraction limit approaches $\theta\approx 0.08\arcsec$ ($\theta\approx
0.14\arcsec$). With the improved spatial resolution of JWST we will now
readdress the question, whether we can find any evidence for strong
gravitational lensing in this ultra-luminous quasar by detailed modeling of
the quasar’s PSF.

###  4.1 PSF modeling and subtraction

Making use of the photutils package (Bradley et al., 2020), we build an
effective PSF by means of the brightest stars in the quasar field in the
respective filter. To this end we first select stars by searching for peaks
($\geq 30~{}\rm MJy\,sr^{-1}$) in the stacked images, while excluding all
peaks that belong to any extended sources rather than stars as well as stars
too close ($\leq 2.5\arcsec$) to the edge of our field of view. This results
in $5$ ($3$, $4$) stars for the F155W (F200W, F356W) filter from which we
build an effective PSF model (Anderson & King, 2000).

Note that since the spectral energy distribution (SED) of the quasar is much
redder than the SED of stars, we expect small differences in the shape of the
PSFs. However, we nevertheless obtain better results when modeling the PSF
when building an effective PSF using stars in the field than when applying
webbpsf555https://webbpsf.readthedocs.io/en/latest/, which allows the user to
input a quasar SED to construct a PSF model.

Leveraging the capabilities of psfMC, which performs 2D surface brightness
modelling on astronomical images using a MCMC algorithm (Mechtley, 2014;
Marshall et al., 2021), we model the quasar as a single point source and find
the best fit PSF model to the quasar images. Note that psfMC requires a noise
model for the effective PSF. To this end we assume a SNR$=100$ of the
effective PSF at each pixel with a fixed noise floor for all pixels below the
$20$th percentile of pixel values.

The results of the PSF modeling and subtraction are shown in Fig. 3 for the
three different NIRCam filters F115W, F200W and F356W. The first column shows
the quasar images, the second column shows our best fit model, i.e. a single
point source convolved with the PSF model, while the last column shows the
residuals of the quasar after subtracting the PSF model. We do not find any
evidence from deviations of a single PSF in any of the three filters.

In order to further test whether we can find any evidence for a second quasar
image, we repeat the modeling procedure and fit the quasar images with two
PSFs instead of a single one. The best fit parameters for the central pixels
of the two models are less than $1$ pixel (i.e. $0.03\arcsec$) apart, when the
magnitudes of both images are approximately equal. Thus we do not find any
evidence for two quasar images and conclude that any possible effects due to
strong gravitational lensing of the source have to result in image separations
of less than a couple of pixels, i.e. below the diffraction limit of
$\theta\approx 0.05\arcsec$.

###  4.2 Searching for a foreground deflector galaxy

We also search for a possible deflector galaxy in the foreground of the
quasar, which could cause the quasar to be gravitationally lensed. To this
end, we search for any extended emission close to the quasar sightline, since
we expect any low-redshift deflector galaxy to be spatially extended given the
size distribution of low-redshift galaxies. In van der Wel et al. (2014), the
authors report a median effective radius of $R\approx 1-4$ kpc ($R\approx 3-6$
kpc) for early-type (late-type) galaxies with a stellar mass of $M_{\star}\sim
5\times 10^{10}\,M_{\odot}$ at redshifts $0<z<3$, which corresponds to spatial
extents of $\sim 0.1\arcsec-0.5\arcsec$ ($\sim 0.4\arcsec-0.8\arcsec$).

We do not see any extended sources close to the quasar sightline in the NIRCam
filters F115W, F200W and F356W down to the $5\sigma$ maximum sensitivity of
28.7, 29.2, 29.0 magnitude, respectively, in the deepest area of the mosaic
around the quasar (Kashino et al., 2022). Furthermore, HST imaging in the
filter F606W at an effective wavelength of $\lambda_{\rm eff}=5776.43$ Å,
where the quasar’s emission is heavily suppressed thanks to the intervening
intergalactic medium acting as a natural band pass filter, does not reveal any
extended sources around the quasar sightline (Yue et al. in prep.).

If a deflecting object would be right along our line-of-sight to the quasar,
we would only be able to see the object in the images after subtracting the
quasar’s PSF, since the quasar would outshine all foreground objects. However,
we do not see any significant extended emission around the quasar or nearby,
even in the PSF subtracted residual images (right panels of Fig. 3)

###  4.3 Probability for strong gravitational lensing

Figure 4: Distribution of the lensing separation $\Delta\theta$ for sources at
the quasar’s redshift $z_{\rm s}=6.327$, predicted by the analytical model
described in Yue et al. (2022) and Oguri & Marshall (2010, OM10). The
diffraction limit obtained from the JWST/NIRCam F115W quasar image is shown as
the dashed line.

We will now estimate the likelihood that the quasar is affected by strong
gravitational lensing given the aforementioned constraints on the Einstein
radius. We will use two models to describe the population of low-redshift
deflector galaxies described in detail in Oguri & Marshall (2010) and Yue et
al. (2022). Both models use singular isothermal spheres (SIS) to describe the
mass profile of the low-redshift deflector galaxies and parameterize the
velocity dispersion function (VDF) of the deflectors with a Schechter function
(Schechter, 1976). However, the VDF applied in Oguri & Marshall (2010) is
based on early-type galaxies derived from SDSS data (Choi et al., 2007), while
the VDF in Yue et al. (2022) is based on both late- and early-type galaxies in
the local universe (Hasan & Crocker, 2019) and assumes a redshift evolution in
the parameters of the Schechter function as suggested by Geng et al. (2021).
Note that the model by Oguri & Marshall (2010) predicts a lower number of less
massive galaxies compared to Yue et al. (2022), since the VDF is based on
early-type galaxies only (see Fig. $1$ in Yue et al. (2022) for a comparison
between the VDFs).

Making use of these VDFs of the low-redshift deflector galaxy population, we
can built an analytic lensing model of the high-redshift lensed quasar
population (Yue et al., 2022). We show the distribution of image separations
of all lensed objects at the quasar’s redshift of $z=6.327$ in Fig. 4. Both
deflector models peak at approximately $\Delta\theta\approx 0.8\arcsec$
indicating the most probable image separation. Based on the PSF modeling of
the quasar in the JWST/NIRCam images described in the previous section we can
constrain the maximum image separation to be below the diffraction limit.
Thus, the fraction of lensed objects with image separation
$\Delta\theta<0.05\arcsec$ among all lensed objects at $z=6.327$ — which
constitutes the probability that J0100+2802 is gravitationally lensed with an
image separation below the diffraction limit — is
$P(\Delta\theta<0.05\arcsec)\approx 2.2\times 10^{-3}$ for the model by Yue et
al. (2022), and $P(\Delta\theta<0.05\arcsec)\approx 6.4\times 10^{-4}$ for the
model by Oguri & Marshall (2010).

##  5 Summary & Discussion

In this letter we present the rest-frame UV as well as the rest-frame optical
spectrum at observed wavelengths between $0.6\mu{\rm m}\leq\lambda_{\rm
obs}\leq 2.3\mu$m and $3.1\mu{\rm m}\leq\lambda_{\rm obs}\leq 4.0\mu$m,
respectively, of the ultra-luminous quasar J0100+2802 at a redshift of
$z=6.327$. The ground-based data consist of a total of $16.8$ hours observed
with the VLT/X-Shooter spectrograph and Magellan/FIRE, while the rest-frame
optical data are observed with JWST/NIRCam in WFSS mode for a total of $9.7$
hours on source.

We measure the mass of the quasar’s SMBH based on the C IV, Mg II and
H$\beta\,$ emission line properties and find estimates that are consistent
between the different line estimators within their systematic uncertainties of
$0.5$ dex, indicating a SMBH mass of $M_{\bullet}\sim 10^{10}\,M_{\odot}$. In
the near future, larger samples of rest-frame optical spectra of high-redshift
quasars will be able to confirm this long-questioned consistency in black hole
mass measurements based on rest-frame UV and rest-frame optical emission line
properties for a wider range of black hole masses and quasar properties.

We revisit the question whether the luminosity of this ultra-luminous quasar
is magnified due to strong gravitational lensing with the increased spatial
resolution of the JWST/NIRCam images compared to previous HST images, and thus
search for evidence for multiple images of the quasars or arcs in three
different NIRCam filters, i.e. F115W, F200W and F356W. We do not find any
evidence from deviations of a single point source at the diffraction limit,
nor do we find any evidence for a foreground deflector galaxy. We thus
estimate the probability of J0100+2802 being affected by strong gravitational
lensing with image separations below the diffraction limit to be $\lesssim
2.2\times 10^{-3}$ using two different models for the foreground galaxy
population (Oguri & Marshall, 2010; Yue et al., 2022).

Our results confirm that J0100+2802 indeed hosts a SMBH with
$M_{\bullet}\approx 10^{10}\,M_{\odot}$ at $z=6.327$, when the universe is
only about $800$ Myr old, which challenges our current understanding of black
hole growth. It requires us to explore models beyond the standard black hole
formation paradigm, such as for instance massive initial black hole seeds in
excess of stellar remnants (e.g. Oh & Haiman, 2002; Begelman et al., 2006), or
radiatively inefficient black hole accretion episodes (e.g. Begelman &
Volonteri, 2017; Davies et al., 2019). JWST has ushered us into a new era for
exploring the early universe. Future observations of high-redshift quasars
will enable us to study the black hole masses for ensembles of quasars,
understand their large-scale environments, and measure their luminosity
functions to gain new insights onto the properties and origins of SMBH seeds,
and their rapid growth at early cosmic times.

The authors would like to thank Paul Schechter for insightful discussions
about gravitational lensing. Furthermore, we would like to Madeline Marshall
for helpful discussions on the PSF modeling, and Gisella De Rosa for
interesting discussions on quasar evolution. This work is based on
observations collected at the European Organisation for Astronomical Research
in the Southern Hemisphere under ESO programme 096.A-0095. This paper includes
data gathered with the 6.5 meter Magellan Telescopes located at Las Campanas
Observatory, Chile. This work is based in part on observations made with the
NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the
Mikulski Archive for Space Telescopes at the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These
observations are associated with program ID $\\#1243$. The specific
observations analyzed can be accessed via https://doi.org/10.17909/xbs2-v060
(catalog 10.17909/xbs2-v060). DK has been supported by JSPS KAKENHI Grant
Number JP21K13956.

## References

  * Agostino et al. (2022) Agostino, C. J., Salim, S., Ellison, S. L., Bickley, R. W., & Faber, S. M. 2022, arXiv e-prints, arXiv:2211.07679.  https://arxiv.org/abs/2211.07679
  * Ahumada et al. (2020) Ahumada, R., Prieto, C. A., Almeida, A., et al. 2020, ApJS, 249, 3, doi: 10.3847/1538-4365/ab929e
  * Anderson & King (2000) Anderson, J., & King, I. R. 2000, PASP, 112, 1360, doi: 10.1086/316632
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bañados et al. (2016) Bañados, E., Venemans, B. P., Decarli, R., et al. 2016, ApJS, 227, 11, doi: 10.3847/0067-0049/227/1/11
  * Bañados et al. (2018) Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Nature, 553, 473, doi: 10.1038/nature25180
  * Baldwin (1977) Baldwin, J. A. 1977, ApJ, 214, 679, doi: 10.1086/155294
  * Begelman & Volonteri (2017) Begelman, M. C., & Volonteri, M. 2017, MNRAS, 464, 1102, doi: 10.1093/mnras/stw2446
  * Begelman et al. (2006) Begelman, M. C., Volonteri, M., & Rees, M. J. 2006, MNRAS, 370, 289, doi: 10.1111/j.1365-2966.2006.10467.x
  * Bochanski et al. (2009) Bochanski, J. J., Hennawi, J. F., Simcoe, R. A., et al. 2009, The Publications of the Astronomical Society of the Pacific, 121, 1409, doi: 10.1086/648597
  * Bradley et al. (2020) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2020, astropy/photutils: 1.0.0, 1.0.0, Zenodo, doi: 10.5281/zenodo.4044744
  * Choi et al. (2007) Choi, Y.-Y., Park, C., & Vogeley, M. S. 2007, ApJ, 658, 884, doi: 10.1086/511060
  * Clough et al. (2005) Clough, S. A., Shephard, M. W., Mlawer, E. J., et al. 2005, J. Quant. Spec. Radiat. Transf., 91, 233, doi: 10.1016/j.jqsrt.2004.05.058
  * Coatman et al. (2016) Coatman, L., Hewett, P. C., Banerji, M., & Richards, G. T. 2016, MNRAS, 461, 647, doi: 10.1093/mnras/stw1360
  * Coatman et al. (2017) Coatman, L., Hewett, P. C., Banerji, M., et al. 2017, MNRAS, 465, 2120, doi: 10.1093/mnras/stw2797
  * Coatman et al. (2019) —. 2019, MNRAS, 486, 5335, doi: 10.1093/mnras/stz1167
  * Connor et al. (2021) Connor, T., Stern, D., Bañados, E., & Mazzucchelli, C. 2021, ApJ, 922, L24, doi: 10.3847/2041-8213/ac37b5
  * Davies et al. (2019) Davies, F. B., Hennawi, J. F., & Eilers, A.-C. 2019, ApJ, 884, L19, doi: 10.3847/2041-8213/ab42e3
  * Davies et al. (2020) Davies, F. B., Wang, F., Eilers, A.-C., & Hennawi, J. F. 2020, ApJ, 904, L32, doi: 10.3847/2041-8213/abc61f
  * De Rosa et al. (2014) De Rosa, G., Venemans, B. P., Decarli, R., et al. 2014, The Astrophysical Journal, 790, 145, doi: 10.1088/0004-637X/790/2/145
  * Eilers et al. (2017) Eilers, A.-C., Davies, F. B., Hennawi, J. F., et al. 2017, ApJ, 840, 24, doi: 10.3847/1538-4357/aa6c60
  * Eilers et al. (2022) Eilers, A.-C., Hogg, D. W., Schölkopf, B., et al. 2022, ApJ, 938, 17, doi: 10.3847/1538-4357/ac8ead
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306, doi: 10.1086/670067
  * Fujimoto et al. (2020) Fujimoto, S., Oguri, M., Nagao, T., Izumi, T., & Ouchi, M. 2020, ApJ, 891, 64, doi: 10.3847/1538-4357/ab718c
  * Geng et al. (2021) Geng, S., Cao, S., Liu, Y., et al. 2021, MNRAS, 503, 1319, doi: 10.1093/mnras/stab519
  * Grier et al. (2017) Grier, C. J., Trump, J. R., Shen, Y., et al. 2017, ApJ, 851, 21, doi: 10.3847/1538-4357/aa98dc
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hasan & Crocker (2019) Hasan, F., & Crocker, A. 2019, arXiv e-prints, arXiv:1904.00486.  https://arxiv.org/abs/1904.00486
  * Horne (1986) Horne, K. 1986, PASP, 98, 609, doi: 10.1086/131801
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jones et al. (2001) Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python.  http://www.scipy.org/
  * Kashino et al. (2022) Kashino, D., Lilly, S. J., Matthee, J., et al. 2022, arXiv e-prints, arXiv:2211.08254.  https://arxiv.org/abs/2211.08254
  * Lyke et al. (2020) Lyke, B. W., Higley, A. N., McLane, J. N., et al. 2020, ApJS, 250, 8, doi: 10.3847/1538-4365/aba623
  * Marshall et al. (2021) Marshall, M. A., Wyithe, J. S. B., Windhorst, R. A., et al. 2021, MNRAS, 506, 1209, doi: 10.1093/mnras/stab1763
  * Matthee et al. (2022) Matthee, J., Mackenzie, R., Simcoe, R. A., et al. 2022, arXiv e-prints, arXiv:2211.08255.  https://arxiv.org/abs/2211.08255
  * Mazzucchelli et al. (2017) Mazzucchelli, C., Bañados, E., Venemans, B. P., et al. 2017, ApJ, 849, 91, doi: 10.3847/1538-4357/aa9185
  * Mechtley (2014) Mechtley, M. 2014, PhD thesis, Arizona State University
  * Meyer et al. (2019) Meyer, R. A., Bosman, S. E. I., & Ellis, R. S. 2019, MNRAS, 487, 3305, doi: 10.1093/mnras/stz1504
  * Morey et al. (2021) Morey, K. A., Eilers, A.-C., Davies, F. B., Hennawi, J. F., & Simcoe, R. A. 2021, ApJ, 921, 88, doi: 10.3847/1538-4357/ac1c70
  * Oguri & Marshall (2010) Oguri, M., & Marshall, P. J. 2010, MNRAS, 405, 2579, doi: 10.1111/j.1365-2966.2010.16639.x
  * Oh & Haiman (2002) Oh, S. P., & Haiman, Z. 2002, ApJ, 569, 558, doi: 10.1086/339393
  * Pacucci & Loeb (2020) Pacucci, F., & Loeb, A. 2020, ApJ, 889, 52, doi: 10.3847/1538-4357/ab6130
  * Park et al. (2022) Park, D., Barth, A. J., Ho, L. C., & Laor, A. 2022, ApJS, 258, 38, doi: 10.3847/1538-4365/ac3f3e
  * Planck Collaboration et al. (2018) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2018, arXiv e-prints.  https://arxiv.org/abs/1807.06209
  * Prochaska et al. (2020) Prochaska, J., Hennawi, J., Westfall, K., et al. 2020, The Journal of Open Source Software, 5, 2308, doi: 10.21105/joss.02308
  * Schechter (1976) Schechter, P. 1976, ApJ, 203, 297, doi: 10.1086/154079
  * Schindler et al. (2020) Schindler, J.-T., Farina, E. P., Bañados, E., et al. 2020, ApJ, 905, 51, doi: 10.3847/1538-4357/abc2d7
  * Schmidt (1965) Schmidt, M. 1965, ApJ, 141, 1295, doi: 10.1086/148217
  * Simcoe et al. (2013) Simcoe, R. A., Burgasser, A. J., Schechter, P. L., et al. 2013, PASP, 125, 270, doi: 10.1086/670241
  * van der Wel et al. (2014) van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28, doi: 10.1088/0004-637X/788/1/28
  * Vernet et al. (2011) Vernet, J., Dekker, H., D’Odorico, S., et al. 2011, A&A, 536, A105, doi: 10.1051/0004-6361/201117752
  * Vestergaard & Osmer (2009) Vestergaard, M., & Osmer, P. S. 2009, ApJ, 699, 800, doi: 10.1088/0004-637X/699/1/800
  * Vestergaard & Peterson (2006) Vestergaard, M., & Peterson, B. M. 2006, ApJ, 641, 689, doi: 10.1086/500572
  * Vestergaard & Wilkes (2001) Vestergaard, M., & Wilkes, B. J. 2001, ApJS, 134, 1, doi: 10.1086/320357
  * Vietri et al. (2018) Vietri, G., Piconcelli, E., Bischetti, M., et al. 2018, A&A, 617, A81, doi: 10.1051/0004-6361/201732335
  * Wang et al. (2019) Wang, F., Wang, R., Fan, X., et al. 2019, ApJ, 880, 2, doi: 10.3847/1538-4357/ab2717
  * Wang et al. (2021) Wang, F., Yang, J., Fan, X., et al. 2021, ApJ, 907, L1, doi: 10.3847/2041-8213/abd8c6
  * Wu & Shen (2022) Wu, Q., & Shen, Y. 2022, arXiv e-prints, arXiv:2209.03987.  https://arxiv.org/abs/2209.03987
  * Wu et al. (2015) Wu, X.-B., Wang, F., Fan, X., et al. 2015, Nature, 518, 512, doi: 10.1038/nature14241
  * Yang et al. (2021) Yang, J., Wang, F., Fan, X., et al. 2021, ApJ, 923, 262, doi: 10.3847/1538-4357/ac2b32
  * Yue et al. (2022) Yue, M., Fan, X., Yang, J., & Wang, F. 2022, ApJ, 925, 169, doi: 10.3847/1538-4357/ac409b


# A Comparison of Two Lattice Boltzmann Models for Electrodynamics

Jorge Isaac Rubiano Murcia<EMAIL_ADDRESS>National University of
Colombia, Faculty of Sciences, Department of Physics, Bogotá, Colombia.
Alejandro Mario Salas Estrada<EMAIL_ADDRESS>National University of
Colombia, Faculty of Sciences, Department of Physics, Bogotá, Colombia.  Jose
David Hernandez Ortega<EMAIL_ADDRESS>National University of
Colombia, Faculty of Sciences, Department of Physics, Bogotá, Colombia.

###### Abstract

In recent years, various Lattice Boltzmann models for electrodynamics have
been developed as alternatives to classical methods such as Finite Difference
Time Domain (FDTD) and Finite Element Methods (FEM). However, there has been a
lack of systematic comparisons between these models. This paper addresses this
gap by comparing two specific Lattice Boltzmann models, published by Mendoza
and Muñoz (MM), and Hauser and Verhey (HV), respectively.

To compare the models, we utilize time and memory as indicators, considering
the same achieved error, in four standard tests: a dielectric pulse traveling
through two interfaces, the skin effect, the Hertz dipole, and a dielectric
pulse traveling through several interfaces.

The results indicate that both methods accurately simulate the tests and
exhibit convergence as the mesh is refined. However, the MM method outperforms
the HV method regarding time, while its memory efficiency was lower. The
modified Hauser-Verhey model demonstrates itself to be a promising alternative
to the Mendoza-Muñoz model. These findings contribute to the ongoing
development and optimization of numerical methods for electromagnetics
simulations.

###### Index Terms:

Lattice Boltzmann, Maxwell’s equations, distribution function, convergence
analysis.

##  I Introduction

The Lattice-Boltzmann method (LBM), although a recent numerical approach
compared to other alternatives, has proven to be a powerful tool for
simulating processes modeled by conservation equations [1]. Among the use
cases, we could highlight applications in computational fluid dynamics where
LBM allows for the simulation of complex geometries and multi-phase flows [2].
Similarly, it has also been successfully deployed to simulate
magnetohydrodynamics [3], the wave equation [4], Poisson’s equation [5], and
linear and non-linear Schrödinger equations [6].

More recently, however, it was further extended by Mendoza and Muñoz to
integrate the electrodynamic Maxwell’s equations [7]. Following their
publication, numerous authors proposed alternative formulations to recover
Maxwell’s equations. For example, Succi and collaborators [8] proposed a
scheme to simulate three-dimensional wave propagation in dispersive media;
Hauser and Verhey [9, 10] similarly treated complex media while significantly
reducing the complexity of the scheme, compared with previous approaches. In
this same fashion, other alternatives were presented for dispersive media [11,
12] and electromagnetic waves in one-dimensional photonic crystals [13].

All these publications perform some cornerstone simulations such as point
dipole antennas, the Skin effect, and medium changes, and some of them carry
out encouraging comparisons with well-established methods such as FDTDM [7,
10]. Nevertheless, there is no systematic comparison between Lattice-Boltzmann
formulations for Maxwell’s equations in the current literature. What’s more
concerning, many previous publications suggest some degree of suitability over
competing alternatives without solid arguments.

In this work, we fill this literature gap by comparing the schemes’ CPU time
and memory requirement in [7] and [9]. For this, we simulated a Gaussian pulse
crossing an interface, the Skin effect, the dipole radiation, and a one-
dimensional highly non-uniform media, for different refinements of the grid.

###  I-A Bhatnagar Gross and Krook (BGK) scheme

The Boltzmann equation governs the evolution of an equilibrium function $f$ in
the phase space, representing a scalar of the probability distribution of
physical quantities in the system. In the Bhatnagar-Gross-Krook (BGK) scheme,
the collision operator $\Omega$ is proportional to the deviation of the
equilibrium function, driving the system towards thermal equilibrium.

The Boltzmann equation can be expressed as follows, see [8, 14].

$\frac{\partial f}{\partial t}+\vec{v}\cdot\nabla f=\Omega(f)+T,$ (1)

where $\Omega(f)$ represents the collision operator and $T$ denotes the source
terms that can be included if necessary.

In numerical computations, the velocity set of vectors is discretized into a
finite set of directions, while the physical space is discretized into a mesh.
The Boltzmann equation can be solved using a two-step process: collision and
advection. This can be achieved through the BGK model in the context of the
lattice Boltzmann method. The discretized versions of the collision and
advection steps can be expressed as follows:

Collision step:

$f^{*}(\vec{x},\vec{v},t+\Delta t)=f(\vec{x},\vec{v},t)-\frac{\Delta
t}{\tau}\left(f(\vec{x},\vec{v},t)-f_{\text{eq}}(\vec{x},\vec{v})\right)+T,$
(2)

Advection step:

$f(\vec{x}+\Delta t\vec{v},\vec{v},t+\Delta t)=f^{*}(\vec{x},\vec{v},t+\Delta
t),$ (3)

where $\Delta t$ is the time step, $\tau$ is the relaxation time, and
$f_{\text{eq}}$ is the equilibrium distribution function. These equations
represent the numerical implementation of the BGK scheme within the lattice
Boltzmann method.

###  I-B Maxwell’s equations

For linear media, Maxwell’s equations can be written in terms of the electric
field $\vec{E}$ and magnetic field $\vec{B}$, as well as the electric
displacement field $\vec{D}$ and magnetic induction field $\vec{H}$. The
equations are as follows, see [15].

  1. 1.

Gauss’s law for electric fields:

$\nabla\cdot\vec{D}=\rho,$ (4)

where $\rho$ is the electric charge density.

  2. 2.

Gauss’s law for magnetic fields:

$\nabla\cdot\vec{B}=0.$ (5)

  3. 3.

Faraday’s law of electromagnetic induction:

$\nabla\times\vec{E}=-\frac{\partial\vec{B}}{\partial t}.$ (6)

  4. 4.

Ampere’s law with Maxwell’s addition:

$\nabla\times\vec{H}=\vec{J}+\frac{\partial\vec{D}}{\partial t},$ (7)

where $\vec{J}$ is the electric current density.

These equations are supplemented by the continuity equation:

$\nabla\cdot\vec{J}=-\frac{\partial\rho}{\partial t},$ (8)

which expresses the conservation of charge.

It’s important to note that the electric displacement field $\vec{D}$ is
related to the electric field $\vec{E}$ by:

$\vec{D}=\epsilon\vec{E},$ (9)

where $\epsilon$ is the permittivity of the medium.

Similarly, the magnetic induction field $\vec{H}$ is related to the magnetic
field $\vec{B}$ by:

$\vec{H}=\frac{1}{\mu}\vec{B},$ (10)

where $\mu$ is the permeability of the medium.
In addition, the Guass’s laws can be obtained from the Faradys’s law, Ampere-
Maxwell law and continuity equation, see [9, 7].

These equations govern the behavior of electromagnetic fields in linear media.
Remark: Throughout the entire text, the permittivity and permeability are
denoted as $\epsilon=\epsilon_{r}\epsilon_{0}$ and $\mu=\mu_{r}\mu_{0}$, where
$\epsilon_{0}$ and $\mu_{0}$ are the permittivity and permeability in vacuum.
Meanwhile, $\epsilon_{r}$ and $\mu_{r}$ represent the relative permittivity
and relative permeability, respectively.

##  II MM model

The model introduced by M. Mendoza and D. Muñoz in [7], henceforth referred to
as the MM model, utilizes a local basis consisting of electric vectors
$\vec{e}^{p}_{ij}$ and velocity vectors $\vec{v}^{p}_{i}$ within a cubic cell
D3Q13. Additionally, magnetic vectors $\vec{b}^{p}_{ij}$ are employed in a
D3Q7 cell. The planes are labeled as $p=1,2,3$, while $i=1,2,3,4$ represents
the discretized velocity directions, with four directions per plane. Each
velocity vector corresponds to two electric vectors denoted by $j=0,1$.

In this framework, the direction of the velocity vectors can be interpreted as
the direction of energy flux, i.e., the pointing vector. Specifically, the
velocity vectors are defined as follows:

$\begin{array}[]{l}\vec{v}_{i}^{0}=\sqrt{2}(\cos((2i-1)\pi/4),\sin((2i-1)\pi/4),0)\,,\\\
\vec{v}_{i}^{1}=\sqrt{2}(\cos((2i-1)\pi/4),0,\sin((2i-1)\pi/4))\,,\\\
\vec{v}_{i}^{2}=\sqrt{2}(0,\cos((2i-1)\pi/4),\sin((2i-1)\pi/4))\,,\\\
\vec{v}_{0}=(0,0,0)\,.\end{array}$ (11)

The electric and magnetic vectors are given by

$\begin{array}[]{c}\vec{e}_{i0}^{p}=\frac{1}{2}\vec{v}^{p}_{[(i+2)\quad\bmod
4]+1}\quad,\quad\vec{e}_{i1}^{p}=\frac{1}{2}\vec{v}^{p}_{[i\bmod 4]+1}\\\
\vec{b}_{ij}^{p}=\vec{v}_{i}^{p}\times\vec{e}_{ij}^{p}.\end{array}$ (12)

These vector sets exhibit appropriate sum relations, as demonstrated in
Equation (11) of reference [7].

For the electric and magnetic fields, there exists a distribution function
$f^{p(r)}_{ij}$ associated with each $p,i,j$. Here, $r=0$ denotes the electric
field, while $r=1$ represents the magnetic field. Moreover, two distribution
functions are assigned to the rest direction for both electric and magnetic
fields. Consequently, the total number of distribution functions amounts to
$2\times 2\times 3\times 4+2=50$.

The macroscopic fields are derived from the distribution functions through the
following equations:

$\begin{array}[]{c}\epsilon_{r}\vec{E}^{\prime}=\vec{D}^{\prime}=\sum_{i=1}^{4}\sum_{p=0}^{2}\sum_{j=0}^{1}f_{ij}^{p(0)}\vec{e}_{ij}^{p}\\\
\vec{B}=\sum_{i=1}^{4}\sum_{p=0}^{2}\sum_{j=0}^{1}f_{ij}^{p(1)}\vec{b}_{ij}^{p}\\\
\rho_{c}=f_{0}^{(0)}+\sum_{i=1}^{4}\sum_{p=0}^{2}\sum_{j=0}^{1}f_{ij}^{p(0)}\end{array}$
(13)

It is remarkable to say that the displacement field $\vec{D}^{\prime}$ is an
auxiliary field. The real field which satisfies Maxwell’s equations is

$\frac{\vec{D}}{\epsilon_{r}}=\vec{E}=\vec{E}^{\prime}-\frac{\mu_{0}}{4\epsilon_{r}}\overrightarrow{J}\,,$
(14)

where $\vec{J}$ is the current source in Maxwell’s equation. In particular,
for ohmic elements $\vec{J}=\sigma\vec{E}$, that is,

$\vec{J}=\frac{\sigma}{1+\frac{\mu_{0}\sigma}{4\epsilon_{r}}}\frac{\vec{D}^{\prime}}{\epsilon_{r}}\,.$
(15)

Furthermore, the $H-$field is given by

$\vec{H}=\frac{\vec{B}}{\mu_{r}}\,.$ (16)

To avoid dissipative effects, the dynamics of the Boltzmann equation adopt the
BGK scheme with a collision time $\tau=1/2$.

The equilibrium functions are defined as follows:

$\displaystyle f_{ij}^{p(0)\mathrm{eq}}(\vec{x},t)$
$\displaystyle=\frac{1}{16}\vec{v}_{i}^{p}\cdot\vec{J}+\frac{\epsilon_{r}}{4}\vec{E}\cdot\vec{e}_{ij}^{p}+\frac{1}{8\mu_{r}}\vec{B}\cdot\vec{b}_{ij}^{p}\,,$
(17) $\displaystyle f_{ij}^{p(1)\mathrm{eq}}(\vec{x},t)$
$\displaystyle=\frac{1}{16}\vec{v}_{i}^{p}\cdot\vec{J}+\frac{1}{4}\vec{E}\cdot\vec{e}_{ij}^{p}+\frac{1}{8}\vec{B}\cdot\vec{b}_{ij}^{p}\,,$
$\displaystyle f_{0}^{(0)\mathrm{eq}}(\vec{x},t)$
$\displaystyle=f_{0}^{(1)\mathrm{eq}}(\vec{x},t)=\rho_{c}\,.$

The advection collision terms are applied conventionally. The BGK collision
step:

$\displaystyle f^{p(r)\prime}_{ij}(\vec{x},t)$
$\displaystyle=2f_{ij}^{p(r)eq}(\vec{x},t)-f_{ij}^{p(r)}(\vec{x},t)\,,$ (18)
$\displaystyle f_{0}^{0\prime}(\vec{x},t)$
$\displaystyle=2f_{0}^{0eq}(\vec{x},t)-f_{0}^{0}(\vec{x},t)\,.$

The advection step:

$\displaystyle f_{0}^{0}\left(\vec{x},t+\Delta t\right)$
$\displaystyle=f^{0\prime}_{0}(\vec{x},t)\,,$ (19) $\displaystyle
f_{ij}^{p(r)}\left(\vec{x}+\vec{v}_{i}\Delta t,t+\Delta t\right)$
$\displaystyle=f^{p(r)\prime}_{ij}(\vec{x},t)\,.$

Interpreting the physical meaning of the equilibrium functions is not
straightforward. However, M. Mendoza and D. Muñoz claim that these functions
can be regarded as perturbations in the energy density.

###  II-A Maxwell’s equations obtained by the MM model

Finally, the MM model successfully reproduces the Maxwell equations for linear
non-dispersive media. Its proof is based on the Chapman-Enskog expansion in
[7].

$\begin{array}[]{c}\frac{\partial\rho_{c}}{\partial
t}+\nabla\cdot\vec{J}=0\,,\\\
\nabla\times\vec{E}=-\frac{\partial\vec{B}}{\partial t}\,,\\\
\nabla\times\vec{H}=\vec{J}+\frac{\partial\vec{D}}{\partial t}.\end{array}$
(20)

to ensure compliance with Gauss’s laws, the equations in the MM model must be
satisfied at time $t=0$. In the case of vacuum, the speed of electromagnetic
waves in the MM model is equal to $1/\sqrt{2}$ in automaton units, where
$\epsilon_{0}=1$ and $\mu_{0}=2$. It is worth noting that exceeding this speed
limit can lead to numerical instabilities.

The authors of the model explain that this phenomenon is attributed to the CFL
(Courant-Friedrichs-Lewy) condition. The CFL condition sets a constraint on
the time step in numerical simulations to maintain stability.

##  III HV model

The model proposed by A. Hauser and L. Verhey in [9], referred to as the HV
model, was introduced based on the model proposed by Y. Liun and G. Yan in
[16]. According to the authors, unlike the MM model, the HV model remains
stable even in the presence of non-smooth transitions at interfaces between
media with different permeability and permittivity.

The discretized set of velocity directions is represented by a D3Q7 cubic
cell, with velocity vectors $\vec{v}_{i}$ for $i=1,...,6$, and $\vec{v}_{0}$
as the rest vector. Similarly, for each $i\geq 1$, there is a pair of
electromagnetic vectors $\vec{e}_{i}$ and $\vec{b}_{i}$. Furthermore, for each
electromagnetic vector, there are three scalar distribution functions
associated with the $x$, $y$, and $z$ components, resulting in the
representation of the distribution function as 3D vectors.

The velocity vectors can be expressed as:

$\displaystyle\vec{v}_{1}$ $\displaystyle=(1,0,0),$ $\displaystyle\vec{v}_{2}$
$\displaystyle=(0,1,0),$ $\displaystyle\vec{v}_{3}$ $\displaystyle=(-1,0,0),$
(21) $\displaystyle\vec{v}_{4}$ $\displaystyle=(0,-1,0),$
$\displaystyle\vec{v}_{5}$ $\displaystyle=(0,0,-1),$
$\displaystyle\vec{v}_{6}$ $\displaystyle=(0,0,1).$

In total, there are $2\times 6\times 3=36$ scalar distribution functions or 12
vector distribution functions.

In the original paper [9], there is an error in equation (B1b). The correct
equation should be:

$\sum_{i=1}^{6}\vec{v}_{\alpha,i}\cdot\vec{v}_{\beta,i}=2\delta_{\alpha\beta}.$
(22)

This error impacts the computations in the paper, including the distribution
function. Therefore, considering the correction, the equilibrium distributions
are given by:

$\displaystyle\vec{e}_{i}^{\mathrm{eq}}$
$\displaystyle=\frac{1}{6}\left(\vec{D}-3\vec{v}_{i}\times\frac{\vec{B}}{\mu}\right),$
(23) $\displaystyle\vec{h}_{i}^{\mathrm{eq}}$
$\displaystyle=\frac{1}{6}\left(\vec{B}+3\vec{v}_{i}\times\frac{\vec{D}}{\epsilon}\right).$

Note that the equilibrium functions bear similarities to Maxwell’s equations
when considering $\nabla$ as a vector that is a scalar multiple of
$\vec{v}_{i}$. This notion is particularly relevant in the context of harmonic
plane waves, where $\nabla$ is parallel to the wave vector (and pointing
vector). This observation suggests that the equilibrium functions can be
interpreted as perturbations of the Maxwell’s equations themselves.

Similar to the MM model, we use $r=0,1$ to denote the electric and magnetic
vectors, respectively, such that $\vec{f}_{i}^{(0)}=\vec{e}_{i}$ and
$\vec{f}_{i}^{(1)}=\vec{b}_{i}$.

The macroscopic fields can be obtained from the distribution functions as
follows:

$\displaystyle\epsilon\vec{E}=\vec{D}(\boldsymbol{r},t)$
$\displaystyle=\sum_{i=1}^{6}\vec{e}_{i}(\boldsymbol{r},t),$ (24)
$\displaystyle\mu\vec{H}=\vec{B}(\boldsymbol{r},t)$
$\displaystyle=\sum_{i=1}^{6}\vec{h}_{i}(\boldsymbol{r},t).$

Similarly to the MM model, the relaxation time $\tau=1/2$ is used, and the BGK
collision step is applied:

$\displaystyle\vec{f}_{i}^{(r)^{\prime}}(\vec{x},t)$
$\displaystyle=2\vec{f}_{i}^{(r)eq}(\vec{x},t)-\vec{f}_{i}^{(r)}(\vec{x},t).$
(25)

The advection step is given by:

$\displaystyle\vec{f}_{i}^{(r)}(\vec{x}+\vec{v}_{i}\Delta t,t+\Delta t)$
$\displaystyle=\vec{f}_{i}^{(r)^{\prime}}(\vec{x},t).$ (26)

###  III-A Maxwell’s equations obtained by the HV model

As demonstrated in [9], the equations obtained through the Chapman-Enskog
expansion are:

$\begin{array}[]{c}\nabla\times\vec{E}=-\frac{\partial\vec{B}}{\partial
t}\,,\\\ \nabla\times\vec{H}=\frac{\partial\vec{D}}{\partial t}\end{array}$
(27)

Following the proposal in [8, 14], we introduce a source term $\vec{T}_{i}$
for each $i$ in the electric distribution function ($r=0$), defined as:

$\vec{T}_{i}=-\frac{1}{2}\left(\vec{J}\cdot\vec{v}_{i}\right)\vec{v}_{i},$
(28)

then the equation (25) is modified for $r=0$ as:

$\displaystyle\vec{f}_{i}^{(0)^{\prime}}(\vec{x},t)$
$\displaystyle=2\vec{f}_{i}^{(0)eq}(\vec{x},t)-\vec{f}_{i}^{(0)}(\vec{x},t)+\vec{T}_{i}.$
(29)

Thus, we obtain the Maxwell’s equations:

$\begin{array}[]{c}\nabla\times\vec{E}=-\frac{\partial\vec{B}}{\partial
t}\,,\\\ \nabla\times\vec{H}=\vec{J}+\frac{\partial\vec{D}}{\partial
t}\end{array}$ (30)

The HV model with the additional source term is referred to as the modified HV
model.

In the HV model, the speed of electromagnetic waves in the vacuum is equal to
$1/3$ in lattice units. Values greater than this limit can result in numerical
instabilities. However, the authors of the model do not provide an explanation
for this specific limit.

##  IV Numerical test and comparison

###  IV-A Gaussian Pulse Crossing Dielectric Interface

For the first comparison, we simulated a Gaussian pulse of the form

$\displaystyle\vec{B}$
$\displaystyle=E_{0}/C\exp(-(z-z_{0})^{2}/(2\alpha^{2}))\,\hat{y}\,.$ (31)
$\displaystyle\vec{E}$
$\displaystyle=E_{0}\exp(-(z-z_{0})^{2}/(2\alpha^{2}))\,\hat{x}\,,$

where $\alpha=0.05\cdot L_{z}/\sqrt{2}$ and $z_{0}=L_{z}/2-L_{z}/6$ .

The pulse travels from the vacuum $\epsilon_{1,r}=1$ to a media with
$\epsilon_{2,r}=\epsilon_{r}=2.5$ and $\mu_{r}=1$. According to the theory,
the reflected and transmitted amplitude are such that

$\displaystyle\frac{A_{ref}}{A_{inc}}=\frac{\sqrt{r}-1}{\sqrt{r}+1}\,,$ (32)
$\displaystyle\frac{A_{trans}}{A_{inc}}=\frac{2}{\sqrt{r}+1}\,,$ (33)

being $r=\epsilon_{2,r}/\epsilon_{1,r}$ and, $A_{ref}$ and $A_{trans}$ the
amplitude of the reflected and transmitted pulse, respectively [15].

Figure 1: Normalized electric amplitudes $A/A_{0}$ of a Gaussian pulse
crossing a dielectric interface simulated using the MM model for different
grid refinements labeled by the number of cells $N_{z}$. The region where
$z/L_{z}>0.5$ corresponds to a dielectric medium with $\epsilon_{r}=2.5$,
while the region where $z/Lz<0.5$ corresponds to the vacuum
($\epsilon_{r}=1.0$). The reflected and transmitted pulse can be observed.

Our simulations consisted of a one-dimensional grid, progressively refined
while keeping the physical size unchanged. In figure (1), the shape of the
transmitted and reflected pulse is depicted for different refinements using
the MM model. For each grid, we measured the CPU time only for the integration
of Maxwell’s equations with the LBM, i.e., our measurement excluded the CPU
time for the initialization of variables or printing. Similarly, we also
gauged the relative error in the amplitude of the simulated reflected and
transmitted pulse, obtaining the curves in figure (2). Importantly, the time
measurement was made several times for each refinement to compute an average
and error bars.

(a)

(b)

Figure 2: CPU time versus simulation error in the reflected electric amplitude
for the dielectric pulse test.

The plots show that the MM scheme requires less CPU time for a given relative
error. Remarkably, the relation CPU time-error for the transmitted pulse is
erratic for both simulations. Furthermore, the error bars for both
formulations are generally small, suggesting that this behavior is intrinsic
in the two schemes.

(a)

(b)

Figure 3: Memory usage of the distribution functions for both models as a
function of the relative errors for the dielectric pulse test.

Regarding the behavior of the relative error for different refinement levels,
it was observed in both schemes that the error decreases following a power law
as the grid refinement is increased. However, it is notable that the MM model
exhibits the most significant reduction in numerical errors with the
implementation of grid refinement.

A quick measurement of the memory used by the distribution functions was also
conducted as a function of the relative errors (figure 3). The memory
requirement was observed to increase for lower relative errors in both
schemes. However, the MM model exhibits a higher memory requirement than the
HV model due to the difference in the number of distribution functions, as
mentioned in Sections II and III.

###  IV-B Skin Effect

The Skin effect describes the exponential decay of a plane wave’s amplitude
after penetrating a conducting material. Theoretically, the amplitude of the
electric field inside the conductor is given by the expression:

$A_{\text{Theo}}=A_{0}\exp(-z/\delta)$ (34)

where $A_{0}$ represents the amplitude outside of the conductor, and $\delta$
denotes the skin thickness, which can be expressed as [17, p. 130]:

$\delta=\sqrt{\frac{2}{\sigma\mu\omega}}\sqrt{\sqrt{1+\left(\frac{\omega\epsilon}{\sigma}\right)^{2}}+\frac{\omega\epsilon}{\sigma}}\,.$
(35)

where $\omega$ corresponds to the angular frequency of the wave.

Figure 4: Normalized electric field amplitude inside a conductor simulated
with the HV model. A plane wave encounters a conductive medium at $z/L_{z}\geq
0.25$ with conductivity $\sigma=0.1\cdot\epsilon/T$ .

We simulate a plane wave in one dimension, imposing the electromagnetic field
at $z=0$ with a wavevector to the right.

$\displaystyle\vec{B}$ $\displaystyle=E_{0}/C\sin(\omega t)\hat{y}\,$ (36)
$\displaystyle\vec{E}$ $\displaystyle=E_{0}\sin(\omega t)\hat{x}\,,$

where $\omega=2\pi/T$, and the period is $T=17.68\cdot 10^{-3}L_{z}/C$. A
conductive medium is placed at $z/L_{z}\geq 0.25$ with conductivity
$\sigma=0.1\cdot\epsilon/T$ .

The behavior of the oscillating wave and its amplitude after penetrating the
conductor is presented in the figure (4) for the HV model.

We proceeded similarly to the Gaussian pulse to compare the errors and
computing time. We utilized a one-dimensional grid with several refinements
and measured the time and error for each grid, computing an average and error
for the time. The only difference lies in the error measurement, which was
modeled with a cost function

$C=\sum_{j}\bigg{(}\frac{A_{\text{sim}}-A_{\text{theo}}}{A_{\text{theo}}}\bigg{)}^{2}$
(37)

where $A_{\text{sim}}$ and $A_{\text{theo}}$ are the simulated and theoretical
electric field amplitudes after penetrating the conductor, and the sum is
carried out over the cells occupied by the conductor. As figure (5) shows,
where both time measurements overlap, the MM scheme produces a smaller cost;
additionally, the polynomial fit shows that for ever smaller costs, the MM
model requires less computation time.

(a)

(b)

Figure 5: CPU time versus cost for the Skin effect.

###  IV-C Electric Dipole Radiation

An oscillating electric or magnetic dipole represents one of the simplest
antenna systems available and is extensively treated in many textbooks [15],
the reason why is an ideal example to test performance for simulation of
radiating systems.

The simulation involves a Hertz dipole with specific parameters. Instead of a
point source, a localized density current associated with a dipole is
simulated. The current density is given by

$J=J_{0}\exp\left(-\alpha\left(\vec{x}-\vec{x_{d}}\right)^{2}\right)\sin(\omega
t)\hat{z}.$ (38)

where $\alpha=0.5$, the amplitude of the current density $J_{0}$ is set to
0.0001 and $\vec{x_{d}}$ are the dipole coordinates. The period $T$ is
calculated as $\frac{17.68}{100.0}\frac{Lz}{C}$, and $\omega$ is set to
$\frac{2\pi}{T}$. The associated amplitude of the dipole is [7][15, Chapter 9]

$p=\frac{J_{0}}{\omega}\left(\frac{\pi}{\alpha}\right)^{1.5}\,.$ (39)

The vacuum impedance $Z_{0}$ is calculated as
$\sqrt{\frac{\mu_{0}}{\epsilon_{0}}}$, and the wavelength $\lambda$ is
determined as $c\cdot T$, being $c$ the speed of light. The maximum simulated
time is $t_{\text{max}}=T\cdot\frac{70}{25.0}$ in the time test comparisons.

To obtain the radiation patterns, measurements are taken at a radius
$R=\lambda\cdot n$, where $n=\frac{Lz}{2\lambda}-2$ (two wavelengths less
before reaching the boundary of the lattice domain). We use
$t_{\text{max}}=\frac{R+2\lambda}{C}$, and measure the maximum energy flux
during a whole period in the time interval $R+\lambda\leq t\leq R+2\lambda$.
The normalized theoretical radiation pattern in spherical coordinates
$(\phi,\theta)$, taking the z-axis in the direction of the dipole, is
$\sin^{2}(\theta)$, see [15, Chapter 9]. In figure (7) is plotted the
theoretical and simulated radiation pattern in a plane of constant $\phi$.

Finally, we gauged the same cost function (37) and CPU time for the error and
time measurement. Specifically, we compared the simulated and theoretical
electric and magnetic field amplitudes along a line perpendicular to the
dipole moment. From figure (8) is clear that in the studied interval, both
models had very similar performance. Even more remarkable, for a given
computation time, the electric and magnetic field costs are not the same for a
given LB; furthermore, in the limit of small costs, MM performs better for the
magnetic field but worse for the electric field.

(a) Contour lines of magnetic field $B_{y}/J_{0}$ in the x-z plane simulated
by the MM model in automaton units.

(b) Theoretical and simulated magnetic field along the x-axis using the MM
model in automaton units.

Figure 6: Comparison of magnetic field simulations using the MM model. Figure
7: Normalized theoretical and simulated slice of the radiation pattern of the
Hertz dipole using the HV model.

(a)

(b)

Figure 8: CPU time versus cost for the radiating dipole.

###  IV-D Gaussian Pulse Crossing Non-uniform media

The setup was similar to the section on the dielectric pulse crossing one
interface. However, in this case, we have four different mediums with relative
permitivities of $1,1.3,2,3$. The corresponding interfaces are located at
$z_{0},z_{0}+d,z_{0}+3d$, respectively, where $z_{0}=\frac{L_{z}}{2}$ and
$d=\frac{L_{z}}{20}$.

The simulation space consisted of $L_{z}$ cells, and the pulse was positioned
at $z_{0}-\frac{L_{z}}{30}$ using an $\alpha=\frac{L_{z}}{100\sqrt{2}}$ in
equation (31). To achieve an error of 0.03 $\%$ between the simulated and
theoretical transmitted pulses using both methods, the MM model requires 18
seconds of CPU time and 4000 cells, which corresponds to 200,000 distribution
functions. On the other hand, the HV model only requires 10 seconds of CPU
time and 3000 cells, totaling 108,000 scalar distribution functions. Figures 6
and 7 depict the pulse after crossing all the interfaces.

Figure 9: Dielectric pulse crossing several interfaces using the HV model.
Figure 10: Dielectric pulse crossing several interfaces using the MM model.

##  V Conclusion

Both models are suitable alternatives for simulating electrodynamics
phenomena. They exhibited similar behavior in terms of computation time and
error. However, the HV model outperformed the MM model in a test involving a
Gaussian pulse crossing a non-uniform medium. Despite this, the MM model
generally performed slightly better than the HV model in most tests.
Nevertheless, the HV model demonstrated lower memory consumption due to its
reduced number of distribution functions. Therefore, the HV model may be
preferred as it requires less memory, while yielding similar errors to the MM
model with equal CPU times. It appears that reducing the number of
distribution functions does not significantly affect overall performance but
improves memory requirements. The trade-off of the HV model is that it
requires more iterations because the speed of light in automaton units is
lower than that of the MM model. Future works could explore more LBMs and
further test the validity of this claim.

## Authorship contribution statement

JR implemented the Lattice-Boltzmann algorithms and proposed the modified HV
model. AS and JH performed the performance tests. AS and JR conceived the
project. All three authors wrote and revised the article.

## Acknowledgment

The authors would like to express their gratitude to Dr. José Daniel Muñoz for
his valuable comments and insights throughout the development of this
research. Additionally, they would like to thank Dr. Rafael Rey for suggesting
the dielectric pulse test conducted at multiple interfaces and to thank
Dr.Andreas Hauser for clarifying some of our doubts. Their contributions
greatly enhanced the quality and depth of this study. In addition, ChatGPT and
Grammarly were utilized for writing style grammar correction.

## APPENDIX A: Source term in the HV model

$\sum_{i}\left(\partial_{t}f_{i}^{\mathrm{eq}}+\sum_{\alpha}v_{\alpha_{i}}\partial_{\alpha}f_{i}^{\mathrm{eq}}\right)\approx\sum_{i}T_{i}\,.$
(40)

Following this equation for the equilibrium functions $\vec{f}_{i}^{(0)}$ we
obtain in the left-hand side

$\frac{\partial\vec{D}}{\partial t}-\nabla\times\vec{H}\,.$ (41)

On the other hand, we get in the right-hand side

$\displaystyle\sum_{i=1}^{6}\vec{T}_{i,\alpha}$
$\displaystyle=-\frac{1}{2}\sum_{i=1}^{6}\left(\vec{J}\cdot\vec{v}_{i}\right)\vec{v}_{i,\alpha}\,,$
(42)
$\displaystyle=-\frac{1}{2}\sum_{i=1}^{6}\vec{J}_{\beta}\vec{v}_{i,\beta}\vec{v}_{i,\alpha}\,,$
$\displaystyle=-\frac{1}{2}\vec{J}_{\beta}\sum_{i=1}^{6}\vec{v}_{i,\beta}\vec{v}_{i,\alpha}\,,$
$\displaystyle=-\frac{1}{2}\vec{J}_{\beta}(2\delta_{\alpha,\beta})\,,$
$\displaystyle=-\vec{J}\,.$

Thus the Ampere-Maxwell equation is obtained.

$\displaystyle\frac{\partial\vec{D}}{\partial t}-\nabla\times\vec{H}=-\vec{J}$

## Appendix B: PSEUDO-CODE

In the following pseudocode, we outline the steps involved in the Lattice
Boltzmann Automata simulation for solving a specific problem using the MM and
modified HV models. The pseudocode describes the main algorithmic steps
involved in the simulation, including initialization, collision, advection,
and analysis. The specific equations and considerations for different models
are also highlighted.

1:At $t=0$, impose all fields $\vec{B},\vec{E},\vec{J}$.

2:Initialize all the distribution functions as the equilibrium functions
evaluated in $\vec{B},\vec{E}$ and $\vec{J}$.

3:Impose fields in the cells, if required.

4:for $t=1$ to $t_{\text{max}}$ do

5: // In the collision, compute the macroscopic fields by summing the
distribution functions, i.e, using eqs. (13) to (16) for MM model and (24) for
HV model

6: Collision

7: Impose fields if required.

8: Advection

9:end for

10:Analyze and plot.

Note: Sometimes, for instance, in the skin effect and in the Hertz dipole, the
current is imposed during the collision step. In the MM model, the current is
computed with equation (15) after computing $\vec{D^{\prime}}$ with equation
(13), and then it is used in the equilibrium functions (17) during the
collision step as in equation (18). On the other hand, in the HV model, the
current is computed (for example as $\vec{J}=\sigma\vec{E}$) and then is
passed to the source term $T$ as in equation (28).

## References

  * [1] J. Zhang and D. Y. Kwok, _Lattice Boltzmann Method (LBM)_. Boston, MA: Springer US, 2008, pp. 981–985. [Online]. Available: https://doi.org/10.1007/978-0-387-48998-8˙800
  * [2] L. Rui, L. Chengcheng, and L. Junjun, “Lattice boltzmann method.”
  * [3] P. J. Dellar, “Electromagnetic waves in lattice boltzmann magnetohydrodynamics,” _Europhysics Letters_ , vol. 90, no. 5, p. 50002, jun 2010. [Online]. Available: https://dx.doi.org/10.1209/0295-5075/90/50002
  * [4] Y. Guangwu, “A lattice boltzmann equation for waves,” _Journal of Computational Physics_ , vol. 161, no. 1, pp. 61–69, 2000. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0021999100964866
  * [5] Z. Chai and B. Shi, “A novel lattice boltzmann model for the poisson equation,” _Applied Mathematical Modelling_ , vol. 32, no. 10, pp. 2050–2058, 2008. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0307904X07001722
  * [6] S. Palpacelli, “Quantum lattice boltzmann methods for the linear and nonlinear schrödinger equation in several dimensions,” Ph.D. dissertation, Roma Tre University, available at http://www.matfis.uniroma3.it/Allegati/Dottorato/TESI/palpacel/tesipalpacelli.pdf.
  * [7] M. Mendoza and J. D. Muñoz, “Three-dimensional lattice boltzmann model for electrodynamics,” _Phys. Rev. E_ , vol. 82, p. 056708, Nov 2010. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevE.82.056708
  * [8] S. M. Hanasoge, S. Succi, and S. A. Orszag, “Lattice boltzmann method for electromagnetic wave propagation,” _Europhysics Letters_ , vol. 96, no. 1, p. 14002, sep 2011. [Online]. Available: https://dx.doi.org/10.1209/0295-5075/96/14002
  * [9] A. Hauser and J. L. Verhey, “Stable lattice boltzmann model for maxwell equations in media,” _Phys. Rev. E_ , vol. 96, p. 063306, Dec 2017. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevE.96.063306
  * [10] ——, “Comparison of the lattice-boltzmann model with the finite-difference time-domain method for electrodynamics,” _Phys. Rev. E_ , vol. 99, p. 033301, Mar 2019. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevE.99.033301
  * [11] J. Hussain, R. Dasgupta, H. N. Dixit, S. P. Thampi, and A. Roy, “A lattice boltzmann method for electromagnetic wave propagation in medium,” in _2020 IEEE International Conference on Computational Electromagnetics (ICCEM)_ , 2020, pp. 299–301.
  * [12] P. Chen, C.-H. Wang, and J.-R. Ho, “A lattice boltzmann model for electromagnetic waves propagating in a one-dimensional dispersive medium,” _Computers and Mathematics with Applications_ , vol. 65, no. 6, pp. 961–973, 2013, mesoscopic Methods in Engineering and Science. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S089812211300031X
  * [13] H. Ma, B. Wu, L. Song, H. Ren, W. Jiang, W. Guo, and M. Tang, “A new lattice boltzmann scheme for photonic bandgap and defect mode simulation in one-dimensional plasma photonic crystals,” _Photonics_ , vol. 9, no. 7, 2022\. [Online]. Available: https://www.mdpi.com/2304-6732/9/7/464
  * [14] Z. Guo, C. Zheng, and B. Shi, “Discrete lattice effects on the forcing term in the lattice boltzmann method,” _Phys. Rev. E_ , vol. 65, p. 046308, Apr 2002\. [Online]. Available: https://link.aps.org/doi/10.1103/PhysRevE.65.046308
  * [15] J. D. Jackson, _Classical electrodynamics_ , 3rd ed. New York, NY: Wiley, 1999. [Online]. Available: http://cdsweb.cern.ch/record/490457
  * [16] Y. Liu and G. Yan, “A lattice boltzmann model for maxwell’s equations,” _Applied Mathematical Modelling_ , vol. 38, no. 5, pp. 1710–1728, 2014. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0307904X13005684
  * [17] E. Jordan and K. Balmain, _Electromagnetic Waves And Radiating Systems 2Nd Ed._ , ser. Prentice-Hall electrical engineering series. Prentice-Hall Of India Pvt. Limited, 1995. [Online]. Available: https://books.google.com.co/books?id=91QsNQAACAAJ