© The Authors 2023

1 Introduction

The Tarantula Nebula (or 30 Doradus) in the Large Magellanic Cloud (LMC) is a vast, intrinsically bright star-forming region (Kennicutt 1984; Doran et al. 2013). With a large population of young, massive stars (M ≳ 8 M_⊙) at a metallicity of about half Solar (see Mokiem et al. 2007, for an overview), this region is reminiscent of giant starbursts observed in distant galaxies (Cardamone et al. 2009; Crowther et al. 2017). The proximity of 30 Doradus and its relatively unobscured nature allow us to resolve the stellar population of this starburst-like environment. This provides a unique opportunity to calibrate integrated quantities required for studies of starburst galaxies in which the stellar populations are unresolved (Thornley et al. 1999; Brandner 2002; Garcia et al. 2021).

At the heart of 30 Doradus lies the cluster R136, which hosts a rich population of (very) massive stars. In the centre of this young (1–2 Myr), massive (M ~ 2–10 × 10⁴ M_⊙; Hunter et al. 1995; Andersen et al. 2009) cluster are the most massive stars observed to date (de Koter et al. 1997; Crowther et al. 2010; Bestenlehner et al. 2020; Brands et al. 2022; Kalari et al. 2022). With its many (very) massive stars, R136 is responsible for 30–50% of the overall ionising luminosity and wind power of 30 Doradus (Crowther 2019; see also Doran et al. 2013), and not surprisingly, the cluster plays an important role in stellar feedback mechanisms throughout this region (e.g. Pellegrini et al. 2011; Lee et al. 2019; Cheng et al. 2021).

In the current work, we focus on the extinction properties of sightlines towards R136. For the study of individual massive stars, knowledge of extinction is required in order to recover the intrinsic luminosity, which is an important factor in the determination of stellar masses; of particular interest are the masses of the three very massive WNh stars in the core of R136. Moreover, knowledge of the extinction provides insight into the dust properties of this starburst-like environment. This, in turn, might help us to understand the distant starburst galaxies, where stellar populations cannot be resolved and a proper characterisation of the dust properties is key for recovering the intrinsic ultraviolet (UV) brightness and identifying star formation (e.g. Salim & Narayanan 2020, and references therein). To properly account for the dust within (unresolved) starburst galaxies, usually an attenuation law is used rather than extinction law. Attenuation laws take into account not only the removal of light from the line of sight (extinction), but also the scattering of light into the line of sight. A widely used starburst attenuation law is that of Calzetti et al. (1994, 2000). For an in depth overview of dust attenuation laws in galaxies, we refer the reader to Salim & Narayanan (2020).

We investigate the dust properties by analysing the wavelength dependence of the extinction, often referred to as the ‘extinction law’. Extinction laws contain valuable information about the interstellar dust, as grain populations of different sizes are thought to account for the extinction at different wavelengths (Draine 2003). At long wavelengths, up to the optical range, extinction increases almost linearly and is associated with silicate and carbonaceous grains with sizes >250 Å. At shorter wavelengths, extinction due to these types of grains flattens, and the continuing increase of extinction in the UV and the peak around 2175 Å are linked to smaller silicates (having sizes <250 Å) and polycyclic aromatic hydrocarbon molecules (PAHs; e.g. Weingartner & Draine 2001; Xiang et al. 2017). The relative fraction of each particle population determines the shape of the extinction curve. While its global shape is rather similar for all sightlines, indicating that the grains are rather uniformly mixed throughout the interstellar medium (ISM), more subtle differences between the curves exist. This shows that environmental factors can influence the grain populations. In this context, processes related to star formation are of particular importance. Examples of such processes are the erosion of molecular clouds under the influence of stellar winds and UV radiation, the shattering of grains by supernova explosions, and the coagulation of dust during cloud collapse (e.g. Galliano et al. 2018).

A parameter often used to characterise the shape of an extinction curve is the total-to-selective extinction, defined as R_V ≡ A_V/E(B – V), with E(B – V) = A_B – A_V, and A_V and A_B being the extinction in the V and B photometric bands, respectively. The sightlines towards R136, and more generally 30 Doradus, display high values of R_V, with average values ranging from R_V = 4.0 to R_V = 4.5 (Doran et al. 2013; Bestenlehner et al. 2014, 2020; Maíz-Apellániz et al. 2014; De Marchi & Panagia 2014; De Marchi et al. 2016; Crowther et al. 2016), similar to sightlines towards star-forming regions in the Milky Way. Within 30 Doradus, large differences between sightlines exist. For example, Maíz-Apellániz et al. (2014) find values in the range R_V = 3.1–6.7. Such variations in R_V within a relatively small physical region are also observed for H II regions in the Galaxy (Maíz-Apellániz & Barbá 2018).

R_V is usually interpreted as a characterisation of the size distribution of dust grains, where sightlines with high values of R_V are associated with a higher fraction of large dust particles, and a lower fraction of small dust particles. For the Milky Way, the dependence of UV extinction on R_V has been investigated intensively, consistently revealing the same relation (e.g. Cardelli et al. 1989; Fitzpatrick 1999; Fitzpatrick et al. 2019). While this might suggest that this behaviour is universal, we stress that this relation is empirical: a priori we would not expect a fixed relation between R_V, which is defined in a narrow optical region, and the extinction in the UV, which is attributed dominantly to different grain populations from those responsible for the extinction in the optical (e.g. Weingartner & Draine 2001; Xiang et al. 2017).

Analyses of extragalactic sightlines indeed show that the Galactic relation between Rv and UV extinction is not universal. For example, Howarth (1983) and Pei (1992) derive average extinction curves for sight lines in the LMC corresponding to R_V ≈ 3.1, the typical Galactic value, while they predict stronger extinction in the UV compared to Galactic laws of the same R_V. Also Gordon et al. (2003), who study two different sets of LMC sightlines, find curves that differ from the Galactic average. De Marchi & Panagia (2019), who investigate three sight lines in 30 Doradus, find that the excess of large grains in this region (associated with the high R_V values) does not seem to come at the expense of small grains. In other words, the high R_V values found for R136 and 30 Doradus are not associated with UV extinction as weak as for Galactic sightlines with similarly high R_V. De Marchi & Panagia (2019) argue that supernova explosions have affected the dust population in the region.

In this paper, we investigate the dust properties in 30 Doradus by studying 50 sightlines in the core of R136. While extinction laws in the near-infrared (NIR) and optical wavelength range of 30 Doradus and R136 have already been evaluated (De Marchi & Panagia 2014; Maíz-Apellániz et al. 2014), a detailed study of the UV extinction is lacking to date. In the present paper, we measure extinction properties towards and around R136, that is, we construct a UV, optical, and NIR extinction law across the region. This yields a spatial map of the extinction in and around the R136 cluster, an (average) extinction law, and insight into the dust properties. We use these results to obtain improved mass estimates of the WNh stars in the core of R136, the most massive stars known.

The remainder of this paper is structured as follows. In Sect. 2, we describe our sample and data, and in Sect. 3 we provide details of the methods used. In Sect. 4, we present our main findings, which we discuss in a broader context in Sect. 5. Finally, we summarise the main conclusions of our study in Sect. 6.

2 Sample and data

Our sample of the R136 core coincides with the sample of Brands et al. (2022), with the omission of six sources. We omit R136a6 because it comprises two sources that cannot be resolved in the UV spectroscopy, and R136a8, H49, H65, H129, and H162 because HST/WFC3 photometry of De Marchi et al. (2011, see below) is lacking for these sources. For all other stars, 50 in total, we compile a spectral energy distribution (SED) that spans from the UV to the NIR. In the UV, we use flux-calibrated spectroscopy, whereas for the optical and NIR we use photometry. An overview of the used observations can be found in Table 1; we describe the data in more detail below.

We use flux-calibrated HST/STIS spectroscopy presented in Crowther et al. (2016) for the far-UV part of the SED (1150–1710 Å). These observations comprise 17 long-slit (52′′ × 0.2′′) contiguous pointings with grating G140L. In each slit, multiple sources are present (see Crowther et al. 2016, and Brands et al. 2022, their Figs. 1). The spectra are extracted using MULTISPEC, a package tailored to extracting spectra from crowded regions (Maiz-Apellaniz 2005). In the extraction process, the position of each source with respect to the slit centre is taken into account for accurate flux calibration. As a consequence of pointing inaccuracies, uncertainties in the calibration remain; we estimate these to be on the order of 10% (see below), but in the most extreme cases they could be as large as a factor two. Small uncertainties in the flux calibration do not pose a problem for our analysis, as long as they are not systematic. To check our flux calibration for systematic uncertainties, we evaluate the integrated flux of the cluster core. To this end, we sum the flux-calibrated spectra of all stars in the cluster core. We then compare this integrated flux to the large aperture (2′ × 2′) HST/GHRS spectrum of the R136 core (Heap et al. 1992), which covers roughly the same area on the sky and for which the flux calibration is assumed to be reliable. Except for wavelengths <1200 Å, which are not included in our SED analysis, we find good agreement between the two. For all wavelengths considered in the fitting, the integrated fluxes match within 25%, and for wavelengths >1350 Å the match is even within 10% (Fig. 1). On average, the sum of the STIS fluxes is 5% lower than that of the GHRS spectrum. A modest difference is not surprising, as there are small discrepancies between the sky coverage and source extraction of the STIS and GHRS observations. All considered, we regard the calibration of the UV fluxes of Crowther et al. (2016) as reliable, and adopt the fluxes at face value.

For the optical part of the SED, we use HST/WFC3 photometry of De Marchi et al. (2011). We have magnitudes in the filters F336W, F438W, F555W, and F814W, which are roughly equivalent to the Johnson U, B, V, and I filters, respectively. For eight stars, the F814W magnitude is missing, and for two stars the F336W or F555W magnitude is missing (see Table 1). We also use the photometry of De Marchi et al. (2011) to estimate the extinction for 1657 stars in the outskirts of R136 (see Sect. 4.2).

The NIR photometry, complete for all stars in the sample, was taken in three bands (H, J, K_S) using VLT/SPHERE (Khorrami et al. 2017, 2021). We note that J and K_S magnitudes are presented in Khorrami et al. (2017); and H and K_S magnitudes in Khorrami et al. (2021). We adopt the H and K_S magnitudes of the most recent analysis, and the J magnitude of the first paper. In order to resolve the individual sources in the crowded R136 core, adaptive optics was employed, making absolute flux calibration challenging. Khorrami et al. (2021) double check the calibration of their H and K_S magnitudes by comparing to the catalogue of Campbell et al. (2010) and report no systematic difference. Bestenlehner et al. (2020) carry out a similar check with the Khorrami et al. (2017) catalogue, but only for the K_S band, as a dataset for cross-checking the J band does not exist. Nonetheless, we adopt the J magnitude calibration and include these magnitudes in our analysis. This is because none of the checks on the absolute flux calibrations on the H and K_S bands give reason for concern, and the J band calibration was carried out by Khorrami et al. (2017) in the same manner as that of the K_S band. We note that, in practice, the J magnitudes have a negligible effect on the outcome of our analysis and removing them from the fitting would not affect our conclusions; this is because their uncertainties are large compared to those of the other bands, which minimises their weight in the fitting process.

Fig. 1 Integrated UV flux of the core of R136. The upper panel shows the GHRS spectrum of the inner 2′ × 2′ (Heap et al. 1992; light blue crosses), as well as the summed STIS spectra used in this work (dark blue circles), with both covering roughly the same inner region of the cluster. The bottom panel shows the ratio of the former two spectra (yellow dots) as well as the central wavelength of the synthetic UV bands used in this work (red triangles, see Sect. 3.1). The dashed lines correspond to ratios of 0.75 and 1.25.

3 Methods

In order to assess the extinction towards the sources in the core of R136, we employ the ‘extinction without standards’ technique (Whiteoak 1966; Fitzpatrick & Massa 2005). This technique requires an observed SED, a model of the intrinsic SED, and an extinction law of which the shape can be parameterised. In addition, the process requires a fitting algorithm to find the optimal extinction curve parameters given the intrinsic and observed SED. We discuss each of these aspects below.

For 1657 stars in the outskirts of R136 (all sources from the De Marchi et al. (2011) catalogue with V < 19), we assess the extinction in a different way, namely by extrapolating the extinction properties we measure for the R136 core stars. To this end, we use the empirical relation between the colour V – I and A_V that we find for the core stars. We describe this process in Sect. 4.2.

3.1 Observed SEDs

We compile observed SEDs spanning from the UV to NIR (see also Sect. 2). The optical and NIR observed magnitudes are converted to fluxes using filter transmission curves taken from the SVO Filter Profile Service¹. For the UV, we bin the observed flux points by taking the flux average over nine different wavelength regions as specified in Table A.1. Taking such a flux average is equivalent to defining a filter with a transmission of 1.0 between the two wavelengths, and 0.0 outside that region. With the 4 optical and 3 NIR photometry points for the R136 core stars, the observed SEDs consist of 16 flux points in total.

The photometric points from the UV to the NIR that we include in our fitting are spread over the SED in such a way that each point represents an approximately equal amount of energy radiated by the star. In other words, if we take an SED and integrate over the wavelength range spanned between the central wavelength of one filter and its neighbour, this amount of energy (per cm² per sec) is roughly equal for each filter. This ensures that the different parts of the SED are as equally represented as possible in the fitting process. We note that the results are robust to changes in the number of UV points that we include in our analysis; a change in the number of points has a negligible effect on the outcome of our fits and would not affect our conclusions. An example of one of our observed SEDs is shown in Fig. 2.

Fig. 2 Example of an observed SED (source: R136a3) covering the UV (dark blue circles), optical (light blue squares), and NIR (turquoise triangles). The UV flux points are derived from flux-calibrated STIS spectroscopy (grating: G140L), and the optical and NIR fluxes come from broadband photometry; the (equivalent) name of each band is indicated.

3.2 Intrinsic SEDs

We compute models of intrinsic SEDs with the model atmosphere code FASTWIND (Santolaya-Rey et al. 1997; Puls et al. 2005; Rivero González et al. 2012; Carneiro et al. 2016; Sundqvist & Puls 2018). This code, which is tailored to hot stars with winds, solves radiative transfer subject to the solution of the non-local thermal equilibrium number-density rate equations and takes into account the effects of line blocking and line blanketing. Only a subset of the spectral lines is explicitly computed and the output SED therefore does not resolve individual spectral lines. For our analysis, this is not directly an issue, as we are not fitting individual lines. Still, the absence of lines in the model SED can result in small differences in flux of the integrated bands. However, the effect of this on the outcome of our analysis is only minor. We mimic the magnitude of the effect by representing the UV absorption-line forest by an overall decrease in the UV model flux of 10%. In this case, the best-fitting UV extinction is only slightly lower (Δc₂ = −0.05) compared to the case where the continuum is unmodified. This change is within the typical range for uncertainties of individual c₂ values. Moreover, the best-fit values for A₅₄₉₅, R₅₄₉₅, and luminosity remain unaffected. We conclude that the effect of using SEDs without explicit spectral lines is only minor.

To further ensure that our models are valid as intrinsic SEDs, we compare FASTWIND models with low mass-loss rates (10⁻¹⁰ M_⊙ yr⁻¹) to hydrostatic models from the TLUSTY grid of Lanz & Hubeny (2003), which have been used for other extinction studies (e.g. MA14). First, we assess how the FASTWIND models behave in the U – B versus B – V plane to make sure that the behaviour around the Balmer jump is correct, as the latter is only modelled in a crude way in the FASTWIND SEDs.

We find that the qualitative behaviour of the two sets of models is similar, and that the absolute differences between the FASTWIND and TLUSTY colours are never more than 2%. Second, we compare the behaviour in the J – H versus H – K plane of the TLUSTY and FASTWIND SEDs. Again, we find that the qualitative behaviour of the two sets of models is similar, and furthermore, that the differences are never larger than 1%. Overall, we conclude that the FASTWIND models are suitable to be used as intrinsic SEDs for this study.

For the stellar parameters, we adopt the values of the optical and UV analysis of Brands et al. (2022) for stars in the core of R136. We note that for computing these models, we need to adopt a stellar radius R_*, while this is one of the parameters that we constrain in this study (see Sect. 3.4). However, for small changes in R_*, the structure of the atmosphere will not change significantly and this means that we can change R_* of the model simply by multiplying the output SED with a factor that represents an increase or decrease in the stellar surface area. Using this approach, we need only one model SED per star, and are still able to fit R_*.

3.3 Extinction law

When considering the SED from the NIR to the UV, none of the existing extinction laws seem suitable for sightlines towards 30 Doradus. The average LMC curves of Howarth (1983) and Pei (1992) correspond to an optical slope of R_V ≈ 3.1, whereas higher values in the range of R_V = 4.0–4.5 are found for 30 Doradus (for references, see Sect. 1). Also, the average LMC curve of Gordon et al. (2003), who find R_V = 3.41, does not match the value of 30 Doradus, nor does their curve towards LMC2 (more commonly referred to as ‘LMC SGS 2’), lying southeast of 30 Doradus. For LMC SGS 2, these latter authors find a rather low value of R_V = 2.76. Misselt et al. (1999) also find low values of R_V, both for stars in LMC SGS 2 (R_V ≤ 3.31) and for stars elsewhere in the LMC (R_V ≤ 2.61)². The law of Maíz-Apellániz et al. (2014) is tailored to 30 Doradus, but is only valid in the optical and NIR.

Also, the Galactic laws, many of which allow for a varying R_V, are not suitable for 30 Doradus when the analysis is extended to the UV. We demonstrate this in Fig. 3, by showing the fits of two stars, with and without UV, while adopting the Galactic law of Fitzpatrick (1999, hereafter F99), which has a typical Galactic R_V-dependent UV extinction. The two stars shown are R136a3 (WNh) and HSH95-55 (or H55, O2 V((f*))z), representative of components within R136. When fitting only the NIR and optical, we reproduce the high values of R_V obtained by previous studies. For example, for R136a3, we find R_V = 4.49 ± 0.57, and for H55, we obtain R_V = 4.37 ± 0.78. Indeed, the lower panels of the figure clearly show that a high R_V better reproduces the observed optical photometry. However, the high values of R_V very poorly reproduce the UV part of the SED. To fit the full SED, that is, including the UV, lower values are needed: R_V = 3.10 ± 0.11 for R136a3 and R_V = 2.94 ± 0.11 for H55. The behaviour shown in Fig. 3 for R136a3 and H55 is observed for all stars in the R136 core, and is also seen when adopting a different Galactic law, such as that of Cardelli et al. (1989) or that of Fitzpatrick et al. (2019).

As no existing extinction law meets our requirements, we need to either adjust an existing extinction law, or derive our own. The data we have available are insufficient to derive a completely new law: in the NIR and optical, we have only a few data points, and in the near-UV (between the U band and the red end of the STIS/G140M grating at 1710 Å), we lack data altogether. Given the complex shape of extinction curves, the number of free parameters we would need to fit would exceed the number of data points we have in these wavelength regions. We therefore decided to modify an existing law.

The parameterisation of our law is similar to that of F99. These latter authors provide an R_V -dependent Galactic law consisting of a cubic spline going through anchor points at fixed wavelengths for the optical and NIR, and a UV part (λ ≲ 3000 Å) that is described by the simple parameterisation of Fitzpatrick & Massa (1990). The strength of extinction at each optical and NIR spline point is tailored to Galactic sightlines and is R_V -dependent. For this work, we do not adopt the optical and NIR spline point values of F99, but rather use the shape of the extinction law of MA14. MA14 provide a family of extinction laws that depend on R₅₄₉₅ ≡ A₅₄₉₅/(A₄₄₀₅ – A₅₄₉₅), the monochromatic equivalent of R_V. We therefore use monochromatic quantities throughout the paper, with the exception of Sect. 5.3, where we compare with other studies that adopted broadband quantities. We note that since the extinction towards R136 is moderate, and the SEDs of the stars we study are fairly similar, the differences between the broadband quantities (A_V and R_V) and their monochromatic equivalents (A₅₄₉₅ and R₅₄₉₅) are small. We refer the reader to Maíz-Apellániz 2013, for a discussion on monochromatic versus broadband quantities.

MA14 use a combination of the seventh-order polynomials of Cardelli et al. (1989) and correction factors on those polynomials to express their law. We convert this functional form to the format of F99, that is, we use a linear relation for expressing the value of each spline point³, so that the dependence on R₅₄₉₅ of each spline point is explicitly provided by a formula for the sake of clarity. Furthermore, we omit the spline points of MA14 for x ≥ 2.7 µm⁻¹ (λ ≥ 3703 Å); for these wavelengths we adopt the UV prescription of Fitzpatrick & Massa (1990, see below). We add an extra spline point at x = 3.0 µm⁻¹ (λ = 3304 Å), and the R₅₄₉₅-dependent extinction at this point we take from the law of MA14. The latter spline point is crucial in order to fit R₅₄₉₅, reflecting the slope of the extinction curve in the optical and NIR, and the slope of the UV extinction as parameterised by Fitzpatrick & Massa (1990), as truly independent parameters. All adopted spline points are listed in Table 2. Our parameteri-sation matches the law of MA14 exactly at the spline points, and within ≤0.1% for all other optical and NIR wavelengths.

Before we continue to describe the UV part of our law, we note that the slope of the NIR part of the extinction law of MA14 is possibly not ideal for extinction towards R136. The NIR part of the MA14 law can be approximated by a power law of the form A_λ ∝ λ^αNIR with α_NIR = 1.61; this value is lower than values obtained by recent studies of Galactic extinction, which find α_NIR = 2.1–2.4 (Stead & Hoare 2009; Nogueras-Lara et al. 2019; Maíz-Apellániz et al. 2020). Upon inspecting the residuals of our fits as a function of A_V, we see clear patterns for all NIR bands, indicating that the law is not completely adequate. However, we adopted it as it is because we have only a few data points in the NIR, and it is beyond the scope of this paper to improve the NIR extinction curve.

For the UV part of the curve, we use the parameterisation of Fitzpatrick & Massa (1990). Expressed in terms of the quantity k(λ – V) ≡ (A_λ – A_V)/(A_b – A_V) ≡ E(λ - V)/E(B – V), the UV part of the law of F99 has the following form: (1)

where x ≡ 1/λ µm⁻¹, a₁ = 0.5392, and a₂ = 0.05644⁴. The parameters c₁ and c₂ relate to a linear background (see below); c₃ indicates the strength of the 2175 Å feature (see below), and c₄ and c₅ indicate the strength and start of the far-UV curvature, respectively. Fitzpatrick & Massa (1990) and F99 do not treat the start of the UV curvature as a free parameter, but adopt a fixed value of c₅ = 5.9; we do the same. We leave c₄ as a free parameter, so that we can constrain the UV curvature in the extinction curve of sightlines towards R136. The 2175 Å feature is described by the Drude profile: (2)

with x₀ being the position and γ the width of the bump. We do not fit the data of the bump as we lack observations here, and instead adopt the values that Gordon et al. (2003) find for LMC SGS 2 (near 30 Doradus), that is, c₃ = 1.463, γ = 0.945, and x₀ = 4.558.

Lastly, and key for our study, we discuss the parameters c₁ and c₂. These parameters do not have a fixed value in the average curve of F99, but instead are expressed in terms of R_V and each other. For Galactic sightlines, F99 derive: (3)

It is the dependence on R_V that links the shape of the extinction in the optical and the NIR to that in the UV. However, while the relation in Eq. (3) is typical for Galactic sightlines, it does not apply to our case (see Fig. 3). Therefore, in the present study, we leave c₂ as a free parameter in the fitting process. In other words, we remove the Galactic dependence of the UV extinction on R_V from our law: the strength of extinction in the UV (c₂) and the slope of extinction in the optical (R₅₄₉₅) will be fitted independently. While we fit the slope of the UV extinction, the equation for c₁ we leave unchanged. This is because c₂ and c₁ are (to some extent) degenerate, and we do not have enough data to break this degeneracy. Experiments where we leave c1 free in the fitting process show that, indeed, it is not possible to disentangle these two parameters with our data.

Briefly, we use the optical and NIR law of Maíz-Apellániz et al. (2014), which we transform to the functional form of Fitzpatrick (1999). For the UV, we adopt the parameterisation of F99 and Fitzpatrick & Massa (1990). The 2175 Å feature is described as in Gordon et al. (2003), and the parameters R₅₄₉₅ (slope in the optical; monochromatic equivalent of R_V), c₂ (strength of UV extinction), and c4 (far-UV curvature) are free parameters.

Fig. 3 Example fits of SEDs of the stars R136a3 (left; WNh) and H55 (right; O2 V((f*))z), adopting the standard Galactic law of F99. The figure shows that this law is not suitable for our analysis: the Galactic dependence of UV extinction on Rv does not hold for sightlines towards R136. In other words, the reddened flux cannot be modelled accurately over the full wavelength range: an RV based on the optical and NIR photometry gives a total mismatch with the observed UV fluxes, whereas a much lower value – that can reproduce the UV fluxes – gives a very poor fit to the observed slope in the optical. The top panels display the flux from UV to NIR; the bottom panels zoom in on the optical regime for which UBVI photometry is available. Circles and squares indicate observed fluxes in the UV, and optical or NIR, respectively. The best-fitting reddened model including the UV is shown in orange; the best fit excluding the UV constraints in red. For this example, we used the law of F99, but the same behaviour is seen with other RV -dependent Galactic laws.

3.4 Fitting

For each individual star, we optimise the free parameters of the adopted extinction law (Sect. 3.3) in order to find the best match between the observed SED and the reddened model SED. We fit five parameters: R₅₄₉₅, affecting the shape of the optical and NIR part of the curve; c₂ and c₄, relating to the UV part of the curve; the extinction in the V band A_V; and the stellar angular radius θ_R= R_*/d, where R_* is the physical radius of the star, and d is the distance. We adopt a fixed LMC distance of d = 49.59 kpc for all stars (Pietrzyński et al. 2019) and therefore in practice the parameter that we vary is R_*, or, equivalently, the stellar luminosity L_*. We assume values for T_eff as described in Sect. 3.2. We stress that changes in the parameters A_V and R_* have different effects and are not degenerate, as changing A_V affects the extinction at each wavelength differently⁵, while changing R_* changes the flux of the intrinsic SED by an equal factor for all wavelengths.

In order to find the best-fitting values of the free parameters A_V, R_*, and those describing the shape of the curve (R₅₄₉₅, c₂, and c₄), we use the minimize function of the Python package lmfit⁶ (Newville et al. 2014). We adopt the least-squares Levenberg-Marquardt method for the minimisation, where we minimise the χ² value: (4)

where N is the number of data points i of the SED that is considered in the fit, ℱ_mod,i the reddened flux of the model, ℱ_obs,i the observed flux of the reddened star, and 𝒪_i is the observational uncertainty on each flux point. For the optical and NIR, we obtain observational uncertainties from the literature; for the UV, we use the standard deviation of the fluxes in each synthetic band as an uncertainty (Sect. 3.1). The reduced-χ², the χ² value per degree of freedom, is defined as , with n_dof = N – n_free the degrees of freedom, where n_free is the number of free parameters.

In order to obtain ℱ_mod,i, we first obtain the distance-corrected reddened model flux, f_λ, by applying the extinction law under consideration to the intrinsic model spectrum: (5)

with F_λ the intrinsic (unreddened), distance-corrected model flux, and A_λ the extinction in magnitudes as a function of wavelength, dictated by the adopted law and A_V. As discussed above, θ_R scales with stellar radius (luminosity) and is a free parameter, as is A_V. After obtaining the reddened model SED, we get the fluxes ℱ_mod,i, by using the transmission curves from filters used for the construction of our observed SEDs.

4 Results

4.1 An R₅₄₉₅-dependent average extinction law for R136

We fit the SEDs of the stars in the core of R136 and for each star we obtain best-fit values and uncertainties for A₅₄₉₅, R₅₄₉₅, c₂, c₄, and luminosity. We find a large range of extinction values throughout the cluster, with A₅₄₉₅ ranging from A₅₄₉₅ = 0.76 ± 0.06 (H108) to A₅₄₉₅ = 2.63 ± 0.05 (H36), and R₅₄₉₅ ranging from R₅₄₉₅ = 1.91 ± 0.26 (H108) to R₅₄₉₅ = 5.79 ± 0.31 (H36), and cluster averages of A₅₄₉₅ = 1.71 ± 0.41 and of R_V = 4.38 ± 0.87. The best-fit parameters of individual stars are listed in Table B.1. Two example fits are shown in Fig. 4.

The of our best fits range from (for H90) to (for H120), with an average of and a median of . These values are relatively high, suggesting that the adopted extinction law does not describe the dust properties sufficiently well, and/or that the adopted uncertainties are underestimated. The H and K_s bands dominate the high χ² values: relative to the observational uncertainties the residuals in H and K_s are on average six times larger than in the other bands. As noted before, the extinction law in the NIR might be imperfect. However, the data that we have do not allow us to improve on this and we therefore accept the values of our fit.

We find a strong gradient in extinction A₅₄₉₅ across the core of R136. We discuss this in more detail in Sects. 4.3 and 5.1; for now it is important to note that the wide range of values in R₅₄₉₅ that we find throughout the core of R136 correspond to similar values of c₂, for which we find values in the range of c₂ = 0.78 – 2.36, with an average of c₂ = 1.30 ± 0.22. The values show no significant trend as a function of R₅₄₉₅, contrary to what is observed for Galactic sightlines (Eq. (3)).

We can now construct an average R5495-dependent extinction law towards R136. We do this by adopting the R₅₄₉₅-dependent optical and NIR law MA14 (using the spline point parameteri-sation of F99, as described in Sect. 3.3), and adopting the UV parameterisation of Fitzpatrick & Massa (1990), assuming average values of c₂ and c₄ that we derive from the fitting, and 2175 Å feature parameters of Gordon et al. (2003). We discuss how this curve compares to other LMC extinction curves in Sect. 5.3. The parameters of the extinction curve towards R136 are summarised in Table 3; a Python implementation of the law is presented in Appendix C, and can also be found on Github⁷.

Fig. 4 Example fits of the stars R136a3 (WN5h; left) and H55 (O2 V((f*))z; right). In the upper part of each panel, we show the photometric points that were used for the fitting; with red solid squares and yellow filled circles indicating points in the UV and optical or NIR, respectively. The thick solid blue line in the background shows the UV spectroscopy from which the UV photometric points were derived. The light blue dotted and dark blue dashed lines in the top panels show the adopted intrinsic SED, and the reddened SED (resulting from the fitting process), respectively. The lower panel shows residuals, and the best-fit parameters are printed in the upper panel. Lastly, each panel contains an inset where the shape of the extinction curve under consideration (solid orange line) is compared to the Galactic curve of Cardelli et al. (1989, RV = 3.1; black dashed line).

4.2 Extinction towards the outskirts of R136

Upon plotting the slope in the optical and NIR, captured by V – I, versus the extinction A₅₄₉₅, we see a strong correlation between the two (Fig. 5): (6)

Under the assumption that the dust properties of sightlines towards the outskirts of R136 are similar to those of sightlines towards the core, we can use this relation to estimate A₅₄₉₅ for 1657 sources with V < 19 in the catalogue of De Marchi et al. (2011) for which both V and I magnitudes are available. The magnitude cut-off ensures that we exclude the vast majority of the pre-main sequence stars included in the De Marchi et al. (2011) catalogue (see their Fig. 8). For the 1657 stars for which we estimate the extinction in this way, we find an average of A₅₄₉₅ = 1.88 ± 0.94. The spatial distribution of the obtained extinction values is presented in Fig. 6, and shows that the trend in A₅₄₉₅ that was observed throughout the cluster core continues east of the cluster core. We discuss the trends in extinction towards R136 relative to the larger field in the following subsection and in Sect. 5.1.

Fig. 5 Extinction A5495 versus observed colour V – I. Each diamond corresponds to a star in the core of R136; light-blue error bars indicate 1σ uncertainties. The red solid line is the best fit through all points, and in yellow we show the bootstrapped uncertainties on the linear fits, represented by linear fits to 1000 randomly chosen samples.

Fig. 6 Extinction map of the R136 cluster and surroundings. The dots indicate stars of the De Marchi et al. (2011) sample for which AV was estimated using V – I (see Sect. 4.2), as well as stars in the core of the R136 cluster, for which AV was determined with a full SED fit. The core of the cluster is indicated with a black circle, and a zoom onto the core of the R136 cluster is shown in the inset in the upper left corner. We note that in order to make more details visible, the dots in the inset plot are of a smaller size relative to the background image than those of the main panel. Background image credits: ESO/R. Fosbury (ST-ECF), R. O'Connell (University of Virginia, Charlottesville), and the Wide Field Camera 3 Science Oversight Committee (colours adapted).

4.3 Spatial trends in extinction

We find a significant variation in extinction across the core of R136. What we refer to as the ‘core’ spans a region of about 4′′ in diameter, corresponding to 1 pc for the LMC distance. Panel a of Fig. 7 shows a spatial map of extinction in the core. A gradient from east to west is clearly visible, where the east side of the cluster is more extincted than the west side. Inspecting the stars of the De Marchi et al. (2011) sample in Fig. 6, we see that this trend extends to outside the core.

The extinction gradient is visualised quantitatively in Fig. 8, where A₅₄₉₅ values of the core stars are plotted as a function of right ascension (a higher value means more to the east) and declination (a lower value means more to the south). We see a particularly strong relation as a function of right ascension, but there seems to be no significant trend as a function of declination. We can express the extinction A₅₄₉₅ across the core of R136 as a function of spatial coordinates as follows: (7)

where RA is the right ascension expressed in decimal degrees; as we do not find a significant trend as a function of declination, Eq. (7) only depends on right ascension.

Panel b of Fig. 7 shows another map of the core but this time with the values for R₅₄₉₅ overplotted. We see a spatial trend for this quantity too, where we find higher values of R₅₄₉₅ towards the east of the cluster, although this trend is less clear than in the case of A₅₄₉₅. Panel c of Fig. 7 shows that there is a strong correlation between R5495 and A_V, and that the relation between the two implies A₄₄₀₅ – A₅₄₉₅ ≈ 0.4. We note that if we refit the SED while adopting a fixed value of R₅₄₉₅ (equal for all stars), we still recover the spatial trend for A₅₄₉₅, although the gradient is slightly less steep. We discuss possible causes of the spatial trend in Sect. 5.1.

4.4 The UV curvature (c₄)

For the parameter c₄, describing the curvature of the extinction curve in the far-UV (x > 5.9 µm⁻¹ or λ < 1695 Å), we find an average value of c₄ = 0.09 ± 0.08. This value is lower than found by Misselt et al. (1999) and Gordon et al. (2003) for the supergiant shell LMC SGS 2 near 30 Doradus; these authors find c₄ = 0.42 ± 0.08 and c₄ = 0.29 ± 0.06, respectively. A low value of c₄ corresponds to a weak far-UV curvature.

We note that the value of c₄ we find for the R136 core stars is possibly even lower than the value quoted above. This is due to possible (relative) flux-calibration issues of the STIS spectra, as revealed in Fig. 1 : upon comparing the integrated flux of the R136 core as recorded by STIS with the flux as recorded by GHRS, the two start diverging at around λ < 1400 Å. In this wavelength regime, which is especially sensitive to the value of c₄, the STIS flux is lower than the GHRS flux, namely by a factor of about 0.88. If the low fluxes are due to an unknown calibration issue of the STIS spectra, and in reality the shape of the spectra is more like that recorded by GHRS, then in our analysis we have overestimated c₄. In order to match the reddened model fluxes to the GHRS flux (i.e., for the reddened model fluxes to be a factor 1/0.88 = 1.14 higher than the STIS spectra, at λ = 1280 Å), the value of c₄ would have to approach zero, that is, no far-UV curvature at all.

5 Discussion

5.1 Extinction due to an extension of the Stapler Nebula

We find a strong gradient in the amount of extinction across the 1 pc core of the R136 cluster, ranging from A₅₄₉₅ ≈ 1.0 in the west, up to A₅₄₉₅ ≈ 2.6 east of the cluster core. The estimated extinction of the stars in the outskirts of R136 suggests that this trend extends outside the cluster core (Fig. 6). We compare the extinction gradient with images of the cluster in order to identify the larger-scale structure of the intervening gas and dust. In the composite optical and UV HST image of the cluster and surroundings (Fig. 9, top row), we see a dark cloud to the northeast of the cluster. The cloud is an extension of the Stapler Nebula; Kalari et al. (2018) study the nebula and its relation to R136 in detail. The fact that there are hardly any stars visible in front of this cloud suggests that it is situated in the foreground (Kalari et al. 2018; Wong et al. 2022). This cloud is even more clearly visible in the IR images captured by NIRCam and MIRI on the James Webb Space Telescope (JWST, Fig. 9, bottom row). These high-resolution images reveal the detailed structure of the cloud, with one arm of the cloud extending all the way to the easternmost star of our sample, for which we find a relatively high extinction. Nonetheless, seen in the NIR, optical, and UV, the west edge of the cloud has a projected distance from the centre of R136 of ≈1–2 pc, and these images therefore do not provide direct evidence that the cloud is responsible for the extinction gradient across the cluster.

However, images at even longer wavelengths reveal that the molecular cloud stretches out farther than can be seen in the HST and JWST images. Wong et al. (2022) observed 30 Doradus with ALMA, and map the strength of the CO(2−1) rotational line, a tracer of cold molecular gas. With high spatial resolution and sensitivity, ALMA (beam size: 1′′75) allows us to distinguish details of the distribution of the cold gas. The cluster core spans about 2.3 × the beam size, allowing us to distinguish between the ¹³CO column density in the east versus the west side of the cluster. Figure 10 shows a map of the column density of¹³ CO, as well as the positions of the R136 core stars, marked by their extinction. It is clear that in regions where¹³ CO is detected (eastern half of the cluster core), the extinction is higher, suggesting that this molecular cloud contains the dust grains that are causing the extinction gradient. Furthermore, we see that the position of the highest ¹³CO column density (≥10¹⁵ cm⁻²) corresponds to the position of the dark cloud in Fig. 9 (yellow-green dashed line). The molecular gas mapped by the CO(2−1) line therefore seems to be associated with the dark cloud of Fig. 9.

Having addressed the spatial trend in A_V, we now turn to the trend we find for R₅₄₉₅. On average, we find higher values of R₅₄₉₅ towards the east of the cluster than towards the west side. Higher R₅₄₉₅ values are interpreted as the result of fewer very small grains (having sizes <250 Å) of either silicate particles (Xiang et al. 2017) or a mixture of graphitic and silicate grains (Weingartner & Draine 2001), relative to larger grains. As R₅₄₉₅ has a strong positive correlation with A_V (see panel c) of Fig. 7), we can see from Fig. 10 that higher values of R₅₄₉₅ map to higher values of¹³ CO column density. This could indicate that grain growth through aggregation is more efficient in the denser environments of the molecular cloud, or alternatively that denser cloud regions shield the grains more efficiently from strong radiation fields that may break apart dust particles.

Very limited information is available on the relation between A₅₄₉₅ and R₅₄₉₅ in other star forming regions in the LMC. Maíz-Apellániz & Barbá (2018) study the extinction towards O-type stars in Galactic HII regions, and find results that are not in line with ours. These authors also find differences in R₅₄₉₅ on relatively small spatial scales (≈1–8 pc, A₅₄₉₅ = 2.0–7.3), but higher values of R₅₄₉₅ do not always coincide with higher values of A₅₄₉₅. On the contrary, combining the measurements of four different HII regions (shown in their Fig. 7), a tentative opposite trend appears, with larger values of R₅₄₉₅ generally corresponding to lower A₅₄₉₅.

Fig. 7 Spatial trends in AV and R5495 in the core of R136. Panels a and b show maps of A5495 and R5495 overplotted on HST/WFC3 V-band (F555W) photometry (O’Connell 2010). The colour of each hexagon corresponds to the average value of the quantity in that region of the core, with red tints corresponding to higher values of A5495 and R5495 and the yellow tints to lower values. A clear spatial trend is visible for A5495, with higher values of A5495 towards the east of the cluster (see also Fig. 8). A similar, albeit weaker trend is visible for R5495. Panel c shows that there is a strong correlation between A5495 and R5495 (colours as in Fig. 5). The dotted lines in the background indicate constant values of A4405 – A5495 ranging from 0.2 to 0.8; the obtained relation implies A4405 – A5495 ≈ 0.4.

Fig. 8 Extinction of the stars in the core of R136 (A5495) plotted against their right ascension (left) and declination (right). Colours as in Fig. 5.

Fig. 9 Multi-wavelength view of the cluster R136 and surroundings. Panels a) and b) show the cluster imaged in optical and UV with WFC3/HST. The wider area around R316 is shown on the left, on the right we zoom in on the core of the cluster. The extinction map is projected onto the zoomed image in yellow to red colours (as in Fig. 7). Imaged in optical and UV, the dark cloud north-east of the cluster, an extension of the Stapler Nebula, appears to have a considerable projected distance from the cluster core (≈1–2 pc). The green-yellow dashed line indicates the contour of a 13CO column density of 1015 cm−2 (see Fig. 10). Panels c)–e) are as panel b, but with different background images: a NIR image captured with NIRCam on James Webb Space Telescope (JWST), a mid-infrared (MIR) image captured with MIRI on JWST, and a composite NIR/(sub-)mm image, captured in NIR with HAWK- I/VLT and VISTA, and in (sub-)mm wavelengths with ALMA. In panel c) the orange-brown colours show cold gas corresponding roughly with the dark cloud in the HST image. Moving to longer wavelengths in panel c), the stars fade and the cool gas consisting of hydrocarbons is lighting up (turquoise). In panel e) light pink regions correspond to relatively hot gas (NIR), and red-yellow areas indicate the presence of cold, dense gas (ALMA). It seems possible that the extinction gradient in R136 may be associated with a lower density fringe of the extension of the Stapler Nebula; see also Fig. 10 for additional evidence. Image credits (from left to right): NASA, ESA, F. Paresce (INAF-IASF, Bologna, Italy), R. O’Connell (University of Virginia, Charlottesville), and the Wide Field Camera 3 Science Oversight Committee; NASA, ESA, CSA, and STScI; IMAGE: NASA, ESA, CSA, STScI, Webb ERO Production Team; ESO, ALMA (ESO/NAOJ/NRAO)/Wong et al., 2022 ESO/M.-R. Cioni/VISTA Magellanic Cloud survey, Acknowledgment: Cambridge Astronomical Survey Unit.

5.2 Revised WNh star masses

When accounting for the new values of the extinction, we find changes in luminosity of the R136 core stars of on average 0.03 dex compared to the luminosities obtained by Brands et al. (2022); the revisions in log L/L_⊙ for individual stars range from −0.13 to +0.04 dex. The changes in luminosity imply small changes of the stellar masses. For the O-stars in the core of R136, we estimate the mass-reduction by eye using their positions on the Hertzsprung–Russell diagram (HRD) and the tracks of Brott et al. (2011) and Köhler et al. (2015)⁸. We estimate that the stellar mass is typically changed by 5–10%.

We performed a detailed analysis of the masses of the most massive stars in the sample, namely R136a1, R136a2, and R136a3. We do this in the same manner as in Brands et al. (2022), only now with our new luminosity values. For this, we use the Bayesian tool BONNSAI⁹ (Schneider et al. 2017), in combination with the evolutionary models of Brott et al. (2011) and Köhler et al. (2015). BONNSAI allows the comparison of observed stellar parameters with stellar evolution models in order to infer posterior distributions of model parameters such as initial and current stellar mass. For our derivations of stellar masses and ages we use temperature, helium abundance, and surface gravity as derived by Brands et al. (2022), and the revised luminosity from this work. As priors, we adopt the Salpeter (1955) initial mass function and the rotation distribution of Ramírez-Agudelo et al. (2013).

We find that the evolutionary and initial masses we derive for R136a1 and R136a3 agree within uncertainties with those of Brands et al. (2022), while for R136a2 our masses are 25% lower. The new age that we derive for R136a3 agrees well with that derived by Brands et al. (2022), and for R136a1 and R136a2 we find downward and upward revisions for the age of about 20%, respectively. We note that Brands et al. (2022) provide two sets of values for the stellar parameters; we adopt the values of their optical + UV fits. The new values can be found in Table 4. Comparing our initial masses with literature values, we find that they are generally higher than those derived by Rubio-Díez et al. (2017) and Kalari et al. (2022), although they do agree with the latter within errors. We note that initial masses of these stars, both ours and those presented in literature, are highly model dependent (see, e.g. Gräfener 2021, Higgins et al. 2022, and Brands et al. 2022, their Fig. 17). Moreover, even within a given set of evolution models, the derived initial mass is sensitive to the combination of observables used for the comparison with the evolutionary models. Lastly, all these analyses are based on single-star evolution scenarios, but it is important to keep in mind that these very massive stars might be merger products (e.g. Banerjee et al. 2012).

Fig. 10 13CO column density map of R136 and surroundings, obtained with ALMA (Wong et al. 2022). Darker regions indicate a higher 13CO column density; in white regions no 13CO was detected. Contours, created with DS9 (smoothing = 4), correspond to lines of constant13CO column density, the value of which is indicated in the plot in units of log[cm−2]. The dashed circle in the lower left corner of the image indicates the beam size. Small filled circles indicate the positions of the stars in the R136 core; the colour of each circle corresponds to the measured extinction, with darker colours corresponding to a higher extinction. Comparing the extinction of the R136 core stars with the 13CO map, we see that in regions where13 CO is detected we measure stronger extinction.

5.3 Anomalous extinction towards R136?

The average R_V value we find for R136 (R_V ≈ R₅₄₉₅ = 4.38 ± 0.87) is in good agreement with previously obtained averages for the region¹⁰ (Crowther et al. 2016; Bestenlehner et al. 2020, for R136, and Doran et al. 2013; Bestenlehner et al. 2014; Maíz-Apellániz et al. 2014; De Marchi & Panagia 2014; De Marchi et al. 2016, for the wider 30 Doradus region). Such high values of R_V are, in the Milky Way, associated with low UV extinction, but this is not what we observe for R136. The UV extinction for R_V = 4.4 in R136 is far higher than the UV extinction in the Galactic case for the same R_V (Fig. 11). From a physical point of view, it is not unexpected that the slope in the optical (R_V) and the strength of UV extinction can vary independently, as these two parts of the curve are thought to be associated with complementary dust particle populations (e.g. Weingartner & Draine 2001; Xiang et al. 2017). Nevertheless, one population can be related to the other, as in the Milky Way their relative contributions do not vary randomly but follow a relatively simple relation (i.e. Eq. (3)). This relation is dictated by the properties of the interstellar dust, which in turn are affected by the environment. In 30 Doradus, or at least in R136, we do not observe this relation between R_V and the UV extinction. Possibly, the dust in and near 30 Doradus is affected by the intense radiation and mechanical energy that is being deposited in the interstellar medium by the large population of hot stars in the region, as well as by previous powerful supernova explosions of local massive stars (De Marchi & Panagia 2019).

In Fig. 11, we compare our Rv-dependent extinction law towards R136 with other extinction curves that are tailored to regions of the LMC. Contrasting the average curve towards R136 with the curves of Gordon et al. (2003), it appears that, similar to Galactic sightlines, higher values of R_V correspond to lower values of UV extinction. This would suggest that a relation between R_V and UV extinction may exist, but that it differs from the Galactic relation. The fact that we do not observe this relation within R136 could be related to the uncertainties in the UV flux calibration, which increase the scatter on our measured c₂ values. On the other hand, Howarth (1983), who also derive an LMC average curve and obtain R_V = 3.1 for their sample, find a UV extinction that is too strong with respect to their R_V to match the tentative relation between R_V and UV extinction that the other curves in Fig. 11 might suggest. Furthermore, the result of Howarth (1983) appears to be in contradiction with the results of Gordon et al. (2003).

The differences between the ‘LMC average’ curves of Gordon et al. (2003) and Howarth (1983) can likely be understood by the fact that their samples consist of different sightlines. Howarth (1983, and also Pei 1992, who find a very similar average curve,) use the samples of Nandy et al. (1980,1981), and Koornneef & Code (1981), of which about half of the sightlines lie around 30 Doradus, with the other half being spread out over other regions in the LMC. For their LMC average curve, Gordon et al. (2003) use sightlines that are more clustered and lower in number. The sightlines of the different LMC samples are shown in Fig. 12. The sample used to derive the curve towards LMC SGS 2 is also indicated, as are R136 and 30 Doradus. We note that while sightlines near 30 Doradus are considered in all samples, no sightline towards 30 Doradus itself is included in any of the samples; the present work on R136 is therefore unique in this respect.

In any case, we cannot draw firm conclusions on the R_V dependence of UV extinction with only four curves. It therefore remains unclear as to whether or not there exists a relation between the different grain populations in the LMC, or what the nature of this relation could be. It would be of value to carry out an LMC-wide study of the R_V dependence of UV extinction in order to investigate this further.

Fig. 11 Comparison of extinction curves tailored to (specific regions within) the LMC. Shown are the average extinction curve of Howarth (1983, green dotted line), which is nearly equal to that of Pei (1992, not shown) and two curves of Gordon et al. (2003), based on samples within LMC SGS 2 (blue dash-dotted line), and a sample in other parts of the LMC (‘LMC avg.’, orange dashed line). We compare these curves to the average curve derived in this work, which we show for R5495 = 4.4 (the average value towards R136; red solid line). For reference, we also show the Galactic curve of Fitzpatrick (1999) for different values of RV (grey lines). The grey areas indicate parts of the SED that we analysed in this work. The white area contains the 2175 Å feature, a wavelength region that was not covered by our data; the extinction curves are also shown in this wavelength range for completeness, but we highlight the fact that the bump parameters of our R136 curve are simply taken from the LMC SGS 2 curve of Gordon et al. (2003).

Fig. 12 Position of stars of different samples that were used to derive extinction curves towards the LMC. The samples of Nandy et al. (1980, 1981) and Koornneef & Code (1981, yellow circles) were used to construct the curves of Howarth (1983) and Pei (1992). Gordon et al. (2003) split their sample into two and derive a separate curve for each, referring to one sample as LMC2 (LMC SGS 2, blue squares) and the other as ‘LMC Average’ (red diamonds). The green triangle denotes the position of R136. This figure was made with use of the Aladin Sky Atlas (Bonnarel et al. 2000); the background consists of DSS2 colour images.

6 Conclusion

Employing the extinction-without-standards method, we infer NIR to UV extinction characteristics towards 50 stars in the core of R136. On average, we find an extinction of A_V ≈ A₅₄₉₅ = 1.70 ± 0.45. However, we infer a strong spatial gradient in extinction properties across the cluster core, where the extinction in the east is about one magnitude higher than the extinction in the west of the cluster. Comparing our extinction map to multi-wavelength observations of the same region, we conclude that the observed extinction gradient is likely caused by material belonging to an extension of a molecular cloud called the Stapler Nebula, which lies to the northeast of the cluster and stretches all the way to the R136 core.

In line with previous studies, we obtain a relatively high average value of R_V ≈ R₅₄₉₅ = 4.38 ± 0.87 towards R136. Moreover, we find that the UV extinction towards R136 is significantly stronger than the canonical Galactic extinction at the same value for the same R_V, implying a relatively large fraction of small particles near R136. The intense radiation field and mechanical energy that is being deposited in the 30 Doradus interstellar medium by the hot stars and their powerful core-collapse supernovae could play a role in this process. A consequence of the stronger UV absorption is that less ionising photons can escape. At A_V = 1.0, the extinction towards R136 at UV wavelengths in the range λ ≈ 1700–1250 Å is about one magnitude higher compared to the canonical Galactic extinction at the same R_V and A_V, implying that the fraction of ionising photons that can escape is a factor 2.5 lower ().

We have now investigated the relation between R_V and extinction in the UV for one region in the LMC, namely R136. Extending this investigation to different environments throughout the Magellanic Clouds would be an interesting topic for future studies: knowledge about the interdependence of different dust populations could provide insights into the environmental factors determining the properties and evolution of the dust therein. This would be of particular interest in the context of starburst galaxies, where stellar populations are unresolved and dust properties need to be known in order to interpret observations. This includes the role of dust in star-bursting regions in absorbing and reprocessing of ionising photons trying to escape to inter-galactic space.

Acknowledgements

We thank the referee Jesús Maíz-Apellániz for providing constructive comments. This publication is part of the project ‘Massive stars in low-metallicity environments: the progenitors of massive black holes’ with project number OND1362707 of the research TOP-programme, which is (partly) financed by the Dutch Research Council (NWO). This research is based in part on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. These observations are associated with programme 12465, and with Early Release Science observations made by the WFC3 Scientific Oversight Committee. This paper makes use of the following ALMA data: ADS/JAO.ALMA 2019.1.00843.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), Observatoire de Paris/LESIA (Paris), and Observatoire de Lyon (OSUL/CRAL). This work is based in part on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with programme 2729. This research made use of DS9, a tool for data visualisation supported by the Chandra X-ray Science Center (CXC) and the High Energy Astrophysics Science Archive Center (HEASARC) with support from the JWST Mission office at the Space Telescope Science Institute for 3D visualisation. This research has made use of “Aladin sky atlas” developed at CDS, Strasbourg Observatory, France (Bonnarel et al. 2000). This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France (Wenger et al. 2000).

Appendix A UV bands

Table A.1 lists the wavelength ranges that we used for binning our UV flux measurements. We use the same ranges for all stars. The average flux of each wavelength interval is the value that was used for the fitting. This is equivalent to constructing a (synthetic) passband with a transmission of 1.0 in the wavelength range corresponding to the filter, and 0.0 outside that range.

Appendix B Best-fit values for individual sources

The best-fit values for individual sources in the core of R136 can be found in Table B.1.

Appendix C The extinction law towards R136

The R₅₄₉₅-dependent spline points that we used for the optical and NIR part of the law can be found in Table 2. A Python function for the computation of the R₅₄₉₅-dependent extinction law towards R136 is presented in Listing 1 and can also be found on Github¹¹. This function contains both the optical and NIR part as derived by MA14, but parameterised in the format of F99, and the modified UV part (Table 3).

Listing 1 A Python function that can be used for the computation of the R5495-dependent extinction law towards R136.

References

All Tables

All Figures

Fig. 1 Integrated UV flux of the core of R136. The upper panel shows the GHRS spectrum of the inner 2′ × 2′ (Heap et al. 1992; light blue crosses), as well as the summed STIS spectra used in this work (dark blue circles), with both covering roughly the same inner region of the cluster. The bottom panel shows the ratio of the former two spectra (yellow dots) as well as the central wavelength of the synthetic UV bands used in this work (red triangles, see Sect. 3.1). The dashed lines correspond to ratios of 0.75 and 1.25.

In the text

	Fig. 2 Example of an observed SED (source: R136a3) covering the UV (dark blue circles), optical (light blue squares), and NIR (turquoise triangles). The UV flux points are derived from flux-calibrated STIS spectroscopy (grating: G140L), and the optical and NIR fluxes come from broadband photometry; the (equivalent) name of each band is indicated.
In the text

Fig. 3 Example fits of SEDs of the stars R136a3 (left; WNh) and H55 (right; O2 V((f*))z), adopting the standard Galactic law of F99. The figure shows that this law is not suitable for our analysis: the Galactic dependence of UV extinction on Rv does not hold for sightlines towards R136. In other words, the reddened flux cannot be modelled accurately over the full wavelength range: an RV based on the optical and NIR photometry gives a total mismatch with the observed UV fluxes, whereas a much lower value – that can reproduce the UV fluxes – gives a very poor fit to the observed slope in the optical. The top panels display the flux from UV to NIR; the bottom panels zoom in on the optical regime for which UBVI photometry is available. Circles and squares indicate observed fluxes in the UV, and optical or NIR, respectively. The best-fitting reddened model including the UV is shown in orange; the best fit excluding the UV constraints in red. For this example, we used the law of F99, but the same behaviour is seen with other RV -dependent Galactic laws.

In the text

Fig. 4 Example fits of the stars R136a3 (WN5h; left) and H55 (O2 V((f*))z; right). In the upper part of each panel, we show the photometric points that were used for the fitting; with red solid squares and yellow filled circles indicating points in the UV and optical or NIR, respectively. The thick solid blue line in the background shows the UV spectroscopy from which the UV photometric points were derived. The light blue dotted and dark blue dashed lines in the top panels show the adopted intrinsic SED, and the reddened SED (resulting from the fitting process), respectively. The lower panel shows residuals, and the best-fit parameters are printed in the upper panel. Lastly, each panel contains an inset where the shape of the extinction curve under consideration (solid orange line) is compared to the Galactic curve of Cardelli et al. (1989, RV = 3.1; black dashed line).

In the text

	Fig. 5 Extinction A5495 versus observed colour V – I. Each diamond corresponds to a star in the core of R136; light-blue error bars indicate 1σ uncertainties. The red solid line is the best fit through all points, and in yellow we show the bootstrapped uncertainties on the linear fits, represented by linear fits to 1000 randomly chosen samples.
In the text

Fig. 6 Extinction map of the R136 cluster and surroundings. The dots indicate stars of the De Marchi et al. (2011) sample for which AV was estimated using V – I (see Sect. 4.2), as well as stars in the core of the R136 cluster, for which AV was determined with a full SED fit. The core of the cluster is indicated with a black circle, and a zoom onto the core of the R136 cluster is shown in the inset in the upper left corner. We note that in order to make more details visible, the dots in the inset plot are of a smaller size relative to the background image than those of the main panel. Background image credits: ESO/R. Fosbury (ST-ECF), R. O'Connell (University of Virginia, Charlottesville), and the Wide Field Camera 3 Science Oversight Committee (colours adapted).

In the text

Fig. 7 Spatial trends in AV and R5495 in the core of R136. Panels a and b show maps of A5495 and R5495 overplotted on HST/WFC3 V-band (F555W) photometry (O’Connell 2010). The colour of each hexagon corresponds to the average value of the quantity in that region of the core, with red tints corresponding to higher values of A5495 and R5495 and the yellow tints to lower values. A clear spatial trend is visible for A5495, with higher values of A5495 towards the east of the cluster (see also Fig. 8). A similar, albeit weaker trend is visible for R5495. Panel c shows that there is a strong correlation between A5495 and R5495 (colours as in Fig. 5). The dotted lines in the background indicate constant values of A4405 – A5495 ranging from 0.2 to 0.8; the obtained relation implies A4405 – A5495 ≈ 0.4.

In the text

	Fig. 8 Extinction of the stars in the core of R136 (A5495) plotted against their right ascension (left) and declination (right). Colours as in Fig. 5.
In the text

Fig. 9 Multi-wavelength view of the cluster R136 and surroundings. Panels a) and b) show the cluster imaged in optical and UV with WFC3/HST. The wider area around R316 is shown on the left, on the right we zoom in on the core of the cluster. The extinction map is projected onto the zoomed image in yellow to red colours (as in Fig. 7). Imaged in optical and UV, the dark cloud north-east of the cluster, an extension of the Stapler Nebula, appears to have a considerable projected distance from the cluster core (≈1–2 pc). The green-yellow dashed line indicates the contour of a 13CO column density of 1015 cm−2 (see Fig. 10). Panels c)–e) are as panel b, but with different background images: a NIR image captured with NIRCam on James Webb Space Telescope (JWST), a mid-infrared (MIR) image captured with MIRI on JWST, and a composite NIR/(sub-)mm image, captured in NIR with HAWK- I/VLT and VISTA, and in (sub-)mm wavelengths with ALMA. In panel c) the orange-brown colours show cold gas corresponding roughly with the dark cloud in the HST image. Moving to longer wavelengths in panel c), the stars fade and the cool gas consisting of hydrocarbons is lighting up (turquoise). In panel e) light pink regions correspond to relatively hot gas (NIR), and red-yellow areas indicate the presence of cold, dense gas (ALMA). It seems possible that the extinction gradient in R136 may be associated with a lower density fringe of the extension of the Stapler Nebula; see also Fig. 10 for additional evidence. Image credits (from left to right): NASA, ESA, F. Paresce (INAF-IASF, Bologna, Italy), R. O’Connell (University of Virginia, Charlottesville), and the Wide Field Camera 3 Science Oversight Committee; NASA, ESA, CSA, and STScI; IMAGE: NASA, ESA, CSA, STScI, Webb ERO Production Team; ESO, ALMA (ESO/NAOJ/NRAO)/Wong et al., 2022 ESO/M.-R. Cioni/VISTA Magellanic Cloud survey, Acknowledgment: Cambridge Astronomical Survey Unit.

In the text

Fig. 10 13CO column density map of R136 and surroundings, obtained with ALMA (Wong et al. 2022). Darker regions indicate a higher 13CO column density; in white regions no 13CO was detected. Contours, created with DS9 (smoothing = 4), correspond to lines of constant13CO column density, the value of which is indicated in the plot in units of log[cm−2]. The dashed circle in the lower left corner of the image indicates the beam size. Small filled circles indicate the positions of the stars in the R136 core; the colour of each circle corresponds to the measured extinction, with darker colours corresponding to a higher extinction. Comparing the extinction of the R136 core stars with the 13CO map, we see that in regions where13 CO is detected we measure stronger extinction.

In the text

Fig. 11 Comparison of extinction curves tailored to (specific regions within) the LMC. Shown are the average extinction curve of Howarth (1983, green dotted line), which is nearly equal to that of Pei (1992, not shown) and two curves of Gordon et al. (2003), based on samples within LMC SGS 2 (blue dash-dotted line), and a sample in other parts of the LMC (‘LMC avg.’, orange dashed line). We compare these curves to the average curve derived in this work, which we show for R5495 = 4.4 (the average value towards R136; red solid line). For reference, we also show the Galactic curve of Fitzpatrick (1999) for different values of RV (grey lines). The grey areas indicate parts of the SED that we analysed in this work. The white area contains the 2175 Å feature, a wavelength region that was not covered by our data; the extinction curves are also shown in this wavelength range for completeness, but we highlight the fact that the bump parameters of our R136 curve are simply taken from the LMC SGS 2 curve of Gordon et al. (2003).

In the text

Fig. 12 Position of stars of different samples that were used to derive extinction curves towards the LMC. The samples of Nandy et al. (1980, 1981) and Koornneef & Code (1981, yellow circles) were used to construct the curves of Howarth (1983) and Pei (1992). Gordon et al. (2003) split their sample into two and derive a separate curve for each, referring to one sample as LMC2 (LMC SGS 2, blue squares) and the other as ‘LMC Average’ (red diamonds). The green triangle denotes the position of R136. This figure was made with use of the Aladin Sky Atlas (Bonnarel et al. 2000); the background consists of DSS2 colour images.

In the text

	Listing 1 A Python function that can be used for the computation of the R5495-dependent extinction law towards R136.
In the text