# The Prevalence of Galaxy Overdensities Around UV-Luminous Lyman
$\mathbf{\alpha}$ Emitters in the Epoch of Reionization

E. Leonova1,2,3, P. A. Oesch1,4, Y. Qin5,6,7, R. P. Naidu8, J. S. B.
Wyithe5,6, S. de Barros1, R. J. Bouwens9,R. S. Ellis10, R. M. Endsley11, A.
Hutter12, G. D. Illingworth13, J. Kerutt1, I. Labbé14, N. Laporte15,16,D.
Magee13, S. J. Mutch5,6, G. W. Roberts-Borsani17, R. Smit18, D. P. Stark11, M.
Stefanon9,S. Tacchella19, A. Zitrin20
1Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290
Versoix, Switzerland
2GRAPPA, Anton Pannekoek Institute for Astronomy and Institute of High-Energy
Physics,
University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
3Nikhef, Science Park 105, 1098 XG Amsterdam, The Netherlands
4Cosmic Dawn Center (DAWN), Niels Bohr Institute, University of Copenhagen,
Jagtvej 128, København N, DK-2200, Denmark
5School of Physics, University of Melbourne, Parkville, VIC 3010, Australia
6ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D)
7Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy
8Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA
9Leiden Observatory, Leiden University, PO Box 9500, 2300 RA Leiden, The
Netherlands
10Department of Physics and Astronomy, University College London, Gower
Street, London WC1E 6BT, UK
11Steward Observatory, University of Arizona, 933 N Cherry Ave, Tucson, AZ
85721, USA
12Kapteyn Astronomical Institute, University of Groningen, PO Box 800, NL-9700
AV Groningen, the Netherlands
13Department of Astronomy and Astrophysics, UCO/Lick Observatory, University
of California, Santa Cruz, CA 95064, USA
14Centre for Astrophysics & Supercomputing, Swinburne University of
Technology, PO Box 218, Hawthorn, VIC 3112, Australia
15Kavli Institute for Cosmology, University of Cambridge, Madingley Road,
Cambridge CB3 0HA, UK
16Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue,
Cambridge CB3 0HE, UK
17Department of Physics and Astronomy, University of California, Los Angeles,
430 Portola Plaza, Los Angeles, CA 90095, USA
18Astrophysics Research Institute, Liverpool John Moores University, 146
Brownlow Hill, Liverpool L3 5RF, UK
19Department of Physics, Ulsan National Institute of Science and Technology
(UNIST), Ulsan 44919, Republic of Korea
20Physics Department, Ben-Gurion University of the Negev, P.O. Box 653, Be’er-
Sheva 84105, Israel E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Before the end of the epoch of reionization, the Hydrogen in the Universe was
predominantly neutral. This leads to a strong attenuation of Ly$\alpha$ lines
of $z\gtrsim 6$ galaxies in the intergalactic medium. Nevertheless, Ly$\alpha$
has been detected up to very high redshifts ($z\sim 9$) for several especially
UV luminous galaxies. Here, we test to what extent the galaxy’s local
environment might impact the Ly$\alpha$ transmission of such sources. We
present an analysis of dedicated Hubble Space Telescope (HST) imaging in the
CANDELS/EGS field to search for fainter neighbours around three of the most UV
luminous and most distant spectroscopically confirmed Ly$\alpha$ emitters:
EGS-zs8-1, EGS-zs8-2 and EGSY-z8p7 at $z_{\mathrm{spec}}=7.73$, 7.48, and
8.68, respectively. We combine the multi-wavelength HST imaging with Spitzer
data to reliably select $z\sim 7-9$ galaxies around the central, UV-luminous
sources. In all cases, we find a clear enhancement of neighbouring galaxies
compared to the expected number in a blank field (by a factor $\sim
3-9\times$). Our analysis thus reveals ubiquitous overdensities around
luminous Ly$\alpha$ emitting sources in the heart of the cosmic reionization
epoch. We show that our results are in excellent agreement with expectations
from the Dragons simulation, confirming the theoretical prediction that the
first ionized bubbles preferentially formed in overdense regions. JWST follow-
up observations of the neighbouring galaxies identified here will be needed to
confirm their physical association and to map out the ionized regions produced
by these sources.

###### keywords:

galaxies: high-redshift, galaxies: abundances, galaxies: formation, cosmology:
reionization, galaxies: groups: general

††pubyear: 2021††pagerange: The Prevalence of Galaxy Overdensities Around UV-
Luminous Lyman $\mathbf{\alpha}$ Emitters in the Epoch of Reionization–A

##  1 Introduction

The Epoch of Reionization (EoR) remains one of the key frontiers of
extragalactic studies. Understanding exactly when and how reionization
occurred within the first Gyr of cosmic history is a major goal of astronomy
over the next decade. While great progress has been made to constrain the
overall time frame of reionization through different measurements, including
the polarization of the Cosmic Microwave Background (CMB; Planck Collaboration
et al., 2020), the onset and duration of the EoR are still uncertain (e.g.,
Greig & Mesinger 2017; Naidu et al. 2020; Bosman et al. 2021; also see
Robertson 2021 for a recent review).

An efficient tool used to further constrain the evolution of the neutral
fraction at $z>6$ is the Ly$\alpha$ line emitted by galaxies (so-called
Ly$\alpha$ Emitters; LAEs). Due to its resonant nature, Ly$\alpha$ is easily
scattered and absorbed in the neutral hydrogen in the IGM making LAEs a
sensitive probe of the progress of cosmic reionization (e.g., Miralda-Escudé,
1998; McQuinn et al., 2007; Dayal et al., 2011; Mason et al., 2018). Indeed,
the fraction of continuum selected galaxies that do show Ly$\alpha$ has been
shown to decrease rapidly when entering the neutral phase of the Universe at
$z>6$ (e.g., Stark et al. 2010; Fontana et al. 2010; Schenker et al. 2012;
Treu et al. 2013; Pentericci et al. 2014; Hoag et al. 2019; also see Ouchi et
al. 2020 for a recent review).

In stark contrast to this overall decrease of LAE fractions in the EoR are the
successful detections and spectroscopic confirmations of particularly
luminous, continuum-selected galaxies at $z>7$. From a UV-luminous sample with
$M_{\mathrm{UV}}\sim-22$ mag identified in the CANDELS fields (Roberts-Borsani
et al., 2016), all four galaxies were subsequently confirmed spectroscopically
through their Ly$\alpha$ line (Oesch et al., 2015; Zitrin et al., 2015;
Roberts-Borsani et al., 2016; Stark et al., 2017; Pentericci et al., 2018),
and in some cases also through other UV emission lines. These galaxies thus
first hinted at a trend of larger Ly$\alpha$ fractions in more luminous
galaxies at $z>7$ compared to fainter sources (see also Ono, 2012; Mason et
al., 2018). This has subsequently been confirmed with larger statistics (e.g.,
Jung et al., 2020; Endsley et al., 2021b; Laporte et al., 2021), albeit not
all luminous sources at $z>7$ revealed detectable Ly$\alpha$.

Several possible explanations for larger detection rates of Ly$\alpha$ in
luminous sources have been discussed in the literature. In particular, the
sample from Roberts-Borsani et al. (2016) was initially selected through an
IRAC excess, which originates from extremely strong [OIII]+H$\beta$ emission
lines, with rest-frame equivalent widths in excess of 800 Å. Such strong lines
seem to be typical at these redshifts (e.g., Labbé et al., 2013; Smit et al.,
2015; De Barros et al., 2019; Endsley et al., 2021a). The strong radiation
field required to power such extreme lines might also impact the local
environment around these galaxies, increasing the transmission of Ly$\alpha$
(Stark et al., 2017; Endsley et al., 2021b). An additional contribution of
ionizing photons could come from an active galactic nucleus (AGN), which might
be present in some luminous sources (Tilvi et al., 2016; Laporte et al., 2017;
Mainali et al., 2018; Endsley et al., 2021b). Additionally, luminous sources
have been observed to show significant velocity offsets between their
Ly$\alpha$ lines and the systemic redshifts, likely due to strong outflows
(e.g., Erb et al., 2014; Willott et al., 2015; Stark et al., 2017; Hashimoto
et al., 2019; Matthee et al., 2020). This would shift the red wing of the
Ly$\alpha$ lines further away from the IGM damping wing, further increasing
transmission. However, it is still not clear whether this effect is large
enough to explain the observed luminosity-dependent evolution of the
Ly$\alpha$ fraction over $z=6-8$ (e.g., Mason et al., 2018).

Another possible explanation for a very high Ly$\alpha$ fraction in luminous
sources could be that they are embedded in the most overdense regions. These
overdense volumes are expected to reionize first thanks to the combined
emission from both the luminous galaxies and the numerous fainter sources
around them (e.g., Barkana & Loeb, 2004; Dayal et al., 2009; Pentericci et
al., 2011; Endsley et al., 2021b). Indeed, the bright galaxies themselves are
most likely unable to emit enough ionizing photons to produce an ionized
bubble large enough to allow Ly$\alpha$ photons to escape (Wyithe & Loeb,
2005).

Observationally, several galaxy overdensities have now been identified during
the EoR both for continuum and narrow-band selected galaxies (e.g., Trenti et
al., 2012; Ishigaki et al., 2016; Hu et al., 2021). In several cases, there is
strong evidence for a connection between galaxy overdensities and overlapping
ionized bubbles that allow Ly$\alpha$ to be transmitted (e.g., Castellano et
al., 2016; Tilvi et al., 2020; Endsley et al., 2021b). In this scenario,
reionization starts in the first overdense regions, implying large spatial
fluctuations in the ionized fraction and an inhomogeneous topology. This is
indeed what is found in cosmological simulations (e.g., Iliev et al., 2006;
Hutter et al., 2017; Geil et al., 2017; Kulkarni et al., 2019; Gronke et al.,
2020; Ocvirk et al., 2020; Qin et al., 2021; Smith et al., 2021).

In this paper, we extend the search for overdensities and possible ionized
bubbles to very high redshifts. In particular, we obtained dedicated HST
imaging around three of the most distant and UV-luminous Ly$\alpha$-detected
galaxies, in order to search for fainter neighbouring galaxies. In total, we
study the environment of four galaxies in the EGS field: three confirmed
Ly$\alpha$ emitting sources at $z=7.5-8.7$, in addition to a photometrically
selected galaxy at $z\sim 9$ (see Sect. 2).

The paper is organized as follows: in Section 2, we present the HST and
ancillary datasets used. Section 3 focuses on the selection of our galaxy
sample, before we quantify the environment around the UV luminous targets in
Section 4 and compare our results with simulations in Section 5. Finally, we
discuss implications for the Lyman-alpha visibility of luminous reionization-
era galaxies in Section 6.

Throughout this paper, we adopt a standard cosmology with $\Omega_{m}=0.3$,
$\Omega_{\Lambda}=0.7$, $\mathrm{h}=0.7$. Magnitudes are given in the AB
system.

Figure 1: Position of all selected sources in the CANDLES/EGS field centred
around our target galaxies. Sources with redshift $7<z_{\mathrm{phot}}<8$ are
shown in the left panel, while galaxies with $z_{\mathrm{phot}}>8$ are
indicated on the right. The fields are shown North up, East left. The light
shaded region shows the CANDELS $H_{160}$ coverage. The darker grey areas show
fields where F105W imaging is available and that have been included in our
search. The three pointings from our own program were centred on EGS-zs8-1
(orange outline), EGS-zs8-2 (yellow) and EGSY-z8p7 (pink). Also highlighted
(in green) is the F105W outline around EGS910-3 presented in Bouwens et al.
(2016) that we include in our analysis. Black squares correspond to the
location of Ly$\alpha$ emitters. The colourbars indicate the photometric
redshifts of the sources. The sizes of the circles represent their $H_{160}$
magnitude. Also indicated is the 3D distance of 9 physical Mpc between the two
Ly$\alpha$ galaxies at $z\sim 7.6$ on the left.

##  2 Imaging data

###  2.1 Ancillary Data

The EGS field is one of the five extragalactic fields observed with the Cosmic
Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS) program
(Koekemoer et al., 2011; Grogin et al., 2011). It is centred at (RA,
DEC)=(214.825000, + 52.825000) and was originally covered with 2-orbit
exposures of Advanced Camera for Surveys (ACS) $V_{606}$, $I_{814}$ and Wide
Field Camera 3 (WFC3/IR) $J_{125}$ and $H_{160}$ imaging (0.67 orbits in
F125W, and 1.33 in F160W). The 3D-HST survey covers the field with $JH_{140}$
(F140W) imaging (Skelton et al., 2014; Momcheva et al., 2016), and a few
individual follow-up programs have added targeted $Y_{105}$ (F105W) images for
a few sources (e.g., Bouwens et al., 2016). Additionally, the field has
extensive multi-wavelength data from ground-based facilities and was covered
with Spitzer’s Infrared Array Camera (IRAC) as part of the S-CANDELS program
(Ashby et al., 2015). The total observed area of the EGS field is $\sim$200
arcmin2. For more information on the ancillary multi-wavelength data in this
field see, e.g., Skelton et al. (2014) and Stefanon et al. (2017).

###  2.2 Dedicated Follow-up of Luminous High-Redshift Sources

The EGS field is especially rich in early galaxies (e.g., Bouwens et al.,
2015). It has revealed several particularly luminous high-redshift candidates
at $z>7$ (e.g. Bouwens et al., 2015, 2019; Roberts-Borsani et al., 2016;
Finkelstein et al., 2021), three of which have subsequently been confirmed
through Ly$\alpha$ emission lines. These are EGS-zs8-1 at
$z_{\mathrm{spec}}=7.7302\pm 0.0006$ (Oesch et al. 2015; see also Stark et al.
2017; Tilvi et al. 2020), EGS-zs8-2 at $z_{\mathrm{spec}}=7.4770\pm 0.0008$
(Roberts-Borsani et al., 2016; Stark et al., 2017), and EGSY-z8p7 at
$z_{\mathrm{spec}}=8.683^{+0.001}_{-0.004}$ (Zitrin et al. 2015; see also
Mainali et al. 2018). All of these sources are expected to exhibit strong
radiation fields due to high-equivalent width [OIII]+H$\beta$ emission lines
identified through their Spitzer/IRAC photometry Roberts-Borsani et al.
(2016). Indeed, highly excited UV lines have been confirmed in the MOSFIRE
spectra of EGS-zs8-1 (C iii; Stark et al. 2017) and EGSY-z8p7 (N v; Mainali et
al. 2018). Especially the detection of the N v line in EGSY-z8p7 points to a
non-negligible, possible contribution from an AGN.

While other particularly distant EoR galaxies have been spectroscopically
confirmed since the discovery of these sources in the EGS field (e.g. Oesch et
al., 2016; Jiang et al., 2021; Hashimoto et al., 2018; Jung et al., 2020;
Laporte et al., 2021), they still rank among the most distant detected
Ly$\alpha$ emission lines (see also Larson et al., in prep; Larson et al.,
2020). For more information on the luminous Ly$\alpha$ sources, see Section
3.1 and Tables 1 and 2.

To explore whether these very early Lyman-alpha emitters commonly reside in
overdense regions, we have obtained dedicated, additional multi-colour HST
WFC3/IR observations around these three luminous galaxies (GO-15103; De
Barros, 2017). In particular, we added the missing F105W filter that is
necessary to reliably identify $z\sim 7.5-9$ galaxies and measure improved
redshifts, in addition to deeper F125W and F160W imaging. The program was
designed to identify galaxies that are up to an order of magnitude less
luminous than the main targets, which have $H_{AB}=25.0-25.3$ mag. In
particular, we obtained 5-orbit exposures in three individual pointings
centred on each of the three luminous LAEs (2 orbits in F105W, 1.5 in F125W,
and 1.5 in F160W). When combined with the previous, CANDELS HST imaging our
data reach 5$\sigma$ limits of 27.9 mag in F105W, 27.5 mag in F125W, and 27.5
mag in F160W, respectively.

In addition to our dedicated HST data, we retrieved all the HST ACS and
WFC3/IR data available around the EGS sources from the HST archive and
combined them into one mosaic as part of the Hubble Legacy Field project
(Illingworth et al, in prep).

###  2.3 Source Detection and Multi-Wavelength Photometry

The HST images in different bands were convolved to the same point-spread
function (PSF; using convolution kernels analogous to the ones made available
by the 3D-HST team; Skelton et al. 2014). We then use SExtractor (Bertin &
Arnouts, 1996) to detect sources in a combined F125W+F160W image and measure
their colours through matched, circular apertures. In particular, we use
apertures with 0$\aas@@fstack{\prime\prime}$5 diameter to measure colours,
which are then corrected to total fluxes using the SExtractor AUTO fluxes in
the F160W image.

In order to obtain IRAC photometry in all four channels, we follow the same
procedure as described in Stefanon et al. (2021a). To overcome confusion, we
subtract contaminating flux of neighbouring galaxies using the Mophongo tool
(see, e.g., Labbé et al., 2013, 2015), before measuring the IRAC fluxes in
1$\aas@@fstack{\prime\prime}$8 diameter apertures and correcting to total
fluxes. When the contamination by neighbour is too large for a reliable
subtraction, the IRAC fluxes are flagged and ignored in the subsequent
spectral energy distribution fits. We use the same procedure to derive
photometry from ground-based K-band images. Overall, we thus obtain 11-band
spectral energy distributions (SED) for all sources (6$\times$ HST, K-band,
and 4$\times$ IRAC), spanning 0.6 to 8 µm.

##  3 Sample Selection

###  3.1 Central LAE Targets

Our HST follow-up observations are centred around three UV luminous galaxies
at $z\geq 7.5$ that have been confirmed through Ly$\alpha$ in the EGS field
(see Section 2.2). One of these, EGS-zs8-1, has subsequently been observed
through narrow-band observations, from which two fainter LAEs have been
identified in its immediate proximity (Tilvi et al., 2020). These latter
sources have also been confirmed through Ly$\alpha$ spectroscopy and they were
named z8_SM at $z=7.767$ and z8_4 at $z=7.748$. They likely sit in a common,
ionized bubble with EGS-zs8-1, since they are only separated by 0.7 to 0.8
physical Mpc (pMpc) from each other and their estimated ionized bubble radii
reach $0.5-1.0$ pMpc.

Additionally, the galaxy group around EGS-zs8-1 is separated by only 9.1 pMpc
from the other UV-luminous, spectroscopically confirmed source EGS-zs8-2 at
z=7.48. This distance is smaller than the average bubble radius of
$\lesssim$10 pMpc expected at the end of cosmic reionization (Wyithe & Loeb,
2005; Geil et al., 2017), and points at a larger scale overdensity in this
field at this redshift. One of the main goals of our follow-up observations
was to further characterize this overdensity through fainter sources. Hence,
we centred two of our HST pointings around EGS-zs8-1 and EGS-zs8-2,
respectively (orange and yellow outlines in Fig 1).

The last pointing was centred around EGSY-z8p7, one of the most distant known
Ly$\alpha$ emitters with $z_{\mathrm{spec}}=8.683$ (pink outline in Fig 1).
This source has a potential companion at $<10$ pMpc: the galaxy EGS910-3
(Bouwens et al. 2016, see also Bouwens et al. 2019). On the sky, these two
sources are separated by only 3.7 arcmin. However, the source only has a
photometric redshift measurement ($z_{\mathrm{phot}}=9.0^{+0.5}_{-0.7}$ as
derived in Bouwens et al. 2016) and has not been confirmed through
spectroscopy so far. Hence, the actual physical separation from EGSY-z8p7 is
uncertain. The source EGS910-3 has also been targeted by F105W follow-up
observations, which we retrieve from the archive and include in our analysis
(see, e.g., Fig 6).

###  3.2 Selection of Neighbouring z=7-9 Candidates

####  3.2.1 Selection Criteria

The main goal of our program was to identify fainter sources around the
central UV-luminous targets. We thus use the HST-detected multi-wavelength
catalogs described in Section 2.3 to select candidate high-redshift sources.
Specifically, we use a combination of detection and non-detection criteria as
expected for Lyman break galaxies together with a photometric redshift
selection. We only consider sources that are significantly detected in both
the F125W and F160W bands with the following criteria:

$S/N(J_{125})>3\land(S/N)_{\mathrm{det}}>5\land S/N(H_{160})>4,$ (1)

where $(S/N)_{\mathrm{det}}$ corresponds to the signal-to-noise in the
combined F125W and F160W image. These criteria allow us to identify sources
down to $H_{160}=27.7$ mag. Additionally, the sources were required to not be
detected in the two optical HST filters with: $S/N(V_{606})<2\land
S/N(I_{814})<2$. In selecting these sources, we have only included areas that
also have F105W imaging coverage, either from our own follow-up program or
from archival observations (darker areas in Fig 1).

Photometric redshifts are then derived for all the sources that pass the above
criteria. We use the python version of the publicly available code
EAZY111https://github.com/gbrammer/eazy-py (Brammer et al., 2008) which
provides the photometric redshift probability distribution function (pdf),
$\mathcal{P}(z)$, based on a $\chi^{2}$ fit of SED models to the observed
photometry. Specifically, we use the template set ‘eazy_v1.1’ which includes
emission lines. We have also tested other template sets, but have not found
large differences. EAZY also accounts for the intergalactic absorption by
neutral hydrogen (Madau, 1995). The redshifts have been searched over $0.1\leq
z\leq 12$.

Based on the EAZY output we then only include sources that have a best-fit
photometric redshift in our range of interest, depending on the central target
sources. For the fields around EGS-zs8-1 and -2, we select galaxies with
$z_{\mathrm{phot}}=7-8$. To ensure that these redshifts are reasonably
accurate, we exclude galaxies with a photometric redshift pdf that is too
wide, i.e., galaxies are required to have $\mathcal{P}(6\leq z\leq 9)\geq
50\%$. For the fields around EGSY-z8p7, we select sources with
$z_{\mathrm{phot}}>8$ that have $\mathcal{P}(7\leq z\leq 10)\geq 50\%$. In
Figure 6, we illustrate the SEDs and redshift pdfs of the central UV-luminous
Ly$\alpha$ emitters and EGS910-3. As can be seen, the photometric redshift
determinations of the Ly$\alpha$ sources are in excellent agreement with their
spectroscopic redshifts, with a very narrow $\mathcal{P}(z)$. For the last
source, we confirm the photometric redshift solution from Bouwens et al.
(2016), as we measure a consistent value of
$z_{\mathrm{phot}}=9.18^{+0.95}_{-0.55}$.

####  3.2.2 Exclusion of Dwarf Stars

Cool dwarf stars can mimic the colours of high-redshift Lyman break galaxies
(e.g. Wilkins et al., 2014). They can remain undetected in optical bands (in
F606W and F814W filters) and peak in the Near-infrared (NIR). To remove such
sources from our catalogs of LBG candidates we use a combination of two
different measurements. First, we identify potential unresolved sources using
the SExtractor parameter $\mathrm{CLASS\\_STAR}$, and second we run EAZY with
a set of dwarf star templates fixed at $z=0$ to compare the $\chi^{2}$ from
both dwarf stars and galaxy templates. The $\mathrm{CLASS\\_STAR}$ parameter
on its own is not reliable, especially at faint magnitudes.

Objects with $\mathrm{CLASS\\_STAR}$ more than 0.6 and $\chi^{2}$ from star
templates less than $\chi^{2}$ from galaxy template are removed from the
catalogs:

$\centering\mathrm{CLASS\\_STAR}>0.6\land(\chi^{2}_{\mathrm{galaxy}}/\chi^{2}_{\mathrm{star}})>1\@add@centering$
(2)

In total, only two sources have been removed based on these criteria. Visual
inspection of these sources indeed confirmed them to be unresolved, and hence
more likely to be stars than high-redshift galaxies.

Table 1: Table with parameters of all sources with redshift $z\sim 7-8$
identified in the F105W footprints around the two LAEs.

ID | RA | DEC | H160 [mag] | zphot | $\int\mathcal{P}(6\leq z\leq 9)$ | zspec | Separationa [′] | References
---|---|---|---|---|---|---|---|---
EGS-zs8-2 group
EGS-zs8-2 | 215.0503352 | 53.0074442 | $25.25\pm 0.03$ | $7.43^{+0.08}_{-0.09}$ | 1.00 | $7.477$ | 0.0 | (1),(3)
EGSY-94562 | 215.0394433 | 53.0207175 | $27.2\pm 0.3$ | $7.90^{+0.04}_{-1.90}$ | 0.55 | $-$ | 0.9 |
EGSY-80113 | 215.056893 | 53.0110575 | $27.4\pm 0.2$ | $7.22^{+0.02}_{-1.22}$ | 0.50 | $-$ | 0.3 |
EGSY-76521 | 215.0697905 | 53.0151445 | $27.5\pm 0.2$ | $7.12^{+0.20}_{-0.93}$ | 0.90 | $-$ | 0.8 |
EGSY-91929 | 215.0269545 | 53.0077343 | $27.7\pm 0.2$ | $7.96^{+0.20}_{-0.95}$ | 0.93 | $-$ | 0.8 |
EGS-zs8-1 group
EGS-zs8-1 | 215.1453377 | 53.0042107 | $25.12\pm 0.04$ | $7.88^{+0.06}_{-0.05}$ | 1.00 | 7.728 | 0.0 | (1),(2),(3),(4)
EGSY-23959 | 215.1508931 | 52.9895618 | $26.1\pm 0.1$ | $7.59^{+0.13}_{-0.36}$ | 1.00 | $-$ | 0.9 |
EGSY-46810 | 215.154764 | 53.0289543 | $26.8\pm 0.1$ | $7.67^{+0.23}_{-1.67}$ | 0.66 | $-$ | 1.5 |
EGSY-46009 | 215.1168772 | 53.0010665 | $26.9\pm 0.2$ | $7.22^{+0.10}_{-1.16}$ | 0.86 | $-$ | 1.0 |
z8_SM | 215.1487577 | 53.0025984 | $27.2\pm 0.1$ | $7.89^{+0.16}_{-0.29}$ | 1.00 | 7.767 | 0.2 | (4)
EGSY-44464 | 215.1314765 | 53.0092156 | $27.4\pm 0.2$ | $7.21^{+0.09}_{-1.21}$ | 0.76 | $-$ | 0.6 |
EGSY-31530 | 215.1334993 | 52.9899563 | $27.5\pm 0.2$ | $7.86^{+0.41}_{-1.86}$ | 0.51 | $-$ | 1.0 |
EGSY-21312 | 215.150041 | 52.9846593 | $27.7\pm 0.2$ | $7.85^{+0.24}_{-0.70}$ | 0.96 | $-$ | 1.2 |
z8_4∗ | 215.1465376 | 52.994614 | $26.8\pm 0.3$ | $6.16^{+0.90}_{-0.48}$ | $-$ | 7.748 | 0.6 | (4)

Notes: a Separation from the central luminous source in arcmin.

∗z8_4 was not selected in our sample due to a lower photometric redshift.

References: (1) Roberts-Borsani et al. (2016), (2) Oesch et al. (2015), (3)
Stark et al. (2017), (4) Tilvi et al. (2020)

Table 2: Table with parameters of all sources with redshift $z\sim 8-9$
identified in the F105W footprints around the two $z>8$ galaxies.

ID | RA | DEC | H160 [mag] | zphot | $\int\mathcal{P}(7\leq z\leq 10)$ | zspec | Separationa [′] | References
---|---|---|---|---|---|---|---|---
EGSY-z8p7 group
EGSY-z8p7 | 215.0354269 | 52.8906918 | $25.19\pm 0.03$ | $8.55^{+0.14}_{-0.15}$ | 1.00 | $8.683$ | 0.0 | (1),(2),(3),(4)
EGSY-6887∗ | 215.0465028 | 52.8895178 | $25.8\pm 0.1$ | $8.05^{+0.22}_{-0.54}$ | 0.95 | $-$ | 0.4 |
EGSY-12580† | 215.0373571 | 52.8918284 | $26.4\pm 0.2$ | $8.07^{+0.19}_{-0.29}$ | 0.99 | $-$ | 0.1 |
EGSY-27727 | 215.031966 | 52.9122172 | $27.3\pm 0.2$ | $8.15^{+0.79}_{-1.78}$ | 0.68 | $-$ | 1.3 |
EGS910-3 group
EGS910-3 | 214.9386673 | 52.91176 | $26.8\pm 0.2$ | $9.18^{+0.95}_{-0.55}$ | 0.63 | $-$ | 0.0 | (2),(3)
EGSY-66990 | 214.968693 | 52.9296537 | $26.4\pm 0.2$ | $9.12^{+0.14}_{-1.08}$ | 0.95 | $-$ | 1.5 |
EGSY-69674 | 214.9675507 | 52.932962 | $26.6\pm 0.1$ | $8.60^{+0.25}_{-0.44}$ | 0.99 | $-$ | 1.6 | (3),(4)
EGSY-56680 | 214.9699133 | 52.9140338 | $27.2\pm 0.1$ | $8.60^{+0.14}_{-2.60}$ | 0.66 | $-$ | 1.1 |

Notes: a Separation from the central luminous source in arcmin.

∗ This source lies very close to a clumpy foreground galaxy, which could
affect its photometric measurements.

† This source has a close companion, which was excluded from our list due to
its photometric redshift just below the cutoff of our selection.

References: (1) Zitrin et al. (2015), (2) Bouwens et al. (2016), (3) Bouwens
et al. (2019), (4) Finkelstein et al. (2021)

####  3.2.3 Final Sample

After excluding two likely dwarf stars, our final sample of high-redshift
galaxies within the four WFC3/IR pointings around the four UV-luminous central
targets consists of 21 sources in total. Their positions are shown in Fig. 1.
In the following we will discuss the individual environment around the four
central sources.

With 8 sources, the pointing around EGS-zs8-1 revealed the largest number of
high-redshift galaxy candidates at $z=7-8$ (including the central source
itself). The EGS-zs8-2 field resulted in 5 sources. Note that the
spectroscopically confirmed source z8_4 is not included in our photometric
sample due to a photometric redshift that was too low. For completeness, we
show the galaxy in Fig. 1.

At higher redshifts, the fields around EGSY-z8p7 and EGS910-3 reveal 4 sources
with $z_{\mathrm{phot}}>8$ each. All these sources are listed in Tables 1 and
2. The photometric redshift distribution functions and stamps for the luminous
sources can be found in the appendix in Fig 6, while the fainter, neighbouring
sources are shown in Fig. 7 to 10.

Figure 2: The cumulative number of the sources within an area of 4.5 arcmin2
at $z=7-8$ (left) and $z>8$ (right). The left panels show the cumulative
numbers as a function of absolute magnitude. The blue line corresponds to the
expected number based on the UV LF and our selection function, while the
coloured lines correspond to the actual numbers of galaxies that we have
identified in the HST data in the different fields (orange for EGS-zs8-1;
yellow for EGS-zs8-2; pink for EGSY-8p7; green for EGS910-3). The shaded
regions correspond to the 16th to 84th percentile uncertainties for a Poisson
distribution convolved with cosmic variance (see text for details). The right
panels show the total number of sources selected down to our magnitude limits
using the same colours as for the cumulative histograms,with the points offset
horizontally for clarity. Additionally, the purple histograms show the
distribution of high-redshift galaxies selected in 82 independent 4.5 arcmin2
cells in the GOODS-S and GOODS-N fields, using the same selection criteria as
for our HST data in the EGS field (see section 4.1.2). The dashed black line
corresponds to the mean in the GOODS fields, which is consistent with the mean
expected number from the UV LF (blue). As is evident from these figures, the
environment around all four UV luminous sources are overdense relative to the
general field. The enhancement ranges from 3$\times$ to 9$\times$.

##  4 The Environment of UV Luminous Galaxies

A principle goal of our work is to quantify the environment and possible
overdensities of UV luminous sources with Ly$\alpha$ emission at $z>7$. To do
this, we first need to estimate the expected number of galaxies in a given
area, which we do based on two different approaches that are described in the
following subsections.

###  4.1 Expected Number of Galaxies

####  4.1.1 Based on the UV Luminosity Function

The first approach to estimate how many galaxies we expect to find within a
given area is based on the UV luminosity function (LF) and the galaxy
selection efficiency. Following previous literature (e.g. Castellano et al.,
2016), we estimate the selection efficiency of our survey by simulating the
photometry for a large number of high-redshift galaxies. Specifically, we
introduce $>60,000$ artificial galaxies with a range of known (input)
magnitudes $24-29$ mag, colour distribution, and redshifts $6<z<10$.

Since high-redshift galaxy selections can be sensitive to the UV colours, we
simulate the artificial sample with a colour distribution that matches
observations. To do this, we start with a Bruzual & Charlot (2003) template
with an age of 100 Myr and a constant star-formation history. Calzetti et al.
(2000) dust is then applied following a normal distribution of E(B-V) with a
mean of 0.1 and a $\sigma$ = 0.15, truncated at zero (following Finkelstein et
al., 2021). Nebular continuum and emission lines have been added empirically
based on the total number of ionizing photons produced by the underlying SED
and using line ratios from Anders & Fritze-v. Alvensleben (2003). Finally, IGM
attenuation is applied according to the prescription of Madau (1995) and set
by the simulated redshift of each galaxy. For $z>7$, this essentially
corresponds to a sharp break at the redshifted Ly$\alpha$ line. The final SEDs
are then convolved with all the HST and Spitzer IRAC transmission curves to
derive the expected, intrinsic fluxes. These are subsequently disturbed using
Gaussian flux uncertainties, based on the median flux $\sigma$ as measured for
real sources in our HST footprints in the respective filters.

We then perform the same analysis as for the real catalog. First, the
detection and optical non-detection criteria are applied and the EAZY
photometric redshifts and pdfs are computed as was done for the real galaxies.
This allows us to compute the completeness in bins of input magnitude and
redshift as the ratio between the number of galaxies that pass our selection
and the number of input galaxies in the respective bin:

${\cal{P}}(M,z)=\frac{N_{\rm{sel}}}{N_{\rm{in}}}.$ (3)

Now that we have the selection efficiency, it is straightforward to compute
the expected number of galaxy detections based on the empirical UV LF. To do
this, we use the Schechter parameter evolution derived in Bouwens et al.
(2021) at $z\sim 2-10$:

$\begin{cases}M_{\rm{UV}}^{*}=-21.03-0.04(z-6)\\\
\phi^{*}/\rm{Mpc}^{-3}=0.4\times 10^{-3}\times
10^{-0.33(z-6)-0.24(z-6)^{2}}\\\ \alpha=-1.94-0.11(z-6)\\\ \end{cases}$ (4)

The expected number of detections can then be derived by multiplying the UV LF
with the probability that a galaxy of absolute magnitude $M$ at redshift $z$
is selected, $\mathcal{P}(M,z)$, within the same area as used for real
sources.

$N_{\mathrm{exp}}=\int\Omega\frac{dV_{\rm{com}}}{dz}dz\int\phi(M)\mathcal{P}(M,z)dM$
(5)

where $\phi(M)$ is the UV LF, $\Omega$ is the survey area for each HST field,
i.e., 4.5 arcmin2 for a WFC3/IR pointing.

For the $z\sim 7-8$ galaxy samples around EGS-zs8-1 and -2, we derive an
expected number of 1.57 galaxies per WFC3/IR pointing down to our detection
limit of 27.7 mag (see Fig. 2). For the $z\sim 8-9$ samples, the expected
number is only 0.45 sources per WFC3/IR field.

####  4.1.2 Empirical Approach

Another way to determine whether a given HST field might be overdense can be
derived empirically using the number counts from another field that is as deep
or deeper than the field in question. We can thus use the CANDELS/GOODS-North
and CANDELS/GOODS-South fields with an area of $\sim$170 arcmin2 each, from
which high redshift galaxy sources have been selected in the past.
Specifically, we use the sample from Bouwens et al. (2015), for which multi-
wavelength catalogs are available from GREATS (Stefanon et al., 2021b). For
consistency, we re-derive the photometric redshifts and pdfs for this sample
using the same EAZY setup as for our main selection. We then select the groups
of sources with best-fit redshifts at $7\leq z\leq 8$ and $8\leq z\leq 9$.
From these, we limit the sources at $H_{160}=27.7$ mag, which corresponds to
the 5$\sigma$ detection limit in our EGS data around the Ly$\alpha$ emitters.

To get a histogram of expected sources in an area equivalent to a single
WFC3/IR pointing, we divide the GOODS-S and GOODS-N fields into circular cells
with radius 1.2 arcmin, corresponding to an area of 4.5 arcmin2 each, and
count the number of the selected sources within each cell. The positioning of
the cells was chosen on a grid allowing for a small ($<10\%$) overlap between
different areas. In total, the GOODS-S field was thus split into 42
independent regions, while the GOODS-N field was divided into 40 separate
cells.

The number count histograms of these 82 regions are also shown in Fig. 2. The
average number of sources per cell is 1.2 at $z\sim 7-8$ and 0.45 at $z\sim
8-9$, which is consistent with the average expected number from the UV LF
derived in the previous section. However, the empirical histograms further
allow us to visualize the effect of cosmic variance. For example, out of 82
independent cells in both GOODS fields, none contain 8 sources, while only 2
contain $\geq$5 $z\sim 7-8$ galaxies. Empirically, finding $\geq$5 galaxies is
thus expected only with a 2.4% probability.

At higher redshift, the largest number of galaxies in a given 4.5 arcmin2
region is 3. Based on this, our finding of 4 (5-8) sources around the UV
luminous LAEs at $z>8$ ($z=7-8$) clearly indicates that these galaxies live in
early overdensities. We quantify this in more detail in the next section.

Figure 3: UV luminosity functions in the EGS field centred around the four UV-
luminous sources (dots) compared to expected average Schechter LFs as derived
in Bouwens et al. (2021, solid lines), and double power-law LFs Bowler et al.
(2020, dashed lines). The purple curve shows the LFs for the $z=7-8$ galaxies,
while the cyan line corresponds to $z=8-9$. This analysis shows that these UV
luminous sources lie in areas with a general density enhancement of
$3-9\times$ that extends to fainter galaxies.

###  4.2 Prevalent Overdensities Around UV Luminous Galaxies

We now ready evaluate the environment of the HST fields around the UV luminous
sources. We start with the two Ly$\alpha$ sources at $z\sim 7.5$, around which
we find 5 and 8 sources in 4.5 arcmin2. This is 3$\times$ and $5\times$
overdense relative to the expected number of 1.57 sources based on the UV LF.

In order to compute more accurate probabilities to find 5 or 8 sources when
1.57 are expected we specifically account for cosmic variance. To do this, we
use the publicly available cosmic variance calculator
galcv222https://github.com/adamtrapp/galcv of Trapp & Furlanetto (2020). For a
4.5 arcmin2 field and a depth of $\Delta z=1$, this predicts a relative
$\sigma_{CV}$ of 32% at $z=7.5$ and 38% at $z=8.5$, down to our magnitude
limits. These numbers are similar (albeit somewhat smaller) than the predicted
cosmic variance from Trenti & Stiavelli (2008).

The full probabilities are then derived from a Monte Carlo simulation in which
we first perturb the expected numbers of galaxies from the UV LF by the
relative cosmic variance and then draw from the Poisson distribution. Doing
this, we find that the probability to find 8 sources, when 1.57 are expected
is indeed small with 0.08%. This corresponds to a 3.2$\sigma$ result.
Similarly, the probability to find $\geq 5$ galaxies is only 3.1%, consistent
with our empirical estimate above. This still represents a significance of
1.9$\sigma$.

For the higher redshift fields, we find 4 sources in each pointing, while only
0.46 were expected on average. This corresponds to an enhancement of
8.7$\times$. Doing the same calculation as above, we find that the probability
of finding 4 (or more) sources is only 0.2%, corresponding to a 2.9$\sigma$
result for each field separately. These calculations show that overdensities
around UV luminous sources have been significantly detected in all four cases,
with overdensity values ranging between $3-9\times$.

This excess can also be seen when looking at the UV LF derived in our target
fields alone. Fig. 3 shows these UV LFs at $z=7-8$ and $z=8-9$ as derived from
the two WFC3/IR pointings centred at the respective redshifts. To compute
this, we use the same selection functions as derived in Section 4.1.1 and
compute the effective volume in a given magnitude bin. Given the small area
(of 9 arcmin2) centred on UV luminous galaxies, the UV LF in these fields is
significantly enhanced in the brightest bins (with $M_{UV}\sim-22$). However,
as can be seen, at the fainter end these UV LFs continue to lie significantly
above the expected field average values. This is consistent with the
hypothesis that these UV luminous galaxies lie in overdensities that are
generally enhanced relative to the rest of the Universe. Our results thus
indicate that the most UV luminous sources can indeed be used to pinpoint the
first overdense regions in the early Universe.

##  5 Expectations from Simulations

In this section we discuss what simulations predict for the environment of UV
luminous galaxies similar to those we have observed.

Specifically, Qin et al. (2021) studied the environment of $z=8$ galaxies in
the Meraxes semi-analytical model (Mutch et al., 2016; Qin et al., 2017) as
part of the Dragons simulation program (Poole et al., 2016).

###  5.1 Number of Neighbouring Galaxies

Qin et al. (2021) predicted that the most luminous galaxies preferentially lie
in overdense regions, in agreement with what we find here from observations.
In particular, Qin et al. (2021) used their simulations to compute the number
of expected sources within a WFC3/IR field-of-view. This can be directly
compared to our observational data.

Figure 4 shows the cumulative number of galaxies expected around the eight
most luminous central sources (having $M_{\mathrm{UV}}<-21.5$) in the
simulations down to the same magnitude limit as used in our observations. As
can be seen, the cumulative number varies between 3 to 13 sources, with a
median of 8.5. This is in excellent agreement with our findings around the two
sources EGS-zs8-1 and EGS-zs8-2. Both the simulations and the observations
thus find a boost in the number counts around luminous sources compared to the
field by $\sim 3-5\times$.

Figure 4: The cumulative number of the sources within an area of 4.5 arcmin2
at $z\sim 8$, both in simulation and in our observations. The left panel shows
the cumulative numbers as a function of absolute magnitude. The coloured lines
correspond to the actual numbers of galaxies that we have identified in the
HST data in different fields (orange for EGS-zs8-1; yellow for EGS-zs8-2). The
grey lines correspond to the numbers from the Dragons simulation (Qin et al.
2021). The right panel shows the total number of sources as a function of the
magnitude of the brightest source selected down to our 27.7 magnitude limit
using the same colours as for the cumulative histograms. The blue line
corresponds to the expected number based on the UV LF and our selection
function (same as Fig 2), with the shaded blue region showing the 16th to 84th
percentile. The horizontal black dashed line shows the average number of
sources in the GOODS fields within the same area. Both the simulations and the
observations thus find a boost of $\sim 3-5\times$ in the number of galaxies
around the most luminous sources.

###  5.2 Ionized Bubble Sizes

The Meraxes simulation further allows us to measure the radii of the ionized
bubbles around galaxies. Qin et al. (2021) show that the most luminous sources
sit in some of the largest ionized bubbles, resulting in a higher total
transmission of the Ly$\alpha$ emission line (see also Wyithe & Loeb, 2005;
Geil et al., 2017; Hutter et al., 2017). The most luminous simulated central
sources lie in ionized regions with measured radii of around 8-9 comoving Mpc,
corresponding to $\sim 1$ physical Mpc at these redshifts (see Fig. 5). We can
thus expect similar ionized region sizes around our observed galaxies. Using
simple estimates, this is indeed what we infer as discussed below.

To convincingly conclude that all of the associated sources we find for our
sample lie within the respective ionised bubbles would require both
spectroscopic confirmation of the association and, ideally, the presence of
Lyman alpha emission. Nonetheless, using the same approach as in Matthee et
al. (2018) we can estimate the sizes of the ionized bubbles that each of our
sources individually could carve out of the neutral IGM. Assuming a spherical
shape of the ionized regions around the galaxies, the expected radii range
from 0.5 pMpc up to 1.1 pMpc. These results are consistent with previous
estimates in the literature for possible bubble sizes at these redshifts (see
e.g., Castellano et al., 2018; Tilvi et al., 2020; Endsley et al., 2021b).

While Tilvi et al. (2020) have confirmed the physical association of two
sources around EGS-zs8-1 through spectroscopy, the new galaxies identified in
this paper only have photometric redshifts available at the moment. The
respective uncertainties in 3D distances are thus too large to determine the
physical associations of these galaxies with the central UV luminous targets.
However, our detection of an enhancement of fainter galaxies around all these
sources points at a larger scale overdensity in the EGS field (see also
Bouwens et al., 2015, 2019; Finkelstein et al., 2021).

Interestingly, the two UV luminous LAEs EGS-zs8-1 and EGS-zs8-2 are only
separated by 9.1 pMpc from each other. It is thus possible that these two
sources lie in a larger, ionized region that spans from $z=7.5$ to $z=7.7$.
Future observations, e.g., with JWST spectroscopy will be needed to confirm
this.

Figure 5: Number of neighboring galaxies brighter than ${\rm H}_{160}{=}27.5$
mag for EGS-zs8-1, EGS-zs8-2 and EGSY-z8p7 (see the star symbols) as a
function of their UV magnitudes. Theoretical expectation using the DRAGONS
simulation (see more in e.g., Fig. 2 of Qin et al. 2021) is shown with circles
for comparison. The size of the corresponding H ii bubbles (in units of cMpc;
indicated by varying circle sizes) is shown with the uncertainties drawn from
500 mock observations.

##  6 Summary and Conclusions

In this paper, we identified faint sources at $z>7$ around three confirmed UV
luminous Lyman $\alpha$ emitters and another, nearby $z_{\mathrm{phot}}\sim 9$
target in the CANDELS/EGS field. We presented dedicated HST data and combined
this with ancillary HST and Spitzer/IRAC imaging in this field.

We find a significant enhancement of fainter galaxies within a WFC3/IR
pointing of 4.5 arcmin2 around each of these UV luminous, central sources. By
comparing the number of detected galaxy candidates with the expected numbers
from the UV LFs at these redshifts and with the numbers in the two GOODS
fields, we estimate these areas to be enhanced by $3-9\times$ compared to the
average field.

Our observational findings are in excellent agreement with the predictions
form the Meraxes simulation (Qin et al., 2021), which also find a boost in
fainter neighbouring sources around UV luminous galaxies. In these
simulations, the combined ionizing photon output of the ensemble of galaxies
around the central sources produce ionized regions with radii of $\sim 1$
pMpc, allowing Ly$\alpha$ photons to be transmitted even when the surrounding
IGM is still highly neutral.

Combined with these simulation results, our observational finding of overdense
regions around these luminous $z=7.5-8.7$ Ly$\alpha$ emitters thus strongly
suggest that the local galaxy environment plays a driving role in the
Ly$\alpha$ transmission. Specifically, the field around EGSY-z8p7 at
$z_{\mathrm{spec}}=8.683$ is enhanced by a factor $\sim 8.7\times$, suggesting
that this source might sit in one of the first overdensities and ionized
bubbles in the Universe, only $\sim 500$ Myr after the Big Bang.

The crucial next step is to obtain spectroscopic redshifts for these
overdensities in order to map out their 3D structure and estimate the size of
the overall ionized regions. The approved JWST cycle 1 program GO-2279 (Naidu
et al., 2021) will achieve exactly this based on NIRCam/grism observations of
EGSY-z8p7 and EGS910-3, which is only separated by 3.7 arcmin. These
observations will provide a sample of [OIII]5007 emission line selected
sources at $z=7-9$, and will thus significantly improve the overdensity
measurement. Other NIRCam/grism and NIRSpec observations over other fields are
poised to identify a large number of ionized bubbles in the EoR in the near
future. This will further clarify the importance of other possible origins for
enhanced Ly$\alpha$ transmission from UV luminous galaxies, such as velocity
offsets or very hard and strong ionization fields, e.g., from an AGN
contribution.

The fact that we find an overdensity around each of the three UV luminous
Ly$\alpha$ emitting galaxies only based on photometric redshifts strongly
suggests that the Ly$\alpha$ transmission is intimately connected to the
overall galaxy environment that helps to create ionized regions. Turning this
around, these UV luminous, early Ly$\alpha$ emitters can be used to pinpoint
the earliest overdense regions in the Universe, where star-formation and
reionization first started. This will become particularly interesting with the
availability of larger samples of UV luminous LBGs and LAEs from future wide-
area surveys with the Euclid or Roman Space Telescopes and follow-up
spectroscopy.

## Acknowledgements

We acknowledge support from: the Swiss National Science Foundation through the
SNSF Professorship grant 190079 (PAO, JK). The Cosmic Dawn Center (DAWN) is
funded by the Danish National Research Foundation under grant No. 140. This
research was supported by the Australian Research Council Centre of Excellence
for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number
CE170100013. Part of this work was performed on the OzSTAR national facility
at Swinburne University of Technology, which is funded by Swinburne University
of Technology and the National Collaborative Research Infrastructure Strategy
(NCRIS). ST is supported by the 2021 Research Fund 1.210134.01 of UNIST (Ulsan
National Institute of Science & Technology). RSE acknowledges funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (grant agreement No 669253). RJB and MS
acknowledge support from TOP grant TOP1.16.057. RS acknowledges support from a
STFC Ernest Rutherford Fellowship (ST/S004831/1). AH acknowledges support from
the European Research Council’s starting grant ERC StG-717001 (“DELPHI”). AZ
acknowledges support by Grant No. 2020750 from the United States-Israel
Binational Science Foundation (BSF) and Grant No. 2109066 from the United
States National Science Foundation (NSF), and by the Ministry of Science &
Technology, Israel.

## Data Availability

The HST data underlying this article will be publicly released as High-Level
Science Product (HLSP) as part of the Hubble Legacy Fields (HLF) project
(Illingworth et al., 2016; Whitaker et al., 2019). A pre-release version can
be shared on reasonable request to the corresponding author.

## References

  * Anders & Fritze-v. Alvensleben (2003) Anders P., Fritze-v. Alvensleben U., 2003, A&A, 401, 1063
  * Ashby et al. (2015) Ashby M. L. N., et al., 2015, ApJS, 218, 33
  * Barkana & Loeb (2004) Barkana R., Loeb A., 2004, ApJ, 609, 474
  * Bertin & Arnouts (1996) Bertin E., Arnouts S., 1996, A&AS, 117, 393
  * Bosman et al. (2021) Bosman S. E. I., et al., 2021, arXiv e-prints, p. arXiv:2108.03699
  * Bouwens et al. (2015) Bouwens R. J., et al., 2015, ApJ, 803, 34
  * Bouwens et al. (2016) Bouwens R. J., et al., 2016, ApJ, 830, 67
  * Bouwens et al. (2019) Bouwens R. J., Stefanon M., Oesch P. A., Illingworth G. D., Nanayakkara T., Roberts-Borsani G., Labbé I., Smit R., 2019, ApJ, 880, 25
  * Bouwens et al. (2021) Bouwens R. J., et al., 2021, arXiv e-prints, p. arXiv:2102.07775
  * Bowler et al. (2020) Bowler R. A. A., Jarvis M. J., Dunlop J. S., McLure R. J., McLeod D. J., Adams N. J., Milvang-Jensen B., McCracken H. J., 2020, MNRAS, 493, 2059
  * Brammer et al. (2008) Brammer G. B., van Dokkum P. G., Coppi P., 2008, ApJ, 686, 1503
  * Bruzual & Charlot (2003) Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
  * Calzetti et al. (2000) Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storchi-Bergmann T., 2000, ApJ, 533, 682
  * Castellano et al. (2016) Castellano M., et al., 2016, ApJ, 818, L3
  * Castellano et al. (2018) Castellano M., et al., 2018, ApJ, 863, L3
  * Dayal et al. (2009) Dayal P., Ferrara A., Saro A., Salvaterra R., Borgani S., Tornatore L., 2009, MNRAS, 400, 2000
  * Dayal et al. (2011) Dayal P., Maselli A., Ferrara A., 2011, MNRAS, 410, 830
  * De Barros (2017) De Barros S., 2017, Beacons in the dark: using the most distant galaxies to probe cosmic reionization, HST Proposal 15103
  * De Barros et al. (2019) De Barros S., Oesch P. A., Labbé I., Stefanon M., González V., Smit R., Bouwens R. J., Illingworth G. D., 2019, MNRAS, 489, 2355
  * Endsley et al. (2021a) Endsley R., Stark D. P., Chevallard J., Charlot S., 2021a, MNRAS, 500, 5229
  * Endsley et al. (2021b) Endsley R., Stark D. P., Charlot S., Chevallard J., Robertson B., Bouwens R. J., Stefanon M., 2021b, MNRAS, 502, 6044
  * Erb et al. (2014) Erb D. K., et al., 2014, ApJ, 795, 33
  * Finkelstein et al. (2021) Finkelstein S. L., et al., 2021, arXiv e-prints, p. arXiv:2106.13813
  * Fontana et al. (2010) Fontana A., et al., 2010, ApJ, 725, L205
  * Geil et al. (2017) Geil P. M., Mutch S. J., Poole G. B., Duffy A. R., Mesinger A., Wyithe J. S. B., 2017, MNRAS, 472, 1324
  * Greig & Mesinger (2017) Greig B., Mesinger A., 2017, MNRAS, 465, 4838
  * Grogin et al. (2011) Grogin N. A., et al., 2011, ApJS, 197, 35
  * Gronke et al. (2020) Gronke M., et al., 2020, arXiv e-prints, p. arXiv:2004.14496
  * Hashimoto et al. (2018) Hashimoto T., et al., 2018, Nature, 557, 392
  * Hashimoto et al. (2019) Hashimoto T., et al., 2019, PASJ, 71, 71
  * Hoag et al. (2019) Hoag A., et al., 2019, ApJ, 878, 12
  * Hu et al. (2021) Hu W., et al., 2021, Nature Astronomy, 5, 485
  * Hutter et al. (2017) Hutter A., Dayal P., Müller V., Trott C. M., 2017, ApJ, 836, 176
  * Iliev et al. (2006) Iliev I. T., Mellema G., Pen U. L., Merz H., Shapiro P. R., Alvarez M. A., 2006, MNRAS, 369, 1625
  * Illingworth et al. (2016) Illingworth G., et al., 2016, arXiv e-prints, p. arXiv:1606.00841
  * Ishigaki et al. (2016) Ishigaki M., Ouchi M., Harikane Y., 2016, ApJ, 822, 5
  * Jiang et al. (2021) Jiang L., et al., 2021, Nature Astronomy, 5, 256
  * Jung et al. (2020) Jung I., et al., 2020, ApJ, 904, 144
  * Koekemoer et al. (2011) Koekemoer A. M., et al., 2011, ApJS, 197, 36
  * Kulkarni et al. (2019) Kulkarni G., Keating L. C., Haehnelt M. G., Bosman S. E. I., Puchwein E., Chardin J., Aubert D., 2019, MNRAS, 485, L24
  * Labbé et al. (2013) Labbé I., et al., 2013, ApJ, 777, L19
  * Labbé et al. (2015) Labbé I., et al., 2015, ApJS, 221, 23
  * Laporte et al. (2017) Laporte N., Nakajima K., Ellis R. S., Zitrin A., Stark D. P., Mainali R., Roberts-Borsani G. W., 2017, ApJ, 851, 40
  * Laporte et al. (2021) Laporte N., Meyer R. A., Ellis R. S., Robertson B. E., Chisholm J., Roberts-Borsani G. W., 2021, arXiv e-prints, p. arXiv:2104.08168
  * Larson et al. (2020) Larson R. L., et al., 2020, in American Astronomical Society Meeting Abstracts #235. p. 207.23
  * Madau (1995) Madau P., 1995, ApJ, 441, 18
  * Mainali et al. (2018) Mainali R., et al., 2018, MNRAS, 479, 1180
  * Mason et al. (2018) Mason C. A., et al., 2018, ApJ, 857, L11
  * Matthee et al. (2018) Matthee J., Sobral D., Gronke M., Paulino-Afonso A., Stefanon M., Röttgering H., 2018, A&A, 619, A136
  * Matthee et al. (2020) Matthee J., Sobral D., Gronke M., Pezzulli G., Cantalupo S., Röttgering H., Darvish B., Santos S., 2020, MNRAS, 492, 1778
  * McQuinn et al. (2007) McQuinn M., Lidz A., Zahn O., Dutta S., Hernquist L., Zaldarriaga M., 2007, MNRAS, 377, 1043
  * Miralda-Escudé (1998) Miralda-Escudé J., 1998, ApJ, 501, 15
  * Momcheva et al. (2016) Momcheva I. G., et al., 2016, ApJS, 225, 27
  * Mutch et al. (2016) Mutch S. J., Geil P. M., Poole G. B., Angel P. W., Duffy A. R., Mesinger A., Wyithe J. S. B., 2016, MNRAS, 462, 250
  * Naidu et al. (2020) Naidu R. P., Tacchella S., Mason C. A., Bose S., Oesch P. A., Conroy C., 2020, ApJ, 892, 109
  * Naidu et al. (2021) Naidu R., et al., 2021, Where Cosmic Dawn Breaks First: Mapping the Primordial Overdensity Powering a z 9 Ionized Bubble, JWST Proposal. Cycle 1
  * Ocvirk et al. (2020) Ocvirk P., et al., 2020, MNRAS, 496, 4087
  * Oesch et al. (2015) Oesch P. A., et al., 2015, ApJ, 804, L30
  * Oesch et al. (2016) Oesch P. A., et al., 2016, ApJ, 819, 129
  * Ono (2012) Ono Y., 2012, in Umemura M., Omukai K., eds, American Institute of Physics Conference Series Vol. 1480, First Stars IV - from Hayashi to the Future -. pp 277–280, doi:10.1063/1.4754367
  * Ouchi et al. (2020) Ouchi M., Ono Y., Shibuya T., 2020, ARA&A, 58, 617
  * Pentericci et al. (2011) Pentericci L., et al., 2011, ApJ, 743, 132
  * Pentericci et al. (2014) Pentericci L., et al., 2014, ApJ, 793, 113
  * Pentericci et al. (2018) Pentericci L., et al., 2018, A&A, 619, A147
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Poole et al. (2016) Poole G. B., Angel P. W., Mutch S. J., Power C., Duffy A. R., Geil P. M., Mesinger A., Wyithe S. B., 2016, MNRAS, 459, 3025
  * Qin et al. (2017) Qin Y., et al., 2017, MNRAS, 472, 2009
  * Qin et al. (2021) Qin Y., Wyithe J. S. B., Oesch P. A., Illingworth G. D., Leonova E., Mutch S. J., Naidu R. P., 2021, arXiv e-prints, p. arXiv:2108.03675
  * Roberts-Borsani et al. (2016) Roberts-Borsani G. W., et al., 2016, ApJ, 823, 143
  * Robertson (2021) Robertson B. E., 2021, arXiv e-prints, p. arXiv:2110.13160
  * Schenker et al. (2012) Schenker M. A., Stark D. P., Ellis R. S., Robertson B. E., Dunlop J. S., McLure R. J., Kneib J.-P., Richard J., 2012, ApJ, 744, 179
  * Skelton et al. (2014) Skelton R. E., et al., 2014, ApJS, 214, 24
  * Smit et al. (2015) Smit R., et al., 2015, ApJ, 801, 122
  * Smith et al. (2021) Smith A., Kannan R., Garaldi E., Vogelsberger M., Pakmor R., Springel V., Hernquist L., 2021, arXiv e-prints, p. arXiv:2110.02966
  * Stark et al. (2010) Stark D. P., Ellis R. S., Chiu K., Ouchi M., Bunker A., 2010, MNRAS, 408, 1628
  * Stark et al. (2017) Stark D. P., et al., 2017, MNRAS, 464, 469
  * Stefanon et al. (2017) Stefanon M., et al., 2017, ApJS, 229, 32
  * Stefanon et al. (2021a) Stefanon M., Bouwens R. J., Labbé I., Illingworth G. D., Gonzalez V., Oesch P. A., 2021a, arXiv e-prints, p. arXiv:2103.16571
  * Stefanon et al. (2021b) Stefanon M., et al., 2021b, arXiv e-prints, p. arXiv:2110.06226
  * Tilvi et al. (2016) Tilvi V., et al., 2016, ApJ, 827, L14
  * Tilvi et al. (2020) Tilvi V., et al., 2020, ApJ, 891, L10
  * Trapp & Furlanetto (2020) Trapp A. C., Furlanetto S. R., 2020, MNRAS, 499, 2401
  * Trenti & Stiavelli (2008) Trenti M., Stiavelli M., 2008, ApJ, 676, 767
  * Trenti et al. (2012) Trenti M., et al., 2012, ApJ, 746, 55
  * Treu et al. (2013) Treu T., Schmidt K. B., Trenti M., Bradley L. D., Stiavelli M., 2013, ApJ, 775, L29
  * Whitaker et al. (2019) Whitaker K. E., et al., 2019, ApJS, 244, 16
  * Wilkins et al. (2014) Wilkins S. M., Stanway E. R., Bremer M. N., 2014, MNRAS, 439, 1038
  * Willott et al. (2015) Willott C. J., Carilli C. L., Wagg J., Wang R., 2015, ApJ, 807, 180
  * Wyithe & Loeb (2005) Wyithe J. S. B., Loeb A., 2005, ApJ, 625, 1
  * Zitrin et al. (2015) Zitrin A., et al., 2015, ApJ, 810, L12

##  Appendix A Stamps and Photometric Redshift Likelihoods for the High-z
Sample

Here, we show the image stamps, SEDs and photometric redshift likelihood
functions for all sources that were selected in our fields. In 6, we show how
our photometric redshifts compare to the spectroscopic redshifts measured for
the UV-luminous central targets. The neighbouring sources around the four
central targets are show in Fig. 7, 8, 9, and 10. While there is clearly a
relatively large spread, the combined likelihood functions for the neighbours
agree well with the redshifts of the central sources, indicating that the
majority of our selected galaxies indeed have a good chance to be physically
associated with the central UV-luminous sources. Spectroscopy will be needed
to confirm this in the future, however.

Figure 6: Multi-wavelength images and SEDs of the UV-luminous central targets:
three Lyman $\alpha$ emitters and EGS910-3. For each source, the top panels
show, from left to right, HST/ACS $V_{606},\ I_{814}$, HST/WFC3 $Y_{105},\
J_{125},\ H_{160}$, and Spitzer/IRAC $3.6\mu m$ and $4.5\mu m$ stamps. The
lower panels show the SED (purple curve with uncertainties), the photometry
(black squares with uncertainties), and the resulting probability distribution
of redshift $\mathcal{P}(z)$ (green curve). The pink vertical lines represent
the spectroscopic redshift of the Ly$\alpha$ emitters, showing that the
photometric redshifts are accurate for these sources. The last source,
EGS910-3, does not have a spectroscopic redshift measurement yet. Our
photometric redshift is in excellent agreement with the original determination
of Bouwens et al. (2016). Figure 7: Redshift probability distribution
functions $\mathcal{P}(z)$ and stamps of fainter $z_{\mathrm{phot}}>8$ sources
in the pointing around EGSY-z8p7. The $\mathcal{P}(z)$ of the UV luminous
central target is also shown (orange line), indicating that the neighbouring
sources have a significant probability to be physically associated to the
central source. The stamps show the same filters as 6 and are labelled on the
top. Figure 8: The same as Figure 7 for the neighbouring sources around
EGS910-3 with $z_{\mathrm{phot}}>8$. Figure 9: The same as Figure 7 for the
neighbouring sources around EGS-zs8-2 with $z_{\mathrm{phot}}=7-8$. Figure 10:
The same as Figure 7 for the neighbouring sources around EGS-zs8-1 with
$z_{\mathrm{phot}}=7-8$.


$\displaystyle(v_{i}^{2})^{-R/2-\kappa_{i}^{(0)}/2-1}\exp\left\\{-\dfrac{1}{2v_{j}^{2}}\sum_{r=1}^{R}\left[(x_{i}^{(r)}-\mu_{i})-\sum_{X_{k}\in
Pa(X_{i})}\beta_{ki}(x_{k}^{(r)}-\mu_{k})\right]^{2}\right\\}$
$\displaystyle\propto$
$\displaystyle(v_{i}^{2})^{-R/2-\kappa_{i}^{(0)}/2-1}\exp\left\\{-\dfrac{1}{2v_{j}^{2}}\sum_{r=1}^{R}u_{i}^{(r)2}-\dfrac{\lambda_{i}^{(0)}}{2v_{j}^{2}}\right\\}=\mbox{Inv-}\Gamma\left(\dfrac{\kappa_{i}^{(R)}}{2},\dfrac{\lambda_{i}^{(R)}}{2}\right),$

where $\kappa_{i}^{(R)}=\kappa_{i}^{(0)}+R$,
$\lambda_{i}^{(R)}=\lambda_{i}^{(0)}+\sum_{r=1}^{R}u_{i}^{(r)2}$ and
$u_{i}^{(r)}=(x_{i}^{(r)}-\mu_{i})-\sum_{X_{k}\in
Pa(X_{i})}\beta_{ki}(x_{k}^{(r)}-\mu_{k}).$

Third, we derive the conditional posterior distribution of mean coefficient
$\mu_{i}$ with $i=1,2,\ldots,m+1$ for any CPP and CQA,

$\displaystyle
p(\mu_{i}|\mathcal{X},\boldsymbol{\mu}_{-i},\boldsymbol{v}^{2},\boldsymbol{\beta})\propto
p(\mu_{i})\prod_{r=1}^{R}\left[p(x_{i}^{(r)}|x_{Pa(X_{i})}^{(r)})\prod_{j\in\mathcal{S}(X_{i})}p(x_{j}^{(r)}|x_{Pa(X_{j})}^{(r)})\right]$
$\displaystyle\propto\exp\Bigg{\\{}-\dfrac{1}{2v_{i}^{2}}\sum_{r=1}^{R}\left[(x_{i}^{(r)}-\mu_{i})-\sum_{X_{k}\in
Pa(X_{i})}\beta_{ki}(x_{k}^{(r)}-\mu_{k})\right]^{2}$
$\displaystyle~{}~{}~{}~{}-\sum_{r=1}^{R}\sum_{X_{j}\in\mathcal{S}(X_{i})}\dfrac{1}{2v_{j}^{2}}\left[(x_{j}^{(r)}-\mu_{j})-\sum_{X_{k}\in
Pa(X_{j})}\beta_{kj}(x_{k}^{(r)}-\mu_{k})\right]^{2}-\dfrac{1}{2\sigma_{i}^{(0)2}}\left(\mu_{i}-\mu_{i}^{(0)}\right)^{2}\Bigg{\\}},$
$\displaystyle\propto\exp\left\\{-\dfrac{1}{2v_{i}^{2}}\sum_{r=1}^{R}\left(\mu_{i}-a_{i}^{(r)}\right)^{2}-\sum_{r=1}^{R}\sum_{X_{j}\in\mathcal{S}(X_{i})}-\dfrac{1}{2v_{j}^{2}}\left(\beta_{ij}\mu_{i}-c_{ij}^{(r)}\right)^{2}\right.$
$\displaystyle~{}~{}~{}~{}\left.-\dfrac{1}{2\sigma_{i}^{(0)2}}\left(\mu_{i}-\mu_{i}^{(0)}\right)^{2}\right\\},$
$\displaystyle\propto\exp\left\\{-\dfrac{\mu_{i}^{2}}{2}\left(\dfrac{R}{v_{i}^{2}}+\sum_{X_{j}\in
S(X_{i})}\dfrac{R\beta_{ij}^{2}}{v_{j}^{2}}+\dfrac{1}{\sigma_{i}^{(0)2}}\right)+\mu_{i}\left(\sum_{r=1}^{R}\dfrac{a_{i}^{(r)}}{v_{i}^{2}}+\sum_{r=1}^{R}\sum_{X_{j}\in\mathcal{S}(X_{i})}\dfrac{\beta_{ij}c_{ij}^{(r)}}{v_{j}^{2}}\right.\right.$
$\displaystyle~{}~{}~{}~{}\left.\left.+\dfrac{\mu_{i}^{(0)}}{\sigma_{i}^{(0)2}}\right)\right\\}=\mathcal{N}(\mu_{i}^{(R)},\sigma_{i}^{(R)2}),$

where
$\mu_{i}^{(R)}=\sigma_{i}^{(R)2}\left[\dfrac{\mu_{i}^{(0)}}{\sigma_{i}^{(0)2}}+\sum_{r=1}^{R}\dfrac{a_{i}^{(r)}}{v_{i}^{2}}+\sum_{r=1}^{R}\sum_{X_{j}\in
S(X_{i})}\dfrac{\beta_{ij}c_{ij}^{(r)}}{v_{j}^{2}}\right]$ and
$\dfrac{1}{\sigma_{i}^{(R)2}}=\dfrac{1}{\sigma_{i}^{(0)2}}+\dfrac{R}{v_{i}^{2}}+\sum_{X_{j}\in
S(X_{i})}\dfrac{R\beta_{ij}^{2}}{v_{j}^{2}},$ with
$a_{i}^{(r)}=x_{i}^{(r)}-\sum_{X_{k}\in
Pa(X_{i})}\beta_{kj}(x_{k}^{(r)}-\mu_{k})$ and
$c_{ij}^{(r)}=\beta_{ij}x_{i}^{(r)}-(x_{j}^{(r)}-\mu_{j})+\sum_{X_{k}\in
Pa(X_{j})/\\{X_{i}\\}}\beta_{kj}(x_{k}^{(r)}-\mu_{k}).$

###  D.2 Knowledge Learning for Cases with Mixing Data

Except the case with complete production data discussed in Section D.1, we
consider the cases with additional incomplete data corresponding to certain
“Top Sub-Graph", denoted by
$G(\mathbf{N}^{\prime}|\boldsymbol{\theta}(\mathbf{N}^{\prime}))$ with
$\mathbf{N}^{\prime}\subseteq\mathbf{N}$, such that any CQA node
$X_{j}\in\mathbf{N}^{\prime}$ has $Pa(X_{j})\subset\mathbf{N}^{\prime}$. Since
batch data collected from biopharmaceutical production process are usually
limited, we want to fully utilize both complete and incomplete data to
estimate the BN model coefficients and improve our knowledge of production
process.

Without loss of generality, we consider the real-world data including two data
sets $\mathcal{X}=\\{\mathcal{X}_{1},\mathcal{X}_{2}\\}$ with the complete
data
$\mathcal{X}_{1}=\\{(x_{1}^{(r_{1})},x_{2}^{(r_{1})},\ldots,x_{m+1}^{(r_{1})})$
for $r_{1}=1,2,\ldots,R_{1}\\}$ and the incomplete data
$\mathcal{X}_{2}=\\{(x_{i}^{(r_{2})}:X_{i}\in\mathbf{N}^{\prime})$ for
$r_{2}=R_{1}+1,R_{1}+2,\ldots,R\\}$, where $R=R_{1}+R_{2}$. Our approach can
be easily extended to cases with multiple incomplete data sets. We use the
same prior distribution
$p(\boldsymbol{\mu},\boldsymbol{v}^{2},\boldsymbol{\beta})$ as shown in
Equation (13). Given the mixing data
$\mathcal{X}=\\{\mathcal{X}_{1},\mathcal{X}_{2}\\}$, we can derive the
posterior distribution of $\boldsymbol{\theta}$,

$p(\boldsymbol{\mu},\boldsymbol{v}^{2},\boldsymbol{\beta}|\mathcal{X})\propto\prod_{r_{1}=1}^{R_{1}}\left[\prod_{i=1}^{m+1}p(x_{i}^{(r_{1})}|x_{Pa(X_{i})}^{(r_{1})})\right]\prod_{r_{2}=R_{1}+1}^{R}\left[\prod_{X_{i}\in\mathbf{N}^{\prime}}p(x_{i}^{(r_{2})}|x_{Pa(X_{i})}^{(r_{2})})\right]p(\boldsymbol{\mu},\boldsymbol{v}^{2},\boldsymbol{\beta}).$

For $\beta_{ij}$ with $X_{j}\notin\mathbf{N}^{\prime}$ or $v_{i}^{2}$ and
$\mu_{i}$ with node $X_{i}\notin\mathbf{N}^{\prime}$, the conditional
posterior is the same as complete data case and we can utilize Equations (15),
(16) and (17) by replacing $\mathcal{X}$ with $\mathcal{X}_{1}$.

Thus, to derive the full Gibbs sampler, we only need to provide the updated
conditional posterior accounting for those nodes included in the incomplete
data set $\mathcal{X}_{2}$. We first derive the conditional posterior
distribution for weight coefficient $\beta_{ij}$ with
$X_{j}\in\mathbf{N}^{\prime}$.

$\displaystyle
p(\beta_{ij}|\mathcal{X},\boldsymbol{\mu},\boldsymbol{v}^{2},\boldsymbol{\beta}_{-ij})\propto\left[\prod_{r=1}^{R_{1}+R_{2}}p(x_{j}^{(r)}|x_{Pa(X_{j})}^{(r)})\right]p(\beta_{ij}),$
(30) $\displaystyle\propto$
$\displaystyle\exp\left\\{-\sum_{r=1}^{R_{1}+R_{2}}\dfrac{1}{2v_{j}^{2}}\left[(x_{j}^{(r)}-\mu_{j})-\beta_{ij}(x_{i}^{(r)}-\mu_{i})-\sum_{k\in
Pa(j)/\\{i\\}}\beta_{kj}(x_{k}^{(r)}-\mu_{k})\right]^{2}\right.$
$\displaystyle\left.-\dfrac{1}{2\tau_{ij}^{(0)2}}\left(\beta_{ij}-\theta_{ij}^{(0)}\right)^{2}\right\\},$
$\displaystyle\propto$
$\displaystyle\exp\left\\{-\dfrac{1}{2v_{j}^{2}}\sum_{r=1}^{R_{1}+R_{2}}\left(\alpha_{i}^{(r)}\beta_{ij}-m_{ij}^{(r)}\right)^{2}-\dfrac{1}{2\tau_{ij}^{(0)2}}\left(\beta_{ij}-\theta_{ij}^{(0)}\right)^{2}\right\\},$
$\displaystyle\propto$
$\displaystyle\exp\left\\{-\dfrac{\beta_{ij}^{2}}{2}\left(\sum_{r=1}^{R_{1}+R_{2}}\dfrac{\alpha_{i}^{(r)2}}{v_{j}^{2}}+\dfrac{1}{\tau_{ij}^{(0)2}}\right)+\beta_{ij}\left(\sum_{r=1}^{R_{1}+R_{2}}\dfrac{\alpha_{i}^{(r)}m_{ij}^{(r)}}{v_{j}^{2}}+\dfrac{\theta_{ij}^{(0)}}{\tau_{ij}^{(0)2}}\right)\right\\}$
$\displaystyle=$
$\displaystyle\mathcal{N}(\theta_{ij}^{(R_{1}+R_{2})},\tau_{ij}^{(R_{1}+R_{2})2}),$

where
$\theta_{ij}^{(R_{1}+R_{2})}=\dfrac{\tau_{ij}^{(0)2}\sum_{r=1}^{R_{1}+R_{2}}\alpha_{i}^{(r)}m_{ij}^{(r)}+v_{j}^{2}\theta_{ij}^{(0)}}{\tau_{ij}^{(0)2}\sum_{r=1}^{R_{1}+R_{2}}\alpha_{i}^{(r)2}+v_{j}^{2}}$
and
$\tau_{ij}^{(R_{1}+R_{2})2}=\dfrac{\tau_{ij}^{(0)2}v_{j}^{2}}{\tau_{ij}^{(0)2}\sum_{r=1}^{R_{1}+R_{2}}\alpha_{i}^{(r)2}+v_{j}^{2}}$
with $\alpha_{i}^{(r)}=x_{i}^{(r)}-\mu_{i}$ and
$m_{ij}^{(r)}=(x_{j}^{(r)}-\mu_{j})-\sum_{X_{k}\in
Pa(X_{j})/\\{X_{i}\\}}\beta_{kj}(x_{k}^{(r)}-\mu_{k})$ for $r=1,2,\ldots,R$.

Then, we derive the conditional posterior distribution for $v_{i}^{2}$ with
$X_{i}\in\mathbf{N}^{\prime}$,

$\displaystyle
p(v_{i}^{2}|\mathcal{X},\boldsymbol{\mu},\boldsymbol{v}_{-i}^{2},\boldsymbol{\beta})\propto\left[\prod_{r=1}^{R_{1}+R_{2}}p(x_{i}^{(r)}|x_{Pa(X_{i})}^{(r)})\right]p(v_{i}^{2}),$
(31) $\displaystyle\propto$
$\displaystyle(v_{i}^{2})^{-(R_{1}+R_{2})/2-\kappa_{i}^{(0)}/2-1}\exp\left\\{-\dfrac{1}{2v_{j}^{2}}\sum_{r=1}^{R_{1}+R_{2}}\left[(x_{i}^{(r)}-\mu_{i})-\sum_{X_{k}\in
Pa(X_{i})}\beta_{ki}(x_{k}^{(r)}-\mu_{k})\right]^{2}\right\\},$
$\displaystyle\propto$
$\displaystyle(v_{i}^{2})^{-(R_{1}+R_{2})/2-\kappa_{i}^{(0)}/2-1}\exp\left\\{-\dfrac{1}{2v_{j}^{2}}\sum_{r=1}^{R_{1}+R_{2}}u_{i}^{(r)2}-\dfrac{\lambda_{i}^{(0)}}{2v_{j}^{2}}\right\\}$
$\displaystyle=$
$\displaystyle\mbox{Inv-}\Gamma\left(\dfrac{\kappa_{i}^{(R_{1}+R_{2})}}{2},\dfrac{\lambda_{i}^{(R_{1}+R_{2})}}{2}\right),$

where $\kappa_{i}^{(R_{1}+R_{2})}=\kappa_{i}^{(0)}+R$ and
$\lambda_{i}^{(R_{1}+R_{2})}=\lambda_{i}^{(0)}+\sum_{r=1}^{R}u_{i}^{(r)2}$
with $u_{i}^{(r)}=(x_{i}^{(r)}-\mu_{i})-\sum_{X_{k}\in
Pa(X_{i})}\beta_{ki}(x_{k}^{(r)}-\mu_{k})$ for $r=1,2,\ldots,R$.

After that, we derive the conditional posterior for mean coefficient $\mu_{i}$
with $X_{i}\in\mathbf{N}^{\prime}$,

$\displaystyle
p(\mu_{i}|\mathcal{X},\boldsymbol{\mu}_{-i},\boldsymbol{v}^{2},\boldsymbol{\beta})\propto
p(\mu_{i})\prod_{r=1}^{R_{1}+R_{2}}p(x_{i}^{(r)}|x_{Pa(X_{i})}^{(r)})\prod_{r_{1}=1}^{R_{1}}\prod_{X_{j}\in\mathcal{S}(X_{i})}p(x_{j}^{(r_{1})}|x_{Pa(X_{j})}^{(r_{1})})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\cdot\prod_{r_{2}=R_{1}+1}^{R_{1}+R_{2}}\prod_{X_{j}\in
S(X_{i})\cap\mathbf{N}^{\prime}}p(x_{j}^{(r_{2})}|x_{Pa(X_{j})}^{(r_{2})}),$
$\displaystyle\propto\exp\Bigg{\\{}-\dfrac{1}{2v_{i}^{2}}\sum_{r=1}^{R_{1}+R_{2}}\bigg{[}(x_{i}^{(r)}-\mu_{i})-\sum_{X_{k}\in
Pa(X_{i})}\beta_{ki}(x_{k}^{(r)}-\mu_{k})\bigg{]}^{2}$
$\displaystyle~{}~{}~{}~{}-\sum_{r_{1}=1}^{R_{1}}\sum_{X_{j}\in\mathcal{S}(X_{i})}\dfrac{1}{2v_{j}^{2}}\bigg{[}(x_{j}^{(r_{1})}-\mu_{j})-\sum_{X_{k}\in
Pa(X_{j})}\beta_{kj}(x_{k}^{(r_{1})}-\mu_{k})\bigg{]}^{2}$
$\displaystyle~{}~{}~{}~{}-\sum_{r_{2}=R_{1}+1}^{R_{1}+R_{2}}\sum_{X_{j}\in
S(X_{i})\cap\mathbf{N}^{\prime}}\dfrac{1}{2v_{j}^{2}}\bigg{[}(x_{j}^{(r_{2})}-\mu_{j})-\sum_{X_{k}\in
Pa(X_{j})}\beta_{kj}(x_{k}^{(r_{2})}-\mu_{k})\bigg{]}^{2}$
$\displaystyle~{}~{}~{}~{}-\dfrac{1}{2\sigma_{i}^{(0)2}}\left(\mu_{i}-\mu_{i}^{(0)}\right)^{2}\Bigg{\\}},$
$\displaystyle\propto\exp\left\\{-\dfrac{\mu_{i}^{2}}{2}\left(\dfrac{R_{1}+R_{2}}{v_{i}^{2}}+\sum_{X_{j}\in\mathcal{S}(X_{i})}\dfrac{R_{1}\beta_{ij}^{2}}{v_{j}^{2}}+\sum_{X_{j}\in
S(X_{i})\cap\mathbf{N}^{\prime}}\dfrac{R_{2}\beta_{ij}^{2}}{v_{j}^{2}}\dfrac{1}{\sigma_{i}^{(0)2}}\right)\right.$
$\displaystyle~{}~{}~{}~{}+\left.\mu_{i}\left(\sum_{r=1}^{R_{1}+R_{2}}\dfrac{a_{i}^{(r)}}{v_{i}^{2}}+\sum_{r_{1}=1}^{R_{1}}\sum_{X_{j}\in\mathcal{S}(X_{i})}\dfrac{\beta_{ij}c_{ij}^{(r_{1})}}{v_{j}^{2}}+\sum_{r_{2}=R_{1}+1}^{R}\sum_{X_{j}\in\mathcal{S}(X_{i})\cap\mathbf{N}^{\prime}}\dfrac{\beta_{ij}c_{ij}^{(r_{2})}}{v_{j}^{2}}+\dfrac{\mu_{i}^{(0)}}{\sigma_{i}^{(0)2}}\right)\right\\},$
$\displaystyle=\mathcal{N}\left(\mu_{i}^{(R_{1}+R_{2})},\sigma_{i}^{(R_{1}+R_{2})2}\right),$
(32)

$\mu_{i}^{(R_{1}+R_{2})}=\sigma_{i}^{(R_{1}+R_{2})2}\left[\dfrac{\mu_{i}^{(0)}}{\sigma_{i}^{(0)2}}+\sum_{r=1}^{R_{1}+R_{2}}\dfrac{a_{i}^{(r)}}{v_{i}^{2}}+\sum_{r_{1}=1}^{R_{1}}\sum_{X_{j}\in\mathcal{S}(X_{i})}\dfrac{\beta_{ij}c_{ij}^{(r_{1})}}{v_{j}^{2}}+\right.\\\
\left.\sum_{r_{2}=R_{1}+1}^{R}\sum_{X_{j}\in\mathcal{S}(X_{i})\cap\mathbf{N}^{\prime}}\dfrac{\beta_{ij}c_{ij}^{(r_{2})}}{v_{j}^{2}}\right]$,
and
$\dfrac{1}{\sigma_{i}^{(R_{1}+R_{2})2}}=\dfrac{1}{\sigma_{i}^{(0)2}}+\dfrac{R_{1}+R_{2}}{v_{i}^{2}}+\sum_{X_{j}\in\mathcal{S}(X_{i})}\dfrac{R_{1}\beta_{ij}^{2}}{v_{j}^{2}}+\sum_{X_{j}\in\mathcal{S}(X_{i})\cap\mathbf{N}^{\prime}}\dfrac{R_{2}\beta_{ij}^{2}}{v_{j}^{2}}$
with $a_{i}^{(r)}=x_{i}^{(r)}-\sum_{X_{k}\in
Pa(X_{i})}\beta_{kj}(x_{k}^{(r)}-\mu_{k})$ and
$c_{ij}^{(r)}=\beta_{ij}x_{i}^{(r)}-(x_{j}^{(r)}-\mu_{j})+\sum_{X_{k}\in
Pa(X_{j})/\\{X_{i}\\}}\beta_{kj}(x_{k}^{(r)}-\mu_{k})$ for $r=1,2,\ldots,R$.
Here for illustration, we have only provided the conditional posteriors with
two datasets $\mathcal{X}_{1}$ and $\mathcal{X}_{2}$. These derivations can be
easily extended to similar cases with multiple datasets collected from
complete graph and different top sub-graphs.

###  D.3 Gibbs Sampling Procedure for BN Model Bayesian Inference

Based on the derived conditional posterior distributions in Sections D.1 and
D.2, we provide the Gibbs sampling procedure in Algorithm 3 to generate
posterior samples $\widetilde{\boldsymbol{\theta}}^{(b)}\sim
p(\boldsymbol{\theta}|\mathcal{X})$ with
$\widetilde{\boldsymbol{\theta}}^{(b)}=(\widetilde{\boldsymbol{\mu}}^{(b)},\widetilde{\boldsymbol{v}}^{(b)2},\widetilde{\boldsymbol{\beta}}^{(b)})$
and $b=1,\ldots,B$. We first set the vague prior
$p(\boldsymbol{\theta})=p(\boldsymbol{\mu},\boldsymbol{v}^{2},\boldsymbol{\beta})$
as Equation (13), and generate the initial point
$\boldsymbol{\theta}^{(0)}=(\boldsymbol{\mu}^{(0)},\boldsymbol{v}^{(0)2},\boldsymbol{\beta}^{(0)})$
by sampling from the prior. Within each $t$-th iteration of Gibbs sampling,
given the previous sample
$\boldsymbol{\theta}^{(t-1)}=(\boldsymbol{\mu}^{(t-1)},\boldsymbol{v}^{(t-1)2},\boldsymbol{\beta}^{(t-1)})$,
we sequentially compute and generate one sample from the conditional posterior
distribution for each coefficient $\beta_{ij}$, $v_{i}^{2}$ and $\mu_{i}$. By
repeating this procedure, we can get samples
$\boldsymbol{\theta}^{(t)}=(\boldsymbol{\mu}^{(t)},\boldsymbol{v}^{(t)2},\boldsymbol{\beta}^{(t)})$
with $t=1,\ldots,T$. To reduce the initial bias and correlations between
consecutive samples, we remove the first $T_{0}$ samples and keep one for
every $h$ samples. Consequently, we obtain the posterior samples
$\widetilde{\boldsymbol{\theta}}^{(b)}\sim p(\boldsymbol{\theta}|\mathcal{X})$
with $b=1,\ldots,B$.

Input: the prior $p(\boldsymbol{\theta})$ and real-world data $\mathcal{X}$.

Output: Posterior samples
$\widetilde{\boldsymbol{\theta}}^{(b)}=(\widetilde{\boldsymbol{\mu}}^{(b)},\widetilde{\boldsymbol{v}}^{(b)2},\widetilde{\boldsymbol{\beta}}^{(b)})\sim
p(\boldsymbol{\theta}|\mathcal{X})$ with $b=1,\ldots,B$.

(1) Set the initial value
$\boldsymbol{\theta}^{(0)}=(\boldsymbol{\mu}^{(0)},\boldsymbol{v}^{(0)2},\boldsymbol{\beta}^{(0)})$
by sampling from prior $p(\boldsymbol{\theta})$;

for  _$t=1,2,\ldots,T$_ do

(2) Given the previous sample
$\boldsymbol{\theta}^{(t-1)}=(\boldsymbol{\mu}^{(t-1)},\boldsymbol{v}^{(t-1)2},\boldsymbol{\beta}^{(t-1)})$;

(3) For each $\beta_{ij}$, generate $\beta_{ij}^{(t)}\sim
p(\beta_{ij}|\mathcal{X},\beta_{12}^{(t)},\ldots,$
$\beta_{i,j-1}^{(t)},\beta_{i,j+1}^{(t-1)},\ldots,\beta_{m,m+1}^{(t-1)},\boldsymbol{\mu}^{(t-1)},\boldsymbol{v}^{(t-1)2})$
through Equation (15) for complete data or (30) for mixing data;

(4) For each $v_{i}^{2}$, generate $v_{i}^{(t)2}\sim
p(v_{i}^{2}|\mathcal{X},\boldsymbol{\beta}^{(t)},v_{1}^{(t)2},$
$\ldots,v_{i-1}^{(t)2},v_{i+1}^{(t-1)2},\ldots,v_{m+1}^{(t-1)2},\boldsymbol{\mu}^{(t-1)}))$
through Equation (16) for complete data or (31) for mixing data;

(5) For each $\mu_{i}$, generate $\mu_{i}^{(t)}\sim
p(\mu_{i}|\mathcal{X},\boldsymbol{\beta}^{(t)},\boldsymbol{v}^{(t),2},$
$\mu_{1}^{(t)},\ldots,\mu_{i-1}^{(t)},\mu_{i+1}^{(t-1)},\ldots,\mu_{n}^{(t-1)})$
through Equation (17) for complete data or (32) for mixing data;

(6) Obtain a new posterior sample
$\boldsymbol{\theta}^{(t)}=(\boldsymbol{\mu}^{(t)},\boldsymbol{v}^{(t)2},\boldsymbol{\beta}^{(t)})$;

(7) Set
$\widetilde{\boldsymbol{\theta}}^{(b)}=\boldsymbol{\theta}^{(T_{0}+(b-1)h+1)}$
with some constant integer $T_{0}$ and $h$, to reduce the initial bias and
correlation between consecutive samples.

Algorithm 3 Gibbs Sampling Procedure for BN Model Uncertainty Quantification

##  Appendix E Simulated Biopharmaceutical Production Data

To study the performance of proposed framework, we generate the simulated
production process data $\mathcal{X}$, which mimics the “real-world data
collection." The BN with coefficients $\boldsymbol{\theta}^{c}$ characterizing
the underlying production process interdependence is used for data generation,
which is built according to the biomanufacturing domain knowledge. The ranges
of CPPs/CQAs are listed Table 9. For each CPP $X_{j}\in\mathbf{X}^{p}$ with
range $[x_{j}^{low},x_{j}^{up}]$, we can specify the marginal distribution
$X_{j}\sim\mathcal{N}(\mu_{j}^{c},(v_{j}^{c})^{2})$ with mean
$\mu^{c}_{j}=(x_{j}^{low}+x_{j}^{up})/2$ and standard deviation
$v_{j}^{c}=(x_{j}^{up}-x_{j}^{low})/4$. For each CQA
$X_{i}\in\\{\mathbf{X}^{a}\cup\mathbf{Y}\\}$ with range
$[x_{i}^{low},x_{i}^{up}]$, we have mean
$\mu_{i}^{c}=(x_{i}^{low}+x_{i}^{up})/2$ and marginal variance
$\mbox{Var}(X_{i})=[(x_{i}^{up}-x_{i}^{low})/4]^{2}$. Based on Equation (12),
the corresponding coefficient $v^{c}_{i}$ can be computed through back-
engineering. For the complex interdependence, Table 10 provides the relative
associations with levels (i.e., high, median, low) between input CPPs/CQAs
with output CQAs in each operation unit, which is built based on the “cause-
and-effect matrix" in [31]. For the high, median and low association between
$X_{i}$ to $X_{j}$, we set the coefficient $\beta^{c}_{ij}=0.9,0.6,0.3$
respectively. Thus, we can specify the underlying true coefficients
$\boldsymbol{\theta}^{c}=(\boldsymbol{\mu}^{c},(\boldsymbol{v}^{2})^{c},\boldsymbol{\beta}^{c})$.
To mimic the “real-world" data collection, we generate the production batch
data $\mathcal{X}$ using the BN model with $\boldsymbol{\theta}^{c}$. Then, to
assess the performance of proposed framework, we assume that the true
coefficient values are unknown.

Table 9: Range of CPPs/CQAs in the production procedure.

Process Unit Operation | CPP | Range | CQA | Range
---|---|---|---|---
Main Fermentation | pH | 6.8-7.2 | impurities | 3-11 pl
| temperature | 20-30 C | protein content | 1-5 g/L
| Oxygen | 2.5-7.5% | bioburden | 5-15 CFU/100mL
| agitation rate | 1.1-2.5 m/s |  |
Centrifuge | temperature | 20 to 30 C | impurities | 3-11 pl
| rotation speed | 3-5K RPM | protein content | 5-15 CFU/100mL
Chromatography | pooling window | 10-30 min | impurities | 3-11 pl
| temperature | 2-10 C | protein content | 1-5 g/L
|  |  | bioburden | 5-15 CFU/100mL
Filtration | size of sieve | 0.1-0.5 um | impurities | 3-11 pl
| flow rate | 25-100 mL/min | protein content | 1-5 g/L

Table 10: Relative association between input CPPs/CQAs with output CQAs in
each process unit operation.

Process Unit Operation | Input CPPs/CQAs | Output CQAs
---|---|---
|  | impurities | protein content | bioburden
Main Fermentation | pH | high | high | low
| temperature | high | high | low
| Oxygen | high | high | low
| agitation rate | high | high | low
Centrifuge | temperature | medium | medium | —
| rotation speed | medium | medium | —
| impurities (main fermentation) | medium | medium | —
| protein content (main fermentation) | medium | medium | —
| bioburden (main fermentation) | medium | medium | —
Chromatography | pooling window | high | medium | high
| temperature | high | medium | high
| impurities (centrifuge) | high | medium | high
| protein content (centrifuge) | high | medium | high
Filtration | size of sieve | low | medium | medium
| flow rate | low | medium | medium
| impurities (chromatography) | low | medium | medium
| protein content (chromatography) | low | medium | medium
| bioburden (chromatography) | low | medium | medium

##  Appendix F Study the Bayesian Learning and Inference

To evaluate the accuracy and efficiency of proposed Bayesian learning, we
empirically study the convergence of BN coefficient inference. In each $k$-th
macro-replication, we first mimic the “real-world" production batch data
collection through generating
$\mathcal{X}^{(k)}=\\{\mathbf{X}^{(k)}_{1},\ldots,\mathbf{X}^{(k)}_{R}\\}$
with $\mathbf{X}_{i}^{(k)}\sim F(\mathbf{X}|\boldsymbol{\theta}^{c})$ for
$i=1,\ldots,R$ and $k=1,\ldots,K$. Then, we generate $B$ posterior samples
$\widetilde{\boldsymbol{\theta}}^{(k,b)}\sim
p(\boldsymbol{\theta}|\mathcal{X}^{(k)})$ with $b=1,2,\ldots,B$. For the Gibbs
sampler in Algorithm 3 provided in online Appendix D.3, we set the initial
warm-up length $T_{0}=500$ and step-size $h=10$. With different size of
complete “real-world" batch data $R=30,100,500$, we compute the mean squared
error (MSE) for each coefficient $\theta_{\ell}\in\boldsymbol{\theta}$:
$\mbox{MSE}(\theta_{\ell})=\iint\left({\theta_{\ell}}-\theta_{\ell}^{c}\right)^{2}dP(\theta_{\ell}|\mathcal{X})dP(\mathcal{X}|\boldsymbol{\theta}^{c}).$
Based on $K=20$ macro-replications and $B=1000$ posterior samples of BN
coefficients, we estimate $\mbox{MSE}(\theta_{\ell})$ with
$\widehat{\mbox{MSE}}(\theta_{\ell})=\dfrac{1}{KB}\sum_{k=1}^{K}\sum_{b=1}^{B}\left(\widetilde{\theta}_{\ell}^{(k,b)}-\theta_{\ell}^{c}\right)^{2}.$
Since the total number of coefficients is large, we further group coefficients
by mean $\boldsymbol{\mu}$, conditional variance $\boldsymbol{v}^{2}$ and
linear coefficients $\boldsymbol{\beta}$, and take average of the sample MSE
respectively:
$\widehat{\mbox{MSE}}(\boldsymbol{\mu})=\frac{1}{|\boldsymbol{\mu}|}\sum_{\theta_{\ell}\in\boldsymbol{\mu}}\widehat{\mbox{MSE}}(\theta_{\ell})$,
$\widehat{\mbox{MSE}}(\boldsymbol{v}^{2})=\frac{1}{|\boldsymbol{v}^{2}|}\sum_{\theta_{\ell}\in\boldsymbol{v}^{2}}\widehat{\mbox{MSE}}(\theta_{\ell})$,
and
$\widehat{\mbox{MSE}}(\boldsymbol{\beta})=\frac{1}{|\boldsymbol{\beta}|}\sum_{\theta_{\ell}\in\boldsymbol{\beta}}\widehat{\mbox{MSE}}(\theta_{\ell})$.
The corresponding results are reported in Table 11. As the size of real-world
data $R$ increases, the average MSE decreases, which implies the posterior
samples obtained by Gibbs sampling procedure can converge to the true
coefficients.

Table 11: The MSE of $\boldsymbol{\mu}$, $\boldsymbol{v}^{2}$ and $\boldsymbol{\beta}$ esimated by using the Gibbs sampling. Batch Data Size | $\widehat{\mbox{MSE}}(\boldsymbol{\mu})$ | $\widehat{\mbox{MSE}}(\boldsymbol{v}^{2})$ | $\widehat{\mbox{MSE}}(\boldsymbol{\beta})$
---|---|---|---
$R=30$ |  0.122$\pm$0.032 |  0.276$\pm$0.029 |  0.0225$\pm$0.0013
$R=100$ |  0.075$\pm$0.023 |  0.061$\pm$0.006 |  0.0063$\pm$0.0004
$R=500$ |  0.009$\pm$0.003 |  0.013$\pm$0.001 |  0.0011$\pm$0.00004