11institutetext:  Astronomical Institute Anton Pannekoek, Amsterdam
University, Science Park 904, 1098 XH, Amsterdam, The Netherlands
11email<EMAIL_ADDRESS>22institutetext: Department of Physics and
Astronomy University of Sheffield, Sheffield, S3 7RH, United Kingdom
33institutetext: Institute of Astrophysics, KU Leuven, Celestijnenlaan 200D,
3001, Leuven, Belgium 44institutetext: Embry-Riddle Aeronautical University,
Department of Physical Science, 1 Aerospace Blvd, Daytona Beach, FL 32114
55institutetext: Argelander-Institut für Astronomie, Universität Bonn, Auf dem
Hügel 71, 53121 Bonn, Germany

# Extinction towards the cluster R136 in the Large Magellanic Cloud

An extinction law from the near-infrared to the ultraviolet

Sarah A. Brands 11 Alex de Koter 1133 Joachim M. Bestenlehner 22 Paul A.
Crowther 22 Lex Kaper 11 Saida M. Caballero-Nieves 44 Götz Gräfener 55

(Received 23 November, 2022 / Accepted 3 March 2023)

###### Abstract

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\pm 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}\sim 4.4$, about three times fewer ultraviolet photons
can escape from the ambient dust environment relative to the canonical
Galactic value of $R_{V}\sim 3.1$ at the same $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.

###### Key Words.:

Stars: early-type – Stars: massive – dust, extinction – Magellanic Clouds –
Galaxies: star clusters: individual: R136

##  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\gtrsim 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\sim 2-10\cdot 10^{4}\leavevmode\nobreak\ \mathrm{M}_{\odot}$;
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\leavevmode\nobreak\ \AA$. 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}\equiv 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}\approx 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.

##  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.

Table 1: Data used to compile SEDs of 50 stars in the core of R136.  Regime | Instrument | Grating or filter
---|---|---
UVa) | HST/STIS |  G140L ($1160-1710$ Å)
Opticalb) | HST/WFC3 |  F336W, F438W, F555W, F814Wc,d)
NIRe) | VLT/SPHERE | $J$
NIRf) | VLT/SPHERE |  $H$, $K_{s}$
a) Crowther et al. (2016). b) De Marchi et al. (2011). c) These filters
correspond roughly to the $U$, $B$, $V$ and $I$ bands, respectively. d) The
stars H35, H71, H73, H86, H116, H121, H135 lack an F814W magnitude. H69 lacks
both an F555W and F814W magnitude. H139 lacks an F336W magnitude. For all
other sources we have magnitudes for all specified bands. e) Khorrami et al.
(2017). f) Khorrami et al. (2021).

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, 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, 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 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.

Figure 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.

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, 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.

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

###  3.1 Observed SEDs

Figure 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.

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 Service111
http://svo2.cab.inta-csic.es/theory/fps/. For the UV, we bin the observed flux
points by taking the flux average over nine different wavelength regions as
specified in Table 5. 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.

###  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 ($\Delta 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_{*}$.

###  3.3 Extinction law

Figure 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 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
$R_{V}$-dependent Galactic laws.

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}\approx 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}\leq 3.31$) and
for stars elsewhere in the LMC ($R_{V}\leq 2.61$)222A note on nomenclature:
Misselt et al. (1999) use ‘30 Doradus’ to refer to a much larger region than
we do. In this study, we use the name 30 Doradus to refer to the inner
$\approx 10^{\prime}$ as shown in for example Fig. 2 of Walborn (1991) and
Fig. 1 of Evans et al. (2020), whereas Misselt et al. (1999), but also for
example Clayton & Martin (1985), consider a larger region when they refer to
30 Doradus; sometimes these authors use 30 Doradus interchangeably with LMC
SGS 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\pm 0.57,$ and for H55, we obtain $R_{V}=4.37\pm
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\pm 0.11$ for
R136a3 and $R_{V}=2.94\pm 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 ($\lambda\lesssim 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}\equiv
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 point333For one spline point,
F99 use a quadratic function; we use linear functions in all cases., 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\geq 2.7\leavevmode\nobreak\ \mu$m-1 ($\lambda\geq 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\leavevmode\nobreak\ \mu$m-1
($\lambda=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 $\leq 0.1\leavevmode\nobreak\ \%$ 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_{\lambda}\propto\lambda^{\alpha_{\mathrm{NIR}}}$ with
$\alpha_{\mathrm{NIR}}=1.61$; this value is lower than values obtained by
recent studies of Galactic extinction, which find
$\alpha_{\mathrm{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.

Table 2: $R_{5495}$-dependent values of optical and NIR spline anchor points based on the law of MA14.  $\lambda$ (Å) |  $\lambda^{-1}$ ($\mu$m-1) |  $A_{\lambda}/(A_{4405}-A_{5495}$)
---|---|---
$\infty$ | 0.000 | 0.0
26500 | 0.377 |  $-0.1097+0.1195\ \times$ $R_{5495}$
18000 | 0.556 |  $-0.2046+0.2228\ \times$ $R_{5495}$
12200 | 0.820 |  $-0.3826+0.4167\ \times$ $R_{5495}$
10000 | 1.000 |  $-0.5270+0.5740\ \times$ $R_{5495}$
8696 | 1.150 |  $-0.6392+0.7147\ \times$ $R_{5495}$
5495 | 1.820 |  $-0.0002+1.0000\ \times$ $R_{5495}$
4670 | 2.141 |  $\phantom{-}0.7455+1.0023\ \times$ $R_{5495}$
4405 | 2.270 |  $\phantom{-}1.0004+1.0000\ \times$ $R_{5495}$
4110 | 2.433 |  $\phantom{-}1.3149+0.9887\ \times$ $R_{5495}$
3704 | 2.700 |  $\phantom{-}1.7931+0.9661\ \times$ $R_{5495}$
3304 | 3.027 |  $\phantom{-}2.2580+0.9689\ \times$ $R_{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}\equiv
A_{5495}/(A_{4405}-A_{5495})$.

For the UV part of the curve, we use the parameterisation of Fitzpatrick &
Massa (1990). Expressed in terms of the quantity $k(\lambda-V)$
$\equiv(A_{\lambda}-A_{V})/(A_{B}-A_{V})\equiv E(\lambda-V)/E(B-V)$, the UV
part of the law of F99 has the following form:

$\frac{E(\lambda-V)}{E(B-V)}=\\\
\left\\{\begin{array}[]{ll}c_{1}+c_{2}x+c_{3}D(x,x_{0},\gamma)&x\leq c_{5}\\\
c_{1}+c_{2}x+c_{3}D(x,x_{0},\gamma)+c_{4}a_{1}(x-c_{5})^{2}&\\\ \hskip
71.13188pt+c_{4}a_{2}(x-c_{5})^{3}&x>c_{5},\end{array}\right.$ (1)

where $x\equiv 1/\lambda$ $\mu$m-1, $a_{1}=0.5392,$ and $a_{2}=0.05644$444As
we express the optical and NIR part of the law in monochromatic quantities,
following MA14, the $V$ and $B$ band quantities in Eq. 1 are replaced by
monochromatic values at $\lambda=5495$ Å and $\lambda=4405$ Å, respectively..
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:

$D(x,x_{0},\gamma)=\frac{x^{2}}{(x^{2}-x_{0}^{2})^{2}+x^{2}\gamma^{2}},$ (2)

with $x_{0}$ being the position and $\gamma$ 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$, $\gamma=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:

$\begin{split}c_{2}&=-0.824+4.717/R_{V}\\\
c_{1}&=\phantom{-}2.030-3.007c_{2}.\end{split}$ (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.

Figure 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).

###  3.4 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 $\theta_{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_{\rm{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 differently555Where the change as a
function of wavelength is dictated by the shape of the adopted extinction
curve., 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 lmfit666We used
lmfit version 1.0.3 in combination with Python version 3.8.5. Furthermore, we
made use of numpy version 1.19.2 (Harris et al., 2020) and scipy version 1.8.1
(Virtanen et al., 2020). (Newville et al., 2014). We adopt the least-squares
Levenberg-Marquardt method for the minimisation, where we minimise the
$\chi^{2}$ value:

$\chi^{2}=\sum_{i=0}^{N}\left(\frac{\mathcal{F}_{{\mathrm{mod}},i}-\mathcal{F}_{{\mathrm{obs}},i}}{\mathcal{O}_{i}}\right)^{2},$
(4)

where $N$ is the number of data points $i$ of the SED that is considered in
the fit, $\mathcal{F}_{{\mathrm{mod}},i}$ the reddened flux of the model,
$\mathcal{F}_{{\mathrm{obs}},i}$ the observed flux of the reddened star, and
$\mathcal{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-$\chi^{2}$, the $\chi^{2}$ value
per degree of freedom, is defined as
$\chi^{2}_{\mathrm{red}}\equiv\chi^{2}/n_{\mathrm{dof}}$, with
$n_{\mathrm{dof}}=N-n_{\mathrm{free}}$ the degrees of freedom, where
$n_{\mathrm{free}}$ is the number of free parameters.

In order to obtain $\mathcal{F}_{\mathrm{mod},i}$, we first obtain the
distance-corrected reddened model flux, $f_{\lambda}$, by applying the
extinction law under consideration to the intrinsic model spectrum:

$f_{\lambda}=\theta^{2}_{R}F_{\lambda}10^{-0.4A_{\lambda}},$ (5)

with $F_{\lambda}$ the intrinsic (unreddened), distance-corrected model flux,
and $A_{\lambda}$ the extinction in magnitudes as a function of wavelength,
dictated by the adopted law and $A_{V}$. As discussed above, $\theta_{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
$\mathcal{F}_{\mathrm{mod},i}$, by using the transmission curves from filters
used for the construction of our observed SEDs.

##  4 Results

###  4.1 An $R_{5495}$-dependent average extinction law for R136

Figure 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\sigma$ 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.  Figure 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).

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\pm 0.06$ (H108) to
$A_{5495}=2.63\pm 0.05$ (H36), and $R_{5495}$ ranging from $R_{5495}=1.91\pm
0.26$ (H108) to $R_{5495}=5.79\pm 0.31$ (H36), and cluster averages of
$A_{5495}=1.71\pm 0.41$ and of $R_{V}=4.38\pm 0.87$. The best-fit parameters
of individual stars are listed in Table 6. Two example fits are shown in Fig.
4.

The $\chi^{2}_{\mathrm{red}}$ of our best fits range from
$\chi^{2}_{\mathrm{red}}=0.52$ (for H90) to $\chi^{2}_{\mathrm{red}}=164.30$
(for H120), with an average of $\chi^{2}_{\mathrm{red}}=18.3$ and a median of
$\chi^{2}_{\mathrm{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 $\chi^{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\pm 0.22$. The values show no significant trend as a function of
$R_{5495}$, contrary to what is observed for Galactic sightlines (Eq. 3).

Figure 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}\approx 0.4$.

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
Github777 https://github.com/sarahbrands/ExtinctionR136/.

Table 3: Average $R_{5495}$-dependent extinction law for R136. Parameter |  Value(s) |  Source
---|---|---
Spline points |  See Table 2 |  MA14, F99, this work
$c_{2}$ |  $1.30$ |  This work
$c_{1}$ |  2.030 – 3.007 $c_{2}$ |  F99
$c_{3}$ |  1.463 |  Gordon et al. (2003)
$c_{4}$ |  $0.09$ |  This work
$c_{5}$ |  5.9 |  Fitzpatrick & Massa (1990)
$x_{0}$ |  4.558 |  Gordon et al. (2003)
$\gamma$ |  0.945 |  Gordon et al. (2003)
Notes. Only the optical/NIR part of the law has a dependence on $R_{5495}$.
The optical/NIR part of the law is as derived by MA14, but we adopt a
functional form similar to that of F99; see Section 3.3. A Python
implementation of the full law can be found in Appendix C and on
Github${}^{\ref{githubfootnote}}$.

###  4.2 Extinction towards the outskirts of R136

Figure 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.

Upon plotting the slope in the optical and NIR, captured by $V-I$, versus the
extinction $A_{5495}$, we see a strong correlation between the two (Fig. 5):

$A_{5495}=[0.38\pm 0.23]+[3.86\pm 0.58](V-I).$ (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_{5495}$ for 1657 sources with
$V\leavevmode\nobreak\ <\leavevmode\nobreak\ 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\pm 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.

###  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_{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:

$A_{5495}=[1.93\pm 0.04]+[1.74\pm 0.23](RA-84.68)/0.006,$ (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_{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}\approx 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.

Figure 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 ($\approx 1-2$ pc). The
green-yellow dashed line indicates the contour of a 13CO 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. Figure 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 constant 13CO column density, the value of which is indicated in the plot
in units of $\log[\mathrm{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 where 13CO is detected we measure stronger extinction.

###  4.4 The UV curvature ($c_{4}$)

For the parameter $c_{4}$, describing the curvature of the extinction curve in
the far-UV ($x>5.9\leavevmode\nobreak\ \mu\mathrm{m}^{-1}$ or
$\lambda<1695\leavevmode\nobreak\ \mathrm{\AA}$), we find an average value of
$c_{4}=0.09\pm 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\pm 0.08$ and $c_{4}=0.29\pm 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 $\lambda<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
$\lambda=1280$ Å), the value of $c_{4}$ 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_{5495}\approx 1.0$ in the west, up to
$A_{5495}\approx 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
$\approx 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$\times$ the
beam size, allowing us to distinguish between the 13CO column density in the
east versus the west side of the cluster. Fig. 10 shows a map of the column
density of 13CO, as well as the positions of the R136 core stars, marked by
their extinction. It is clear that in regions where 13CO 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 13CO column
density ($\geq 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\leavevmode\nobreak\ \AA$) 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 13CO 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_{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 ($\approx 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}$.

###  5.2 Revised WNh star masses

Table 4: Revised ages and masses of the WNh stars with 1$\sigma$ errors based on the stellar parameters derived by Brands et al. (2022) and accounting for the extinction values obtained in this work.  |  R136a1 |  R136a2 |  R136a3
---|---|---|---
${M}_{\rm act}$ (M⊙) | $237\pm^{28}_{22}$ | $140\pm^{10}_{9}$ | $174\pm^{18}_{10}$
${M}_{\rm ini}$ (M⊙) | $277\pm^{24}_{27}$ | $171\pm^{11}_{7}$ | $209\pm^{10}_{9}$
Age (Myr) | $0.94\pm^{0.17}_{0.16}$ | $1.64\pm^{0.13}_{0.12}$ | $1.32\pm^{0.19}_{0.25}$
log $L/L_{\odot}$ | $6.82\pm^{0.06}_{0.04}$ | $6.59\pm^{0.02}_{0.03}$ | $6.69\pm^{0.02}_{0.02}$
Figure 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 $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).  Figure 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.

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_{\odot}$ 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)888The HRD for the core of R136 with tracks
of Brott et al. (2011) and Köhler et al. (2015) can be found in Fig. 10 of
Brands et al. (2022).. 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 Bonnsai999The Bonnsai web-service is available at
https://www.astro.uni-bonn.de/stars/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).

###  5.3 Anomalous extinction towards R136?

The average $R_{V}$ value we find for R136 ($R_{V}\approx R_{5495}=4.38\pm
0.87$) is in good agreement with previously obtained averages for the
region101010In this section, we use the broadband and the monochromatic
definitions of the total-to-relative extinction, $R_{V}$ and $R_{5495}$,
interchangeably: as the sources that we discuss are not heavily reddened, and
their SEDs are fairly similar, the differences between them are assumed to be
sufficiently small. (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 $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, and also , 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.

##  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}\approx A_{5495}=1.70\pm 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}\approx R_{5495}=4.38\pm 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 $\lambda\approx
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^{\tau}\approx
e^{A_{V}}\approx 0.4$).

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). Observations were taken with the NASA/ESA HST,
obtained from the data archive at the Space Telescope Institute. This research
has made use of the SIMBAD database, operated at CDS, Strasbourg, France
(Wenger et al., 2000).

## References

  * Andersen et al. (2009) Andersen, M., Zinnecker, H., Moneti, A., et al. 2009, ApJ, 707, 1347
  * Banerjee et al. (2012) Banerjee, S., Kroupa, P., & Oh, S. 2012, MNRAS, 426, 1416
  * Bestenlehner et al. (2020) Bestenlehner, J. M., Crowther, P. A., Caballero-Nieves, S. M., et al. 2020, MNRAS, 499, 1918
  * Bestenlehner et al. (2014) Bestenlehner, J. M., Gräfener, G., Vink, J. S., et al. 2014, A&A, 570, A38
  * Bonnarel et al. (2000) Bonnarel, F., Fernique, P., Bienaymé, O., et al. 2000, A&AS, 143, 33
  * Brandner (2002) Brandner, W. 2002, in Astronomical Society of the Pacific Conference Series, Vol. 285, Modes of Star Formation and the Origin of Field Populations, ed. E. K. Grebel & W. Brandner, 105
  * Brands et al. (2022) Brands, S. A., de Koter, A., Bestenlehner, J. M., et al. 2022, A&A, 663, A36
  * Brott et al. (2011) Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682
  * Calzetti et al. (1994) Calzetti, D., Kinney, A. L., & Storchi-Bergmann, T. 1994, ApJ, 429, 582
  * Campbell et al. (2010) Campbell, M. A., Evans, C. J., Mackey, A. D., et al. 2010, MNRAS, 405, 421
  * Cardamone et al. (2009) Cardamone, C., Schawinski, K., Sarzi, M., et al. 2009, MNRAS, 399, 1191
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
  * Carneiro et al. (2016) Carneiro, L. P., Puls, J., Sundqvist, J. O., & Hoffmann, T. L. 2016, A&A, 590, A88
  * Cheng et al. (2021) Cheng, Y., Wang, Q. D., & Lim, S. 2021, MNRAS, 504, 1627
  * Clayton & Martin (1985) Clayton, G. C. & Martin, P. G. 1985, ApJ, 288, 558
  * Crowther (2019) Crowther, P. A. 2019, Galaxies, 7, 88
  * Crowther et al. (2016) Crowther, P. A., Caballero-Nieves, S. M., Bostroem, K. A., et al. 2016, MNRAS, 458, 624
  * Crowther et al. (2017) Crowther, P. A., Castro, N., Evans, C. J., et al. 2017, The Messenger, 170, 40
  * Crowther et al. (2010) Crowther, P. A., Schnurr, O., Hirschi, R., et al. 2010, MNRAS, 408, 731
  * de Koter et al. (1997) de Koter, A., Heap, S. R., & Hubeny, I. 1997, ApJ, 477, 792
  * De Marchi & Panagia (2014) De Marchi, G. & Panagia, N. 2014, MNRAS, 445, 93
  * De Marchi & Panagia (2019) De Marchi, G. & Panagia, N. 2019, ApJ, 878, 31
  * De Marchi et al. (2016) De Marchi, G., Panagia, N., Sabbi, E., et al. 2016, MNRAS, 455, 4373
  * De Marchi et al. (2011) De Marchi, G., Paresce, F., Panagia, N., et al. 2011, ApJ, 739, 27
  * Doran et al. (2013) Doran, E. I., Crowther, P. A., de Koter, A., et al. 2013, A&A, 558, A134
  * Draine (2003) Draine, B. T. 2003, ARA&A, 41, 241
  * Evans et al. (2020) Evans, C., Lennon, D., Langer, N., et al. 2020, The Messenger, 181, 22
  * Fitzpatrick (1999) Fitzpatrick, E. L. 1999, PASP, 111, 63
  * Fitzpatrick & Massa (1990) Fitzpatrick, E. L. & Massa, D. 1990, ApJS, 72, 163
  * Fitzpatrick & Massa (2005) Fitzpatrick, E. L. & Massa, D. 2005, AJ, 130, 1127
  * Fitzpatrick et al. (2019) Fitzpatrick, E. L., Massa, D., Gordon, K. D., Bohlin, R., & Clayton, G. C. 2019, ApJ, 886, 108
  * Galliano et al. (2018) Galliano, F., Galametz, M., & Jones, A. P. 2018, ARA&A, 56, 673
  * Garcia et al. (2021) Garcia, M., Evans, C. J., Bestenlehner, J. M., et al. 2021, Experimental Astronomy, 51, 887
  * Gordon et al. (2003) Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, ApJ, 594, 279
  * Gräfener (2021) Gräfener, G. 2021, A&A, 647, A13
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Heap et al. (1992) Heap, S. R., Ebbets, D., & Malumuth, E. 1992, in European Southern Observatory Conference and Workshop Proceedings, Vol. 44, European Southern Observatory Conference and Workshop Proceedings, 347
  * Higgins et al. (2022) Higgins, E. R., Vink, J. S., Sabhahit, G. N., & Sander, A. A. C. 2022, MNRAS, 516, 4052
  * Howarth (1983) Howarth, I. D. 1983, MNRAS, 203, 301
  * Hunter et al. (1995) Hunter, D. A., Shaya, E. J., Holtzman, J. A., et al. 1995, ApJ, 448, 179
  * Kalari et al. (2022) Kalari, V. M., Horch, E. P., Salinas, R., et al. 2022, ApJ, 935, 162
  * Kalari et al. (2018) Kalari, V. M., Rubio, M., Elmegreen, B. G., et al. 2018, ApJ, 852, 71
  * Kennicutt (1984) Kennicutt, R. C., J. 1984, ApJ, 287, 116
  * Khorrami et al. (2021) Khorrami, Z., Langlois, M., Clark, P. C., et al. 2021, MNRAS, 503, 292
  * Khorrami et al. (2017) Khorrami, Z., Vakili, F., Lanz, T., et al. 2017, A&A, 602, A56
  * Köhler et al. (2015) Köhler, K., Langer, N., de Koter, A., et al. 2015, A&A, 573, A71
  * Koornneef & Code (1981) Koornneef, J. & Code, A. D. 1981, ApJ, 247, 860
  * Lanz & Hubeny (2003) Lanz, T. & Hubeny, I. 2003, ApJS, 146, 417
  * Lee et al. (2019) Lee, M. Y., Madden, S. C., Le Petit, F., et al. 2019, A&A, 628, A113
  * Maiz-Apellaniz (2005) Maiz-Apellaniz, J. 2005, Instrument Science Report STIS 2005-02 (Baltimore:STSc, Space Telescope STIS Instrument Science Report
  * Maiz-Apellaniz (2007) Maiz-Apellaniz, J. 2007, MULTISPEC: A Code for the Extraction of Slitless Spectra in Crowded Fields
  * Maíz Apellániz (2013) Maíz Apellániz, J. 2013, in Highlights of Spanish Astrophysics VII, ed. J. C. Guirado, L. M. Lara, V. Quilis, & J. Gorgas, 583–589
  * Maíz Apellániz & Barbá (2018) Maíz Apellániz, J. & Barbá, R. H. 2018, A&A, 613, A9
  * Maíz Apellániz et al. (2014) Maíz Apellániz, J., Evans, C. J., Barbá, R. H., et al. 2014, A&A, 564, A63
  * Maíz Apellániz et al. (2020) Maíz Apellániz, J., Pantaleoni González, M., Barbá, R. H., García-Lario, P., & Nogueras-Lara, F. 2020, MNRAS, 496, 4951
  * Misselt et al. (1999) Misselt, K. A., Clayton, G. C., & Gordon, K. D. 1999, ApJ, 515, 128
  * Mokiem et al. (2007) Mokiem, M. R., de Koter, A., Vink, J. S., et al. 2007, A&A, 473, 603
  * Nandy et al. (1981) Nandy, K., Morgan, D. H., Willis, A. J., Wilson, R., & Gondhalekar, P. M. 1981, MNRAS, 196, 955
  * Nandy et al. (1980) Nandy, K., Morgan, D. H., Willis, A. J., et al. 1980, Nature, 283, 725
  * Newville et al. (2014) Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python, Zenodo
  * Nogueras-Lara et al. (2019) Nogueras-Lara, F., Schödel, R., Najarro, F., et al. 2019, A&A, 630, L3
  * O’Connell (2010) O’Connell, R. W. 2010, in American Astronomical Society Meeting Abstracts, Vol. 215, American Astronomical Society Meeting Abstracts #215, 222.05
  * Pei (1992) Pei, Y. C. 1992, ApJ, 395, 130
  * Pellegrini et al. (2011) Pellegrini, E. W., Baldwin, J. A., & Ferland, G. J. 2011, ApJ, 738, 34
  * Pietrzyński et al. (2019) Pietrzyński, G., Graczyk, D., Gallenne, A., et al. 2019, Nature, 567, 200
  * Puls et al. (2005) Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669
  * Ramírez-Agudelo et al. (2013) Ramírez-Agudelo, O. H., Simón-Díaz, S., Sana, H., et al. 2013, A&A, 560, A29
  * Rivero González et al. (2012) Rivero González, J. G., Puls, J., Najarro, F., & Brott, I. 2012, A&A, 537, A79
  * Rubio-Díez et al. (2017) Rubio-Díez, M. M., Najarro, F., García, M., & Sundqvist, J. O. 2017, in The Lives and Death-Throes of Massive Stars, ed. J. J. Eldridge, J. C. Bray, L. A. S. McClelland, & L. Xiao, Vol. 329, 131–135
  * Salim & Narayanan (2020) Salim, S. & Narayanan, D. 2020, ARA&A, 58, 529
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161
  * Santolaya-Rey et al. (1997) Santolaya-Rey, A. E., Puls, J., & Herrero, A. 1997, A&A, 323, 488
  * Schneider et al. (2017) Schneider, F. R. N., Castro, N., Fossati, L., Langer, N., & de Koter, A. 2017, A&A, 598, A60
  * Stead & Hoare (2009) Stead, J. J. & Hoare, M. G. 2009, MNRAS, 400, 731
  * Sundqvist & Puls (2018) Sundqvist, J. O. & Puls, J. 2018, A&A, 619, A59
  * Thornley et al. (1999) Thornley, M. D., Schreiber, N. M. F., Spoon, H. W. W., et al. 1999, in New Views of the Magellanic Clouds, ed. Y. H. Chu, N. Suntzeff, J. Hesser, & D. Bohlender, Vol. 190, 247
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Walborn (1991) Walborn, N. R. 1991, in The Magellanic Clouds, ed. R. Haynes & D. Milne, Vol. 148, 145
  * Weingartner & Draine (2001) Weingartner, J. C. & Draine, B. T. 2001, ApJ, 548, 296
  * Wenger et al. (2000) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9
  * Whiteoak (1966) Whiteoak, J. B. 1966, ApJ, 144, 305
  * Wong et al. (2022) Wong, T., Oudshoorn, L., Sofovich, E., et al. 2022, ApJ, 932, 47
  * Xiang et al. (2017) Xiang, F. Y., Li, A., & Zhong, J. X. 2017, ApJ, 835, 107

##  Appendix A UV bands

Table 5 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.

Table 5: Minimum and maximum wavelengths ($\lambda_{\rm min}$ and $\lambda_{\rm max}$) of synthetic passbands used for the UV, named by their (rounded) central wavelength. Passband |  $\lambda_{\rm min}$ (Å) |  $\lambda_{\rm max}$ (Å)
---|---|---
$\lambda$1280 | 1265.00 | 1295.00
$\lambda$1319 | 1307.00 | 1330.00
$\lambda$1351 | 1347.00 | 1355.00
$\lambda$1420 | 1405.00 | 1435.00
$\lambda$1455 | 1440.00 | 1470.00
$\lambda$1498 | 1475.00 | 1520.00
$\lambda$1585 | 1570.00 | 1600.00
$\lambda$1618 | 1612.00 | 1624.00
$\lambda$1688 | 1675.00 | 1700.00

##  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 6.

Table 6: Best-fit values of the R136 core stars$\dagger$. Source | $R_{V}$ | $A_{V}$ | $A_{K}$ | $\log L/L_{\odot}$ | $c_{2}$ | $c_{4}$ | $\chi^{2}_{\mathrm{red}}$
---|---|---|---|---|---|---|---
R136a1 | 4.25$\pm$ | 0.10 | 1.85$\pm$ | 0.02 | 0.24 | $6.82\pm 0.01$ | 1.28$\pm$ | 0.05 | 0.00$\pm$ | ${}^{0.12}_{0.00}$ | 4.44
R136a2 | 2.93$\pm$ | 0.15 | 1.25$\pm$ | 0.04 | 0.14 | $6.59\pm 0.01$ | 1.62$\pm$ | 0.11 | 0.00$\pm$ | ${}^{0.34}_{0.00}$ | 14.32
R136a3 | 4.63$\pm$ | 0.08 | 1.99$\pm$ | 0.02 | 0.26 | $6.69\pm 0.01$ | 1.14$\pm$ | 0.03 | 0.07$\pm$ | ${}^{0.12}_{0.07}$ | 1.81
R136a4 | 3.30$\pm$ | 0.14 | 1.38$\pm$ | 0.03 | 0.17 | $6.23\pm 0.01$ | 1.43$\pm$ | 0.07 | 0.00$\pm$ | ${}^{0.42}_{0.00}$ | 4.83
R136a5 | 4.67$\pm$ | 0.31 | 2.00$\pm$ | 0.05 | 0.27 | $6.34\pm 0.01$ | 1.33$\pm$ | 0.13 | 0.00$\pm$ | ${}^{0.02}_{0.00}$ | 13.18
R136a7 | 5.04$\pm$ | 0.61 | 1.82$\pm$ | 0.07 | 0.25 | $6.29\pm 0.02$ | 1.87$\pm$ | 0.30 | 0.00$\pm$ | ${}^{0.69}_{0.00}$ | 24.60
R136b | 5.66$\pm$ | 0.21 | 2.57$\pm$ | 0.04 | 0.36 | $6.39\pm 0.01$ | 1.18$\pm$ | 0.08 | 0.06$\pm$ | ${}^{0.20}_{0.06}$ | 8.72
H30 | 4.13$\pm$ | 0.13 | 1.63$\pm$ | 0.03 | 0.21 | $5.72\pm 0.01$ | 1.15$\pm$ | 0.06 | 0.06$\pm$ | ${}^{0.15}_{0.06}$ | 3.90
H31 | 3.94$\pm$ | 0.24 | 1.53$\pm$ | 0.05 | 0.20 | $5.93\pm 0.02$ | 1.24$\pm$ | 0.13 | 0.06$\pm$ | ${}^{0.39}_{0.06}$ | 19.94
H35 | 4.66$\pm$ | 0.29 | 1.57$\pm$ | 0.04 | 0.21 | $5.81\pm 0.01$ | 1.34$\pm$ | 0.09 | 0.05$\pm$ | ${}^{0.20}_{0.05}$ | 5.69
H36 | 5.79$\pm$ | 0.31 | 2.63$\pm$ | 0.05 | 0.37 | $6.31\pm 0.02$ | 1.15$\pm$ | 0.12 | 0.03$\pm$ | ${}^{0.35}_{0.03}$ | 17.16
H40 | 5.55$\pm$ | 0.41 | 1.95$\pm$ | 0.06 | 0.27 | $5.93\pm 0.02$ | 1.30$\pm$ | 0.15 | 0.11$\pm$ | ${}^{0.41}_{0.11}$ | 20.02
H45 | 5.38$\pm$ | 0.13 | 2.20$\pm$ | 0.02 | 0.30 | $5.82\pm 0.01$ | 1.25$\pm$ | 0.06 | 0.13$\pm$ | ${}^{0.19}_{0.13}$ | 3.39
H46 | 5.04$\pm$ | 0.32 | 2.38$\pm$ | 0.06 | 0.32 | $6.12\pm 0.02$ | 1.18$\pm$ | 0.19 | 0.13$\pm$ | ${}^{0.64}_{0.13}$ | 23.42
H47 | 5.47$\pm$ | 0.52 | 2.50$\pm$ | 0.09 | 0.35 | $5.99\pm 0.03$ | 1.15$\pm$ | 0.31 | 0.22$\pm$ | ${}^{1.01}_{0.22}$ | 39.86
H48 | 5.42$\pm$ | 0.37 | 2.31$\pm$ | 0.06 | 0.32 | $5.98\pm 0.01$ | 1.21$\pm$ | 0.13 | 0.06$\pm$ | ${}^{0.35}_{0.06}$ | 13.09
H50 | 4.17$\pm$ | 0.10 | 1.83$\pm$ | 0.02 | 0.24 | $5.84\pm 0.00$ | 1.11$\pm$ | 0.04 | 0.09$\pm$ | ${}^{0.11}_{0.09}$ | 1.59
H52 | 4.80$\pm$ | 0.23 | 1.67$\pm$ | 0.03 | 0.23 | $5.69\pm 0.01$ | 1.29$\pm$ | 0.10 | 0.14$\pm$ | ${}^{0.24}_{0.14}$ | 4.87
H55 | 4.66$\pm$ | 0.38 | 1.73$\pm$ | 0.06 | 0.23 | $5.76\pm 0.02$ | 1.27$\pm$ | 0.16 | 0.07$\pm$ | ${}^{0.44}_{0.07}$ | 18.26
H58 | 4.96$\pm$ | 0.41 | 2.06$\pm$ | 0.06 | 0.28 | $5.87\pm 0.01$ | 1.31$\pm$ | 0.18 | 0.15$\pm$ | ${}^{0.48}_{0.15}$ | 13.36
H62 | 4.54$\pm$ | 0.17 | 1.68$\pm$ | 0.02 | 0.22 | $5.61\pm 0.01$ | 1.18$\pm$ | 0.06 | 0.21$\pm$ | 0.14 | 1.86
H64 | 5.23$\pm$ | 0.45 | 2.35$\pm$ | 0.08 | 0.32 | $5.85\pm 0.02$ | 1.29$\pm$ | 0.21 | 0.04$\pm$ | ${}^{0.60}_{0.04}$ | 35.26
H66 | 4.06$\pm$ | 0.17 | 1.78$\pm$ | 0.03 | 0.23 | $5.65\pm 0.01$ | 1.12$\pm$ | 0.07 | 0.25$\pm$ | 0.19 | 2.27
H68 | 5.21$\pm$ | 0.13 | 2.34$\pm$ | 0.02 | 0.32 | $5.69\pm 0.01$ | 1.20$\pm$ | 0.07 | 0.03$\pm$ | ${}^{0.22}_{0.03}$ | 4.19
H69 | 4.92$\pm$ | 0.59 | 1.75$\pm$ | 0.07 | 0.24 | $5.45\pm 0.01$ | 1.46$\pm$ | 0.14 | 0.08$\pm$ | ${}^{0.14}_{0.08}$ | 1.19
H70 | 4.18$\pm$ | 0.45 | 1.85$\pm$ | 0.08 | 0.24 | $5.68\pm 0.02$ | 1.34$\pm$ | 0.28 | 0.19$\pm$ | ${}^{0.88}_{0.19}$ | 17.69
H71 | 4.45$\pm$ | 0.82 | 1.55$\pm$ | 0.11 | 0.20 | $5.44\pm 0.03$ | 1.35$\pm$ | 0.33 | 0.07$\pm$ | ${}^{0.76}_{0.07}$ | 43.68
H73 | 3.68$\pm$ | 0.23 | 1.65$\pm$ | 0.04 | 0.21 | $5.15\pm 0.01$ | 0.96$\pm$ | 0.07 | 0.16$\pm$ | 0.16 | 3.09
H75 | 3.72$\pm$ | 0.12 | 1.41$\pm$ | 0.02 | 0.18 | $5.45\pm 0.00$ | 1.40$\pm$ | 0.06 | 0.17$\pm$ | 0.14 | 1.34
H78 | 4.59$\pm$ | 0.24 | 1.71$\pm$ | 0.03 | 0.23 | $5.50\pm 0.01$ | 1.45$\pm$ | 0.09 | 0.07$\pm$ | ${}^{0.20}_{0.07}$ | 2.17
H80 | 3.91$\pm$ | 0.21 | 1.40$\pm$ | 0.03 | 0.18 | $5.09\pm 0.01$ | 1.19$\pm$ | 0.08 | 0.00$\pm$ | ${}^{0.19}_{0.00}$ | 2.30
H86 | 2.36$\pm$ | 0.60 | 1.06$\pm$ | 0.14 | 0.11 | $5.42\pm 0.03$ | 1.22$\pm$ | 0.35 | 0.00$\pm$ | ${}^{0.02}_{0.00}$ | 56.89
H90 | 4.85$\pm$ | 0.15 | 1.79$\pm$ | 0.02 | 0.24 | $5.34\pm 0.00$ | 1.31$\pm$ | 0.05 | 0.03$\pm$ | ${}^{0.09}_{0.03}$ | 0.52
H92 | 4.66$\pm$ | 0.16 | 1.64$\pm$ | 0.02 | 0.22 | $5.21\pm 0.00$ | 1.27$\pm$ | 0.05 | 0.01$\pm$ | ${}^{0.10}_{0.01}$ | 0.63
H94 | 5.15$\pm$ | 0.26 | 1.87$\pm$ | 0.03 | 0.26 | $5.39\pm 0.01$ | 1.39$\pm$ | 0.08 | 0.04$\pm$ | ${}^{0.17}_{0.04}$ | 1.44
H108 | 1.91$\pm$ | 0.26 | 0.77$\pm$ | 0.06 | 0.07 | $4.88\pm 0.01$ | 1.21$\pm$ | 0.15 | 0.08$\pm$ | ${}^{0.33}_{0.08}$ | 9.05
H112 | 4.19$\pm$ | 0.98 | 1.85$\pm$ | 0.15 | 0.24 | $5.08\pm 0.03$ | 1.16$\pm$ | 0.35 | 0.10$\pm$ | ${}^{0.72}_{0.10}$ | 24.59
H114 | 4.64$\pm$ | 0.69 | 1.55$\pm$ | 0.08 | 0.21 | $5.21\pm 0.02$ | 1.40$\pm$ | 0.24 | 0.00$\pm$ | ${}^{0.70}_{0.00}$ | 11.09
H116 | 4.85$\pm$ | 1.57 | 1.45$\pm$ | 0.15 | 0.20 | $4.86\pm 0.03$ | 1.54$\pm$ | 0.54 | 0.07$\pm$ | ${}^{1.01}_{0.07}$ | 37.99
H120 | 2.36$\pm$ | 0.88 | 0.93$\pm$ | 0.20 | 0.09 | $4.78\pm 0.06$ | 1.24$\pm$ | 0.57 | 0.00$\pm$ | ${}^{1.19}_{0.00}$ | 164.30
H121 | 4.33$\pm$ | 1.39 | 1.46$\pm$ | 0.17 | 0.19 | $4.77\pm 0.03$ | 1.24$\pm$ | 0.48 | 0.23$\pm$ | ${}^{1.04}_{0.23}$ | 48.33
H123 | 3.39$\pm$ | 0.87 | 1.22$\pm$ | 0.14 | 0.15 | $4.87\pm 0.03$ | 1.37$\pm$ | 0.44 | 0.03$\pm$ | ${}^{1.06}_{0.03}$ | 59.44
H132 | 3.79$\pm$ | 0.51 | 1.49$\pm$ | 0.08 | 0.19 | $5.01\pm 0.01$ | 1.22$\pm$ | 0.19 | 0.10$\pm$ | ${}^{0.37}_{0.10}$ | 7.05
H134 | 4.16$\pm$ | 1.38 | 1.27$\pm$ | 0.16 | 0.16 | $4.72\pm 0.04$ | 1.46$\pm$ | 0.62 | 0.09$\pm$ | ${}^{1.43}_{0.09}$ | 62.54
H135 | 2.98$\pm$ | 0.69 | 1.36$\pm$ | 0.13 | 0.16 | $4.65\pm 0.02$ | 0.78$\pm$ | 0.20 | 0.32$\pm$ | ${}^{0.33}_{0.32}$ | 10.99
H139 | 4.67$\pm$ | 1.13 | 1.47$\pm$ | 0.13 | 0.20 | $4.90\pm 0.02$ | 1.43$\pm$ | 0.38 | 0.00$\pm$ | ${}^{1.08}_{0.00}$ | 21.19
H141 | 5.03$\pm$ | 0.53 | 1.43$\pm$ | 0.05 | 0.19 | $4.67\pm 0.01$ | 2.36$\pm$ | 0.30 | 0.14$\pm$ | ${}^{0.69}_{0.14}$ | 2.21
H143 | 3.75$\pm$ | 0.62 | 1.75$\pm$ | 0.11 | 0.22 | $5.11\pm 0.02$ | 1.06$\pm$ | 0.21 | 0.22$\pm$ | ${}^{0.43}_{0.22}$ | 10.19
H159 | 5.36$\pm$ | 0.92 | 1.72$\pm$ | 0.09 | 0.24 | $4.69\pm 0.01$ | 1.47$\pm$ | 0.32 | 0.29$\pm$ | ${}^{0.70}_{0.29}$ | 7.74
H173 | 3.54$\pm$ | 0.54 | 1.18$\pm$ | 0.07 | 0.15 | $4.42\pm 0.01$ | 1.31$\pm$ | 0.22 | 0.19$\pm$ | ${}^{0.42}_{0.19}$ | 4.91
Average | 4.38$\pm$ | 0.87 | 1.71$\pm$ | 0.41 | 0.22 |  |  $1.30$$\pm$ | $0.22$ |  $0.09$$\pm$ | $0.08$ |
$\dagger$ The parameters $c_{2}$ and $c_{4}$ correspond to the slope of the UV
extinction, and the far-UV curvature, respectively. For all stars, we adopted
$c_{1}=2.030-3.007c_{2}$; $c_{3}=1.463$; $c_{5}=5.9$, $x_{0}=4.558$, and
$\gamma=0.945$. The uncertainties quoted here are statistical uncertainties
resulting from the extinction fits; in the case of luminosity, the
uncertainties associated with the derivation of stellar parameters (see Brands
et al. (2022), their Table A.1.) are not included.

##  Appendix C The extinction law towards R136

The $R_{5495}$-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_{5495}$-dependent extinction law towards R136 is presented in
Listing 1 and can also be found on
Github111111https://github.com/sarahbrands/ExtinctionR136/. 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
$R_{5495}$-dependent extinction law towards R136.

⬇

import numpy as np

import scipy.interpolate as interpolate

def exctinction_R136(waves,Rv):

"""

Extinction law tailored to the cluster R136 in the Large Magellanic Cloud,

where a strong gradient in R5495 values is observed.

- Shape of the optical and NIR part as in Maiz-Apellaniz et al. (2014),

but parameterised in a way similar to that of Fitzpatrick (1999).

- UV part of the curve tailored to R136.

[Input]

- waves [numpy array]: wavelength range Angstrom

- Rv [float]: monochromatic total-to-relative extinction R_5495:

R5495 = A(lam=5495)/(A(lam=4405)-A(lam=5495))

[Output]:

- Alam_Av [numpy array]: normalised extinction as a function of

wavelength, that is, A(lambda)/A_5495.

Note that (in broadband equivalents): A(lambda)/A(V) = curve / Rv, with

curve = k(lambda-V) + Rv = E(lambda-V)/E(B-V) + Rv = Alam / E(B-V)

This expression is used in e.g. Fitzpatrick et al. (2007).

"""

# Parameters of the UV part of the extinction curve

c2 = 1.30 # R136 average

c1 = 2.030 - 3.007*c2 # As in Fitzpatrick (1999)

c3 = 1.463 # Bump parameters from Gordon et al. (2003).

c4 = 0.09 # R136 average

c5 = 5.9 # As in Fitzpatrick (1999)

x0 = 4.558 # Bump parameters from Gordon et al. (2003).

gamma = 0.945 # Bump parameters from Gordon et al. (2003).

# Defining inverse wavelength range and the curve components

xx = 10000./ np.array(waves)

curve = xx*0.

xcutuv = 10000.0/2700.0

xspluv = 10000.0/np.array([2700.0,2600.0])

iuv = np.where(xx >= xcutuv)[0]

iopir = np.where(xx < xcutuv)[0]

if (len(iuv) > 0):

xuv = np.concatenate((xspluv,xx[iuv]))

else:

xuv = xspluv

# UV part of the curve

yuv = c1 + c2*xuv

yuv = yuv + c3*xuv**2/((xuv**2-x0**2)**2 +(xuv*gamma)**2)

yuv = yuv + c4*(0.5392*(np.maximum(xuv,c5)-c5)**2 +

0.05644*(np.maximum(xuv,c5)-c5)**3) + Rv

yspluv = yuv[0:2]

if (len(iuv) > 0):

curve[iuv] = yuv[2::]

# Optical and NIR spline points

xsplopir = np.concatenate(([0],10000.0/np.array([26500.0, 18000.0, 12200.0,

10000.0, 8696.0, 5495.0, 4670.0, 4405.0, 4110.0, 3704.0, 3304.0])))

ysplopir = np.array((np.polyval([-0.1097 ,0.1195][::-1], Rv),

np.polyval([-0.2046 ,0.2228][::-1], Rv),

np.polyval([-0.3826 ,0.4167][::-1], Rv),

np.polyval([-0.5270 ,0.5740][::-1], Rv),

np.polyval([-0.6392 ,0.7147][::-1], Rv),

np.polyval([-0.0002 ,1.0000][::-1], Rv),

np.polyval([0.7455 ,1.0023][::-1], Rv),

np.polyval([1.0004 ,1.0000][::-1], Rv),

np.polyval([1.3149 ,0.9887][::-1], Rv),

np.polyval([1.7931 ,0.9661][::-1], Rv),

np.polyval([2.2580 ,0.9689][::-1], Rv)))

ysplopir = np.concatenate((np.array([0]),ysplopir))

# If the wavelength range covers optical/NIR, interpolate from UV to optical

if (len(iopir) > 0):

tck=interpolate.splrep(np.concatenate((xsplopir,xspluv)),

np.concatenate((ysplopir,yspluv)),k=3)

curve[iopir] = interpolate.splev(xx[iopir], tck)

Alam_Av = curve/Rv

return Alam_Av