UV&Ly$\alpha$ halos of Ly$\alpha$ emitters across environments at z=2.84

We present UV and Ly$\alpha$ radial surface brightness (SB) profiles of Ly$\alpha$ emitters (LAEs) at $z=2.84$ detected with the Hyper Suprime-Cam (HSC) on the Subaru Telescope. The depth of our data, together with the wide field coverage including a protocluster, enable us to study the dependence of Ly$\alpha$ halos (LAHs) on various galaxy properties, including Mpc-scale environments. UV and Ly$\alpha$ images of 3490 LAEs are extracted, and stacking the images yields SB sensitivity of $\sim1\times10^{-20}\mathrm{~erg~s^{-1}~cm^{-2}~arcsec^{-2}}$ in Ly$\alpha$, reaching the expected level of optically thick gas illuminated by the UV background at $z\sim3$. Fitting of the two-component exponential function gives the scale-lengths of $1.56\pm0.01$ and $10.4\pm0.3$ pkpc. Dividing the sample according to their photometric properties, we find that while the dependence of halo scale-length on environment outside of the protocluster core is not clear, LAEs in the central regions of protoclusters appear to have very large LAHs which could be caused by combined effects of source overlapping and diffuse Ly$\alpha$ emission from cool intergalactic gas permeating the forming protocluster core irradiated by active members. For the first time, we identify ``UV halos'' around bright LAEs which are probably due to a few lower-mass satellite galaxies. Through comparison with recent numerical simulations, we conclude that, while scattered Ly$\alpha$ photons from the host galaxies are dominant, star formation in satellites evidently contributes to LAHs, and that fluorescent Ly$\alpha$ emission may be boosted within protocluster cores at cosmic noon and/or near bright QSOs.


INTRODUCTION
Corresponding author: Satoshi Kikuta satoshi.kikuta@nao.ac.jp Gas inflow from the cosmic web fuels star formation and active galactic nucleus (AGN) activity in galaxies. Following these activities, outflows play a significant role in expelling or heating a large amount of gas, thereby reducing further star formation and supermassive black hole (SMBH) growth. The circumgalactic medium (CGM) contains vital information on these flow components and thus holds an important key to revealing galaxy evolution (see Tumlinson et al. 2017, for a recent review). The CGM of z 2 star-forming galaxies is now routinely detected as diffuse Lyα nebulae, or Lyα halos (LAHs), around star-forming galaxies such as Lyα emitters (LAEs) and Lyman break galaxies (LBGs) at high-redshift both individually (Rauch et al. 2008;Wisotzki et al. 2016;Leclercq et al. 2017;Erb et al. 2018;Bacon et al. 2021;Kusakabe et al. 2022) and through a stacking technique (Hayashino et al. 2004;Steidel et al. 2011;Matsuda et al. 2012;Feldmeier et al. 2013;Momose et al. 2014Momose et al. , 2016Xue et al. 2017;Lujan Niemeyer et al. 2022a,b). Together with a technique based on absorption lines in the spectra of neighboring background sources (Adelberger et al. 2005;Steidel et al. 2010;Rudie et al. 2012;Chen et al. 2020;Muzahid et al. 2021), LAHs have provided a crucial empirical window into the CGM of distant galaxies.
To extract useful information on the CGM from LAHs, the physical origins of the Lyα emission should be identified. Lyα surface brightness (SB) profiles of LAHs hold the key since they are determined by the distribution and kinematics of gas and the relative importance of various Lyα production mechanisms such as scattering of Lyα photons from host galaxies, star formation in neighboring galaxies, collisional excitation of inflow gas powered by gravitational energy (sometimes called gravitational cooling radiation), and recombination following photoionization by external sources, often referred to as "fluorescence." Theoretical studies have attempted to reproduce and predict observed LAHs by considering these mechanisms (Haiman & Rees 2001;Dijkstra & Loeb 2009;Goerdt et al. 2010;Kollmeier et al. 2010;Faucher-Giguère et al. 2010;Zheng et al. 2011;Dijkstra & Kramer 2012;Rosdahl & Blaizot 2012;Yajima et al. 2013;Cen & Zheng 2013;Cantalupo et al. 2014;Lake et al. 2015;Mas-Ribas & Dijkstra 2016;Gronke & Bird 2017;Mitchell et al. 2021;Byrohl et al. 2021). Powerful outflows (so-called "superwind", Taniguchi & Shioya 2000;Mori et al. 2004) have also been proposed to excite gas, but often for more energetic/massive counterparts such as Lyα blobs (LABs) and nebulae around QSOs and radio galaxies 1 .
Dependence of LAH shapes on e.g., their hosts' halo mass and large-scale overdensity is naturally expected because both gas and sources of Lyα and ionizing photons are more abundant in massive halos and/or denser environments (Zheng et al. 2011;Mas-Ribas & Dijkstra 2016;Kakiichi & Dijkstra 2018). Current simulations cannot treat all relevant physics with sufficient accuracy and it is only very recently that such predictions are reported in the literature with a statistical number of simulated galaxies (e.g., Byrohl et al. 2021). Observations of LAHs can help theorists pin down which Lyα production processes are at work by revealing the dependence of LAH SB profile shapes on various properties such as UV and Lyα luminosity, Lyα equivalent width (EW Lyα,0 , Momose et al. 2014Momose et al. , 2016Wisotzki et al. 2016Wisotzki et al. , 2018Leclercq et al. 2017), and the large-scale environment Xue et al. 2017). The results in the literature are, however, far from converging (see e.g., Figure 12 of Leclercq et al. (2017)). It is often parametrized by an exponential function ∝ exp(−r/r h ) with a scale-length r h , which is fit to observed profiles. The reported scale-lengths of individual LAEs as a function of UV magnitude show a large scatter (from < 1 physical kpc (pkpc hereafter) to ∼ 10 pkpc) and the relation for stacked LAEs shows large differences as well.
In the case of large-scale environment, relevant observations of LAHs of LAEs are still scarce. First, Steidel et al. (2011) found very large LAHs with a scalelength of 25 pkpc around LBGs in three protoclusters at z = 2.3-3.1. Following this result, Matsuda et al. (2012) suggested that the scale-length of LAHs of LAEs are proportional to galaxy overdensity squared δ 2 . On the other hand, Xue et al. (2017) found no such dependence with LAEs in two overdense regions at z = 2.66 and z = 3.78.
A major problem with some previous work is poor sensitivity.
While only a few studies investigated LAHs with deep images of a fairly large sample of LAEs (N LAE > 2000; Matsuda et al. 2012;Momose et al. 2014), others used the insufficient number of LAEs (N LAE ∼ a few×100-1000) and/or images taken with 4m telescopes (e.g., Feldmeier et al. 2013;Xue et al. 2017).
Because LAHs beyond the virial radii of LAEs have extremely low SB (< 10 −19 erg s −1 cm −2 arcsec −2 ), interpretation of relations derived from shallow data would not be straightforward. Even with sufficiently deep data, the extent of LAHs is difficult to measure. Disagreements among different studies can be attributed in part to different fitting methods, fitting range, radial bin size of SB profiles, and sensitivity of observational data among each study. To alter this situation and to provide a firm observational basis for theorists, a well-controlled statistical sample of LAEs drawn from a wide dynamic range of environments, with sufficiently deep images, is required.
In this paper, we present a new LAH study with deep narrow-band (NB468) data obtained using the Hyper Suprime-Cam (HSC; Miyazaki et al. 2012) on the Subaru Telescope toward the HS1549 protocluster at z = 2.84 Mostardi et al. 2013;Kikuta et al. 2019) to probe what shapes LAHs. Thanks to the HSC's large field of view (φ ∼ 1.5 deg, corresponding to 160 comoving Mpc at z = 2.84), we can construct a large LAE sample across environments from a protocluster to surrounding lower density fields at the same time. The sample size of our study of N = 3490 is one of the largest to date, giving robust UV and Lyα SB profiles to be compared with simulations. As a result, we for the first time detect "UV halos" which directly prove the contribution of star formation in satellite galaxies. Moreover, we detect very extended LAHs for the protocluster LAEs which suggest an important role of locally enhanced ionizing radiation fields for LAHs. This paper is structured as follows. In Section 2, we describe our LAE sample, followed by how we divide them for the stacking analyses described in Section 3. The results of the analyses are shown in Section 4. Based on these, we present discussion in Section 5 and summarize the work in Section 6. Throughout this paper, we use the AB magnitude system and assume a cosmology with Ω m = 0.3, Ω Λ = 0.7, and H 0 = 70 km s −1 Mpc −1 , unless otherwise noted.

DATA AND SAMPLE
We used the LAE sample described by Kikuta et al. (2019). Here, we highlight key properties of the data and refer readers to the aforementioned paper for details. The target protocluster contains a hyperluminous QSO, namely HS1549+1919 (L 1450Å = 1.5 × 10 14 L , Mostardi et al. 2013) at its center. The field was observed in the g-band (central wavelength λ c = 4712Å, FWHM = 1479Å) and NB468 (λ c = 4683Å, FWHM = 88Å) narrow-band filters. The global sky subtraction method 2 was used to estimate and subtract the sky on scales larger than that of individual CCDs in the mosaic, with a grid size of 6000 pixels (17 ) not to subtract diffuse emission. The FWHMs of stellar sources in the final images are 0. 77 (0. 65) for g-band (NB468). The NB468 image was smoothed with a circular Gaussian function to match the FWHM of stellar sources in the g-band image. A 5σ limiting magnitude measured with 1. 5 diameter aperture is 27.4 (26.6) mag for the g-band (NB468) image. Our criterion for NB excess, g − NB468 > 0.5, corresponds to EW obs > 45Å after considering the 0.1 mag offset. SExtractor (Bertin & Arnouts 1996) was used to perform 1. 5 aperture photometry with double-image mode, using the NB468 image as the detection band and a background mesh size of 64 pixels (= 11 ) for local sky estimation. This small mesh size is only used for LAE detection, since this is optimal for detecting compact sources such as distant galaxies. As a result, we detected 3490 LAEs within r < 36 from the QSO position. Their sky distribution is shown in Figure 1.
To study the dependence of LAHs on various photometric galaxy properties, we divided the sample into several groups according to the following five quantities: UV magnitude, Lyα luminosity, rest-frame Lyα equivalent width (EW 0 ), environment, and distance from the HLQSO, as summarized in Table 1. The first three quantities, which are derived with 1. 5 aperture (= 9 pixels = 2 × [PSF FWHM]), are obviously not independent of one another. UV continuum and Lyα luminosity are both good proxies for star formation rate (SFR), but the latter's resonant nature and a different impact of dust extinction thereby make interpretation difficult (e.g., Scarlata et al. 2009). Moreover, UV slope or hardness of UV emission is a strong function of age and metallicity, and thus the Lyα equivalent width changes accordingly (Schaerer 2003;Hashimoto et al. 2017a). There is a known observational relation between the EW distribution and UV limiting magnitude, known as "the Ando effect" (a fainter UV threshold tends to include more high-EW LAEs; Ando et al. 2006) 3 . Since Lyα emission can be powered by mechanisms other than star formation in galaxies, UV magnitude is the most robust to use here as a tracer for SFR. Binning in UV magnitude, Lyα luminosity, and Lyα equivalent width (as well as distance from the HLQSO) were done so that each subsample has approximately the same number of LAEs (N LAE ∼ 700).
The projected distance from the HLQSO is used to test whether QSO radiation affects the LAHs of surrounding LAEs. HS1549 is so luminous that the entire field covered by the HSC could experience a higher ionizing radiation field than the cosmic average at z ∼ 3 (Haardt & Madau 2012) if QSO radiation has had time to propagate, and additional ionization induced by the QSO could increase Lyα luminosity of LAEs in the field (see Section 5.2.1 for discussion). Here, we use a projected distance from the QSO. Note, however, that the NB468 filter has an FWHM of ∆λ = 88Å or ∆z = 0.075, corresponding to 19 pMpc width when centered at z = 2.84. This brings uncertainty in a line-ofsight distance and therefore also in real (3D) distance.
3 Note that this effect could be purely attributed to the selection bias and not intrinsic (Nilsson et al. 2009a;Zheng et al. 2010).
The boundaries defining the LAE subsamples are indicated by concentric circles in Figure 1. Lastly, grouping based on environment is done using projected (surface) LAE overdensity δ ≡ (n −n)/n measured locally with an aperture radius of 1. 8 (= 0.83 pMpc) to be consistent with the measurement of Matsuda et al. (2012). Here, n andn are, respectively, the number of LAEs within a circle with a 1. 8 radius centered at the position of interest, and its average over the entire field. This division is visually illustrated in Figure 1 by gray contours (see also Figures 1 and 2 of Kikuta et al. 2019). The boundary of δ = 2.5 is set to only include protocluster members in the densest subgroup by visual inspection. The next boundary δ = 1 is set by a trade-off between tracing sufficiently dense regions and including adequate numbers to allow sufficient S/N in stacks. The remainder are set so as to roughly equalize the number of LAEs in each bin. Figure 2 shows cumulative distributions of UV and Lyα luminosity, rest-frame Lyα equivalent width, distance from the HLQSO, and overdensity δ for each division, illustrating correlations between these quantities. We note that the odd behavior of the thin blue and purple curves in the second panel from left in the third row of Figure 2 is likely to be artificial; EWs of LAEs not detected in g band (above 2σ) are just lower limits. The samples of faint low-EW LAEs are also incomplete due to the lack of dynamic range in the measurement of g − NB468, distorting the distribution. The protocluster subsample (LAEs with δ > 2.5; thick red curves in panels shown with δ) evidently stands out among others, while the projected distance from the HLQSO seems not to make a significant difference except for δ. These differences should be kept in mind when interpreting the result in Section 4. We visualize these quantities also in Figure 1 which roughly includes all of above information.  (1)  3. ANALYSES

Image Stacking
Before stacking, additional sky subtraction needs to be made, as the global sky subtraction (see Section 2) alone is generally insufficient to eliminate the influence of artifacts such as halos of bright stars. Since we are interested in diffuse and extended components, we use a background mesh size of 176 pixels (= 30 arcsec) for additional background evaluation of g-band and NB images using SExtractor and then subtract the sky. Then a Lyα (continuum) image was created by subtracting the g-band (NB468) image from the NB468 (g-band) image after scaling by their relative zero points and considering the difference in the filters' transmission curves assuming flat continuum (see Appendix B of Mawatari et al. 2012, for details). A segmentation image of the continuum image was used for masking. The segmentation image is an output of SExtractor and it specifies which pixel is detected as a source. If detected, that pixel has a non-zero integer value that corresponds to the ID number in the output catalog. We set DETECT MINAREA and DETECT THRESH parameters to 5 and 2.0. We confirmed that the overall results do not depend significantly on the choice of the threshold value.
We created cutout Lyα and UV continuum images centered on each LAE. The centers of the LAEs are identified as their centroids in the NB image. The mask is applied to each cutout image. If the LAE at the center of the cutout image is detected in the continuum image, masking for the object is turned off so as not to underestimate Lyα and UV continuum emission near the center. Stacking is executed using the IRAF task "imcombine" with median/average and with/without sigma clipping to further eliminate unrelated signals. Lyα SB profiles of stacked images are measured in a series of annuli with a width of 2 pixels.

Uncertainties and Limitations
To estimate the noise level of the stacked image, we created cutouts of the Lyα and continuum image centered on randomly selected points in the field. After applying the continuum source mask, these "sky cutouts" are stacked to make stacked sky images and their radial SB profiles are measured in the same manner as stacked LAE images. This time the turning off of the masking of some sources is not employed. The 1σ noise level of each annulus is estimated by repeating this 1000 times and deriving the standard deviation of the distribution of to-tal count in the annulus 4 . As shown in Figure 3 (Top), we confirmed that the noise level decreases almost as ∝ N −1/2 . The result also demonstrates that the choice of stacking method does not make a significant difference in the noise level except for the case of average stacking without sigma clipping (in this case, too many artifacts remain; green points for "ave no clip" except for N = 3000 are far above the graph's upper bound). Thus, we decided to extrapolate this relation between the noise and N (the number of images for stacking) to estimate the noise level of stacked images with any N rather than to iterate 1000 times for every possible N in the column (4) of Table 1. The same method is used to estimate the noise level of the continuum image, which behaves almost ∝ N −1/2 as well. In figure 3 (Bottom), the noise level of the sky stack as a function of radius is shown for the N stack = 700 case. Due to the pixel spatial correlation, they do not perfectly decrease as (the number of pixels in each annulus) −0.5 . Stacking methods again do not change the result.
At the same time, the average values of the sum of the counts in each annulus were measured. Due to systematic errors and sky residuals, the average sky counts are not exactly equal to zero. To correct this effect, we subtract the average sky value when we derive the radial profile. Typical sky value of the Lyα and continuum images are ∼ −5 × 10 −21 erg s −1 cm −2 arcsec −2 and ∼ 1.2 × 10 −32 erg s −1 cm −2 Hz −1 arcsec −2 .
Since the Lyα image was created by subtracting the gband image from the NB image, any difference between the PSFs of the images could produce spurious patterns around sources in the Lyα image. Even if the simple Gaussian smoothing done in Section 2 can match the FWHM of stellar sources, it cannot exactly match the shape of the PSFs of the two images. Moreover, the shape of the PSF at a large radius may introduce additional errors. To examine the detailed shapes of the PSF in the two images, we first select bright unsaturated sources from a source catalog using the SExtractor output CLASS STAR, which is a parameter characterizing the stellarity of sources. CLASS STAR is 1 if an object is a point source and drops to 0 if extended. Here we use following criteria: CLASS STAR> 0.95 and 18 < g < 22. In total, 3980 sources are stacked to de-   . The x-axis shows how many sky cutouts were stacked, and the y-axis shows the estimated noise level evaluated in 50th annulus from the center (thus r = 100 pixel or 17 arcsec) containing 1270 pixels. Different points indicate the different stacking methods (average VS. median, with 3σ or 5σ clipping VS. no sigma clipping), which are well converged except for the case with average without sigma clipping. The red curve is a fitting function with a form a × N b to the blue circles, which is consistent with inverse square root proportionality (b = −0.493 ± 0.0074). Symbols for "ave no clip" are mostly out of the upper boundary. (Bottom) The estimated noise level of each annulus for the N stack = 700 case. Their behavior with respect to the stacking methods is the same as the above figure. termine the central part of the PSFs in the images 5 . To determine the much fainter outer part of the PSFs, we extracted stars with magnitude 13 < g SDSS < 15 from the SDSS DR14 catalog (Abolfathi et al. 2018). After excluding stars with bright nearby companions and/or obviously extended sources when seen with our deep images, images of 113 bright stars are stacked. Since point sources in this magnitude range start to saturate, the PSF of the brighter sources is connected at r = 20 pixels or 3.4 arcsec with that of fainter sources, following a method described in Infante-Sainz et al. (2019). Derived PSFs from 0.17 arcsec to 40 arcsec are shown in Figure  4. The PSF of the NB image is slightly smaller than that of the g band at r = a few arcsec. The PSFs of both bands beyond several arcsec agree very well. They are not Gaussian-like and have power-law tails with a slope of ∼ −2.8.
To check whether the slight difference between the broadband and NB PSFs affects our surface brightness measurement, we created a stacked "non-LAE" image, following a method described in Momose et al. (2014). "Non-LAE" sources are defined as objects not selected as LAEs which have almost the same distribution in the FWHM NB468 vs. NB468 magnitude plane as the real LAEs ( Figure 5). Since the majority of LAEs are distributed in the range 0.75 arcsec < FWHM < 3.25 arcsec and 24 < NB468 < 26.5, we select non-LAE sources from this range for stacking. Any signal detected in the stacked Lyα image of non-LAEs can be used to estimate the effect not only of the PSF difference but also of other unknown systematics such as errors associated with flatfielding and sky subtraction as discussed in Feldmeier et al. (2013).   Figure 6 shows the median-stacked Lyα and continuum images of all LAEs and non-LAEs without sigma clipping. The SB profiles of them are shown in Figure  7. We confirmed that the profiles do not depend on the stacking methods (except for average stacking without sigma clipping; see Figure 3); they show > 1σ deviation only at very low-S/N regime near r ∼ 100 pkpc. Hereafter, we present the results for median stacking without sigma clipping. The non-LAE has a negative ring-like structure around the center. This probably arises from the slight differences in their colors and PSFs of g-band and NB images. Still, the absolute value of the Lyα SB profile of the non-LAE is about an order of magnitude smaller than that of LAEs in Figure 7. Beyond 2 arcsec, the SB profile of the non-LAE stack is almost consistent with the sky value and thus we conclude the effect of the PSF difference is negligible, in particular at large radii of r > 2 arcsec of most interest to the present work. The PSF, which is shown with the gray curve in Figure 7, drops much more rapidly than the Lyα profile of LAEs. From the above arguments, we conclude that LAHs around LAEs at z = 2.84 are robustly detected down to ∼ 1 × 10 −20 erg s −1 cm −2 arcsec −2 and the effects of systematic errors cannot have a significant impact on the derived Lyα SB profiles out to ∼ 100 pkpc.

Stacked profiles and effects of systematics
On the other hand, the UV SB profile of LAEs seems to be negative beyond 20 pkpc. Similar patterns can be also seen in some previous studies performing stacking analyses Momose et al. 2016) but the exact reasons were never identified in the literature. When estimating the sky value with SExtractor, pixels with counts above some threshold are masked. However, since our LAEs are selected with the NB image as a detection band, some LAEs in our sample are too faint in the UV continuum image to be masked. In addition, when creating the continuum image, the Lyα contribution is subtracted even if the source is not detected or significantly affected by the sky noise in the g-band image for UV-faint sources. These effects cause oversubtraction in the continuum image of UV-faint LAEs, affecting the UV SB profiles of subsamples that contain UV-faint LAEs. In our sample, there are 430 (1391) LAEs with < 2σ (< 5σ) detection in the g-band image, respectively. The numbers of LAEs with < 2σ detection in all subsamples are also shown in Table 1. Thus, UV SB profiles of subsamples that contain many UV non-detected LAEs should be interpreted with caution 6 .
To check whether or not the mesh size for sky estimation matters, we derived a stacked Lyα profile of all LAEs with sky mesh sizes different from 30 arcsec. A larger mesh size enables us to probe possible largescale emission around LAEs, at the same time increasing   errors due to residual non-astrophysical signals (which stem from e.g. halos around bright stars). A smaller mesh size may lead to oversubtraction of the real signal while reducing the errors described above. To find a better compromise, we tested sky mesh sizes of 1 arcmin, 2 arcmin, and 11 arcsec (64 pixels, the default mesh size used in LAE selection in Section 2). The 1σ errors and residual sky emission to be subtracted were derived in the same way as described in Section 3.2. In Figure 8, we showed the results of this test. Except for the case of 11 arcsec (blue curve), the derived Lyα SB profiles are all consistent with each other within uncertainty. Also, no systematic trend is evident with increasing mesh size in the slight offset seen in the outer part. This suggests the effect of oversubtraction of diffuse Lyα emission is minor at this sensitivity, on this scale. However, as the larger sky mesh sizes are utilized, residual artificial emission around bright stars in the Lyα images becomes clearer as well. On the other hand, a mesh size of 64 pixels = 11 arcsec = 85 pkpc clearly oversubtract halo emission. Considering that LAHs are detected out to ∼ 100 pkpc, a mesh size of 85 pkpc, which is comparable to the extent of LAHs, should not be used. The same trend is seen in the continuum image as well. Thus we decided to use 0.5 arcmin = 30 arcsec sky mesh.
To quantify the extent of SB profiles, we performed fits to both UV and Lyα SB profiles using the following exponential function(s) and power-law function: where "PSF * " means convolution with the measured PSF of NB468. While exponential functions have commonly been used in previous observational work, a power-law function is motivated by an analytical model by Kakiichi & Dijkstra (2018). C 2 is set to zero for 1-component exponential fitting and let r 1 < r 2 if otherwise, thus r 1 is scale-length for the core component and r 2 is for the halo component. Unlike most previous work, we do not assume that the scale-lengths for the core component of Lyα and UV SB profile are the same, and they were fitted separately. The result shown in Figure 9 clearly demonstrates the need for non-zero C 2 or a power-law function to fit the Lyα SB profile, while a single exponential function will do for the UV SB profile fitting. The 2-component exponential fit to the UV SB profile does not converge, and the power-law fit clearly deviates from the observed UV profile.

Subsamples
The stacked Lyα and UV continuum images of all subsamples are shown in Figure 10; their corresponding Lyα and UV continuum SB profiles are shown in Figures  11 and 12, respectively with fitting curves. Rest-frame Lyα equivalent widths calculated in each annulus are also plotted with orange dots in Figure 12. The resulting fit parameters for the Lyα SB profiles are given in Table 1, and those of UV SB profiles in Table 2. In the right panel of Figure 10, the effect of oversubtraction discussed in Section 4.1 is clearly manifested by the black ring-like structures around the central emission in the stacked UV continuum images of the UV/Lyα faintest and highest-EW subsamples (rightmost three panels from top to middle). Again, the UV SB (and  the EW) profiles of UV-faint subsamples should be interpreted with caution.
We detected Lyα emission more extended than UV (stellar) emission in all subsamples (Figure 11), while most UV SB profiles can be fitted with the onecomponent exponential function ( Figure 12). However, UV SB profiles of the UV/Lyα brightest and lowest-EW subsamples clearly require two-component exponential functions and indeed can be fitted remarkably well. To our knowledge, this is the first detection of "UV halos" in high-redshift LAEs. We confirmed that this detection is robust against the choice of used stacking methods, the sky mesh size (as long as it is not too small), and the masking threshold. In addition, these subsamples are securely detected in the g-band and thus the effect of oversubtraction (Section 4.1) should be minor. For subsamples for which two-component exponential fitting is well converged, we show the resultant fitting parameters in Table 3. Figures 19 and 20 in Appendix A show results in a different manner for a clearer comparison within each photometric property. In Figure 19, clear systematic differences can be seen in bins of UV, L Lyα , and EW 0,Lyα in a way that UV/Lyα-bright LAEs and low EW 0,Lyα LAEs have larger LAHs. We see hints of UV halos also in the second and third UV brightest subsamples in the upper right panel of Figure 19. On the other hand, the difference of profiles for the projected distance and local environment subsamples is not obvious except for the protocluster subsample (those with δ > 2.5). There is no significant difference in UV SB profiles in both cases except for the protocluster subsample, which contains more UV bright galaxies as seen in Figure 2.
The scale-lengths and power-law index of fitting functions can be used for more quantitative discussion. In Figure 13, we show the results of fitting for each photometric property. r 1,UV plotted in Figures 13 and 14 are those obtained from a one-component fit (including subsamples with UV-halo), while r 1,Lyα are from a twocomponent fit. In the left three columns of the top row, both r 1,Lyα and r 1,UV show nearly monotonic behavior. While the power-law index α also shows consistent behavior, the power-law functions often deviate from the real data beyond a few tens of pkpc in Figure 11. On the other hand, r 2,Lyα behaves not that simply. This would result from both astrophysical and observational reasons as we discuss in Section 5.2. Figure 14 compares scale-lengths of UV/Lyα core/halo components. In the left panel, the scalelengths for the core components are compared. The centroids of UV continuum emission of LAEs are known to show some offset from those of Lyα emission which are defined as the image centers in this work. The vast majority have offset lower than 0.2 arcsec (Shibuya et al. 2014;Leclercq et al. 2017) which is comparable to HSC's pixel scale of 0.17 arcsec. This should not affect the measurement of r 2,UV but leads to an overestimation of r 1,UV . In addition, some subsamples contain g-band nondetected sources (see Table 1) and could be affected by oversubtraction. Although one should be aware of these potential issues, r 1,Lyα seems to be almost always larger than r 1,UV and they correlate with Spearman's rank correlation coefficient of 0.74 and pvalue of 8.2 × 10 −5 . On the other hand, similarly to the trend seen in the panels in the middle row of Figure 13, r 2,Lyα do not have a clear correlation with r 1,UV nor r 1,Lyα (p-value 0.61 and 0.65, respectively).

Sources of Differences from Previous Observational Studies
While we clearly detect LAHs more extended than UV continuum for all subsets of LAEs, some previous works report non-detections of such components (e.g., Bond et al. 2010;Feldmeier et al. 2013). We describe a number of possible reasons for such discrepancies between this study and others. The Lyα morphology cannot be properly captured by just comparing simple quantities such as their FWHMs or half-light radii without enough sensitivity or taking a very large aperture for total luminosity estimate (Nilsson et al. 2009b). Detailed analyses on SB profiles are desirable but with enough sensitivity of (at least) ∼ 10 −19 erg s −1 cm −2 arcsec −2 , and even higher sensitivity is required for safer arguments beyond a mere detection. The scale-lengths of an exponential function(s) are most widely used in the literature for such analyses (e.g., Steidel et al. 2011;Matsuda et al. 2012;Momose et al. 2014Momose et al. , 2016Wisotzki et al. 2016;Leclercq et al. 2017;Xue et al. 2017). However, resulting scale-lengths vary a lot depending not only on the data quality but also on the details of analyses: results depend on sample selection criteria, sky subtraction, masking, and stacking methods, the range of radius used for fitting, the radial binning size, fitting functions, whether to assume r 1,UV = r 1,Lyα or not, etc. First, we investigated the impact of the sensitivity using randomly chosen LAEs with a smaller sample size. We randomly took 100 LAEs and obtained r 2,Lyα of stacked Lyα images, and repeated this process 1000 times. While the medians of obtained distribution of r 2,Lyα did not show a systematic trend, the distribution had a 3.4 times larger standard deviation compared to that obtained for the 700 LAE case (see Section 5.2.1), although this number may just represent diversity in our sample. Results from lower-sensitivity data could be thus more uncertain. Secondly, most of the observed Lyα SB profiles are downwardly convex ( Figure 11) due to the flattening at ∼ 15 pkpc and thus the scale-length becomes smaller when the outer (inner) boundary of the fitting range is smaller (larger). Indeed, if we limit our fitting range up to < 30 pkpc, obtained r 2,Lyα is underestimated by 35%. The outer boundary is also affected by the sensitivity; without deep data, one has no choice other than to set it to smaller values where the signal is detected. Thirdly, as we showed in Figure 8, an insufficiently small mesh size for local sky background estimate leads to underestimation of the scale-length of LAHs, but in many cases the mesh size (or even a brief summary about sky subtraction) has not been presented in the text. As for fitting functions, a two-component fit can more robustly capture the shape of the Lyα SB profiles (see Appendix C of Xue et al. 2017 Figure 11. Radial Lyα SB profiles of all subsamples (specified in Table 1, solid orange curve with errorbars) with fitting curves. Green dashed, red dot-dashed, and gray dotted curves show the result of fitting with two-and one-component exponential functions and a power-law function. Downward triangles show 1σ error levels after residual sky subtraction. Thin gray curve shows the normalized PSF shape.
For example, Xue et al. (2017) reported a halo scalelength of LAEs of ∼ 5-9 pkpc and found no evidence for environmental dependence based on NB surveys of two protoclusters, at z = 3.78 and z = 2.66. These observations are an order of magnitude shallower than the present study; Xue et al. (2017) used the Mayall 4m telescope for the z = 3.78 data, and the Subaru telescope for the z = 2.66 data but used an intermediate band filter (IA445 on Suprime-Cam, ∆λ = 201Å) which has lower line sensitivity. The number of LAEs used to examine environmental dependence was at most 139 (for an intermediate density sample). The criteria used to select LAEs picks up relatively high EW LAEs, with EW 0,Lyα > 50Å, in the z = 2.66 protocluster field. Our study suggests that this could have biased the results toward smaller LAHs (Figure 19). Lower ionizing radia-tion field strength and/or lower abundance of cool gas in their protoclusters may also produce smaller LAHs (see Section 5.2.1 and 5.2.2). As for observations with integral field spectrographs, Wisotzki et al. (2016); Leclercq et al. (2017) probed r < 30 pkpc of LAEs at z > 3 and obtained r 2,Lyα of 15 pkpc with the majority with r 2,Lyα < 5 pkpc. Analyses on an individual basis would suffer from greater noise and sample variance than ours. On the other hand, Chen et al. (2021) probed 59 star-forming galaxies (including non-LAEs) at z = 2-3 with sufficiently deep Keck/KCWI observational data and conducted stacking. They got r 1,Lyα = 3.71 +0.06 −0.04 pkpc and r 2,Lyα = 15.6 +0.5 −0.4 pkpc for their stacked Lyα SB profiles. Their sample is typically about an order of magnitude more massive and star-forming than our sample and this could explain the size difference. To    Table 1, solid blue curve with errorbars) with fitting curves. Green dashed and red dot-dashed curves show the result of fitting with two-and one-component exponential functions. The gray dashed curve is Lyα SB profiles converted from FLyα units to fν units by simply dividing by the FWHM of the NB filter. Orange dot indicates rest-frame Lyα equivalent width inÅcalculated in each annulus, with its value on the right axis. Downward and upward triangles show 1σ limits of UV emission and equivalent width (some datapoints in the rightmost panels are above the upper boundary) after residual sky subtraction. Thin gray curve shows the normalized PSF shape. Note the difference in the range of the x-axis from that of Figure 11. summarize, comparing results obtained using inhomogeneous analyses in the literature without consistent reanalysis is difficult, and thus we will not attempt a detailed comparison here.

Dependence of Scale-length on Galaxy Properties
In Figure 13, UV/Lyα-bright or low-EW LAEs tend to have larger r 1,Lyα and r 1,U V than UV/Lyα-faint or high-EW LAEs. This is qualitatively consistent with the HST-based results of Leclercq et al. (2017) although our results from a ground-based telescope tend to show larger values (0.1-1 pkpc vs > 1 pkpc). Considering that UV/Lyα luminous and/or low-EW LAEs tend to be more massive, the trend is also consistent with the known trend between M UV and effective radius in the UV (e.g., Shibuya et al. 2019), or more generally the so-called size-luminosity or size-mass relation. A larger scale-length in massive LAEs also results from suppression of UV/Lyα light due to more abundant dust especially in the central region (Laursen et al. 2009) than in less massive LAEs. Lyα photons are then further affected by resonant scattering and differential extinction caused thereby, leading to r 1,Lyα > r 1,UV .
On the other hand, r 2,Lyα do not show a simple trend with respect to any photometric properties. This is again almost consistent with Leclercq et al. (2017). This fact can be attributed to both observational and astrophysical reasons. First, we may simply lack the sensitivity to reveal a real trend. Even with our deep images, the Lyα SB profiles in Figure 11 at r > 50 pkpc are not well constrained. In addition, the astrophysics involved in diffuse emission in LAHs is notoriously complicated as we discuss in Section 5.4. Dominant mechanisms for Lyα production might be different over different mass/luminosity ranges, making a simple trend (if any) difficult to discern. For example, as compiled in Kusakabe et al. (2019, their Figure 7), the total Lyα luminosity of the halo component may depend in a different way on halo mass with respect to production mechanisms (e.g., collisional excitation in cold streams vs. scattering). Lastly, both observations (Wisotzki et al. 2016;Leclercq et al. 2017) and numerical studies (Lake et al. 2015;Byrohl et al. 2021) have shown that the Lyα SB profiles of individual LAEs are very diverse, even within galaxies with similar integrated properties. To conclude whether there is any trend, even larger samples and/or deeper observations are needed.
With better data, we could more easily select whether two-component exponential functions or the power-law functions are preferred.

Curious Behavior of Distance Subsamples: QSO
Radiative History Imprinted?  Subsamples based on the projected distance from the HLQSO (d Q ) show a significant variation with a minimum at d Q ∼ 10 pMpc ( Figure 13, the second panel from right in the second row). This could just be due to the stochasticity discussed above, but if real, it could be related to the QSO's radiative history. r 2,Lyα of the two largest (r 2,Lyα > 13 pkpc, d Q < 6.2 pMpc and 14.8 < d Q < 16.9 pMpc subsample) and the smallest (r 2,Lyα < 7 pkpc, 9.5 < d Q < 12 pMpc subsample) differ significantly; when we repeatedly select 700 LAEs at random from the whole sample, stacking their Lyα images, and measuring r 2,Lyα 1000 times, we get r 2,Lyα < 7 pkpc 0.8% of the time and r 2,Lyα > 13 pkpc 12% of the time (with the median value of r 2,Lyα = 10.5 pkpc) 7 . If the HLQSO was active ∼ 50 Myrs ago, followed by ∼ 30 Myrs of inactivity, and was re-ignited 20 Myrs ago, the ionizing photons emitted by the QSO would have traveled a distance of > 15 pMpc and < 6 pMpc from the QSO 8 . These photons can ionize the envelopes of the LAEs and boost their Lyα luminosity, explaining the observed behavior 9 . Assuming the HLQSO 50 Myrs ago had the same luminosity as it has today ( luminosity near a rest-frame wavelength 1450Å, νL ν,1450 = 5.7 × 10 47 erg s −1 , Trainor & Steidel 2012) and isotropic radiation with escape fraction of unity, the ionizing radiation at 16 pMpc from the QSO can still dominate over the cosmic average UV background at z ∼ 3, Γ z=3 bkg = 1.0 × 10 −12 s −1 (Becker & Bolton 2013) by a factor of a few. Cantalupo et al. (2005) calculated the fluorescent Lyα emission due to QSOs in addition to the cosmic background and gave a fitting formula for an effective boost factor b eff (their Equation 14-16) which can be used to estimate SB of illuminated gas clouds. In our case, the resulting SB would be SB = (0.74 + 0.50(11.5(r/16 pMpc) −2 (1.6 α /α)) 0.89 )SB HM , where r is the distance from the HLQSO in pMpc, α is the QSO's spectral slope (L ν ∝ ν −α ), and SB HM = 3.67 × 10 −20 erg s −1 cm −2 arcsec −2 is the expected SB without QSO boost. Assuming α ∼ 1 gives SB = 2.7 × 7 If we exclude the 55 protocluster LAEs from the d Q < 6.2 pMpc subsample, r 2,Lyα becomes 11.3 pkpc (the ∼ 30 percentile). 8 These time estimates are lower limits calculated with projected distance and the speed of light. Propagation of ionization fronts could be delayed in some situations (Shapiro et al. 2006). 9 see also Trainor & Steidel (2013); Borisova et al. (2016) where the authors used QSOs associated with spectroscopic high-EW(> 240Å) LAEs to place limits on QSO lifetime. We see consistent behavior also in the fraction of high EW LAEs as a function of distance from the HLQSO (Kida et al. 2019). Those EW are derived with 1. 5 aperture. However, boosting EW of such a central part of LAEs with L Lyα > 10 41 erg s −1 at 16 pMpc distance would be energetically not feasible with the current HLQSO luminosity.
10 −19 erg s −1 cm −2 arcsec −2 . Thus, QSO-induced fluorescence is energetically possible to have caused variation in r 2,Lyα in the projected distance subsamples, at least in the optimistic case. In reality, our narrow-band selection picks up LAEs with line-of-sight distance uncertainty of ∼ 19 pMpc, which is the same level as the radius of the FoV of our images, and this randomizes light-travel time from the QSO to each LAE. Such effects further complicate the situation, but we have shown that under some circumstances with appropriate QSO light curve and line of sight distribution of LAEs, the observed trend of r 2,Lyα might be explained. Upcoming instruments such as Prime Focus Spectrograph (PFS; a wide-field multi-fiber spectrograph, Takada et al. 2014) on the Subaru Telescope can test this fluorescence scenario by obtaining systemic redshifts of LAEs and thus reducing uncertainties on their 3D distances from the QSO.

Origin of the Large LAH of Protocluster LAEs
We showed that the dependence of LAHs on environment is not large where δ < 2.5, but in the protocluster (δ > 2.5) the LAHs show elevated flux out to > 100 pkpc (Figure 13,20). If LAHs trace the gas distribution of the CGM, the former indicates that the large-scale environment (except for protocluster environments) does not have a large impact on the matter distribution out to ∼ 100 pkpc and it is determined by individual halo mass or other internal processes. Alternatively, LAEs could simply be poor tracers of large-scale environments. Recently, Momose et al. (2021) found that LAEs behind the foreground large-scale structure tend to be missed due to absorption by the foreground structure (see also Shimakawa et al. 2017). It may also be the case that even our large field of view of 1.2 deg diameter may be insufficient to capture diverse environments including voids while targeting a single protocluster. Other line emitters or continuum-selected galaxies would be ideal, though would be expensive for current facilities. On the other hand, the HS1549 protocluster is confirmed by overdensity of LAEs and continuum-selected galaxies, with ∼ 200 member galaxies spectroscopically identified Mostardi et al. 2013;C. Steidel et al. in prep.). Enlarging the size of cutout images, we confirmed that flux higher than > 1σ level continues out to ∼ 500 pkpc. In the following, we investigate the possible cause of this emission in the δ > 2.5 subsample.
First, we further divide the protocluster sample into a "core" group (LAEs within a projected distance of < 500 pkpc from the HLQSO but excluding the HLQSO itself, N = 25) and an "outskirt" group (the remainder, N = 29), and stacked them separately. This time, pix-els with SB> 10 −17 erg s −1 cm −2 arcsec −2 around the HLQSO and the associated bright nebula ) are masked before stacking in each cutout image to exclude their contribution. As shown in Figure 15, the core sample (the yellow curve) clearly shows excess emission which does not decrease toward the edge above the original orange curve, while the outskirt stack (the green curve) shows just a mild bump around ∼ 15 pkpc and no evidence of excess. This suggests the extended emission of δ > 2.5 sample is solely produced by the LAEs at the core of the protocluster.
Recalling that the protocluster LAEs tend to have higher M UV , L Lyα and lower EW ( Figure 2) and that such galaxies generally have more extended Lyα SB profiles ( Figure 13 and Figure 19), the overabundance of such LAEs leads to larger LAHs. Moreover, the core region is sufficiently crowded that Lyα emission from neighboring LAEs may overlap, thereby leading to an overestimation of the SB profile. We evaluated this effect using the stacked Lyα images of subsamples based on UV magnitude (those shown in the top row of Figure  10) as follows. We embedded their images in a blank image, mimicking the observed spatial distribution of LAEs (including the HLQSO) using IRAF task "mkobjects". When embedding LAEs with a certain M UV , we used the stacked Lyα image of the appropriate M UV subsample scaled to match the observed UV magnitude. For example, we embedded the scaled stacked image of the −19.2 < M UV < −18.6 subsample at the location of LAEs in the same UV magnitude range. Cutout images of the simulated Lyα image were then created at the locations of the embedded LAEs and stacked in the same manner as the real "core" LAEs. The result is also shown in Figure 15 (the blue curve). The observed large LAH (the orange curve) is not reproduced. Thus, we conclude that the combined effect of the overabundance of bright LAEs and overlapping do contribute to the large LAHs, but cannot fully explain the extent of the LAH of the core LAEs. Kikuta et al. (2019) showed that there is a Mpc-scale diffuse Lyα emitting structure around the HLQSO. Such diffuse emission would explain the remaining excess in the core region of the protocluster. The excess comes not only from gas directly associated with LAEs but also from gas out of the LAEs; the pixel value distribution of the core region (< 500 pkpc from the HLQSO) after masking < 10 arcsec (∼ 80 pkpc) regions around the detected LAEs and a ∼ 200 pkpc ×260 pkpc box covering the bright QSO nebula clearly shows excess at 1 × 10 −18 <SB< 4 × 10 −18 erg s −1 cm −2 arcsec −2 compared to that of the outer region (with regions around LAEs masked in the same manner as the core case).
Previously, two studies reported similar large diffuse LAHs in protoclusters at z = 2-3; Steidel et al. (2011) observed HS1549 (z = 2.84, same as this study 10 ), HS1700 (z = 2.30; see also Erb et al. 2011), and SSA22 (z = 3.09; also observed by Matsuda et al. 2012). For the HS1549 and SSA22, direct evidence of Mpc-scale diffuse Lyα emission has been observed Umehata et al. 2019, respectively), and the identification of filamentary structure traced by six Lyα blobs in the HS1700 protocluster by Erb et al. (2011) suggests that this protocluster also harbors such diffuse emission. In the forming protocluster core at z = 2-3, a large amount of cold gas can be accreted through the filamentary structure (the cosmic web) penetrating the cores (Kereš et al. 2005(Kereš et al. , 2009Dekel et al. 2009). In addition to abundant gas, Umehata et al. (2019) showed that an enhanced ionizing UV background due to a local overdensity of star-forming galaxies and AGNs may play a crucial role in boosting fluorescent Lyα emission to a detectable level. The HS1549 protocluster also has tens of active sources (i.e., AGN, Lyα blobs (Kikuta et al. 2019, C. Steidel et al. in prep.), and SMGs (Lacaille et al. 2019)) within a few arcmins from the HLQSO, producing sufficient UV radiation to power the diffuse emission; in Supplementary Material S9 of Umehata et al. (2019), they calculated the required UVB strength to boost SB of optically thick gas. Fifty times stronger UVB than the cosmic average at z = 3, which is easily realized by the HLQSO alone in the area within several pMpc from it, would boost SB to ∼ 10 −18 erg s −1 cm −2 arcsec −2 level. To summarize, the very large LAHs of stacked Lyα profiles reported previously and in this work can be attributed to an overlap of crowded LAEs and diffuse fluorescent Lyα emission within the forming protocluster core. The δ 2 dependence of LAH scale-length claimed in Matsuda et al. (2012) should be revisited using new data targeting more protocluster fields at similar redshift together with appropriate analysis methods as discussed in Section 5.1.

Discovery of "UV Halos" and Its Implication to Low-Mass Galaxy Evolution
As seen in Figure 12 and Figure 13 in Section 4.2, we have discovered "UV halos" around UV/Lyα-bright and/or low-EW LAEs. This has a significant impact on our understanding of the origin of LAHs, and also of galaxy evolution, because it provides direct evidence of 10 Although their LBG sample is an order of magnitude brighter in UV than our LAE sample, we did a consistency check with Steidel et al. (2011) results; we confirmed that our stack of LBGs used also in their sample gives a consistent Lyα SB profile as their work. 2.5 < δ core w/o QSO outskirt emb All Figure 15. Stacked Lyα SB profiles of the protocluster (δ > 2.5) sample (orange), the core sample (yellow), the outskirt sample (green), embedded and stacked image (blue, see text), and all LAEs (purple, same as Figure 7). Downward triangles show 1σ error levels after residual sky subtraction.
star formation activity in the outskirts of high-redshift low-mass galaxies.
To gain more insight on the latter point, we used the data products from the TNG100 run of the IllustrisTNG simulation (e.g., Nelson et al. 2019;Pillepich et al. 2018). We make median stacked SFR surface density profiles of FOF (friends-of-friends) halos, which roughly represent collections of gravitationally bound DM particles, at z = 3 for 4 SFR bins; 0.1 < SFR < 1 M yr −1 , 1 < SFR < 10 M yr −1 , 10 < SFR < 100 M yr −1 , and 100 M yr −1 < SFR and compare them with those of our three UV-brightest subsamples after converting the simulation data to UV flux density using the SFR-UV luminosity density conversion of Murphy et al. (2011);Kennicutt & Evans (2012) relation (Figure 16). The SFR surface density profiles were convolved with a Gaussian with 0. 77 FWHM 11 to make a fair comparison. Here, SFR of TNG galaxies denotes a total of all particles which belong to one FOF halo. We also plot the prediction of Lake et al. (2015, convolved with a Gaussian with 1. 32 FWHM or 10.3 pkpc at z = 3.1) for discussion in Section 5.4. The three UV-bright subsamples have median UV absolute magnitude of −19.62, −18.88, −18.31 (Table 1), but these were derived using 1. 5 diameter apertures and thus may be underestimated. The total magnitudes derived by integrating the UV SB profiles down to the radii where emission is detected at more than 1σ significance are −21.34, −20.49, −20.10, respectively, corresponding to SFR of 17, 7.7, 5.4 M yr −1 . The SFR of simulated galaxies whose UV SB profiles match those of our LAEs is higher than that of our LAEs (10-100 vs. 5.4-17), but there is considerable uncertainty in the conversion between UV luminosity and SFR, as it depends on stellar age, dust attenuation, and metal abundance, and not all galaxies in simulations would be selected as LAEs. Reconciling this mismatch is beyond the scope of our paper. Rather, to compare the profile shapes we normalized each curve in the bottom panel of Figure 16. While the UV SB profiles of the two fainter subsamples seem to be slightly more compact than the SFR density profiles of TNG galaxies, the UV-brightest subsample has a remarkably similar shape as the SFR surface density profiles of TNG galaxies with 1 < SFR < 10 and 10 < SFR < 100 subsamples.
We further decompose the simulated FOF halos into the main halo and subhalos. The decomposed SFR profiles demonstrate that the flattened outer part is dominated by the contribution of satellites. Given the similarity of the profiles, we suggest that the UV halo of the UV-brightest LAEs is also due to such satellites. To characterize the subhalos around the central galaxy, we extracted dark matter halo mass, gas mass, stellar mass, and SFR of those within 50 pkpc (∼ 6 arcsec, 2D projected distance) from their main halos. We only handle those halos with M DM > 7.5 × 10 7 M (= 10×DM particle mass), M stellar > 7.0 × 10 6 M (= 5×stellar particle mass), and non-zero SFR to avoid spurious objects. The distribution of DM mass, stellar mass, gas mass, and SFR of satellites are plotted in Figure 17. Satellites responsible for the UV halo of our LAEs would be similar to those around central galaxies with 1 < SFR < 10 and 10 < SFR < 100. On average, they have 1.9 and 2.3 satellites, respectively, with median DM halo masses of 3.3×10 9 M and 4.4×10 9 M , and mean total halo SFR of 0.30 and 2.6 M yr −1 (i.e., ∼ 10% of central galaxies). This suggests that the UV halo is comprised of a few satellite galaxies around the main halo and not by an intrinsically diffuse halo, unlike optical stellar halos of local galaxies (e.g., D'Souza et al. 2014). Under this hypothesis, individual LAEs would not have smooth UV SB profiles like those presented in Figure 12, but would be more likely to exhibit stochastic shapes made by ∼ 2 discrete satellites; thus the term UV "halo" may not be appropriate, if the quoted simulations represent the reality. The smooth profile of stacked UV SB is simply a reflection of the radial and SFR distribution of satellite galaxies.
Are such satellites observable? The SFR distributions of satellites around galaxies with 0.1 < SFR < 10 in the simulations have medians in the range 0.01 < SFR/[M yr −1 ] < 0.1, although some very low-mass objects may be affected by resolution effects (Figure 17). This translates into UV absolute magnitude of −13.3 and −15.8, and apparent g-band magnitude of 32.1 to 29.6 mag at z = 2.84 (assuming no K-correction, or equivalently, flat UV continua). Brighter satellites are thus well within the reach of HST and JWST sensitivity in some deep fields or with the aid of gravitational lensing (e.g., Alavi et al. 2016;Bouwens et al. 2022 Lastly, we infer the origin of LAHs of LAEs through comparison with recent numerical simulations. As introduced in Section 1, Lyα photons in the halo regions are generated either by ex-situ (mostly from the host galaxy) and transported by resonant scattering in neutral gas, or in-situ via photoionization followed by recombination or collisional excitation. In-situ photoionization is maintained by ionizing photons from star formation and/or AGN activity in satellites, central galaxies, or other nearby sources of ionizing UV such as QSOs.
In-situ collisional excitation would be driven by shocks due to fast outflows due to feedback or gravitational energy of inflowing gas, but the former is considered to be effective in more massive and energetic sources such as radio galaxies and Lyα blobs (e.g., Mori et al. 2004). In-situ Lyα photons may also experience resonant scattering, but its effect on the redistribution of photons is likely to be relatively weak due to lower HI column densities than the central part by more than a few dex (Hummels et al. 2019;van de Voort et al. 2019).
To predict Lyα SB profiles around galaxies, one needs to know physical parameters such as hydrogen density, neutral fraction, gas kinematics, and temperature which depend on SFR and AGN activity around a point of interest, and ionizing/Lyα photon escape fraction, etc. Solving all these quantities is practically impossible, but recent simulations are beginning to reproduce observed Lyα SB profiles reasonably well. Among such studies is Byrohl et al. (2021, hereafter B21); B21 presents full radiative transfer calculations via post-processing of thousands of galaxies in the stellar mass range 8.0 < log(M * /M ) < 10.5 drawn from the TNG50 simulations of the IllustrisTNG project. One of the advantages of B21 is the sample size, which is far larger than those of previous studies. For example, Lake et al. (2015, hereafter L15) performed radiative transfer modeling of 9 LAEs, obtaining significantly diverse SB profiles. Those with massive neighbors have elevated SB profiles both in UV and Lyα, significantly boosting the average profile. But the 9 galaxies modeled in the simulation would not be representative of star-forming galaxies at z ∼ 3 if not carefully selected. Mitchell et al. (2021) calculated Lyα SB profiles of a single galaxy at z = 3-4 using the RAMSES-RT code, a radiation hydrodynamics extension of the RAMSES code, and succeeded in reproducing SB profiles similar to MUSE observations. The small sample size is somewhat mitigated by using all available outputs between z = 4 and z = 3, but still there may remain biases with respect to environment or evolutionary phase. For these reasons, we compare our results primarily with those of B21.

Satellite Galaxies
In Figure 18, we plot B21 (convolved with a Gaussian with 0. 7 FWHM) and L15 (convolved with a Gaussian with 1. 32 FWHM) results with our observations of subsamples based on UV magnitude. Although B21 probed only up to < 50 pkpc, the overall shapes of the 9.5 < log(M * /M ) < 10.0 stack and our UV-brightest subsample stack match remarkably well, and the 9.0 < log(M * /M ) < 9.5 stack and the other UV subsamples stack also show fairly good agreement. We also highlight the diversity of Lyα SB profiles here by drawing all the profiles of galaxies with 9.0 < log(M * /M ) < 9.5 in B21 in Figure 18 with thin curves. Similarly to L15, bumps in some curves are caused by companion galaxies. This demonstrates the difficulty of studying halo origins with a small sample as discussed in Section 5.2. Future (observational and theoretical) studies should keep this in mind before discussing the halo dependence on physical properties. B21 concluded that scattering of Lyα photons originate from star formation in the central and satellite galaxies is almost always dominant within 50 pkpc from the center (∼ 50% at > 20 pkpc), with insitu collisions and recombination explaining the remaining 30% and 20%, respectively. The contribution from satellite star formation dominates over that from the central galaxy beyond ∼ 40 pkpc, and they show that halos that have more massive neighbors within 500 pkpc can have very extended LAHs compared to those resid-ing in normal environments (Figure 12 of B21). Similar conclusions are reached by L15 and Mitchell et al. (2021) as to the importance of satellites. The dashed curve shows the L15 result which extends to 100 pkpc. As we saw in Figure 16 (see also Figure 4 of L15), their galaxies (whose mean stellar mass is ∼ 2.9 × 10 10 M ) have more star formation outside of the host halos and have an enhanced SB profile at outer regions. A comparably massive galaxy sample is required to confirm their prediction. Our first detection of satellites (Section 5.3) and reasonable agreement of both UV ( Figure 16) and Lyα ( Figure 18) SB profiles with simulations suggest that star formation in central and satellite galaxies are important Lyα sources contributing to LAHs.
Rest-frame Lyα EW of each annulus of observed LAEs is plotted with the right axis of Figure 12. Considering the outward diffusion of Lyα from the central galaxies, this EW is always an upper limit for the EW of in-situ Lyα emission. If low-mass satellites are responsible for the outer LAHs as simulations suggest, then their expected dark matter halo masses are about 10 9−10 M . EW 0,Lyα of 200Å observed at r ∼ 30 pkpc of the UV/Lyα-brightest or EW-lowest subsamples can be explained by halo star formation alone if these low-mass galaxies have average EW 0,Lyα of ≥ 200Å, even without scattered Lyα from central galaxies. For other subsamples, the EW 0,Lyα is > 240Å due to the faintness of UV SB and the extended Lyα SB profiles at outer regions. Such high EW 0,Lyα is hard to explain by star formation alone (Schaerer 2003); scattering of Lyα photons produced elsewhere and in-situ recombination/collisional excitation should dominate the Lyα photon budget.

Non-Negligible Contribution from In-Situ Lyα Production
Processes other than star formation, e.g. QSOboosted fluorescence or collisional excitation via gravitational cooling, can still be important at large radii, since they are predicted to make non-negligible contributions and there are situations under which such processes become more important. In Section 5.2.1 and Section 5.2.2, we specified conditions where fluorescence could be dominant, namely regions near bright QSOs and/or near protocluster cores. A number of simulations have suggested that stronger ionizing radiation fields boost fluorescent Lyα emission from the CGM/IGM (Cantalupo et al. 2005;Kollmeier et al. 2010, see also Appendix A5 of B21). A major problem is that both processes are very hard to accurately predict even with state-of-the-art simulations; recombination emissivity could be significantly boosted without changing total hydrogen column density if there are many tiny ( 1 pkpc) clumps with locally increased density in the CGM/IGM regions unresolved in current standard simulations but suggested from observations (Rauch et al. 1999;Cantalupo et al. 2014Cantalupo et al. , 2019Mc-Court et al. 2018;Hummels et al. 2019;van de Voort et al. 2019). The total Lyα luminosity from gravitational cooling should have a strong dependence on halo mass (Goerdt et al. 2010). In addition, the emissivity of the collisional process depends extremely sensitively on temperature exponentially in the range T = 10 4 -10 5 K characteristic of cold accretion, and treatments of the effect of self-shielding against the UVB may have a critical impact on results (Faucher-Giguère et al. 2010;Kollmeier et al. 2010;Rosdahl & Blaizot 2012).
There remains a possibility that our main conclusion about the dominance of central/satellite star formation may apply only to relatively massive LAEs, since lowermass halos would have less scattering media and less satellites. For example, high-EW LAEs are efficient producers of Lyα photons. But because they are on average less massive and should have lower HI gas (Rakic et al. 2013;Turner et al. 2017), they have more com-pact LAHs despite efficient Lyα production. Kakiichi & Dijkstra (2018) showed that scattering of Lyα photons produced by central galaxies with realistic HI distributions constrained by Lyα forest observations of LBGs results in power-law-like Lyα SB profiles (see also Steidel et al. 2011). In Figure 11, UV/Lyα-faint and high-EW LAEs seem to deviate from power-law fits at r > 25 kpkc. This could be a hint of the dominance of other processes. Similar conclusion of the dominance of scattered Lyα from central galaxies and possible contribution from the other processes is reached by recent observational studies (Lujan Niemeyer et al. 2022a,b). In this way, much information is buried beyond the "flattening radius" around 20 pkpc, outside of which contributions from central/satellites appear insufficient to explain the observations. With a larger sample and deeper data, we can further probe the behavior of LAHs e.g., by making EW-based subsamples with matched UV magnitude, etc.

SUMMARY
We have investigated the rest-frame UV continuum (λ ∼ 1225Å) and Lyα radial surface brightness (SB) profiles of LAEs at z = 2.84 through stacking analyses of UV and Lyα images created from Subaru/HSC gband and the NB468 narrow-band images. The depth and wide field coverage, including a known protocluster, enable us to study both SB profiles with unprecedented depth because of the large sample (N = 3490) at z ∼ 3. Our major findings are as follows: 1. Stacking of 3490 LAEs yields a SB sensitivity of ∼ 1 × 10 −20 erg s −1 cm −2 arcsec −2 in Lyα and ∼ 1 × 10 −33 erg s −1 cm −2 Hz −1 arcsec −2 (Figure 9). Our analyses reveal that systematic errors should be at the same levels at most (Section 4.1 and that the choice of mesh size for local sky estimation could have a large impact on the results (Figures  7 and 8).
2. By dividing the LAEs into subsamples according to various photometric properties, UV magnitude, Lyα luminosity, Lyα EW, projected distance from a hyperluminous (HL) QSO residing at the center of the protocluster (as a proxy for radiation field strength boosted by the HLQSO), and LAE overdensity δ on a 1. 8 (∼ 840 pkpc) scale, we study the dependence of the SB profiles on these quantities (Figures 11 and 12). To quantify the radial dependence of SB profiles, we fit 2-component exponential functions (Equation 1) to observed profiles. For Lyα SB profiles, we consistently obtain r 2,Lyα , the scale-length of the more extended component, of ∼ 10 pkpc for all subsamples. However, we do not observe any clear trend of r 2,Lyα with any property probed here (Figure 13), whereas the scale-length of the compact components (both of r 1,Lyα and r 1,UV ) varies monotonically with respect to UV magnitude, Lyα luminosity, and Lyα EW.
3. We find an exceptionally large exponential scalelength r 2,Lyα for LAEs in the inner core (those within 500 pkpc from the HLQSO) of the protocluster, and a significant variation in r 2,Lyα with respect to the projected distance from the HLQSO. These findings could be explained by enhanced ionizing background radiation due to abundant active sources and cool gas at the forming protocluster core and the past activity of HLQSO, respectively (Section 5.2.1 and Section 5.2.2).
4. We for the first time identify extended UV components (i.e., r 2,UV inconsistent with zero), or "UV halos" around some bright LAE subsamples, which provides direct evidence for the contribution of star formation in halo regions and/or satellite galaxies to Lyα halos. Comparison with cosmological hydrodynamical simulations suggests that UV halos could be composed of 1-2 low-mass (M DM ∼ 10 9.5 M ) galaxies (Section 5.3) with total SFR of ∼ 10% of that of their central galaxies.
5. Combining our results with predictions of recent numerical simulations, we conclude that star formation in both the central galaxy and in satellites, together with resonant scattering of Lyα photons are the dominant factors determining the Lyα SB profiles at least within a few tens pkpc. In outer regions (projected distances 30 pkpc) other mechanisms such as fluorescence can also play a role especially in certain situations like dense regions of the Universe and near zones of bright QSOs (Section 5.4).
The low-mass satellite galaxies suggested by our deep UV stacked images will be very important targets for revealing the role of minor mergers in galaxy evolution and cosmic reionization, as they are believed to be promising analogues of the main galaxy contributors to ionizing photon budget at z > 6 ( Robertson et al. 2013;Mason et al. 2015). In a year or so, deep surveys with JWST will detect these galaxies within several arcsecs from LAE-class galaxies, which are to be observed in coming programs. Finally, Hα observations of LAHs open up a new pathway to study star formation in the CGM and fluorescent clumps without the blurring effect of resonant scattering. Observations of Hα from z = 2.84 has not been possible with groundbased telescopes due to heavy atmospheric absorption and extremely bright backgrounds, but now this is also becoming feasible thanks to JWST. In the future, we will combine cross-analyses with a larger LAE sample and new constraints from JWST to further constrain the origin of LAHs.

ACKNOWLEDGMENTS
We are grateful to the anonymous referee for their careful reading and constructive comments and suggestions, Chris Byrohl for providing simulation data, and Haibin Zhang, Masafumi Yagi, Masayuki Umemura, Tadafumi Takata, Kazuhiro Shimasaku, Yusei Koyama, and Kazuhiro Hada for fruitful discussions. We thank Yukie Oishi and the HSC pipeline team for their helpful comments on HSC data analyses. We would like to acknowledge all who supported our observations at the Subaru Telescope, including the staff of the National Astronomical Observatory of Japan, Maunakea Observatories, and the local Hawaiian people who have been making efforts to preserve and share the beautiful dark sky of Maunakea with us. We are honored and grateful for the opportunity of observing the Universe from Maunakea, which has the cultural, historical, and natural significance in Hawaii. Data analysis was carried out on the Multi-wavelength Data Analysis System operated by the Astronomy Data Center (ADC), National Astronomical Observatory of Japan and on analysis servers at Center for Computational Astrophysics, National Astronomical Observatory of Japan. SK acknowledges supports from the Japan Society for the Promotion of Science (JSPS) KAKENHI grant Nos. 18J11477, 19H00697 Figures 19 and 20, we show Lyα (left) and UV (right) SB profiles of each subsample in each row for easier comparison.
In the top-left panel of Figure 20, the curves appear to deviate from each other beyond 2 arcsec, although the sensitivity at these angular separations is not very high. To check whether this difference is significant, we stacked Lyα images of 700 randomly selected LAEs (roughly corresponding to the number of LAEs in each bin of projected distance subsample) from the whole (N = 3490) LAE sample and derived its Lyα SB profile. We repeated this 500 times and derived the 5th and 95th percentile of the SB distribution in each annulus. These are shown as a gray shaded region in Figure 20; in the bottom panel, we also show the SB distribution with 1000 randomly selected LAEs to see the difference between δ < 1.0 subsamples. The curves are almost within the shaded regions, suggesting that the difference is apparently marginal.  EW Lyα,0 < 30A 30A < EW Lyα,0 < 55A 55A < EW Lyα,0 < 90A 90A < EW Lyα,0 < 160A 160A < EW Lyα,0 Figure 19. Lyα SB profiles (Left) and UV continuum SB profiles (Right) derived by stacking analysis. From top to bottom, the LAE sample is grouped by their UV magnitude, Lyα luminosity, and Lyα equivalent width in the manner specified in Table 1x10  Dist < 6.2 pMpc 6.2 pMpc < Dist < 9.5 pMpc 9.5 pMpc < Dist < 12.0 pMpc 12.0 pMpc < Dist < 14. 2.5 < δ 1.0 < δ < 2.5 0.3 < δ < 1.0 -0.15 < δ < 0.3 -1.0 < δ < -0.15 Figure 20. Same as Figure 19, but here the LAEs are grouped by their projected distances from the HLQSO (Top) and local environment (Bottom) in the manner specified in Table 1. The gray shaded regions in the Top (Bottom) panel shows 5th and 95th percentile of the Lyα SB distribution of stacked images created with randomly selected 700 (1000) LAEs.