© The Authors 2025

1. Introduction

The evolution of galaxies is impacted by the gas reservoirs around them, known as the circumgalactic medium (CGM; Tumlinson et al. 2017). The CGM plays a central role in the exchange of gas, dust, and metals between the galaxy and its surroundings. It funnels toward galaxies the gas needed to form stars, and it is the site where powerful galactic outflows end up and where gas recycling happens (Putman et al. 2012; Fox & Davé 2017; Péroux & Howk 2020; Veilleux et al. 2020). The CGM presents a complex density and kinematic structure, which has been studied in the past decades thanks to multi-wavelength observations and simulations. They have revealed that the CGM is “multiphase”, as it consists of different components ranging over wide intervals of density, temperature, and ionization state (e.g., Péroux et al. 2019; Wakker et al. 2012; Ford et al. 2013; Anderson et al. 2013; Suresh et al. 2017; Weng et al. 2022). This complexity makes the modeling of the geometry and kinematics of the CGM challenging, in particular at high redshift, where it is usually modeled with simplified spherically or cylindrically symmetric static or expanding gas geometries (such as expanding ellipsoids) or with bipolar outflows (e.g., Eide et al. 2018), but the currently available data do not allow us to distinguish between them.

The Lyman-α (Lyα) line is one of the best observational signatures used in such studies, in particular at high-z, given its dependence on the structure, ionization, and kinematics of the H I gas where its photons propagate (Osterbrock 1962; Dijkstra et al. 2016; Gronke & Dijkstra 2016). It is the strongest tracer of recombining ionized hydrogen in young star-forming galaxies (Partridge & Peebles 1967) and is ubiquitously detected at high-z (e.g., Zitrin et al. 2015; Vanzella et al. 2017; Caminha et al. 2023; Bunker et al. 2023; Nakane et al. 2024), but its interpretation is still debated. In fact, different from other lines such as the Hα, whose photons propagate undisturbed to reach us, Lyα has a resonant nature, and thus a Lyα photon can undergo a great number of scatterings after its emission. The number of scatterings that it experiences before being able to leave its emission site depends on the H I column density, geometry, and kinematics (Adams 1972; Dijkstra 2014); on quantum mechanical probabilities (Stenflo 1980); and on the properties of the region where it originated. For instance, centrally emitted Lyα photons, for example, created as nebular emission powered by star formation, significantly scatter before escaping, potentially giving rise to an observed spatially extended Lyα emission. However, spatially extended Lyα emission can also be produced by cooling gas (Haiman et al. 2000), gas that has been shock-heated by supernova explosions (Mori et al. 2004) and galactic winds (Taniguchi & Shioya 2000), fluorescent radiation from an external ionizing field (Hogan & Weymann 1987; Cantalupo et al. 2005), or extended star formation (Momose et al. 2016; Mas-Ribas et al. 2017). These features are encoded in the intensity spectrum of the source (Ahn et al. 2002, 2003; Verhamme et al. 2006; Dijkstra & Loeb 2008; Gronke et al. 2015), causing the broadening and shifting of its Lyα line profile (Neufeld 1990; Dijkstra et al. 2006). Moreover, the Lyα spectrum is also affected by radiative transfer effects at the interstellar medium (ISM) and CGM scales and is dependent on the inclination and evolutionary stage of the galaxy (see Blaizot et al. 2023, for a study on a simulated Lyα emitter galaxy at z ∼ 3 − 4). The spectrum and the Lyα surface brightness profile, which reveals the spatial distribution of the Lyα emission and the diffusion process of Lyα photons, are the observables most frequently used to investigate the nature of the Lyα emission. But thanks to its resonant nature, we can additionally leverage the Lyα degree and direction of linear polarization typically represented by the Stokes parameters Q and U.

The Lyα degree of polarization increases if photons are scattered in a preferential direction, and the resulting value mainly depends on two fundamental factors: the production mechanism and site where the Lyα photons are created, and the geometry of the gas where they scatter before escaping, as it is particularly sensitive to the isotropy and homogeneity of the emission and gas distribution (e.g., Lee & Ahn 1998; Ahn et al. 2002; Eide et al. 2018). Theoretical studies have shown that different models may present similar spectra but different degrees of polarization of the Lyα line (e.g., Dijkstra & Loeb 2008; Gronke et al. 2015; Eide et al. 2018).

Due to the observational difficulty of measuring the Lyα degree of polarization of distant sources and the limited number of available (spectro-) polarimeters intended for extragalactic use (Hayes & Scarlata 2011; Beck et al. 2016), only a few studies that map the polarization of Lyα are currently present, and they have mostly made use of the narrow-band imaging mode. Moreover, they have all targeted bright (L_Lyα > 10⁴³ erg s⁻¹) and extended (up to ∼150 kpc) Lyα emissions at z ∼ 2 − 3 in dense environments, with no clear central powering sources, the so-called Lyα blobs (LAB) or around extreme environments, such as high-z quasars, overdense regions (clusters or protoclusters), and Active Galactic Nuclei (AGN) or radio galaxies. The interpretation of the results remains questionable due to a lack of general consensus. Koratkar et al. (1995) made use of the Hubble Space Telescope (HST) Faint Object Spectrograph (FOS) spectropolarimeter to observe three quasars at z = 0.5 − 1.6, including PG 1630+377, which presents an increase of the degree of polarization up to ∼20% blueward the Lyman break and a Lyα line polarization of (7.3 ± 1.6)% that can be explained with the presence of two sources, one completely obscured and producing the scattered polarized signal and a secondary unpolarized redder source. Vernet & Cimatti (2001) presented low resolution spectropolarimetry of the submillimeter-selected galaxy SMM J02399−0136, finding a Lyα degree of polarization of $2.1_{-0.5}^{+0.9}\%$. Vernet et al. (2001) measured Lyα degrees of polarization < 2% in nine radio galaxies at z ∼ 2.5 with Keck II/Low Resolution Imaging Spectrometer (LRIS). Prescott et al. (2011) observed the LAB LABd05, containing an obscured AGN at z = 2.656. Due to the coarse spatial resolution, they could only put an upper limit of ≈5% on the Lyα degree of polarization within an aperture of 65 kpc. With deeper data, Kim et al. (2020) found a consistent polarization fraction of (6.2 ± 0.9)% within the same aperture, but could also detect a spatially resolved polarization varying from 5% at the Lyα peak to 20% 45 kpc away, consistent with Lyα photons not being scattered in the central region between the AGN and the Lyα peak, but only in the outer gas surrounding the nebula. Similarly, Humphrey et al. (2013) measured a low (< 5%) polarization fraction in the center of the LAB hosting the TXS 0211−122 radio galaxy at z = 2.34, increasing to (16.4 ± 4.6)% in the eastern section. You et al. (2017) found no polarized signal in the center of the LAB B3, which surrounds a radio-loud AGN at z = 3.09, and a degree of polarization of 3% (17%) at 10 (25) kpc, with an asymmetric distribution. Hayes et al. (2011) observed LAB1, in the SSA22 protocluster at z = 3.09 and did not detect a polarized signal in the center, but a 20%-polarized ring at approximately 45 kpc. The results were confirmed by using spectropolarimetric observations by Beck et al. (2016), who found increasing polarization toward the wings of the Lyα spectral profile, which can be explained by the presence of outflows. Making use of integral field spectroscopy observations for this system, Herenz et al. (2020) found that regions with a larger degree of polarization also have high velocity shifts and narrow line profiles, and they associated this evidence with Lyα scattering from a central source. Finally, North et al. (2024) did not detect a polarized signal in the LAB hosting the radio-quiet quasar SDSS J1240+145, at z = 3.11.

Additionally, state-of-the-art models and simulations are needed to interpret the observations. In the last few years, significative steps forward have been done in such studies, with the implementation of Lyα polarization in advanced Lyα radiative transfer codes (e.g., Ahn et al. 2000; Dijkstra & Loeb 2008; Ahn & Lee 2015; Chang et al. 2017, 2023; Eide et al. 2018; Seon et al. 2022; Chang & Gronke 2024), that allow us to predict the polarization behavior in different geometries (spherically symmetric, as with the expanding shell, and non spherically symmetric, as with the expanding ellipsoids or biconical outflows) and with different physical properties (density and clumpiness of the gas) and kinematics.

In this paper we investigate the Lyα origin and the geometry of the scattering CGM through the first study of the degree of polarization of the Lyα emission line for a typical clumpy, star-forming galaxy (Abell 2895a in Livermore et al. 2015) at z ∼ 3.4, strongly lensed by the cluster of galaxies Abell 2895 (A2895) into three multiple images. Thanks to the image multiplicity and the lensing magnification, we are able to study the properties of this source in great detail using the available multi-wavelength dataset, which includes ancillary observations taken with HST, the Multi Unit Spectroscopic Explorer (MUSE), the Enhanced Resolution Imager and Spectrograph (ERIS), and the Spectrograph for INtegral Field Observations in the Near Infrared (SINFONI) at the Very Large Telescope (VLT), the Atacama Large Millimeter/submillimeter Array (ALMA), and new observations taken with VLT/FOcal Reducer/low dispersion Spectrograph 2 (FORS2) instrument, with its Polarimetric Multi Object Spectroscopy (PMOS) mode.

The paper is organized as follows. In Section 2, we present the data available for this system, summarizing the ancillary archival data and then focusing on the new FORS2 PMOS observations. We also describe the pipeline and the step followed in the data reduction. In Section 3, we describe the 1D spectra extraction, the dilution correction, the measurement of the Stokes parameters and polarization fraction as a function of the wavelength, and the assumed bins in wavelength. In Section 4 we expose the radiative transfer models we developed to interpret the results, and we compare observations and models in Section 4.3. In Section 5 we discuss our results and the future perspectives. In Section 6, we summarize our results.

Throughout this paper, we adopt a flat ΛCDM cosmology with Ω_Λ = 0.7, Ω_m = 0.3, and H₀ = 70 km s⁻¹ Mpc⁻¹. We report all the measurements corrected for lensing effects.

2. Data

2.1. Our target galaxy and its ancillary data

Our target, Abell 2895a, is a clumpy star-forming (SFR = 10 ± 0.3 M_⊙ yr⁻¹) galaxy at z ∼ 3.4, characterized by an extended Lyα emission with a total flux of (1.41 ± 0.04)×10⁻¹⁷ erg s⁻¹ cm⁻², offset by 1.2 ± 0.2 kpc with respect to the UV continuum (Iani et al. 2021; Zanella et al. 2024). Abell 2895a is strongly lensed by the galaxy cluster A2895 into three multiple images with coordinates (M1; M2; M3) (RA, Dec) = (01:18:11.19, −26:58:04.4; 01:18:10.89, −26:58:07.5; 01:18:10.57, −26:58:20.5) (Livermore et al. 2015; Iani et al. 2021). They are located in the inner region of the A2895 cluster, angularly close to the brightest cluster galaxy (BCG; at z = 0.227). In this work, we focus on the two most magnified images, M1 and M2, whose average lensing magnifications are μ = 5.5 ± 0.7 and μ = 4.5 ± 0.3, respectively (Iani et al. 2021).

Abell 2895a has a redshift of z_opt = 3.39535 ± 0.00025, estimated from optical emission lines (Iani et al. 2021), which is consistent with the z_[C II] = 3.39548 ± 0.00007 value estimated from the [C II] far-infrared line (Zanella et al. 2024). It presents a clumpy morphology in the UV continuum and in [C II] emission line observations. At least four star-forming clumps within a diffuse emission are detected in the UV continuum, one of which is also observed in the [C II] data. Two additional clumps, without a detected UV counterpart, are identified in [C II], whereas no dust continuum is detected down to F_cont < 34 μJy (Zanella et al. 2024), in agreement with the measured blue UV-continuum β slope of −2.53 ± 0.15 and low reddening of E(B − V) < 0.16 mag (Iani et al. 2021).

Abell 2895a offers a suite of ancillary data, presented in Iani et al. (2021) and Zanella et al. (2024). The rest-frame UV imaging was observed with HST (SNAP program 10881, PI: G. Smith) with a FWHM ∼ 0.13″ resolution in the F606W filter. The Hβ and [O III] lines have been studied with VLT/SINFONI spectra (Livermore et al. 2015), and new VLT/ERIS observations (program IDs 110.2576, 112.25HA, 114.273Y, PI: A. Zanella), aimed at spatially-resolving them, are ongoing. The Lyα emission was targeted by VLT/MUSE observations (program IDs 60.A-9195(A), 0102.B-0741(A), PI: A. Zanella) in the Adaptive Optics (AO) Wide Field Mode (WFM), with a resolution of FWHM ∼ 0.4″. It appears spatially offset with respect to the clumps visible in the UV continuum and the optical emission lines. The HST rest-frame UV and the Lyα contours are shown in Fig. 1. The Lyα line also presents an asymmetric profile that is redshifted with a relative velocity Δv = 403 ± 4 km s⁻¹ with respect to the systemic redshift (Iani et al. 2021). Finally, ALMA Band 8 observed Abell 2895a (program ID 2019.1.01676.S, PI: E. Iani) to detect the [C II] emission line and the underlying continuum (Zanella et al. 2024). The beam size (assuming a natural weighting) is FWHM = 0.31″ × 0.26″.

Fig. 1. Left: HST F606W image of the inner part of A2895, where the three multiple images of Abell 2895a (M1, M2, and M3) appear. It represents the rest-frame UV at the redshift of Abell 2895a, z ∼ 3.4. In orange, we show the 2, 3, and 5σ contours of the Lyα emission, detected with MUSE. The green box represents the FORS2 1.4″ × 22″ adopted slit, that includes M1 and M2, and one image of another source at z ∼ 3.7 (Iani et al. 2023), whose Lyα 2, 3, and 5σ contours are shown in red. Thanks to the PMOS mode, the signal included in the slit is split in an ordinary (o) and extraordinary (e) ray, with orthogonal polarizations, shown with different tones of green. Right: Stacked 2D spectrum with φ = 0.0° obtained after calibration, cosmic ray rejection, and sky subtraction. In particular, the sky is subtracted only in the relevant slits which contain Abell 2895a, and appear darker. Each slit has a spatial (vertical) width of 22″, and the sky level is estimated and then subtracted from the closest ordinary and extraordinary slits (see the e and o labels on the right), respectively, as shown by the green arrows on the right. In the slits including the target, we detect the extended Lyα emission associated with M1 and M2, at approximately 5340 Å, the continuum from the nearby BCG, in the bottom part, and the Lyα emission of the source at z ∼ 3.7 (red squares, close to the edge of the slitlets). We show the 2D spectrum taken with φ = 0.0° on top, and smaller cutouts around the Lyα line for φ = 22.5°, 45.0° and 67.5° on the bottom.

A2895 also benefits from a robust strong lensing model, introduced by Iani et al. (2021). The 2D-projected total mass distribution of the cluster is modeled as a combination of an extended cluster-scale halo and multiple galaxy-scale double pseudo-isothermal elliptical components (dPIE; Elíasdóttir et al. 2007), whose centers and shapes are constrained by the respective surface brightness centroids, ellipticities, and position angles from the HST F606W image. The cluster members are selected through the color-magnitude diagram method (e.g., Richard et al. 2014), and the total mass associated to each member is computed from its luminosity, through the Faber-Jackson relation for elliptical galaxies (Faber & Jackson 1976), as is usual in strong lensing modeling on cluster scales (e.g., Caminha et al. 2023; Bergamini et al. 2023). The model is constrained by using the locations of the multiple images of Abell 2895a and those of another triply imaged system with spectroscopic redshift z = 3.721 (Livermore et al. 2015; Iani et al. 2023), shown in Fig. 1. The best-fit model reproduces the location of the multiple images with a root-mean-square (rms) displacement of 0.09″ between the observed and modeled multiple images positions. The convergence maps derived from this strong lensing model are exploited to estimate the magnification factors in the locations of M1 and M2, which we use to derive the magnification-corrected quantities adopted in the following.

2.2. FORS2 PMOS observations and data reduction

Abell2895a was observed with VLT/FORS2 (Appenzeller et al. 1998) between September 2021 and August 2022 (program ID 108.2260, PI: A. Zanella), for a total of 18.1 hours in PMOS mode. The polarization optics are composed of a superachromatic half-wave plate mosaic followed by a Wollaston prism, which separates the light into two beams with orthogonal polarization (the “ordinary” (o) and “extraordinary” (e) rays). Half of the MOS mask slitlets in front of the polarization optics are fully closed to avoid the overlap of the o and the e beams, leaving eight 22″ high slitlets for science targets.

Observations were executed with seeing < 0.9″, clear sky conditions, fraction of lunar illumination < 0.4, and airmass ≲1.6. The run was divided into twelve sets of four 1200 s exposures with the half wave plate position angles (φ) set successively to 0°, 22.5°, 45°, and 67.5°. The target acquisition was performed, with a typical precision of < 0.1″, through a blind offset from a bright star, distant ∼40″ from the target. We used the MIT red CCD together with the 1400V grism and a slit width of 1.4″ for all observations providing an effective spectral resolution of about 3.6 Å FWHM covering the wavelength range from 4560 to 5860 Å. The slit was oriented at 40° North to East so that both M1 and M2 fit in a single slitlet (see Fig. 1), which also includes one multiple image of another source at z ∼ 3.7 (Iani et al. 2023). However, this galaxy is too close to the edge of the slit, and the slit losses are too important to analyze the polarization of this second target too, so it is not considered in the following analysis.

We reduced the data with the standard FORS2 PMOS pipeline v5.14, making use of the ESO Recipe Execution Tool (EsoRex; ESO CPL Development Team 2015) pipeline. We reduced separately the observations of each of the twelve OBs and combined them as the last step. We focused on the data taken with the CHIP1 (Norma) CCD, which contains the spectrum of Abell 2895a in the bottom part. We ran the calibration recipe to correct all the raw exposures using the associated BIAS frames, identify the slitlets limits, the dispersion relation, and the spatial distortion, and correct for the FLAT fields. These products were used as inputs in the science recipe, that produces as output wavelength-calibrated optical distortion-corrected 2D spectra. Given the faintness of our target in the single OB, we disabled the sky subtraction automatically performed by the pipeline, as it may affect the resulting signal-to-noise ratio (S/N) of the target. We detected and rejected the cosmic ray traces with the Astro-SCRAPPY (McCully et al. 2018) Python package, based on L.A. Cosmic (van Dokkum 2001). To properly subtract the sky background contribution for the o and e beams, we estimated the median flux at each wavelength in the two respective closest slits, located 22″ apart (Fig. 1, right panel, arrows on the right). We did not use the slits containing the target themselves to estimate the sky as they are dominated by the emission of M1, M2, and the BCG contribution. We checked that, after subtracting the median sky flux from the target slits, the residuals did not show systematics or gradients. Three OBs presented a strongly polarized background contamination, likely due to the presence of the moon, and thus with a not robust sky subtraction. We decided to exclude them so as not to bias our analysis. These OBs were the only ones taken with the moon above the horizon, and had the lowest angular separation between the target and the Moon (∼90° instead of the 140° −150° of the rest of the OBs). We stacked the remaining nine OBs by spatially matching the position of the Lyα peak. We computed, for each OB, a profile in the spatial y-direction around the Lyα line and identified the peaks of the two Lyα glows, relative to M1 and M2. We noticed that, in the selected OBs, the peaks do not show significant shifts (≲0.5″, much smaller than the extraction aperture used for 1D spectra in the following), and thus we directly stacked them, obtaining four 2D spectra, one for each of the four φ values (shown in Fig. 1). We measured the 2D variance spectra by converting the observed spectra from ADU/s to total counts, and propagating the uncertainties associated to the object, the sky, and the readout noise assuming Poissonian statistics.

The four resulting 2D spectra of the o and e channels including Abell 2895a consist of 22″ slits, with a spatial pixel sampling of ∼0.25″ pix⁻¹ and spanning in wavelength from 4560 Å to 5860 Å, with an effective spectral resolution of about 3.6 Å FWHM, and with a 0.64 Å pix⁻¹ dispersion.

3. Analysis

3.1. Intensity spectrum of Abell 2895a and dilution correction

To maximize the S/N, we extracted the 1D spectra by summing, for each wavelength, the signal included in an aperture of 7.5″, designed to encompass both the M1 and M2 Lyα emission. After ray-tracing it to the source plane, this aperture corresponds to a physical size of approximately 10 kpc at z = 3.4. The eight resulting 1D spectra, with intensities $I_{\phi }^{o,e}(\lambda )$, relative to the o and e channels with φ = 0.0° ,22.5° ,45.0°, and 67.5°, showed a significant contribution from the flux of the angularly close BCG. This signal is not polarized and it does not directly affect the polarization measurements of the Lyα line of Abell 2895a, but dilutes it. Thus, the maximum polarization fraction we could measure was limited by the dilution factor, f_d. We first estimated the total intensity spectrum I(λ) within the extraction aperture, which is equivalent to the total intensity Stokes parameter spectrum, as the sum of the eight 1D spectra. Then, we evaluated the total intensity spectrum associated with the BCG, I_BCG(λ), by applying the same procedure to eight spectra extracted within 4″-apertures centered on the bottom part of each slit, where there is no contribution from the Lyα blobs. We normalized the I(λ) and I_BCG(λ) spectra, and smoothed I_BCG(λ) to mitigate noise features. The I(λ) spectrum, before and after subtracting I_BCG(λ), is shown in the left panels of Fig. 2. We derived the dilution factor, defined as the fraction between the total intensity spectrum from the underlying continuum, mainly due to the BCG, and that of Abell 2895a only, as

Fig. 2. Left panels: Total intensity spectra (green solid line) and 1σ uncertainties (light green shaded region) in the spectral region around the Lyα line. The top (bottom) panel shows the normalized spectrum before (after) the subtraction of the BCG contribution. The black filled squares, with 1σ uncertainties, represent the binned data. Right panels: Polarization (P) measurements obtained before (top) and after (bottom) the dilution correction described in Eq. (1). The black open circles represent the 1σ upper limits obtained by applying the correction of Simmons & Stewart (1985). In the bottom right panel, only one datapoint is visible, as the others are outside the plotted range (upper limits for P ∼ 80%−95%).

$$\begin{array}{c}\hfill f_{\mathrm{d}}(\lambda )=\frac{I_{\mathrm{BCG}}(\lambda )}{I(\lambda )-I_{\mathrm{BCG}}(\lambda )}·\end{array}$$(1)

We measured f_d(λ) values ranging from approximately 2 at the peak of the Lyα to 10 in the tails.

3.2. The reduced Stokes parameters q, u, and P spectra and variance spectra

We combined the measured $I_{\phi }^o(\lambda )$ and $I_{\phi }^e(\lambda )$ intensities to measure the reduced Stokes parameters q and u¹ as a function of the wavelength. Combining observations taken with four different position angles is crucial to correctly handle the different gain factors in the o and e channels, and to extract reliable q(λ) and u(λ) spectra (Cohen et al. 1997). The degree of polarization is then measured as

$$\begin{array}{c}\hfill P(\lambda )=\sqrt{q^2(\lambda )+u^2(\lambda )},\end{array}$$(2)

and the uncertainties on q(λ), u(λ), and P(λ) are estimated through the propagation of the uncertainties associated with each of the eight I_φ^o(λ) and I_φ^e(λ) 1D intensity spectra. In principle, q(λ) and u(λ) allow us to also measure the polarization angle, but given the low S/N regime, we could not obtain significant measurements, and we do not consider the polarization angle in the following.

The degree of polarization measured through Eq. (2) is, by definition, a positive quantity. In low S/N regimes, the uncertainties on q(λ) and u(λ) lead to an increase of the biased measured value of P(λ), that will differ from the true unbiased value P₀(λ). We applied the correction by Simmons & Stewart (1985) that proposes, for different S/N regimes, four different possible methods to estimate P₀(λ) (the average estimator from Serkowski 1958, the Wardle & Kronberg 1974 estimator, the maximum likelihood estimator, and the median estimator) and their uncertainties. In high S/N regimes all the four methods predict consistent results, while they are particularly effective in estimating P₀(λ) in low S/N regimes, such as those presented in this work. In this regime, the quality of the measured q(λ), u(λ), and P(λ) can be moreover enhanced by binning the spectra in a number of bins that depends on the target S/N and on the spectral resolution needed for a proper interpretation of the data. We adopted two different approaches: we included the entire Lyα line in a single bin, integrating from 1215.5 Å to 1219.6 Å, or we divided it into three bins, one including the blue tail (1215.5−1216.3 Å), a central one including the peak (1216.3−1217.9 Å) and one including the red tail (1217.9−1219.6 Å). We also included two bins to sample the continuum blueward (one narrower and closer to the Lyα, from 1210.0 Å to 1215.5 Å, and one broader, from 1037 Å to 1210 Å) and redward (similarly, from 1219.6 Å to 1224.0 Å, and from 1224 Å to 1260 Å) the Lyα emission. These are partially visible in Fig. 2 as the horizontal black lines. In the following, we will refer to the case with a single bin for the Lyα, while the results for the other case are shown in Appendix A, as they are equivalent and bring to the same conclusions.

We applied Eq. (2), and measured low P values, consistent with zero. After the positive bias correction and the binning, we obtained observational upper limits on P₀ in the considered bins. We multiplied these upper limits by the dilution factor f_d(λ), correspondingly binned within the same chosen spectral windows (e.g., f_d ≈ 3 in the Lyα bin). The measured polarization fractions before and after the dilution correction are shown in the right panels of Fig. 2. The P₀ values revealed that we can put tighter constraints on the polarization fraction at the peak of the Lyα line, where the S/N is at its maximum, and increasingly shallower ones moving toward the tails. Far from the line, in the continuum where I(λ)≈I_BCG(λ), it is not possible to put informative constraints. After the dilution correction, we measure for the bin including the Lyα 1σ, 2σ and 3σ unbiased upper limits on the degree of polarization, P_Lyα, of 4.6%, 5.8%, 6.5%, respectively. Due to the large dilution factor of approximately 35 (100), the degree of polarization is barely (not) constrained for the narrow (broad) continua bins.

We checked whether the P_Lyα measurements could be affected by our choice of the aperture adopted to extract the spectra, including two multiple images of Abell 2895a. In fact, the strong lensing critical lines at the redshift of Abell 2895a pass between the two images, that appear reasonably mirrored (as can be seen from the UV morphology and the Lyα contours on the left panel of Fig. 1). We extracted the spectra separately from the extended Lyα of M1 and M2, visible on the right panel of Fig. 1, with apertures of 4″ centered on the peaks (against the 7.5″ aperture that includes both peaks). By adopting the same Lyα bin and corrections, we obtained consistent 1σP_Lyα upper limits of 7.4% and 5.8%, which are less stringent due to the lower S/N.

Both in the case of including the Lyα emission from the M1 and M2 multiple images together or separately, the adopted extraction apertures globally contain the Lyα emission from the entire galaxy. It has been observed that this approach can lower the measured degree of polarization, as the result of the cancellation of opposite contributions from different sides of the emission (e.g., Humphrey et al. 2013; You et al. 2017). In order to check whether this effect could have affected the low polarization we measured, we extracted the spectra from two regions, including separately the two different spatial halves of the Lyα emission. The results of this test are described in detail in Appendix B. We found, for the two different regions, upper limits on P_Lyα of 5.1% and 7.5% at the 1σ level, consistent with those found in our reference case. The low S/N did not allow us to include the polarization spatial information in our analysis and comparison with the models, which we discuss in Section 5.2.

4. Radiative transfer models for Lyα

Iani et al. (2021) modeled the spectral profile of Lyα through Lyα radiative transfer modeling (Dijkstra 2019), using the fitting pipeline of Gronke et al. (2015), adopting the shell model (e.g., Ahn et al. 2000; Verhamme et al. 2006), which is composed of a thin H I shell with a single constant radial velocity and a central Lyα source. However, due to the symmetry of the shell model, the integrated Lyα spectrum is always unpolarized. If the scattering medium is not symmetric, such as a bipolar wind or ellipsoidal halo, the integrated Lyα can be polarized (Dijkstra & Loeb 2008; Eide et al. 2018). Thus, to explore polarized Lyα, we adopt a new asymmetric model using the radiative transfer code RT-scat (Chang et al. 2023; Chang & Gronke 2024). We describe the geometry of our wind model and the resulting Lyα spectrum in Section 4.1 and polarization in Section 4.2. We extensively compare observations and models in Section 4.3.

4.1. Lyα spectrum in the wind model

Our new model is composed of a bipolar wind and a point Lyα source. We show a schematic illustration of the model in Fig. 3. The wind model is characterized by four main parameters: the H I column density, N_H I, the expansion velocity, v_exp, the half opening angle of the bipolar wind, θ_o, Wind, and the width of intrinsic Lyα, σ_Src. The bipolar wind outflows expand radially with a constant velocity v_exp, similar to the shell model. The wind’s temperature is fixed at 10⁴ K. The wind’s H I number density is constant, and its inner radius is fixed to 10% of the outer radius. We adopt this geometry as the bipolar wind is intended to represent the outflows outside of the galaxy. In our model, we assumed that there is no H I outside the wind outflows, although some H I from inflowing gas or satellites might be present. Given that different H I properties can have opposite effects on the resulting polarization and that we do not have any evidence to constrain these properties, we excluded the presence of additional H I outside the wind outflows, leaving it for future spatially resolved studies (see the discussion in Section 5.2).

Fig. 3. Schematic illustration of the wind model composed of a central point source (orange) and a bipolar outflow with the radius Ro (red). The central source emits Lyα photons, following a Gaussian profile with a width σSrc. The bipolar outflow is characterized by the H I column density NH I, the opening angle θo, Wind (increasing from the +z-axis, as the blue solid arrow), and the expansion velocity vexp parameters. The inner radius of the outflow Ri is fixed at 0.1 Ro. As the wind model is symmetric about the z-axis, the line of sight angle θLOS is the angle from the +z-axis following the orange arrow, and thus with θLOS = 0° meaning observing in the direction of the outflow, and θLOS = 90° representing the equatorial view.

The central source emits Lyα with a Gaussian profile with a width of σ_Src. This assumption is chosen as it represents the most general case, and is supported by the observational evidence that the galaxy is dust poor (both from the study of optical lines, the blue UV-continuum β slope, and the low E(B − V) reddening in Iani et al. 2021 and from the non detection of dust continuum from ALMA observations in Zanella et al. 2024), that could suppress or modify the Lyα shape (Laursen et al. 2009). Additionally, we consider the angle of the line of sight θ_LOS, the azimuthal angle from the +z-axis, as the bipolar wind is symmetric about the z-axis. In the simulations, we consider 10⁶ photons and extract the escaping Lyα spectrum for various θ_LOS.

In our new model we assumed a simplified geometry as it allows us to focus on the physical processes that originate polarization. This approach is analogous to that commonly employed to analyze the spectra, where the shell model is usually adopted as the standard model and, even if it does not reflect reality and can be affected by many degeneracies (e.g., Gronke et al. 2016; Li & Gronke 2022), it can help us to decrypt the information in the Lyα line. Currently, there is not a similar standard in the joint study of polarization and spectra together, and we decided to adopt the bipolar wind model as it allows us to explore a large variety of scenarios in a well-motivated physical frame.

Fig. 4 shows the observed I(λ) spectrum around the Lyα line, the simulated spectrum from the shell model, and the simulated spectra of the wind model. Iani et al. (2021) estimated the physical properties of the shell model (N_H I ∼ 10²⁰ cm⁻², v_exp ∼ 200 km s⁻¹, and σ_Src ∼ 100 km s⁻¹), that we adopted as the starting point of our wind model. We found that all the simulated spectra of the wind model for different θ_o, Wind values (e.g., we show θ_o, Wind = 30°, 60°, and 75° in the left panel of Fig. 4) do not match the observed spectrum, especially near the red peak of the Lyα. This is because scattered photons escape outside the outflow opening angle, in the equatorial direction, unlike in the shell model.

Fig. 4. Comparisons between observed and simulated Lyα spectra. The simulated spectra of the wind models are normalized by setting the same peak height as the spectrum of the previous fit with the shell model. The black dashed line is the observed Lyα spectrum, which is continuum subtracted. The gray dashed line is the simulated Lyα spectrum of the shell model from Iani et al. (2021). In the left and middle panels line colors represent, for σSrc equal to, respectively, 100 km s−1 and 150 km s−1, θo, Wind = 30° (red), 60° (orange), and 75° (blue), with θLOS fixed to 0°. In the right panel, θo, Wind is fixed to 60°, and line colors represent θLOS = 0° (orange), 40° (blue), and 90° (green).

To address this discrepancy, we assumed a broader intrinsic Lyα (σ_Src = 150 km s⁻¹) to better reproduce the observations. This wider intrinsic spectrum can stem from radiative transfer effects within the inner ISM (e.g., Gronke et al. 2018). We explored the effect of this broadening in Appendix C, finding that the intrinsic Lyα profile can broaden from 50 km s⁻¹ to 150 km s⁻¹ due to radiative transfer effects occurring within the inner ISM before the photons penetrate into the wind. In addition, we showed in Appendix D that a higher random motion of the bipolar wind cannot reproduce the observed spectrum. In conclusion, a broader intrinsic Lyα is required for our modeling with asymmetric geometry.

With this assumption, simulated spectra for θ_o, Wind = 60° and 75° are similar to those of the previous shell modeling and to the observations, as can be seen in the middle panel of Fig. 4.

The right panel of Fig. 4 shows the simulated spectra for three θ_LOS at θ_o, Wind = 60°. The spectra at θ_LOS = 0° and 40° match the observed spectrum well. On the contrary, the spectrum at θ_LOS = 90° does not resemble the observed spectrum and has enhanced flux in the vicinity of the systemic velocity. In the case of a large θ_LOS, as this one, intrinsic photons are observed directly, without scattering, resulting in a central peak and enhanced red wing in the simulated spectrum. In summary, to reproduce the red peak and the suppressed blue peak of the observed Lyα, we need a larger σ_Src = 150 km s⁻¹ and θ_LOS to be smaller than θ_o, Wind. The larger σ_Src allows Lyα photons to be emitted over a wider velocity range, broadening the red peak, while the smaller θ_LOS ensures that most intrinsic photons undergo scattering within the wind cone, suppressing the blue Lyα peak emission.

4.2. Polarization of Lyα in the wind model

Dijkstra & Loeb (2008) showed that different models, that assume different production mechanisms of the Lyα photons and different geometries of the gas surrounding the emitting source, give rise to different polarization levels. From an observational point of view, we detect the signal from an ensemble of photons and not from individual ones, hence the main factor influencing the detected polarization is the geometry of the scattering medium and the presence of a preferential polarization direction, as the signal from highly polarized individual photons may result in an average not polarized signal if their polarization angles are not aligned. Thus, to understand the polarization behavior of scattered Lyα, two fundamental mechanisms to develop polarization are important.

First, the polarization of the integrated Lyα strongly depends on the symmetry of the scattering geometry (Eide et al. 2018). Thus, in Fig. 5, at θ_LOS = 0°, the total degree of polarization P becomes zero due to symmetry, regardless of θ_o, Wind. Similarly, in the right panel, the total P at θ_o, Wind = 90° is zero for any line of sight. The overall P at θ_LOS = 0° and those at θ_o, Wind = 90° are approximately zero. Additionally, as θ_o, Wind increases, the simulated Lyα halo becomes more symmetric, leading to a decrease in P (panels from left to right in Fig. 5).

Fig. 5. Degree of polarization P for models with θLOS from 0° to 90°, with steps of 10°. Given its low significance, we do not display the degree of polarization when the flux in the simulated intensity spectrum is smaller than 5% of the peak flux of the Lyα line. The left, middle, and right panels show the results for θo, Wind = 30°, 60°, and 90°, respectively. The overall degree of polarization increases with increasing θLOS and decreasing θo, Wind.

Second, the degree of polarization P of scattered photons increases when the scattering angle, defined as the angle between the incident and scattered directions, approaches 90° (e.g., Chandrasekhar 1960; Chang et al. 2017; Seon et al. 2022). Scattering at angles close to 0° (forward scattering) or 180° (backward scattering) causes the P of the scattered photon to be identical to that of the incident photon. Based on this mechanism, when θ_LOS is close to 90°, the fraction of photons undergoing perpendicular scattering increases, leading to a higher P. This behavior is evident in the left and middle panels of Fig. 5, where the overall P increases with increasing θ_LOS. As a result, P decreases with increasing θ_o, Wind and decreasing θ_LOS, making the information on geometrical properties of the H I medium imprinted in the polarization of the Lyα line.

Another key factor determining the degree of polarization in symmetric geometries is the H I column density. When Lyα emission is polarized and spatially diffused through scatterings, its polarization strongly depends on N_H I. However, this dependence is non-monotonic. Dijkstra & Loeb (2008) found that the polarization at N_H I = 10²⁰ cm⁻² is weaker than at N_H I = 10¹⁹ cm⁻², due to a higher number of scatterings. In other words, the polarization increases as N_H I decreases.

Conversely, in the lower N_H I regime (< 10¹⁸ cm⁻²), the polarization degree decreases with decreasing N_H I due to the dominance of different types of Lyα scatterings (Seon et al. 2022; Chang et al. 2023). Core (resonance) scatterings, which occur at the line center in the rest frame of the H I atom, generally produce weaker polarization compared to wing (Rayleigh) scatterings that occur far from the line center (see also, Ahn & Lee 2015). Consequently, metal resonance lines exhibit lower polarization compared to Lyα, as discussed further in Section 5.3.

However, as discussed above, symmetric geometries inherently result in unpolarized spectropolarimetric data due to their symmetry, making them insufficient for reproducing the observed polarization properties of Lyα. Thus, in this work we focus on modeling observed data with asymmetric H I geometries.

4.3. Comparing observations and models

We compared our observations of Abell 2895a with the simulated results by using both their total intensity spectra I(λ), and the polarization fraction of the Lyα, P_Lyα. To properly compare the I(λ), we convolved those from simulations with a Gaussian function, to take into account the instrumental spectral resolution of ∼3.6 Å (R ∼ 1500). We adopted the following physical properties of the wind model (also summarized in Table 1): N_H I = 10²⁰ cm⁻², v_exp = 200 km s⁻¹, and σ_Src = 150 km s⁻¹. We considered a large variety of geometries for the modeled biconical outflows, varying the wind opening angles θ_o, Wind from 0° to 90° with steps of 15° and the line of sight angles θ_LOS from 0° to 90° with steps of 10°.

We evaluated the goodness of the matching of their total intensities I by adopting a reduced χ_ν² metric over N = 136 wavelength elements around the Lyα emission line, defined as

$$\begin{array}{c}\hfill {\chi }_{\nu }^2=\frac{1}{N}{\displaystyle \sum _{\lambda =1206\mathrm{\AA }}^{1225\mathrm{\AA }}}{\left[\frac{I(\lambda )-I_{\mathrm{model}}(\lambda )}{{\sigma }_I(\lambda )}\right]}^2,\end{array}$$(3)

where I(λ), σ_I(λ), and I_model(λ) are, respectively, the observed, its uncertainty, and the modeled normalized total intensity spectra. Thus, models with lower χ_ν² values were preferred, as they indicate that the observed and modeled spectra are more in agreement. We remark that, given that the flux measurements are correlated between several pixels, the χ_ν² values must not be interpreted in an absolute way, but rather qualitatively. We decided to adopt this metric because it can offer a straightforward and clear visualization of the goodness of the agreement, and allows one to directly compare different models, and evaluate those ruled out by I(λ).

The observed Lyα shows an asymmetric profile that is redshifted with a relative velocity Δv = 403 ± 4 km s⁻¹ with respect to the systemic redshift (Iani et al. 2021), as it is typical for outflows. In Fig. 6, we compared observed data and simulated results from the wind models for various θ_LOS and θ_o, Wind. As we discussed in Section 4.1, only the models with θ_LOS < θ_o, Wind are able to reproduce these spectral features (see the top panels of Fig. 6). At θ_LOS > θ_o, Wind, the peak of the simulated spectrum is centered on Δv = 0 km s⁻¹, due to directly escaping photons from the central source. However, at θ_o, Wind ≤ 15°, even if θ_LOS < θ_o, Wind, the spectrum does not reproduce the observed red wing since a small θ_o, Wind induces less scattering.

Fig. 6. Top and third rows: Observed (black, with 1σ uncertainties in gray) and modeled (with different colors denoting different line-of-sight angles, θLOS, reported in the legend) normalized total intensity spectra I, assuming a H I column density NH I of 1020 cm−2, an outflow velocity vexp of 200 km s−1, and a Gaussian width of intrinsic LyασSrc of 150 km s−1, in the case of θo, Wind = 0° ,15° ,30° ,45° ,60° ,75° ,90°, from left to right, as indicated in the gray labels. Second and bottom rows: Polarization fraction relative to the model described in the corresponding row above. The black open circles represent the observational 1σ upper limit of PLyα = 4.6%. The uncertainties for the models are smaller than the linewidths.

The cases with θ_o, Wind = 45° and 60° are those better representing the observations, for θ_LOS < θ_o, Wind, and have the lowest χ_ν² values. They present an asymmetric profile and can well reproduce the observed red tail, in particular for θ_LOS ≲ 40°. The case with θ_o, Wind = 30° (75°) has a slightly larger χ_ν² value, because the peak is less (more) redshifted than the observations, but consistent within 1σ. The trend continues to the θ_o, Wind = 90° case, the one resembling the expanding ellipsoid geometry of the gas, that is only marginally consistent with the observed total intensity spectrum. The observed and modeled I(λ) spectra for all the models can be seen in the top panels of Fig. 6, while their χ_ν² values are shown in Fig. 7 in the grayscale, where lighter tones represent better agreement.

Fig. 7. Comparison between observations and wind models at NH I = 1020 cm−2, vexp = 200 km s−1, σSrc = 150 km s−1. The sketches along the axes give a visual representation of the models, with increasing θo, Wind along the y-axis and increasing θLOS along the x-axis, represented with arrows colored according to Fig. 6. The gray scale indicates the agreement between the observed and modeled total intensities I, evaluated through the χν2 value (colorbar on the right), with lighter tones representing better agreement (see the top row of Fig. 6). Red hatched regions rule out models whose degree of polarization PLyα is larger than (and thus not consistent with) the observational 1σ upper limit PLyα (see the bottom row of Fig. 6).

We binned the degree of polarization from the models with the same binning adopted for the observations, and considered the polarization of the bin including the observed Lyα peak, from 1215.5 Å to 1219.6 Å. Given that we only had upper limits on P_Lyα from the observations, we considered as consistent those models whose Lyα polarization fraction, within 1σ, is lower than the observational upper limit. We show the P_Lyα from the models and from the observations in the bottom row of Fig. 6. As discussed in Section 4.2, the low P_Lyα value suggests that the observed system is fairly symmetric, and thus the perfectly symmetric θ_o, Wind = 0° and 90° cases with P_Lyα ≃ 0 are fully consistent with observations. The case with θ_o, Wind = 15° (30°) presents a P_Lyα value that increases with θ_LOS, always consistent (consistent for θ_LOS < 60°) with the observational constraint. In the cases with θ_o, Wind = 45°, 60°, and 75°, the binned P_Lyα values do not increase with increasing θ_LOS, due to the non-vanishing polarization close to the systemic redshift (Δv ≈ 0), which can also be seen in the center panel of Fig. 5 (in particular for the θ_LOS = 50° and 60° curves). It results in larger degrees of polarization that can exceed the observational upper limit for the models with θ_o, Wind ∼ θ_LOS. The models ruled out by the observational constraints (at the 1σ level) on the polarization are in hatched red in Fig. 7. We highlight that the I(λ) and P_Lyα constraints rule out complementary regions in the θ_LOS − θ_o, Wind plane, proving the effectiveness of combining these different tracers to investigate the geometry of the scattering H I gas around star-forming galaxies at high redshift and the mechanism of production of the Lyα photons.

5. Discussion

5.1. Origin mechanisms of the Lyα

As introduced in Section 4.2, there are two main mechanisms that originate the Lyα photons that we observe from distant sources. The first, usually referred to as “in situ”, includes both the recombination (in photo-ionized gas) and the collisional excitation scenarios. Recombination happens when an electron is captured by a proton, resulting in a hydrogen atom in an excited state that can eventually decade to the ground state, emitting Lyα (Haiman & Rees 2001; Cantalupo et al. 2005; Arrigoni Battaia et al. 2019). The probability of emitting Lyα mainly depends on the temperature and density of the medium and, in the usually adopted B-case recombination regime, the probability is as large as approximately 68% at T = 10⁴ K (e.g., Dijkstra 2014), also explaining why the Lyα is the intrinsically stronger emission line in star-forming galaxies. The collisional process instead involves an electron and a hydrogen atom that is left in an excited state after the close encounter. This process converts thermal energy of the electrons into radiation, and is thus also called “Lyα emission by cooling radiation”. The second mechanism is related to the scattering, as the Lyα photons that are created by a central source (e.g., a star-forming galaxy or AGN), can then escape the production site after numerous scatters in the surrounding H I cloud (Hayes et al. 2011; Beck et al. 2016), causing them to propagate far from their production site. Likely these two mechanisms act simultanesouly (e.g., Kim et al. 2020). Polarization can be detected in the photoionization case, if enough neutral hydrogen is present in extended regions, producing a considerable scattering probability. Moreover, after many scatterings, the emergent Lyα line can be depolarized, especially in environments where the medium is isotropic. We can thus conclude that, if the Lyα is polarized, there is photon scattering, while if the Lyα is not (or low) polarized, there is no scattering, the gas has a specific geometry, or there are many scatterings.

The Lyα intensity spectrum of Abell 2895a requires a scattering contribution to explain the asymmetric profile and the redshifted peak of the line with respect to the systemic redshift. On the other hand, the low polarization upper limit, P_Lyα, is also consistent with photoionization only and no scatter, or with scatter in a fairly symmetric medium that, given the lack of a preferential scattering direction, would result in a measured low polarization signal, even if the individual photons might be highly polarized. For example, as described in Section 4.3 for Abell 2895a, centrally emitted photons that experience scattering can give rise to an integrated P_Lyα ∼ 0, consistent mostly with the large θ_o, Wind and low θ_LOS cases, that are the most symmetric ones. These preferred values of θ_o, Wind and θ_LOS are also consistent with the observed displacement of 1.2 ± 0.2 kpc (on the source plane, Iani et al. 2021) between the UV stellar continuum and the peak of the Lyα, if we assume that the Lyα photons are emitted by the stars and that only one side of the wind is illuminated (it can be explained with different scenarios, e.g., an initial transient feature before breaking through the disk, an enhanced dust obscuration of the far wind, or an effect of tidal stripping). Deep spatially resolved polarization observations, which we discuss in the next section, are needed to reconstruct in detail the production site of Lyα photons and the geometry of the scattering medium, allowing us to disentangle between the different possible production mechanisms of the Lyα.

5.2. Spatially resolved observations and future perspectives

Scattering of Lyα photons through a geometrically asymmetric medium gives rise to polarized signal (Angel 1969; Lee & Ahn 1998), detectable also without resolving the system (Eide et al. 2018), as we described in Sections 3 and 4 for Abell 2895a. However, as discussed in Section 5.1, this leaves the question about the origin of the Lyα photons (and thus the region from where they are emitted) and the geometry and orientation of the scattering medium. This degeneracy can be broken by learning the geometry of the galaxy from other independent probes, or by using spatially resolved and deeper spectropolarimetric observations. Spatially resolved observations can put constraints on the polarization over different regions that we can compare with models predicting different degrees of polarization ranging, for example in a biconical wind geometry, from P_Lyα ∼ 0% in the center to P_Lyα up to 80% in the outflows (Eide et al. 2018). Additionally, deeper observations can help putting tighter constrains to the measured polarization levels, and achieving observational S/N high enough to measure also the polarization angle, that is predicted to significantly change between different models, giving hints on, for example, the direction of the outflow in a biconical wind geometry (Eide et al. 2018). As described in Section 3 and in Appendix B, we tentatively measured the polarization fraction by considering only half of the Lyα emission, with the goal of including these information in the comparison with the models. But, due to the low S/N and angular resolution of our observations and the moderate tangential stretch from lensing, we could only measure Lyα polarization upper limits of approximately 7.5% at 1σ.

Currently, Abell 2895a is the only distant star-forming galaxy with Lyα spectropolarimetric observations. The main reasons for the lack of such observations are the typical low luminosity of high-z star-forming galaxies and the fact that the light beam has to be split in two components (its o and e components, with orthogonal polarizations) and rotated, significantly increasing the observational exposure time necessary to reach sufficient S/N values, even with the most advanced spectropolarimeters mounted on cutting-edge telescopes. This issue can be mitigated by observing galaxies at z ∼ 2 − 4 (with redshifted Lyα line included in the instrument’s wavelength coverage) that are strongly lensed (i.e., magnified and distorted) with a large magnification factor (μ ≳ 10), that would sensibly lower the needed exposure time and increase the spatial resolution. The observation of the most magnified and distorted sources lying across the lensing critical lines would allow us to zoom-in and study the spatial variation of the polarization across the Lyα emission on sub-kpc scales and put more stringent limits to the origin of the Lyα emission and the scattering H I geometry. Additionally, extending such observations to a small sample of galaxies will also allow us to put the results obtained for Abell 2895a in a broader context.

5.3. Polarization of other resonant lines

In principle, the study presented in this paper can be executed by employing different resonant line emissions than the Lyα, such as the Mg II, C IV, O VI, N V, Si IV doublets, commonly used in astrophysics (e.g., Prochaska et al. 2011; Hayes et al. 2016; Henry et al. 2018; Berg et al. 2019; Katz et al. 2022; Dutta et al. 2023). The scattering processes of such resonance doublets are very similar to those of the Lyα line, because they present similar atomic properties with one electron in the outer orbit and have the same atomic structure composed of two transitions, called “K” (S​_1/2 − P_3/2) and “H” (S​_1/2 − P_1/2). Due to their resonance nature, metal resonance doublets are also spatially extended and polarized via scattering. The inclusion of such lines would allow us to expand the usable redshift range (atomic line center wavelengths from ∼1032 Å for the O VI to ∼2800 Å for the Mg II) and to explore the multi-phase nature of the CGM, exploring gas temperatures from 10^3 − 4 K (ionization energies of 13.6 eV for the Lyα and 15 eV for the Mg II) to 10^5 − 6 K (138 eV of the O VI). Moreover, the two K and H transitions that are not possible to observationally disentangle for the Lyα line, because of the small energy difference of ∼10⁻⁴ eV ≈1.5 km s⁻¹, are resolved for other lines, such as the Mg II (transition K at 2796.4 Å and H at 2803.5 Å, with a separation of 760 km s⁻¹) and the Si IV (K at 1393.8 Å and H at 1402.8 Å, with the largest separation of 1926 km s⁻¹ between the above-mentioned resonant lines).

Recently, radiative transfer models that include the Mg II doublet have been developed, also including polarization (Seon 2024; Chang & Gronke 2024). In particular, Chang & Gronke (2024) considered three-dimensional shell, sphere, and clumpy sphere geometries, and studied the joint Lyα and Mg II escape using RT-scat. They reveal that, despite being driven by similar atomic processes, the emerging Mg II and Lyα spectra are very different, because of the different Mg II and H I column densities, with N_H I ≫ N_Mg II. Moreover, they confirmed a correlation between the escape of LyC radiation and the Mg II doublet ratio, defined as the ratio of the K and H lines (see also Henry et al. 2018; Chisholm et al. 2020; Izotov et al. 2022; Katz et al. 2022; Xu et al. 2023). They also found that the Mg II degree of polarization decreases with increasing column density N_Mg II, because the multiple scatterings without a preferential direction decrease the resulting polarization. This results in a low (< 5%) degree of polarization in the center, growing up to 10% (25%) with N_Mg II = 10¹⁴ cm⁻² (N_Mg II = 10¹³ cm⁻²). Consequently, they suggested that the Mg II polarization can be used to estimate the doublet ratio of the extended halo. Unfortunately, the Mg II emission is much fainter than the Lyα, making it challenging to resolve and measure its degree of polarization, and requiring extremely large amount of observational time, especially at high redshift. Again, strongly lensed candidates would allow us to obtain spatially extended and spectropolarimetric observations and to pursue this kind of studies.

6. Summary and conclusions

In this paper we put novel constraints on the geometry of the H I region surrounding a clumpy star-forming galaxy at z ≈ 3.4, Abell 2895a (Livermore et al. 2015; Iani et al. 2021; Zanella et al. 2024), strongly lensed into three multiple images by the cluster of galaxies Abell 2895. We made use of new VLT/FORS2 observations taken with the Polarimetric Multi Object Spectroscopy (PMOS) mode to extract the spectra relative to the ordinary (o) and extraordinary (e) beams in four different half-wave plate angles, from an aperture including the M1 and M2 multiple images of Abell 2895a. We combined the different spectra to measure the total intensity 1D spectrum, I(λ). We focused on the Lyα emission line, and we computed the Stokes parameters q and u, used to measure the polarization fraction P_Lyα, corrected both for the positive bias of polarization and for the dilution effect due to the contamination of the angularly close, unpolarized, BCG of the Abell 2895 cluster. We obtained P_Lyα upper limits of 4.6%, 5.8%, and 6.5% at the 1σ, 2σ, and 3σ level.

To interpret the observational constraints, we developed a Lyα radiative transfer model including a bipolar wind geometry (characterized by the half opening angle θ_o, Wind, the radial expansion velocity v_exp, and the H I column density N_H I parameters) and a central source that emits Lyα photons with a Gaussian profile with a width of σ_Src. We assumed N_H I ∼ 10²⁰ cm⁻², v_exp ∼ 200 km s⁻¹, and σ_Src ∼ 150 km s⁻¹ and tried different θ_o, Wind values (from 0° to 90° with steps of 15°) and simulated different line of sight angles (θ_LOS from 0° to 90° with steps of 10°).

In order to reproduce the spectral profile and shift with respect to the systemic velocity of the Lyα line, we needed Lyα photons that are created in the inner regions and that undergo several scattering events before escaping, and θ_LOS < θ_o, Wind and θ_o, Wind > 15°. Models with θ_o, Wind of 45°, 60°, or 75° with θ_LOS of 30° −60° can well reproduce the Lyα spectral profile, but the combination of θ_o, Wind and θ_LOS resulting in an asymmetric geometry and the scattering due to the larger θ_o, Wind, make the predicted polarization much larger than the low observed P_Lyα, that is consistent with no scattering or with scattering in a fairly symmetric medium, and thus are excluded.

Summarizing, the models that satisfy both the Lyα spectral profile and polarization requirements are those with θ_o, Wind ∼ 30° for θ_LOS ≤ 20°, θ_o, Wind ∼ 45° for θ_LOS ≤ 20°, θ_o, Wind ∼ 60° for θ_LOS ≤ 20°, θ_o, Wind ∼ 75° for θ_LOS ≤ 40°, and θ_o, Wind ∼ 90° for any θ_LOS, where θ_LOS = 0° means observing in the direction of the outflow. We showed that polarization constraints are complementary to those usually achieved by using only the total Lyα intensity and line profile, demonstrating the effectiveness of this novel technique. These results will pave the way to future spatially resolved spectropolarimetric observations, needed to discriminate between the different production mechanisms of the Lyα and the geometry of the scattering medium.

Acknowledgments

We are grateful to the referee, Matthew Hayes, for the insightful comments and constructive suggestions, which have significantly enhanced the quality of this paper. This research is based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 108.2260, It is also based on observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESA) and the Canadian Astronomy Data Centre (CADC/NRC/CSA). The research activities described in this paper have been co-funded by the European Union – NextGeneration EU within PRIN 2022 project n.20229YBSAN – Globular clusters in cosmological simulations and in lensed fields: from their birth to the present epoch. The authors thank Sabine Moehler for the helpful discussions regarding the FORS2 PMOS data reduction. A.B. and A.Z. acknowledge support from the INAF minigrant 1.05.23.04.01 “Clumps at cosmological distance: revealing their formation, nature, and evolution”. M.G. thanks the Max Planck Society for support through the Max Planck Research Group. E.I. acknowledges funding from the Netherlands Research School for Astronomy (NOVA). We acknowledge the use of the numpy (Harris et al. 2020), matplotlib (Hunter 2007), astropy (Astropy Collaboration 2022), and pandas (The pandas development team 2023) packages.

References

Appendix A: Lyα polarization measured in the three-bin case

We show here the results obtained by dividing the Lyα line over three bins, one including the blue tail (1215.5-1216.3 Å), a central one including the peak (1216.3-1217.9 Å) and one including the red tail (1217.9-1219.6 Å). The binning is completed by sampling the blue and red continua with two bins: one narrower and closer to the line and one broader, from 1210.0 Å to 1215.5 Å, and from 1037 Å to 1210 Å for the blue, respectively, and from 1219.6 Å to 1224.0 Å, and from 1224 Å to 1260 Å for the blue, respectively. The bins are shown in the left panels of Fig. A.1. For the blue tail, central, and red tail bins, we measure 1σ upper limits on P_Lyα of 6.6%, 4.2%, and 13.1%, respectively. Due to the large dilution factor, the polarization degree is barely (not) constrained for the narrow (broad) continua bins.

Fig. A.1. Left panels: Total intensity spectra (green solid line) and 1σ uncertainties (light green shaded region) in the spectral region around the Lyα line, similar to those presented in Fig. 2, but in the three-bin case. The top (bottom) panel shows the normalized spectrum before (after) the subtraction of the BCG contribution. The black filled squares, with 1σ uncertainties, represent the binned data. Right panels: Polarization (P) measurements obtained before (top) and after (bottom) the dilution correction described in Eq. 1. The black open circles represent the 1σ upper limits obtained by applying the correction of Simmons & Stewart (1985). In the bottom right panel, only two datapoints are visible, as the others are outside the plotted range (upper limits for P from 13 to 95%).

Appendix B: Tentative spatially resolved Lyα polarization measurement

We show here the results obtained by dividing the spatial extraction aperture over two different regions, in order to verify whether the low polarization fraction measured might be the result of the cancellation of opposite contributions from different sides. Unfortunately, the S/N is not sufficient to perform such test on the individual M1 and M2 multiple images but, thanks to their mirrored nature, we designed two apertures to include the same spatial half of each Lyα emission. In particular, the first aperture has a width of ≈4″ and goes from the M1 Lyα peak to the M2 Lyα peak.

The second aperture has a total width of 6″, divided into two apertures of 3″, from each Lyα peak toward the edge of the slit. Considering the case with a single bin for the Lyα line, we obtain upper limits of 5.1%, 7.1%, and 8.3% at the 1σ, 2σ, and 3σ level for the first aperture, while 7.5%, 10.2%, and 11.9%, respectively, for the second aperture. The results are shown in Fig. B.1 and B.2. The results show that low polarization levels are measured in different spatial regions, and thus that the golbal low polarization is not the result of cancellation of the Stokes vectors from two different sides. Unfortunately, due to the low S/N and spatial resolution, it is not possible to include the spatial information in our analysis and in a more detailed comparison with the model.

Fig. B.1. Left panels: Total intensity spectra (green solid line) and 1σ uncertainties (light green shaded region) in the spectral region around the Lyα line, similar to those presented in Fig. 2, but extracted from a region designed to include a spatial half of the Lyα blob. The aperture, represented in orange in the inset, has a width of ≈4″ and goes from the M1 Lyα peak to the M2 Lyα peak. The top (bottom) panel shows the normalized spectrum before (after) the subtraction of the BCG contribution. The black filled squares, with 1σ uncertainties, represent the binned data. Right panels: Polarization (P) measurements obtained before (top) and after (bottom) the dilution correction described in Eq. 1. The black open circles represent the 1σ upper limits obtained by applying the correction of Simmons & Stewart (1985).

Fig. B.2. As in Fig. B.1, but with the spectra extracted from the other spatial half of the Lyα emission. In particular, the extraction aperture, represented in orange in the inset, has a total width of 6″, divided into two apertures of 3″, from each Lyα peak toward the edge of the slit.

Appendix C: Lyα Line Broadening

In Section 4.1, our modeling with the wind geometry required the intrinsic Lyα emission with σ_Src = 150 km s⁻¹ to reproduce the observed spectrum. However, this is three times broader than the observed width of ∼50 km s⁻¹ of the Hβ line (Iani et al. 2021). This is puzzling as the intrinsic Lyα width should be similar to the width of the observed Balmer lines. Similar discrepancies have been noted previously in the literature (Yang et al. 2016; Orlitová et al. 2018). In this section, we aim to understand the broad intrinsic Lyα through the effect of Lyα radiative transfer in a small-scale H I gas in ISM or star-forming regions since its resonance nature induces the broadening of the escaping spectrum and can potentially explain the differences between the line widths (Neufeld 1990; Gronke et al. 2018).

To test this broadening effect, we assumed a spherical geometry composed of a central Lyα point source surrounded by a static spherical H I halo. The line width of the central Lyα emission was fixed at 50 km s⁻¹, corresponding to the Hβ width. Two types of H I gas distributions were considered: a smooth medium with uniform H I density and a clumpy medium composed of small H I clumps in an empty inter-clump medium, where the clump size is 100 times smaller than the halo size. The smooth medium is characterized by the H I column density, N_H I, and the random motion, σ_ran. The clumpy medium is described by the total H I column density, N_H I,total, the clump random motion, σ_cl, and the covering factor, f_c, which represents the mean number of clumps along the line of sight from the central source (see, e.g., Hansen & Oh 2006). The total column density N_H I,total is f_cN_H I,cl, where N_H I,cl is the H I column density of a single clump. Details of this geometric setup are described in Chang & Gronke (2024). Fig. C.1 illustrates simulated spectra (dotted lines) compared to a Gaussian profile with σ_Src = 150 km s⁻¹ (black solid lines).

Fig. C.1. Simulated spectra in the spherical geometry considering smooth (top) and clumpy (bottom) media. The black solid lines represent Gaussian functions with a width of σ = 150 km s−1, corresponding to the intrinsic Lyα profile estimated using the wind model in Section 4.1.

For the smooth medium, we fixed σ_ran at 50 km s⁻¹, corresponding to the intrinsic Lyα width. However, in the top panel of Fig. C.1, the simulated spectrum exhibits a double-peak profile that does not match the Gaussian shape. This outcome arises because optically thick H I medium at the Lyα line center generally produces double-peaked profiles (Neufeld 1990). To resolve this, we decreased N_H I for optically thin gas at the line center and increase σ_ran to enhance spectral broadening. When σ_ran reaches 150 km s⁻¹, the simulated spectrum closely matches the Gaussian profile.

In the clumpy medium, we considered small values of f_c since low covering factors enable single-peak profiles via surface scattering (Neufeld 1991; see also Hansen & Oh 2006; Chang et al. 2023) when f_c is below the critical covering factor (Gronke et al. 2016, 2017). In the bottom panel of Fig. C.1, the simulated spectra successfully reproduce the Gaussian profile with σ_Src = 150 km s⁻¹.

In summary, the intrinsic Lyα profile can broaden easily from 50 km s⁻¹ to 150 km s⁻¹ due to radiative transfer effects occurring within the inner ISM before the photons penetrate into the wind. We can simulate this effect within the context of a clumpy medium in which a relatively small number of clumps (f_c ≲ f_c, crit) can lead to a Gaussian line broadening, or in a smooth medium using a large turbulent motion σ_ran and a low N_H I. These findings underscore the distinct mechanisms by which Lyα line formation is influenced by the physical properties of small-scale cold gas.

Appendix D: Impact of turbulent motion within a bipolar wind on the formation of Lyα

Galactic winds, as most astrophysical media, are highly turbulent which is imprinted, for instance, in the large random velocities seen in both simulations and observations orthogonal to the bulk flow direction (e.g., Schneider et al. 2020; Veilleux et al. 2020). In this section, we want to investigate the effect of this turbulent component on the emergent Lyα spectrum.

Fig. D.1 shows Lyα spectra for three different random motions of the bipolar wind: σ_ran = 12.8, 100, and 200 km s⁻¹.² The spectrum at σ_ran = 12.8 km s⁻¹ corresponds to the model at θ_o, Wind = 60° in the left panel of Fig. 4. As σ_ran increases, the spectrum becomes more redshifted and deviates from the observed spectrum. This behavior arises because higher random motion induces a more significant peak shift in the Lyα line.

Fig. D.1. Simulated Lyα spectra for various random motions of the H I bipolar wind (σran = 12.8, 100, and 200 km s−1) at θo, Wind = 60° and θLOS = 0°. The intrinsic line width σSrc is fixed at 100 km s−1, based on the shell model in Iani et al. (2021). The black and gray dashed lines represent the observed spectrum and the shell model spectrum, respectively. The orange line corresponds to the spectrum at θo, Wind = 60° from Fig. 4. The navy and violet lines are spectra for σran = 100 and 200 km s−1, respectively. Higher random motion induces redshifting and deviation from the observed spectrum. The simulated spectra are normalized by setting the same peak height as the spectrum of the previous fit with the shell model.

Consequently, accounting for increased random motion within this set of other parameters, the bipolar wind model fails to reproduce the observed spectrum. However, naturally other degeneracies exist which are, however, beyond the scope of this study to explore. Interestingly, the shell-model fit of Iani et al. (2021) produces an effective temperature of log(T_eff/K) = 5.3 ± 0.2, that is, assuming a temperature of T = 10⁴ K corresponding to a random motion of ∼50 km s⁻¹.

All Tables

All Figures

Fig. 1. Left: HST F606W image of the inner part of A2895, where the three multiple images of Abell 2895a (M1, M2, and M3) appear. It represents the rest-frame UV at the redshift of Abell 2895a, z ∼ 3.4. In orange, we show the 2, 3, and 5σ contours of the Lyα emission, detected with MUSE. The green box represents the FORS2 1.4″ × 22″ adopted slit, that includes M1 and M2, and one image of another source at z ∼ 3.7 (Iani et al. 2023), whose Lyα 2, 3, and 5σ contours are shown in red. Thanks to the PMOS mode, the signal included in the slit is split in an ordinary (o) and extraordinary (e) ray, with orthogonal polarizations, shown with different tones of green. Right: Stacked 2D spectrum with φ = 0.0° obtained after calibration, cosmic ray rejection, and sky subtraction. In particular, the sky is subtracted only in the relevant slits which contain Abell 2895a, and appear darker. Each slit has a spatial (vertical) width of 22″, and the sky level is estimated and then subtracted from the closest ordinary and extraordinary slits (see the e and o labels on the right), respectively, as shown by the green arrows on the right. In the slits including the target, we detect the extended Lyα emission associated with M1 and M2, at approximately 5340 Å, the continuum from the nearby BCG, in the bottom part, and the Lyα emission of the source at z ∼ 3.7 (red squares, close to the edge of the slitlets). We show the 2D spectrum taken with φ = 0.0° on top, and smaller cutouts around the Lyα line for φ = 22.5°, 45.0° and 67.5° on the bottom.

In the text

Fig. 2. Left panels: Total intensity spectra (green solid line) and 1σ uncertainties (light green shaded region) in the spectral region around the Lyα line. The top (bottom) panel shows the normalized spectrum before (after) the subtraction of the BCG contribution. The black filled squares, with 1σ uncertainties, represent the binned data. Right panels: Polarization (P) measurements obtained before (top) and after (bottom) the dilution correction described in Eq. (1). The black open circles represent the 1σ upper limits obtained by applying the correction of Simmons & Stewart (1985). In the bottom right panel, only one datapoint is visible, as the others are outside the plotted range (upper limits for P ∼ 80%−95%).

In the text

Fig. 3. Schematic illustration of the wind model composed of a central point source (orange) and a bipolar outflow with the radius Ro (red). The central source emits Lyα photons, following a Gaussian profile with a width σSrc. The bipolar outflow is characterized by the H I column density NH I, the opening angle θo, Wind (increasing from the +z-axis, as the blue solid arrow), and the expansion velocity vexp parameters. The inner radius of the outflow Ri is fixed at 0.1 Ro. As the wind model is symmetric about the z-axis, the line of sight angle θLOS is the angle from the +z-axis following the orange arrow, and thus with θLOS = 0° meaning observing in the direction of the outflow, and θLOS = 90° representing the equatorial view.

In the text

Fig. 4. Comparisons between observed and simulated Lyα spectra. The simulated spectra of the wind models are normalized by setting the same peak height as the spectrum of the previous fit with the shell model. The black dashed line is the observed Lyα spectrum, which is continuum subtracted. The gray dashed line is the simulated Lyα spectrum of the shell model from Iani et al. (2021). In the left and middle panels line colors represent, for σSrc equal to, respectively, 100 km s−1 and 150 km s−1, θo, Wind = 30° (red), 60° (orange), and 75° (blue), with θLOS fixed to 0°. In the right panel, θo, Wind is fixed to 60°, and line colors represent θLOS = 0° (orange), 40° (blue), and 90° (green).

In the text

Fig. 5. Degree of polarization P for models with θLOS from 0° to 90°, with steps of 10°. Given its low significance, we do not display the degree of polarization when the flux in the simulated intensity spectrum is smaller than 5% of the peak flux of the Lyα line. The left, middle, and right panels show the results for θo, Wind = 30°, 60°, and 90°, respectively. The overall degree of polarization increases with increasing θLOS and decreasing θo, Wind.

In the text

Fig. 6. Top and third rows: Observed (black, with 1σ uncertainties in gray) and modeled (with different colors denoting different line-of-sight angles, θLOS, reported in the legend) normalized total intensity spectra I, assuming a H I column density NH I of 1020 cm−2, an outflow velocity vexp of 200 km s−1, and a Gaussian width of intrinsic LyασSrc of 150 km s−1, in the case of θo, Wind = 0° ,15° ,30° ,45° ,60° ,75° ,90°, from left to right, as indicated in the gray labels. Second and bottom rows: Polarization fraction relative to the model described in the corresponding row above. The black open circles represent the observational 1σ upper limit of PLyα = 4.6%. The uncertainties for the models are smaller than the linewidths.

In the text

Fig. 7. Comparison between observations and wind models at NH I = 1020 cm−2, vexp = 200 km s−1, σSrc = 150 km s−1. The sketches along the axes give a visual representation of the models, with increasing θo, Wind along the y-axis and increasing θLOS along the x-axis, represented with arrows colored according to Fig. 6. The gray scale indicates the agreement between the observed and modeled total intensities I, evaluated through the χν2 value (colorbar on the right), with lighter tones representing better agreement (see the top row of Fig. 6). Red hatched regions rule out models whose degree of polarization PLyα is larger than (and thus not consistent with) the observational 1σ upper limit PLyα (see the bottom row of Fig. 6).

In the text

Fig. A.1. Left panels: Total intensity spectra (green solid line) and 1σ uncertainties (light green shaded region) in the spectral region around the Lyα line, similar to those presented in Fig. 2, but in the three-bin case. The top (bottom) panel shows the normalized spectrum before (after) the subtraction of the BCG contribution. The black filled squares, with 1σ uncertainties, represent the binned data. Right panels: Polarization (P) measurements obtained before (top) and after (bottom) the dilution correction described in Eq. 1. The black open circles represent the 1σ upper limits obtained by applying the correction of Simmons & Stewart (1985). In the bottom right panel, only two datapoints are visible, as the others are outside the plotted range (upper limits for P from 13 to 95%).

In the text

Fig. B.1. Left panels: Total intensity spectra (green solid line) and 1σ uncertainties (light green shaded region) in the spectral region around the Lyα line, similar to those presented in Fig. 2, but extracted from a region designed to include a spatial half of the Lyα blob. The aperture, represented in orange in the inset, has a width of ≈4″ and goes from the M1 Lyα peak to the M2 Lyα peak. The top (bottom) panel shows the normalized spectrum before (after) the subtraction of the BCG contribution. The black filled squares, with 1σ uncertainties, represent the binned data. Right panels: Polarization (P) measurements obtained before (top) and after (bottom) the dilution correction described in Eq. 1. The black open circles represent the 1σ upper limits obtained by applying the correction of Simmons & Stewart (1985).

In the text

	Fig. B.2. As in Fig. B.1, but with the spectra extracted from the other spatial half of the Lyα emission. In particular, the extraction aperture, represented in orange in the inset, has a total width of 6″, divided into two apertures of 3″, from each Lyα peak toward the edge of the slit.
In the text

	Fig. C.1. Simulated spectra in the spherical geometry considering smooth (top) and clumpy (bottom) media. The black solid lines represent Gaussian functions with a width of σ = 150 km s−1, corresponding to the intrinsic Lyα profile estimated using the wind model in Section 4.1.
In the text

Fig. D.1. Simulated Lyα spectra for various random motions of the H I bipolar wind (σran = 12.8, 100, and 200 km s−1) at θo, Wind = 60° and θLOS = 0°. The intrinsic line width σSrc is fixed at 100 km s−1, based on the shell model in Iani et al. (2021). The black and gray dashed lines represent the observed spectrum and the shell model spectrum, respectively. The orange line corresponds to the spectrum at θo, Wind = 60° from Fig. 4. The navy and violet lines are spectra for σran = 100 and 200 km s−1, respectively. Higher random motion induces redshifting and deviation from the observed spectrum. The simulated spectra are normalized by setting the same peak height as the spectrum of the previous fit with the shell model.

In the text