# Observations of a Magellanic Corona

This work uses archival HST/COS and FUSE spectra to show evidence for the Magellanic Corona. Here we describe our data reduction, Voigt profile fitting and ionization modelling methods. Throughout this work, all reported values and uncertainties are medians and 68% confidence intervals, unless otherwise specified.

### Impact parameters and projection effects

Unlike CGM studies of extragalactic systems, this work focuses on the CGM surrounding the LMC at a distance of only D = 50 kpc (ref. 31). This proximity means that background quasars at large angular separations from the LMC, θ, correspond to relatively small physical separations. The impact parameter, ρLMC, is found using ρLMC = Dsin(θ). At θ > 45°, this assumption no longer results in realistic impact parameter estimates and we would require a true 3D model of the location of gas absorbers to calculate a physical separation between gas absorbers and the LMC. Furthermore, for sightlines at large θ, it is harder to kinematically distinguish absorption lines from the LMC and the Milky Way. A larger-scale understanding of the Magellanic Corona and multiphase CGM would only be possible with more reliance on models and simulations to identify the 3D locations of gas absorbers. To keep this work more focused on the observationally derived results, our analysis is thus strictly limited to sightlines within 45° of the LMC, corresponding to 35 kpc in impact parameter.

### HST/COS observations

We design our sample to consist of HST/COS far UV observations of background quasars using both the G130M and G160M gratings, covering the wavelength ranges of about 1,150–1,450 Å and about 1,405–1,775 Å, respectively. The combination of these gratings enable us to examine the following absorption lines: O I λ 1302, N I λλ 1199, 1200, 1200.7, C II λ 1334, Al II λ 1670, Si II λλ 1260, 1193, 1190, 1526, 1304, Si II λλ 1250, 1253, 1259, Fe II λλ 1608, 1144, Si III λ 1206, C IV λλ 1548, 1550 and Si IV λλ 1393, 1402. The COS spectra are processed following previously developed custom reduction and wavelength-calibration methods10,32 based on the raw products from the calcos33 data-reduction pipeline. To remove geocoronal airglow contamination in O I λ 1302 and Si II λ 1304, we use a second calcos reduction of the data, using only observations taken during orbital night-time.

The COS far UV observations have a native pixel size of 2.5 km s−1 and a spectral resolution (full width at half maximum) of about 20 km s−1 and about 15 km s−1 for G130M and G160M spectra, respectively. We bin all spectra such that the resulting spectra are Nyquist sampled with two pixels per resolution element.

### FUSE observations

For 15 sightlines in our sample, archival FUSE spectra are also available and analysed to search for O VI λ 1031, 1037 absorption. However, only six sightlines had high enough S/N ratio to make a measurement. These wavelengths fall on the FUSE LIF1 channel with a spectral resolution of about 20 km s−1 and native pixel size of 2 km s−1, which we bin to Nyquist sample with two pixels per resolution element. These FUSE data are reduced and aligned following the customized methods similar to those used for the HST/COS spectra34,35. The O VI λ 1031 may have contamination from molecular H2 absorption at λ 1032.356, which corresponds to roughly 130 km s−1 in the O VI frame. However, the expected contribution from this contamination is very small because of the high Galactic latitude of the sightlines and is, in most cases, negligible.

### Absorption-line measurements

We use the open-source Python software, VoigtFit36, to perform Voigt profile fitting of the absorption in several ions observed with HST/COS with the G130M and G160M gratings. This process uses a least-squares optimizer37 with recent atomic data38,39,40 and convolves the Voigt profile with an approximate instrumental profile of a Gaussian with full width at half maximum corresponding to the observed grating resolution. Although this Gaussian approximation of the instrumental profile is not an exact representation, it has been shown to have a nearly negligible effect on fit results for weak high-velocity components21. For all ion fits, we normalize the spectra using a third-order polynomial fit to the continuum surrounding absorption lines of interest. Regions of the absorption spectra that are contaminated by high-redshift absorption components are then flagged to avoid fitting.

We first fit the absorption in all low and intermediate ions (O I, N I, C II, C II*, Si II, Si III, Al II, Fe II) simultaneously, allowing component-line centres to be tied across ions when they show general agreement. The C II* line always contaminates the measurement of absorption in C II at +250 km s−1. When there are blended C II components at this velocity, we fix the C II* column density to a constant value of 1013.8 cm−2, based on average measurements from previous work41, but—in these cases—the measured C II columns near +250 km s−1 are not used in our analysis. The Si III λ 1206 transition is frequently saturated, requiring the linewidths to be tied to match the fit Si II linewidths. A minimum allowed linewidth of 9 km s−1 is applied on the basis of the instrumental resolution and maximum linewidths are only added as a constraint for highly blended components if they are needed to converge to a best fit.

C IV and Si IV are then fit simultaneously following the same procedure, but independent of the low-ion results to avoid biasing our analysis, because the high-ion component structure may be different. O VI absorption from FUSE is also fit independently when data are available and a reasonable continuum can be determined. If the component structure of the low and high ions match, they are flagged after the fitting process so that their column densities, linewidths and line centres can be compared in the subsequent steps. Furthermore, we calculate upper limits of any transitions in which absorption is not seen on the basis of the S/N ratio of the observed spectra42,43. Last, fit components attributed to the Milky Way or known intermediate-velocity or high-velocity clouds are flagged to avoid contaminating our analysis. We note that some contamination from fixed pattern noise persists in our reduced spectra, which may affect our measured column densities and is not accounted for in our estimated errors.

In total, across 28 sightlines, we initially identify 112 unique velocity components that may be attributed to the Magellanic system. We then impose a velocity threshold and only consider absorbers at vLSR > 150 km s−1 to avoid contamination from absorbers associated with the Milky Way44. The velocity threshold of 150 km s−1 was determined using a combination of the observed component velocities and simulations of the Magellanic system7; it represents the value that best separates the Galactic and Magellanic components and is consistent with previous kinematic studies of Magellanic absorption12,14. Furthermore, this velocity threshold is supported by dynamical arguments: given the LMC mass, the virial theorem predicts that Magellanic gas has a velocity dispersion of 50 km s−1 centred on the LMC velocity of 280 km s−1, implying that 95% of Magellanic gas should be within 180 km s−1 and 380 km s−1. As a result, our final sample has 52 unique Magellanic velocity components that are further analysed on the basis of their kinematics and photoionization modelling. The Voigt profile model parameters for these 52 C IV and Si IV components are given in Extended Data Table 1 and the ten unique Magellanic O VI absorption components are shown in Extended Data Table 2. Extended Data Fig. 1 shows our measured C IV λ 1548 and O VI λ 1031 absorption-line spectra for our sample. Extended Data Fig. 2a shows the total measured HST/COS column densities in several low and high ions from the Magellanic absorbers at vLSR > 150 km s−1 as a function of the LMC impact parameter. All low ions show a declining radial profile, similar to the relation shown in the high ions (Fig. 2).

A comparison of our observed radial profile with that seen in the COS-Dwarfs survey23 and M31 (ref. 24) is shown in Extended Data Fig. 2b. We normalize impact parameter measurements across these surveys based on the radius enclosing a mean overdensity of 200 times the critical density, R200, which is often used as a measure of the virial radius in CGM studies. In the radial region of overlap between these surveys and our work, the declining profile of the LMC is more concentrated, with a possibly truncated profile. Because the LMC halo is already within the virial radius of the Milky Way, it is expected to be tidally truncated, hence such a truncated profile is expected. However, the uncertainties in estimates of R200 are estimated to be 50% in the COS-Dwarfs and M31 surveys, with the LMC value we use at R200 = 115 ± 15 kpc.

Our spectra can be accessed publicly on the Barbara A. Mikulski Archive for Space Telescopes (MAST). A full table of our fit parameters, along with summary plots of our best fits, can be accessed on GitHub at https://github.com/Deech08/HST_MagellanicCorona.

### Ionization models

We use 1D Cloudy45 radiative transfer models to simulate the physical conditions of the absorbing gas. Our Cloudy models require four key inputs to run: (1) an external radiation field, (2) the observed column density measurements, (3) a specified stopping condition to reach for convergence and (4) a gas-phase metallicity. All models assume a plane-parallel geometry and constant gas density.

Incident radiation fields in Cloudy require a shape and intensity. We adopt the Milky Way escaping radiation field model to set the shape of the radiation field, assuming that the radiation fields from the LMC and the SMC have the same spectral shape10,46,47. The intensity of the radiation field towards each sightline is set by a hydrogen-ionizing photon flux ΦH determined from published ionization models, which includes contributions from the LMC, the SMC and the Milky Way20. We reconstruct this model in 3D space to interpolate an initial value for ΦH for any specified location. In our model, we allow ΦH to be a free parameter, because a precise distance to the absorbing material is not known. We also include a constant contribution from an extragalactic UV background48 and cosmic ray background49.

We use Cloudy’s built-in ‘optimize’ command to vary our free parameters and find optimal parameters to explain our observed column densities and upper limits45,50. The optimize models use up to three possible free parameters: (1) the hydrogen-ionizing photon flux, ΦH, described above, (2) the total hydrogen number density, nH, which is the sum of the ionic, atomic and molecular hydrogen densities of the plasma that is to be modelled, and (3) the neutral hydrogen column density (NH I) stopping condition. For sightlines with an H I or O I detection, the observed H I or O I column-density measurement serves as the stopping condition and the model only uses the first two free parameters (ΦH and nH). For sightlines without an H I or O I, we use all three free parameters (ΦH, nH and NH I). Once Cloudy’s optimize method has found a possible solution of parameters, we run one final Cloudy model at the specified optimal parameters to produce predictions of ion column densities and gas temperatures, including predictions for high-ion (Si IV, C IV, O VI) column densities. To ensure Cloudy does not settle at local minima in the optimization process, we use a broad range of initial densities from log10(nH/cm−3) = −3 to 1 and ionizing fluxes ΦH in a range of 3 dex around the model prediction at D = 50 kpc, but still find a resulting narrow range of total hydrogen densities (nH), ionized gas temperatures (Te), neutral atomic hydrogen columns (NH I) and ionized-to-neutral atomic hydrogen ratios N(H II)/N(H I) across all sightlines and velocity components. Furthermore, we have also run a coarse grid at a larger range of free parameters to help confirm that our solutions are indeed optimal and not local minima.

Although interstellar medium gas-phase metallicities have been measured in the LMC, SMC, Magellanic Bridge51 and Magellanic Stream52,53,54,55, the metallicity of the Magellanic CGM is highly uncertain. To estimate the gas metallicity, we use a sightline in our sample towards HE 0226–4110 that overlaps with recently published analysis of FUSE spectra to measure neutral hydrogen column densities56. Two absorption components towards this sightline may belong to the Magellanic Corona at vLSR = +174 km s−1 and +202 km s−1, providing a measured neutral hydrogen column density to set as a stopping condition in Cloudy. Unfortunately, there is no detected O I absorption in either the COS or the FUSE data, so a metallicity is calculated using a Cloudy optimize model (described above), allowing the total hydrogen density, hydrogen-ionizing photon flux and metallicity to vary. The Cloudy models are optimized on the basis of the measured COS column densities across all available metal ions and any upper limits when absorption is not detected. The results for these two components are log10(ΦH/photons s−1) = 5.06, log10(nH/cm−3) = −1.58 and [Z/H] = −0.72 for the vLSR = +174 km s−1 component and log10(ΦH/photons s−1) = 4.95, log10(nH/cm−3) = −1.91 and [Z/H] = −0.62 for the vLSR = +202 km s−1 component. On the basis of these results, we adopt the average [Z/H] = −0.67 as the gas-phase metallicity for photoionized gas. For hotter gas in interfaces and the corona, we assume a gas-phase metallicity of [Z/H] = −1, because we expect this more primordial gas to be at lower metallicity.

Our optimal set of Cloudy models provides predictions for the expected column densities of the high ions Si IV, C IV and O VI for a single-phase photoionized gas. However, the observed high-ion columns are much greater (by orders of magnitude) than the photoionization predictions. Across all sightlines and absorption components that may be associated with the Magellanic system, we find that 72% of Si IV and 84% of C IV absorption components are less than 10% photoionized. We use this 10% (1-dex) threshold to define our sample of Magellanic absorbers that are not photoionized (see shaded components in Fig. 2). These C IV and Si IV absorbers probably arise in interfaces in the range T = 104.3–4.9 K.

The observed triply ionized Magellanic absorption is well described using either equilibrium or time-dependent non-equilibrium collisional ionization models16. In both cases, we can infer an electron temperature based on the ratio of C IV and Si IV column densities, because the close similarity of the C IV and Si IV line profiles indicates that the two ions are co-spatial. The modelled relation of this column-density ratio with temperature for the equilibrium model and for isobaric and isochoric time-dependent models is shown in Extended Data Fig. 3 for a range of metallicities. The inferred temperature is then used to determine a C IV ionization fraction, from which the total ionized hydrogen (H II) column density can be calculated, resulting in the measurements shown in Fig. 3b. In total, the temperature distributions of the photoionized and collisionally ionized gas are shown in Extended Data Fig. 4b. In the sightlines in which we have measured O VI absorption, we find that the O VI absorbing gas requires a higher temperature than the C IV and Si IV absorbing gas, indicating that the O VI arises in a separate, hotter phase. Whereas at high metallicity lower-temperature solutions for our observed column-density ratios are possible, this is not the case at the lower metallicities (below 0.1 solar) expected for Magellanic coronal gas.

We also consider more recent collisional ionization models that include photoionization from an extragalactic background57. However, these models do not include the non-isotropic radiation fields necessary for modelling clouds near the Milky Way and the LMC, and only offer approximate predictions using a general background radiation field. Instead, we only consider the two cases of entirely photoionized or entirely collisionally ionized in this work, but note that a full picture will require considering collisional ionization and photoionization from the Milky Way and the Magellanic Clouds together.

### Statistical significance of results

Here we describe the statistical tests we used to support our claims of significance. Throughout this work, we adopt a significance threshold P-value of P = 0.05.

#### Velocity structure

In our Voigt profile fitting process, individual components are initially paired across low ions and high ions based on their approximate centroid velocities. This pairing process is inherently biased, as it assumes that components across ions are physically tied and results in the lowest possible differences in velocity centroids for our analysis. However, for the low and high ions, the velocity structure was qualitatively well matched to one another, with absorption components at similar velocities for both cases. This correspondence is less clear for the O VI absorption line centroids, so matching O VI components in the same manner is much more uncertain. Combined with the relatively low S/N ratio (≈10) and moderate velocity resolution (20 km s−1) of our spectra, we are unable to fully resolve all absorption components. We therefore find that comparisons of the kinematic properties of low and high ions are generally inconclusive. However, the kinematics are still consistent with our primary conclusion that C IV and Si IV arise in the interfaces between cool clouds and a Magellanic Corona, because in an interface model the velocity structure of the low ions and the high ions should be linked. When considering O VI, we calculate the velocity offset from the closest absorption component in other ions (Si III or C IV) and find that the widths of the velocity-offset distributions have standard deviations of $${\sigma }_{{\rm{O}}{\rm{VI}}-{\rm{Si}}{\rm{III}}}=2{2}_{-4}^{+7}\,{\rm{km}}\,{{\rm{s}}}^{-1}$$ and $${\sigma }_{{\rm{O}}{\rm{VI}}-{\rm{C}}{\rm{IV}}}=2{2}_{-6}^{+10}\,{\rm{km}}\,{{\rm{s}}}^{-1}$$, respectively. This is $${7}_{-8}^{+12}\,{\rm{km}}\,{{\rm{s}}}^{-1}$$ greater in width of the distribution of velocity differences between the low ions and C IV matched in the same manner, supporting the result that O VI exists in a different phase.

#### Linewidths

We show the paired (matched on the basis of their velocities during the Voigt profile fitting process) differences of component linewidths in Extended Data Fig. 5a. Differences in paired linewidths do not show statistical significance. However, when considering our populations of linewidth measurements, we do find a statistically significant difference between the linewidth distributions of singly ionized C and Si in comparison with triply ionized C and Si (see their distributions in Extended Data Fig. 4c,d). The Anderson–Darling statistical test of the null hypothesis that the singly and triply ionized linewidths are drawn from the same underlying population can generally be rejected at the P-value threshold of 0.05 for both C and Si. We perform the test on 10, 000 bootstrap samples to account for measurement errors of linewidths. The C IV and C II linewidths return a P-value (with 68% confidence intervals) of $${P}_{{\rm{C}}}=0.00{8}_{-0.007}^{+0.08}$$, with 78% of bootstrap samples below our P-value threshold of 0.05. Similarly, the Si IV and Si II linewidths return P-values of $${P}_{{\rm{Si}}}=0.00{1}_{-0.0}^{+0.01}$$, with 93% of bootstrap samples below our significance threshold.

We test the statistical significance of the anti-correlation between the C IV and Si IV with LMC impact parameter using Kendall’s τ rank correlation coefficient with censoring, which provides a robust measure of the monotonic relationship between two variables58,59. The Magellanic Corona shows a distribution of coefficients that are negative for both C IV and Si IV, with mean values of τ = −0.4 ± 0.1 and τ = −0.3 ± 0.1, respectively, as shown in Extended Data Fig. 6. The P-values for C IV allow the null hypothesis of no correlation to be rejected at the 0.05 level for 97% of bootstrap samples, whereas the P-values for Si IV can only be rejected for 73% when considering all of our data. When only considering the absorbers at ρLMC > 7 kpc, the significance of the Si IV anti-correlation becomes stronger, with P < 0.05 for 89% of 10,000 bootstrap samples and a mean value of τ = −0.4 ± 0.1, but the change for C IV is negligible. The best-fit lines for the anti-correlation are found using a Markov chain Monte Carlo analysis with censoring to account for upper limits and measurement errors60. For the O VI measurements, Kendall’s τ rank correlation coefficient is less reliable, as we only have six data points, and is not conclusive.

### Magellanic Corona versus tidally stripped stream with interfaces

Previous simulations have been able to explain much of the ionized gas associated with the Magellanic Stream by tidal stripping, without the presence of a corona61. If this were the case, and the Stream were the dominant source of ionized gas, we would expect to see a stronger correlation of C IV column density as a function of distance from the Magellanic Stream (absolute Magellanic Stream latitude) than as a function of the LMC impact parameter. We use the partial Spearman rank-order correlation test to assess the strength of the correlation between our measured ion column densities and either the LMC impact parameter or the absolute Magellanic Stream latitude, while removing the effects of the other. We note that, for this test, we are only considering the collisionally ionized C IV and Si IV columns, but considering all the observed columns for low ions. The correlation coefficients and P-values of the test with a null hypothesis of no correlation are given in Extended Data Table 3. For most ions, the correlation is substantially stronger with the LMC impact parameter, after removing the effects of the absolute Magellanic Stream latitude. However, the partial correlation test for Fe II is inconclusive and the test for O I suggests a stronger correlation with the absolute Magellanic Stream latitude. These tests are consistent with a Magellanic Corona and CGM origin to the gas absorbers we have measured, with the exception of O I, which may be more biased towards tracing cooler, tidally stripped gas in the Magellanic Stream.

In Extended Data Fig. 7, we show our measurements of collisionally ionized C IV columns on a map of the Magellanic system in Magellanic coordinates, alongside measurements of all C IV absorption from a previous survey of the Magellanic Stream10. When considering all C IV, the surface density profile is much more extended along the direction of the Magellanic Stream, but with our adopted velocity threshold and removal of photoionized gas, the radial profile centred on the LMC is apparent, especially when considering sightlines that overlap on our sample and the previous sample.

In this previous work, much of the observed C IV absorption was interpreted to arise from interfaces around the tidally stripped, cooler gas from the LMC with a hot, approximately 106 K Milky Way corona. The basic premise of this conclusion is still valid in our sample, but the strong radial profile centred on the LMC suggests that the hotter gas interacting to form the interfaces should also be centred on the LMC, not the Milky Way. Therefore, a Magellanic Corona at roughly 105.5 K can explain our observed radial profile and the observed C IV absorption.

### Mass estimates

Our estimates of the mass for each phase of the Magellanic CGM are derived from the relation between the ionized hydrogen column density and the LMC impact parameter. For each phase (approximately 104 K, approximately 104.9 K and approximately 105.5 K), a best-fit linear regression model is fit to the ionized hydrogen column as a function of ρLMC. Then the ionized hydrogen mass in each phase is calculated using

$${M}_{{\rm{H}}{\rm{II}}}={\int }_{0\,{\rm{kpc}}}^{35\,{\rm{kpc}}}\,{N}_{{\rm{H}}{\rm{II}}}\left({\rho }_{{\rm{LMC}}}\right){m}_{{\rm{p}}}\,2\pi {\rho }_{{\rm{LMC}}}\,{f}_{{\rm{cov}}}\,{\rm{d}}{\rho }_{{\rm{LMC}}}$$

(1)

in which mp is the proton mass and fcov is the covering fraction.

For the approximately 104 K gas, the ionized hydrogen column density in each direction is derived directly from the Cloudy models, with a covering fraction fcov = 0.82, as low ions are detected at Magellanic velocities in 23/28 directions in our sample. However, we note that the covering fraction of low ions tends to decrease as a function of the LMC impact parameter, but use a constant covering fraction as an approximation.

For the 104.9 K gas, the total ionized hydrogen column density in each sightline is derived on the basis of the C IV column density and best-fit temperature from the collisional ionization models16 using

$${N}_{{\rm{H}}{\rm{II}}}=\frac{{N}_{{\rm{C}}{\rm{IV}}}}{{f}_{{\rm{C}}{\rm{IV}}}\left[{\rm{Z/H}}\right]},$$

(2)

in which fC IV ≡ C3+/C is the fraction of triply ionized carbon at the best-fit temperature and [Z/H] = 0.21[Z/H] is the metallicity. For C IV, the covering fraction is set to fcov = 0.78 for ρLMC < 30 kpc and fcov = 0.3 for ρLMC ≥ 30 kpc based on the observed detection rate of C IV absorption in our sample. The relation between the derived column densities and the LMC impact parameter allow for our mass calculations for the approximately 104 K and approximately 104.9 K gas to converge, changing by at most 0.1 dex if instead we integrate out to 500 kpc.

For the approximately 105.5 K gas, the mass is again found on the basis of the O VI absorption columns in the collisional models, using

$${N}_{{\rm{H}}{\rm{II}}}=\frac{{N}_{{\rm{O}}{\rm{VI}}}}{{f}_{{\rm{O}}{\rm{VI}}}\left[{\rm{Z/H}}\right]},$$

(3)

and using the same covering fraction correction used for the approximately 104.9 K gas. Here we use the maximal fO VI value for each of the collisional models, which peak near 105.5 K at fO VI ≈ 0.2. The best-fit line for this phase does not converge, so integrating the radial profile depends highly on the radial range considered. Instead, we only integrate between the bounds of our observations (6.7 kpc < ρLMC < 32.5 kpc) and present an approximate corona mass for this region only.