Synergistic Radiative Transfer Modeling of $\rm Mg\,{\textsc{II}}$ and Ly$\alpha$ Emission in Multiphase, Clumpy Galactic Environments: Application to Low-Redshift Lyman Continuum Leakers

In general, there are four primary mechanisms responsible for Mg ii doublet emission in the ISM / CGM (e.g., Burchett et al., 2021; Chang & Gronke, 2024; Seon, 2024):

(1) Recombination in central H ii regions: High-energy UV photons ($\gtrsim$ 15 eV) from massive stars (e.g., O and B types) can doubly ionize magnesium atoms, creating Mg²⁺ ions. These ions can then recombine to produce Mg ii doublet emission, exhibiting a doublet ratio ($F_{\rm 2796}/F_{\rm 2803}$) of approximately 2:1.

(2) Collisional excitation: Mg⁺ ions can be collisionally excited by electrons near the H ii regions, followed by de-excitation, which results in Mg ii doublet emission. This process also yields a doublet ratio close to 2:1.

(3) Continuum pumping: UV continuum photons near 2800 Å can excite Mg⁺ ions, leading to de-excitation and subsequent Mg ii emission. In this scenario, the spectrum often displays a P-Cygni profile, characterized by an absorption trough and an emission peak atop the continuum.

(4) Ionization by diffuse FUV radiation: FUV continuum radiation in the diffuse warm neutral medium (WNM) at around 1620 Å can ionize magnesium atoms to form Mg⁺ ions. These ions, when collisionally excited and then de-excited by electrons, produce Mg ii doublet emission, again with a doublet ratio close to 2:1.

In this work, we focus primarily on the first three mechanisms that produce Mg ii doublet emission from the central galaxy within the ISM, which either generate line emission with a 2:1 ratio or continuum emission. Our subsequent modeling suggests that such a choice is sufficient to reproduce nearly all of the observed Mg ii emission spectra. Therefore, in line with Occam’s Razor, we believe that while additional mechanisms may contribute, their impact is likely minor.

We model a multiphase, clumpy gaseous medium that resembles the physical structure of the CGM of low-$z$ LyC leakers, as illustrated in Figure 1. We assume that Mg ii and Ly$\alpha$ emission are produced by the central galaxy and propagate outward through such a two-phase CGM, which consists of cool ($\sim 10^{4}\,\rm K$), outflowing clumps and a hot ($\sim 10^{6}\,\rm K$), diffuse, outflowing inter-clump medium (ICM). The Mg⁺ ions are assumed to reside solely in the cool clumps, as they are likely to be fully ionized in the hot medium. In contrast, H i atoms are assumed to exist both in the clumps, with high column densities ($\gtrsim 10^{17}\,\rm cm^{-2}$), and in the ICM, with much lower column densities ($\sim 10^{15}\,\rm cm^{-2}$). We note that although in reality, Mg⁺ and H i may be relatively well-mixed in the cool clumps, their properties (such as outflow velocities and column densities) will be modeled independently – we will focus on Mg ii RT in Section 4 and Ly$\alpha$ RT in Section 5, respectively.

In each model, we further assume that both the clumps and the CGM halo are spherical, and set the clump radius to $R_{\rm cl}=100$ pc and halo radius to $R_{\rm halo}=10$ kpc¹¹1Note that these physical sizes are generally rescalable in idealized RT simulations., based on the typical physical extent of the observed Ly$\alpha$ and Mg ii emission for this sample (Henry et al., 2015, 2018). The clumps are distributed following a power-law $n_{\rm cl}\propto r^{-2}$, with the total number determined by the clump volume filling factor, $F_{\rm V}$. This volume filling factor can be converted into the clump covering factor (i.e., the average number of clumps per sightline), $f_{\rm cl}$, using the relation $f_{\rm cl}=\frac{3}{4}F_{\rm V}R_{\rm halo}/{R_{\rm cl}}$ (Li et al., 2024). In our standard model, the clumps are assumed to have a constant Mg ii column density, $N_{\rm MgII,\,cl}$, and the average total Mg ii column density along a given sightline is $N_{\rm MgII,\,tot}=\frac{4}{3}f_{\rm cl}N_{\rm MgII,\,cl}$ (Dijkstra & Kramer, 2012; Gronke & Dijkstra, 2016).

The motion of the clumps is characterized by two modes: a microscopic turbulent motion²²2Such internal turbulence can result from momentum transfer from the hot, ambient medium (see e.g., Nikolis & Gronke 2024)., represented by the clump Doppler parameter $b_{\rm D,\,cl}$, and a macroscopic radial outflow. For the latter mode, instead of using simplistic velocity functions (e.g., linear or power-law), we adopt a realistic clump velocity profile where the clumps are accelerated by an $r^{-\alpha}$ force and decelerated by the gravitational pull from an NFW-type dark matter halo. Specifically, the kinematic equation governing an outflowing clump is given by:

The solution to this equation is (see Li et al. 2024 for details):

where $M_{\rm vir}$ is the virial mass of the dark matter halo, $r_{\rm vir}$ is the virial radius, $c$ is the halo concentration parameter, $r_{\rm s}=r_{\rm vir}/c$ is the halo scale radius, $r_{\rm min}$ is the clump launch radius (we have fixed it at 1 kpc), and $\mathcal{V_{\infty}}$ is the asymptotic outflow velocity in the absence of gravitational deceleration. For simplicity, in our modeling, we adopt a typical dark matter halo mass of $M_{\rm vir}=10^{11}M_{\odot}$ and an average redshift of $z=0.2$ for our sample. Consequently, the outflow velocity profile depends on only two free parameters: $\mathcal{V_{\infty}}$ and $\alpha$. Varying combinations of $\mathcal{V_{\infty}}$ and $\alpha$ yield different velocity profiles and clump maximum outflow velocities ($v_{\rm cl,\,max}$). We present several example $v_{\rm cl,\,out}(r)$ profiles in Appendix A.

In addition to the resonant scattering of Mg ii, we account for the effect of dust absorption and scattering within the clumps. The dust scattering albedo near $\sim 2800$ Å is approximately $\alpha_{\rm d}$ = 0.57, with an asymmetry factor of $g=0.55$ in the Henyey-Greenstein scattering phase function (Draine, 2003a). We use the dust absorption optical depth in each clump, $\tau_{\rm d,\,cl}$, to characterize the amount of dust absorption³³3Note that the difference in the dust absorption optical depth of the K and H transitions are negligible. We have also accounted for the fact that the ratio of $\tau_{\rm d,\,cl}$ for Ly$\alpha$ and Mg ii is 9.26 in the SMC dust model (this ratio varies across different dust models; see Chang & Gronke 2024).. We also employ a fiducial Small Magellanic Cloud (SMC) dust model, as described by Draine (2003a, b).

Before fitting the observed Mg ii spectra with the RT models, we need to introduce several additional parameters. First of all, since most of the observed Mg ii spectra exhibit positive equivalent widths (EW) that suggest net emission, we introduce the ratio of line photons to continuum photons, $R_{\rm line}=N_{\rm line}/N_{\rm continuum}$⁴⁴4$R_{\rm line}$ can be converted to an intrinsic emission EW using the following relation: ${\rm EW}_{\rm int,\,K+H}=R_{\rm line}2\Delta v{\lambda_{0}}/c\simeq 28R_{\rm line}\rm(\AA)$., to characterize the source function of the Mg ii emission. For each model configuration, two separate RT simulations are performed – one with 10,000 line photons (in a 2:1 ratio for the K and H transitions) and another with 10,000 continuum photons near 2800 Å. In the first case, the intrinsic line emission is modeled as two Gaussian functions centered at the K and H line centers with $\sigma=100\,\rm km\,s^{-1}$. In the second case, the intrinsic emission is assumed to be a flat continuum within $\Delta v=\pm 1500\,\rm km\,s^{-1}$ from the K line center⁵⁵5Here “intrinsic” refers to the line profiles emerging from the ISM, meaning that we essentially assume that all RT effects come from the outflowing gas in the CGM, while the ISM only contributes to the initial line broadening.. Since the RT of each photon is independent, we can construct a composite model spectrum by combining the continuum photons with a fraction of the line photons, characterized by $R_{\rm line}$⁶⁶6Our subsequent modeling shows that all galaxies, except J1648+4957, require $R_{\rm line}<1$. . We also introduce an aperture correction factor to account for aperture losses, which may arise due to the limited slit size of spectrographs used in Mg ii observations (see e.g., Scarlata & Panagia 2015). To incorporate this effect, we define a parameter, $b_{\rm max}$, representing the maximum impact parameter within which photons are included in constructing the model spectra⁷⁷7Note that our treatment of aperture loss is somewhat idealized, as the apertures of slit spectrographs are typically non-circular in reality..

In addition, we incorporate a velocity shift parameter, $\Delta v$, to account for any potential offset between the systemic redshift of Mg ii emission in the model and the observed systemic redshift derived from nebular emission lines. We also include a scaling factor, $f_{\rm scale}$, to provide the model flexibility in matching the continuum level of the data. Our results show that $\Delta v$ is typically small and falls within the observational uncertainties, while $f_{\rm scale}$ is always close to 1, as both the models and the data are normalized prior to comparison. As a result, our modeling has 10 free parameters in total: the clump volume filling factor $F_{\rm V}$, the clump Mg ii column density $N_{\rm MgII,\,cl}$, the clump Doppler parameter $b_{\rm D,\,cl}$, the clump asymptotic outflow velocity $\mathcal{V_{\infty}}$, the clump acceleration power-law index $\alpha$, the ratio of line photons to continuum photons $R_{\rm line}$, the photon maximum impact parameter $b_{\rm max}$, the clump dust absorption optical depth $\tau_{\rm d,\,cl}$, the velocity shift $\Delta v$, and the continuum scaling factor $f_{\rm scale}$.

Our fitting pipeline utilizes the python nested sampling package dynesty (Skilling, 2004, 2006; Speagle, 2020). To avoid the computational expense of real-time RT calculations, the fitting process relies on a pre-computed grid of Mg ii RT models. For each parameter space point visited, the model spectrum is computed through a parameter-weighted multi-dimensional linear interpolation the grid of Mg ii RT model spectra. The resulting spectrum is then convolved with a Gaussian function (FWHM = $50\,\rm km\,s^{-1}$) to simulate the finite instrumental resolution before being compared to the observed Mg ii spectrum⁸⁸8For the four objects observed with HET, the spectral resolution is relatively low (Xu et al., 2023). We therefore apply a convolution with FWHM = $125\,\rm km\,s^{-1}$ in these cases..

For each object, we start by performing a test fit using a fiducial model grid. After examining the posterior distribution, we expand the grid along different dimensions as needed for individual objects. This approach allows us to account for the wide range of best-fit parameters across the parameter space. By customizing the model grid for each object, we ensure that the best-fit parameter values remain within the initial prior bounds. We present the parameter values (i.e., prior ranges) for the fiducial model grid in Table 1.

We have identified an intriguing trend between the best-fit maximum clump radial outflow velocity $v_{\rm MgII,\,max}$ and the total Mg ii column density $N_{\rm MgII,\,tot}$ inferred from the RT modeling. Here $v_{\rm MgII,\,max}$ is defined as the maximum velocity of the clump radial velocity profile (see Appendix A for examples), whereas $N_{\rm MgII,\,tot}$ is defined as the average total Mg ii column density along a sightline, which is the product of the average number of Mg ii clumps along a sightline and the Mg ii column density in each clump.

In Figure 2, we plot all 33 modeled LyC leakers on the $v_{\rm MgII,\,max}$ – $N_{\rm MgII,\,tot}$ plane, categorized into four distinct colors based on their LyC escape types. We find a necessary (though not sufficient; see discussion below) condition for a LyC leaker to be a strong leaker: a high maximum clump radial outflow velocity ($v_{\rm MgII,\,max}\gtrsim 390\,\rm km\,s^{-1}$) and a low total $\rm Mg\,{\textsc{II}}$ column density ($N_{\rm MgII,\,tot}\lesssim 10^{14.3}\,\rm cm^{-2}$).

We now discuss the behavior of each of the four types of LyC leakers and explore the connection between their best-fit parameters and the morphology of their Mg ii spectra individually (best-fits and parameters presented in Figures 3 - 6):

• •

  Strong leakers: The nine strong leakers are all located in the upper left corner of Figure 2, defined by $v_{\rm MgII,\,max}\gtrsim 390\,\rm km\,s^{-1}$ and $N_{\rm MgII,\,tot}\lesssim 10^{14.3}\,\rm cm^{-2}$. In terms of spectral morphology, all strong leakers exhibit an absorption trough extending beyond $v\approx-400\,{\rm km\,s}^{-1}$, indicating their high outflow velocities. However, the Mg ii column density does not appear to be directly associated with any specific spectral feature, making RT modeling necessary for accurate extraction of its value. We note that J1517+3705 was observed with HET/LRS2, which has a particularly low resolution (FWHM $\simeq 125\,\rm km\,s^{-1}$), resulting in a $N_{\rm MgII,\,tot}$ of $\sim 10^{14.3}\,\rm cm^{-2}$ with large uncertainties. Aside from this object, all other strong LyC leakers have $N_{\rm MgII,\,tot}\lesssim 10^{13.7}\rm\,cm^{-2}$.

• •

  Moderate leakers: Among the six moderate leakers, two (J1301+5104 and J1038+4527) are located in the same regime as the strong leakers. This suggests that the distinction between moderate and strong leakers may not be clear-cut; instead, they could represent the same galaxy population with similar gas properties, with only slight differences in LyC leakage due to subtle environmental variations and observational uncertainties. The remaining four moderate leakers, however, all exhibit high Mg ii total column densities with $N_{\rm MgII,\,tot}\gtrsim 10^{14.1}\,\rm cm^{-2}$.

• •

  Potential leakers: The five potential leakers exhibit a broad range of total Mg ii column densities but share a common characteristic: they all have very low maximum clump outflow velocities ($v_{\rm MgII,\,max}\lesssim 320\,\rm km\,s^{-1}$). This is evident from their spectral morphology: the absorption trough at $v<0$ is either negligible (as seen in J1648+4957 and J1133+6513) or is situated relatively close to the line center ($|v_{\rm trough}|\lesssim 300\,\rm km\,s^{-1}$).

  

  

• •

  Non-leakers: Among the thirteen non-leakers, one (J0723+4146) is actually situated in the same regime as the strong leakers. It remains unclear whether this non-leaker does share certain characteristics with the strong leakers or this is an artifact caused by its relatively noisy spectrum and/or the low-resolution observation of J1517+3705 (see discussion above). The rest of the non-leakers appear relatively scattered in the $v_{\rm MgII,\,max}$ – $N_{\rm MgII,\,tot}$ plane but can be categorized into two subgroups. One subgroup has fairly low maximum clump outflow velocities ($v_{\rm MgII,\,max}\lesssim 360\,\rm km\,s^{-1}$), while the other exhibits high total Mg ii column densities Mg ii column densities $N_{\rm MgII,\,tot}\gtrsim 10^{14.5}\,\rm cm^{-2}$. Morphologically, the absorption trough is either insignificant or situated relatively close to the line center ($|v_{\rm trough}|\lesssim 300\,\rm km\,s^{-1}$), or it is very deep and extended due to large total Mg ii column densities (e.g. J1346+1129 and J1314+1048).

Among all 33 modeled objects, we note that the spectrum of J0826+1820 is of poor quality, resulting in a less reliable model fit. We therefore do not fully trust the best-fit parameters for this object and use a hollow black circle in Figure 2 to represent it. In addition, two objects – J1235+0635 (a moderate leaker) and J1244+0215 (a non-leaker) – exhibit an intriguing “quadruple peak” spectral morphology (see Figures 4 and 6) due to strong absorption at the line center of both the Mg ii H and K transitions. This phenomenon can be attributed to their inferred combination of very high Mg ii total column density ($N_{\rm MgII,\,tot}\gtrsim 10^{14.4}\,\rm cm^{-2}$) and extremely low clump outflow velocity ($v_{\rm MgII,\,max}\lesssim 30\,\rm km\,s^{-1}$). We will explore this point in more detail in the next section.

To gain a deeper understanding of the results from Mg ii RT modeling, this section explores the effect of varying individual model parameters on the Mg ii spectra using a suite of RT simulations. In addition to analyzing the shape of Mg ii line profiles, we calculate the equivalent width ratio, $\mathcal{R}=\rm EW_{2796}/EW_{2803}$, for each profile to characterize the relative strength of the K and H lines. In the limit where the spectrum is purely emission, $\mathcal{R}$ converges to the doublet flux ratio $F_{\rm 2796}/F_{\rm 2803}$, which has been frequently used in prior studies (e.g., Chisholm et al. 2020; Chang & Gronke 2024; Seon 2024).

In the following, we vary each key model parameter individually while keeping the others fixed to examine their effects on the Mg ii model spectra:

• •

  Clump maximum outflow velocity $v_{\rm MgII,\,max}$: We explore the effect of this parameter in two different Mg ii total column density regimes. The first example has $N_{\rm MgII,\,tot}=10^{13.0}\,\rm cm^{-2}$, representing strong leakers such as J1158+3125, J1033+6353, J1410+4345, and J1327+4218. The second example has $N_{\rm MgII,\,tot}=10^{14.0}\,\rm cm^{-2}$, representing non-leakers like J1129+4935, J0129+1459, J0723+4146, and J1314+1048. As shown in Figure 7, the major difference between the two cases is that the absorption troughs are generally much shallower and narrower in the low column density scenario, where $v_{\rm MgII,\,max}$ is primarily constrained by the relative intensity of the two peaks (see Appendix B).

  Overall, increasing $v_{\rm MgII,\,max}$ tends to suppress both emission peaks (to a greater extent for the main peak at the K transition) and reduces the equivalent width ratio $\mathcal{R}$. It also shifts the absorption trough away from the line centers of the K and H transitions, with $v_{\rm trough}\sim-v_{\rm MgII,\,max}$, since that is where resonant absorption happens and optical depth is the greatest. In addition, an interesting phenomenon occurs when $v_{\rm MgII,\,max}$ approaches the velocity separation between the K and H transitions ($\sim 770\,\rm km\,s^{-1}$). In this case, the H peak gets enhanced as some K photons are scattered into the H transition, causing a sharp drop in $\mathcal{R}$. Moreover, the absorption trough associated with the H transition becomes blurred and ill-defined (see also Chang & Gronke 2024).

• •

  Clump Mg ii total column density: This parameter is defined as the average total Mg ii along a sightline: $N_{\rm MgII,\,tot}=\frac{4}{3}f_{\rm cl}N_{\rm MgII,\,cl}$, where $f_{\rm cl}=\frac{3}{4}F_{\rm V}R_{\rm halo}/{R_{\rm cl}}$ is the average number of clumps per sightline and $N_{\rm MgII,\,cl}$ is the Mg ii column density of each clump. Here we also consider two examples with different $v_{\rm MgII,\,max}$ values in Figure 8: 350 and 740 $\rm km\,s^{-1}$, respectively. In both cases, increasing $N_{\rm MgII,\,tot}$ deepens the absorption troughs without shifting their locations in velocity space, and also decreases the doublet EW ratio $\mathcal{R}$.

  Nevertheless, we have observed a very interesting phenomenon: in the low $v_{\rm MgII,\,max}$ example, increasing $N_{\rm MgII,\,tot}$ enhances both emission peaks, as no photon exchange occurs between the K and H transitions, and the deepening of the absorption trough is compensated by the increase of the associated emission peak. In contrast, in the high $v_{\rm MgII,\,max}$ example, increasing $N_{\rm MgII,\,cl}$ actually suppresses the K peak while enhancing the H peak, due to more K photons being scattered into the H transition.

• •

  The Doppler parameter of each clump $b_{\rm D,\,cl}$: This parameter represents the turbulent velocity within each clump. As shown in Figure 9, increasing $b_{\rm D,\,cl}$ tends to broaden and strengthen both the emission peaks and absorption troughs, as the perturbed velocity field of the Mg ii gas causes photons across a wider frequency range to participate in scattering. The effect on the doublet EW ratio $\mathcal{R}$ is relatively minor – interestingly, $\mathcal{R}$ slightly increases as $b_{\rm D,\,cl}$ rises from 8 $\rm km\,s^{-1}$ to 26 $\rm km\,s^{-1}$ due to a more significant decrease in the EW for the H transition. However, as $b_{\rm D,\,cl}$ continues to increase, $\mathcal{R}$ decreases again.

• •

  The relative strength of line emission versus continuum emission $R_{\rm line}$: This parameter is defined as the ratio of line photons to continuum photons in an RT simulation, reflecting the relative contribution of the two emission mechanisms. For example, $R_{\rm line}=0$ corresponds to pure continuum emission, while $R_{\rm line}=1$ indicates equal numbers of line and continuum photons are emitted. As $R_{\rm line}$ increases, the spectrum gradually shifts from being continuum-dominated to line-dominated – this is evident as the continuum level decreases and the emission peaks become more prominent. The EWs of both the K and H transitions increase, with the K transition experiencing a greater enhancement, leading to an overall increase in the EW ratio $\mathcal{R}$.

• •

  The clump dust absorption optical depth $\tau_{\rm d,\,cl}$: In our modeling, we have assumed the presence of dust within the clumps, which influences the RT of Mg ii photons through dust scattering and absorption. As $\tau_{\rm d,\,cl}$ increases, the emission peaks diminish due to dust extinction, and the absorption troughs become shallower, as dust extinction weakens the Mg ii resonant scattering. In general, the effect of dust is more pronounced for the K line than for the H line, leading to a decrease in the doublet EW ratio $\mathcal{R}$.

• •

  The aperture correction factor $b_{\rm max}/R_{\rm halo}$: This parameter is designed to mimic the aperture loss of Mg ii photons in real observations, where $b_{\rm max}$ represents the maximum impact parameter of the scattered photons included in the model spectra. The effect of this parameter on the Mg ii spectra is generally minor, except when the total Mg ii column density is very high and the clump’s maximum outflow velocity is very low. In this regime, the clump’s optical depth is high enough for photons to travel to large impact parameters.

  Here we use a “quadruple peak” example to illustrate the effect of $b_{\rm max}/R_{\rm halo}$. In this case, photons that are originally near the line center are scattered to large impact parameters, and when $b_{\rm max}/R_{\rm halo}$ is small, these photons are excluded, yielding a deep trough at the line center (even below the continuum level). As $b_{\rm max}/R_{\rm halo}$ increases, the photons with frequencies near the line center are gradually included in the spectra, and the flux at the line center rises above the continuum level again.

Combining the analysis from this section and the previous section, we find that the two most important parameters shaping the Mg ii line profiles are the clump outflow velocity and the total Mg ii column density. The other four parameters analyzed above ($b_{\rm D,\,cl},R_{\rm line},\tau_{\rm d,\,cl},$ and $b_{\rm max}/R_{\rm halo}$) have only subdominant effects on the Mg ii spectra, and there are no significant correlations between their best-fit values inferred from RT modeling and the amount of LyC leakage. More importantly, it is evident that the effects of different model parameters are highly intertwined, suggesting any empirical indices characterizing specific features of a Mg ii spectrum may not be sufficient to deduce the underlying gas properties. For instance, a higher doublet EW ratio $\mathcal{R}$ could result from either a lower clump outflow velocity, a lower clump column density, a higher turbulent velocity within the clumps, a higher fraction of line emission, a lower clump dust absorption optical depth, or simply a greater aperture loss. To accurately determine the physical parameters of the scattering gas in the ISM / CGM, systematic RT modeling of the full Mg ii doublet profile is indispensable, as it is virtually impossible to make any reliable inferences about the underlying gas properties without such simulations.

Although modeling the Mg ii emission can provide insights into the properties of clumpy gas in the ISM / CGM, these gas properties cannot, in principle, be directly used to infer a galaxy’s LyC leakage. This is because Mg ii emission only constrains the properties of Mg ii gas, and the relationship between Mg ii gas and H i gas (which is ultimately responsible for LyC leakage) remains unclear. It is uncertain whether Mg ii is consistently associated with H i in each gas clump; even if they do co-exist, it is unknown whether they share the same kinematics or spatial extent, or if they have a fixed Mg ii-to-H i ratio. Therefore, rather than relying on empirically calibrated Mg ii-H i relations, we conduct RT modeling of Ly$\alpha$ profiles for certain objects in our sample to directly constrain the H i properties in their CGM and further infer their LyC leakage.

We searched for archival HST COS/G160M observations for our sample and found that six galaxies have high-resolution Ly$\alpha$ profiles published (Henry et al. 2015; Yang et al. 2017; Henry et al. 2018, HST Program ID 12928, 14201, 15865). Among these six galaxies, two are strong leakers (J0917+3152 and J0911+1831), three are potential leakers (J1133+6513, J1248+1234, and J1246+4449), and one is a non-leaker (J1244+0215). The continuum-subtracted Ly$\alpha$ profiles of these six objects are shown in Figure 10. At first glance, if we follow the conventional wisdom on inferring LyC leakage from Ly$\alpha$ profiles (e.g., based on shell model results, Verhamme et al. 2015; see also Naidu et al. 2022), one might expect that the galaxy with the most significant flux near the line center (J1248+1234) would be the strongest LyC leaker, with J1133+6513 being the second strongest. However, this is not the case: both are actually potential leakers with $F_{\rm\lambda LyC}/F_{\rm\lambda 1100}<0.05$ detected at $\lesssim 2\sigma$ significance. Meanwhile, the two strong leakers, J0917+3152 and J0911+1831, exhibit virtually no flux at the line center. This phenomenon poses a challenge to the conventional Ly$\alpha$-LyC connection and prompts us to conduct RT modeling to investigate further.

We use a similar setup for Ly$\alpha$ modeling as in our Mg ii modeling, but we also incorporate a hot, inter-clump medium (ICM) that may contain some residual H i capable of scattering Ly$\alpha$ photons (Li & Gronke, 2022; Erb et al., 2023). This hot gas component, assumed to have a temperature of $\sim 10^{6}$ K, is a volume-filling medium characterized by two parameters: the H i number density, $n_{\rm HI,\,ICM}$, and the radial outflow velocity, $v_{\rm ICM}$ (assumed to be constant). We refer the readers to Figure 1 for a schematic representation of the model. The fitting procedure is similar to what we described in Section 2; however, since we do not need $R_{\rm line}$ and $f_{\rm scale}$ for Ly$\alpha$ modeling, each fitting run still contains a total of 10 free parameters. We present the parameter values for the fiducial Ly$\alpha$ model grid in Table 2.

We present the best-fit results using the multiphase, clumpy RT model in Figure 10. While the RT models have successfully reproduced all six Ly$\alpha$ line profiles, several best-fit parameter values are both intriguing and somewhat puzzling, warranting further discussion. The key findings are as follows:

• •

  Clump maximum outflow velocity: We do not observe a clear correlation between the maximum clump outflow velocities inferred from Ly$\alpha$ and Mg ii modeling. For J0917+3152, J0911+1831, and J1133+6513, the velocities inferred by both methods are comparable, suggesting the Mg⁺ and H i gas may track each other kinematically; however, for the other three objects, the outflow velocities inferred from Ly$\alpha$ are at least twice as large as those inferred from Mg ii, and the Mg ii-inferred clump maximum outflow velocity lies between the ICM outflow velocity and the Ly$\alpha$-inferred clump maximum outflow velocity (i.e., $v_{\rm ICM}<v_{\rm cl,\,max,\,MgII}<v_{\rm cl,\,max,\,Ly\alpha}$).

• •

  Clump H i column density: There does not seem to be a clear correlation between the inferred clump H i and Mg ii column densities either. For the two strong leakers, J0917+3152 and J0911+1831, the ratio of the total column densities of Mg ii to H i is approximately $10^{-5.5}$, similar to the Mg abundance [Mg/H] in the warm neutral medium of the Milky Way (Jenkins, 2009). In contrast, J1133+6513 exhibits a much lower Mg-to-H i ratio of $10^{-6.7}$, while J1248+1234, J1246+4449 and J1244+0215 show significantly higher ratios¹⁰¹⁰10Note that a higher Mg-to-H i ratio does not necessarily indicate a higher [Mg/H] abundance, as the column density of H⁺ in the CGM is unknown from RT modeling (Xu et al., 2022b, see e.g., Section 5.3.2 in). of $10^{-4.0}$ to $10^{-2.7}$.

• •

  H i in the hot, inter-clump medium: This component primarily contributes to additional Ly$\alpha$ scattering near the line center. Unsurprisingly, the two objects with the most significant flux near the line center, J1133+6513 and J1248+1234, have the lowest H i column density in the diffuse inter-clump medium ($N_{\rm HI,\,ICM}\lesssim 10^{15}\,\text{cm}^{-2}$). In particular, J1248+1234 has an ICM outflow velocity of $\sim 120\,\text{km}\,\text{s}^{-1}$, which further reduces the H i optical depth at the line center.

• •

  Clump volume filling factor (and therefore clump covering factor): For J0917+3152, J0911+1831, J1248+1234 and J1246+4449, the $F_{\rm V}$ values inferred from both Ly$\alpha$ and Mg ii modeling are consistent, around 2 – 4%. However, for J1133+6513 and J1244+0215, we find large $F_{\rm V}$ values inferred from Mg ii modeling (3% and 4%), but much smaller values inferred from Ly$\alpha$ modeling (0.4% and 0.7%). This discrepancy suggests that the H i gas and Mg⁺ may not be fully co-spatial; the H i gas may reside in smaller clumps, thereby occupying a smaller volume fraction of the halo.

Having completed the modeling of the observed Ly$\alpha$ profiles using the multiphase, clumpy RT models, it is now possible to theoretically predict the LyC escape fraction from the best-fit models, as we have determined the spatial distribution and column densities of the H i clumps. To do this, we inject 10,000 photons at wavelengths far from the Ly$\alpha$ line center (to simulate LyC photons) and record the H i column densities encountered by the photons before they either escape or are absorbed by dust. Since the H i cross-section for LyC photons is $\sim 6.3\times 10^{-18}\,\rm cm^{-2}$, we classify photons that encounter a total H i column density lower than $\sim 1.6\times 10^{17}\,\rm cm^{-2}$ and are not absorbed by dust as successfully escaping. The theoretical LyC escape fraction is then calculated by dividing the number of escaped LyC photons by the total number of injected LyC photons. We find that, given the high best-fit clump H i column densities (all greater than $10^{17}\rm cm^{-2}$), the dominant factor for LyC escape is the clump covering factor, $f_{\rm cl}$. Using the 2-$\sigma$ ranges for the best-fit $F_{\rm V}$ values, along with the corresponding parameter values from the high-likelihood region of the posteriors derived from Ly$\alpha$ modeling and the relation $f_{\rm cl}=\frac{3}{4}F_{\rm V}R_{\rm halo}/{R_{\rm cl}}$, we infer the value and uncertainties of $f_{\rm cl}$ for all six objects. We then run RT simulations using the derived $f_{\rm cl}$ values to estimate the values and uncertainties of the theoretical LyC escape fractions.

We plot the RT-inferred LyC escape fractions, $f_{\rm esc,\,RT}$, against the observed escape fraction determined using three different metrics (observed flux near $\lambda_{\rm rest}=912\,\rm\AA$, H$\rm\beta$, and the UV SED, Flury et al. 2022) in Figure 11. While the two strong leakers do, on average, exhibit higher $f_{\rm esc,\,RT}$ values compared to two of the potential leakers, it is puzzling to observe that the non-leaker J1244+0215 and the potential leaker J1133+6513 actually exhibit the highest RT-inferred LyC escape fractions. We double-checked these two objects and found that their $f_{\rm cl}$ inference is quite robust, consistently remaining below one. This result indicates that the average number of clumps per sightline is less than one, suggesting a low overall gas covering fraction and the existence of channels that are free of high-$N_{\rm HI}$ clumps. Interestingly, at face value, there seems to be no clear connection between their high RT-inferred LyC escape fractions and their spectral morphology – J1133+6513 shows significant flux near the line center, whereas J1244+0215 displays almost no flux at the line center.

We find that the reason J1133+6513 and J1244+0215 favor low $f_{\rm cl}$ values is due to the relatively strong blue peaks in their Ly$\alpha$ spectra. From an RT perspective, the relative intensity of the blue peak compared to the red peak is primarily determined by two parameters: the clump covering factor $f_{\rm cl}$ and the clump maximum outflow velocity $v_{\rm cl,\,max}$. Increasing $f_{\rm cl}$ tends to suppress the blue peak relative to the red peak, as it increases the likelihood of scattering, which, in an outflowing medium, favors the escape of red photons over blue photons. Similarly, increasing $v_{\rm cl,\,max}$ also suppresses the blue peak because, in the reference frame of the outflowing clumps, more blue photons are redshifted toward the line center and require additional scattering to escape. In the case of J1133+6513 and J1244+0215, where the red wings extend to around $1000\rm\,km\,s^{-1}$, their $v_{\rm cl,\,max}$ must be reasonably large, requiring $f_{\rm cl}$ to be low in order to produce the strong blue peak. We present several examples to illustrate this point in Figure 12, where we vary $f_{\rm cl}$ and $v_{\rm cl,\,max}$ individually while keeping all other model parameters fixed.

At this point, it remains unclear why J1133+6513 and J1244+0215 exhibit particularly strong blue Ly$\alpha$ peaks, which result in their notably low gas covering factors. Nevertheless, our experiments suggest that modeling spatially integrated Ly$\alpha$ spectra alone may not be sufficient to accurately infer the LyC escape fractions of LyC leakers. In the future, spatially resolved observations using high-resolution IFU spectrographs may help resolve the degeneracy between model parameters by providing additional spatial information (e.g., Erb et al. 2023). Additionally, certain assumptions in the RT model may contribute to this puzzle, such as the assumption of angular isotropy, whereas in reality, the escape of Ly$\alpha$ and LyC photons may be anisotropic (e.g., Almada Monter & Gronke 2024). We plan to explore these possibilities in future work.

As discussed in Section 4.1, we identified a necessary condition for a LyC leaker to be a strong leaker: a high maximum clump radial outflow velocity ($v_{\rm MgII,\,max}\gtrsim 390\,\rm km\,s^{-1}$) and a low total Mg ii column density ($N_{\rm MgII,\,tot}\lesssim 10^{14.3}\,\rm cm^{-2}$). We now explore the physical reasoning behind this condition.

Our Ly$\alpha$ RT modeling of six objects in our sample reveals that the only two galaxies with large ($>10\%$) theoretical LyC escape fractions exhibit particularly low volume filling factors $F_{\rm V}$ and correspondingly low gas covering factors $f_{\rm cl}$. These objects achieve significant LyC escape not because of lower clump H i column densities (which must exceed $\sim 10^{17}\rm cm^{-2}$ to match the broad peaks of their Ly$\alpha$ profiles), but because certain sightlines encounter only zero or one clump. In other words, our modeling results support a “picket fence” geometry for the CGM (e.g., Heckman et al., 2011; Rivera-Thorsen et al., 2017; Gazagnes et al., 2018, 2020; Saldana-Lopez et al., 2022) rather than a “density-bounded” scenario, where the clump H i column density is not high enough to prevent the penetration of LyC photons. We illustrate this point with a schematic in Figure 13.

The conditions we identified, namely a combination of high $v_{\rm MgII,\,max}$ and low $N_{\rm MgII,\,tot}$, may therefore favor the formation of such a picket-fence-like CGM that facilitates LyC escape. A high Mg ii outflow velocity likely indicates strong supernova feedback, which can blow out the gas clumps and create low-density H i channels (e.g., Li et al., 2015; Fielding et al., 2017; Smith et al., 2018; Sarbadhicary et al., 2022). Moreover, statistically speaking, a lower total Mg ii column density is generally associated with lower H i column densities, further suggesting the presence of these low-density channels. Therefore, galaxies exhibiting both characteristics are understandably more likely to be strong LyC leakers.

While the $v_{\rm MgII,\,max}$ – $N_{\rm MgII,\,tot}$ criterion we present is relatively straightforward and clear, we recommend exercising caution when using it as an indirect indicator of LyC leakage. A larger sample that has both Mg ii and Ly$\alpha$ observations is needed to fully test this criterion. Furthermore, additional high-resolution numerical simulations are essential to deepen our understanding of the physical processes driving this correlation.

A recent study by Carr et al. (2024) employed the semi-analytical line transfer (SALT) model to analyze Mg ii line profiles for a similar sample. The SALT model aims to analytically solve the radiative transfer equation using the Sobolev approximation, which simplifies the treatment of radiative transfer. As such, SALT does not simulate the resonant scattering of photons; its setup also differs from the clumpy RT model used in this work – it assumes a bi-conical, continuous outflowing wind with power-law radial density and velocity profiles, without accounting for turbulent gas motion in the CGM. Given these differences in the modeling approach, it is expected that the SALT model yields different results, such as gas outflow velocities and Mg ii column densities, compared to our findings.

In their study, Carr et al. (2024) modeled 29 galaxies from the sample used in this work, achieving satisfactory fits for 20 galaxies with an outflowing SALT model and for 6 galaxies with a non-outflowing double-Gaussian ISM model. For these six galaxies (J0826+1820, J1033+6353, J1133+6513, J1158+3125, J1235+0635, and J1410+4345), except for J0826+1820 due to its noisy spectrum, our RT model provides a good fit for the remaining five with inferred maximum clump outflow velocities of 533, 272, 432, 30, and 452 $\rm km\,s^{-1}$, respectively. In addition, for two profiles that were not reproduced by the SALT model, our model gives a good fit for J1248+1234 and a decent fit for J1310+2148¹¹¹¹11Note that J1310+2148 exhibits an unusual redshifted absorption feature, possibly caused by an additional inflowing absorber along our sightline; nevertheless, our modeling has successfully captured the overall P-Cyngi-like shape of the line profile.. We find that differences in our data processing procedures may explain the different modeling results – Carr et al. (2024) used a finer binning than the instrumental resolution for their spectra, whereas in this work, we used the resolution-based rebinned version of the Mg ii spectra as presented in Xu et al. (2023). These differences in data processing may result in different levels of significance in the blueshifted absorption troughs, particularly for J1033+6353, J1158+3125, and J1410+4345 – the three galaxies categorized as strong leakers in both studies. This may also explain why no evidence of outflows was found for these galaxies in Carr et al. (2024).

To do a more detailed comparison, we examined the best-fit total Mg ii column densities¹²¹²12We refrain from making a direct comparison of the maximum gas outflow velocities between the two studies, as we have made very different assumptions about the radial velocity and density profiles of the Mg⁺ gas (e.g., Carr et al. 2024 noted that their terminal velocity $v_{\rm\infty}$ may become unconstrained when the density field drops to an undetectable level before the velocity field reaches its maximum). derived by Carr et al. (2024) for the 20 galaxies successfully modeled with their outflowing SALT framework, as shown in Figure 14. We find that the Mg ii column densities reported in their work are systematically higher, with larger uncertainties, compared to those derived from our RT model. However, we note that this difference does not necessarily imply that one model outperforms the other; instead, the different $N_{\rm MgII,\,tot}$ values may result from the different model geometries (i.e., spherically symmetric vs. bi-conic), and the different uncertainties may simply reflect the different volumes of the parameter spaces explored.

Given the more comprehensive treatment of resonant scattering physics in our full RT modeling, we suggest that our multiphase, clumpy model may provide valuable constraints on the underlying gas parameters of the CGM. Nonetheless, the SALT model, with its simpler, semi-analytical framework, can also yield useful results, particularly under the assumption of an anisotropic geometric configuration. Future work will be essential to further assess and compare the strengths and limitations of these two models in different contexts. Further Mg ii observations using high-resolution IFUs, along with comparative studies against cosmological simulations, will also be essential (e.g., Zheng et al. 2010, 2011; Gronke & Dijkstra 2016; Gronke et al. 2018; Blaizot et al. 2023, Jennings et al., in prep; Carr et al., in prep).