Extinction towards the cluster R136 in the Large Magellanic Cloud An extinction law from the near-infrared to the ultraviolet

Context. The cluster R136 in the giant star-forming region 30 Doradus in the Large Magellanic Cloud (LMC) offers a unique opportunity to resolve a stellar population in a starburst-like environment. Knowledge of the extinction towards this region is key for the accurate determination of stellar masses, and for the correct interpretation of observations of distant, unresolved starburst galaxies. Aims. Our aims are to construct an extinction law towards R136, and to measure the extinction towards individual sources inside the cluster. This will allow us to map the spatial distribution of the dust, to learn about dust properties, and to improve mass measurements of the very massive WNh stars inside the cluster. Methods. We obtain the near-infrared to ultraviolet extinction towards 50 stars in the core of R136, employing the ‘extinction without standards’ method. To assure good fits over the full wavelength range, we combine and modify existing extinction laws. Results. We detect a strong spatial gradient in the extinction properties across the core of R136, coinciding with a gradient in density of cold gas that is part of an extension of the Stapler Nebula, a molecular cloud lying northeast of the cluster. In line with previous measurements of R136 and the 30 Doradus region, we obtain a high total-to-relative extinction ( R V = 4 . 38 ± 0 . 87). However, the high values of R V are accompanied by relatively strong extinction in the ultraviolet, contrary to what is observed for Galactic sightlines. Conclusions. The relatively strong ultraviolet extinction towards R136 suggests that the properties of the dust towards R136 differ from those in the Milky Way. For R V ∼ 4 . 4, about three times fewer ultraviolet photons can escape from the ambient dust environment relative to the canonical Galactic extinction at the same R V and A V . Therefore, if dust in the R136 star-bursting environment is characteristic for cosmologically distant star-bursting regions, the escape fraction of ultraviolet photons from such regions is overestimated by a factor of three relative to the standard Milky Way assumption for the total-to-selective extinction. Furthermore, a comparison with average curves tailored to other regions of the LMC shows that large differences in ultraviolet extinction exist within this galaxy. Further investigation is required in order to decipher whether or not there is a relation between R V and ultraviolet extinction in the LMC.


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 4 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 stel-lar 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. (1994Calzetti et al. ( , 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;Maíz Apellániz et al. 2014;De Marchi & Panagia 2014;De Marchi et al. 2016;Crowther et al. 2016;Bestenlehner et al. 2020), 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 R V 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 Section 2, we describe our sample and data, and in Section 3 we provide details of the methods used.In Section 4, we present our main findings, which we discuss in a broader context in Section 5. Finally, we summarise the main conclusions of our study in Section 6.

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" x 0.2") contiguous pointings with grating G140L.In each slit, multiple sources are present (see Crowther et al. 2016, andBrands et al. 2022, their Figs. 1).The spectra are extracted using multispec, a package tailored to extracting spectra from crowded regions (Maiz-Apellaniz 2005, 2007).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 HST/STIS G140L (1160  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 x 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 Section 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(Khorrami et al. , 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.

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 Section 4.2.

Observed SEDs
We compile observed SEDs spanning from the UV to NIR (see also Section 2).The optical and NIR observed magnitudes are converted to fluxes using filter transmission curves taken from the SVO Filter Profile Service1 .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 cm2 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.

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 2 = −0.05)compared to the case where the continuum is unmodified.This change is within the typical range for uncertainties of individual c 2 values.Moreover, the best-fit values for A 5495 , R 5495 , 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 −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 Section 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 * .

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 Section 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) 2 .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, The figure shows that this law is not suitable for our analysis: the Galactic dependence of UV extinction on R V does not hold for sightlines towards R136.In other words, the reddened flux cannot be modelled accurately over the full wavelength range: an R V 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 bestfitting 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 R V -dependent Galactic laws.
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 5495 ≡ A 5495 /(A 4405 − A 5495 ), the monochromatic equivalent of R V .We therefore use monochromatic quantities throughout the paper, with the exception of Section 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 5495 and R 5495 ) 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 point3 , so that the dependence on R 5495 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 −1 (λ ≥ 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 −1 (λ = 3304 Å), and the R 5495 -dependent extinction at this point we take from the law of MA14.The latter spline point is crucial in order to fit R 5495 , 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 parameterisation 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 , the UV part of the law of F99 has the following form: where x ≡ 1/λ µm −1 , a 1 = 0.5392, and a 2 = 0.056444 .The parameters c 1 and c 2 relate to a linear background (see below); c 3 indicates the strength of the 2175 Å feature (see below), and c 4 and c 5 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 = 5.9; we do the same.We leave c 4 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: with x 0 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 3 = 1.463, γ = 0.945, and x 0 = 4.558.Lastly, and key for our study, we discuss the parameters c 1 and c 2 .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 2 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 2 ) and the slope of extinction in the optical (R 5495 ) will be fitted independently.While we fit the slope of the UV extinction, the equation for c 1 we leave unchanged.This is because c 2 and c 1 are (to some extent) degenerate, and we do not have enough data to break this degeneracy.Experiments where we leave c 1 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 5495 (slope in the optical; monochromatic equivalent of R V ), c 2 (strength of UV extinction), and c 4 (far-UV curvature) are free parameters.

Fitting
For each individual star, we optimise the free parameters of the adopted extinction law (Section 3.3) in order to find the best match between the observed SED and the reddened model SED.We fit five parameters: R 5495 , affecting the shape of the optical and NIR part of the curve; c 2 and c 4 , 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.59kpc 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 Section 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 5 , 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 5495 , c 2 , and c 4 ), we use the minimize function of the Python package lmfit 6 (Newville et al. 2014).We adopt the least-squares Levenberg-Marquardt method for the minimisation, where we minimise the χ 2 value: where N is the number of data points i of the SED that is considered in the fit, F mod,i the reddened flux of the model, F obs,i the observed flux of the reddened star, and O 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 (Section 3.1).The reduced-χ 2 , the χ 2 value per degree of freedom, is defined as χ 2 red ≡ χ 2 /n dof , with n dof = N − n free the degrees of freedom, where n free is the number of free parameters.
In order to obtain F mod,i , we first obtain the distancecorrected reddened model flux, f λ , by applying the extinction law under consideration to the intrinsic model spectrum: 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 F mod,i , by using the transmission curves from filters used for the construction of our observed SEDs.

An R 5495 -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 5495 , R 5495 , c 2 , c 4 , and luminosity.We find a large range of extinction values throughout the cluster, with A 5495 ranging from A 5495 = 0.76 ± 0. The χ 2 red of our best fits range from χ 2 red = 0.52 (for H90) to χ 2 red = 164.30(for H120), with an average of χ 2 red = 18.3 and a median of χ 2 red = 9.6.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 χ 2 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 5495 across the core of R136.We discuss this in more detail in Section 4.3 and Section 5.1; for now it is important to note that the wide range of values in R 5495 that we find throughout the core of R136 correspond to similar values of c 2 , for which we find values in the range of c 2 = 0.78 − 2.36, with an average of c 2 = 1.30 ± 0.22.The values show no significant trend as a function of R 5495 , contrary to what is observed for Galactic sightlines (Eq.( 3)).
We can now construct an average R 5495 -dependent extinction law towards R136.We do this by adopting the R 5495 -dependent optical and NIR law MA14 (using the spline point parameterisation of F99, as described in Section 3.3), and adopting the UV parameterisation of Fitzpatrick & Massa (1990), assuming average values of c 2 and c 4 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 Section 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 Github7 .

Extinction towards the outskirts of R136
Upon plotting the slope in the optical and NIR, captured by V −I, versus the extinction A 5495 , we see a strong correlation between  2011) sample for which A V was estimated using V − I (see Section 4.2), as well as stars in the core of the R136 cluster, for which A V 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).the two (Fig. 5): 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 5495 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 5495 = 1.88 ± 0.94.The spatial distribution of the obtained extinction values is presented in Fig. 6, and shows that the trend in A 5495 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 Section 5.1.

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 5495 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 5495 across the core of R136 as a function of spatial coordinates as follows: 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 5495 overplotted.We see a spatial trend for this quantity too, where we find higher values of R 5495 towards the east of the cluster, although this trend is less clear than in the case of A 5495 .Panel (c) of Fig. 7 shows that there is a strong correlation between R 5495 and A V , and that the relation between the two implies A 4405 − A 5495 ≈ 0.4.We note that if we refit the SED while adopting a fixed value of R 5495 (equal for all stars), we still recover the spatial trend for A 5495 , although the gradient is slightly less steep.We discuss possible causes of the spatial trend in Section 5.1.

The UV curvature (c 4 )
For the parameter c 4 , describing the curvature of the extinction curve in the far-UV (x > 5.9 µm −1 or λ < 1695 Å), we find an average value of c 4 = 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 4 = 0.42 ± 0.08 and c 4 = 0.29 ± 0.06, respectively.A low value of c 4 corresponds to a weak far-UV curvature.
We note that the value of c 4 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 4 , 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 4 .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 4 would have to approach zero, that is, no far-UV at all.

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 5495 ≈ 1.0 in the west, up to A 5495 ≈ 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 13 CO column density in the east versus the west side of the cluster.Fig. 10 shows a map of the column density of 13 CO, as well as the positions of the R136 core stars, marked by their extinction.It is clear that in regions where 13 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 13 CO column density (≥ 10 15 cm −2 ) 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 5495 .On average, we find higher values of R 5495 towards the east of the cluster than towards the west side.Higher R 5495 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 5495 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 5495 map to higher values of 13 CO column density.This could indi- 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 13 CO column density of 10 15 cm −2 (see Fig. 10).Panels c), d) and e) are as panel b), but with different background images: a NIR image captured with NIRCam on James Web 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  cate 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 5495 and R 5495 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 5495 on relatively small spatial scales (≈ 1 − 8 pc, A 5495 = 2.0 − 7.3), but higher values of R 5495 do not always coincide with higher values of A 5495 .On the contrary, combining the measurements of four different Hii regions (shown in their Figure 7), a tentative opposite trend appears, with larger values of R 5495 generally corresponding to lower A 5495 .

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 Fig. 10. 13 CO column density map of R136 and surroundings, obtained with ALMA (Wong et al. 2022).Darker regions indicate a higher 13 CO column density; in white regions no 13 CO was detected.Contours, created with DS9 (smoothing = 4), correspond to lines of constant 13 CO 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 13 CO map, we see that in regions where 13 CO is detected we measure stronger extinction.
Table 4. Revised ages and masses of the WNh stars with 1σ errors based on the stellar parameters derived by Brands et al. (2022)  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 Bonnsai9 (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).

Anomalous extinction towards R136?
The average R V value we find for R136 (R V ≈ R 5495 = 4.38 ± 0.87) is in good agreement with previously obtained averages for the region10 (Crowther et al. 2016;Bestenlehner et al. 2020, for R136, andDoran 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 previ-  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 R 5495 = 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 R V (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).Nandy et al. (1980Nandy et al. ( , 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 Fig. 11, we compare our R V -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 2 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, andalso Pei 1992, who find a very similar average curve,) use the samples of Nandy et al. (1980Nandy et al. ( , 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.

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 5495 = 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 multiwavelength 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 5495 = 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 (e τ ≈ e A V ≈ 0.4).

Fig. 1 .
Fig. 1.Integrated UV flux of the core of R136.The upper panel shows the GHRS spectrum of the inner 2 x 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 Section 3.1).The dashed lines correspond to ratios of 0.75 and 1.25.

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.

Fig. 3 .
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 R V does not hold for sightlines towards R136.In other words, the reddened flux cannot be modelled accurately over the full wavelength range: an R V 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 bestfitting 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 R V -dependent Galactic laws.

Fig. 4 .
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, R V = 3.1; black dashed line).

Fig. 5 .
Fig. 5. Extinction A 5495 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 .
Fig. 6.Extinction map of the R136 cluster and surroundings.The dots indicate stars of the De Marchi et al. (2011) sample for which A V was estimated using V − I (see Section 4.2), as well as stars in the core of the R136 cluster, for which A V 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).

Fig. 7 .
Fig. 7. Spatial trends in A V and R 5495 in the core of R136.Panels (a) and (b) show maps of A 5495 and R 5495 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 A 5495 and R 5495 and the yellow tints to lower values.A clear spatial trend is visible for A 5495 , with higher values of A 5495 towards the east of the cluster (see also Fig. 8).A similar, albeit weaker trend is visible for R 5495 .Panel (c) shows that there is a strong correlation between A 5495 and R 5495 (colours as in Fig. 5).The dotted lines in the background indicate constant values of A 4405 − A 5495 ranging from 0.2 to 0.8; the obtained relation implies A 4405 − A 5495 ≈ 0.4.

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

Fig. 9 .
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 13 CO column density of 10 15 cm −2 (see Fig. 10).Panels c), d) and e) are as panel b), but with different background images: a NIR image captured with NIRCam on James Web 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., ESO/M.-R.Cioni/VISTA Magellanic Cloud survey, Acknowledgment: Cambridge Astronomical Survey Unit.
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 13 CO column density of 10 15 cm −2 (see Fig. 10).Panels c), d) and e) are as panel b), but with different background images: a NIR image captured with NIRCam on James Web 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., ESO/M.-R.Cioni/VISTA Magellanic Cloud survey, Acknowledgment: Cambridge Astronomical Survey Unit.

Fig. 11 .
Fig. 11.Comparison of extinction curves tailored to (specific regions within) the LMC.Shown are the average extinction curve ofHowarth (1983,  green dotted line), which is nearly equal to that ofPei (1992, not shown)  and two curves ofGordon 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 R 5495 = 4.4 (the average value towards R136; red solid line).For reference, we also show the Galactic curve ofFitzpatrick (1999) for different values of R V (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 ofGordon et al. (2003).

Fig. 12 .
Fig.12.Position of stars of different samples that were used to derive extinction curves towards the LMC.The samples ofNandy et al. (1980Nandy et al. ( , 1981) ) andKoornneef & Code (1981, yellow circles)  were used to construct the curves ofHowarth (1983) andPei (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.

Table 1 .
Data used to compile SEDs of 50 stars in the core of R136.

Table 2 .
(Virtanen et al. 2020ues of optical and NIR spline anchor points based on the law of MA14.λ(Å) λ −1 (µm −1 ) A λ /(A 4405 − A 5495 )Notes.For the interpolation between these points we use a cubic spline and the function interpolate.splrep of the Python package scipy(Virtanen et al. 2020).R 5495 is the monochromatic equivalent of R V : R 5495 ≡ A 5495 /(A 4405 − A 5495 ).
and accounting for the extinction values obtained in this work.

Table B .
1. Best-fit values of the R136 core stars †.