2.1. Observations and data reduction

The MAGPI surveyFootnote ^a is an ongoing Large Program on the VLT/MUSE, targeting 56 fields from the Galaxy and Mass Assembly (GAMA; Driver et al. Reference Driver2011) G12, G15 and G23 fields. MAGPI also includes archival observations of legacy fields Abell 370 and Abell 2744. The survey targets a total of 60 primary galaxies with stellar masses $M_{*} \gt \! 7 \times10^{10}$ M $_\odot$ and $\sim$ 100 satellite galaxies with $M_{*} \gt \! 10^{9}$ M $_\odot$ . The primary objective of MAGPI is to perform a detailed spatially resolved spectroscopic analysis of stars and ionised gas within $0.25\lt z \lt0.35$ galaxies (see Foster et al. Reference Foster2021). Data are taken using the MUSE Wide Field Mode ( $1'\times1'$ ) with a spatial sampling rate of $0.2''$ /pixel and the median Full Width at Half Maximum (FWHM) is $0.64''$ , $0.6''$ and $0.55''$ in the g, r and i bands, respectively. Each field is observed in six observing blocks, each comprising $2\times 1\,320$ s exposures, resulting in a total integration time of $4.4$ h. The survey primarily employs the nominal mode (data are taken with the blue cut-off filter in place), providing a wavelength coverage ranging from $4\,700$ to $9\,350$ Å, with a dispersion of $1.25$ Å. Ground-layer adaptive optics (GLAO) is used to correct atmospheric seeing effects, resulting in a gap between $5\,805$ and $5\,965$ Å due to the GALACSI laser notch filter. The MUSE data reduction process is thoroughly described in Foster et al. (Reference Foster2021).

2.2. LAE identification

The depth of the MAGPI data enables the detection of both foreground sources within the Local Universe and distant background sources, including LAEs in the redshift range $2.9 \lesssim z \lesssim 6.6$ . To identify emission lines, we employ LSDCat1.0 Footnote ^b (Herenz & Wisotzki Reference Herenz and Wisotzki2017), a Python-based tool designed to detect faint emission lines in integral-field spectroscopic datacubes by correlating a matched filter with the data. A detection threshold of $\mathrm{SN}_{\mathrm{thresh}} = 7.0$ was adopted to minimise false positives. This choice is in line with the thresholds used in other MUSE-based LAE searches; for instance, the MUSE-WIDE survey applies a threshold of 5–6.4 for LAE detection (see Kerutt et al. Reference Kerutt2022). The initial list of emission-line candidates is then passed through the LAE_Scanner Footnote ^c tool, which filters out spurious and low-significance detections by enforcing criteria on spectral shape, substantially reducing the false-positive rate. For the remaining candidates, segmentation maps were created using the PROFOUND R package (Robotham et al. Reference Robotham2018) and then global spectra are extracted. Redshifts are estimated using the MARZ Footnote ^d software (Hinton et al. Reference Hinton, Davis, Lidman, Glazebrook and Lewis2016), which provides automated redshift identification based on spectral templates. During this step, all candidates undergo visual inspection to confirm the presence and shape of the Ly $\alpha$ line, and to exclude contaminants such as low-redshift [O ii] emitters by checking whether the observed line profile matches the expected [O ii] $\lambda3726,3729$ doublet separation at lower redshifts. Following these steps, we identified 417 LAEs in the first 35 MAGPI fields. A circular aperture of $1.4''$ radius is used to compute spectroscopic properties, chosen to match the typical extent of LAEs while minimising contamination from noisy spaxels. The redshift distribution of Ly $\alpha$ luminosities is presented in Figure 1.

Figure 1. Ly $\alpha$ luminosity as a function of spectroscopic redshift for 417 LAEs identified in first 35 MAGPI fields, spanning the redshift range $2.9 \lesssim z \lesssim 6.6$ .

5.1. Spectral measurements

We measure the following spectroscopic quantities from the intrinsic model parameters, reflecting true profiles before instrumental broadening:

• • Peak velocity separation ( $\Delta_{\mathrm{peak}}$ ): Defined as the velocity difference between the red and blue peaks, $\Delta_{\mathrm{peak}} = V_{\mathrm{red}} - V_{\mathrm{blue}}$ , where both velocities are measured with respect to the systemic redshift, in km s^-1.

• • Total Ly $\alpha$ flux and luminosity: The total Ly $\alpha$ flux is obtained by integrating the flux density under the area of the DAG. Flux uncertainties are drawn from the posterior distributions of the MCMC chains, providing parameter uncertainties that better reflect the full likelihood landscape. The Ly $\alpha$ luminosity, $L_{\mathrm{Ly} \alpha}$ , is derived from the total flux, and its corresponding errors.

• • Integrated signal-to-noise ratio, (S/N) $_{\mathrm{total}}$ : Calculated as the total integrated Ly $\alpha$ flux divided by the flux uncertainty.

• • Blue-to-total flux ratio ( $F_{\mathrm{blue}}\, / F_{\mathrm{total}}$ ): The relative strength of the blue peak is quantified as $F_{\mathrm{blue}} / F_{\mathrm{total}} = F_{\mathrm{blue}} / (F_{\mathrm{blue}} + F_{\mathrm{red}})$ , where $F_{\mathrm{blue}}$ and $F_{\mathrm{red}}$ are the integrated fluxes of the blue and red peaks, respectively.

• • FWHM of each peak: The line width of each peak is calculated using the fitting parameters (see Mukherjee et al. Reference Mukherjee2024):

  (2) \begin{equation} \mathrm{FWHM}\, (\mathrm{km}\, \mathrm{s}^{-1}) = \frac{2\sqrt{2\,\ln \! (2)} \, w}{1 - 2\,\ln \! (2)\, a_{\mathrm{asym}}^2} \end{equation}
  where w and $a_{\mathrm{asym}}$ are the width and asymmetry parameters of the fit, respectively. Uncertainties are derived from errors in the fit parameter returned by emcee.

• • Flux density in the absorption trough ( $F_{\mathrm{trough}}$ ): measured directly from the intrinsic DAG profile as the flux density at the minimum between the two peaks, with uncertainties estimated from MCMC sampling.

• • Ly $\alpha$ red peak asymmetry parameter, $A_{f}$ : This is defined as follows:

  (3) \begin{equation} A_{f} = \frac{\int_{\lambda_{\mathrm{red}}}^{\infty} f_{\lambda} d\lambda}{\int_{\lambda_{\mathrm{trough}}}^{\lambda_{\mathrm{red}}} f_{\lambda} d\lambda} \end{equation}
  where $\lambda_{\mathrm{red}}$ and $\lambda_{\mathrm{trough}}$ are the wavelengths of the red peak and absorption trough, respectively (Rhoads et al. Reference Rhoads2003; Kakiichi & Gronke Reference Kakiichi and Gronke2021).

5.2. Global Ly $\alpha$ profiles

This section describes the global properties of Ly $\alpha$ emission measured from spatially integrated spectra. Full spectroscopic measurements are listed in Table A1 in the Appendix.

5.2.1. Fraction of blue and red-dominated profiles

In our parent double-peaked LAE sample, the total integrated S/N of the Ly $\alpha$ line ranges from 3.52 to 92.11. We find $18/108 (\sim 17$ %) double-peaked profiles are blue-dominated ( $F_{\mathrm{blue}}\, / F_{\mathrm{total}}$ $ \gt 0.5$ ). Three of these with MAGPI IDs 1534253170, 1207170349, and 2302106291 have systemic redshifts that confirm inflows or accretion-dominated sources. Of which two bright cases showing extended halos are reported in Mukherjee et al. (Reference Mukherjee2023). The remaining 90 sources are red-dominated ( $F_{\mathrm{blue}}\, / F_{\mathrm{total}}$ $ \lt 0.5$ ), indicating gas outflows. This blue-to-red fraction is in agreement with the cosmological zoom-in simulation (Blaizot et al. Reference Blaizot2023) of a galaxy, where $\lt$ 20% of the profiles are found to be blue-dominated. However, a recent study of LAEs from the MUSE Extremely Deep Field found a significantly higher fraction ( $\gt$ 40%) (Vitte et al. Reference Vitte2025). However, we note that numerous blue-dominated profiles lacking systemic redshift measurements could also arise from backscattering events (Verhamme et al. Reference Verhamme, Schaerer and Maselli2006), where the observed smaller red bump originates from redshifted photons scattered off the far side of the outflowing medium. In particular, we observe several such cases at $z\gt4$ . Without systemic redshifts, it remains uncertain whether these cases truly represent gas inflows or are instead manifestations of such backscattering effects. Figure 5 shows the blue-to-toal flux ratio versus (S/N) $_{\mathrm{total}}$ , including data from Vitte et al. (Reference Vitte2025). We find that blue peak flux increases with decreasing (S/N) $_{\mathrm{total}}$ , while extreme blue-dominated LAEs with $F_{\mathrm{blue}}\, / F_{\mathrm{total}}$ $\gt 0.8$ are only detected for (S/N) $_{\mathrm{total}} \gtrsim 7$ .

Table 1. Summary of the double-peaked LAE sample used in this work.

Figure 5. Blue-to-total flux ratio as a function of the total S/N of the Ly $\alpha$ line in our parent sample. Pink triangles are taken from Vitte et al. (Reference Vitte2025). The horizontal black dashed line marks the division between blue- and red-dominated regions.

For robust measurements, we adopt a conservative threshold of $(\mathrm{S}/\mathrm{N})_{\mathrm{blue}} \geq 3$ and $(\mathrm{S}/\mathrm{N})_{\mathrm{red}} \geq 3$ throughout the rest of the paper, resulting in a ‘Conservative sample’ of 76 sources (see Table 1). Our conservative sample includes 9 blue-dominated LAEs.

5.2.2. Ly $\alpha$ luminosity versus flux ratio

In Figure 6, we study the relationship between Ly $\alpha$ luminosity and $F_{\mathrm{blue}}\, / F_{\mathrm{total}}$ for our conservative sample. For this, we take sources only up to $z \lt 4.5$ . Sources at higher redshift ( $z \gt 4.5$ ) are excluded, as they are strongly affected by IGM absorption and tend to be more luminous, which could bias the relation. We find a strong anticorrelation (Spearman $\rho = - 0.48$ , p-value $= 4.95 \times 10^{-5}$ ) between Ly $\alpha$ luminosity and $F_{\mathrm{blue}}\, / F_{\mathrm{total}}$ (Figure 6): Brighter Ly $\alpha$ profiles show stronger red peaks, while fainter ones have more flux in the blue peak. We perform a quantile-based binning, taking the same number of sources within each $F_{\mathrm{blue}}\, / F_{\mathrm{total}}$ bin. We find that in the red-dominated region, the bin median luminosity decreases as $F_{\mathrm{blue}}\, / F_{\mathrm{total}}$ increases (see Figure 6). This trend is consistent with cosmological simulations (Blaizot et al. Reference Blaizot2023) using mock Ly $\alpha$ spectra. Blaizot et al. (Reference Blaizot2023) identify this relation in a single simulated galaxy with a relatively constant stellar mass and star formation rate (SFR) during its evolutionary phase. Our results extend this trend to a broader population of double-peaked LAEs across a wide redshift range.

Figure 6. Observed Ly $\alpha$ luminosity versus $F_{\mathrm{blue}}\, / F_{\mathrm{total}}$ for 64 LAEs in our conservative sample at $z \lt 4.5$ , where the IGM has less effect. The vertical black dashed line separates the blue-dominated and red-dominated regions. Orange stars represent the median luminosity within each of the 7 equal-population bins (9 sources per bin), with the scatter in each bin indicated by black error bars.

The relation flattens at lower luminosities and higher ratios. The flattening is likely a selection effect rather than astrophysical. Since each Ly $\alpha$ peak must have S/N $\geq 3$ in our conservative sample, galaxies with low total S/N (i.e. (S/N) $_{\mathrm{total}} \sim 5 - 6)$ can only be included if both peaks are of similar strength, driving $F_{\mathrm{blue}}\, / F_{\mathrm{total}}$ $\to 0.5$ . This bias emerges near the luminosity where (S/N) $_{\mathrm{total}}$ drops below $\sim$ 6 and may also contribute to the apparent slope at higher luminosities and lower ratios.

Figure 7. Redshift distribution of Ly $\alpha$ peak line widths in our conservative sample. The left panel shows the FWHM of the red peak ( $\mathrm{FWHM}_{\mathrm{red}}$ ), while the right panel presents the FWHM of the blue peak ( $\mathrm{FWHM}_{\mathrm{blue}}$ ), both obtained from DAG fitting. In the left panel, the pink triangle markers represent LAEs from the MUSE-Wide and MUSE-Deep surveys (Kerutt et al. Reference Kerutt2022).

5.2.4. Correlations among velocities, peak separation and line width

Our conservative sample spans Ly $\alpha$ peak separations of $\Delta_{\mathrm{peak}} = 223$ – $1\,008\, \mathrm{km}\,\mathrm{s}^{-1}$ , with a mean of $498.98\,\pm\,73.12 \, \mathrm{km}\,\mathrm{s}^{-1} $ . This is consistent with previous findings from MUSE-Wide and MUSE-Deep LAEs ( $481 \,\pm\, 244 \, \mathrm{km}\,\mathrm{s}^{-1} $ ; Kerutt et al. Reference Kerutt2022), the MUSE Extremely Deep Field ( $534 \pm 28 \,\mathrm{km}\,\mathrm{s}^{-1}$ ; Vitte et al. Reference Vitte2025), and $z \sim 2.2$ LAEs ( $500 \pm 56 \,\mathrm{km}\,\mathrm{s}^{-1}$ ; Hashimoto et al. Reference Hashimoto2015). Our sample remains above the minimum measurable threshold of $150\, \mathrm{km}\,\mathrm{s}^{-1}$ reported by Vitte et al. (Reference Vitte2025).

For LAEs with known systemic redshifts, we display the blue and red peak velocities as a function of half the peak separation ( $\Delta_{\mathrm{peak}}/2$ ) in Figure 8, and we notice a strong correlation close to the one-to-one relation. This indicates that double-peaked profiles are typically symmetrically centered around the systemic velocity. Pahl et al. (2024) report that the ionising photon escape fraction, $f^{\mathrm{LyC}}_{\mathrm{esc}}$ correlates more strongly with both $\Delta_{\mathrm{peak}}$ and $V_{\mathrm{blue}}$ than with $V_{\mathrm{red}}$ at $z \sim 0.3$ . This suggests that the visibility of the blue peak in LAEs is closely linked to the escape of ionising photons, potentially indicating that similar neutral-phase gas dynamics govern both the escape of the blue Ly $\alpha$ peak and the ionising photons.

Figure 8. Velocities of Ly $\alpha$ blue and red peak are plotted against half of the peak separation ( $\Delta_{\mathrm{peak}}/2$ ) for LAEs with known systemic redshifts in our sample. The black dashed line shows the one-to-one relation.

We find a strong correlation between $\mathrm{FWHM}_{\mathrm{red}}$ and $\Delta_{\mathrm{peak}}$ (see Figure 9), indicating that broader red components are associated with larger peak separations. This trend is consistent with both observations (Yang et al. Reference Yang2017; Verhamme et al. Reference Verhamme2018; Kerutt et al. Reference Kerutt2022) and radiative transfer models (Verhamme et al. Reference Verhamme, Orlitová, Schaerer and Hayes2015). To carry out a linear regression, we merge our dataset with those of Kerutt et al. (Reference Kerutt2022), Yang et al. (Reference Yang2017), and Verhamme et al. (Reference Verhamme2018), and use LtsFit Footnote ^h (Cappellari et al. Reference Cappellari2013), which accounts for uncertainties in both variables and incorporates intrinsic scatter in the data. The combined data show Spearman correlation coefficients of p-value = $6.4\, \times 10^{-33}, \rho = 0.48$ . The intrinsic scatter estimated from the fit is $29.6 \, \mathrm{km}\,\mathrm{s}^{-1}$ . The resulting best-fit, shown by the solid black line in Figure 9, is given by

\begin{equation} \mathrm{FWHM}_{\mathrm{red}} = (0.336 \pm 0.017) \Delta_{\mathrm{peak}} + (69.1 \pm 8.6)\end{equation}

Figure 9. FWHM of the main (red) peak as a function of double-peak separation of Ly $\alpha$ . Red circles represent $z \sim 0.3$ Green Pea galaxies compiled by Yang et al. (Reference Yang2017). Yellow diamonds denote $z = 3$ –4 LAEs from Verhamme et al. (Reference Verhamme2018), using objects with peak separations listed in their Table 1, and pink triangles correspond to LAEs at $z \sim 2.9$ – $6.6$ from the MUSE-Wide and MUSE-Deep surveys (Kerutt et al. Reference Kerutt2022). The best-fit to the combined data is shown as a black solid line, with the corresponding $1\sigma$ confidence bounds indicated by blue dashed lines. The Spearman correlation coefficients and the equation of the best-fit line are displayed in the top left corner.

This result implies that the redshifted Ly $\alpha$ emission broadens as it escapes further from the systemic velocity. This behaviour has been interpreted via radiative transfer model using homogeneous shell models (Verhamme et al. Reference Verhamme2018), simulations considering anisotropic or bipolar outflows and inhomogeneous gas distributions (Zheng & Wallace Reference Zheng and Wallace2014), and multiphase clumpy media models (Gronke & Dijkstra Reference Gronke and Dijkstra2016). Together, these models predict a robust $V_{\mathrm{red}}$ – $\mathrm{FWHM}_{\mathrm{red}}$ correlation independent of outflow geometry or kinematics.

5.2.5. Trough flux density distribution

We find that several LAEs in our sample show positive residual flux in the absorption trough, $F_{\mathrm{trough}}$ , between the two peaks. The ratio of $F_{\mathrm{trough}}$ and the mean flux of the peak maxima, $\bar{F}_{\mathrm{peak}}$ , serves as a diagnostic to distinguish clumpy from homogeneous shell models, showing different distributions in each (Gronke & Dijkstra Reference Gronke and Dijkstra2016). Throughout the rest of this paper, we refer to this ratioFootnote ⁱ ( $F_{\mathrm{trough}}/ \bar{F}_{\mathrm{peak}}$ ) as the normalised trough flux density. Our result shows an anticorrelation between this ratio and peak separation (Figure 10, left panel), consistent with Orlitová et al. (Reference Orlitová2018) for the green pea galaxies at $z \sim 0.2$ – $0.3$ . Such a correlation is absent in the clumpy model (Gronke & Dijkstra Reference Gronke and Dijkstra2016). While the clumpy model predicts $F_{\mathrm{trough}}/ \bar{F}_{\mathrm{peak}}$ up to $\sim 0.8$ , our sample spans only between 0 and $0.4$ . We also find large scatter in peak separation for $F_{\mathrm{trough}}/ \bar{F}_{\mathrm{peak}} \sim 0 $ , and a decrease to $\lesssim$ $500\, \mathrm{km}\, \mathrm{s}^{-1}$ for $F_{\mathrm{trough}}/ \bar{F}_{\mathrm{peak}} \gtrsim 0.2 $ , supporting the predictions from the homogeneous shell model (Orlitová et al. Reference Orlitová2018).

Figure 10. Ly $\alpha$ double peak separation (left panel) and blue-to-total flux ratio (right panel) versus normalised trough flux density $F_{\mathrm{trough}}/ \bar{F}_{\mathrm{peak}}$ for the conservative sample. Red diamonds are data points from Orlitová et al. (Reference Orlitová2018). Only sources with $F_{\mathrm{trough}}/ \bar{F}_{\mathrm{peak}} \gt 0.009$ , corresponding to the lowest value reported in Orlitová et al. (Reference Orlitová2018), are presented. In the right panel, the horizontal black dashed line indicates the boundary between blue- and red-dominated regions.

In Figure 10 (right panel), we plot $F_{\mathrm{blue}}\, / F_{\mathrm{total}}$ against $F_{\mathrm{trough}}/ \bar{F}_{\mathrm{peak}}$ , revealing distinct distributions in blue- and red-dominated regions. In red-dominated regions, the trough flux rises with increasing blue peak flux. In contrast, strong blue-dominated LAEs ( $F_{\mathrm{blue}}\, / F_{\mathrm{total}}$ $\gt 0.8$ ) show low trough flux ( $F_{\mathrm{trough}}/ \bar{F}_{\mathrm{peak}} \lt 0.2$ ). This is likely due to the S/N of the Ly $\alpha$ line, given that the blue-to-total flux ratio shows a strong dependence on (S/N) $_{\mathrm{total}}$ , as discussed in Section 5.2.1.

5.3. Spatially resolved Ly $\alpha$ profiles

Integral field spectroscopy has enabled spatially resolved studies of LAEs. To investigate general trends and analyse the spectral properties of double-peaked LAEs as a function of radius, we employ a method similar to Erb et al. (Reference Erb2023). Among our 108 LAEs, we select 10 bright LAEs that show double-peaked profiles throughout their spatial extent with good S/N and exhibit clear spatial variations in their spectral profiles. They show extended Ly $\alpha$ halos in the continuum-subtracted narrowband images. While our selection ensures robust spatially resolved measurements, we note that double-peaked morphology does not necessarily persist at all radii in general. For example, Weldon et al. (Reference Weldon, Song, Reddy and Shapley2024) report a halo with a single peak in the central region but double peaks in the inner and outer halo.

Figure 11. Spatially resolved properties of one extended Ly $\alpha$ halo (ID: 1201250296) in our sample. 1st panel: Ly $\alpha$ pseudo narrow-band image; Coordinate is centred on the brightest pixel. Shown contours in white represent 2 and 4 $\sigma$ significance levels. The black circle denotes the location of stellar continuum detected in the MUSE white light image. The regions of spectral extraction using circular annular binning (as as described in Section 5.3) are shown as cyan circles. 2nd panel: Spectra from annular bins. 3rd and 4th panels: Pixel-by-pixel maps of the blue-to-total flux ratio ( $F_{\mathrm{blue}}\, / F_{\mathrm{total}}$ ) and peak separation, respectively, shown only for regions with $\geq 2 \, \sigma$ significance and zoomed to this area, and therefore not on the same spatial scale as the narrow-band image. These two maps are showing azimuthal variations in flux ratio and peak separation across the halo.

Annular spectra are created by binning spaxels radially. The high S/N of the Ly $\alpha$ halos enables reliable spatial binning. Starting from the central spaxel with the highest surface brightness, we include all spaxels with S/N $\gt 2$ in the continuum-subtracted Ly $\alpha$ narrowband images. Each halo is binned in radial increments corresponding to a single spaxel ( $0.2''$ ). The spectra within each radial bin are summed to construct Ly $\alpha$ profiles for different radii. We normalise these spectra by the total Ly $\alpha$ flux to clearly witness the variations in both peaks. We then fit a double-asymmetric Gaussian to these spectra by forward modeling, as in Section 4. In addition to analysing spatially averaged profiles, we also create pixel-by-pixel maps of the blue-to-total flux ratio and peak separation across the halos to investigate their azimuthal variations. Figure 11 shows an example of an extended Ly $\alpha$ halo in our sample, where the pseudo-NB image displays the extended emission. The spatially-resolved properties for all ten halos at $z = 3.2 - 3.7$ are shown in Figure A3 in the Appendix.

Figure 12 presents the spatial variations of the spectroscopic quantities as a function of radius in kpc. Halo radii span $4.6$ –12 kpc. In 9 of 10 halos, $F_{\mathrm{blue}} / F_{\mathrm{total}}$ increases with radius, consistent with $z \sim 2$ double-peaked Ly $\alpha$ studies (Erb et al. Reference Erb2023; Weldon et al. Reference Weldon, Song, Reddy and Shapley2024), likely due to line-of-sight effects in the outer halo showing more flux in the blue peak. In addition, peak separation decreases with radius, in agreement with the findings of Erb et al. (Reference Erb2023) and Weldon et al. (Reference Weldon, Song, Reddy and Shapley2024). However, we find that these spatial variations are not always circularly symmetric. Instead, significant azimuthal differences are present in some cases (see the example given in Figure 11, and also see Figure A3 in the Appendix), with certain regions within the halo exhibiting enhanced blue flux and reduced peak separation, rather than these features being uniformly distributed from the core to the outskirts. As a result, this azimuthal structure also somewhat influences the interpretations drawn from spatially averaged profiles based on annular regions. Azimuthal variations in the peak ratio suggest velocity asymmetries and nonradial gas motions at large radii. We also observe that the normalised trough flux density generally shows an increasing trend with radius, while Ly $\alpha$ red peak asymmetry shows a mild anticorrelation for some of them, but no strong trend is observed. Generally, beyond $\sim$ 4 kpc, the blue flux and the flux in the absorption trough rise, while the peak separation falls from center to outward, though the small sample limits physical constraints on this scale.

We also see some exceptions in these trends. Interestingly, one halo (ID: 1529235291) shows the opposite trend: a decreasing flux ratio with radius, while peak separation still decreases. Similar behaviour is seen in blue-dominated halos (Mukherjee et al. Reference Mukherjee2023; Bolda et al. Reference Bolda, Li, Erb, Steidel and Chen2024), where the central regions with the highest surface brightness are blue-dominated and the outskirts show a stronger red peak. We also notice a moderate decrease in normalised trough flux density with radius for this source. Additionally, the halo with ID: 1503243260 shows a slight increase in peak separation with radius. We explore these trends in more detail in the next section in the context of spatially resolved radiative transfer modeling.

Figure 12. Shown are results from line profile measurements of the spatially resolved spectra extracted from the annular regions described in Section 5.3. (top row) Left panel: blue-to-total flux ratio ( $F_{\mathrm{blue}}\, / F_{\mathrm{total}}$ ) as a function of radius (kpc), Right panel: Ly $\alpha$ peak separation versus radius; (bottom row) Left panel: normalised trough flux density ( $F_{\mathrm{trough}}/ \bar{F}_{\mathrm{peak}}$ ) versus radius, Right panel: Ly $\alpha$ red peak asymmetry ( $A_{f}$ ) versus radius. MAGPI IDs corresponding to each galaxy are indicated in the Right panel of bottom row.

6.1. Radiative transfer modeling

Studies on Ly $\alpha$ radiative transfer through clumpy shells have shown that the Ly $\alpha$ spectral profiles can be more complex than those from homogeneous shells. In this work, we adopt the simple and widely used homogeneous shell model. We use the publicly available radiative transfer code zELDA Footnote ^j (Gurung-López et al. Reference Gurung-López, Gronke, Saito, Bonoli and Orsi2022) to model both spatially integrated and resolved double-peaked Ly $\alpha$ profiles. zELDA is based on radiative transfer code FLaREON (Gurung-López, Orsi, & Bonoli Reference Gurung-López, Orsi and Bonoli2019). This code employs a thin-shell geometry for the gas distribution surrounding the Ly $\alpha$ emitting region, assuming an isothermal, homogeneous, spherical shell of neutral hydrogen at $10^4$ K, characterised by a uniform radial bulk velocity, $V_{\mathrm{shell}}$ . Importantly, zELDA accommodates both inflow ( $V_{\mathrm{shell}} \lt 0$ ) and outflow ( $V_{\mathrm{shell}} \gt 0$ ) scenarios, with emergent spectra from inflow is identical to those of a shell with an outflow velocity, but mirrored around the zero Ly $\alpha$ velocity (systemic). This framework enables an effective analysis of the gas dynamics within the ISM and CGM of the studied LAEs. The model incorporates five free parameters: shell velocity ( $V_{\mathrm{shell}}$ in km $\mathrm{s}^{-1}$ ), neutral hydrogen column density ( $\log(N_{HI})$ in cm $^{-2}$ ), dust optical depth ( $\tau$ ), intrinsic equivalent width ( $\log(EW)_{\mathrm{int}}$ in Å), and intrinsic line width ( $W_{\mathrm{int}}$ in Å). Ly $\alpha$ and continuum photons are generated from the central source with an intrinsic width of $W_{\mathrm{int}}$ and an intrinsic equivalent width of $\log(EW)_{\mathrm{int}}$ and propagate through the H i shell, undergoing absorption or resonant scattering while traveling.

zELDA performs forward modeling that accounts for the instrument’s resolution and pixelisation, and employs an MCMC approach to fit the observed line profiles. The S/N of the spectrum is also incorporated. To estimate uncertainties, the shell parameters and redshift are derived from the observed spectrum after applying noise perturbations 1 000 times. The best-fit models effectively reproduce key features of global double-peaked Ly $\alpha$ profiles. The best-fit H i column densities range from log $(N_{\mathrm{HI}}) = 18.36$ – $21.12 \, \mathrm{cm}^{-2}$ , which is consistent with the typical range observed in previous LAE studies.

6.2. Global spectroscopic properties and $N_{\mathrm{HI}}$ dependence

Radiative transfer modeling of Ly $\alpha$ emission, using both simplified shell geometries and more complex clumpy configurations, has successfully reproduced observed correlations between spectral features and physical model parameters. For example, a strong correlation has been found between the separation of the Ly $\alpha$ peaks and $N_{\mathrm{HI}}$ in the ISM and CGM of galaxies (see Verhamme et al. Reference Verhamme, Orlitová, Schaerer and Hayes2015; Verhamme et al. Reference Verhamme2018; Li & Gronke Reference Li and Gronke2022). A lower $N_{\mathrm{HI}}$ allows Ly $\alpha$ photons to escape with fewer scatterings, resulting in narrower lines and leading to a smaller peak separation. Conversely, at very high redshifts ( $z \gt 6$ ), LAEs within larger ionised bubbles exhibit broader red peaks (e.g. Mukherjee et al. Reference Mukherjee2024). Additionally, the removal of neutral gas at zero velocity due to early photoionisation facilitates the escape of Ly $\alpha$ photons on the blue side (Hayes et al. Reference Hayes, Runnholm, Scarlata, Gronke and Rivera-Thorsen2023). On the other hand, the strength of the blue peak tends to decrease with increasing shell outflow velocity (e.g. Verhamme et al. Reference Verhamme, Orlitová, Schaerer and Hayes2015; Li & Gronke Reference Li and Gronke2022; Erb et al. Reference Erb2023).

The IGM predominantly absorbs the blue peak, diminishing its line flux, while the red peak is shifted out of the resonance frequency. Stacked Ly $\alpha$ spectra reveal a decreasing fraction of blue-side flux with increasing redshift (Hayes et al. Reference Hayes, Runnholm, Gronke and Scarlata2021; Kramarenko et al. Reference Kramarenko2024), making double-peaked structures challenging to observe at high redshifts (Mason & Gronke Reference Mason and Gronke2020). Consequently, the number of double-peaked LAEs decreases rapidly with redshift (see Figure 2), as the increasingly neutral IGM attenuates the blue peak, leading to single-peak detections. Ionised bubbles around high-redshift LAEs may allow blue peaks with small separations to shift out of resonance before encountering the IGM. Therefore, we observe minimal redshift effects on $\mathrm{FWHM}_{\mathrm{blue}}$ . Kerutt et al. (Reference Kerutt2022) suggest that smaller peak separations correspond to smaller ionised bubbles necessary for the blue bump to shift sufficiently to survive IGM attenuation. More realistic radiative transfer models are required to further constrain these properties.

Figure 13. Normalised trough flux density $F_{\mathrm{trough}}/ \bar{F}_{\mathrm{peak}}$ versus global neutral hydrogen column density ( $N_{\mathrm{HI}}$ ) for sources with $F_{\mathrm{trough}}/ \bar{F}_{\mathrm{peak}} \gt 0.009$ (similar to Figure 10).

Figure 14. Examples of spatially resolved radiative transfer modeling of Ly $\alpha$ across four bright halos are shown. Observed spectra (black) and best-fit zELDA models (red) are presented for the ‘Core’ and ‘Outer halo’, with best-fit shell outflow velocities and H i column densities indicated in the top-left corner of each panel.

In Figure 13, the residual trough flux density is shown as a function of the global H i column density. In general, we find that LAEs with relatively lower column densities typically show a higher trough flux. The higher the column density, the more scattering and absorption Ly $\alpha$ photons experience near the line center. This pushes Ly $\alpha$ photons further into the red and blue wings before they can escape. This results in greater suppression of the central flux, causing the trough between the peaks to deepen as $N_{\mathrm{HI}}$ increases. This trend is in agreement with the well-known correlation between peak separation and $N_{\mathrm{HI}}$ , given the anticorrelation between peak separation and $F_{\mathrm{trough}}/ \bar{F}_{\mathrm{peak}}$ in Figure 10. Therefore, a smaller peak separation and a higher trough flux strongly suggest the presence of low- $N_{\mathrm{HI}}$ channels, which we will explore further in Section 6.4.

6.3. Spatially resolved modeling and variations across halos

To characterise the gas properties across the Ly $\alpha$ halos, we perform spatially resolved spectral modeling using zELDA. We fit spectra from the central ‘core’ region (within a $0.2''$ radius aperture) and the ‘outer halo’ region, selecting either the outermost annular bin or the preceding bin based on spectral quality. The modeling results for four bright halos are presented in Figure 14. zELDA successfully reproduces the observed double-peaked profiles in both the core and outer halo regions. The presence of a double-peaked profile in the cores indicates that Ly $\alpha$ photons experience multiple scatterings within the neutral gas shell. This contrasts with the halo studied by Weldon et al. (Reference Weldon, Song, Reddy and Shapley2024), where the central region appears single-peaked. They attribute this to backscattering, implying strong outflows that allow photons to escape without undergoing additional scatterings. In our case, the persistent double-peaked profile suggests that Ly $\alpha$ photons are more likely to undergo multiple scatterings before escaping.

For all ten halos, best-fit shell velocities and $N_{\mathrm{HI}}$ values for the core and outer halo regions are shown in Figure 15, revealing a clear decline in both parameters from the core to the outskirts. These findings align with previous studies (e.g. Erb et al. Reference Erb2023; Weldon et al. Reference Weldon, Song, Reddy and Shapley2024), supporting the notion that neutral gas velocity and $N_{\mathrm{HI}}$ regulate Ly $\alpha$ spectral variations across halos. The Ly $\alpha$ blue-to-total flux ratio is influenced by the line-of-sight velocity, which is largely a geometric effect arising from the decreasing projected component of the radial outflow velocity with increasing radius (Li et al. Reference Li, Steidel, Gronke, Chen and Matsuda2022). Although the decrease in shell expansion velocity in the outermost regions of the Ly $\alpha$ halo can be attributed to the gravitational deceleration of galactic winds at large radii (Weldon et al. Reference Weldon, Song, Reddy and Shapley2024), consistent with several observational and simulation studies (Barai et al. Reference Barai2013; Thompson et al. Reference Thompson, Quataert, Zhang and Weinberg2016; Chen et al. Reference Chen2020), the broader theoretical picture is more nuanced. Recent high-resolution simulations (e.g. Schneider et al. Reference Schneider, Ostriker, Robertson and Thompson2020) demonstrate that momentum transfer from the hot to the cool phase can efficiently accelerate cool gas out to kiloparsec scales. Additionally, the enhanced blue peak observed in the outer regions may also result from the preferential scattering of blueshifted Ly $\alpha$ photons to larger radii (Zheng et al. Reference Zheng, Cen, Trac and Miralda-Escudé2010). Peak separation trends are primarily driven by line-of-sight variations in $N_{\mathrm{HI}}$ , with lower column densities at larger radii allowing photons to escape closer to the line center without suffering much scattering, reducing overall peak separation. Additionally, the observed trend of increasing residual flux in the absorption troughs with radius is mainly governed by the amount of residual neutral hydrogen in the ISM. A high $N_{\mathrm{HI}}$ indicates a larger amount of H i along the line of sight, which tend to absorb more Ly $\alpha$ photons near the line center, leading to a deeper trough between the peaks.

Notably, MAGPI1529235291 shows an increase in shell outflow velocity from the core to outer halo (panel d of Figures 14 and 15), alongside a decreasing flux ratio with radius (Figure 12). Although its global and central Ly $\alpha$ spectra are red-dominated, the radial trend resembles that of a blue-dominated (accretion-dominated) system, suggesting active gas accretion in the central region and outflows in the outskirts. This LAE, with an integrated blue-to-total flux ratio of $0.42$ , exhibits a higher Ly $\alpha$ luminosity ( $6.92 \times 10^{42}\, \mathrm{erg}\, \mathrm{s}^{-1}$ ) compared to others within the same ratio bin (see Figure 6), highlighting its distinct behaviour relative to typical outflow-dominated systems. In the future, spatially resolved studies of sources with comparable blue-to-total flux ratios could constrain the turnover point, beyond which the opposite trend is expected.

Figure 15. Results from the spatially resolved modeling showing the shell velocity ( $V_{\mathrm{shell}}$ ) and H i column density for two distinct regions: the core and the outer halo. MAGPI IDs are indicated in the legend, following the same order as in Figure 12.

6.4. Connection to global LyC escape

Previous studies suggest that significant Ly $\alpha$ flux near systemic velocity may be a potential indicator of LyC escape (e.g. Rivera-Thorsen et al. Reference Rivera-Thorsen2017; Naidu et al. Reference Naidu2022). Given that many LAEs in our sample display positive residual flux at the absorption trough or close to the systemic velocity, we revisit the spatially integrated Ly $\alpha$ profiles to assess their potential as LyC leaker candidates. In the local universe, double peak separation correlates with LyC escape, with narrower peak separations generally associated with higher escape fractions (Verhamme et al. Reference Verhamme2017; Izotov et al. Reference Izotov2018). Peak separation is a well-established parameter in the context of Ly $\alpha$ radiative transfer, and is directly linked to $N_{\mathrm{HI}}$ (Verhamme et al. Reference Verhamme, Orlitová, Schaerer and Hayes2015; Verhamme et al. Reference Verhamme2017; Li & Gronke Reference Li and Gronke2022). However, peak separation alone does not uniquely reveal Ly $\alpha$ escape routes (Kakiichi & Gronke Reference Kakiichi and Gronke2021), as narrow separations can result from multiple small openings in a turbulent medium or a single large, wind-driven cavity (Hu et al. Reference Hu2023). Instead, it can distinguish between different LyC leakage paths. Low LyC leakage, where only a few pathways allow LyC to escape, scattered Ly $\alpha$ photons traverse optically thick channels with high $N_{\mathrm{HI}}$ , resulting in a broad peak separation. In contrast, high leakage corresponds to narrow separations in density-bounded channels. However, at high redshifts ( $z \gt 2$ ), no clear correlation is observed between peak separation and LyC escape fraction (Kerutt et al. Reference Kerutt2024).

We present the red peak asymmetry parameter, $A_{f}$ , originally introduced by Rhoads et al. (Reference Rhoads2003). This parameter has recently been used to quantify the multiphase nature of turbulent H ii regions and their connection to LyC escape (Kakiichi & Gronke Reference Kakiichi and Gronke2021; Hu et al. Reference Hu2023).

Figure 16 illustrates the relationship between Ly $\alpha$ peak separation and red peak asymmetry, revealing an anticorrelation: profiles with strong asymmetry exhibit small peak separations, while a wider range of separations is observed for profiles with low asymmetry. When Ly $\alpha$ photons have direct escape routes through ionised, low-density channels, they experience fewer scatterings, leading to more symmetric profiles. Similarly, if the ISM is largely neutral but has a few clear, direct paths, the result is also relatively symmetric. However, in media with both density-bounded (low $N_{\mathrm{HI}}$ ) and optically thick neutral regions, some photons escape directly, while others scatter extensively within the dense, neutral gas, resulting in asymmetric profiles. In such a scenario, the scattered photons take longer paths, potentially emerging at different velocities than those escaping directly.

Figure 16. Ly $\alpha$ double peak separation versus Ly $\alpha$ red peak asymmetry. The grey-, blue-, and pink-shaded regions correspond to distinct LyC leakage regimes: minimal leakage (grey), escape through low $N_{\mathrm{HI}}$ channels (pink), and escape through large, density-bounded holes (blue). The black dashed line corresponds to the empirical boundary set by Hu et al. (Reference Hu2023). The brown squares mark 5 LAEs in our sample with $F_{\mathrm{trough}}/ \bar{F}_{\mathrm{peak}} \gtrsim 0.2$ , indicating they may be potential LyC leakers.

Following Hu et al. (Reference Hu2023), we classify the diagram (Figure 16) into three regions: (grey) minimal LyC leakage corresponding to systems with no or very few escape paths; (pink) significant leakage through ionisation-bounded channels with low column density ( $f^{\mathrm{LyC}}_{\mathrm{esc}}$ $\gt 10\%$ ); and (blue) substantial leakage through large holes or density-bounded channels ( $f^{\mathrm{LyC}}_{\mathrm{esc}}$ $\gt 10\%$ ). Based on this classification, we identify five LAEs as potential LyC leakers: four in the blue region ( $A_{f} \lt 3$ ), suggesting escape through density-bounded channels, and one in the pink region ( $A_{f} \gt 3$ ), indicating ionisation-bounded leakage in a multiphase medium. Using the empirical relation between peak separation and $f^{\mathrm{LyC}}_{\mathrm{esc}}$ from Izotov et al. (Reference Izotov2018), these five LAEs have $f^{\mathrm{LyC}}_{\mathrm{esc}}$ $\sim 10$ –27%. Many LAEs reside in the grey region, suggesting lower LyC escape fractions compared to those five LAEs in the blue and pink zones. Additionally, as discussed in Section 5.2.5, an anticorrelation exists between the peak separation and $F_{\mathrm{trough}}/ \bar{F}_{\mathrm{peak}}$ . Three LAEs in the blue zone exhibit higher net Ly $\alpha$ flux in the trough (i.e. $F_{\mathrm{trough}}/\bar{F}_{\mathrm{peak}} \gt 0.2$ ), supporting the conclusion that these galaxies have LyC-thin sightlines (Verhamme et al. Reference Verhamme, Orlitová, Schaerer and Hayes2015; Gazagnes et al. Reference Gazagnes, Chisholm, Schaerer, Verhamme and Izotov2020). Recent studies show that some high-redshift strong leakers instead exhibit large peak separations ( $\Delta_{\mathrm{peak}} \gt 350\, \mathrm{km} \,\mathrm{s}^{-1}$ ; Fletcher et al. Reference Fletcher2019; Naidu et al. Reference Naidu2022; Kerutt et al. Reference Kerutt2024). Hu et al. (Reference Hu2023) has established an empirical boundary between the blue and grey zones, located near a peak separation of 400 km s $^{-1}$ . Applying this boundary, along with the trough flux criterion of $F_{\mathrm{trough}}/\bar{F}_{\mathrm{peak}} \gt 0.2$ , we identify an additional two LAEs, resulting in a total of five LAEs as possible LyC leakers.

Cosmological hydrodynamic simulations highlight that the central escape fraction (fraction of the flux near the systemic velocity or absorption trough) is a more robust diagnostic of LyC leakage, whereas peak separation or asymmetry separately are insufficient indicators (Choustikov et al. Reference Choustikov2024). Given the diagnostic potential of double-peaked Ly $\alpha$ profiles in constraining CGM kinematics, column density, and LyC escape, expanding double-peaked sample will provide a more statistically significant dataset to further explore the physical mechanisms driving Ly $\alpha$ and LyC radiative transfer.

6.5. Caveats in modeling and analysis

Finally, we note several caveats, particularly concerning the modeling of Ly $\alpha$ profiles and the interpretation of the derived results. Our analysis employs a spherically symmetric radiative transfer model applied to annularly averaged spectra, assuming that all Ly $\alpha$ photons originate at the galaxy center and then scatter through an expanding or infalling shell before escaping into the CGM. This geometry is highly simplified, and the assumptions of central photon production and spherical symmetry are not fully consistent with the setup of annular averaging. In addition, the best-fit models are mainly driven by the high-S/N spectra from the innermost bins, while the outermost bins, though critical for tracing gradients in spectroscopic properties, provide weaker constraints on model parameters. Extracting physical insights from different halo regions will require more realistic anisotropic radiative transfer models with non-uniform H i distributions. Thus, the results presented here should be regarded as average parameters representative of the modeled regions rather than as precise local measurements.