Clustering dependence on Lyman-$\alpha$ luminosity from MUSE surveys at $3

[Abbreviated] We investigate the dependence of Lyman-$\alpha$ emitter (LAE) clustering on Lyman-$\alpha$ luminosity. We use 1030 LAEs from the MUSE-Wide survey, 679 LAEs from MUSE-Deep, and 367 LAEs from the to-date deepest ever spectroscopic survey, the MUSE Extremely Deep Field. All objects have spectroscopic redshifts of $3<z<6$ and cover a large dynamic range of Ly$\alpha$ luminosities: $40.15<\log (L_{\rm{Ly}\alpha}/[\rm{erg \:s}^{-1}])<43.35$. We apply the Adelberger et al. K-estimator as the clustering statistic and fit the measurements with state-of-the-art halo occupation distribution (HOD) models. From the three main data sets, we find that the large-scale bias factor, the minimum mass to host one central LAE, $M_{\rm{min}}$, and (on average) one satellite LAE, $M_1$, increase weakly with an increasing line luminosity. The satellite fractions are $\lesssim10$% ($\lesssim20$%) at $1\sigma$ ($3\sigma$) confidence level, supporting a scenario in which DMHs typically host one single LAE. We next bisected the three main samples into disjoint subsets to thoroughly explore the dependence of the clustering properties on $L_{\rm{Ly}\alpha}$. We report a strong ($8\sigma$) clustering dependence on $L_{\rm{Ly}\alpha}$, where the highest luminosity LAE subsample ($\log(L_{\rm{Ly}\alpha}/[\rm{erg \:s}^{-1}])\approx42.53$) clusters more strongly ($b_{\rm{high}}=3.13^{+0.08}_{-0.15}$) and resides in more massive DMHs ($\log(M_{\rm{h}}/[h^{-1}\rm{M}_{\odot}])=11.43^{+0.04}_{-0.10}$) than the lowest luminosity one ($\log(L_{\rm{Ly}\alpha}/[\rm{erg \:s}^{-1}])\approx40.97$), which presents a bias of $b_{\rm{low}}=1.79^{+0.08}_{-0.06}$ and occupies $\log(M_{\rm{h}}/[h^{-1}\rm{M}_{\odot}])=10.00^{+0.12}_{-0.09}$ halos. We discuss the implications of these results for evolving Ly$\alpha$ luminosity functions, halo mass dependent Ly$\alpha$ escape fractions, and incomplete reionization signatures.


Introduction
Dark matter halos (DMHs) serve as sites of galaxy formation but their co-evolution is still a matter of investigation. Observations deliver snapshots of the luminosities of galaxies at given redshifts, while numerical analyses succeed at simulating the evolution and copiousness of DMHs. Linking these two constituents is not straightforward but, because the spatial distribution of baryonic matter is biased against that of dark matter (DM), the former indirectly traces the latter. The evolutionary stage of the two distributions depends on both the epoch of galaxy formation and the physical properties of galaxies (see Wechsler & Tinker 2018 for a review). Thus, studying the dependence of the baryonic-DM relation on galaxy properties is essential for better understanding the evolution of the two components.
Exploring the spatial distribution of high-redshift (z > 2) galaxies and its dependence on physical properties provides an insight into the early formation and evolution of the galaxies we observe today. Clustering statistics yield observational constraints on the relationship between galaxies and DMHs, as well as on their evolution. Traditional studies of high-z galaxies (Steidel et al. 1996;Hu et al. 1998;Ouchi et al. 2003;Gawiser et al. 2007;Ouchi et al. 2010;Khostovan et al. 2019) model the large-scale (R 1 − 2 h −1 cMpc) clustering statistics with a two parameter power-law correlation function that takes the form ξ = (r/r 0 ) −γ (Davis & Peebles 1983) to derive the large-scale linear galaxy bias and the associated typical DMH mass. To make full use of the clustering measurements, the smaller separations of the nonlinear regime (R 1 − 2 h −1 cMpc) are modeled by relating galaxies to DMHs within the nonlinear framework of halo occupation distribution (HOD) modeling. In this context, the mean number of galaxies in the DMH is modeled as a function of DMH mass, further assessing whether these galaxies occupy the centers of the DMHs or whether they are satellite galaxies.

Data
The MUSE spectroscopic surveys are based on a wedding cake design, namely: a first spatially wide region (bottom of the cake) is observed with a short exposure time (1 hour), while deeper observations (10 hours exposure) are carried out within the first surveyed area (middle tier of the cake). Contained in the last observed region, an even deeper survey (140 h) is then built up (at the top of the cake). These three surveys are known as: MUSE-Wide (Herenz et al. 2017;Urrutia et al. 2019), MUSE-Deep (Bacon et al. 2017;Inami et al. 2017;Bacon et al. 2022), and MUSE Extremely Deep Field (MXDF; Bacon et al. 2022). Each of them can be seen as a different layer of a wedding cake, where higher layers become spatially smaller and correspond to deeper observations. In what follows, we give further details on survey and galaxy sample construction.

MUSE-Wide
The spectroscopic MUSE-Wide survey (Herenz et al. 2017;Urrutia et al. 2019) comprises 100 MUSE fields distributed in the CANDELS/GOODS-S, CANDELS/COSMOS and the Hubble Ultra Deep Field (HUDF) parallel field regions. Each MUSE field covers 1 arcmin 2 . While 91 fields were observed with an exposure time of one hour, nine correspond to shallow (1.6 hours) reduced subsets of the MUSE-Deep data (see next section; as well as Bacon et al. 2017), located within the HUDF in the CANDELS/GOODS-S region. However, we do not include the objects from this region since they overlap with the MUSE-Deep sample (see next section and gap in the left panel of Fig. 1). The slight overlap between adjacent fields leads to a total spatial coverage of 83.52 arcmin 2 . The red circles in Fig. 1 display the spatial distribution of the LAEs from the MUSE-Wide survey. We refer to Urrutia et al. (2019) for further details on the survey build up, reduction and flux calibration of the MUSE data cubes.
In this paper, we extend (x2 spatially, 50% more LAEs) the sample used in Herrero Alonso et al. (2021) and include all the 1 h exposure fields from the MUSE-Wide survey. Despite the somewhat worse seeing (generally) in the COSMOS region (right panel of Fig.1), we demonstrate in Appendix A that adding these fields does not significantly impact our clustering results but helps in minimizing the effects of cosmic sample variance.
We also expanded the redshift range of the sample. While MUSE spectra cover 4750-9350 Å, implying a Lyα redshift interval of 2.9 < ∼ z < ∼ 6.7, we limited the redshift range to 3 < z < 6 (differing from the more conservative range of Herrero Alonso et al. 2021; 3.3 < z < 6) as the details of the selection function near the extremes are still being investigated. Section 2 of Herrero Alonso et al. (2021) describes the aspects relevant to our analysis on the construction of a sample of LAEs, as well as the strategy to measure line fluxes and redshifts. The redshift distribution of the sample is shown in red in the top panel of Fig. 2. Systematic uncertainties introduced in the redshift-derived 3D positions of the LAEs have negligible consequences for our clustering approach (see Sect. 2.2 in Herrero Alonso et al. 2021).
Within 83.52 arcmin 2 and in the selected redshift interval, we detected a total of 1030 LAEs. This implies a LAE density of more than 13 objects per arcmin 2 or n ≈ 1 · 10 −3 h 3 Mpc −3 (for 3 < z < 6). At the median redshift of the sample z = 4.0, the transverse extent of the footprint is ≈ 43 h −1 Mpc. The range of Lyα luminosities is 40.92 < log(L Lyα /[erg s −1 ]) < 43.35 (see red circles in Fig. 2), with a median value of log L Lyα /[erg s −1 ] = 42.34 (or ≈ L in terms of characteristic luminosity L ; Herenz et al. 2019), which makes this sample the highest luminosity data set of our three considered surveys. The Lyα luminosity distribution is shown in red in the right panel of Fig. 2. The main properties of the MUSE-Wide LAEs are summarized in Table 1.   spatial coverage is 9.92 arcmin 2 . We represent the spatial distribution of the survey in green in Fig. 1. We did, however, remove the MUSE-Deep objects that are selected in the deepest survey, described in the next section. We refer to Bacon et al. (2017Bacon et al. ( , 2022 for a detailed description on survey construction and data reduction. The sources in MUSE-Deep were blindly detected and extracted using ORIGIN (Mary et al. 2020), based on a matched filtering approach and developed to detect faint emission lines in MUSE datacubes. While the redshift measurements and line classifications were carried out with pyMarZ, a python version of the redshift fitting software MarZ (Hinton et al. 2016), the line flux determination was conducted with pyPlatefit, which is a python module optimized to fit emission lines of high-redshift spectra. The redshift distribution of the sample is shown in green in the top panel of Fig. 2, also within 3 < z < 6.

MUSE Extremely Deep
The MUSE Extremely Deep Field (Bacon et al. 2022) is situated in the CANDELS/GOODS-S region and overlaps with MUSE-Deep and MUSE-Wide. It is composed of a single quasi circular field with inner and outer radii of 31" and 41", respectively. While a 140 hour exposure was employed to observe the totality of the field, the inner field is 135 hours deep, decreasing to 10 hours depth at the outer radius. This makes MXDF the deepest spectroscopic survey to date. For further details see Bacon et al. (2022) and the blue data points in Fig.1, where the MXDF field is overplotted on the previous surveys.
The survey assembly and data reduction is described in Bacon et al. (2022) and is similar to the one applied to MUSE-Deep . The source extraction in MXDF and the redshift and flux measurements are conducted following the same procedure as was done for MUSE-Deep. The redshift distribution of the sample is shown in blue in the top panel of Fig. 2.
Contained within ≈1.47 arcmin 2 and over the same redshift range as for the previous catalogues, we detected 367 LAEs, corresponding to a LAE density of n ≈ 3 · 10 −2 h 3 Mpc −3 (432 LAEs per arcmin 2 at 3 < z < 6). With a median redshift of z = 4.2, the footprint covers ≈ 2.8 h −1 Mpc (transversely). The Lyα luminosities span 40.15 < log(L Lyα /[erg s −1 ]) < 43.09 (see blue stars in Fig. 2 and its distribution in the right panel). The median Lyα luminosity is log L Lyα /[erg s −1 ] = 41.22 (or ≈ 0.08L ), more than one order of magnitude fainter than for MUSE-Wide. This makes MXDF the faintest ever observed sample of non-lensed LAEs. The main properties are listed in Table 1.

LAE subsamples
We bisected the main samples into disjoint subsets based on their median Lyα luminosity to investigate the clustering dependence on this quantity. We did not merge the main LAE datasets because their distinct Lyα luminosities, together with their slightly different location on the sky, might introduce systematics in the clustering measurements. The subsample properties are summarized in Table 2 and described in the following.
We split the MUSE-Wide sample at the median Lyα luminosity log L Lyα /[erg s −1 ] = 42.34. The two subsamples consist of 515 LAEs each. The low-luminosity subset has a median redshift and Lyα luminosity of z low = 3.7 and log L Lyαlow /[erg s −1 ] = 42.06, while the high-luminosity subsample has z high = 4.1 and log L Lyαhigh /[erg s −1 ] = 42.53. The median redshift of the number of galaxy pairs for the lowluminosity subset is z pair ≈ 3.4, and that for the high-luminosity one is z pair ≈ 4.1.
We next bisected the MUSE-Deep set at the median Lyα luminosity log L Lyα /[erg s −1 ] = 41.64. The low-luminosity subsample has 340 LAEs and presents a median redshift and Lyα luminosity of z low = 3.7 and log L Lyαlow /[erg s −1 ] = 41.46. The high-luminosity subset is formed by 339 LAEs with z high = 4.5 and log L Lyαhigh /[erg s −1 ] = 41.89. While for the lowluminosity subsample z pair ≈ 3.5, for the high-luminosity one z pair ≈ 4.4.
We also divide the sample with the largest dynamic range of Lyα luminosities (MXDF) at the median Lyα luminosity log L Lyα /[erg s −1 ] = 41.22. While the lower luminosity subset contains 183 LAEs with z low = 4.0 and log L Lyαlow /[erg s −1 ] = 40.97, the higher luminosity subsample consists of 184 LAEs with z high = 4.5 and log L Lyαhigh /[erg s −1 ] = 41.54. For the low-luminosity subset, we have z pair ≈ 3.9, and for the high-luminosity one, we have z pair ≈ 4.8.
The redshift distribution of each subsample is shown in Fig. 3. The corresponding median redshifts are represented with a vertical dashed line. Despite the similar median redshifts between the subsample pairs, the redshift distributions are significantly different, with a higher amount of spike-trough contrasts in the high-luminosity subsets.

K-estimator
Galaxy clustering is commonly measured by two-point correlation function (2pcf) statistics. Samples investigated by this method typically span several square degrees on the sky. With MUSE, we encounter the opposite scenario. By design, MUSE surveys cover small spatial extensions on the sky and provide a broad redshift range. Although the MUSE-Wide survey is the largest footprint of all MUSE samples, its nature is still that of a pencil-beam survey. Its transverse scales are of the order of 40 h −1 Mpc, while in redshift space it reaches almost 1500 h −1 Mpc. If we consider the deeper samples, the difference is even more prominent: 8.7 vs 1500 h −1 Mpc for MUSE-Deep and 2.8 versus 1500 h −1 Mpc for MXDF. It is thus paramount to exploit the radial scales and utilize alternative methods to the traditional 2pcf.
In Herrero Alonso et al. (2021) we applied the so-called Kestimator, introduced by Adelberger et al. (2005), to a subset of our current sample. Here, we build on our previous work by extending the dataset and measuring the small-scale clustering required to perform full HOD modeling. The details of the Kestimator are given in Sect. 3.1 of Herrero Alonso et al. (2021). In the following, we provide a brief description of the method.
The K-estimator measures the radial clustering along line-ofsight distances, Z i j , by counting galaxy pairs (formed by galaxy i and galaxy j) in redshift space at fixed transverse separations, R i j . Although the K-estimator does not need a random sample to carry out the clustering measurements, its nature is very similar to that of the projected two-point correlation function. We bin by R i j , shown with distinct radii in the cylinders of Fig. 4, and count the number of pairs within individual transverse bins, for two different ranges of Z i j , represented in red and blue in Fig. 4. The K-estimator as a function of R i j is then defined as the ratio of galaxy pairs within the first Z i j interval (blue cylinder) and the total Z i j range (red and blue cylinder), quantifying the excess of galaxy pairs in the first Z i j bin with respect to the total one. We optimize the choice of the Z i j ranges, and thus the K-estimator, by seeking out the estimator that delivers the best sensitivity for the clustering signal (i.e., the highest signal-tonoise ratio, S/N; see Sect. 3.1.2 in Herrero Alonso et al. 2021). Although slightly different than in Herrero Alonso et al. (2021), we find nearly identical K-estimators for each of the current samples (K 0,7 7,45 for MUSE-Wide, K 0,7 7,45 for MUSE-Deep, and K 0,7 7,40 for MXDF), whose clustering signals only differ in their S/N. We chose the same K-estimator for the three data sets, K 0,7 7,45 . The K-estimator is directly related to the average underlying correlation function (see Eq. 2 in Herrero Alonso et al. 2021). In fact, its definition is proportional to a combination of projected two-point correlation functions corresponding to the blue and red cylinders of Fig. 4. While the traditional 2pcf method integrates the correlation function ξ(R i j , Z i j ) over line-of-sight separations up to a maximum line-of-sight distance π max , the Kestimator integrates up to a 2 and a 3 . The correlation function ξ(R i j , Z i j ) can be approximated with a power-law following Limber (1953) equations as we did in Herrero Alonso et al. (2021), or modeled with a halo occupation distribution (HOD) model (see Sect. 3.3). For reference, randomly distributed galaxies in space (ξ(R i j , Z i j ) = 0) provide K 0,7 7,45 (R i j ) values equal to 7/45 (see Eq. 2 in Herrero Alonso et al. 2021). Samples with data points significantly above 7/45 dispense clustering signals.

Error estimation for the MUSE-Wide survey
Applying clustering statistics delivers correlated data points. One single galaxy might be part of more than one galaxy pair and can therefore contribute to several R i j bins, especially if they are adjacent. In order to quantify the actual correlation between data points, we applied the jackknife resampling technique, followed by the computation of the covariance matrix (see e.g., Krumpe et al. 2010;Miyaji et al. 2011). For the MUSE-Wide sample, we employed ten logarithmic bins in the range 0.16 < R i j < 27.5 h −1 Mpc, discarding lower R i j scales since they host very few galaxy pairs.
We then found a compromise between the number of independent regions (jackknife zones) and the size of the jackknife zones and divide the sky coverage into N jack = 10 regions, each of which extends ≈ 4 h −1 Mpc in both RA and Dec directions (see Appendix B for a visual representation of the sky division). The limited spatial extent of the survey does not allow for a higher number of jackknife zones. We then constructed N jack jackknife subsamples, excluding one jackknife zone at a time, and computed the K-estimator for each of the subsets. The Kestimator measurements are then used to derive the covariance Fig. 4: Sketch of the K-estimator, representing the relative geometry that probe the one-and two-halo term scales. The empty blue and filled red cylinders, delimited by |a 2 | = 7 h −1 Mpc and |a 3 | = 45 h −1 Mpc respectively, illustrate the line-of-sight distance Z i j intervals within which we count galaxy pairs at fixed transverse separations R i j , represented by nested cylinders. Pairs of LAEs connected with green lines within the same DMH (filled gray circle) contribute to the one-halo term (small R i j scales), while pairs belonging to two different DMHs (yellow lines) probe the two-halo term (larger R i j separations). matrix M i j , which quantifies the correlation between bins i and j. The matrix is expressed as where K k (R i ), K k (R j ) are the K-estimators from the k-th jackknife samples and K(R i ) , K(R j ) are the averages over all jackknife samples in the i, j bins, respectively. The error bar for the K-estimator at the ith bin comes from the square root of the diagonal element ( √ M ii ) of the covariance matrix, our so-called "jackknife uncertainty." This approach could not be followed in Herrero Alonso et al. (2021) because of the smaller sky coverage. Instead, we used a galaxy bootstrapping approach. In Appendix C, we compare the two techniques and show that bootstrapping uncertainties are ≈ 50% larger than the jackknife error bars, in agreement with Norberg et al. (2009), who found that boostrapping overestimates the uncertainties. We next search for the best-fit parameters by minimizing the correlated χ 2 values according to where N bins = 10 is the number of R i j bins, K(R i ), K(R j ) are the measured K-estimators and K(R i ) HOD , K(R j ) HOD are the Kestimators predicted by the HOD model for each i, j bin, respectively.
Regardless of the larger sample considered in this work, we are still limited by the spatial size of the survey, which only permits a small number of jackknife zones. The insufficient statistics naturally lead to a higher noise contribution in the covariance matrix, which cause the χ 2 minimization to mathematically fail (i.e., cases of χ 2 < 0) when the full covariance matrix is included. Hence, we only incorporated the main diagonal of the matrix and its two contiguous diagonals. In Appendix B, we discuss the high level of noise in the matrix elements corresponding to bins that are significantly apart from each other. We also verify the robustness of our approach and show that our clustering results are not altered (within 1σ) by this choice.

Error estimation for the deeper surveys
The small sky coverage of the deeper surveys does not allow us to follow the same error estimation approach as for the MUSE-Wide survey. In Appendix C, we not only compare the bootstrapping technique applied in Herrero Alonso et al. (2021) to the jackknife approach performed in MUSE-Wide, but we also consider the Poisson uncertainties. We demonstrate that Poisson and jackknife errors are comparable in our sample. In fact, we show that while bootstrapping uncertainties are ≈ 50% larger than jackknife errors, Poisson uncertainties are only ≈ 7% higher. Thus, and similarly to Adelberger et al. (2005); Diener et al. (2017); Khostovan et al. (2018), we stick to Poisson uncertainties for the MUSE-Deep and MXDF samples. For these datasets, we measure the K-estimator in eight and six logarithmic bins in the ranges 0.09 < R i j /[h −1 Mpc] < 4.75 and 0.09 < R i j /[h −1 Mpc] < 1.45, respectively, constrained by the spatial extent of the surveys.
We then perform a standard χ 2 minimization to find the best fitting parameters to the K-estimator measurements. Namely, where K(R i ), K(R i ) HOD , and σ i denote the measured Kestimator, the HOD modeled K-estimator and the Poisson uncertainty in the ith bin, respectively. We note that the standard χ 2 minimization does not account for the correlation between bins. Although in Appendix B we show that only contiguous bins are moderately correlated, we should take the resulting fit uncertainties with caution.

Halo occupation distribution modeling
The clustering statistics can be approximated with a powerlaw or modeled with state-of-the-art HOD modeling. Traditional clustering studies make use of power laws to derive the correlation length and slope, from which they infer large-scale bias factors and typical DMH masses. This simple approach deviates from the actual shape of the clustering statistic curve, even in the linear regime, and its inferred DMH masses suffer from systematic errors (e.g., Jenkins et al. 1998 and references therein). To overcome these concerns, physically motivated HOD models do not treat the linear and non-linear regime alike but differentiate between the clustering contribution from galaxy pairs that reside in the same DMH and pairs that occupy different DMHs.
In Herrero Alonso et al. (2021) we only modeled the twohalo term of the K-estimator with HOD modeling, which only delivered the large-scale bias factor and the typical DMH mass of the sample. We now extend into the non-linear regime (i.e., R i j < 0.6 h −1 Mpc) of the one-halo term. We can then model the clustering measured by the K-estimator with a full HOD model, combining the separate contributions from the one-(1h, i.e., galaxy pairs residing in the same DMH) and the two-halo (2h, i.e., galaxy pairs residing in different DMHs) terms: where ξ is the correlation function.
The HOD model we used is the same as in Herrero Alonso et al. (2021), an improved version of that described by Miyaji et al. (2011); Krumpe et al. (2012Krumpe et al. ( , 2015Krumpe et al. ( , 2018. We assumed that LAEs are associated with DMHs, linked by the bias-halo mass relation from Tinker et al. (2005). From Tinker et al. (2005), we also included the effects of halo-halo collisions and scaledependent bias. The mass function of DMHs, which is denoted by φ(M h )dM h , is based on Sheth et al. (2001), and the DMH profile is taken from Navarro, Frenk & White (1997). We use the concentration parameter from Zheng et al. (2007), and the weakly redshift-dependent collapse overdensity from Navarro, Frenk & White (1997); van den Bosch et al. (2013). We further incorporated redshift space distortions (RSDs) in the twohalo term using linear theory (Kaiser infall;Kaiser 1987 andvan den Bosch et al. 2013). We did not model RSDs in the onehalo term because the peculiar velocity has negligible effects to our K-estimator as demonstrated in the following. The velocity dispersion (σ v ) of satellites in a M h halo can be estimated by σ 2 v ≈ GM h /(2R vir ), where R vir is the virial radius (Tinker 2007). Its effect on the line-of-sight physical distance estimate is then σ v /H(z). For 10 11−12 h −1 M DMH masses, which are typical for our sample, with virial radii of ≈ 0.02 − 0.05 (physical) h −1 Mpc, the line-of-sight distance estimation is deviated by ≈ 0.15 − 0.30 h −1 Mpc, corresponding to a peculiar velocity dispersion of σ v ≈ 80 − 170 km s −1 . This is significantly small compared to our a 2 = 7 h −1 Mpc. We thus assume that the onehalo term contributes only to the Z i j = 0 − 7 h −1 Mpc bin. We evaluated the HOD model at the median redshift of N(z) 2 , where N(z) is the redshift distribution of the sampled galaxy pairs. For our three main datasets, z pair ≈ 3.8.
The mean halo occupation function is a simplified version of the five parameter model by Zheng et al. (2007). We fixed the halo mass at which the satellite occupation becomes zero to M 0 = 0 and the smoothing scale of the central halo occupation lower mass cutoff to σ log M = 0, due to sample size limitations. We define the mean occupation distribution of the central galaxy and that of satellite galaxies N s (M h ) as where M min is the minimum halo mass required to host a central galaxy, M 1 is the halo mass threshold to host (on average) one satellite galaxy, and α is the high-mass power-law slope of the satellite galaxy mean occupation function. The total halo occupation is given by the sum of central and satellite galaxy halo The dependencies of the HOD parameters on the shape of the K-estimator are detailed in Appendix D. In short, for the HOD parameters there selected, the clustering amplitude of the twohalo term is ascertained by the hosting DMHs and is thus very sensitive to their mass, M min , and to the fraction of galaxies in massive halos with respect to lower-mass halos, linked to α. The clustering in the one-halo term regime, however, is affected by the three parameters in a complex manner; roughly M min and α vary the amplitude, and α as well as (moderately) M 1 modify the slope.
To find the best-fit HOD model, we construct a 3D parameter grid for M min , M 1 , and α. We vary log(M min /[h −1 M ]) in the range 9.5 − 11.2, log(M 1 /M min ) from 0.5 to 2.5, and α within 0.2 − 4.3, all in steps of 0.1. For each parameter combination, we computed ξ (Eq. 4), converted it to the K-estimator using Eq. 2 in Herrero Alonso et al. (2021), and computed a χ 2 value (Eqs. 2 or 3). We then used the resulting 3D χ 2 grid to estimate the confidence intervals for the HOD parameters. For each point on a 2D plane, we search for the minimum χ 2 for the contouring along the remaining parameter. The contours we plot are at ∆χ 2 = 3.53 and 8.02, which correspond to Gaussian 68% (1σ) and 95% (2σ) confidence levels, respectively, applying the χ 2 distribution for three degrees of freedom. The projections of the 68% probability contours on the three interesting parameters are then used to compute the uncertainty of each HOD parameter.
For each point in the three parameter grid, we also computed the large-scale galaxy bias factor, b, and the fraction of satellite galaxies per halo, f sat , as follows: where b h (M h ) denotes the large-scale halo bias. The typical DMH mass is determined by the large-scale galaxy bias factor. We ultimately compute the bias and f sat distributions from the HOD models that fall within the 68% confidence (for the threeparameter space) contours. These distributions are then used to assess the uncertainties in the bias and f sat .

Fit results from the MUSE-Wide survey
Using the K-estimator K 0,7 7,45 , we compute the clustering of our LAE sample in ten logarithmic bins in the range 0.16 < R i j /[h −1 Mpc] < 27.5, with error bars calculated following the jackknife resampling technique described in Sect. 3.2.1. In the left top panel of Figure 5, we show the measured clustering signal, with all MUSE-Wide data points significantly above the 7/45 baseline, which represents the expected clustering of an unclustered population.
Following the procedure laid out in Sect. 3.3, we obtain constraints on the HOD parameters. From the grid search and the χ 2 minimization, we find the best HOD fit to the K-estimator, colored in black in the same figure and dissected into the oneand two-halo term contributions. It can be seen from the residuals (bottom) that the model is in remarkable agreement with the measurements.
A somewhat intriguing feature, at least at first sight, is the kink in the two-halo term profile at 0.2 < R i j /[h −1 Mpc] < 0.4. This reflects the effect of the halo-halo collision introduced in the HOD model formalism by Tinker et al. (2005), where the galaxy pairs within the same DMH cannot contribute to the two-halo term.
Our fitting allows us to find the best-fit HOD from Eqs. 5 and 6. In the right top panel of Fig. 5, we represent the best HODs for the central, satellite, and total LAEs from the MUSE-Wide survey. While the halo mass needed to host one (central) LAE is log(M h /[h −1 M ]) > 10.6, satellite galaxies are only present if the DMHs are at least one order of magnitude more massive (log(M h /[h −1 M ]) > 11.6).
As described in Sect. 3.3, we also compute the confidence regions for the HOD parameters. We show the probability contours (red) in Fig. 6. The wobbliness of the curves, especially those involving α, is caused by making use of a discrete grid. We list the best-fit HOD parameters in Table 3 Seeking robust information about the number of satellite galaxies, we compute the satellite fraction f sat (Eq. 8) for each parameter combination. We find f sat 0.10 at the 3σ confidence level, being f sat = 0.012 +0.018 −0.009 . That is, ≈ 3% (1σ upper limit) of the LAEs in the MUSE-Wide survey are satellites. In other words, at most ≈ 2 out of ≈ 65 DMHs in our sample host one satellite LAE.

Fit results from MUSE-Deep
We measure the clustering of the MUSE-Deep LAE sample with the same K-estimator in eight logarithmic bins within 0.09 < R i j /[h −1 Mpc] < 4.75. We compute Poisson uncertainties as laid out in Sect. 3.2.2 and display the result in the middle left panel of Fig. 5. Overplotted on the clustering signal, we show the best HOD fit, split into the one-and two-halo term contributions. The good quality of the fit is quantified with the residuals in the bottom panel of the figure.
Following the procedure described in Sect. 3.2.2, we compute the confidence intervals for the HOD parameters and list them in Table 3. We plot the probability contours (green) in satellite fraction is f sat = 0.004 +0.009 −0.002 , consistent with that from the MUSE-Wide LAE sample.
We then compute the best-fit HOD for central, satellite and total LAEs (middle right panel of Fig. 5). In line with the bestfit HOD parameters and somewhat lower than the values found for the MUSE-Wide survey, the smallest DMH that can host a central LAE has a mass of log(M h /[h −1 M ]) > 10.4, more than one order of magnitude lower than that required to host one additional LAE (satellite).

Fit results from the MUSE Extremely Deep Field
We make use of six logarithmic bins in the range 0.09 <  Notes: z is the median redshift of the sample. M min , M 1 are the threshold DMH masses to host a central and a satellite LAE, respectively. α is the high-mass power-law slope of the number of satellite galaxies, f sat is the satellite fraction, b is the large-scale bias factor and M h is the typical DMH mass of the galaxy sample. It is worth pointing out that the three HOD parameters have some degree of degeneracy, printed out in the diagonally elongated probability contours in log M 1 /M min -α space in the bottom right panel of Fig. 6. This can be understood as follows: a higher α in the models causes an increase of satellites at high mass halos, but this can be compensated by producing less satellites by increasing log M 1 /M min . While this correlation is clearly visible for the MUSE-Wide and MUSE-Deep samples, the MXDF dataset only seems to be affected in the 95% confidence contour. We did not observe clear correlations between other parameters with any of our samples. Appendix A shows how our K-estimator varies with the parameters. The causes of parameter degeneracies are also noticeable in Fig. D.1. We note however that while the correlation between the HOD parameters leads to the perturbed shape of the probability contours, the lowest (MXDF) and highest luminosity (MUSE-Wide) sample contours are detached from each other. Thus, for the purposes of this study, simultaneously fitting the three HOD parameters and showing their correlations is preferable over, for instance, fixing α to a dubious value.

Clustering dependence on Lyα luminosity
The complex radiative transfer processes that the Lyα photons are subject to make the search for correlations between Lyα luminosity and other physical properties a difficult task. Despite this complication, Yajima et al. (2018) predicted a correlation between simulated L Lyα and halo mass based on halo merger trees and Lyα radiative transfer calculations. Khostovan et al. (2019) is, however, the only study so far that has reported a clear (5σ) relation between these quantities using observational data. Motivated by these results, we exploited the large dynamic range of Lyα luminosities that we cover to investigate the relation between Lyα luminosity and DMH mass. As a first step, we compare the K-estimator measurements in the MUSE-Wide survey (highest luminosity LAE sample: L Lyα ≈ 10 42.34 erg s −1 , but still fainter than those in Khostovan et al. (2019)) and in MXDF (faintest LAE sample; L Lyα ≈ 10 41.22 erg s −1 ) and show the outcome of this comparison in the left panel of Fig. 7.
The relatively luminous LAEs from the MUSE-Wide survey cluster slightly more strongly (b Wide = 2.65 +0.13 −0.11 ) than the lowluminosity LAEs from MXDF (b MXDF = 2.43 +0.15 −0.15 ). The clustering measurements and bias factor (b = 2.42 +0.10 −0.09 ) in MUSE-Deep (log(L Lyα /[erg s −1 ]) = 41.64) fall between those from MUSE-Wide and MXDF. We convert the bias factors from the three main samples of this study into typical DMH masses and plot them as a function of their median Lyα luminosity with colored symbols in Fig. 8.
Although the three main datasets sample the same region of the sky, their transverse coverage is limited and somewhat differs. Therefore, our results are affected by cosmic sample variance. Ideally, this uncertainty is estimated from the variance of clustering measurements from simulated mocks in different lines of sight. Inferring cosmic variance from a large set of mocks that are able to reproduce the observed clustering of our LAEs is however beyond the scope of this paper.
We further investigate the possible dependence on L Lyα by splitting the main LAE samples into disjoint subsets (see Table 2). We compute the K-estimator in each L Lyα subsample, find the best HOD fit and list the large-scale bias factors and the typical DMH masses in Table 4. We also plot the typical DMH masses in Fig. 8 (empty symbols) as a function of the median L Lyα of the subsamples. We find that typical halo mass increases from 10 10.00 to 10 11.43 M between 10 40.97 and 10 42.53 erg s −1 in line luminosity. For each subsample pair, the high-luminosity subset always clusters more strongly than the low-luminosity one and, in this case, cosmic sample variance effects can be completely neglected because subset pairs span the exact same area on the sky. The most pronounced difference is found when splitting the MXDF sample, the dataset with the largest dynamic range of Lyα luminosity. The best HOD fits deliver b low = 1.79 +0.08 −0.06 and b high = 3.10 +0.24 −0.22 (3.9σ significant). Despite its higher luminosity, we infer a less massive DMH for the MUSE-Deep high-luminosity subsample than for the main dataset. This is due to the higher z pair of the subset (see Sect. 2.4 and 3.3). Because we evaluate the HOD model at z pair , a higher redshift corresponds to HOD models in which the halo mass function presents a lower number density of massive halos and, thus, deliver less massive typical DMHs. The same reasoning applies when comparing the high-luminosity MXDF and low-luminosity MUSE-Deep subsamples and the high-luminosity MUSE-Deep and low-luminosity MUSE-Wide subsets. While each subsample pair presents similar median luminosities, the former also has similar z pair , unlike the latter one (see Sect. 2.4). This translates into similar DMH masses for the first pair but significantly distinct masses for the second.
We last consider the most extreme cases, the low-luminosity subset from MXDF and the high-luminosity one from the MUSE-Wide survey. We show the measured clustering in the two subsamples in the right panel of Fig. 7. The highluminosity LAEs cluster 8σ more strongly than the lowluminosity LAEs, without accounting for cosmic variance. We −0.09 . These results fit well within the assumed framework in which star-forming galaxies that reside in more massive halos present higher star formation rates and thus show more luminous nebular emission lines (Kusakabe et al. 2018). This dependence can then be weakened by low Lyα escape fractions in high mass halos.
Following Sect. 5.4.1 of Herrero Alonso et al. (2021), we matched the redshift distributions of the three main samples and of each subsample pair to verify that the difference in clustering amplitude is not driven by the different redshift distribution of the datasets. For each main sample, we compare individual bins between their corresponding z-distributions and select the one that contains a higher number of objects. We then randomly remove LAEs until we match the number counts of the non-selected samples in that bin. Once all bins have been inspected, we obtain "matched" z-distributions (i.e., equivalent), but with still different Lyα luminosity distributions. We ran the K-estimator in the three "matched" datasets and find consistent results with the original ones. We follow the same approach for the subsamples such that the low-and high-luminosity subsets have exactly the same z-distribution. We find that the clustering difference between the "matched" and original subsamples varies within 1σ. Besides, as we did for L Lyα , we also searched for a possible clustering dependence on redshift and found no trend. Thus, we discarded the possibility of a possible clustering dependence on Lyα luminosity driven by z.
Our results are not driven by AGN or low-redshift emission line contamination either. The Lyα-emitting AGN fraction for L Lyα < 10 43 erg s −1 is close to zero (Spinoso et al. 2020 and references therein) and the four known X-ray detected AGNs (Luo et al. 2017), which only affect MUSE-Wide and MUSE-Deep, were not included in our datasets. Besides, Urrutia et al. (2019) performed a stacking experiment of X-ray images centered on MUSE-Wide LAEs, yielding no signal. The presence of lowredshift interlopers in our spectroscopic samples is also unlikely. More significant is the dependence found in Khostovan et al. (2019) andHerrero Alonso et al. (2021). While the latter measured a 2σ difference in bias factors or DMH masses between two subsets of 349 and 346 LAEs at z ≈ 4 with log(L Lyα /[erg s −1 ]) ≈ 42.14 and log(L Lyα /[erg s −1 ]) ≈ 42.57, the former used various surveys with discrete redshift slices between 2.5 < z < 6 and 42.0 < log(L Lyα /[erg s −1 ]) < 43.6 to find that halo mass clearly (5σ) increases with increasing line luminosity. For a direct comparison, we plot in Fig. 8 (gray crosses) the DMH masses computed by Khostovan et al. (2019) from samples with similar redshifts (z ≈ 3) and Lyα luminosities (log(L Lyα /[erg s −1 ]) ≈ 42) to our current LAE samples. Our results are in good agreement and extend to much fainter Lyα luminosities.
Our results, along with those from the literature, demonstrate that having a broad dynamic range of L Lyα (nearly extending two orders of magnitude) and a large number of LAEs in the samples is crucial to detect the clustering dependence on L Lyα . Alonso et al. (2021) In this section we compare our results with the findings of our previous study (Herrero Alonso et al. 2021, hereafter HA21), where we measured the clustering of a subset (68 fields of the MUSE-Wide survey) of our current sample (91 fields of the MUSE-Wide survey) and fitted the corresponding signal with a two-halo term only HOD modeling. In order to envisage the methodological and statistical improvement of our new investigation, we applied our K 0,7 7,45 estimator to the sample considered in HA21 (695 LAEs at 3.3 < z < 6). We compare the outcome to our current clustering measurement in Fig. 9.

Comparison to Herrero
The two datasets show good agreement within the uncertainties, with smaller errors for the current sample. Besides the higher number of LAEs and larger spatial coverage, the error estimation was carried out following different procedures. While the spatial coverage of the full MUSE-Wide survey allows us to compute the covariance matrix from the jackknife resampling technique, the smaller transverse extent covered by the 68 fields did not allow the split of the surveyed area into a significant number of jackknife zones. Thus, in HA21, we chose bootstrapping error bars as our next most conservative and realistic approach.
The slightly puzzling hump seen in Sect. 4 of HA21 at 4 R i j /[h −1 Mpc] 7 is no longer visible in our new dataset. This confirms the judgement in HA21 that the feature was consistent with a statistical fluctuation resulting from the correlation between datapoints.
In HA21, we limited the range of transverse separations to R i j > 0.6 h −1 Mpc, excluding the smallest scales of the one-halo term. Thus, we fitted the signal with a two-halo term only HOD model (red dotted curve in Fig. 9) in contrast to the full HOD modeling performed in this work (blue dotted curve). While the former only constrained the large-scale bias factor and the typical DMH mass of LAEs, the latter further determines the number of central and satellite galaxies, as well as the required DMH mass to host each type of galaxy. Despite these dissimilarities, the two fits are in good agreement: the bias factor (b = 2.80 +0.38 −0.38 ) and the typical DMH mass of LAEs (log(M DMH / [h −1 M ]) = 11.34 +0.23 −0.27 ) from HA21 are consistent with those derived in this work (b = 2.65 +0.13 −0.11 and log(M DMH / [h −1 M ]) = 11.09 +0.10 −0.09 ). The higher accuracy of our current measurements originates from the larger sample, the availability of more realistic error bars, and constraints from the one-halo term.

Comparison to the literature
A common way to infer the host DMH masses of LAEs is to quantify the galaxy clustering of the detected population through clustering statistics, which is then traditionally approximated with power-laws or fit with physically motivated HOD models.
Following the traditional approach, Gawiser et al. (2007), Ouchi et al. (2010) and Bielby et al. (2016) focused on the clustering of a few hundred LAEs at z = 3.1 − 6.6 to obtain typical DMH masses in the range 10 10 − 10 11 M . Similar masses were found by Khostovan et al. (2018) in a much larger sample (≈ 5000 LAEs) in discrete redshift slices within 2.5 < z < 6, adopting the same procedure. A major improvement in terms of methodology was presented in Lee et al. (2006); Durkalec et al. (2014); Ouchi et al. (2018); Durkalec et al. (2018), who considered samples of high-z galaxies (2000-3000 mainly LAEs and Lyman-break galaxies, LBGs) and quantified the clustering with HOD modeling. While Ouchi et al. (2018) found that their LAEs at z = 5.7 (6.6) are hosted by DMHs with typical masses Fig. 9: Clustering of the full MUSE-Wide sample (blue; this work) compared to the subset considered in HA21 (red). The former measurements show jackknife uncertainties (see Sect. 3.2.1) and the latter bootstrapping errors (see Sect. 3.1.3 in HA21). The blue dotted curve represents our best-fit from full HOD modeling. The red dotted curve displays the two-halo term only best HOD fit found in Sect. 4.3 of HA21. The black straight line shows the expected K value of an unclustered sample. Lee et al. (2006) and Durkalec et al. (2014Durkalec et al. ( , 2018 computed log(M h /h −1 M ) ≈ 11.7 for their sample of galaxies at z = 4 − 5 and z = 3, respectively. Considering that we have performed a full HOD modeling at the median redshift of our number of galaxy pairs (z pair = 3.8) and that the DMH masses are predicted to evolve with cosmic time, our derived typical DMH masses log(M h /h −1 M ) ≈ 10.77 − 11.09 are in good agreement with the literature.
Besides the computation of typical DMH masses, modeling the one-halo term of the clustering statistics with HOD models delivers the minimum DMH mass required to host a central galaxy, M min , that is needed for a satellite galaxy, M 1 , and the power-law slope of number of satellites, α. These three parameters constrain the satellite fraction, f sat . Ouchi et al. (2018) partially exploited the power of HOD models in a sample of ≈ 2000 LAEs to obtain log(M min /M ) = 9.5 +0.5 −1.2 (9.1 +0.7 −1.9 ) at z = 5.7 (6.6). Our derived minimum masses to host a central galaxy at z pair = 3.8 are considerably larger (log(M min /M ) ≈ 10.3−10.7), which can be explained by the different Lyα luminosities covered in the two studies, and by the fact that several HOD parameters were fixed in Ouchi et al. (2018), namely, σ log M = 0.2, log M 0 = 0.76M 1 +2.3, log M 1 = 1.18 log M min −1.28, and α = 1, which are not compatible with ours. This was the only previous study that performed HOD modeling in a sample of LAEs. Lee et al. (2006) and Durkalec et al. (2014) made use of the full potential of HOD models to reproduce the clustering of their LBG population at z = 4 − 5 and 2.9 < z < 5, respectively. Although it is still under debate whether LBGs and LAEs are the same galaxy population (Garel et al. 2015 and references therein), Lee et al. (2006) computed a minimum DMH mass to host a central LBG of log(M min /M ) ≈ 10.8, to host a satellite LBG of log(M 1 /M ) ≈ 12.0, and a power-law slope α for the number of satellites of α ≈ 0.7, with considerable uncertainties. Similarly, Durkalec et al. (2014) found log(M min /M ) = 11.18 +0.56 −0.70 , log(M 1 /M ) = 12.55 +0.85 −0.88 , and α = 0.73 +0.23 −0.30 . While their halo masses are in agreement with our findings, their slope is somewhat shallower. This is partially expected given the dis-similarities in the galaxy populations (i.e., disparate observational selection techniques detect distinct galaxy populations).

Satellite fraction
In the above discussions on HOD modeling, we limit ourselves to the HOD model form expressed by Eqs. 5 and 6, which is rather restrictive. The underlying assumption of the model is that the center of the halo with mass M h > M min is always occupied by one galaxy in the sample (or at least at a M h -independent constant probability). This form may be appropriate for instance, for luminosity or stellar mass thresholding samples, but there is no reason that this has to be the case for samples selected by other criteria.
We note that the inferred value of f sat is sensitive to the form of the parameterized model of the central and satellite HODs. In this work and in the literature, a power-law form of the satellite HOD is customarily assumed. In this case, a lower α would increase the model N s (M h ) at the lower M h end, near M min , and yield fewer satellites in higher mass halos. Since the halo mass function drops with increasing mass, f sat is mainly determined by the HOD behavior around M h ∼ M min ∼ 10 10.5 h −1 M , where the halo mass function is large and the virial radius is r vir ≈ 0.08 h −1 Mpc at z ∼ 3.8 (Zheng et al. 2007). These scales are too small to be well constrained by our observations. Our observed one-halo term mainly constrains the satellite fraction at larger mass halos (M h ∼ M min ∼ 10 13 h −1 M , where r vir ≈ 0.5 h −1 Mpc at the same redshift). Thus, the f sat values from the HOD modeling should be viewed with caution and may well reflect the artefacts of the assumed form of the model. On the other hand, the sheer presence of a significant one-halo term indicates the existence of some satellites at higher halo masses. The extent of the one-halo term up to R i j ≈ 0.5 h −1 Mpc shows that there are indeed satellites up to M h ∼ 10 13 h −1 M .
In spite of the above caveats, the small satellite fraction of the LAEs is likely to be robust. The small f sat values for the assumed HOD model indicate that not only central-satellite pairs are rare, but also satellite-satellite pairs are as well, suggesting that only a small fraction of halos contain multiple LAEs. The small M min values themselves are also an indication that a large majority of the halos (at the low mass end) that contain a LAE are indeed dominated by one galaxy and in this case, the LAE is probably the central galaxy.

Implications
The clustering results of this study do not only have implications on the baryonic-DM relation, but also on evolving Lyα luminosity functions, signatures of incomplete reionization, and halo mass-dependent Lyα escape fractions. We address these aspects in the following.
The relation between halo mass (or clustering strength) and Lyα luminosity (Table 4 and Fig. 8) demonstrates that highluminosity LAEs tend to reside in higher density environments than lower luminosity ones. As a result, overdense regions contain a larger fraction of high-luminosity sources (and a lower fraction of less luminous ones) than environments of lower density. These inferences affect the Lyα LF measurements at 3 < z < 6. While we expect a shallower faint-end slope of the Lyα LF in overdense regions, the slope should steepen in average or low density environments. As a consequence, surveys for relatively high-luminosity (L Lyα ≈ 10 42 erg s −1 ) LAEs are implicitly biased against the lowest density regions and thus gives a biased shape for the LF, which should not be extrapolated towards lower Lyα luminosities.
Assuming that our L Lyα − M h relation still holds at higher redshifts, the Lyα LF at z ≥ 6 would be even more affected, not only because of the above discussion but also because higher redshift bins are mainly populated by high-luminosity sources, contrary to lower redshift bins (typical case for telescopes with higher sensitivity at bluer wavelengths). Thus, it is important to be careful when interpreting Lyα LFs, especially near the epoch of reionization (EoR), where a shallow to steep variation in the slope of the LF from higher (z ≈ 7) to lower redshifts (z ≈ 5.7) is commonly interpreted as a sign of incomplete reionization (Konno et al. 2014;Matthee et al. 2015;Santos et al. 2016).
Simulations at those higher redshifts also tend to find that high-luminosity LAEs are more likely to be observed than lowluminosity ones because they are able to ionize their surroundings and form H ii regions around them (i.e., ionized bubbles; Matthee et al. 2015;Hutter et al. 2015;Yoshioka et al. 2022). These allow Lyα photons to redshift out of the resonance wavelength and escape the region. Lower luminosity LAEs are then observed if they reside within the ionized bubbles of higher luminosity LAEs or if they are able to transmit enough flux through the IGM (Matthee et al. 2015). If our L Lyα − M h relation is still valid at these redshifts, our results would support this simulation paradigm since high-luminosity LAEs (situated in overdense regions) could form large ionized bubbles more efficiently than low-luminosity sources which tend to be located in lower density environments (Tilvi et al. 2020).
Theoretical studies (e.g., Furlanetto et al. 2006;McQuinn et al. 2007) have modeled the size distribution of these H ii regions and predicted an increase in the apparent clustering signal of LAEs towards the epoch of reionization (i.e., towards a more neutral IGM). Large ionized bubbles become rarer as the ionizing fraction declines. This patchy distribution of H ii regions, which mostly surrounds large galaxy overdensities, boosts the apparent clustering of LAEs. This is commonly interpreted as another sign of incomplete reionization (e.g., Matthee et al. 2015;Hutter et al. 2015). Comparisons between observed intrinsic LAE clustering and model predictions have therefore been used to infer the fraction of neutral hydrogen at the EoR (e.g., Ouchi et al. 2018). Nevertheless, if the clustering dependence on Lyα luminosity continues to z ≈ 6, this comparison should be performed with caution. Because the observed high redshift bins (z ≥ 6) mainly contain high-luminosity LAEs, a strong clustering signal at z ≈ 6 may be wrongly interpreted as incomplete reionization when, in fact, it may only reflect the natural relation between Lyα luminosity and clustering strength.
We speculate that our results also play a role in the amount of escaping Lyα photons (Lyα f esc ). Durkalec et al. (2018) observed a dependence between halo mass and absolute UV magnitude (M UV ). The interpretation of their relation goes as follows: M UV traces star formation rate (SFR; e.g., Walter et al. 2012), which, in turn, tracks stellar mass (M * ; e.g., Salmon et al. 2015), which correlates with halo mass (e.g., Moster et al. 2010). Because we observe a similar relation of M h with L Lyα , L Lyα is presumably also a tracer of star formation. If this is correct, the object-toobject variations in Lyα escape fraction cannot be so large that they obscure the trend of SFR -M * -M h . Given the typical Lyα luminosities of our sample, this is in agreement with the model suggestions of Schaerer et al. (2011a); Garel et al. (2015), where the Lyα f esc is of the order of unity for sources with SFR ≈ 1 M yr −1 . The Lyα luminosity would then be a good tracer of the SFR for less luminous LAEs.

A. Effect of different fields on the clustering measurements
In this work, we have analyzed the clustering of LAEs in the full MUSE-Wide sample, including the CANDELS/COSMOS fields and the HUDF parallel fields. Here, we explore the possible effects on the MUSE-Wide clustering results when including or excluding various sets of fields. In appendix A of Herrero Alonso et al. (2021), we showed that the HUDF parallel fields did not alter the clustering results, their exclusion or inclusion mainly affected the clustering uncertainties. We therefore explore the effect of including the CANDELS/COSMOS region by comparing the clustering of the full MUSE-Wide survey with that present in a subsample without the CANDELS/COSMOS fields. The number of LAEs in the CANDELS/COSMOS region is 250. It is clear from Fig. A.1 that the clustering in both samples is in good agreement. The large-scales bias factors derived from the two curves are indistinguishable (within 1σ). The uncertainties corresponding to the smaller sample are (on average) 20% larger than in the full MUSE-Wide sample. We conclude that the inclusion of these fields has no notable effect on our clustering results but helps in reducing cosmic sample variance uncertainties.

B. Covariance matrix
A common approach to quantify the correlation of the clustering data points is to resample the set of galaxies with the jackknife technique, followed by the calculation of the covariance matrix. To apply the jackknife method, we find a compromise between the number and the size of the jackknife zones. Thus, we split the sky area into ten independent regions (see Fig. B.1) with a spatial extent of ≈ 4 h −1 Mpc in both RA and Dec directions. We then construct ten different subsamples, each of them excluding one jackknife zone, and compute the K-estimator in each subset. These measurements are then used to build up the covariance matrix using Eq. 1 (see Sect. 3.2.1). Considering that the probability of one galaxy pair to contribute to various adjacent bins is higher than that to contribute to several distant bins, one would naively expect a higher correlation in the former case. This is indeed what the (normalized) covariance matrix reflects in the left panel of Fig. B.2. In fact, the noise in the matrix elements corresponding to notably separate bins is substantial. In the right panel of Fig. B.2, we plot the normalized matrix elements as a function of bin i for each bin j to better illustrate the high level of noise in the matrix, especially for bins i > 6, where most curves become negative. This is likely due to the limited spatial size of the survey, which does not allow neither for a higher number of jackknife zones nor for spatially larger zones.
As a result of the considerable noise in the matrix on account of barely correlated bins significantly apart from each other, the minimization of the χ 2 values (Eq. 2) including the full covariance fails (i.e., various χ 2 values become negative). We therefore limit the use of the covariance matrix to its main diagonal and two adjacent diagonals (see red section in the left panel of Fig. B.2; our so-called reduced covariance matrix). This means we set the negative part of the curves in the right panel of Fig. B.2 to zero (i.e., no correlation between those bins), in an attempt to smooth out the noise. While incorporating more diagonals results mathematically problematic for the χ 2 minimization, we have verified that the number of adjacent diagonals (one or two) slightly modifies the χ 2 values but the probability contours represented in Fig. 6 remain unaltered. Thus, so do the best-fit HOD parameters.   Despite current limitations, jackknife is still the most robust method to compute the K-estimator uncertainties. While galaxy bootstrapping or Poisson error bars do not account for bin to bin correlations, our reduced covariance matrix only neglects the correlation between bins remarkably separated (expected to be minimal), but accounts for the correlation between nearby bins.

C. Error estimation comparison
In order to quantify the correlation between the K-estimator bins, the covariance matrix must be computed. By splitting the sky area into independent regions, following the jackknife resampling technique, we create as many subsamples from the MUSE-Wide sample as jackknife zones (see Sect. 3.2.1). The Kestimator is then computed in each subset and the measurements are used to quantify the covariance matrix, whose diagonal provides the variance of each clustering data point. The square root of the diagonal represents the 1σ uncertainties and are represented in blue in Fig. C.1 (same along the main paper).
The jackknife resampling method requires a division of the sky area into several independent regions, each of which should ideally be large enough to cover the full range of scales under consideration. Out of the three samples examined in this study, this can only be partially achieved in the MUSE-Wide dataset. We find that Poisson (bootstrapping) errors are, on average, 7% (46%) larger than those computed with the jackknife technique. These findings corroborate the results from Norberg et al. (2009), who found that the bootstrapping approach overestimates the uncertainties.
Similarly as for the MUSE-Wide survey, we find that bootstrapping uncertainties are ≈ 40% (on average) larger than Poisson in both MUSE-Deep and MXDF. We thus decide to use Poisson errors for the deeper samples in an attempt to least overvalue the uncertainties.
We verified that the error estimation method does not significantly affect our clustering results. The best-fit parameters from MUSE-Deep and MXDF using bootstrapping error bars and the χ 2 minimization described in Sect. 3.1.3 of Herrero Alonso et al. (2021) are consistent with those delivered from Poisson statistics. Although in agreement, bootstrapping delivers ≈ 45% larger uncertainties than Poisson for the best-fit HOD parameters.
We last perform the same experiment in MUSE-Deep and MXDF but using scaled Poisson error bars. We decreased the Poisson in 7% (excess found in MUSE-Wide) and find that the best-fit parameters are ≈ 10% less uncertain than if Poisson errors are directly applied.

D. Dependence of HOD parameters on the shape of the K-estimator
Here we visualize and qualitatively describe the effect of the HOD parameters on the K-estimator. Figure D.1 shows the Kestimator for numerous HOD models. Each panel represents the result of varying one HOD parameter with the other two parameters fixed. Before detailing the major effects, it should be pointed out that the exact change in the shape of the K-estimator does not only depend on the varied parameter but also on the specific choice of the other two. Hence, these panels should merely be seen as illustrative examples. The left panel of Fig. D.1 shows the dependence of the Kestimator on M min . Higher values of log M min (i.e., more massive halos) raise the expected K-estimator at all R i j scales (one-and two-halo terms). At large scales, this occurs because more massive halos present larger bias factors, whereas at small scales, this is due to the decline in the contribution from less massive DMHs.
The middle panel of Fig. D.1 shows the dependence of the Kestimator on M 1 /M min . Larger log(M 1 /M min ) values (i.e., more massive halos) reduce the one-halo term clustering amplitude because of the decrease in the contribution from less massive DMHs. The clustering in the two-halo term does not depend on M 1 .
The right panel of Fig. D.1 shows the dependence of the Kestimator on α. Higher values of α increase the fraction of galaxies in massive DMHs with respect to smaller mass DMHs. Given that more massive halos are more strongly biased, the amplitude of the two-halo term increases. The change observed in the one-halo term is explained because galaxies hosted by massive DMHs can contribute to the one-halo term on its largest scales, while galaxies residing in less massive halos can only contribute to the one-halo term at smaller R i j scales. Since α modifies the fraction of galaxies in massive DMHs to less mass DMHs, the corresponding fraction of the clustering contribution also varies. This alters the slope of the one-halo term.