© The Authors 2025

1. Introduction

Studying the early stages of galaxy formation is fundamental for understanding their evolution and diversity in general. This requires observations of the most distant galaxies. Although the field has recently been revolutionised with the launch of the James Webb Space Telescope (JWST), galaxies in these early phases still appear as very faint and close-to-unresolved sources. While JWST has provided rest-frame optical emission line detections out to redshifts z ∼ 7 − 11 (e.g. Bunker et al. 2023; Arrabal Haro et al. 2023; Heintz et al. 2024; Saxena et al. 2023; Tang et al. 2023; Boyett et al. 2024), the rest-frame ultraviolett Lyman-α line (Lyα; λ_Lyα = 1215.68 Å) remains of interest for galaxies in the early universe (see reviews by Hayes 2015; Dijkstra 2019; Ouchi et al. 2020). Originating mainly from the recombining gas in the vicinity of hot O- and B-stars, Lyα is associated in particular with young galaxies that are in the process of rapidly converting gas into stars, as famously envisioned by Partridge & Peebles (1967).

An observational census of Lyα-emitting galaxies (LAEs) has been established for galaxies with Lyα luminosities of L_Lyα ≳ 10^41.5 erg s⁻¹ throughout most of the history of the Universe (e.g., Ouchi et al. 2008; Blanc et al. 2011; Sobral et al. 2018; Wold et al. 2017; Herenz et al. 2019; Liu & Jiang 2023; Thai et al. 2023) up until the epoch of reionisation (EoR; e.g., Ouchi et al. 2010; Matthee et al. 2015; Santos et al. 2016; Konno et al. 2018; Ning et al. 2022; Wold et al. 2022) by means of the LAE luminosity function. These censuses ascertain that LAEs are an abundant population of high-redshift (z ≳ 2) galaxies that can serve as a tracer for addressing fundamental questions of galaxy formation and cosmology. For example, the HETDEX survey will measure redshifts of more than a million LAEs from z ∼ 1.8 to z ∼ 2.5 to investigate the nature of dark energy (Gebhardt et al. 2021; Mentuch Cooper et al. 2023). Moreover, LAEs also serve as an important probe to constrain the timeline and the topology of the intergalactic medium during the EoR (e.g. reviews by Dijkstra 2014 and Choudhury 2022; see also recent results by Nakane et al. 2024, Witstok et al. 2024, and Witten et al. 2024). Finally, large-area cosmological line-intensity mapping experiments that attempt to measure baryonic acoustic oscillations as a function of the age of the Universe will also rely on the Lyα signal from the unresolved population of LAEs (e.g. Bernal & Kovetz 2022; Lujan Niemeyer et al. 2023).

Given the high Lyα absorption cross-section of neutral hydrogen, the resulting Lyα optical depths in typical interstellar media (ISM) are extremely high (τ_Lyα > 10⁶ for neutral columns of N_HI > 10²⁰ cm⁻²). Hence the Lyα line is resonant. The resulting complex radiative transfer within the ISM redistributes the initial Lyα radiation field spatially and in frequency space. Moreover, the radiative transfer increases a Lyα path length of the photon and makes it prone to absorption by dust. Theoretical works demonstrate the resulting intricate relationships between Lyα observables (e.g. Lyα spectral morphology, Lyα surface brightness distribution, and total Lyα line intensity) and the density distribution, dust distribution, and kinematics of the interstellar scattering medium (see, e.g. recent papers by Remolina-Gutiérrez & Forero-Romero 2019; Smith et al. 2019; Lao & Smith 2020; Song et al. 2020; Gronke et al. 2021; Li et al. 2022; Blaizot et al. 2023; Munirov & Kaurov 2023, and references therein).

A key quantity in the context of LAEs is the Lyα escape fraction, $f_{\text{esc}}^{\text{Ly}\alpha }$, which is defined as the ratio between the observed Lyα luminosity and the intrinsic Lyα luminosity. An understanding of the galactic characteristics that regulate $f_{\text{esc}}^{\text{Ly}\alpha }$ is required to answer the question of what makes a star-forming galaxy an LAE. Many observational (e.g., Giavalisco et al. 1996; Keel 2005; Finkelstein et al. 2009; Hayes et al. 2011; Malhotra et al. 2012; Atek et al. 2014; Hathi et al. 2016; Hagen et al. 2016; Oyarzún et al. 2016; Paulino-Afonso et al. 2018; Matthee et al. 2021; McCarron et al.2022; Chávez Ortiz et al. 2023; Napolitano et al. 2023; Hayes et al. 2023; Lin et al. 2024) and modelling (e.g. Schaerer & Verhamme 2008; Verhamme et al. 2008; Forero-Romero et al. 2011; Dijkstra & Wyithe 2012; Dayal & Libeskind 2012; Garel et al. 2015, 2016; Gurung-López et al. 2019) efforts have been conducted with the intention of answering this question. It is now established that the physical and morphological characteristics of LAEs are significantly different from star-forming galaxies where $f_{\text{esc}}^{\text{Ly}\alpha }$ ∼ 0. As reviewed by Ouchi et al. (2020), high-z LAEs are commonly compact (r_e ∼ 1 kpc) star-forming galaxies (SFR ∼ 1…10 M_⊙yr⁻¹) with low-dust content (E(B − V)≲0.2) that harbour young stellar populations (age ∼ 10 Myr) of low stellar mass (∼10⁸…10⁹ M_⊙). Moreover, the ISM of those galaxies is often metal poor, with metallicities ranging from Z ∼ 0.1 Z_⊙ to Z ∼ 0.3 Z_⊙. This differentiates LAEs from non-LAEs and/or continuum-selected high-z galaxies, which are typically more evolved and more massive systems, but there is certainly overlap between the two populations.

The defining criterion for LAEs in surveys is often based on a threshold in the equivalent width of the Lyα line, EW_Lyα. Interestingly, EW_Lyα is linearly correlated with $f_{\text{esc}}^{\text{Ly}\alpha }$ (Sobral & Matthee 2019; Melinder et al. 2023). Sobral & Matthee (2019) explain the measured slope of this relation in a simple scenario where, on average, galaxies with higher EW_Lyα and $f_{\text{esc}}^{\text{Ly}\alpha }$ are dominated by stellar populations with a harder ionising spectrum, as characterised by the ionising photon production efficiency ξ_ion, and lower E(B − V). This idea is consistent with results that show how the overall dust-poor high-EW_Lyα LAEs are indeed characterised by higher ξ_ion (e.g., Trainor et al. 2016; Nakajima et al. 2016, 2018; Maseda et al. 2020; Hayes et al. 2023; Kramarenko et al. 2024). Certainly, ξ_ion and E(B − V) are galaxy characteristics that we expect to directly influence the radiative transfer. The fact that these characteristics occur predominantly in compact star-forming systems of low metallicity and low mass informs us more about the processes of star formation and, more generally, galaxy formation physics in such systems. However, the causal connection between Lyα radiative transport physics and stellar mass or star-formation rate appears to involve some other process. In this respect we reiterate that galaxy and ISM kinematics are causally connected with star formation and stellar mass (see e.g. review by Glazebrook 2013). Thus we may also expect also causal relationships between galaxy kinematics and/or $f_{\text{esc}}^{\text{Ly}\alpha }$ or EW_Lyα.

Significant effort has been devoted to understanding the effect of kinematics and gas distribution on the Lyα line profile morphology (e.g., Blaizot et al. 2023, and references therein). The kinematics of the scattering medium are also frequently studied in absorption, especially using low-ionisation-state absorption lines of Si⁺ and C⁺ ions in the UV (Kunth et al. 1998; Wofford et al. 2013; Henry et al. 2015; Rivera-Thorsen et al. 2015; Hayes 2023; Hayes et al. 2023). These line-profile analyses and absorption-line studies reveal how low-neutral-H I covering fractions and outflow kinematics appear to promote Lyα escape. In principle, the neutral phase can also be probed by the 21 cm line. However, the spatial resolution of the radio beam is usually significantly lower then the spatial scales that can be probed by imaging and spectroscopy (Cannon et al. 2004; Pardy et al. 2014). Nevertheless, such analyses reveal that mergers offset large amounts of H I via tidal interactions and it appears likely that the thereby reduced neutral column towards the star-forming regions promotes Lyα escape (Purkayastha et al. 2022; Le Reste et al. 2024). A promising way forward comes in the form of highly spatially resolved aperture-synthesis observations in combination with other observational probes of the emitting and scattering medium. Especially recent work by Le Reste et al. (2022) informs us about the role of the interplay between H II and H I gas on ISM scales in facilitating Lyα escape.

To date, few studies have sought to understand the role of gas kinematics in shaping the global Lyα characteristics $f_{\text{esc}}^{\text{Ly}\alpha }$, EW_Lyα, and Lyα luminosity (L_Lyα). To this aim we must leverage galaxy-averaged kinematical statistics. Importantly, this averaging should be done on spatial scales most relevant for the radiative-transfer problem. Such an analysis was attempted by Herenz et al. (2016), hereafter H16, who used integral-field spectroscopic observations of the Hα line of the Lyman-α reference sample (LARS; Östlin et al. 2014). This study found that, galaxies whose ionised gas kinematics are dominated by unordered (turbulent) motions can appear as Lyα emitters, whereas galaxies that are dominated by large-scale coherent motions never appear to show significant Lyα escape. Evidence for turbulence-driven Lyα escape was also presented by Puschnig et al. (2020) from a study of the dust and the molecular gas content of LARS. Moreover, a connection between global galaxy kinematics and Lyα escape was reported very recently by Foran et al. (2024) for the first time in high-z galaxies. These authors analysed the global kinematical characteristics of z ∼ 2 and z ∼ 3 galaxies for which consistent EW_Lyα measurements were available. Objects that showed no significant Lyα emission or even Lyα in absorption (i.e. low or negative EW_Lyα) were found to be a mixture of dispersion- and rotation-dominated systems, whereas galaxies with significant Lyα emission were always found to be dispersion dominated.

The goal of the present study is to revisit the relationship between galaxy kinematics and Lyα characteristics observationally, with the main aim being to remove the shortcomings of the rather small sample used in the kinematical study of H16. To this aim, we analysed integral-field spectroscopic observations of the Hα line for all 42 galaxies of the LARS + extended LARS programmes (eLARS; Melinder et al. 2023, herafter M23). The enlarged dataset now allows us to discuss the importance of galaxy kinematics for Lyα observables with respect to other galaxy characteristics in a quantitative manner.

This paper is structured as follows: Section 2 describes the sample and data, especially the newly obtained Calar-Alto 3.5m PMAS observations. Section 3 presents the kinematic analysis and our relations between galaxy kinematics and Lyα observables. Section 4 then discusses the importance of kinematic parameters with respect to the “soup of quantities” that have been found to correlate with Lyα properties of galaxies. An astrophysical interpretation of the results is also presented in Sect. 4. Finally, Sect. 5 concludes this paper.

2. Data

2.1. The (Extended) Lyman-α Reference Sample

The sample analysed here consists of 42 star-forming galaxies selected from the SDSS spectroscopic database. The main selection criteria are cuts in the Hα equivalent width, EW_Hα (as a proxy for stellar age), and FUV luminosity (at λ ≈ 1500 Å from GALEX), L_UV (as a proxy for star-formation rate). Fourteen of the 42 galaxies had to satisfy a cut in EW_Hα > 100 Å. These 14 galaxies (LARS1 – LARS14) comprised the initial Lyα Reference Sample (LARS; Östlin et al. 2014; Hayes et al. 2014). As detailed in Östlin et al. (2014), additional considerations for the sample selection were that the galaxies cover a range in FUV luminosities typical of star-forming galaxies in the early universe, and the availability of relevant Hubble Space Telescope (HST) archival data to design the HST imaging and spectroscopic campaign as economically as possible.

The initial sample was extended with 28 star-forming galaxies (eLARS 1–28) that had to satisfy a cut of EW_Hα > 40 Å, and again the sample was picked to include a dynamical range in L_UV (M23). Only 7 out of those 28 galaxies show EW_Hα > 100 Å. The final sample, LARS+eLARS, now spans a range in log₁₀(L_FUV [L_⊙]) from 9 to 11 and in EW_Hα from 41 Å to 578 Å. Stellar mass and Hα star-formation rate distributions are presented in Sect. 3.5. All galaxies show emission line ratios consistent with being powered by a young stellar population according to the extreme starburst classification demarcation line in the [O III]λ5007]/Hβ vs. [N II]λ6548/Hα diagram provided by Kewley et al. (2001). However, five galaxies could be classified as harbouring an AGN according to the criterion of Kauffmann et al. (2003). The HST imaging observations of all 42 galaxies were used to synthesise images in Lyα, Hα, and Hβ. We use for the present analysis global measurements derived from these images (see Sect. 2.4 below).

2.2. PMAS observations of eLARS

We observed all 28 eLARS galaxies with the Potsdam Multi Aperture Spectrophotometer (PMAS; Roth et al. 2005) at the Calar Alto 3.5 m telescope. An overview of the campaign is provided Table A.1 of Appendix A. We used PMAS in its Lens Array configuration with double magnification – in this configuration the 256 spectral pixels (spaxels) of 1″×1″ size sample contiguously a 16″×16″ field of view. Within this field of view the bulk of the high surface-brightness Hα emission of the galaxies in the sample could be covered with a single pointing. Exceptions are eLARS 3, 5, and 26, which required two pointings (pointings are differentiated by the suffixes A and B in Table A.1). We used the R1200 backward-blazed grating as in H16, but different to H16, here the grating angle could be held fixed so to disperse the first spectral order around Hα close to the centre of the detector for all galaxies. According to the PMAS grating table¹ the expected nominal resolving power for the adopted position (GROTPOS = 109) is R ∼ 5000. In order to sample the spectroscopic line spread function adequately at this resolving power (see Robertson 2017) we did not bin the read out of the PMAS 4k×4k CCD along the dispersion direction (XBIN = 1). This set up was chosen to obtain data that is optimal for a kinematic study of the ionised gas.

The observations took place during the course of two visitor mode observing runs in spring 2016 (run350) and spring 2017 (run362). We observed all our targets at air-masses ≲1.3. The observing conditions were not always optimal, with mediocre seeing and occasional thick cloud cover. Nevertheless, by visually monitoring the incoming read outs of the PMAS 4k×4k CCD after each exposure, we ensured that there is sufficient signal in the Hα lines that could be used for our analysis. Observations that did not contain enough signal due to cloud cover were repeated. If the conditions permitted, we took an additional 400s blank sky exposure close to the target. We ensured that such blank sky exposures were always taken when the redshifted Hα line was in spectral proximity to a bright telluric line. Sky flat exposures were taken at the beginning or end of each night. On-sky exposures were always flanked by a continuum lamp exposure and an arc-lamp exposure (HgNe).

Due to strong ambient temperature fluctuations the spectrograph was not always optimally focused. Thus, especially during the 2016 run, sometimes lower values of R are measured (see Table A.1 – the resolving power determination and correction is described in Appendix B). Fortunately, the galaxy’s Hα profiles are always resolved significantly, also in the slightly suboptimal R ∼ 3000 datasets.

We reduced the observational raw data with the general reduction pipeline for fibre-fed integral-field spectrographs p3d (Sandin et al. 2010, 2012). Our reduction strategy is detailed in Sect. 3 of H16, and for completeness we provide a summary in Appendix A.

2.3. The original LARS PMAS data

In this study, we combine our new observational results with the results from our study of Hα kinematics in the original LARS sample H16. The LARS datacubes are available via the CDS². In the following we reanalyse these LARS datacubes with the methods used for the eLARS observations presented in Sect. 3.1 and described in detail in Appendix B. While there are slight quantitative differences in the recovered velocity fields and velocity dispersion maps (see Sect. 3.3), their visual appearance is not significantly different from what was presented in H16. Thus, we streamline the presentation here by not presenting those maps again; the interested reader is referred to Figs. 5–7 in H16.

2.4. Global measurements from HST imaging

For relating the ionised gas kinematics with the galaxies physical parameters (stellar mass M_⋆, star-formation rate SFR_Hα, starburst age) and their Lyα observables (equivalent width EW_Lyα, escape fraction $f_{\text{esc}}^{\text{Ly}\alpha }$, and observed Lyα luminosity, L_Lyα) we rely on the integrated measurements provided in M23. As described there in detail, the required global quantities were measured within circular apertures on Voronoi-binned maps that were created from UV and optical PSF-matched HST imaging data that was processed with an updated version LaXs pipeline (Hayes et al. 2009). The radii of these apertures, denoted r_tot in M23, are defined on the Lyα maps, as Lyα is often more extended then Hα or the FUV continuum. The EW_Lyα is defined as the ratio between the aperture integrated Lyα luminosity and the UV continuum luminosity at λ_Lyα, and $f_{\text{esc}}^{\text{Ly}\alpha }$ is defined as

$$\begin{array}{c}\hfill f_{\mathrm{esc}}^{\mathrm{Ly}\alpha }=\frac{L_{\mathrm{Ly}\alpha }/L_{\mathrm{H}}\alpha }{{(L_{\mathrm{Ly}\alpha }/L_{\mathrm{H}}\alpha )}_{\mathrm{int}}·{10}^{0.4·k({\lambda }_{\mathrm{H}}\alpha )·E{(B-V)}_n}},\end{array}$$(1)

where k(λ_Hα) is the extinction coefficient at the wavelength of Hα, and E(B − V) is the aperture integrated colour excess that attenuates nebular emission and is derived from the Balmer-decrement,

$$\begin{array}{c}\hfill E(B-V)=\frac{2.5}{k({\lambda }_{\mathrm{H}\beta })-k({\lambda }_{\mathrm{H}}\alpha )}·{log}_{10}\frac{L_{\mathrm{H}}\alpha /L_{\mathrm{H}}\beta }{{(L_{\mathrm{H}}\alpha /L_{\mathrm{H}}\beta )}_{\mathrm{int}}},\end{array}$$(2)

using the integrated luminosities derived from the reconstructed Hα and Hβ images. The intrinsic ratios of Lyα / Hα and Hα/Hβ luminosities are given by the canonical Case-B recombination values (L_Lyα/L_Hα)_int = 8.7 (see, e.g. Fig. 4.2 in Hayes 2015) and (L_Hα/L_Hβ)_int = 2.86, respectively. Moreover, M23 use the Cardelli et al. (1989) extinction curve for k(λ)³ in Eq. (1) and Eq. (2). Nebular extinction is thus assumed to be caused by a homogeneous dust screen in front of the line emitting gas.

We note that no significant Lyα emission could be detected for six galaxies in the sample (LARS 6, LARS 10, LARS13, ELARS 12, ELARS 14, and ELARS 16). Hence, these galaxies only have upper limits in their Lyα observables. Furthermore, we note that M23 performed a refined and more accurate reprocessing of the HST imaging data also for the 14 galaxies constituting the original LARS sample. Hence, the global Lyα luminosities, equivalent widths, and escape fractions are different and supersede the values presented in Hayes et al. (2014) and Runnholm et al. (2020). Adopting the a cut of EW_Lyα ≥ 20 Å to define a galaxy as LAE, the sample contains 12 LAEs.

3. Analysis

3.1. Line-of-sight velocity and velocity dispersion maps

We project the kinematical properties traced by Hα in the IFS cuboids into 2D maps of line-of-sight velocity v_los and velocity dispersion σ_v. Our procedures for creating the S/N, v_los, and σ_v maps closely follow the well established convention for this task that involves fitting single-component Gaussian profiles to the continuum subtracted data (e.g. Alonso-Herrero et al. 2009). We describe the creation of the maps in more detail in Appendix B.

We here note that only for the galaxies eLARS 7 and eLARS 24 numerous spaxels are not well described by a single Gaussian component. The spectral profile in those spaxels show strong bimodality indicating the projection of two distinct kinematic components along the line of sight⁴. Thus the width and position of the Gaussian fits becomes unreliable parameters. A similar situation is also encountered in LARS 9 and LARS 13 from the H16 sample. All four galaxies with double component profiles are excluded in the following statistical analyses involving velocity dispersion measurements.

We show the v_los and σ_v maps alongside their signal-to-noise maps (S/N) in Hα for four eLARS galaxies in Fig. 1; Fig. B.1 in Appendix B provides those maps for the full sample. An atlas that allows for a visual comparison between the HST imaging data with the kinematical maps is presented in Appendix A.1 of Schaible (2023). To the best of our knowledge, kinematic maps of the ionised phase of the galaxies analysed in this study have not been presented elsewhere, with the exception being eLARS 1, for which Hα kinematics from scanning Fabry-Perot observations were analysed in Sardaneta et al. (2020).

Fig. 1. Hα kinematics of four eLARS galaxies from PMAS (16″ × 16″ field of view). For each galaxy we show HαS/N (left panels; colour coded in log-scale from 1 to 1000), vlos [km s−1] (centre panels; velocities are colour coded linearly from blue (approaching) to red (receding) with the centre at vlos = 0 in white and scaled symmetrically to the absolute maximum), and σv [km s−1] (right panels; colour coded linearly from the 2nd percentile to the 98th percentile of the observed σv distribution per galaxy). The name of each galaxy is coloured according to the visual kinematic classification (Sect. 3.2): rotating discs in purple, perturbed rotators in green, and galaxies with complex kinematics in orange. The full sample is shown in Appendix B (Fig. B.1).

3.2. Visual classification of the velocity fields

3.2.1. Classification scheme and results

We visually characterise the v_los maps shown in Fig. B.1 according to three qualitative kinematical categories (see review by Glazebrook 2013): “rotating disc” (RD), “perturbed rotator” (PR), and “complex kinematics” (CK).

For RD galaxies the bulk of the Hα emitting gas needs to be dominated by orbital motion. The resulting velocity field shows a smooth gradient which steepens in the centre and this gradient is aligned parallel to the major axis of the optical image. Projection of a galaxy’s differential rotation lead to a classical symmetric “spider-leg” morphology of observed iso-velocity contours (not shown here; see, e.g. Fig. 8 in Gnerucci et al. 2011 or Fig. 6 in Epinat et al. 2010 for examples). Importantly, the receding and approaching velocities mirror each other. Moreover, for RD galaxies the velocity dispersion map often show elevated σ_v values near the kinematical centre (the reasons for this are provided in Appendix C). The velocity fields of PR galaxies are still dominated by orbital motions, but they show significant deviations from a perfect disc. Lastly, galaxies characterised by CK exhibit no clear gradient indicative of simple orbital motions. For CK galaxies the gradient also appears less steep in comparison to RDs or PRs.

This visual classification is certainly subjective as there is a not a well defined boundary between the classes. Hence, the classification was done by six authors individually (ECH, AS, PL, MH, ALR, and GÖ) and then finalised in a consolidation session. Moreover, as described in detail in Appendix C.1 the classification of the RDs is confirmed by successfully modelling the RD systems by simple disc models with GalPak^3D (Bouché et al. 2015) and by the symmetric appearance of the line-of-sight velocity profiles along the kinematical major axis (Appendix C.2; Fig. C.3). The classifications for the combined LARS and eLARS sample are listed given in Appendix D (Table D.1). Moreover, we colour the label of each galaxy in Fig. B.1 according to its kinematical class.

According to our visual classification, the eLARS sample contains 4 CK galaxies (14% of eLARS), 8 PR galaxies (29% of eLARS), and 16 RD galaxies (57% of eLARS). This is different compared to the original LARS sample, where 7 galaxies have complex kinematics (50% of LARS), 5 are a PR (36% of LARS), and only 2 are a RD (14% of LARS). This difference is a consequence of the sample selection. Whereas LARS galaxies had to satisfy an Hα equivalent width cut of EW_Hα ≥ 100 Å, this requirement was lowered to EW_Hα ≥ 40 Å for eLARS to remove the bias towards irregular and merging systems M23, which are characterised by CK. The combined sample of all 42 eLARS and LARS galaxies consists of 18 (43%) RD galaxies, 13 (31%) PR galaxies, while 11 (26%) galaxies of the sample exhibit CK.

This almost three-way split amongst the kinematical categories bears resemblance to the split found in early IFU studies of star-forming galaxies at similar masses in the high-z (z ≳ 1) universe (Glazebrook 2013, and references therein). More recent surveys, which use stricter criteria to define rotating discs, report disc fractions in excess of 50 % (Wisnioski et al. 2015, 2019), but that could be an overestimate due to mergers acting as disc impostors (Simons et al. 2019).

3.2.2. Comparing Lyα observables between kinematical classes

We may imagine that interstellar medium conditions in galaxies with complex velocity fields are favourable for Lyα escape. For example, irregular and merging systems are often characterised by high specific SFRs, which lead to spatially concentrated injection of momentum and energy from stellar feedback processes. The resulting winds and outflows then may lead to a more effective clearing of escape channels for Lyα (and also Lyman continuum) radiation. The 42 galaxies of the combined LARS and eLARS sample allow for a statistical exploration of such a scenario.

We begin by showing stacked histograms of the Lyα observables colour-coded by their kinematical class in our sample in Fig. 2. The maximum EW_Lyα (56.7 Å; LARS 2), maximum $f_{\text{esc}}^{\text{Ly}\alpha }$ (0.3; LARS 2), and maximum L_Lyα (5.6 × 10⁴² erg s⁻¹; LARS 14) values are found indeed among the CK objects. Moreover, also the averages of those observables are larger in the CK class (⟨EW_Lyα⟩ = 23.15 Å: ⟨$f_{\text{esc}}^{\text{Ly}\alpha }$⟩ = 0.102; ⟨L_Lyα⟩ = 1.02 × 10⁴² erg s⁻¹) than the averages found for the PR (⟨EW_Lyα⟩ = 17.47 Å; ⟨$f_{\text{esc}}^{\text{Ly}\alpha }$⟩ = 0.067; ⟨L_Lyα⟩ = 2.75 × 10⁴¹ erg s⁻¹) and for the RD (⟨EW_Lyα⟩ = 16.36 Å; ⟨$f_{\text{esc}}^{\text{Ly}\alpha }$⟩ = 0.057; ⟨L_Lyα⟩ = 3.10 × 10⁴¹ erg s⁻¹) galaxies, respectively. However, also non-detections of Lyα emission are found amongst objects in all three classes.

Fig. 2. Stacked histograms, colour coded according to kinematical class (rotating discs: purple; perturbed rotators: green; complex kinematics: orange) of the Lyα observables EWLyα (left panel), $f_{\text{esc}}^{\text{Ly}\alpha }$ (middle panel), and LLyα (right panel).

In order to quantify the difference of the observed distributions of Lyα observables amongst the different classes we use the two-sample Kolmogorov-Smirnov (KS) test. Here the six galaxies with upper limits in their Lyα observables are excluded from calculating the KS-test statistic. As mentioned in Sect. 3.2.1, the boundary between our three classes is not well defined, and some PRs could also be classified either as CKs or RDs. For the KS-test we thus regroup our sample into two groups, 𝒜 and ℬ. First, we combine the RDs and PRs into group 𝒜 (27 galaxies) and the CKs into group ℬ (9 galaxies), then we associate only the RDs with group 𝒜 (15 galaxies) and combine the PRs with the CKs in group ℬ (21 galaxies). For both parings and for all three Lyα observables the KS-test is consistent with the null-hypothesis that the samples in the kinematic class are drawn from the same parent population (p-values in the range p_KS ∼ 0.2…0.3). This exercise shows that there is significant overlap in the Lyα observables amongst the kinematical classes. The qualitative nature of a galaxy’s line-of-sight ionised gas velocity field in our sample does not reflect whether it is as a strong Lyα emitter or not.

3.3. Global kinematical measures: v_shear, σ₀^obs, and σ_tot

We now turn to a more quantitative analysis of the observed line-of-sight velocity- and velocity dispersion maps. Therefore we derive two main summarising measures from those maps: the apparent shearing velocity, v_shear, and the intrinsic velocity dispersion, σ₀ (see review by Glazebrook 2013).

For the calculation of v_shear we account for extreme outliers from the distribution of observed velocities. To this aim we calculate $v_{\mathrm{shear}}=\frac{1}{2}(v_{\max }-v_{\min })$ for each galaxy, with v_max and v_min being the values of the 95th percentile and the 5th percentile of the v_los values, respectively. The upper- and lower errors on v_los are estimated from the difference of the 95th and 5th percentiles, respectively, to the mean of the values outside of each percentile.

For σ₀ we adopt the notion of a galaxy’s average line-of-sight velocity dispersion. This simplistic definition implicitly assumes that the σ₀ is isotropic. For measuring this quantity from the observed dispersion fields we need to account for inflated values in individual spaxels due to the convolution of the galaxies’ intrinsic Hα spatial- and spectral- morphology with the atmospheres point spread function and the coarse sampling by 1″×1″ spaxels. In H16 we made no such correction, but the large fraction of RD galaxies in the combined LARS+eLARS sample demands us to be more cautious in this respect. Various methods that attempt to correct for this “beam smearing” or, more correctly, “point-spread function smearing” effect have been presented in the literature (e.g. Davies et al. 2011; Varidel et al. 2016). We compare several such methods against results from morphokinematical modelling of our observations with the GalPaK^3D software (Bouché et al. 2015) in Appendix C. This modelling accounts for the spectral and spatial convolution of the observed Hα kinematics and it adopts also the notion that σ₀ is isotropic. However, GalPak^3D can only be applied to rotating systems. Given the large fraction of more complex systems in our sample we thus require an empirical method for estimating σ₀ directly from the data. In Appendix C we find that the model-based velocity dispersions can be reliably recovered by averaging the σ_v maps, provided that certain spaxels are masked according to a local gradient-measure in the v_los map. The uncertainty on ${\sigma }_0^{\text{obs}}$ is calculated by propagating the Δσ₀’s of the individual spaxels (Sect. 3.1) in quadrature.

The v_shear and ${\sigma }_0^{\text{obs}}$ values as determined above are listed in Table D.1 of Appendix D. We comment that there is a small difference between the v_shear values reported for the LARS galaxies in H16 and the values listed here. This difference is because H16 applied a Voronoi-binning technique to the datacube, whereas here we use each spaxel as is. Moreover, we here used the intervals outside the 95th and 5th percentiles individually to estimate the uncertainty, whereas H16 symmetrised the error-bars. That being so, the here reported updated v_shear measurements and error estimates agree within the error-bars of H16. In contrast, the here tabulated ${\sigma }_0^{\text{obs}}$ values were computed more thoroughly as explained above (see also Appendix C) and are thus slightly different.

We find maximum shearing velocities in eLARS 3 ($v_{\mathrm{shear}}=181.8_{-2.8}^{+1.1}$ km s⁻¹) and LARS 13 ($168.5_{-3.3}^{+18.7}$ km s⁻¹), with the former being a rotating disc and the latter being a merger that exhibits complex kinematics. Minimum shearing velocities are observed in eLARS 28 ($v_{\mathrm{shear}}=14.7_{-2.4}^{+0.9}$ km s⁻¹) and LARS 2 ($v_{\mathrm{shear}}=15.4_{-2.9}^{+8.0}$ km s⁻¹), a perturbed rotator and a complex system, respectively. In Fig. 3 (left panel) we show a stacked histogram of v_shear colour-coded by kinematical class. We note that the distribution of the RDs differs markedly in comparison to the PRs and CKs, with the latter two classes showing overall lower v_shear values. Despite the eLARS extension enlarging the fraction of discs compared to LARS, only one new v_shear extremum is introduced into the combined sample (eLARS 03). Thus, we now sample the v_shear range at hand more densely, as reflected also by the mean and the median of v_shear (57.96 km s⁻¹ and 73.13 km s⁻¹, respectively) being similar or identical to the respective values of the original LARS sample (52.4 km s⁻¹ and 73.13 km s⁻¹). Given that the selection of our sample is limited to young star-bursts, the shearing velocities in our sample are unsurprisingly lower than what is observed for more massive disc galaxies, even when not corrected for inclination.

Fig. 3. Stacked histograms of the global kinematical measures vshear (left panel), ${\sigma }_0^{\text{obs}}$ (middle panel), and vshear/σ0obs (right panel) colour coded according to kinematical class, as in Fig. 2.

For the velocity dispersions the maxima are observed in eLARS 24 (σ₀^obs = 78.0 ± 2.9 km s⁻¹) and LARS 3 (σ₀^obs = 77.3 ± 2.9 km s⁻¹) and the minima are seen in eLARS 16 (σ₀^obs = 23.9 ± 0.9 km s⁻¹) and eLARS 23 (σ₀^obs = 23.9 ± 0.7 km s⁻¹). We show a stacked histogram of ${\sigma }_0^{\text{obs}}$, colour coded by kinematical class, in Fig. 3 (middle panel). There appears not much difference regarding the distributions of ${\sigma }_0^{\text{obs}}$ in the different kinematical classes, except that most of the lowest velocity dispersions are found amongst the RDs and PRs. The eLARS extension introduces a larger fraction of lower velocity dispersion (σ₀^obs ∼ 35 km s⁻¹) galaxies into the sample; this is apparent from Table D.1, and can also be seen by comparing mean and median ${\sigma }_0^{\text{obs}}$ of the full sample (41.1 km s⁻¹ and 37.3 km s⁻¹) to the original sample (55.6 km s⁻¹ and 51.9 km s⁻¹). The on average lower ${\sigma }_0^{\text{obs}}$ values are again a consequence of the lower EW_Hα cut used to define the extended sample (cf. Sect. 3.2.1). This is due to the fact that EW_Hα traces specific SFR (sSFR = SFR per galaxy mass), which tightly correlates with SFR (τ = 0.51, p_τ ∼ 10⁻⁶ for LARS+eLARS; see also analysis by Law et al. 2022), and since the latter, as discussed already in the introduction, tightly correlates with SFR. This is analysed further Sect. 3.5 below.

We also use the ratio v_shear/σ₀^obs that combines the above measurements to a summarising statistic of a galaxy’s kinematical state. These ratios are also listed in Table D.1. The asymmetric error-bars on v_shear were propagated into v_shear/σ₀^obs after moving into log-space and then using method presented in the appendix of Laursen et al. (2019). Our sample is characterised by low v_shear/σ₀^obs values, ranging from 0.35 to 4.46 (median: 1.6, average: 1.8). We show the corresponding stacked-histogram with colour coding by kinematical class in the right panel of Fig. 3. The distributions of v_shear/σ₀^obs differ markedly when comparing the RDs with the CKs and PRs, with the later clearly dominating the low- v_shear/σ₀^obs range. Still, the overall low ratios in the whole sample are below that of typical disc galaxies (v_rot/σ₀ ≳ 8, where v_rot is the maximum rotational velocity; e.g. Epinat et al. 2008) and are more akin to low-mass star-forming galaxies at low- and high-redshifts. Galaxies with $v_{\mathrm{rot}}/{\sigma }_0<\sqrt{3.36}\approx 1.86$ can be classified as dispersion dominated galaxies, as here the turbulent motions exceed the value that would correspond to equal contributions of random motions and rotation to the dynamical support of a turbulent disc (Förster Schreiber & Wuyts 2020; Bik et al. 2022). According to this criterion 24 of the 42 galaxies in LARS+eLARS are dispersion dominated; most of these are PRs or CKs.

We again note that v_shear is not inclination corrected, because such a correction is often ill-defined for the non-rotating complex systems. Thus, strictly speaking, v_shear/σ₀^obs characterises the projected kinematical state. However, we argue that this projected kinematical state, rather than the inclination corrected ratio, is the relevant parameter for relating kinematics with Lyα observables. Especially for disc galaxies the Lyα escape fraction is theoretically expected to be strongly inclination dependent, with edge-on systems having absorbed most of the Lyα photons along the line of sight whereas a face-on system may show higher escape fractions (Laursen & Sommer-Larsen 2007; Laursen et al. 2009b; Behrens & Braun 2014; Verhamme et al. 2012; Smith et al. 2022). At comparable masses, this could lead to an anti-correlation between v_shear/σ₀^obs and $f_{\text{esc}}^{\text{Ly}\alpha }$ for disc galaxies, since face-on discs should have preferentially lower v_shear/σ₀^obs. A deprojected ratio could not trace such types of relations that are the subject of this study. We analyse the inclination dependence on the Lyα observables further in Sect. 3.4.

3.4. Inclination dependence on Lyα escape

To date the most detailed Lyα radiative transfer simulations were carried out for massive disc galaxies (Behrens & Braun 2014; Verhamme et al. 2012; Smith et al. 2022). These single galaxy simulations found a strong inclination dependence of EW_Lyα and $f_{\text{esc}}^{\text{Ly}\alpha }$ (see also Laursen & Sommer-Larsen 2007). The highest EW_Lyα and $f_{\text{esc}}^{\text{Ly}\alpha }$ are found when the simulated galaxies are seen face-on, whereas Lyα is significantly extinguished when the galaxies are viewed edge-on. A plausible conjecture from those simulation results is that high-z samples selected on Lyα emission might be biased towards face-on systems. However, observational results on high-z LAEs, which use the photometric major- to minor axis ratio, b/a, or the flattening, 1 − b/a, as proxy for the inclination (since cos(i) = b/a for infinitely thin discs), are not supportive of this idea (Gronwall et al. 2011; Shibuya et al. 2014; Paulino-Afonso et al. 2018).

In the present contribution, we are equipped with kinematical inclinations, i, from the GalPak^3D disc modelling of the 15 RDs that was performed to verify the empirical masking technique for the calculation of ${\sigma }_0^{\text{obs}}$ (Appendix C.1). We list cos(i) in Table D.1. We compare these measurements to the photometric ratio determined in the I-band by Rasekh et al. (2022), as this quantity minimises the mean difference between photometric axis ratios and kinematic inclinations (⟨cos(i)−(b/a)⟩ = 0.07). Of course, no perfect match is expected, given that the simple approximation cos(i)≈b/a does not account for the finite thickness of the disc. Thus the photometric inclination is biased low for systems that are observed close to edge-on; LARS 11 and eLARS 18 are such cases in our sample. Moreover, Rasekh et al. (2022) calculate the axis ratio from the image moments. However, these measures are sensitive to substructures (i.e. clumps and spiral arms) that are not homogeneously distributed throughout the disc. The moment-based axis ratios can then measure the galaxy as a more elongated structure and thus overestimate the actual inclination; eLARS 4, 17, 19, and 27 are such cases in our sample. Kinematical modelling provides a more accurate measure of the disc inclination.

The inclination dependence of the Lyα observables in our sample is plotted in Fig. 4. It can be seen that there is no preference for larger Lyα observables towards more face-on systems in the LARS+eLARS sample. This is also confirmed by a correlation analysis⁵ using Kendall’s τ (e.g. Puka 2011), with p-values significantly in excess of 0.05. Moreover, no trend is observed when we consider the photometric axis ratios as a proxy for the inclination. Here it needs to be kept in mind that the concept of inclination cannot be applied to CK systems, which is why CK systems are omitted from the plot. Our data obtained for the LARS+eLARS sample does not support the idea that LAE samples are biased towards face-on discs.

Fig. 4. Lyα observables (top panel: LLyα; middle panel: EWLyα; bottom panel: $f_{\text{esc}}^{\text{Ly}\alpha }$ ) vs. galaxy inclination. Filled purple symbols are based on the kinematical inclination, i, from the disc modelling of the RDs with GalPak3D. Unfilled symbols are based on photometric inclinations (purple: RD; green: PR) according to cos(i)≈(b/a)I, where (b/a)I is the I-band axis ratio from Rasekh et al. (2022). Thus, each RD is represented by a filled and an open symbol in these plots. Circles with error bars indicate measurements of the Lyα observables, whereas downward-pointing triangles indicate upper limits. Correlation coefficients (Kendall’s τ) and corresponding p-values for the kinematical inclinations are provided in the top-right corner of each panel.

3.5. Relations between galaxy parameters, M_⋆ and SFR_Hα, and global kinematical statistics

It is established that LAEs at high-z are preferentially low-mass (M_⋆ ≲ 10¹⁰ M_⊙) galaxies with star-formation rates on and slightly above the star-forming main sequence (SFR ∼ 10 M_⊙ yr⁻¹; e.g. Rhoads et al. 2014; Oyarzún et al. 2017; Kusakabe et al. 2018; Ouchi 2019; Pucha et al. 2022; Chávez Ortiz et al. 2023). Kinematically, the SFR is known to tightly correlate with the intrinsic velocity dispersion (e.g. Green et al. 2010, 2014; Law et al. 2022). The rotation velocity, v_max, directly traces mass. For rotating discs v_shear ≤ v_max ⋅ sin(i) holds, where i is the inclination and v_max is the maximum value of the rotation curve; equality holds when the sampled velocity field reaches into the flat part of the rotation curve. Thus, before establishing and interpreting any relation between ionised gas kinematics and Lyα observables, we first need to investigate whether such correlations between physical galaxy parameters and galaxy kinematics exist in LARS+eLARS.

We begin by showing in Fig. 5 histograms of the distribution of M_⋆ and SFR_Hα. These global quantities are measured from SED fitting to the HST images of the galaxies and the relevant details can be found in M23. The sample spans a M_⋆ (SFR_Hα) range from 3.5 × 10⁹ M_⊙ (0.15 M_⊙yr⁻¹) to 1.52 × 10¹¹ M_⊙ (67 M_⊙yr⁻¹), with the mean and median being 3.2 × 10¹⁰ M_⊙ (8.6 M_⊙yr⁻¹) and 2 × 10¹⁰ M_⊙ (2.4 M_⊙yr⁻¹), respectively. The masses sample rather uniformly a decade above and half a decade below the characteristic stellar mass, M_⋆^*, of the stellar mass function of z ∼ 3 LAEs (M_⋆^* ≈ 4 × 10¹⁰ M_⊙; Santos et al. 2021). LAE samples at high-z typically probe objects down to stellar masses of M_⋆ ∼ 10⁷ (Ouchi 2019, his Table 3.3), and the fraction of galaxies at lower masses than what is probed here is dominating in those samples. Nevertheless, here the different kinematic classes appear rather uniformly distributed over the masses and star-formation rates of our sample.

Fig. 5. Stacked histograms, colour-coded according to kinematical class (as in Fig. 2) of M⋆ (left panel) and SFRHα (right panel) from M23.

In Fig. 6 we compare the relations between the kinematic parameters and M_⋆ (left panel) or SFR_Hα (right panel). We also use Kendall’s rank correlation coefficient, τ, and the associated p-value, p_τ, to reject the null-hypothesis of no correlation. The calculation of the uncertainty on τ and the assessment of the robustness of our correlation analyses against the error-bars of the parameters is determined via a Monte-Carlo simulation; further details are given in Appendix E. We recover both the σ₀ vs. SFR_Hα relation and also the v_shear vs. M_⋆ relation at high significance (p_τ < 10⁻⁴). The LARS+eLARS IFS sample thus shows the same correlations as larger samples that relate galaxy kinematics with M_⋆ and SFR (e.g. Green et al. 2014; Moiseev et al. 2015; Barat et al. 2020; Law et al. 2022). Moreover, we also recover a relation between SFR_Hα and v_shear. While this correlation shows a larger scatter and decreased |τ| value than the SFR_Hα vs. ${\sigma }_0^{\text{obs}}$ relation it is also robust with respect to the error-bars. The correlation between SFR_Hα and v_shear may be interpreted as the kinematical imprint of the star-formation main sequence. In Fig. 6 we can also appreciate that RDs have, unsurprisingly, overall higher v_shear/σ₀^obs and higher v_shear values than CKs and PRs.

Fig. 6. Global kinematical measures vs. stellar mass (left panels) and star-formation rate (from Hα; right panels) for the LARS (squares) and eLARS galaxies (circles). The left-hand top-, middle-, and bottom panels show vshear vs. M⋆, σ vs. M⋆, and vshear/σ0obs vs. M⋆, respectively. The right-hand top-, middle-, and bottom panels, show vshear vs. SFRHα, ${\sigma }_0^{\text{obs}}$ vs. SFRHα, and vshear/σ0obs vs. SFRHα. Rotating discs, perturbed rotators, and systems with complex kinematics are shown with purple, green, and orange symbols, respectively. Systems with double-component Hα profiles in some spaxels are shown with open symbols; these galaxies are not included in the calculation of Kendal’s τ and pτ for relations involving ${\sigma }_0^{\text{obs}}$ and vshear/σobs.

Summarising, we find that the eLARS galaxies shows similar star-formation rates as found in high-z LAE samples, while their masses are found to be higher than the average masses such samples. Nevertheless, we recover the known trends between those parameters and galaxy kinematics of star-forming galaxies. We thus regard the eLARS sample as a representative probe for studying the relation between galaxy kinematics and Lyα observables, but being aware of potential caveats that arise from the bias towards higher masses (see Sect. 4).

3.6. Relations between global kinematical statistics and Lyα observables

When combining the observed trends between M_⋆ or SFR_Hα and kinematical trends from Sect. 3.5 with the fact that LAEs are predominantly found amongst low-mass galaxies on and above the star-formation main sequence, it is expected that strong LAEs exhibit high velocity dispersions and small shearing amplitudes. This expectation is met by the LARS galaxies (H16) and was recently also reported for z ∼ 2 and z ∼ 3 galaxies by Foran et al. (2024). With the now enlarged sample and more carefully determined intrinsic velocity dispersions at hand (Appendix C) we can put these trends again to the test.

Using the canonical EW_Lyα ≥ 20 Å boundary to separate the sample into LAEs and non-LAEs, we find that the LAEs are characterised by a on average higher ${\sigma }_0^{\text{obs}}$ than the non-LAEs: ${⟨{\sigma }_0^{\mathrm{obs}}⟩}_{({\mathrm{EW}}_{\mathrm{Ly}\alpha }\ge 20\mathrm{\AA })}=53\pm 13$ km s⁻¹ vs. ${⟨{\sigma }_0^{\mathrm{obs}}⟩}_{({\mathrm{EW}}_{\mathrm{Ly}\alpha }<20\mathrm{\AA })}=34\pm 9$ km s⁻¹. Moreover, while the v_shear values of LAEs and non-LAEs cover similar ranges (${⟨v_{\mathrm{shear}}⟩}_{({\mathrm{EW}}_{\mathrm{Ly}\alpha }\ge 20\mathrm{\AA })}=67\pm 51$ km s⁻¹ vs. ${⟨v_{\mathrm{shear}}⟩}_{({\mathrm{EW}}_{\mathrm{Ly}\alpha }<20\mathrm{\AA })}=67\pm 39$ km s⁻¹), LAEs tend to exhibit lower v_shear/σ₀^obs ratios (${⟨v_{\mathrm{shear}}/{\sigma }_0^{\mathrm{obs}}⟩}_{({\mathrm{EW}}_{\mathrm{Ly}\alpha }\ge 20\mathrm{\AA })}=1.38\pm 1.18$) than the non-LAEs (${⟨v_{\mathrm{shear}}/{\sigma }_0^{\mathrm{obs}}⟩}_{({\mathrm{EW}}_{\mathrm{Ly}\alpha }<20\mathrm{\AA })}=2.00\pm 1.05$). This indicates that the LAEs and non-LAEs are likely drawn from different underlying distributions in ${\sigma }_0^{\text{obs}}$ and v_shear/σ₀^obs. This hypothesis is confirmed by a KS-test on the subsamples resulting from a split of LARS+eLARS at EW_Lyα = 20 Å with p-values of p_KS < 10⁻³ for ${\sigma }_0^{\text{obs}}$ and p_KS = 0.04 for v_shear/σ₀^obs. However, in v_shear the subsamples appear indistinguishable (p_KS = 0.51). This result indicates that LAEs are more turbulent and show lower rotational support than non-LAEs.

We now investigate trends between the Lyα observables and the global kinematical characteristics. We do this graphically in Fig. 7, where we show a 3×3 scatter-plot matrix; the columns trace EW_Lyα, $f_{\text{esc}}^{\text{Ly}\alpha }$, and L_Lyα on the abscissa (from left to right), whereas the rows trace v_shear, ${\sigma }_0^{\text{obs}}$, and v_shear/σ₀^obs on the ordinate (from top to bottom). For each of the nine analysed relations we also compute Kendall’s τ and the corresponding p_τ. The resulting τ and p_τ values are provided in each subpanel of Fig. 7. We adopt the canonical p_τ < 0.05 threshold for considering the scatter not being random. This corresponds to |τ|> 0.22 (|τ|> 0.21) for N = 38 (N = 42, for analysis involving only v_shear; cf. Sect. 3.1) data points (Appendix E). We recover correlations for ${\sigma }_0^{\text{obs}}$ vs. EW_Lyα and ${\sigma }_0^{\text{obs}}$ vs. L_Lyα and we recover anti-correlations for v_shear vs. $f_{\text{esc}}^{\text{Ly}\alpha }$ and v_shear/σ₀^obs vs. $f_{\text{esc}}^{\text{Ly}\alpha }$ by adopting this threshold. All of these correlations are robust against perturbations due uncertainties on the involved measurements by the metric defined in Appendix E. Moreover, the data is suggestive of a anti-correlation between v_shear/σ₀^obs and L_Lyα, although the obtained p_τ is slightly above the canonical boundary of 0.05. We discuss the implications of our findings further in Sect. 4. For now we conclude that the data suggests that the kinematical state of a galaxy are in a causal relation to all three Lyα observables.

Fig. 7. Relations between ionised gas kinematics (vshear, ${\sigma }_0^{\text{obs}}$, vshear/σ0obs in the top-, middle-, and bottom row, respectively) and Lyα observables (EWLyα, $f_{\text{esc}}^{\text{Ly}\alpha }$, and LLyα, in the left-, centre-, and right column, respectively). Symbols are the same as in Fig. 6; additionally here upper limits in the Lyα observables are shown as leftward-pointing triangles. Correlation coefficients (Kendall’s τ) and corresponding p-values are provided in the top-right corner of each panel.

As the correlation analysis tests the data for monotonicity, it fails to capture trends that cannot be described by a monotonic function. In this respect the panels for v_shear vs. EW_Lyα or v_shear/σ₀^obs vs. EW_Lyα in Fig. 7 appear noteworthy. Here we observe that low-EW_Lyα galaxies occupy almost the whole dynamical range sampled in v_shear and v_shear/σ₀^obs, whereas high-EW_Lyα galaxies are found predominantly at lower v_shear and v_shear/σ₀^obs values. More quantitatively, adopting the thresholds of EW_Lyα ≥ 20 Å for defining LAEs and v_shear/σ₀^obs < 1.86 for defining dispersion dominated systems (cf. Sect. 3.3), we find that 9 out of 12 LAEs are dispersion dominated (75%). These 9 LAEs are also characterised by v_shear < 100 km s⁻¹. In this regard it appears also of interest to note that all galaxies classified as RDs in Sect. 3.2.1 do not qualify as dispersion dominated. If the trends found in the LARS+eLARS sample can be generalised, then we may suspect that LAEs are preferentially dispersion dominated systems with low velocity shearing and velocity fields that are, moreover, not commensurate with an unperturbed rotating disc.

4. Discussion

4.1. Assessing the importance of the kinematical parameters for the Lyα observables

Our analysis in Sect. 3.6 shows the potential importance of the global kinematical state of the ionised gas in a galaxy regarding its Lyα observables. We found significant correlations of ${\sigma }_0^{\text{obs}}$ with EW_Lyα and L_Lyα, as well as the anti-correlations of v_shear and v_shear/σ₀^obs with f_esc (Fig. 7). We reasoned that such relations are expected as they are the kinematical imprints of the fact that strong Lyα emitters are preferentially lower-mass systems with higher star-formation rates in comparison to weaker Lyα emitters or even absorbers. Of course, higher star-formation rates result in a larger amount of intrinsically produced Lyα photons, but we may further reason that feedback from star formation is creating an environment that facilitates easy escape of those photons. This may be especially the case in low-mass galaxies, where stellar winds and supernovae explosions can remove gas more efficiently because of their shallower gravitational potentials.

We recall that other galaxy characteristics are known to significantly affect the Lyα observability. In particular, the size and the dust content are of high relevance. Studies both at high-z and low-z show that LAEs are preferentially compact systems with only small amounts of dust (e.g. Guaita et al. 2010; Paulino-Afonso et al. 2018; Marchi et al. 2019; Paswan et al. 2022; Napolitano et al. 2023). The low dust content of LAEs is also reflected in their young ages and low metallicities. Moreover, the outflow kinematics and covering fraction of the scattering medium, probed by low-ionisation absorption lines in the rest-frame UV, are important; Lyα emitting galaxies are found to exhibit lower covering fractions and velocity offsets indicative of outflows (e.g. Wofford et al. 2013; Rivera-Thorsen et al. 2015; Trainor et al. 2015; Östlin et al. 2021; Reddy et al. 2022; Hayes 2023). Lastly, Lyα emitters exhibit higher degrees of ionisation and overall a harder UV radiation field (Trainor et al. 2016; Hayes et al. 2023). Clearly, Lyα production and the subsequent escape depend on a multitude of factors, and it is not at all intuitive which are the most important. Thus, we now want to quantify the importance of the global ionised gas kinematics regarding $f_{\text{esc}}^{\text{Ly}\alpha }$, EW_Lyα, and L_Lyα (hereafter responses) relative to the other global observational and observationally inferred physical galaxy characteristics (hereafter parameters) of LARS+eLARS.

Assessing the relative parameter importance from a set of parameters with respect to a response is a difficult statistical problem. It is an area of active and controversial research and there exists no general agreed upon solution (see review by Grömping 2015). A main problem are intra-correlations between parameters, which lead to fundamental problems in defining the concept of importance. This hampers the robustness of the interpretation of heuristic approaches. Moreover, the more robust approaches implicitly assume that measurement errors are negligible, a situation unfortunately not encountered for the here considered responses and parameters. We thus opt for using rather simple metrics that were already applied in the context of understanding Lyα emission from galaxies (Runnholm et al. 2020, hereafter R20; Napolitano et al. 2023; Hayes 2023). Nevertheless, we also need to be aware of their shortcomings.

First, we use the absolute values of Kendall’s τ for each parameter on its own (Sect. 4.1.1). This so-called “marginal perspective” was recently used by Napolitano et al. (2023) alongside a Random Forrest classifier for identifying the important parameters that regulate EW_Lyα in high-z galaxies. The downside of the marginal perspective is that seemingly unimportant parameters may help each other in influencing the response variable (e.g. Guyon & Elisseeff 2003). For example, it may not be sufficient that a galaxy shows high O₃₂ for having a large EW_Lyα, but the galaxy is also required to have a low F_cov – both parameters together may thus be of high-importance. Thus, we also need to analyse the problem from a “conditional perspective” (Grömping 2015). To this aim we first analyse all intra-parameter correlations in our sample before we follow R20 and Hayes et al. (2023) and adopt a procedure known as stepwise-regression in a multi-variable linear regression framework (Sect. 4.1.2).

As in R20 we define a purely observational parameter set and a physical parameter that is derived from the observational parameters. The observational parameter set consists of luminosities in the FUV, U, B, and I-band (expressed as luminosity densities L_X, with X ∈ {FUV, U, B, I}), emission line luminosities of the lines Hα, Hβ, [O III] λ5007, [O II] λ3727 + λ3729, [N II] λ6584 (L_Y, with Y ∈ {Hα,Hβ,[O II],[O III][N II]}). The physical parameter set consists of stellar mass M_⋆, the flux ratio between the [O III] and [O II] lines, O₃₂ = F_{[O III],λ5007}/(F_{[O II],λ3727} + [O II],λ3729) as proxy for the degree of ionisation, the nebular dust attenuation, E(B − V), the star-formation rate, SFR_Hα, the nebular oxygen abundance, 12 + log(O/H), the UV size s_UV, as well as three parameters from low-ionisation absorption line measurements from HST/COS spectroscopy, namely the covering fraction F_cov, the line width, w₉₀, and the velocity offset, v₉₅, with respect to the systemic redshift determined from the nebular emission lines. All parameters are taken from R20 and M23; these publications detail the measurements of most of those parameters and we also refer to Rivera-Thorsen et al. (2015) and Hayes et al. (2023) for a description of the HST/COS based parameters. Both parameter sets consist of 9 parameters, which are again summarised in Table 1.

Lastly, in the following analyses we remove all six galaxies without detected Lyα emission from the statistical analyses presented below, since methods of the conditional perspective cannot deal with upper limits. As mentioned in Sect. 3.1, galaxies with double component Hα profiles are also not considered in our statistical analyses (i.e. we also have to remove additionally the three galaxies LARS 9, eLARS 7, and eLARS 24; LARS 13 has both double components and upper limits in the Lyα observables). Thus, the following analyses consider a sample of N = 33 galaxies.

4.1.1. Importance of kinematics: Marginal perspective

The marginal perspective regarding the importance of observational and physical parameters with respect to the responses EW_Lyα, $f_{\text{esc}}^{\text{Ly}\alpha }$, and L_Lyα in LARS+eLARS is presented graphically in Fig. 8. We include the global kinematical statistics v_shear, ${\sigma }_0^{\text{obs}}$, and v_shear/σ₀^obs in each of the two analysed parameter sets to asses their relative importance regarding the responses in comparison to the other parameters. In Fig. 8 we rank the parameters by the absolute values of Kendall’s τ. The rationale here is that |τ| is a direct measure of the frequency at which the change of a parameter leads to a parallel or anti-parallel change of the response. The higher the frequency of changes in the same direction, the higher the value of |τ|. Thus, an ordering of the parameters by |τ| can be interpreted as a ranking of their importance in influencing the response. As before, we deem a correlation significant if p_τ < 0.05, which corresponds to |τ|> 0.239 (indicated by a vertical dashed line in Fig. 8) for N = 33. Moreover, the correlation is deemed robust if the inter-quartile range of τ’s from a Monte-Carlo simulation, that takes errors on the parameters and observables into account, does not overlap with this threshold (Appendix E).

Fig. 8. Importance of observational and physical parameters for $f_{\text{esc}}^{\text{Ly}\alpha }$, EWLyα, and LLyα using the absolute value of Kendall’s τ (marginal perspective; Sect. 4.1.1). In each panel, the parameters are ranked according to |τ| from top to bottom. The length of each bar indicates the absolute value of τ and the error bars indicate the upper- and lower quartile from the robustness analysis presented in Appendix E. Each bar is coloured in red or blue depending on whether τ > 0 (correlated trend between parameter and response) or τ < 0 (anti-correlated trend between parameter and response), respectively. The vertical dashed red line indicates the threshold |τ|> 0.239, which corresponds to pτ < 0.05 for N = 33 (Appendix E). The left, centre, and right panels show the results for the response EWLyα, $f_{\text{esc}}^{\text{Ly}\alpha }$, and LLyα, respectively, where the left and right subpanels show the rankings for the observational and physical parameters, respectively.

Focusing first on the response EW_Lyα we find that ${\sigma }_0^{\text{obs}}$ dominates in importance over all parameters in the observational- and physical parameter sets. On the contrary, v_shear is found to be least influential and, moreover, v_shear/σ₀^obs appears also irrelevant. Formally, the emission line luminosities are ranked directly below ${\sigma }_0^{\text{obs}}$, followed by the broad-band luminosity densities. However, here only the correlation with L_Hβ is statistically significant and robust. Lastly, consistent with the literature (e.g. Law et al. 2012; Paulino-Afonso et al. 2018), we find an anti-correlation between EW_Lyα and s_UV. Taken together, these trends suggest that compact galaxies that are bright in the rest-frame optical emission lines, and which are characterised by a highly turbulent interstellar medium, are more likely to show higher EW_Lyα.

Next we focus on the response $f_{\text{esc}}^{\text{Ly}\alpha }$. The kinematical parameters v_shear and v_shear/σ₀^obs, while not being highly ranked for EW_Lyα, emerge as highly influential parameters for f_esc. Both anti-correlations dominate over all other observational parameters, albeit the correlation with v_shear/σ₀^obs is not robust (Appendix E. Notably, except for L_[O III] and L_[N II], all observational parameters show even |τ|< 0.1 and appear thus individually completely irrelevant with respect to $f_{\text{esc}}^{\text{Ly}\alpha }$. This is also the case for ${\sigma }_0^{\text{obs}}$. Regarding relations between physical parameters and $f_{\text{esc}}^{\text{Ly}\alpha }$ we find that the anti-correlations with E(B − V) and F_cov are ranked higher than the anti-correlations with v_shear. Moreover, the value of |τ| for v_shear appears on par with the coefficient recovered for M_⋆.

The high importance of ${\sigma }_0^{\text{obs}}$ and simultaneous irrelevance of v_shear or v_shear/σ₀ for regulating EW_Lyα appears noteworthy, given that the situation is reverse for $f_{\text{esc}}^{\text{Ly}\alpha }$. This situation is also evident in Fig. 7. Taken at face value, this is suggestive of different mechanisms regulating EW_Lyα and $f_{\text{esc}}^{\text{Ly}\alpha }$. This appears intriguing, because a tight linear relation between $f_{\text{esc}}^{\text{Ly}\alpha }$ and EW_Lyα is found in our sample (M23), as well as in samples at higher redshifts (Sobral & Matthee 2019; Roy et al. 2023; Tang et al. 2023). However, in the following Sect. 4.1.2 we show that ${\sigma }_0^{\text{obs}}$ also emerges as an important parameter for regulating $f_{\text{esc}}^{\text{Ly}\alpha }$ from the conditional perspective, but only when taken in concert with other parameters that are known to influence the Lyα radiative transfer.

Here we recover another intriguing discrepancy regarding the influence of individual parameters with respect to the responses $f_{\text{esc}}^{\text{Ly}\alpha }$ and EW_Lyα. Namely, we find that E(B − V) and M_⋆ are not correlated with EW_Lyα, whereas both parameters anti-correlate with $f_{\text{esc}}^{\text{Ly}\alpha }$. Kinematically, the presence or absence of an anti-correlation between M_⋆ and $f_{\text{esc}}^{\text{Ly}\alpha }$ or EW_Lyα, respectively, is reflected in our sample by the different behaviour of both parameters with respect to v_shear, which we just discussed. Moreover, this discrepancy does not vanish under the conditional perspective in Sect. 4.1.2 and it is appears also in contrast to some literature results (Marchi et al. 2019; Napolitano et al. 2023). Nevertheless, other studies also do not report monotonic relations between EW_Lyα and E(B − V) or M_⋆ (Atek et al. 2014; Bolan et al. 2024).

Finally, we focus on the response L_Lyα. Here we find that all the other observational parameters are ranked above the kinematical parameters. However, ${\sigma }_0^{\text{obs}}$ is ranked in second place amongst the physical parameters below SFR_Hα; v_shear and v_shear/σ₀^obs appear in the bottom half of the ranks and are thus not important in influencing L_Lyα. The strong relation between the emission line luminosities of Hβ and Hα for L_Lyα indicates that the production of intrinsic Lyα is a dominating cause in regulating the total Lyα photon output. This normalisation of Lyα luminosity by recombination rate is also reflected in the physical parameter set, where SFR_Hα (∝L_Hα) is ranked first.

As already mentioned in Sect. 4.1, we need to refrain from over-interpreting the rankings based on the marginal perspective for getting a grip on the physics that regulate the escape and observability of Lyα emission. For now we only conclude that our findings here are not inconsistent with the hypothesis that the global gas kinematics play a very important role in regulating the Lyα observables of a galaxy.

4.1.2. Importance of kinematics: Conditional perspective

In order to move away from the marginal perspective, we first analyse intra-parameter correlations in our parameter sets from Table 1. To this aim we use heat-map visualisations in Fig. 9. Using Kendall’s τ as a correlation measure Fig. 9 visualise the strength and directions of possible monotonic relations between all possible parameter combinations. The aim of this analysis is primarily to illustrate possible caveats when analysing only single parameter correlations with respect to the Lyα observables. We therefore also refrained from carrying out a robustness analysis on τ, instead we choose a more stringent p-value of p_τ = 0.03 or, correspondingly, |τ|> 0.265 for N = 33 (Eq. (E.2)). Intra-parameter relations resulting in a smaller |τ| are crossed out in Fig. 9.

Fig. 9. Correlation matrix heat map visualisation (Kendall’s τ) for observational parameters (left panel) and physical parameters (right panel), including vshear, ${\sigma }_0^{\text{obs}}$, and vshear/σ0obs, as well as the Lyα observables $f_{\mathrm{esc}}^{\mathrm{Ly}\alpha }$, EWLyα, and LLyα. Statistically insignificant trends (i.e. where p0 > 0.05, or correspondingly |τ|> 0.239 for N = 33) are crossed out.

Considering the observational parameters in Fig. 9 (left panel) it becomes apparent that galaxies brighter in the continuum (irrespective of the band) are also almost always brighter in the rest-frame optical emission lines (τ ≳ 0.5, p_τ ≲ 10⁻³) and also in the Lyα line. This “more is more” effect regarding L_Lyα might be specific to the here analysed sample that excludes Lyα absorbers. Regarding our kinematical parameters, only correlations between ${\sigma }_0^{\text{obs}}$ and the emission lines, including Lyα, appear significant⁶, and references therein).

Focusing on the physical parameters in Fig. 9, we see that M_⋆ shows the highest number of intra-correlations with other parameters; only lg(O/H) and E(B − V) are found not to be correlated with M_⋆. Numerous other significant intra-correlations exist amongst the physical parameters. Amongst the kinematical parameters v_shear/σ₀^obs correlates with v_shear and anti-correlates with ${\sigma }_0^{\text{obs}}$. We also note that all Lyα observables correlate with each other. While some of the uncovered monotonic relations are certainly imprints of astrophysical scaling relations, other correlations or non-correlations might be specific to the sample at hand. We not discuss this further here, as this is not the main aim of our study. Of partilcular interest in our case, the high number of intra-correlations highlights the complex interdependence between the individual parameters in the observational and in the physical parameter sets. This complex network of relationships amongst the analysed parameters complicates the problem of assessing parameter importance, but it also highlights the need for a more sophisticated approach than just looking at parameter pairings.

We opt for a simple heuristic approach based on multi-variate regression,

$$\begin{array}{c}\hfill y_i={\displaystyle \sum _{i=1}^{N_p}}{c}_i{x}_i,\end{array}$$(3)

where the x_i are the N_p = 12 parameters from the observational- and physical parameter sets, respectively, including the kinematical parameters. The y_i in Eq. (3) are the responses, i.e. here the Lyα observables, and the c_i denote the coefficients maximising the coefficient of determination,

$$\begin{array}{c}\hfill R^2=1-\frac{{\displaystyle {\sum }_{i=1}^N}{(y_i-f_i)}^2}{{\displaystyle {\sum }_{i=1}^N}{(y_i-⟨y_i⟩)}^2},\end{array}$$(4)

where the f_i are the predictions of the observable y_i and ⟨y_i⟩ is the mean over all y_i. Formally, R² encapsulates the variance explained by the empirical model of Eq. (3). If R² ∼ 0, the linear relation in x_i is not usable to describe and predict the y_i, whereas for R² ∼ 1 Eq. (3) is a nearly perfect empirical model.

In above framework we apply a forward and backward selection technique (also known as step-wise regression) to establish importance rankings for the parameters. Here we iteratively include (forward selection) or exclude (backward selection) parameters from Eq. (3), such that the included or excluded parameter leads to the largest increase or smallest decrease of R², respectively. Moreover, this procedure requires that the parameters are standardised to a comparable scale. As in Hayes et al. (2023), here we use Z-score normalisation, that is, we subtract the mean over all values for each parameter and then we divide the mean shifted values by their standard deviation. Finally, we account for the statistical errors on the parameters and responses by performing N_MC = 10³ Monte-Carlo realisations of the above procedure. In each realisation we perturb both the response and the parameters according to their statistical uncertainties. This procedure provides us then with N_MC ranks for each parameter x_i under forward- and backward selection; we write for these j ranks r_j^FW(x_i) and r_j^BW(x_i), respectively, where j = 1, …, N_MC.

We find that in almost all realisations the rankings obtained under forward selection are not completely consistent with the rankings obtained under backward selection. This so called “violation of anonymity” of the parameters regarding their positions in the regression is driven by the parameter intra-correlations (Fig. 9). It is a fundamental flaw of step-wise regression (Grömping 2015). This flaw also remains in bi-directional step-wise regression, where parameters are included and excluded in each iteration (Smith 2018). A potential alternative are more computational intensive methods related to game theory, which assign relative importance to a parameter based on all possible permutations that a parameter can have with all other parameters (Grömping 2015). However, given that we also require to include the measurement uncertainties via Monte Carlo simulations in our importance estimate, this approach is computationally expensive and will be the subject of future work. Here we instead opt for an ad hoc definition of a final ranking score. With

$$\begin{array}{c}\hfill r(x_i)=\frac{1}{2N_{\mathrm{MC}}}{\displaystyle \sum _{j=1}^{N_{\mathrm{MC}}}}{r}_j^{\mathrm{FW}}(x_i)+r_j^{\mathrm{BW}}(x_i)\end{array}$$(5)

we define the rank importance r(x_i) of parameter x_i as the average over all ranks obtained under forward and backward selection for a particular Monte Carlo realisation. The smaller r(x_i), the more important is x_i. The so obtained order of importance for the responses EW_Lyα, $f_{\text{esc}}^{\text{Ly}\alpha }$, and L_Lyα are provided in Table 2.

Finally, following Hayes et al. (2023), we visualise in Fig. 10 the probability distributions of all ranks r(x_i), p(r(x_i)) by means of a Gaussian kernel-density estimator (kernel size 0.5). We highlight the four curves with the smallest r(x_i) in Fig. 10, as these are formally the top ranked parameters according to Eq. (5). Moreover, we deem a high ranking of a parameter as robust if the p(r(x_i)) curve for a given parameter peaks at a low r. Contrariwise, the ranking of a parameter is not considered as robust if the curve shows a flat distribution that encompasses a larger number of possible ranks. As can be seen from Fig. 10, and as we discuss below, for most of the parameters no clear order of importance can be established, even if a parameter is formally ranked high according to Eq. (5). However, for some parameters a clear picture emerges.

Fig. 10. Visualisation of the parameter importance for responses EWLyα (top panels), fescLyα (middle panels), and LLyα (bottom panels) from a conditional perspective using step-wise regression in a multivariate linear regression framework (see Sect. 4.1.2 for details). The panels on the left- and right-hand sides visualise the distribution of rankings p(r(xi)) for the parameters xi from the observational and physical parameter sets, respectively, where each set includes the kinematical parameters that have been derived in this work. Low values of r(xi) indicate a high importance of xi, whereas high values indicate a low importance. All but two of the curves were obtained using a kernel-density estimator (Gaussian kernel, width = 0.5) on the ranks from NMC = 103 Monte Carlo realisations. The exceptions are the ${\sigma }_0^{\text{obs}}$ curves for the response EWLyα (see footnote 7). The four most important parameters according to Eq. (5) are shown with bold lines, whereas the remaining ranks are shown with thin lines and in subdued colours. The legend at the top of each panel is sorted according to the ranks from Eq. (5), which are also listed in Table 2.

The most notable and robust result from this statistical experiment is that σ₀ is found as the most important parameter in determining EW_Lyα (top panels of Fig. 10). In fact, in almost all Monte-Carlo realisations ${\sigma }_0^{\text{obs}}$ emerges as the top ranked parameter, both under forward and backward selection⁷. The dominating importance of σ₀ in regulating EW_Lyα is thus found both from the marginal perspective (Sect. 4.1.1) and from the conditional perspective. We discuss this result in Sect. 4.2 below. Moreover, the distribution for the parameters ranked formally below ${\sigma }_0^{\text{obs}}$ do not show significant peaks. We thus regard the order of those parameters obtained by Eq. (5) not as robust indications of their relative importance with respect to EW_Lyα. The average R² for all 12 parameters regarding the response EW_Lyα is 0.55 ± 0.04 and 0.50 ± 0.04 for the observational- and physical parameter set, respectively.

The distributions that we recover for the observational parameters with respect to $f_{\text{esc}}^{\text{Ly}\alpha }$ appear completely nondescript (centre left panel in Fig. 10). However, here a clearer picture emerges regarding the role of the physical parameter with respect to $f_{\text{esc}}^{\text{Ly}\alpha }$ (centre right panel in Fig. 10). Here the four top ranked parameters – E(B − V), O₃₂, F_cov, and σ_obs – are characterised by single peaked distributions in the left side of the diagram, whereas the remaining parameters show flattened out distributions towards the right side. For both parameter sets we recover an average R² = 0.60 ± 0.05.

Attending the bottom panels of Fig. 10 we see that for L_Lyα the I-band luminosity arises as a top ranked observational parameter, with the other parameters not yielding conclusive rankings. Here R² = 0.966 ± 0.003, but this high R² is a result of the strong positive correlations of almost all parameters with L_Lyα. In the physical parameter set the top three ranked parameters – O₃₂, M_⋆, and F_cov – separate themselves distinctly from the remaining parameters. The R² for the physical parameter set is 0.64 ± 0.05.

4.2. Astrophysical interpretation

Our analysis above leads us to the conclusion that ${\sigma }_0^{\text{obs}}$ appears as a dominating parameter in regulating EW_Lyα when compared to other observational and physical parameters. This result is found both from the marginal – (Sect. 4.1.1) and from the conditional perspective (Sect. 4.1.2). A slightly different picture emerges regarding the role of ${\sigma }_0^{\text{obs}}$ with respect to $f_{\text{esc}}^{\text{Ly}\alpha }$. Here no direct correlation is found (Sect. 3.6 and Sect. 4.1.1), i.e. ${\sigma }_0^{\text{obs}}$ appears unimportant for $f_{\text{esc}}^{\text{Ly}\alpha }$ when taken on its own, however it is found as an important parameter in regulating $f_{\text{esc}}^{\text{Ly}\alpha }$ when considered in concert with other parameters (Sect. 4.1.2). These parameters are the neutral gas covering along the line of sight traced by F_cov, the ionisation state traced by O₃₂, and the dust extinction traced by E(B − V). That these three other parameters are also expected and known to influence the escape of Lyα photons provides credibility to our method and result.

Unmistakably, Eq. (3) from Sect. 4.1.2 cannot replace an astrophysical model. This statistical method just helps to uncover trends in the data beyond the conventional two-parameter correlation analyses. The result from our statistical analysis requires us now to understand the dominant role of ${\sigma }_0^{\text{obs}}$ for EW_Lyα and the role of ${\sigma }_0^{\text{obs}}$ in modulating $f_{\text{esc}}^{\text{Ly}\alpha }$.

4.2.1. The role of the age of the stellar population

The ${\sigma }_0^{\text{obs}}$ values in our sample are all significantly above the expected combined velocity dispersions due to thermal and turbulent motions of the individual H II regions (${\sigma }_{\mathrm{therm}}=\sqrt{2k_{\mathrm{B}}T/m_{\mathrm{H}}}=12.85$ km s⁻¹ for T = 10⁴ K). Therefore, the observed large values of ${\sigma }_0^{\text{obs}}$ are dominated by the turbulence of the cold and/or warm neutral medium in which the unresolved H II regions are embedded. The monotonic relation between SFR and σ₀ (Sect. 3.5) thus indicates a causal connection between star-formation and turbulence of the interstellar medium.

There are competing scenarios regarding the main physical mechanism that are responsible for the SFR vs. σ₀ relation. One the one hand, it is suggested that gravitational instabilities, which are required to form stars, may stir-up the interstellar medium and decay as turbulence (e.g. Krumholz & Burkert 2010). On the other hand, feedback from star-formation (i.e. stellar winds and especially supernovae) may act as the main source of energy driving this turbulence (e.g. Bacchini et al. 2020). Observational results are often used to argue in favour of latter scenario (Moiseev et al. 2015; Law et al. 2022), but there are also results that favour the former scenario (e.g. Krumholz & Burkhart 2016; Yu et al. 2021). Recent modelling by Ejdetjärn et al. (2022) indicates that the warm-ionised phase traced by Hα may actually also effectively couple to kinematic energy injected by feedback.

Adopting feedback as the driving force beyond σ₀ as a working hypothesis we recall that for a simple stellar population the amount of energy injected by stellar winds and supernovae is related to the duration since the star-burst event (burst age, e.g. Leitherer et al. 1999). Furthermore, we recall that the intrinsic equivalent width, EW_Lyα^int, is controlled by the ionising photon production efficiency, ξ_ion, and both ξ_ion and EW_Lyα are controlled by the burst age and metallicity (Charlot & Fall 1993; Schaerer 2003; Raiter et al. 2010; Inoue 2011). Thus, the evolutionary state of the stellar population is inherently linked to Lyα production, and given our working hypothesis and the results from Sect. 3.6 and Sect. 4.1.2, should also affect the escape of Lyα photons and modulate the observed equivalent width. Indeed, this idea finds support by the recent work of Hayes et al. (2023), who report a strong anti-correlation between EW_Lyα and burst age in their sample of 87 low-z galaxies observed with HST/COS. These authors also report an equally strong anti-correlation between burst age and $f_{\text{esc}}^{\text{Ly}\alpha }$. If our hypothesis is true, then the effect that burst age has on EW_Lyα and $f_{\text{esc}}^{\text{Ly}\alpha }$ should also be identifiable as an effect on ${\sigma }_0^{\text{obs}}$.

In order to analyse the dependence of ${\sigma }_0^{\text{obs}}$, EW_Lyα, and $f_{\text{esc}}^{\text{Ly}\alpha }$ on stellar age, we require estimates of the latter. The work by M23 provided age estimates for each Voronoi cell from the HST imaging data. We here use an average of those ages by weighing with the extinction corrected UV continuum. We note that our estimate differs from the age averages tabulated in M23, as these were only weighted by the extinction uncorrected UV continuum. As shown in Appendix F, where we also provide a table of the extinction corrected ages, for the majority of galaxies we now obtain significantly younger ages since the dust-uncorrected ages did not account for the very young populations enshrouded in dust.

We visualise the dependence of ${\sigma }_0^{\text{obs}}$, EW_Lyα, and $f_{\text{esc}}^{\text{Ly}\alpha }$ on the so derived burst ages in Fig. 11. None of the three parameters are found to correlate with burst age. Hence, the anti-correlations between burst age and EW_Lyα or $f_{\text{esc}}^{\text{Ly}\alpha }$ that were recovered in the sample of Hayes et al. (2023) are not found in our sample. However, a comparison with this sample is not on an equal footing. The crucial methodological difference is that Hayes et al. (2023) estimated stellar ages directly from fitting stellar templates to UV and optical spectra, whereas our ages are derived photometrically. Especially the UV spectra contain many age sensitive features (e.g. Chisholm et al. 2019). Nevertheless, the lack of a dirrect correlation between starburst age and ${\sigma }_0^{\text{obs}}$ may indicate that kinetic energy released from the stellar population and the conversion to turbulent energy follows more complex patterns.

Fig. 11. Relations between ${\sigma }_0^{\text{obs}}$ (left panel), EWLyα (middle panel), and $f_{\text{esc}}^{\text{Ly}\alpha }$ (bottom panel) with the dust-corrected burst age. Kendall’s τ and the respective p-value are provided in the bottom left of each panel. Different symbols are used for rotating discs (circles), perturbed rotators (squares), and galaxies exhibiting complex kinematics (stars). We colour code the points according to EWLyα (left and right panel) or $f_{\text{esc}}^{\text{Ly}\alpha }$ (middle panel). We also highlight points with red outlines that show LLyα/LHα above the maximum dust screen expectation from Eq. (6). As visual aid for the interpretation we show model calculations for the mechanical luminosity, Lmech, and the intrinsic Lyα equivalent width, EWLyαint, as a function of burst age with grey dotted lines. We note that the ticks on the right side of each panel provide the scale for these model calculations. In particular, we show the Lmech of a Z = 0.008 starburst normalised to 1 M⊙ in the left panel and we indicate the temporal evolution of EWLyαint for a Z = 0.008 burst and the EWLyαint range for constant star-formation at Z = 0.001 − 0.040 in the middle panel (see text for details).

Fig. 12. Lyα equivalent width, EWLyα, vs. nebular extinction, E(B − V). Points are coloured according to ${\sigma }_0^{\text{obs}}$ as indicated by the colour bar on the right. Different symbols indicate the kinematical classes of the galaxies introduced in Sect. 3.2. The dotted blue line indicates the zone of EWLyα ≤ 20 Å and E(B − V)≤0.4 discussed in the text. There is no correlation between both quantities (τ = −0.049/pτ = 0.7).

To analyse this further we plot in Fig. 11 (upper panel) the evolution of mechanical luminosity (normalised to 1 M_⊙) according to the STARBURST99 stellar population model (Leitherer et al. 1999) for a Z = 0.008 population (with a power law initial mass function of slope α = −2.35 and upper mass cut-off at 100 M_⊙). To exemplify how EW_Lyα^int depends on the burst-age we plot in the in the middle panel of Fig. 11 the EW_Lyα^int evolution with time for a star-burst with Z = 0.008 Z_⊙ according to the calculations by Inoue (2011). We also show the range of EW_Lyα^int for similar types of low-Z populations with constant star formation from the calculations of Schaerer (2003). Considering the idealised star-burst, where all stars form instantaneously, initally winds from the highest mass stars (O stars) dominate the mechanical energy output. During these first few 10⁶ years after the burst, the intrinsic Lyα equivalent widths are expected to be very high, with EW_Lyα^int ∼ 300 Å at the typical metallicities (Z ∼ 10⁻²) of our sample. After ≈5 × 10⁶ years the injection of kinematic energy will increase substantially due to the most massive stars going supernova. The absence of ionising photons from these stars in otherwise very blue continua will lead to a gradual decrease of EW_Lyα^int. At ≈10⁷ years the O-stars are gone, the ionising luminosity drops by a factor of ∼100 compared to the initial phase of the burst, and the EW_Lyα^int decreases even further. At this phase the now the injected kinetic energy stems only from stellar winds of lower mass stars.

Of course, the galaxies here consist of multiple clusters of different ages. Moreover, turbulent energy dissipates over varying timescales (∼10 − 100 Myr; e.g. Bacchini et al. 2020), so the link between mechanical luminosity or EW_Lyα^int with burst age in Fig. 11 serves only as a schematic illustration. Still, in such a simple scenario, we would expect the lowest ${\sigma }_0^{\text{obs}}$ values for the oldest population, but this is not what we observe. Also the trend for the EWs appears more nuanced. Except for one outlier, this plot is characterised by a zone of avoidance of old objects with high-EWs. Moreover, the highest EWs indeed appear at the youngest ages. Lastly, no discernible trend appears in the $f_{\text{esc}}^{\text{Ly}\alpha }$ plot. We conclude that we do not recover the anticipated anti-correlation of burst-age on EW_Lyα^int and $f_{\text{esc}}^{\text{Ly}\alpha }$, nor can we establish a causal relation between burst-age and ${\sigma }_0^{\text{obs}}$.

4.2.2. Gas kinematics can shape dust geometry

We analyse in Fig. 12 the distribution of our sample in the EW_Lyα vs. E(B − V) plane. In this diagram we find a zone of avoidance: the five objects with the highest EW_Lyα measurements are found at low extinction (E(B − V) < 0.4), whereas the objects with lower EW_Lyα’s occupy the full range in E(B − V). The distribution of galaxies in the EW_Lyα vs. E(B − V) plane appears similar as in the compilation of low-z starbursts presented by Atek et al. (2014). We colour code the points in Fig. 12 by ${\sigma }_0^{\text{obs}}$ to visualise the influence of this parameter on the observed EW_Lyα vs. E(B − V) distribution. While EW_Lyα is strongly correlated with ${\sigma }_0^{\text{obs}}$ (Sect. 3.6 and Sect. 4.1.1), we formally have to accept the null-hypothesis that ${\sigma }_0^{\text{obs}}$ and E(B − V) are not correlated (see also Fig. 9, numerically τ = 0.223 / p₀ = 0.07). Still, in the E(B − V) vs. EW_Lyα plane it is apparent, how the corner of low extinction objects (0 ≤ E(B − V)≤0.4) with low Lyα equivalent widths (EW_Lyα ≤ 20 Å) appears to be dominated by galaxies with lower ${\sigma }_0^{\text{obs}}$, whereas objects outside of this corner appear to be characterised by higher ${\sigma }_0^{\text{obs}}$. These boundaries may appear somewhat arbitrary, but EW_Lyα ≥ 20 Å represents the canonical cut to define LAEs, while E(B − V) = 0.4 halves the E(B − V) range probed by our sample. We find that 18 objects with EW_Lyα ≤ 20 Å and 0 ≤ E(B − V)≤0.4 show a mean (median) ${\sigma }_0^{\text{obs}}$ of 32.8 km s⁻¹ (31.2 km s⁻¹), whereas the other 15 galaxies are characterised by a mean (median) ${\sigma }_0^{\text{obs}}$ of 50.3 km s⁻¹ (51.2 km s⁻¹). Thus, excluding objects in this low-σ₀^obs and low-E(B − V) region leaves us with objects that show simultaneously EW_Lyα > 20 Å and high velocity dispersions, but we also find galaxies with higher dust extinction that still show appreciable Lyα equivalent widths (EW_Lyα ≈ 20 Å) and, intriguingly, their intrinsic velocity dispersions are also high.

Next, we analyse in Fig. 13 the distribution of our sample in the $f_{\text{esc}}^{\text{Ly}\alpha }$ vs. EW_Lyα plane. There we also show the linear relation $f_{\text{esc}}^{\text{Ly}\alpha }$ ∝ 4 × 10⁻³ EW_Lyα from M23⁸. As in Fig. 12, we colour code the points by ${\sigma }_0^{\text{obs}}$. We furthermore now highlight objects outside of the low-EW_Lyα – low-E(B − V) region defined in Fig. 12. We see that the dominant fraction of objects that are below the $f_{\text{esc}}^{\text{Ly}\alpha }$ ∝ 4 × 10⁻³ EW_Lyα line are outside this E(B − V)≤0.4 and EW_Lyα ≤ 20 Å box of Fig. 12 (9 out of 12 objects). Thus, galaxies that have low $f_{\text{esc}}^{\text{Ly}\alpha }$ but higher EW_Lyα are characterised by high ${\sigma }_0^{\text{obs}}$.

Fig. 13. Lyα escape fraction, $f_{\text{esc}}^{\text{Ly}\alpha }$ vs. Lyα equivalent width, EWLyα for the in Sect. 4 analysed subsample of 33 galaxies as in Fig. 12. Both quantities are highly correlated (τ = 0.453/pτ ∼ 10−3), even in a linear sense (Pearson’s r = 0.78/pr ≃ 10−7). The dotted line shows the linear fit to the whole sample from M23. Points are coloured according to ${\sigma }_0^{\text{obs}}$ as indicated by the colour bar on the right. Symbols are the same as in Fig. 12, but here we highlight by blue outlines the objects outside of the EWLyα and E(B − V)≤0.4 zone in Fig. 12.

Interestingly, almost all of the outliers below the f_esc ∝ 4 × 10⁻³ EW_Lyα relation appear incompatible with extinction by a homogeneous dust screen. This can be appreciated in Fig. 14 where we plot the Lyα/Hα luminosity ratio against the Balmer decrement, Hα/Hβ, for our sample. In Fig. 14 we also indicate the expectation of Lyα/Hα as a function of Hα/Hβ for extinction by a homogeneous dust-screen in the absence of radiative transfer effects:

$$\begin{array}{c}\hfill \frac{L_{\mathrm{Ly}\alpha }}{L_{\mathrm{H}}\alpha }={\left(\frac{L_{\mathrm{Ly}\alpha }}{L_{\mathrm{H}}\alpha }\right)}_{\mathrm{int}}·\frac{e^{-{\tau }_{d,\mathrm{Ly}\alpha }}}{e^{-{\tau }_{d,\mathrm{H}\alpha }}}.\end{array}$$(6)

Here

$$\begin{array}{c}\hfill {\tau }_{d,\mathrm{Ly}\alpha }=0.4·ln(10)·k({\lambda }_{\mathrm{Ly}\alpha })·E(B-V)\end{array}$$(7)

and

$$\begin{array}{c}\hfill {\tau }_{d,\mathrm{H}\alpha }={\tau }_{\mathrm{Ly}\alpha }·(k({\lambda }_{\mathrm{H}}\alpha )/k({\lambda }_{\mathrm{Ly}\alpha }))\end{array}$$(8)

are the dust optical depths at λ_Lyα and λ_Hα, respectively; the definitions for (L_Lyα/L_Hα)_int, k(λ), and E(B − V) were given in Sect. 2.4.

As discussed extensively in the literature, the intuition that the observed Lyα/Hα is always below the ratio given in Eq. (6) due to the loss of Lyα photons by radiative transfer effects proves often to be wrong (Scarlata et al. 2009; Atek et al. 2014; Hayes et al. 2014; Bridge et al. 2018; M23). As can be seen from Fig. 14, here we have 9 objects in our sample of 33 galaxies that are significantly above the pure dust attenuation expectation of Eq. (6). We also see in this figure that these 9 galaxies all are outside the EW ≤ 20 Å and E(B − V)≤0.4 box from Fig. 12. Moreover, 8 of those 9 galaxies are below the $f_{\text{esc}}^{\text{Ly}\alpha }$ ∝ 10⁻³EW_Lyα line from Fig. 13. Thus, galaxies that are not compatible with extinction by a homogeneous dust screen are characterised by high intrinsic velocity dispersions.

Fig. 14. Lyα/Hα ratio against the Hα/Hβ luminosity ratio for the in Sect. 4 analysed subsample of 33 galaxies as. there is no correlation between both quantities (τ = −0.09/pτ = 0.47). Symbols, colouring, and highlighting of the symbols is the same as in Fig. 13. The expectation of Lyα/Hα as a function of Hα/Hβ for a homogeneous dust-screen without radiative transfer effects in front of the line emission is shown as a solid line (Eq. (6)). The curves with a colour gradient indicate possible values for Lyα/Hα and Hα/Hβ for a clumpy dust screen, where on average NC = {3, 7, 9} clumps block the sight lines (Eq. (9)). These lines are drawn with a colour gradient that indicates the dust optical depth τC at 1216 Å of a single clump (colour bar at the top of the plot).

Scarlata et al. (2009) point out that galaxies where the measured Lyα/Hα is higher than the expectation of Eq. (6), can be explained by obscuration by a system of dusty clumps, where on average N_C clumps obscure the lines of sight. For this geometrical scenario Natta & Panagia (1984) derive the effective dust optical depth

$$\begin{array}{c}\hfill {\tau }_{\mathrm{eff}}(\lambda )=N_C·(1-e^{-{\tau }_c(\lambda )})\text{,}\end{array}$$(9)

where τ_c(λ) is the dust optical depth of a single clump. In Fig. 14 we plot the curves from this scenario parameterised by τ_c(λ_Lyα) for N_C = {3, 5, 7}. Here we assume that the dust-extinction of the individual clumps obey also the Cardelli et al. (1989) law that was assumed for homogeneous screen. The notable feature of these curves is that they describe loops: With increasing dust optical depth the individual clumps start to significantly block emission over all involved lines and while the overall intensities decrease, with lim_{τc → ∞}τ_eff(λ) = N_C, the ratio starts to become dominated by the sightlines that are not obscured by the dusty clumps. As evident from Fig. 14 all 9 galaxies above the expectation from Eq. (6) show Lyα/Hα ratios that are on or close to the loci associated with clumpy dust. Our analysis thus implies a connection between the small scale distribution of dust and ionised gas kinematics as traced by ${\sigma }_0^{\text{obs}}$.

4.2.3. The influence of turbulence on f_esc^Lyα

We close by discussing the anticipated effects of turbulence on Lyα radiative transfer in the very idealised set-up of a homogeneous static medium (e.g. Neufeld 1990; Verhamme et al. 2006; Laursen et al. 2009a). This idealised scattering medium is characterised by the density of neutral hydrogen, n_H, or a corresponding neutral column, N_H = n_Hs, along the path length s. Without turbulence the velocities of the atoms follow a Maxwellian distribution with the corresponding thermal velocity dispersion σ_therm = 12.849 km s${}^{-1}\times \sqrt{T_4}$ for a given temperature T₄ in 10⁴ K. Turbulent motions, characterised by a velocity dispersion σ_turb, broaden the velocity dispersion. In the framework of micro-turbulence this broadening is expressed by a convolution of the velocity distribution of the thermal motions with the velocity distribution of the turbulent motions (e.g. Verhamme et al. 2006; Smith et al. 2022). The intrinsic velocity dispersion is thus given by ${\sigma }_0=\sqrt{{\sigma }_{\mathrm{therm}}^2+{\sigma }_{\mathrm{turb}}^2}$. Adopting Eq. (9) from Verhamme et al. (2006) we can express the monochromatic optical depth of this medium at the Lyα line centre as

$$\begin{array}{c}\hfill {\tau }_{{\lambda }_{\mathrm{Ly}\alpha }}=4.253\times {10}^{-13}\times \frac{N_{\mathrm{H}}[{\mathrm{cm}}^{-2}]}{{\sigma }_0[\mathrm{km}{\mathrm{s}}^{-1}]}.\end{array}$$(10)

Hence, a highly turbulent medium will naturally become less opaque for resonant absorption of Lyα photons. The reduced optical depth will decrease the number of scatterings before a Lyα photon can escape the medium, and thereby reduce the probability that it will be absorbed by dust.

In fact, there exists an analytically derived approximate solution for $f_{\text{esc}}^{\text{Ly}\alpha }$ for the homogeneous static medium, if we assume a monochromatic Lyα source at the centre of a plane-parallel slab and if we assume that this medium is mixed with dust of optical depth τ_d, Lyα (Neufeld 1990). Following Verhamme et al. (2006) this solution can be written as

$$\begin{array}{c}\hfill f_{\mathrm{esc}}^{\mathrm{Ly}\alpha }\approx 1/cosh\left(\frac{\sqrt{3}}{{\pi }^{5/12}\zeta }\left[{\left({(a{\tau }_{{\lambda }_{\mathrm{Ly}\alpha }})}^{1/3}{\tau }_{d,\mathrm{Ly}\alpha }\right)}^{1/2}\right]\right)\text{.}\end{array}$$(11)

Here a = A₂₁/4πΔν_D, with A₂₁ = 6.265 × 10⁸ s⁻¹ being the Einstein coefficient of the Lyα transition and Δν_D = (σ₀/c)×ν_Lyα being the Doppler frequency broadening due to the thermal and turbulent motions (ν_Lyα = 2.466 × 10¹⁵ Hz). Moreover, ζ in Eq. (11) is a fitting parameter that Neufeld (1990) derive to be ζ ≈ 0.525. Using Eq. (7) and Eq. (10) we can rewrite Eq. (11) as

$$\begin{array}{c}\hfill f_{\mathrm{esc}}^{\mathrm{Ly}\alpha }\approx 1/cosh(2.4\times {10}^{-2}\frac{{(N_{\mathrm{H}}[{\mathrm{cm}}^{-2}])}^{1/6}E{(B-V)}^{1/2}}{{({\sigma }_0[\mathrm{k}m{s}^{-1}])}^{2/6}})\text{.}\end{array}$$(12)

This analysis shows that already in a very idealised set-up no simple monotonic relation between $f_{\text{esc}}^{\text{Ly}\alpha }$ and ${\sigma }_0^{\text{obs}}$ is expected since additionally variations in neutral column and dust content will factor into the solution in a non-trivial manner. This result is in line with our finding that ${\sigma }_0^{\text{obs}}$ on its own does not affect $f_{\text{esc}}^{\text{Ly}\alpha }$ (Sect. 4.1.1), but that it emerges as an important parameter in concert with E(B − V) and two other parameters (F_cov and O₃₂, Sect. 4.1.2) that are related to the neutral column.

In Fig. 15 we show curves according to Eq. (12) in the $f_{\text{esc}}^{\text{Ly}\alpha }$ vs. σ₀ parameter space for $f_{\text{esc}}^{\text{Ly}\alpha }$, ${\sigma }_0^{\text{obs}}$, and E(B − V) values that are encountered in our sample. There we also plot the $f_{\text{esc}}^{\text{Ly}\alpha }$ vs. ${\sigma }_0^{\text{obs}}$ pairings of our sample colour coded by E(B − V). We see that the majority of our low-E(B − V) objects are commensurate with Eq. (12) for N_H ∼ 10¹⁹ cm⁻² and that high-E(B − V) objects would require lower neutral columns. However, most of these low-$f_{\text{esc}}^{\text{Ly}\alpha }$ high-E(B − V) objects show Lyα/Hα ratios indicative of a clumpy dust distribution (Sect. 4.2.2).

Fig. 15. Lyα escape fraction, $f_{\text{esc}}^{\text{Ly}\alpha }$, vs. intrinsic velocity dispersion σ0. Solid, dashed, and dotted lines indicate the analytically derived approximate solution, our Eq. (12), for monochromatic Lyα radiation escaping a plane parallel slab of NH = 1018 cm−2 (solid lines), 1019 cm−2 (dashed lines), and 1020 cm−2 (dotted line), respectively. For NH = 1018 cm−2 we show curves with E(B − V) = {0.1, 0.4, 0.7}, for NH = 1019 cm−2 curves with E(B − V) = {0.1, 0.4}, and for NH = 1020 cm−2 the curve for E(B − V) = 0.1. We also show the $f_{\text{esc}}^{\text{Ly}\alpha }$ vs. ${\sigma }_0^{\text{obs}}$ pairings of the discussed sample. Scatter points and curves are colour coded by E(B − V) as indicated by the colour bar on the right. As in Fig. 11 we highlight points by red outlines that show LLyα/LHα ratios above the maximum dust screen expectation from Eq. (6).

We caution not to over overinterpret the comparison to the here presented analytic approximation, as Eq. (12) requires a highly idealised set-up. Our values on ${\sigma }_0^{\text{obs}}$ trace the frequency distribution of the intrinsic Lyα radiation field, which is thus far from monochromatic, and the static plane-parallel slab oversimplifies the density and velocity distribution of the scattering medium encountered in reality. Moreover, Munirov & Kaurov (2023) recently highlighted the limited applicability of the micro-turbulence framework for Lyα radiative transfer. These authors show how the macroscopic nature of turbulence may significantly alter τ_Lyα and thus also the number of scatterings and the spectral distribution Lyα photons emerging from a turbulent medium. This is because turbulence leads to zones of coherent velocities whose size- and velocity distributions are statistically described. Hence the simple convolution of the thermal- and turbulent velocity distributions used in the micro-turbulence framework may not be justified in reality. Munirov & Kaurov (2023) show especially that macroscopic turbulence can decrease the number of scatterings by orders of magnitude before Lyα photons escape the medium. While this underpins the importance of turbulence in regulating $f_{\text{esc}}^{\text{Ly}\alpha }$, quantitative inferences will require observational estimates of the characteristic scale lengths of the turbulent flows.

5. Conclusions

Here we present an analysis of the ionised gas kinematics obtained from all 42 galaxies that comprise the combined LARS+eLARS samples. The whole sample was observed in medium spectral resolution with the PMAS integral field unit at the Calar-Alto 3.5 m telescope. We constructed line-of-sight velocity fields and velocity-dispersion maps for each galaxy (Sect. 3.1) and present the newly obtained kinematic maps for the 28 eLARS galaxies (Fig. B.1). We combined this sample with the 14 galaxies that were analysed in H16. To date, this sample is the largest kinematical sample of galaxies that are also observed in Lyα emission and, thanks to the comprehensive HST imaging analyses of the sample presented in M23, now allows more robust statistical analyses regarding the possible trends found in H16. The reduced datacubes and the analysed velocity fields are released to the community with this publication. Our aim here is to check whether the Lyα observables EW_Lyα, $f_{\text{esc}}^{\text{Ly}\alpha }$, and L_Lyα are affected by galaxy kinematics.

Our analysis of this sample revealed two null results. First, the Lyα observables are not affected by the morpho-kinematical nature of the velocity fields (Sect. 3.2.2). Second, the inclination of a system in our sample also does not affect the Lyα observables (Sect. 3.4). The first null result is new, and requires verification with larger samples at high redshift. The second result is in agreement with the null result found in photometric inclination estimates in high-z studies. Still, given the discrepancies between photometric and kinematic inclinations, verification with kinematical samples at high-z appears desirable.

We proceeded to quantify the global kinematical state of the galaxies via measurements of their intrinsic velocity dispersion, ${\sigma }_0^{\text{obs}}$ (corrected for beam-smearing effects; Appendix C), and their projected maximum velocity amplitude, v_shear (dubbed shearing velocity, Sect. 3.3). The gas kinematics of the sample are characterised by highly turbulent motions, with 30 km s⁻¹ ≲ σ₀^obs ≲ 80 km s⁻¹. Moreover, the majority of the shearing velocities in the sample are significantly below the values for disc galaxies. More than half of the sources of the sample (24 galaxies) show dispersion-dominated kinematics.

Our main result is that clear trends emerge between the kinematical measurements and the Lyα observables (Sect. 3.6). In particular, we find that EW_Lyα and L_Lyα correlate with ${\sigma }_0^{\text{obs}}$, whereas $f_{\text{esc}}^{\text{Ly}\alpha }$ anti-correlates with v_shear and v_shear/σ₀^obs. Moreover, adopting the canonical criterion of EW_Lyα ≥ 20 Å to define a LAE, 75 % of the LAEs in our sample are dispersion-dominated systems; that is, they are characterised by v_shear/σ₀^obs < 1.86. This result is in agreement with the recent study of high-z galaxies by Foran et al. (2024). We also find that ${\sigma }_0^{\text{obs}}$ strongly correlates with L_Lyα. However, this correlation is driven by the normalisation of L_Lyα with L_Hα in our sample. The trends uncovered here are commensurate with our initial hypothesis that the kinematics play an important role in regulating the Lyα observability of a galaxy.

We then tried to quantify the relative importance of the kinematical parameters with respect to other known observational and physical parameters of the galaxy sample. Our statistical analyses in Sect. 4.1.1 and Sect. 4.1.2 reveal that ${\sigma }_0^{\text{obs}}$ emerges as a dominant variable in regulating EW_Lyα. We also find that ${\sigma }_0^{\text{obs}}$ is of significant importance regarding $f_{\text{esc}}^{\text{Ly}\alpha }$ and L_Lyα. However, no clear trends regarding the importance of v_shear or v_shear/σ₀^obs emerged.

We note that the metrics used here to discuss the parameter importance are relatively simplistic. Moreover, our step-wise approach and also the marginal approach are actually criticised and disfavoured in the modern statistical literature (Grömping 2015; Smith 2018). Future analyses of the problem would benefit from more refined importance metrics, such as those reviewed by Grömping (2015). Adopting these methods is beyond the scope of the present analysis. Nevertheless, we suggest that the analysis of simple one-to-one correlations in particular might lead to erroneous conclusions or over-simplifications regarding the astrophysics of Lyα radiative transport.

While the empirical results are clear, we have not yet identified the dominating mechanism injecting turbulence and thereby regulating ${\sigma }_0^{\text{obs}}$. This complicates any interpretation of the fact that ${\sigma }_0^{\text{obs}}$ is found to be of utmost importance in regulating EW_Lyα. Motivated by results from stellar-population analysis, we discussed a simple scenario in which the age of the starburst simultaneously regulates ${\sigma }_0^{\text{obs}}$, $f_{\text{esc}}^{\text{Ly}\alpha }$, and EW_Lyα. However, this scenario does not agree with our data.

That ${\sigma }_0^{\text{obs}}$ is of dominating importance for regulating EW_Lyα while not being of similar importance for regulating $f_{\text{esc}}^{\text{Ly}\alpha }$ deserves a closer look, especially because of the highly significant linear relationship between EW_Lyα and $f_{\text{esc}}^{\text{Ly}\alpha }$. In Sect. 4.2, we discuss the findings that galaxies that do not follow this linear relationship are often better described by a clumpy dust screen. However, those galaxies were also found to exhibit high ${\sigma }_0^{\text{obs}}$. This suggests a connection between ionised gas kinematics and small-scale dust distribution.

Lastly, we discuss how, even in a highly idealised medium, a complicated interplay between N_H, E(B − V), and σ₀ is expected to influence $f_{\text{esc}}^{\text{Ly}\alpha }$. We also highlight that the conventional treatment of turbulence as a microscopic effect may lead to underestimation of the actual quantitative influence of σ₀ on $f_{\text{esc}}^{\text{Ly}\alpha }$. The results presented here suggest that further insight could be gained by investigating the trends we have uncovered in future analyses from a modelling and radiative transfer perspective, especially the role of turbulence in facilitating Lyα escape. Moreover, future spatially resolved kinematical studies of high-z galaxies with and without Lyα in emission could help us to understand whether or not our results are indeed applicable to the early Universe.

Data availability

The fully reduced PMAS cuboids are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/693/A252

Acknowledgments

Based on observations collected at the Centro Astronómico Hispano-Alemán (CAHA) at Calar Alto, operated jointly by Junta de Andalucía and Consejo Superior de Investigaciones Científicas (IAA-CSIC). The authors acknowledge financial support by the Centre National D’Etudes Spatiales (CNES). We acknowledge the use of the following software: ds9 (Joye & Mandel 2003), QFitsView (Ott 2012), Astropy (Astropy Collaboration 2022), Numpy (Harris et al. 2020), SciPy (Virtanen et al. 2020), matplotlib Hunter (2007), seaborn (Waskom 2021), cmasher (van der Velden 2020), GalPak^3D (Bouché et al. 2015), dust_extinction (Gordon 2024), and p3d (Sandin et al. 2010). MH is fellow of the Knut & Alice Wallenberg Foundation. AS acknowledges the support by Jörg Main and the hospitality at Institute for Theoretical Physics 1 (Stuttgart University) to carry out this research. AS’s contribution to this project was partly made possible by funding from the Carl-Zeiss-Stiftung. We thank especially Daniel Kunth for organising meetings at Institut Astrophysique de Paris and for continuous support with the project. ECH and AS acknowledge travel funds provided by the IAP for attending these meetings. We thank N. Bouché for support regarding the use of GalPak^3D. We thank Rikke Saust as well as the motivated staff at Calar-Alto observatory for help with the preparation and execution of the visitor mode observations. ECH acknowledges the generous hospitality of Leiden University for carrying out parts of this research. ECH acknowledges funding of an internship for AS within the DAAD RISE scheme. ECH dedicates this work in loving memory to a cat named Pixel.

References

Appendix A: PMAS data reduction

The raw observational data (Table A.1) were reduced with p3d⁹ (Sandin et al. 2010, 2012) with post-processing by custom developed Python routines. We outline the steps of the data reduction in the following, noting that a more comprehensive description can be found in H16.

First, we create a master bias from the individual bias frames for each night. Next, we create trace masks and dispersion masks for all continuum- and arc-lamp exposures. The trace masks are then used to extract 1D spectra for each spectral trace of the target exposures. For this extraction of the spectral traces from the detector we use the optimal weighting scheme from Horne (1986). For the spectral re-binning in the extraction we opt for the “Drizzle” algorithm by Fruchter & Hook (2002). The dispersion masks that control the relevant mapping from detector pixel to wavelength space during the extraction are found by fitting 4^th order polynomials to the HgNe lines from the arc-lamp exposures. Cosmic-rays hits on the detector are rejected with the L.A.Cosmic algorithm of van Dokkum (2001) using the same parameters as in H16. In the final step the extracted data are flat-fielded with the sky-flats that were extracted in the same manner as the science data.

Different to H16 here we do not flux calibrate the eLARS PMAS data, since an absolute flux calibration is not needed for the kinematic analysis. For the same reason we also do not subtract telluric emission from the extracted spectra, except for galaxies where the redshifted Hα line is in close proximity or even overlapping with a sky line (eLARS 2, and eLARS 22 to 28).

The data for targets where open-shutter time was accumulated over multiple exposures is stacked via simple addition. The resulting data products from above reduction chain are a row-stacked representation of the 256 spaxels for each science target exposure in analogue-digital-units. p3d also propagates errors through each step of the reduction chain and stores them in a separate row-stacked spectra array. These arrays were then reformatted to 3D cuboids using the lookup table that comes with p3d. The unvignetted spectral range covers ∼ 6100 Å to 7400 Å, and the spectral sampling is 0.46 Å per wavelength bin.

We registered the world coordinate system (WCS; Greisen & Calabretta 2002) of all eLARS PMAS cuboids against the HST Hα images from Melinder et al. (2023). We did this visually¹⁰, taking inspiration from the alignment procedure outlined in Ferruit (2010). For the visual alignment we relied on two images. First, we created surface brightness contours that outline important morphological features visible in Hα. These contours were overlaid over a Hα narrow-band image that was extracted from the PMAS datacube at a given field of view position. Secondly, we created a simulated PMAS Hα flux map from the HST Hα image at a given field of view position. The field of view position was then adjusted manually until the contours and the simulated observations matched satisfactorily with the observations. The final field of view positions with respect to the HST imaging data are shown in Appendix A.1 of Schaible (2023). We estimate that the relative offset between HST imaging WCS and the so determined PMAS cuboid WCS is better than 0.25″.

For the analysis in the present article we only use this WCS information to resample the datacuboids of the three galaxies that required two pointings onto a common grid (eLARS 3, 5, and 26). This is achieved with the reproject python package (Robitaille et al. 2020). However, the presence of a WCS in the PMAS datacuboids of LARS and eLARS, which is relatively accurate with respect to the HST imaging data, may be exploited in future analysis that aim at investigating local relations between kinematics and Lyα observables.

Appendix B: Kinematic map creation

In order to create the kinematic maps shown in Fig. B.1 we first subtract a running median in spectral direction (width = 151 wavelength bins) from each spaxel of the datacube. This removes telluric- and stellar continuum and only the emission line profiles from the galaxies remain. Next, we model the Hα profile in each spaxel with a simple Gaussian profile

$$\begin{array}{c}\hfill f_i(\lambda |A_i,{\lambda }_{0,i},{\sigma }_i)=A_{i}exp(-{[(\lambda -{\lambda }_{0,i})/{\sigma }_i]}^2).\end{array}$$(B.1)

Fig. B.1. Hα kinematics of eLARS from PMAS (16″×16″ field of view). For each galaxy we show HαS/N (left panels; colour coded in log-scale from 1 to 1000), vlos [km s−1] (centre panels; approaching to receding velocities are colour coded linearly from blue to red with the centre at vlos = 0 in white and scaled symmetrically to the absolute maximum), and σv [km s−1] (right panels; colour coded linearly from the 2th-percentile, median, and 98-th percentile of the observed σv distribution per galaxy). The name of each galaxy is coloured according to the visual kinematic classification (Sect. 3.2): rotating discs in purple, perturbed rotators in green, and galaxies with complex kinematics in orange. Galaxies with two PMAS pointings are assembled at the bottom of this two-page figure.

Fig. B.2. continued.

The best-fitting A_i, λ_0, i, and σ_i for each spaxel i are found with the Levenberg-Marquardt algorithm. For the fitting we use the framework provided by astropy.modelling (Astropy Collaboration 2022). The initial parameter guesses for the algorithm are derived from the data itself – for λ_0, i and σ_i via the first- and second moment and for A_i via the peak value of the observed line profile. For each spaxel the best-fit λ_0, i and σ_i are then converted into v_los, i and σ_obs, i in km s⁻¹ via

$$\begin{array}{c}\hfill v_{\mathrm{los},i}=c({\lambda }_{0,i}-\overline{\lambda })/{\lambda }_{0,i}\end{array}$$(B.2)

and

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{obs},i}=c\sqrt{({\sigma }_i^2-{\sigma }_{\mathrm{LSF},i}^2)/{\lambda }_{0,i}^2},\end{array}$$(B.3)

where $\overline{\lambda }$ is the average of all reliable λ_0, i fits and ${\sigma }_{\mathrm{LSF},i}={\mathrm{FWHM}}_{\mathrm{LSF},i}/(2\sqrt{2ln2})$ relates to resolving power R_i = λ_0, i/FWHM_LSF, i at λ_0, i for each spaxel i. We note that cosmic ray hits severely hampered the line fits for some galaxies that were observed only with a single exposure. We mask the affected spaxels manually. This results in a few gaps in the velocity field (see, e.g. eLARS 4 and eLARS 11 in Fig. B.1 for the worst cases).