of the sources within the bounding box (e.g., the yellow star in the top left
panel of Figure 1). We match our field with the Gaia catalog (Gaia
Collaboration et al., 2016, 2018) and mask out stars outside of the bounding
box. As we have already detected peaks throughout the whole cutout image, we
also generate a mask for objects outside the bounding box to reduce the impact
of scattered light from bright galaxies on the modeling of the flat-sky
source.

The optimization process uses the adaptive proximal gradient method (Melchior
et al., 2019), which is a robust method for optimization with constraints. The
model is considered to be converged when the relative changes of parameters
are smaller than e_rel = 2e-4. Typically, convergence is achieved after $\sim
50$ steps of optimization, and the whole modeling process takes about 40s for
a typical LSBG. It takes longer when the target galaxy has a large angular
size as more peaks are modeled. We note that scarlet only finds the maximum
likelihood estimation of the model instead of deriving the full posteriors,
thus a good initialization becomes especially important for fast convergence.

##  Appendix C Spergel profiles

Figure 8: Correspondence between the Sérsic profile (C3) and the Spergel
profile (C1). We fit a Spergel profile to each Sérsic while fixing the total
luminosity and half-light radius. The black line in the left panel shows the
best-fit Spergel index $\nu$ as a function of the Sérsic index $n$. We also
provide a symbolic approximation to this relation in (C5), shown as the gray
dashed-dotted line. In the right panel, two Sérsic profiles (solid) and their
best-fit Spergel profiles (dashed-dotted) are shown. The Spergel profile
approximates the Sérsic profile well at the small Sérsic index.

Figure 9: Distribution of UPGs (blue) and mock galaxies (red) in $g-i$ color,
ellipticity, and Sérsic index $n$. We translate the Spergel index to the
Sérsic index using Equation (C5). The mock galaxy population represents the
real UDG/UPG population well.

In this appendix, we demonstrate that the Spergel profile can approximate the
Sérsic profile and provide a lookup table for the correspondence between the
Sérsic index $n$ and Spergel index $\nu$.

The surface brightness of a Spergel profile has the form (Spergel, 2010):

$I_{\nu}(r)=\frac{c_{\nu}^{2}L_{0}}{2\pi
r_{0}^{2}}f_{\nu}\left(\frac{c_{\nu}r}{r_{0}}\right),$ (C1)

where

$f_{\nu}(u)=\left(\frac{u}{2}\right)^{\nu}\frac{K_{\nu}(u)}{\Gamma(\nu+1)},$
(C2)

and $K_{\nu}(u)$ is the Modified Bessel function of the second kind. The half-
light radius is $r_{0}$, the total luminosity is $L_{0}$, and $c_{\nu}$
satisfies the equation $(1+\nu)f_{\nu+1}(c_{\nu})=1/4$. The Spergel profile
has a simple analytical expression in Fourier space, making it easy to
convolve with a PSF.

The surface brightness of a Sérsic profile follows (Sérsic, 1963; Graham &
Driver, 2005):

$I(r)=I_{\mathrm{e}}\exp\left\\{-b_{n}\left[\left(\frac{r}{r_{\mathrm{e}}}\right)^{1/n}-1\right]\right\\},$
(C3)

where $r_{e}$ is the half-light radius, $I_{e}$ is the surface brightness at
$r=r_{e}$, $n$ is the Sérsic index. The value of $b_{n}$ satisfies
$\Gamma(2n)=2\gamma\left(2n,b_{n}\right)$, where $\gamma(a,x)$ is the
incomplete gamma function. According to Graham & Driver (2005), the total
luminosity of a Sérsic profile is given by

$L_{0}=I_{e}r_{e}^{2}\,2\pi n\,e^{b_{n}}\left(b_{n}\right)^{-2n}\Gamma(2n).$
(C4)

To study the correspondence between Sérsic and Spergel profiles, we generate
Sérsic 1D profiles with different Sérsic indices using
astropy181818https://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Sersic1D.html,
and try to fit the Sérsic profiles using the Spergel model. The Sérsic index
ranges from $n=0.5$ to $n=4.5$, and both profiles are normalized at the half-
light radius (see the right panel in Figure 8). For each Sérsic profile, we
calculate the total luminosity $L_{0}$ according to (C4) and plug it into the
Spergel profile (C1) as a fixed value. Therefore, only the Spergel index $\nu$
is allowed to vary during the fitting. The fitting uses the least squares
fitting with Levenberg–Marquardt algorithm implemented in
astropy.modeling.fitting191919https://docs.astropy.org/en/stable/api/astropy.modeling.fitting.LevMarLSQFitter.html#astropy.modeling.fitting.LevMarLSQFitter,
which is a commonly used robust algorithm to solve nonlinear least square
fitting problems. The best-fit Spergel index $\nu$ as a function of the Sérsic
index $n$ is shown in the left panel of Figure 8. A Spergel profile with
$\nu=0.5$ is exactly an exponential profile with $n=1$. The de Vaucouleurs
profile (de Vaucouleurs, 1948) with $n=4$ can be approximated by a Spergel
profile with $\nu=-0.74$. As two examples, we show two Sérsic profiles (solid)
and their best-fit Spergel profiles (dashed-dotted) in the right panel.
Overall, the Spergel profile provides a good approximation for the Sérsic
profile.

As shown in Figure 8, a Sérsic profile with a small Sérsic index ($0.5<n<1.5$)
can be well-approximated by a Spergel profile, although the Spergel profile
seems to be more extended than Sérsic in the outskirts. For a Sérsic profile
with a high Sérsic index, the approximation gets worse at both small and large
radii. Overall, the Spergel profile is a good approximation to Sérsic from
$n\approx 0.5$ to $n\approx 4.5$. Using the symbolic regression package PySR
(Cranmer, 2023), the best-fit relation between $n$ and $\nu$ for $0.5<n<4.5$
is given as follows:

$\nu=1.06\,n^{-2.5}\cos(\sin(3.74/n))^{1/2}-0.365\,n^{1/2}$ (C5)

It is well-known that the light profiles of low-mass galaxies are quite flat
and can be described using Sérsic profiles with $0.5<n<1.5$ (e.g., van Dokkum
et al., 2015; Lange et al., 2015; Greco et al., 2018; Zaritsky et al., 2021;
Carlsten et al., 2021). It is thus reasonable to use Spergel profiles to model
the LSBGs (§3.3) and enjoy its convenience in Fourier space.

##  Appendix D Completeness and Measurement Uncertainty

###  D.1 Completeness

The total completeness is a multiplication between the detection completeness
and deblending completeness. Below we describe how we generate mock galaxies
and derive the recovered fraction.

To estimate the detection completeness, we inject $\sim 700,000$ mock galaxies
with single-Sérsic light profiles into the co-added images. The mock galaxies
are generated using GalSim (Rowe et al., 2015) and are injected after being
convolved with HSC PSF generated by hscPipe. We directly inject the mock
galaxies into the co-added images without adding extra noise. The galaxies are
placed randomly on the image but we require all mock galaxies to be separated
by at least 80″ from each other to minimize the impact of mock galaxies on the
image noise properties. We have also done extensive tests on injecting mock
galaxies into the raw images, adding Poisson noise, and going through the
entire data reduction pipeline using SynPipe (Huang et al., 2018b). This is
very expensive in terms of both CPU time and disk space as we must run the
full hscPipe. However, we find no noticeable difference in the completeness
between this method and the direct injection to co-added images.

To resemble the real LSBG population in G18, we generate mock galaxies
following uniform distributions in size ($2\arcsec\leqslant r_{e}\leqslant
21\arcsec$), surface brightness ($23\leqslant\overline{\mu}_{\rm
eff}(g)\leqslant 28.5\ \mathrm{mag\ arcsec^{-2}}$), Sérsic index
($0.8<n<1.2$), and ellipticity ($0<\varepsilon<0.6$). They are randomly
assigned to have a blue ($g-i=0.47,\ g-r=0.32$), medium ($g-i=0.64,\
g-r=0.43$), and red ($g-i=0.82,\ g-r=0.56$) color with equal chance. The
ranges spanned by these colors and Sérsic indices cover most of the observed
LSBGs in G18. Then we run the detection step described in §3.1 and cross-match
the detection catalog with the input mock galaxy catalog to calculate
completeness. We split the size and surface brightness range into 15 bins with
$\Delta r_{e}=0\farcs 86$, $\Delta\overline{\mu}_{\mathrm{eff}}(g)=0.33\
\mathrm{mag\ arcsec}^{-2}$, and interpolate over bins using an isotropic
Gaussian kernel with $\sigma=0.5$. We note that the size and surface
brightness shown in Figure 3 are all the intrinsic values for the mock
galaxies, which are different from the measured ones. We find negligible
dependence of detection completeness on the Sérsic index, color, and
ellipticity. Therefore, we neglect the dependence of detection completeness on
parameters other than size and surface brightness hereafter.202020We performed
a smaller set of image simulations where we extended the ellipticity range and
also simulated galaxies with additional structures (such as star-forming
clumps). We find the completeness declines for $\varepsilon>0.6$, suggesting
that edge-on disk galaxies may be missing from our sample. Please see Greene
et al. (2022) for more details.

For the deblending completeness, we inject 5000 mock Sérsic galaxies into the
co-added images and run vanilla scarlet on them. The mock galaxies follow the
same uniform distribution in size, surface brightness, ellipticity, Sérsic
index as for deriving the detection completeness, but follow a Gaussian
distribution in color: $g-i\sim\mathcal{N}(0.6,0.2^{2}),\
g-r=0.7\cdot(g-i)+\mathcal{N}(0,0.03^{2})$. We show the distribution of mock
galaxies and UPGs in Figure 9. This mock galaxy population represents the real
UDG/UPG population well. To be specific, the Gaussian color prior covers 86%
of the UPGs within $1.5\sigma$, the ellipticity prior covers 98% of them, and
the Sérsic index prior covers 60% of them. We measure the size and surface
brightness on the scarlet models of the mock galaxies, then we apply the
deblending cuts. We also find that the deblending cut mainly depends on the
size and surface brightness.

###  D.2 Measurement Bias and Uncertainty

Due to the low-surface-brightness nature of LSBGs, their size, magnitude,
surface brightness, and shape are hard to characterize and are typically
associated with large uncertainties. Häussler et al. (2007) used mock single-
Sérsic galaxies and parametric fitting codes to demonstrate that the estimated
size and total magnitude are sensitive to local sky estimation and masking of
neighboring objects. In the low-surface-brightness regime, the measured size
and total magnitude can be quite biased and have large uncertainties.
Therefore, measurement bias and uncertainty must be considered when studying
the properties of LSBGs. Zaritsky et al. (2021, 2022) characterized their
measurement bias and uncertainty by injecting mock Sérsic galaxies into co-
added images and comparing the recovered properties with the truth. Tanoglidis
et al. (2022) recently proposed a new method to estimate the measurement error
using a Bayesian neural network. In this paper, we simply take the former
method because of its simplicity.

In order to test how well we recover the photometric and structural parameters
in our measurements (Sec 3.3), we take the 5000 mock Sérsic galaxies used for
computing the deblending completeness and model them using the Spergel light
profile. We model the bias and uncertainty in size, surface brightness, total
magnitude, and color as a function of other parameters including size, surface
brightness, and shape. However, we do not find any significant dependence of
the bias on color, ellipticity, and the Spergel index. Thus we just model the
bias and uncertainty as a function of the measured angular size $r_{e}$ and
surface brightness $\overline{\mu}_{\mathrm{eff}}(g)$.

Given the size and surface brightness of our simulated galaxies, we set the
range of the observed size and surface brightness to be
$r_{e}\in[1\arcsec,15\arcsec],\
\overline{\mu}_{\mathrm{eff}}(g)\in[23,29]\,\mathrm{mag\ arcsec}^{-2}$. We
then split the observed $r_{e}-\overline{\mu}_{\rm eff}(g)$ plane using an
$8\times 8$ grid, and calculate the mean bias $\Delta X=X_{\rm truth}-X_{\rm
meas}$ within each bin, where $X=\\{\overline{\mu}_{\mathrm{eff}}(g),\ g-i,\
g-r\\}$. For $r_{e}$, we calculate the relative bias $\Delta r_{e}/r_{e}$
instead because it is less dependent on the angular size. Then we interpolate
over the grid using a multi-quadratic
kernel212121https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RBFInterpolator.html
with epsilon=0.5, smooth=1. Unlike Zaritsky et al. (2021) where they fit
models in 4-D space, we find that interpolation works well enough in 2D and we
do not use polynomial fitting to avoid meaningless results outside of the
fitting range. We emphasize that the bias and uncertainty are modeled to be
functions of the measured properties, not intrinsic ones, such that we can
correct for bias based on our measurements.

For LSBGs in our sample, we apply corrections for the bias using the
interpolated bias terms. We first correct for the bias in size, $g$-band
average surface brightness, and $g-r$ and $g-i$ colors. Then we calculate the
$g$-band total magnitude following $m_{g}=\overline{\mu}_{\rm
eff}(g)-2.5\log(2\pi r_{e}^{2})$. The magnitudes and surface brightnesses in
other bands are derived using $g$-band magnitude, surface brightness, and
colors. In this way, we apply a self-consistent correction for the measurement
biases to the data. The measurement uncertainty consists of a statistical
uncertainty, which is determined by the shape of the likelihood (posterior)
surface, and a systematic uncertainty, which is related to various factors
including sky subtraction, neighbor contamination, etc. Unlike other
parametric modeling codes such as imfit (Erwin, 2015) and the Tractor (Lang et
al., 2016), scarlet does not explore the full posterior space but rather finds
one optimal solution. Thus we have no access to the statistical error on the
measured properties from scarlet. Fortunately, by comparing the recovered
properties of mock galaxies with the truth, we can empirically estimate the
measurement uncertainty without knowing the impact of each factor. Following
the same method as that of calculating the bias, in each bin we compute the
$1\sigma$ standard deviation of the difference between the truth and the bias-
corrected measurement, then we interpolate over the grid. The measurement
errors $\sigma(X)$ are shown in Figure 4 as contours, and they have the same
units as the biases. We set minimum uncertainties to be
$\sigma(r_{e})\geqslant 0.3,\ \sigma(\overline{\mu}_{\rm eff})\geqslant 0.05,\
\sigma(g-i)\geqslant 0.05$ to avoid meaningless uncertainty due to small
statistics.