<figure><img src="image_1.png" /></figure>

# Evidence of Dust Grain Evolution from Extinction Mapping in the IC 63 Photodissociation Region*

Dries Van De Putte<fnref n="1" />, Karl D. Gordon<fnref n="1" />,<fnref n="2" />, Julia Roman-Duval<fnref n="2" />, Benjamin F. Williams<fnref n="3" />, Maarten Baes<fnref n="1" />, Kirill Tchernyshyov<fnref n="4" />, Brandon L. Lawton<fnref n="2" />, and Heddy Arab<fnref n="5" />

Received 2019 July 19; revised 2019 October 31; accepted 2019 November 6; published 2019 December 30

## Abstract

Photodissociation regions (PDRs) are parts of the ISM consisting of predominantly neutral gas, located at the interface between H II regions and molecular clouds. The physical conditions within these regions show variations on very short spatial scales, and therefore PDRs constitute ideal laboratories for investigating the properties and evolution of dust grains. We have mapped IC 63 at high resolution from the UV to the NIR (275 nm to 1.6 µm), using the *Hubble Space Telescope* WFC3. Using a Bayesian SED fitting tool, we simultaneously derive a set of stellar ($T_{eff}$, log(g), distance) and extinction ($A_V$, $R_V$) parameters for 520 background stars. We present maps of $A_V$ and $R_V$ with a resolution of 25 arcsec based on these results. The extinction properties vary across the PDR, with values for $A_V$ between 0.5 and 1.4 mag, and a decreasing trend in $R_V$, going from 3.7 at the front of the nebula to values as low as 2.5 further in. This provides evidence for evolution of the dust optical properties. We fit two modified blackbodies to the MIR and FIR SED, obtained by combining the $A_V$ map with data from *Spitzer* and *Herschel*. We derive effective temperatures (30 and 227 K) and the ratio of opacities at 160 µm to V band $\kappa_{160}/\kappa_V$ (7.0 × 10⁻<fnref n="4" /> and 2.9 × 10⁻⁹) for the two dust populations. Similar fits to individual pixels show spatial variations of $\kappa_{160}/\kappa_V$. The analysis of our *HST* data, combined with these *Spitzer* and *Herschel* data, provides the first panchromatic view of dust within a PDR.

*Unified Astronomy Thesaurus concepts: Interstellar medium (847); Photodissociation regions (1223); Interstellar dust extinction (837)*

## 1. Introduction

Photodissociation regions (PDRs) are regions of the interstellar medium (ISM) where the physical properties of the gas are mainly determined by the radiation field of nearby O or B stars. They form a boundary layer between ionized (H II) regions, and the rest of the molecular cloud where they reside. The FUV radiation from ionizing stars is quickly attenuated by the opacity of the relatively high-density medium shielding the rest of the gas (several 10 s up to ~10⁶ cm⁻<fnref n="3" />). As a consequence, there are typically sharp transitions between different regimes (e.g., H II, H I, and H₂), resulting in a layered structure (Hollenbach & Tielens 1997, 1999). This makes PDRs the perfect laboratory for studying the evolution of the ISM, including the dust over different physical and excitation conditions.

Interstellar dust grains make a very significant contribution to the total opacity of the ISM, and modify any impinging radiation field through the effects of absorption, scattering, and re-emission (Draine 2003; Steinacker et al. 2013). Aside from attenuating the UV radiation field that regulates the physics inside a PDR, they also couple to the gas through other processes. Dust grains play a major role in the formation of H₂, through reactions that take place on grain surfaces (Hollenbach & Salpeter 1971; Wakelam et al. 2017). They provide a major contribution to the

heating of the gas through the photoelectric effect, but can also have a cooling effect when gas-grain collisions occur (Bakes & Tielens 1994; Weingartner & Draine 2001).

To model these effects, knowledge about the composition and size distribution of the dust grains is necessary. While there are many dust models that can explain the observed extinction curves and emission spectra (Desert et al. 1990; Zubko et al. 2004; Draine & Li 2007; Compiègne et al. 2011; Jones et al. 2017), it remains difficult to accurately constrain the exact properties of the grains. Moreover, the effects of dust evolution processes can change these properties depending on the time and the environment. Some models, such as THEMIS (Jones et al. 2017), have built-in ways to follow the changes in dust properties. Providing better constraints on these models is crucial for understanding not only the dust itself, but also the structure and evolution of PDRs and the ISM in general.

Evidence for dust evolution in PDRs has been found through observations in the mid- and far-infrared. By applying Blind Signal Separation methods to *Spitzer* IRS data of several PDRs (Ced 201, NGC 7023 east and northwest, ρ Oph), Berné et al. (2007) identified two spectral shapes. One mainly contains the aromatic infrared bands (AIBs), and is linked to polycyclic aromatic hydrocarbons (PAHs). The other exhibits a combination of broad AIBs and MIR continuum emission, and was found to correspond to very small grains (VGSs). Later work by Boersma et al. (2014), uses k-means spectral clustering to identify zones in NGC 7023 that have similar spectral

<footnote n="1">¹Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281 S9, B-9000 Gent, Belgium; drvdputt.vandeputte@ugent.be</footnote>

<footnote n="2">²Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="3">³Department of Astronomy, Box 351580, University of Washington, Seattle, WA 98195, USA</footnote>

<footnote n="4">⁴Department of Physics and Astronomy, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA</footnote>

<footnote n="5">⁵Observatoire Astronomique de Strasbourg, Université de Strasbourg, 11 rue de l'Université, F-67000 Strasbourg, France</footnote>

<footnote>* Based on observations made with the NASA/ESA *Hubble Space Telescope*, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. These observations are associated with program GO-14186.</footnote>
shapes, and find a spatial evolution in the PAH band strength ratios.

Radiative transfer modeling by Compiègne et al. (2008) shows that differences in excitation conditions are not enough to explain the observed variations of the AIBs and MIR continuum in NGC 2023 and the Horsehead nebula, compared to the diffuse ISM. Changes in the relative abundances of the PAHs and VSGs are suspected. In the Orion bar, a radiative transfer model which uses the dust abundances of the diffuse ISM is sufficient to explain the dust variations derived from *Herschel*/PACS/SPIRE and *Spitzer*/IRAC maps, but overestimates the PAH emission at 3.6 µm (Arab et al. 2012). Using theoretical dust emission models which include the coagulation of grains, Köhler et al. (2015) were able to reproduce the changes in temperature, spectral index, opacity, and MIR emission, that are observed when transitioning from the diffuse ISM to high-density regions.

IC 63 is a nearby nebula which is illuminated by the B IV star $\gamma$ Cas. Based on the *Hipparcos* parallax (Perryman et al. 1997) of 5.32 mas, the distance to $\gamma$ Cas is roughly 190 pc. A comprehensive study of the physics and chemistry of IC 63 can be found in Jansen et al. (1994, 1995, 1996). A study using *Infrared Space Observatory (ISO)* data has characterized the fine structure and H$_2$ lines in IC 63 (Thi et al. 2009), and the distribution of PAHs has been derived from *Spitzer* data (Fleming et al. 2010). Recently, *Herschel* FIR maps and FIR spectroscopy data were combined with [C II] 157 µm velocity maps from the GREAT instrument on board the SOFIA observatory to revisit IC 63 in greater detail (Andrews et al. 2018).

For the radiation field at the tip of IC 63, the model of Jansen et al. (1995) used a value of 650 Draine field units (Draine 1978), or $G_0 = 1100$ in Habing field units (Habing 1968). Here, an edge-on orientation was assumed, meaning that the distance between the tip of IC 63 and $\gamma$ Cas is taken to be equal to the projected distance of 1.3 pc. In the work of Andrews et al. (2018), a value of $G_0 \sim 150$ is obtained instead, based on measurements of the FIR emission. This implies that IC 63 might be several times further away from $\gamma$ Cas than the projected distance, and that the orientation is not truly edge-on. The whole cloud has a projected size of $\sim 0.5$ pc, and the area studied in this work, the tip, is $\sim 0.1$ pc wide. IC 63 is still a dense PDR ($1.2 \times 10^4 \text{ cm}^{-3}$), but it has a relatively low column density ($2.3 \times 10^{20} \text{ cm}^{-2}$; Andrews et al. 2018). Compared to the Horsehead nebula or the Orion bar, this is about an order of magnitude lower in both number and column density, making IC 63 sufficiently transparent to detect background stars.

In this work we aim to study the spatial variations of the dust properties observed through extinction, in particular through the $A_V$ and $R_V$ parameters (Cardelli et al. 1989). IC 63 is a most suitable target for this, as it has many observable background stars, each of which provides information about the medium along a specific line of sight. Grouping these stars into spatial bins makes it possible to measure the average extinction for a number of regions on the sky. To measure the extinction for each star, we use the same approach as in Gordon et al. (2016) for the Panchromatic Hubble Andromeda Treasury (PHAT) survey data (Dalcanton et al. 2012). A catalog of point sources is generated from broadband *Hubble Space Telescope (HST)* observations, to which our Bayesian Extinction And Stellar fitting Tool (BEAST, Gordon et al. 2016) is applied.

In Section 2, we describe our seven-band photometric observations with *HST*, and how we extract photometric

measurements for 520 background sources. Our Bayesian extinction fitting tool is introduced, and some necessary modifications to it are explained. Section 3 presents the individual fit results for the sources, and how these were processed to create $A_V$ and $R_V$ maps. In Section 4, we compare our findings to earlier studies of IC 63, and use data from *Herschel* and *Spitzer* combined with our maps and modified blackbody fits to derive the $A_V$-normalized dust surface brightness and dust optical depth. To conclude, we present some simple per-pixel fits to the FIR SED, and compare some of the results for IC 63 with the Horsehead nebula and NGC 7023.

## 2. Data and Analysis

### 2.1. Hubble Observations

The photometric data for the various sources behind IC 63 were obtained through observations with the Wide Field Camera 3 (WFC3) of *HST*, using both the UVIS (UV and visual) and IR (infrared) channels. We obtained photometric images in the F275W, F336W, F475W, F625W, and F814W bands with the UVIS chip, and in the F110W and F160W bands with the IR chip. This set of seven broadband filters was chosen to cover the stellar SED from the UV to the near-infrared, and to optimize the extraction of the dust extinction parameters through SED fitting. The IR measurements provide a way to resolve the degeneracy between the stellar surface temperature and the reddening, while the UV filters constrain the type of extinction. A similar observing strategy, and a more detailed reasoning for choosing each filter can be found in Dalcanton et al. (2012).

The UVIS images span a field of view (FOV) of 162 by 162 arcsec, while the FOV for the IR images is somewhat smaller, at 123 by 136 arcsec (Figure 1). The observing program consisted of two visits, each two orbits long. During the first visit (2016 August 26), the F275W and F336W exposures were taken, followed by F814W and F160W. During the second visit (2016 August 31), the F625W, F110W, F475W, and additional F275W exposures took place. See Table 1 for the exposure times. The time difference of several days between the visits is intentional, to prevent persistence for the IR exposures.

Guard exposures for bright stars were taken for 5 s in F475W and F814W between the two orbits of the first and second visit, respectively. To help deal with macroscopic features of the detector and cosmic rays, a gap line dither pattern was used for all exposures with UVIS except the guard exposures. For the IR, a small line pattern was used instead to minimize blurring due to IR persistence. A post-flash was used for all UVIS exposures, except F625W.

For all figures that make use of this *HST* data (such as Figure 2), the *drz* images were used. These images are created by the standard “Drizzle” algorithm of the *HST* pipeline Fruchter & Hook (2002), which aligns and combines the dithered exposures, removes cosmic rays, and artifacts, and corrects for geometric distortion. The photometry described in Section 2.3 works with the *flt* files (for IR) and *flc* files (for UVIS). These contain the individual calibrated exposures. The *flc* files in particular have been corrected for the charge transfer efficiency (CTE) of the detector. See the WFC3 data handbook.<fnref n="6" />

<footnote n="6">⁶ http://www.stsci.edu/hst/instrumentation/wfc3/documentation</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Footprints of the observed data and retrieved archival data. The image (credit: Ken Crawford) shows IC 63 and its environment, including the nearby nebula IC 59, using a combination of exposures in blue, green, red, Hα, and SII filters. Between IC 63 and IC 59, a structure of dark dusty clouds can be observed, which appears more clearly in FIR images. The gray arrow indicates the direction of the radiation from γ Cas. The physical length scale shown assumes a distance of 190 pc.</figcaption></figure>

<table><caption>**Table 1**
Filters and Exposure Times as Listed in the Observing Program (Arab 2015)</caption><thead><tr><td>Channel and Filter</td><td>λ<sub>eff</sub></td><td>Exposure Time</td></tr></thead><tbody><tr><td>UVIS F275W</td><td>275 nm</td><td>2 × 600 s + 2 × 349 s</td></tr><tr><td>UVIS F336W</td><td>335 nm</td><td>2 × 650 s</td></tr><tr><td>UVIS F475W</td><td>475 nm</td><td>4 × 450 s</td></tr><tr><td>UVIS F625W</td><td>625 nm</td><td>4 × 420 s</td></tr><tr><td>UVIS F814W</td><td>814 nm</td><td>4 × 449 s</td></tr><tr><td>IR F110W</td><td>1.10 μm</td><td>2 × 350 s</td></tr><tr><td>IR F160W</td><td>1.60 μm</td><td>2 × 400 s</td></tr></tbody></table>

## 2.2. Ancillary Data

The FIR PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 2010) imaging data were obtained from the Herschel Science Archive. They were originally part of a key program for the Herschel Space Observatory (Pilbratt et al. 2010), “Evolution of Interstellar Dust” (Abergel et al. 2010). We also retrieved IRAC (Fazio et al. 2004) and MIPS (Rieke et al. 2004) photometric data from the Spitzer heritage archive, taken under the “Star Formation in Bright Rimmed Clouds” program (ID 202). An overview of the images can be found in Figure 2, where these data were reprojected onto the same coordinate frame as the drizzled HST UVIS images.

## 2.3. Point-source Photometry

### 2.3.1. Source Extraction

We employ a technique analogous to the one used to extract similar photometric catalogs from the PHAT (Dalcanton et al. 2012; Williams et al. 2014), SMIDGE (Yanchulova Merica-Jones et al. 2017), and METAL (Roman-Duval et al. 2019) data. Following a precise astrometry and alignment step, the photometry software DOLPHOT (Dolphin 2016) is used to automatically detect the point sources, and to determine their

positions and fluxes, based on point-spread function (PSF) fitting. The details about this photometric extraction routine are described in Williams et al. (2014). However, we do use the same tweaks as Roman-Duval et al. (2019). The Tiny Tim PSF libraries (Krist et al. 2011) are used instead of the Anderson libraries (Anderson & King 2000). Additionally, for UVIS, we make use of the flc images, which have already been corrected for CTE. Therefore, a separate CTE correction step after the PSF fitting is no longer necessary.

### 2.3.2. Removal of Spurious Detections

Due to the relatively bright extended emission in IC 63 (stellar light scattered by dust, gas emission lines due to recombination), the majority of the sources listed by the automatic source detection algorithm mentioned above are false positives. This can be seen in Figure 3, where most of the detections clearly coincide with the extended emission. Fortunately, we found some simple criteria to separate the bulk of these spurious detections from the real point sources, making use of some of the quantities provided by the PSF fitting routine. These quantities are given in certain columns of the photometric catalog.

The first criterion is a cut on the relative flux error in the F814W band. We only keep the sources for which the error to flux ratio F814W_ERR < 0.08 (S/N > 12.5). This excludes almost all of the entries due to the extended emission and some diffraction spikes of bright stars.

Second, an extra cut on the crowding is performed, specifically F814W_CROWD < 0.25. This removes a handful of overlapping sources, as well as several remaining detections located at diffraction spikes. The thresholds for these first two cuts were obtained through some trial and error. By visually inspecting plots such as Figure 3, we found that these values remove most of the spurious detections, while keeping detections of obvious point sources intact.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Overview of the stacked and drizzled *HST* images, and a selection of archived imaging data at longer wavelengths. The IRAC 4.5 and 5.8 µm images are not shown. The archived data has been reprojected onto the coordinate frame of the *HST* UVIS images. The color scales show the flux in MJy sr⁻¹.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Illustration of the cleaning step, using a zoom-in on the tip of the nebula. The orange boxes indicate the positions of detected sources in our catalog, before (top) and after (bottom) applying the criteria described in Section 2.3.2. All of the detections due to the extended emission are removed, while only a couple of the stars are missed.</figcaption></figure>

Lastly, we remove all sources that are outside the region where observations exist in all bands, to provide a homogeneous data set for analysis. The right panels of Figure 3 show the positions of the sources that are left after applying the cuts. There are 520 sources remaining, over the area covered by both the UVIS and IR chip.

## 2.4. Physics Model

To model the stellar SED attenuated by dust for each source, we employ the Bayesian Extinction And Stellar Tool (BEAST; Gordon et al. 2016). This tool works with a 7D grid of SED models, one-dimension for each parameter ($M_{ini}$, t, Z, d, $A_V$, $R_V$, $f_A$; see the description below and Table 2). For each set of parameters, the physics model makes a prediction for the observed SED of the source. This is done by first constructing models for the stellar spectrum and the extinction curve. The model spectrum is then extinguished according to this curve at each of its wavelength points, and integrated over the transmission curve of each filter.

### 2.4.1. Stellar Parameterization

In the current implementation, the models for the stellar spectra are calculated based on four parameters, specifying the star’s birth mass $M_{ini}$, age t, metallicity Z, and distance d.

Starting from the first three parameters, the luminosity of the star and the shape of the spectrum are calculated using a combination of publicly available stellar atmosphere grids (Castelli & Kurucz 2003; Lanz & Hubeny 2003, 2007) and evolutionary tracks (Marigo et al. 2008; Girardi et al. 2010; Bressan et al. 2012, 2013). The reasoning behind this choice of parameters can be found in Gordon et al. (2016). The predicted fluxes generated by this model are purely synthetic.

### 2.4.2. New BEAST Feature: Distance as an Extra Stellar Parameter

Prior to this work only a single value for the distance was supported, as this tool had only been applied in the extragalactic context where a single distance to all stellar sources can be assumed. For IC 63, setting a single distance is of course not possible, as all the background stars are located within the Galaxy, and hence they have a wide variety of distances. Therefore, in the version of the BEAST used for this work, the distance was implemented as an extra parameter for the stellar model.

### 2.4.3. Extinction Parameterization

The other three parameters describe the shape and magnitude of the extinction curve. The BEAST features a two-component dust extinction model, mixing two different shapes (wavelength dependencies); the type-$A$ curve, which models the average extinction as measured in the Milky Way (MW), and the type-$B$ curve, which models that of the Small Magellanic Cloud (SMC) Bar (Gordon et al. 2003). For the MW-like type, the shape of the extinction curve $A_λ/A_V$ depends on $R_V$

$$ \frac{A_λ}{A_V} = k_{λ,A}(R_V) \qquad (1) $$

where $A_V$ is the total Johnson V-band extinction by the dust, and $R_V = A_V/E(B - V)$. This $R_V$ dependence is modeled according to Fitzpatrick (1999). The SMC-like extinction model $k_{λ,B}$ does not depend on $R_V$, and its wavelength dependence is modeled using an average of the measurements given by Gordon et al. (2003). Within this model, the value of $A_V$ is a way to express the total dust column, while an increase in $R_V$ is believed to reflect a shift of the grain size distribution toward larger sizes (Cardelli et al. 1989). A third parameter, $f_A$, takes a linear combination of shapes $A$ and $B$:

$$ k_{λ}(R_V, f_A) = f_A k_{λ,A}(R_V) + (1 - f_A) k_{λ,B}. \qquad (2) $$

Note that this extra degree of freedom is necessary because the type-$A$ and type-$B$ curves only describe the average MW and SMC extinction. For individual sightlines, both components are needed to model the scatter on the $R_V$-dependent relationship of Cardelli et al. (1989). There are certain sightlines in the MW that exhibit SMC-like extinction in the UV (Valencic et al. 2003), as well as sightlines in the SMC with MW-like extinction (Gordon et al. 2003).

This extinction model is then applied to the stellar spectra generated from the first four parameters.

### 2.4.4. Model Details for This Work

The details about the parameter ranges and their spacing in the BEAST physics model grid are summarized in Table 2. Suitable ranges and resolutions were obtained through a combination of previous experience (such as Table 1 in Gordon et al. 2016), some
<table><caption>**Table 2**
Specification of the Parameter Grid for the Physics Model</caption><thead><tr><th>Parameter</th><th>Description</th><th>Min</th><th>Max</th><th>Resolution</th><th>Prior</th></tr></thead><tbody><tr><td>log(t) (yr)</td><td>stellar age</td><td>6.0</td><td>10.13</td><td>0.3</td><td>constant SFR</td></tr><tr><td>log(M<sub>ini</sub>) (M<sub>☉</sub>)</td><td>stellar mass</td><td>-1.0</td><td>2.5</td><td>variable<sup>a</sup></td><td>Kroupa IMF</td></tr><tr><td>log(Z) (mass fraction)</td><td>stellar metallicity</td><td>-3</td><td>-1.2</td><td>0.15</td><td>flat in Z</td></tr><tr><td>log(d) (pc)</td><td>distance</td><td>2.7</td><td>4.2</td><td>0.037<sup>b</sup></td><td>flat in log(d)</td></tr><tr><td>A<sub>V</sub> (mag)</td><td>dust column</td><td>0</td><td>10.05</td><td>0.08</td><td>flat</td></tr><tr><td>R<sub>V</sub></td><td>dust average grain size</td><td>1.5</td><td>5.5</td><td>0.5</td><td>flat</td></tr><tr><td>f<sub>A</sub></td><td>dust mixing parameter</td><td>0</td><td>1</td><td>0.2</td><td>flat</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Supplied by isochrone model grid.

<sup>b</sup> This corresponds to 40 logarithmic steps from 500 to 15000 pc.

trial and error, and some compromises considering computational resources. Assuming a distance to IC 63 of about 190 pc, we initially considered a lower bound of 150 pc for the new distance parameter. A test run for the fitting showed that only a handful of sources were closer than 500 pc. Therefore, we decided to change the lower bound for the distance to 500 pc, providing a slight increase in resolution. Analogously, we confirmed that an upper bound of 15,000 pc was sufficient to fit all observed stars.

The parameters shown in the table give rise to about $5.3 \times 10^9$ models, and without compression, about 1.6 TB of disk space is needed to store this grid. The resolution for most of the parameters is relatively low, because a relatively large range and number of bins is needed for the distance. For our purposes, however, a relatively rough estimate of $A_V$ and $R_V$ suffices, since we only need the average values over relatively large pixels. Within such a pixel, we found the range of $A_V$ and $R_V$ values to be broader than the chosen resolution. This spread can either be caused by the different stellar distances resulting in a different contribution by the diffuse Galactic ISM, or by small-scale variations in the cloud. The spread on $A_V$ and $R_V$ within each pixel will be used to estimate the error of the mean values. In any case, the precision on the $A_V$ and $R_V$ maps is mostly limited by the number of sources and the chosen pixel size, and not by the precision of individual $A_V$ or $R_V$ measurements.

The priors are shown in the rightmost column of Table 2. For the new distance parameter, the prior is chosen to be flat as a function of log(d). As a consequence, for some stars a degeneracy exists between the estimated distance and mass (or luminosity): an observed star can either be far away and of a luminous type, or nearby and of a fainter type. The multiband photometry partly resolves this degeneracy, and we did not find this effect to be problematic for constraining the extinction.

In the future, more physically motivated priors could be implemented for the distance, such as an exponentially decreasing space density (Bailer-Jones et al. 2018). Another option would be to make direct use of the parallaxes provided by Gaia DR2 (Gaia Collaboration et al. 2016, 2018; Luri et al. 2018), and use them as source-specific priors in a post-processing step. However, since there are only 59 out of ~500 stars in our data set with Gaia data available, this would have a limited impact on the results of this work. Moreover, we did cross-check the Gaia data with our fit results, and found that the parallaxes derived from the distances were consistent with the fairly uncertain Gaia parallaxes (Figure 4).

## 2.5. Noise Model

### 2.5.1. Artificial Star Tests

By applying the physics model to all combinations of the seven parameters, a grid of theoretical SEDs is constructed, representing the observations in the relevant filters under perfect circumstances. But in a realistic setting, when the flux is measured from imaging data, the flux values that are extracted for a certain point source are not only affected by shot (Poisson) noise and the PSF, but also by nearby objects, various sources of background, instrumental artifacts, and the extraction algorithm itself. Nearby point sources can create crowding effects (Gordon et al. 2016), which increase the noise and cause an over-estimation of the flux. However, since the density of observed stars in IC 63 is very low compared to the size of the PSF, these effects are negligible. For IC 63, the main expected contributor to the observation model is the presence of a foreground and background of extended emission, such as H$\alpha$ emission of the gas or stellar light scattered by dust grains.

To quantify how significantly one of the model SEDs differs from the observed SED of a certain star, we create a model describing the uncertainty (error $\sigma$) and systematic deviation (bias $\mu$) for each flux. While the photometric catalog produced by the PSF fitting does contain values for the uncertainty on each flux (which we used to perform the input cleaning step), these do not cover all possible error sources. Instead, we treat the errors as a property of each theoretical model: for each grid point $i$, we aim to calculate the error or uncertainty $\sigma_i$ and a bias $\mu_i$ in every band.

The starting point for this calculation is a set of artificial star tests (ASTs). The input for a single AST consists of a theoretical SED and a position. These parameters are used as input for the fake star mode of DOLPHOT. This routine will insert a fake star with the given SED into each of our observations, simulate the PSF and photon noise, and then perform the exact same photometric extraction routine. The result is a set of output fluxes, representing a mock-observation of a point source with that specific theoretical SED. By putting the same source at many different positions, statistical information is obtained about the deviations that occur when observing that specific SED $i$, from which appropriate values for $\sigma_i$ and $\mu_i$ can be derived.

Ideally, such a set of tests would be performed for each SED $i$ in the physics grid, leading to an individual measurement of the noise for each model. However, performing a single AST takes about 2 minutes of computing time on a modern
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Comparison of the Gaia DR2 parallaxes and the parallax derived from the distance expectation values $d_{\text{exp}}$, for the 59 sources present in both catalogs.</figcaption></figure>

processor, so doing this for the billions of models contained in the physics grid is practically impossible. Therefore, it is customary to produce a sample of SEDs that is representative of all the models relevant for fitting the observed sources.

### 2.5.2. New BEAST Feature: Uniform SED Sampling by Magnitude

Since the distance parameter can strongly scale the model SEDs, faint stars at long distances and luminous stars at short distances will have model fluxes that are far away from the minimum and maximum of the observed catalog. Their noise parameters are not useful for fitting our observed SEDs. Therefore we only consider models that have fluxes that fall between the minimum and maximum of the observed SEDs, with 1 magnitude of leeway to allow for a range of uncertainty. The magnitude ranges of the remaining models are then discretized into 25 bins for each filter, and SEDs are chosen randomly from the physics grid and added to the AST input list. For each SED that is chosen, the corresponding magnitude is determined in each band, so we can keep track of the number of ASTs that will cover each filter-mag bin. The algorithm continues until we have 50 samples in each bin, and the resulting set of AST input SEDs will have a magnitude distribution that is more or less flat in each filter.

Figure 5 shows the difference between the original method which evenly samples by stellar age, and the new method which evenly samples by magnitude. To show the magnitudes that need to be covered, the distribution of the observed catalog is also shown. The original method does not sufficiently sample the low brightness range, while simultaneously producing many samples that are several magnitudes brighter than the maximum of the observed catalog. As designed, the new method produces a certain minimum number of samples for each flux bin, within a more suitable range. The peaks that appear in the distribution depend on the contents of the physics grid and the brightness cutoffs. The extra stars that populate those peaks were picked to fill the rarer flux ranges in other filters. Note that the new method still produces many samples outside of the observed flux range, because we only require that

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Red: magnitude distribution produced by the new sampling method, which aims to provide a minimum number of samples for each magnitude interval. Blue: the old sampling method, which simply picks a fixed number of SED models from each stellar age bin. Black outline: magnitude distribution of the observed catalog (only fluxes > 0 shown).</figcaption></figure>

at least three bands fall within that range. This allows for a larger variety of SED shapes to be picked for the ASTs.

Using this approach, a list of about 2500 SEDs was generated. This list is duplicated multiple times while giving a random position to each SED, leading to 180,000 unique ASTs. The list of ASTs was split up into a set of jobs, each of which have a manageable runtime. These jobs were executed across three machines with the help of GNU parallel (Tange 2011), using 78 of the available cores over the course of 3 days. The results were then merged and processed into a single fake star catalog, of the same format as the one containing the observed photometry.

### 2.5.3. From ASTs to Noise Model

Before the fake star catalog is further processed by the BEAST to generate the noise model, the same selection function used in Section 2.3.2 is applied. This way it is ensured that the observation model generates a set of noise parameters that is representative of our cleaned ensemble of observed stars. Note that this also removes ASTs which are missing a flux in one or more bands, which includes those positioned outside the IR chip's FOV. After this step, about 80,000 fake stars remain. Although a large fraction (100,000 out of 180,000) of the fake stars ended up either undetected or rejected, we confirmed that
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Flux-dependent model of the bias $\mu$ and uncertainty $\sigma$ for the shortest and longest wavelengths. The vertical axis is a symmetrical log scale (linear between $10^{-3}$ and $10^{-3}$), which shows the sign change of $\mu$. The absolute value of $\mu$ is also plotted, showing that the bias correction is within $1\sigma$ almost everywhere.</figcaption></figure>

the remaining catalog still sufficiently covers the required magnitude range, as in Figure 5.

To generate the noise model for each SED in the physics grid, the BEAST’s most conservative method is used, which considers each flux individually and as such ignores any correlations between the bands. This method minimizes the number of ASTs needed, which is already quite large due to the broad range of fluxes needed to model all the distances. Similarly to how the AST input SEDs were chosen, the flux ranges of the model SEDs are divided into 30 logarithmic bins. For each bin, the input and output SEDs of the relevant ASTs are compared, producing statistical values for the error and the bias in that flux range, for a specific filter. By interpolating between these bin-averaged values, each filter now has a model that predicts $\sigma$ and $\mu$ as a function of the theoretical flux $F$. This behavior is shown in Figure 6.

One of the main features of the noise model is that the bias seems to be negative for all but the brightest fluxes, in all bands. This means that the fluxes extracted by the photometry routine are on average smaller than the input fluxes of the fake stars. To take this effect into account, we forward model the observed fluxes by adding the value of $\mu$ to the theoretical model $F$, as shown in Equation (6) of the next section. In this case specifically, the negative value for $\mu$ makes it possible to fit our model to the many stars with negative flux measurements in the UV. The theoretical model cannot produce negative flux values, but adding the bias shifts these values down to the correct level.

For weak fluxes, the constant slope of the relative $\sigma$ and $\mu$ implies that the noise parameters are dominated by photon noise and/or constant features in the image. For stronger fluxes, the relative error levels off while the bias becomes positive (but much smaller than the error). This means that fractional deviations dominate instead, which are most likely caused by the PSF or imperfections in the photometry routine. For all bands except the UV band, the absolute value of $\mu$ seems to be smaller than $\sigma$. This means that the exact value of $\mu$ will have less of an effect on the results, as it only causes deviations within $1\sigma$.

Since we suspect that the background/foreground extended emission provides the main contribution to the noise model, we experimented with a noise model that is stratified according to the intensity of this emission. In the Appendix, we explain how we determine 4 regions of similar background/foreground intensity, and create a separate noise model for each of them. Ultimately, the fitting results obtained using these individual noise models were not found to be significantly different, so the rest of the analysis was performed using the single noise model presented in this section. However, these tools integrate well with and improve upon the existing code that deals with varying levels of crowding, and have therefore been incorporated in the main BEAST branch.

## 2.6. Fitting Procedure and Output

Once the necessary models have been set up, the fitting step comes down to calculating a 7D posterior probability distribution

$$P(\theta|F_{\text{obs}}) \propto P(F_{\text{obs}}|F(\theta), \Sigma(\theta))P(\theta) \quad (3)$$

for each of the observed point sources. Here, $F_{\text{obs}}$ is the observed SED, and $F(\theta)$ and $\Sigma(\theta)$ represent the SED and the noise parameters for one of the models, evaluated at a specific point $\theta$ of the 7D parameter grid. The different priors listed in Table 2 are multiplied to form the 7D prior $P(\theta)$. The grid point $\theta$ at the global maximum of the posterior distribution is considered to be the best-fitting model. The likelihood function is a 7D multivariate Gaussian:

$$P(F_{\text{obs}}|F(\theta), \Sigma(\theta)) = \frac{1}{N} e^{-\chi^2/2} \quad (4)$$

$$\text{with } \chi^2 = \Delta(\theta)^T \mathcal{C}(\theta)^{-1} \Delta(\theta) \quad (5)$$

$$\text{and } \Delta(\theta) = F_{\text{obs}} - \mu(\theta) - F(\theta) \quad (6)$$

where the dimension $N=7$ equals the number of fluxes, and $\mathcal{C}(\theta) = \text{diag}(\sigma(\theta))$ and $\mu(\theta)$ are the covariance matrix and the bias vector. The values for $\sigma(\theta)$ and $\mu(\theta)$ are calculated as a function of $F(\theta)$, as described in Section 2.5.3. Note that $\mathcal{C}$ is diagonal for this work, since this calculation is done for each band individually.

In principle, the whole 7D log-likelihood can be written to disk for each star. In practice, only a local area around the peak of the posterior distribution is written out. Additionally, the 1D marginal probability distributions for each parameter are also available. For our purposes, however, some statistics derived from the posterior are sufficient. For all the parameters, the best fit, the expectation values, and the 16th, 50th, and 84th percentiles are calculated and stored in a table that contains one row per source. The $\chi^2$ of the least-squares model and $\ln P$ of the best-fitting model are also stored. The former gives us some
insight about the quality of the grid and the noise model, and the source itself. A high $\chi^2$ usually means that the grid range or resolution is not sufficient, that there might be something wrong with the noise model, or that a source in particular has an SED that is very much unlike that of an extinguished star.

## 3. Results

### 3.1. Post Processing

Following the workflow described in the previous section, we obtained a catalog for all 520 sources (Figure 3) that were left after applying the criteria from Section 2.3.2. This catalog contains statistics for the seven main parameters, as well as a set of derived stellar quantities: current mass $M_{\text{act}}$ (derived from age and initial mass $M_{\text{ini}}$), radius $R$, luminosity $\log(L)$, surface gravity $\log(g)$, temperature $\log(T)$, and absolute bolometric magnitude $m_{\text{bol}}$. The distribution of the expectation values of the seven parameters and the correlations between them are shown in Figure 7, and the $\chi^2$ distribution is shown in Figure 8. The latter seems to be slightly peaked around 7, which could naively be attributed to the fact that we have measurements in seven bands. But in practice, the distribution of this $\chi^2$ is not as simple to interpret, because there are significant correlations between the observed fluxes in different bands. The value can only be used as a rough check that the fits are reasonable.

The two bottom panels of Figure 8 show that for almost all sources with a higher $\chi^2$ (30–200), the width of their 1D posterior $A_V$ and distance distribution is near the minimum cutoff. This means that the resulting fits are acceptable despite their high $\chi^2$ values; their accuracy is simply limited by the resolution of the grid. Decreasing the $A_V$ and distance spacing will likely improve the $\chi^2$ and make the $A_V$ distributions narrower, but we consider the precision achieved with the current parameters sufficient for the analysis that follows. Since the distance is a nuisance parameter which is integrated over to obtain the $A_V$ and $R_V$ expectation values, increasing the distance resolution will arguably have a negligible effect.

Based on the corner plot of Figure 7, there seem to be no correlations or degeneracies that could cause problems. The $R_V-f_A$ relationship is to be expected, as the type-B dust extinction component has a fixed $R_V = 2.74$ (Gordon et al. 2016); the results with low $f_A$ have an $R_V$ which is closer to this value. It is however remarkable that the population of stars seems to have a bimodal distribution of the fitted initial mass $M_{\text{ini}}$. Figure 9 shows a color-magnitude diagram, and it can be seen that the models without extinction seem to have a gap along the color axis. Upon closer inspection of the derived parameters, the high-mass group centered at around $0.8 M_\odot$ has temperatures ranging from 4300 to 6200 K, which corresponds to main-sequence stars. For the low-mass group at $0.3 M_\odot$, the temperatures range from 2900 to 3600 K, which means it entirely consists of red dwarfs. Looking at the fitted distances in Figure 7, the red dwarfs are only observed from 1000 up to 7300 pc, while the more massive, hotter stars are distributed across the entire range from 1500 to 15,000.

To check whether this feature could be produced by two different populations, we created a simulated catalog of stars around the line of sight toward IC 63, using the Besançon Galaxy model (Robin et al. 2003). When no magnitude limit is used (AbsMag_max = 99, all observed band limits at 99), the resulting log-mass distribution forms a broad peak between roughly 0.05 and $1.4 M_\odot$. At the center of this peak ($0.4 M_\odot$), a

small dip does appear, of about 10% of the height of the distribution. We experimented with some magnitude cuts for the generated catalog, and these generally caused the dip in the log-mass distribution to become less pronounced. Concluding, it is not out of the question that the observed mass bimodality could be caused by the intrinsic mass distribution along the line of sight. While we suspect the observational selection function and possible fitting artifacts to have a significant effect, it remains unclear why exactly there is such a prominent gap in the distribution of the fitted masses.

We perform a second cleaning step, which removes sources from the output catalog for which the obtained parameters are likely to be incorrect. For certain sources, the shape of the observed SED does simply not resemble that of a star with a certain extinction. The reason for this can be physical, because it is a type of star not covered by the BEAST model, or because there is some other phenomenon at play which is not common enough across the image to be captured by the noise model (e.g., unusually strong scattered light, background galaxy, diffraction spikes). Most of these are easy to spot, with $\chi^2$ usually being several hundreds or thousands. Based on the histogram shown in Figure 8, we apply a cutoff of $\chi^2 < 100$ to remove the worst offenders. They have been highlighted in red in Figure 7.

There also seems to be a group of stars for which the $\chi^2$ is not necessarily large, yet their $A_V$ is unusually high ($\gtrsim 7$ mag). When inspecting the other parameters for these stars (Figure 10), they all seem to have high temperatures ($\log T \gtrsim 3.8$), luminosities ($\log L/L_\odot \gtrsim 1$), masses ($\log M_{\text{ini}} \gtrsim 1.3$), and distances ($d \gtrsim 10,000$ pc). This likely points to a type of star that the BEAST fails to fit correctly, due to the necessary stellar models not being present. The resulting parameters might provide an acceptable SED, but their values are likely unphysical. In the second panel of Figure 10, a gap can be seen at a temperature of $\log T = 3.87$. We choose to place a cut on the temperature there; all stars that have an expectation value or best-fit value of $\log T > 3.87$ are deemed suspicious, and are removed from the catalog. The stars that were removed using this cut are highlighted in gray in Figure 7.

This phenomenon has some similarities with a previously documented failure case for the BEAST, which occurs when an attempt is made to fit the SED model to a Thermally Pulsing Asymptotic Giant Branch (TP-AGB) star.<fnref n="7" />

### 3.2. $A_V$ and $R_V$ Maps

The main objective of this work is to look for signs of dust evolution, by mapping the spatial variations of the $A_V$ and $R_V$ parameters. The maps, shown in Figure 11, are constructed by creating a rectangular grid, aligned with the R.A. and decl. axes. The angular size of each map pixel is roughly 25", which provides us with about 10–20 sources per pixel (except at the edges, as indicated with a green n on the figure). The $A_V$ and $R_V$ values for these pixels are calculated by considering all the stars that fall within their boundaries, and taking the median of the expectation values. Some of the pixels have no sources in them at all, mainly due to geometry effects, as the grid does not line up with the edges of the IR data. These pixels are indicated with a red X. The error for each pixel is estimated by computing the standard deviation for the same set of values, and dividing it by the square root of the number of sources $N$.

<footnote n="7">⁷ https://beast.readthedocs.io/en/latest/beast_issues.html</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Distribution and correlation of the expectation values for the fit parameters, of all 520 sources. The density of the points is represented by the color scale. Red: stars with $\chi^2 > 100$. Gray: stars with $\log T > 3.78$. After applying these cuts, 462 stars remain.</figcaption></figure>

The typical error is ~0.25 mag for $A_V$, and ~0.15 for $R_V$. To show the size of the errors versus the size of the variations, cuts through these maps are shown in Figure 12. The numeric values of the maps are given in Tables 5 and 6 in the Appendix.

To highlight the variations of $A_V$ ($R_V$) across the nebula, some contours are shown that have been calculated directly from the maps. The contour levels for the maps were chosen by taking four linearly spaced values between the 16th and 84th percentiles of the total range of $A_V$ ($R_V$) values. The same contours are plotted onto the SPIRE 350 µm image, to indicate the spatial correspondence between features of the $A_V$ and $R_V$ maps and those of the FIR emission.

It should be noted that $A_V$ does not drop to zero outside of the visible gas of the nebula. In the area at the bottom of the map, just outside the visible edge of the gas, an $A_V > 2$ is still observed while the 350 µm emission reaches a minimum (near the background level of ~22 MJy/ sr). A selection of 50 stars in this area, has an $A_V$ distribution centered at ~2.5 mag, with a spread of 0.70 mag. They have a mean $A_V$ of 2.51 ± 0.10, and

we assume that this is the extinction produced by material in the background (or foreground) of IC 63.

This simplistic method of accounting for the diffuse Galactic extinction comes with an important caveat, especially since this value is larger than the residual $A_V$, i.e., the extinction by the nebula. Since the stars are at widely varying distances (1–15 kpc), the extinction between each star and the cloud could be very different. To investigate whether this could contaminate the structure of the $A_V$ map, we examined the distance distribution in each pixel by making maps of the 16th, 50th, and 84th percentiles of the observed set of distances: $d_{16}$, $d_{50}$, and $d_{84}$. These maps did not show any obvious structure that could drastically contaminate the observed trends in $A_V$. This was further confirmed by examining scatter plots (Figure 13) of the $A_V$ pixels versus the corresponding distance distribution percentiles. Furthermore, the range of $d_{50}$ values for the pixels is quite constrained, with almost all the values being between 4 and 6 kpc, with more or less a flat distribution. Similarly $d_{16}$ lies between 2000 and 4000 for most pixels, and
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Top panel: $\chi^2$ histogram for the fit results and the theoretical distribution for $\log(\chi^2)$ that would naively be expected for 7 degrees of freedom. Bottom two panels: relationship between $\chi^2$ and the width of the 1D posterior distributions for $A_V$ and $\log(d)$ (defined by their 84th and 16th percentiles). The color scale visualizes the density of overlapping points on the plot. The dashed line indicates $\chi^2 = 30$. Almost all sources that have a $\chi^2$ higher than 30 have narrow posterior distributions for $A_V$, $\log(d)$, or both.</figcaption></figure>

$d_{84}$ lies between 5000 and 8000. Therefore, the diffuse Galactic extinction introduces an extra source of noise caused by variations in the average distance, but we do not expect it to produce false trends.

# 4. Discussion

## 4.1. $A_V$

The $A_V$ map shows that the extinction peaks at $00^{\prime\prime}59^{\prime\prime}02^{\prime\prime}$; $60^{\prime\prime}$ 53′30″, with $A_V = (3.91 \pm 0.23)$ mag. This is near the tip of the nebula. Given the background $A_V$ estimate of the previous section, the average amount of extinction added by the tip of IC 63 is hence

$$ \Delta A_{V,\text{tip}} = (1.41 \pm 0.25)\text{mag.} \quad (7) $$

The higher $A_V$ values seem to lie along a ridge-like feature containing this maximum, more or less parallel to the direction of the radiation (see arrow in Figure 1). The extinction of this part in the nebula typically ranges from 1.0 to 1.4 mag. With the exception of the bottom row of the map which we used to

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Color–magnitude diagram of the (nonnegative) raw data and the expectation values of the theoretical model fluxes. The dashed lines indicate the median magnitude and color of each set. Bright sources seem to have colors and magnitudes that are closely matched, as can be seen from the nearly overlapping crosses and diamonds. For weak sources, this overlap disappears due to noise.</figcaption></figure>

measure the $A_V$ background, we find a minimum $A_V$ of 3.0 mag, or a net extinction of at least 0.5 mag for the rest of the values.

In the top right of the map, a feature with a much higher extinction ($A_V \sim 5$ mag) can be seen, 2.5 mag when subtracting the background. The cuts in Figure 12 show that this peak is significantly more pronounced than the other features. On the full version of the SPIRE map, this feature is a bright part of a larger structure in the FIR, which does not have optical emission. In the NIR (e.g., F160W in Figure 2), it can be seen very faintly. We do not consider this structure part of the PDR, and we ignore its contribution in the analysis that follows.

## 4.2. $R_V$

The $R_V$ measurements over individual stars have a spread of 0.5, and a mean of $R_V = 3.556 \pm 0.001$ (the error on the mean is this small due to the number of sources in our sample, but excludes systematics). The contours in Figure 11 seem to indicate a slight drop in $R_V$, going from ~3.6 at the front, to ~3.4 at the back of the PDR. Given the errors on the map of about 0.15, this drop is not significant when comparing individual pixels. Instead, we first examine the diagonal cuts in Figure 12. For both the orange and the blue curves, a slope is observed, and the difference between the second and fourth points is larger than the error bars. Lastly, to quantify the significance of the $R_V$ decrease, we focus on the two magenta regions shown in the same figure. We gather the individual $R_V$ results for the stars in each of the two regions, and calculate the means. For the front region (bottom right rectangle in Figure 12), we find $R_V^{\text{front}} = 3.63 \pm 0.08$, and for the back region $R_V^{\text{back}} = 3.42 \pm 0.05$. A standard two-sample t-test comparing the two sets of $R_V$ values yields $t = 2.25$, with a
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Illustration of the suspicious, high-A_V fit results described in Section 3.1. Every dot corresponds to the fit result for a single star. Red: stars removed by applying cuts on the temperature (log T_exp and log T_best < 3.87). Most of these high temperature stars also have high luminosities (middle panel), as well as high distances and masses (bottom panel). This is suspicious, and despite their seemingly normal χ² values (top panel), these fit results are likely to be incorrect.</figcaption></figure>

two-tailed p-value of 0.02, indicating that the mean R_V differs significantly between the front and the back region.

One also has to take into account that the observed R_V is a mix of that of the nebula and the background. Since the A_V of the background (~2.5 mag) is larger than that of the nebula (0.5–1.4 mag), the drop in R_V of the dust in IC 63 will be stronger than what is observed. We can write the measured R_V in terms of the extinction parameters of the two parts: (A_V^B, R_V^B) for the background, and (A_V^N, R_V^N) for the nebula

$$R_V = \frac{A_V}{E(B - V)} = \frac{A_V}{\frac{A_V^N}{R_V^N} + \frac{A_V^B}{R_V^B}} \quad (8)$$

where we use the fact that the measured E(B-V) is the sum of the two contributions. Isolating R_V^N from this equation, we obtain

$$R_V^N = \frac{A_V^N}{\frac{A_V}{R_V} - \frac{A_V^B}{R_V^B}} = \frac{A_V - A_V^B}{\frac{A_V}{R_V} - \frac{A_V^B}{R_V^B}} \quad (9)$$

Knowing A_V and R_V from our maps, and using our value of A_V^B = 2.51 mag obtained earlier, we just need to assume a value for R_V^B. It should be noted, however, that Equation (9) is quite sensitive to the assumed value of R_V^B. For example, picking a value near the Galactic average of ~3.0 often leads to very high or even negative values for R_V^N. If we take the average R_V over the same region that we used to determine A_V^B, we obtain R_V^B = 3.6.

The galactic R_V map by Schlafly et al. (2017) shows values between 3.0 and 3.7, with variations at kiloparsec scales. While there might still be variations at smaller scales, we assume that A_V^B and R_V^B stay constant over our FOV. With these assumptions, we find that R_V^N ranges from 3.7 at the bottom right of the map (tip of the nebula), to values between 2.5 and 3.3 near the top left. Since R_V correlates with the average grain size, this points to the existence of processes making the average grain size of the dust population larger at the front of the nebula. Possible candidates are coagulation and accretion at the tip of the nebula due to higher gas densities (Stepnik et al. 2003; Köhler et al. 2015).

Andersson et al. (2013) determined A_V and R_V for 14 observed background stars of IC 63. This was done using a spectral classification of the stars, followed by a calculation of R_V using the relation R_V = 1.12 × E(V-K)/E(B-V) + 0.02 (Fitzpatrick 1999). Via the total-to-selective extinction and the color excess E (B-V), A_V was then obtained. Unfortunately these stars are not present in our sample, as they were either too bright (and removed per Section 2.3.2), or simply outside of the FOV. They obtained an average value 〈R_V〉 = 2.2 ± 0.5, which is significantly lower than the values listed by our R_V map. We do not believe that this is due to different physical conditions, as the positions of the stars used in Andersson et al. (2013) probe different areas of the cloud: some are at the front edge, and others are behind or inside the cloud. The low R_V values listed for all of these stars, and the fact that there is no overlap between the samples, make it hard to compare these results.

### 4.3. Column Density to A_V Ratio

Our results have an average A_V across IC 63 of 3.3 mag, or 0.8 mag with the background subtracted. As a check on the validity of this measurement, we make an estimate of the N_H (column density of hydrogen nuclei) to A_V ratio.

The relationship between the extinction A_V and the column density N_H of hydrogen has been found to be linear, on average in the MW. Different values for the ratio 〈N_H/A_V〉 have been found, by using far-UV extinction observations (Bohlin et al. 1978; Diplas & Savage 1994) and observations of X-ray sources (Reina & Tarenghi 1973; Gorenstein 1975; Predehl & Schmitt 1995; Güver & Özel 2009). The most recent measurements by Zhu et al. (2017) show that 〈A_V/N_H〉 is more or less invariant across the Galaxy, and N_H = (2.08 ± 0.02) × 10²¹ A_V cm⁻² mag⁻¹.
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Left: maps of $A_V$ and $R_V$ and their contours. The contour levels are indicated on the respective colorbars. Middle: errors on the map values. Right: the same contours overplotted onto a SPIRE 350 µm image. Red X: No sources. Green n: low ($\le 5$) number of sources.</figcaption></figure>

The HI4PI 21 cm survey (HI4PI Collaboration et al. 2016) contains two sightlines which cover our area of interest; the listed H I column densities are $4.47 \times 10^{21} \text{ cm}^{-2}$ and $5.53 \times 10^{21} \text{ cm}^{-2}$. We assume that the average atomic hydrogen column density for IC 63 is somewhere between these two values: $N_{H1} = (4.50 \pm 0.03) \times 10^{21} \text{ cm}^{-2}$. Since we cannot distinguish the 21 cm contribution of IC 63 from that of the background, we assume that the $N_{H1}/A_V$ ratio is the same for IC 63 and its background, and we divide by the average $A_V^{\text{total}} = 3.3 \text{ mag}$.

$$ \frac{N_{H1}}{A_V} = \frac{N_{H1}^{\text{total}}}{A_V^{\text{total}}} = 1.36 \times 10^{21} \text{ cm}^{-2} \text{ mag}^{-1}. \quad (10) $$

Andrews et al. (2018) used H$_2$ excitation diagrams based on IRS spectra from *Spitzer* to determine the gas temperature and column density for the molecular component. They found that it consists of two components, of which the cooler one has a temperature of $T = (207 \pm 30) \text{ K}$ and a column density of $N_{H2} = 2.3 \times 10^{20} \text{ cm}^{-2}$, while the warm component has $T = (740 \pm 47) \text{ K}$ and $N_{H2} = 9.3 \times 10^{17} \text{ cm}^{-2}$. Given the dominant (cold component) value from Andrews et al. (2018) for the column density of H$_2$, and our background subtracted value for the average $A_V = 0.8 \text{ mag}$, we find a $2N_{H2}/A_V$ ratio of

$$ \frac{2N_{H2}}{A_V} = 5.75 \times 10^{20} \text{ cm}^{-2} \text{ mag}^{-1}. \quad (11) $$

Since we are looking at the PDR, we expect the medium to be mostly neutral. IPHAS data (Drew et al. 2005;

Barentsen et al. 2014) show that most of the H$\alpha$ emission is emitted at the bright, irradiated edge of the cloud. Additionally, the PDR model of Jansen et al. (1995) predicts that the H$^+$ abundance is at least six orders of magnitude lower than that of H or H$_2$. Therefore, we will not take into account the column density of ionized hydrogen here. The combined value is

$$ \frac{N_H}{A_V} = \frac{N_{H1}}{A_V} + \frac{2N_{H2}}{A_V} = 1.94 \times 10^{21} \text{ cm}^{-2} \text{ mag}^{-1} \quad (12) $$

This value is quite close to the value from the literature mentioned above. Looking for a separate velocity component in the H I data, and combining our measurements with better MIR spectroscopy, would provide a way to map the molecular fraction and the gas to dust ratio for IC 63.

### 4.4. $A_V$ to Surface Brightness Relation

We investigate how our $A_V$ and $R_V$ maps compare to the dust emission measured using *Spitzer* and *Herschel* data. To look for a relation between $A_V$ and the surface brightness in the different bands, the IRAC, MIPS, PACS, and SPIRE maps are first reprojected onto the same coordinate grid, using the *reproject* package for Python.<fnref n="8" /> For each reprojected image, the pixels are matched with the ones from the $A_V$ map, and the reprojected fluxes are plotted against the $A_V$ map values, as shown in Figures 14 and 15. For all bands, we find a positive

<footnote n="8">⁸ https://reproject.readthedocs.io/en/stable/</footnote>
<figure><img src="image_13.png" /><figcaption>Figure 12. Cuts through the $A_V$ and $R_V$ maps. Arrows: starting point and direction of each cut. Curves and error bars: value and error along each arrow, color coded. Different offsets have been added on the x-axis for legibility reasons. The magenta rectangles indicate the two regions used to determine the significance of the $R_V$ decrease.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13. Histogram and scatter plot for the distance percentiles in each pixel. Each point corresponds to one of the pixels of the $A_V$ map. The plot for $d_{84}$ is not shown here, but shows a similar scatter.</figcaption></figure>

correlation between $A_V$ and the observed flux. Upon applying the same analysis to the $R_V$ map, no significant correlation between any of the fluxes and $R_V$ was found.

For each panel in Figures 14 and 15, we determine a slope and an uncertainty using a linear fit, which is weighted according to the errors on the $A_V$ pixels. The inverse of this slope gives us the flux produced per $A_V$-unit of material present, which can be interpreted as an emission coefficient for the dust. Figure 17 shows these inverse slopes as a function of

wavelength. We applied the same technique to our HST data (Figure 16), showing that this relation extends into the visual and UV, albeit less significant and with greater scatter. As a consequence, the error bars for the emission coefficients are larger. This is to be expected, as the observed flux is more complex for these wavelengths. There are contributions from the gas emission, as well as scattered light from $\gamma$ Cas.

These measurements can serve as useful observational constraints for future SED modeling efforts of IC 63 or similar objects, especially those that model the dust composition and emissivity. We therefore present the values and their uncertainties in Table 3. In the next section, we will use these results to fit a simple dust emission model.

### 4.5. Model for the Average Dust SED

The surface brightness of dust with a single equilibrium temperature $T_d$ is given by

$$S_{\nu}(\lambda) = \kappa(\lambda)\Sigma_{d}B_{\nu}(\lambda; T_{d}) \quad (13)$$

where $\kappa$ is the grain absorption cross section per unit dust mass as a function of $\lambda$, $\Sigma_d$ is the dust surface mass density and $B_\nu$ is the Planck function. The subscript $\nu$ indicates that we are using quantities per unit frequency (e.g., MJy). Within this model with $\kappa$ and $\Sigma_d$ as parameters, $\kappa$ and $\Sigma_d$ are degenerate. But the quantity that we have measured, that is displayed in Figure 17 is actually the surface brightness per unit of $A_V$:

$$\frac{S_{\nu}(\lambda)}{A_V} = \kappa(\lambda) \frac{\Sigma_d}{A_V} B_{\nu}(\lambda; T_d) \quad (14)$$

So if a value for the ratio $\Sigma_d/A_V$ is known, it becomes possible to measure $\kappa$. In practice, we use a modified blackbody model
<figure><img src="image_15.png" /><figcaption>Figure 14. Correlation of $A_V$ with the reprojected Herschel fluxes. The black line shows the best-fit slope. The Pearson correlation coefficient $\rho$ is shown, together with the $p$-value for a two sided hypothesis test, where a slope of 0 is the null hypothesis.</figcaption></figure>

of the form

$$ \frac{S_{\nu}(\lambda)}{A_{V}} = \kappa_{\text{eff},160} \frac{\Sigma_{d}}{A_{V}} \left( \frac{\lambda}{160 \, \mu\text{m}} \right)^{-\beta_{\text{eff}}} B_{\nu}(\lambda; T_{\text{eff}}) \quad (15) $$

as in Gordon et al. (2014). The scale factor $\kappa_{\text{eff},160}$ determines the effective value of $\kappa$ at 160 $\mu$m. The second and third parameters are the effective spectral index $\beta_{\text{eff}}$ and the effective temperature $T_{\text{eff}}$. These parameters are called effective because the observed surface brightness per $A_V$ is a mix of measurements along different sightlines, and each sightline is a mix of different dust populations. By using the well known relationship

$$ A_V = 1.086 \tau_V \text{ mag} = 1.086 \kappa_{\text{eff},V} \Sigma_d \text{ mag} \quad (16) $$

we can write the model as

$$ \frac{S_{\nu}(\lambda)}{A_{V}} = \frac{\tau_{160}}{A_{V}} \left( \frac{\lambda}{160 \, \mu\text{m}} \right)^{-\beta} B_{\nu}(\lambda; T) \quad (17) $$

$$ = \frac{\kappa_{160}}{1.086\kappa_V} \left( \frac{\lambda}{160 \, \mu\text{m}} \right)^{-\beta} B_{\nu}(\lambda; T) \text{mag}^{-1} \quad (18) $$

where we drop the “eff” subscript for conciseness. Therefore, by fitting this model to our measurements we can determine the average 160 $\mu$m to V band opacity ratio. We treat the ratio $\tau_{160}/A_V = \kappa_{160}/1.086\kappa_V \text{mag}^{-1}$ as a single fit parameter.

We found that the sum of two modified blackbody models works well to fit the slopes calculated from the Herschel and

Spitzer maps simultaneously. We use a least-squares fitting approach using the standard options of Astropy's modeling module. Three fit results are displayed in Table 4, and the corresponding curves are shown in Figure 17.

The first modified blackbody was fitted to only the PACS and SPIRE points, with a total $\chi^2$ of 0.61. The second model fit included the IRAC points and the MIPS 24 $\mu$m point, but not the IRAC 3.6 $\mu$m point. The latter was excluded because it is systematically higher than the IRAC 4.5 point, preventing the model from properly fitting the other data points. For this model, the $\kappa$-ratio and $\beta$ parameters are not well constrained, and $\chi^2$ in this case is 3.55. The estimated Pearson correlation coefficients between the three parameters are all over 90%, which is a known problem for this parameterization (Gordon et al. 2014).

Fixing $\beta = 2.06$ for both models gives stronger constraints on the other two parameters, especially on the temperatures. Lastly, we simultaneously fitted the sum of the two modified blackbodies with the same fixed value for $\beta$, leading to the gray, dashed curve in Figure 17. For this two-component modified blackbody model, $\chi^2$ is 4.18 when fitting the four free parameters to our nine data points, and the resulting parameters are shown on the third line of Table 4.

Before we discuss these fit results, some remarks should be made about the meaning of these two models and their temperatures. The use of two modified blackbody components does not relate to the existence of a warm and a cool component of the gas in the cloud (Andrews et al. 2018). Instead, they each model a different wavelength range in the
<figure><img src="image_16.png" /><figcaption>Figure 15. Correlation of $A_V$ with the reprojected Spitzer fluxes, analogous to Figure 14.</figcaption></figure>

SED, for which the observed flux is emitted by a different subpopulation of the grains.

The first component and its temperature parameter $T_1$ can be interpreted rather straightforwardly, as the observed flux at the Herschel wavelengths consists of thermal continuum emission from big grains. These grains are in thermal equilibrium with the local radiation field, and $T_1$ can be used as a measure of their average temperature. This temperature is expected to be higher for grains residing in the low density PDR region where the radiation field is stronger, and lower for those in the cold dense parts of the cloud.

On the contrary, the flux in the Spitzer bands originates from VSGs and PAHs. The VSGs are heated stochastically by photons, and the resulting SED shows a complex behavior with several peaks in the MIR (Draine & Li 2001). Despite this, the MIR SED (except the IRAC 3.6 µm point) seems to be well described by a modified blackbody with effective temperature $T_2$.

Using the same PACS 70 and 160 µm images, and assuming a fixed $\beta = 1.8$, Andrews et al. (2018) derive a color dust temperature of 30 K in IC 63. They also perform single-component modified blackbody fits with the same fixed $\beta$ at the tip of the nebula, and again obtain a dust temperature of ~30 K. This matches well with our value for $T_1$. For the gas temperature, they report two values: 207 ± 30 K from H$_2$ excitation diagrams, and 130 K from a PDR model. Note that incidentally, these values for the temperature are of the same order of magnitude as $T_2$. We find it unlikely, however, that these two values would be coupled, as the physics driving the

MIR emission by VSGs and PAHs is very different from the processes driving the H$_2$ excitation (Habart et al. 2011).

As for the optical depth, Andrews et al. (2018) find that $\tau_{160} = 8 \times 10^{-4}$. Using our value of $A_{V,tip} = 1.41$ mag yields $\tau_{160}/A_V = 5.7 \times 10^{-4}$ mag$^{-1}$, which is compatible with our value for $\tau_{160;1}/A_V$. Multiplying the last two values for $\tau_{160}/A_V$ in Table 4 with 1.086 (Equation (16)) gives us the result

$$ \frac{\kappa_{160;1}}{\kappa_V} = (6.98 \pm 0.88) \times 10^{-4} \qquad (19) $$

$$ \frac{\kappa_{160;2}}{\kappa_V} = (2.91 \pm 0.53) \times 10^{-9}. \qquad (20) $$

In the above equations $\kappa_V$ is the combined V-band cross section of the two components. Since the first component dominates, the value of $\kappa_{160;1}/\kappa_V$ will depend on the optical properties of the big grains. The second component makes only a small contribution to the total value, and therefore $\kappa_{160;2}/\kappa_V$ will also depend on the relative mass density of the VSGs and PAHs. We can compare the value of $\kappa_{160;1}/\kappa_V$ to theoretical dust models, using $\kappa_{0.54}$ (the cross section at 540 nm) as an approximation for $\kappa_V$. With their standard grain size distributions, we find $\kappa_{160}/\kappa_{0.54} = 1.4 \times 10^{-3}$ for the model of Draine & Li (2007), and $\kappa_{160}/\kappa_{0.54} = 7.4 \times 10^{-4}$ for the THEMIS model (Jones et al. 2017).

## 4.6. Per-pixel Modified Blackbody Fits

Instead of making a fit to the average $\Delta S_i(\lambda)/\Delta A_V$ SED obtained from the slopes, we also perform individual fits for
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Analogous to Figures 14 and 15, for the HST data. Some of the correlations are somewhat less significant, as can be seen from the *ρ* and *p* values, but we still obtain usable values for the slopes, albeit with larger error bars.</figcaption></figure>

each pixel of the maps, based on the same reprojected *Spitzer* and *Herschel* images. In this case, the scale parameter for the modified blackbodies is τ₁₆₀ instead of τ₁₆₀/A<sub>V</sub>, since we are fitting directly to the observed fluxes, instead of average flux versus A<sub>V</sub> slopes. We are not fitting a differential measurement this time, and unlike the previous section we need to determine and subtract the background in each band. Since we have a measurement of the A<sub>V</sub> background of 2.51 mag, we use the slope and intercept of the lines in Figures 15 and 14 to calculate an appropriate level for the background in each band. By doing these per-pixel fits, we obtain maps of the same resolution for each parameter, shown in Figure 18. We do not show the parameters for the MIR model, as T₂ and β₂ are quite noisy,

and the τ₁₆₀;₂ map looks very similar to that of τ₁₆₀;₁, besides being of a different order of magnitude.

Figure 18 shows a gradient in the temperature, perpendicular to the direction of the incoming radiation field. The typical uncertainty on T₁ is about 2–5 K. A feature of the same shape exists for the β₁ and τ₁₆₀;₁ parameters too, which shows an anticorrelation with the temperature. The β–T relationship is well known and can be explained by variations in the dust properties (Dupac et al. 2003; Ysard et al. 2015). However, because of noise and temperature variations along the line of sight, least-squares fitting naturally leads to this type of correlation (Shetty et al. 2009a, 2009b; Juvela & Ysard 2012). Therefore, we performed the same fits with β fixed to 2.06. A very notable difference is that the
<table><caption>**Table 3**
Measurements and Uncertainties of the Average Surface Brightness per $A_V$ Unit</caption><thead><tr><th>Band</th><th>ΔS<sub>ν</sub>(λ)/ΔA<sub>V</sub><br>(MJy sr<sup>-1</sup> mag<sup>-1</sup>)</th><th>σ<br>(MJy sr<sup>-1</sup> mag<sup>-1</sup>)</th></tr></thead><tbody><tr><td>F275W</td><td>0.183</td><td>0.025</td></tr><tr><td>F336W</td><td>0.250</td><td>0.035</td></tr><tr><td>F475W</td><td>0.669</td><td>0.350</td></tr><tr><td>F625W</td><td>0.520</td><td>0.094</td></tr><tr><td>F814W</td><td>0.408</td><td>0.069</td></tr><tr><td>F110W</td><td>1.007</td><td>0.239</td></tr><tr><td>F160W</td><td>1.115</td><td>0.242</td></tr><tr><td>IRAC 3.6</td><td>2.472</td><td>0.332</td></tr><tr><td>IRAC 4.5</td><td>1.216</td><td>0.164</td></tr><tr><td>IRAC 5.8</td><td>11.770</td><td>1.563</td></tr><tr><td>IRAC 8.0</td><td>33.961</td><td>4.387</td></tr><tr><td>MIPS 24</td><td>26.741</td><td>3.496</td></tr><tr><td>PACS 70</td><td>399.678</td><td>54.112</td></tr><tr><td>PACS 160</td><td>345.684</td><td>45.670</td></tr><tr><td>SPIRE 250</td><td>103.251</td><td>12.661</td></tr><tr><td>SPIRE 350</td><td>38.976</td><td>4.734</td></tr><tr><td>SPIRE 500</td><td>13.184</td><td>1.731</td></tr></tbody></table>

trend in *T* becomes inverted, and that the temperature range becomes narrower. Because the $β-T$ trend appears partly due to noise-induced correlations and partly due to a physical relationship, these results are not easily interpreted. On the one hand, fixing $β$ might elminate any exaggerated variations of *T* that occur due to the correlations in the fit. On the other hand, choosing a constant $β$ value introduces a bias in *T*, which will depend on the real value of $β$. A more in-depth analysis (e.g., Galliano 2018) might be able to disentangle these degeneracies, but would go beyond the scope of this work.

Comparing the $τ_{160:1}$ map to the $A_V$ map in Figure 11, we find a similar structure. By dividing the $τ_{160:1}$ values by the $A_V$ value in each pixel (minus the background of 2.51 mag), and using $A_V = 1.086τ_V$, we obtain maps of $κ_{160:1}/κ_V$. We see that this ratio is quite noisy, and comparing the top and bottom rows of Figure 18, fixing $β$ seems to affect the pixels in the bottom left corner more than others. The general trend is that $κ_{160:1}/κ_V$ seems to be higher toward the lower edge of the image. So while the $R_V$ map points to an increase in average grain size at the front of the nebula, this analysis of the Herschel data and the $A_V$ map shows a less pronounced trend which more or less aligns with its south edge.

We did not find any significant spatial correlation between $A_V$ or $R_V$ and any of the fitted parameters, except $A_V$ and $τ$. This holds regardless of whether a free or fixed $β$ was used. For some of the physical conditions, Fleming et al. (2010) provide maps based on Spitzer IRS data: total PAH emission and band ratios; H$_2$ emission, column density, temperature, and ortho-to-para ratio; and the radiation field $G_0/n_H$. These maps cover a 40 by 40 arcsec area at the very tip of the nebula, which corresponds to about 4 pixels of our $R_V$ map. Figure 7 of Fleming et al. (2010) shows an arc shaped feature in the PAH emission, where the H$_2$ column density is also the highest. This feature might be related to the higher $R_V$ area observed at the front of the nebula, but the resolution of our $R_V$ measurements is insufficient to confirm this.

## 4.7. Comparison with Other Clouds

In the Orion bar, the dust temperature shows a much stronger gradient, going from 71 K in front of the bar, to 49 K inside,

and 37 K behind the bar (Arab et al. 2012). The distance between the first and last point for which these temperatures were measured, is of the order 0.1 pc. The size of the area we study in IC 63 is ~2' or also ~0.1 pc, assuming a distance of 190 pc. Over the same spatial scale, there is at best a subtle temperature gradient in IC 63. Note that the radiation field and the density are more than an order of magnitude higher in the Orion bar: $G_0 = 4 \times 10^4$ and $n_H = 2 \times 10^5$ cm$^{-3}$. Note that the conditions in IC 63 are $G_0 = 1100$ (150 in Andrews et al. 2018), and $n_H = 1.2 \times 10^4$ cm$^{-3}$. Arab et al. (2012) measure the values $β = 1.2, 1.6$, and 2.2 for the same three points, respectively. These differences likely originate from dust evolution in the Orion bar. For IC 63 we find it hard to determine the significance of the variations in the fitted $β$, and we observe no correlations between the changes in $β$ and the observed differences in $R_V$.

In Salgado et al. (2016), both the bar and the H II region of the Orion nebula were studied using MIR observations with the FORCAST instrument on board SOFIA, and Herschel/PACS maps. A two-component modified blackbody model was necessary to properly fit the thermal dust emission. The cool component (~40 K) was associated with the molecular cloud, and the warm component (~80 K) with the dust in the PDR/H II region. For IC 63, only one component is needed to fit the thermal dust emission; we added a second component because it was able to describe the MIR emission of the PAHs and VGSs rather well.

For NGC 7023, Köhler et al. (2014) performed similar, per-pixel modified blackbody fits (with $β$ fixed to 2). The dust temperature ranges from 20.1 to 30.3 K in NGC 7023 NW, and from 17.2 to 20.5 K in NGC 7023 E. The temperature differences are smaller than in the Orion bar, and closer to what we observe in IC 63. Through the $τ$ factor of their model (in this case at 250 µm), they derive line-of-sight column densities ranging from 2.8 to $9.8 \times 10^{22}$ cm$^{-2}$, more than an order of

<figure><img src="image_18.png" /><figcaption>**Figure 17.** Surface brightness per $A_V$ unit for each of the MIR and FIR observations provided by Spitzer and Herschel, and modified blackbody fits to these data. The points for the HST data are added for completeness, but are not used in the blackbody fits. The data points are calculated by fitting slopes to the scatter plots in Figures 14–16. Red: fit to the IRAC and MIPS 24 points, ignoring the IRAC 3.6 point. Blue: fit to the PACS and SPIRE points. Gray dashed: combined model, fitting both sets of points.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Results of the per-pixel modified blackbody fitting with free $β_1$ (top row) and fixed $β_1$ = 2.06 (bottom row). Red X: bad value. This is mainly due to missing $A_V$ or $A_V$ < 2.51 (the assumed background). The top and bottom panels for each parameter use the exact same color scale.</figcaption></figure>

<table><caption>**Table 4**

Overview of the Modified Blackbody Fits to Our Measurement of $\Delta S_i(\lambda)/\Delta A_V$</caption><thead><tr><th>Fitting method</th><th>τ<sub>160;1</sub>/A<sub>V</sub><br>(10<sup>-4</sup> mag<sup>-1</sup>)</th><th>β<sub>1</sub></th><th>T<sub>1</sub><br>(K)</th><th>τ<sub>160;2</sub>/A<sub>V</sub><br>(10<sup>-9</sup> mag<sup>-1</sup>)</th><th>β<sub>2</sub></th><th>T<sub>2</sub><br>(K)</th></tr></thead><tbody><tr><td>Individual,<sup>a</sup> free β</td><td>6.4 ± 1.9</td><td>2.06 ± 0.22</td><td>29.9 ± 2.2</td><td>1.4 ± 2.4</td><td>2.4 ± 1.0</td><td>227 ± 6</td></tr><tr><td>Individual, fixed β</td><td>6.43 ± 0.40</td><td>2.06</td><td>29.9 ± 0.5</td><td>2.70 ± 0.71</td><td>2.06</td><td>227 ± 6</td></tr><tr><td>Simultaneous,<sup>b</sup> fixed β</td><td>6.43 ± 0.81</td><td>2.06</td><td>29.9 ± 1.1</td><td>2.68 ± 0.49</td><td>2.06</td><td>227 ± 4.2</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Model 1 is fit to the *Herschel* points, model 2 to the *Spitzer* points.

<sup>b</sup> The two models and sets of points are combined.

magnitude larger than in IC 63. This means we can expect $A_V$ values along the line of sight of several tens of magnitudes. For the radiation field, Pilleri et al. (2012) provide the values $G_0 = 2600$ for NGC 7023 NW, and $G_0 = 250$ for NGC 7023 E. The density $n_H$ of these two PDRs is about $2 \times 10^4$ to $3 \times 10^4$ cm$^{-3}$. Despite the differences in radiation field and optical depth, the dust temperatures in NGC 7023 NW, NGC 7023 E, and IC 63 do not show as much variation as those in the Orion bar. Concluding, the observed dust variations are likely related to the density $n_H$.

Unfortunately, we cannot compare our $R_V$ measurements to other PDRs, because direct extinction measurements in other PDRs are rarely possible, especially in a spatially resolved way. Reiterating, the number of background stars and the low optical depth make IC 63 one of the only PDRs for which this is possible, at least to our knowledge.

**5. Conclusion**

We have expanded the BEAST code to better support the fitting of background stars of nebulae and other nearby, transparent objects. We have performed individual SED fits of a combined stellar and extinction model to the numerous background stars of the IC 63 PDR, as observed by *HST* in seven broadband filters from the UV to the NIR.

From these fit results, we were able to generate maps of the extinction parameters $A_V$ and $R_V$ across the nebula. We find that

$A_V$ varies between 0.5 and 1.4 mag when we assume a background of 2.5 mag. There is a decreasing trend in $R_V$, from 3.7 at the tip of IC 63 to ~3.4 when moving further away from $\gamma$ Cas. With a correction for background and/or foreground material, the latter value can go as low as 2.5.

It is the first time that this type of measurement has been made for a PDR. This approach was possible for IC 63 because of its low optical depth, making this PDR transparent for a sufficient number of detectable background sources. This same technique would therefore also be applicable to IC 59. Most other PDRs, such as the Horsehead nebula, the Orion bar, and NGC 7023, have much higher optical depths and would not allow a sufficient amount of background stars to be measured.

By combining the $A_V$ map with FIR maps from *Spitzer* and *Herschel*, we were able to measure the shape of the average, $A_V$-normalized surface brightness of the dust. A dual modified blackbody model fits the shape of this curve well, and provides two temperatures (30 and 227 K) and measurements of $\kappa_{160}/\kappa_V$ ((6.98 ± 0.88) × 10$^{-4}$ and (2.91 ± 0.53) × 10$^{-9}$). By performing fits of the same model to individual pixels of the $A_V$ map, we derived a map of $\kappa_{160}/\kappa_V$. The values of this map are of the same order of magnitude, but there are local variations that differ significantly (a factor of 2 or 3) from the average. Both the $R_V$ map and the $\kappa_{160}/\kappa_V$ show variations, but no correlations were found between the two maps, and it is unclear if there are relationships with other physical quantities.
We thank the referee for a very thorough report, leading to significant improvements to the quality and completeness of this work. This research was made possible thanks to a BOF (bijzonder onderzoeksfonds) fellowship, a special research fund provided by Ghent University.

*Software:* DOLPHOT (Dolphin 2000, 2016). Tiny Tim (Krist et al. 2011). BEAST (Gordon et al. 2016). This research made use of Astropy,<fnref n="9" /> a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013, 2018).

# Appendix A
## New BEAST Feature: Parallelization Using Subgrids

By introducing the stellar distance as an extra parameter, the dimensionality of the whole SED model grid is increased from 6D to 7D. Additionally the range of distances is potentially very large (see Table 2), and because the distance rescales the whole SED, a granularity on par with the $A_V$ grid is necessary. Since we chose to use ~40 (log spaced) distance bins, the amount of memory needed to work with the grid is also increased 40-fold. Originally we would split the catalog into many pieces, and then let every process work on a small part of the catalog. This approach is problematic because for every process launched, the whole grid has to be loaded into memory. Therefore, the total memory of the machine still imposes a limit on the number of fitting processes that can be run simultaneously.

To deal with these large grids, a memory-saving parallelization strategy was developed, which we call the *subgridding* approach. This approach entails that each process will work on only a small portion of the grid or *subgrid*. This works for the physics model setup, the noise model generation, as well as the fitting step. For each subgrid, a separate set of statistics and 1D posterior distributions is obtained. By properly taking into account the weight (integrated posterior) of each subgrid and combining these individual 1D posterior distributions, the statistics can be recalculated to reflect those of the whole grid.

To facilitate this parallelization scheme, a set of functions for splitting and/or merging the model grids and output statistics was developed, and an example workflow is given in the documentation.<fnref n="10" /> If the size of the RAM proves problematic for running the BEAST, then the number of subgrids can simply be increased. Arguably, this makes the BEAST scalable to much larger grid sizes.

# Appendix B
## New BEAST Feature: Spatial Variation of Noise Model

Since the extended emission varies for different locations in the nebula, the noise model should ideally depend on the position of the source that is being fit. To investigate the spatial dependence of the noise model, a set of tools was developed and added to the BEAST. The concept is that a number of regions on the sky that have similar intensities of background/ foreground emission are determined (background regions). Within each of these background regions, a separate set of ASTs is performed, and an individual noise model is created. If the resulting noise models differ significantly, an individual fitting run can be done for each background region, using the corresponding noise model to fit only the sources that fall within that region.

### B.1. Background Map Tool

The first tool consists of a stand-alone script which can generate a map, with pixels of a user-specified size along the R. A. and decl. axes. To measure the background for a single pixel, it is first determined which of the sources of the input catalog fall into this bin, according to their positions. Then, a reference image is used (the F625W image in our case) to measure the background for each source. The sources themselves are masked out, and their background is determined using an annulus of a certain radius and thickness. Determining an ideal inner and outer radius for the annuli may require some trial and error by the user. The background for this pixel is then set to be the median of these individual background measurements. For IC 63, the resulting background map is shown in the left panel of Figure 19.

As a side note, this same tool can also be used to construct a map of point-source densities in the exact same format. This is useful for catalogs that have a very high but variable density of stars, leading to different crowding effects depending on the position. The source density is derived directly from the positions listed in a given catalog. Since the rest of the code is agnostic of the quantity described by either map, a common approach can be used to model either effects of extended emission or those of crowding.

### B.2. Selection of AST Positions

One of the applications that uses this background map is the routine for selecting the positions of the ASTs. We want to make sure that the areas corresponding to different levels of background intensity are all sufficiently sampled. We achieve this by first choosing a set of bins (usually a small number, due to the available hard disk space and computing time), each of which represents a certain range in background intensity. Then, the pixels of the map described above can be distributed among these bins, leading to groups of pixels which constitute the desired background regions (Figure 19). A fixed set of SEDs is reused for each region, and each SED is assigned a random position within that region. This way, each background region has the same number of samples regardless of its total surface area.

### B.3. Creating Individual Noise Models

Once the selected ASTs have been run through our photometry routine, and a catalog of fake stars has been produced, both the observed and the fake catalog are split up according to the background regions (see again Figure 19, right panel). By applying the same spatial split to both the observed catalog and the fake catalog, a series of catalog pairs is created. The BEAST can then use each pair consisting of one fake and one observed catalog to generate the individual noise models, and perform the individual fitting runs if the resulting noise models are found to be significantly different.

Using the regions shown in the right panel of Figure 19, we created four noise models (Figure 20). Despite the differences between the noise parameters for the background bins, we found that the predicted $A_V$ and $R_V$ are not significantly influenced by using these extra noise models. In each bin there are at most four stars that have an $A_V$ or $R_V$ that differs more than 10%. For some of the other parameters the differences can be larger, but the fit quality is not necessarily better. Therefore, the results presented in Section 3 and the

<footnote n="9">⁹ http://www.astropy.org</footnote>

<footnote n="10">¹⁰ https://beast.readthedocs.io/en/latest/subgrid_parallelism.html</footnote>
<figure><img src="image_20.png" /><figcaption>**Figure 19.** Left: illustration of the background map tool. In this case, the F625W image was used as an input, so the background intensity mainly consists of Hα emission. The point sources have been masked using the positions listed in the cleaned up catalog described in Section 2.3.2. This can be seen from the empty disks on the image. The resolution of the map is chosen so that there's at least a handful of sources in each pixel. Right: example of a set of spatial regions, derived from the background map on the left. The colors indicate which background intensity regime each spatial bin and, by extension, each source (circles) belongs to.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 20.** Four different noise models, one for each background region. The numbers 1, 2, 3, and 4 point to the regions shown in the right panel of Figure 19.</figcaption></figure>

rest of the analysis are based on fits using the single noise model that was shown in Figure 6. Note that this noise model still benefits from the background-based AST position picking, because it makes sure that the tip of the nebula is sufficiently sampled.

## Appendix C
### The $A_V$ and $R_V$ Maps in Table Format

We provide the values and standard deviation for $A_V$ and $R_V$ of the maps in Tables 5 and 6. The R.A. and decl. values in the table are the coordinates of the bottom left corner of each pixel,
<table><caption>**Table 5**
Av Map and Uncertainty (Same Orientation as Figure 11)</caption><thead><tr><th>Decl./R.A. (deg)</th><th>14.8021</th><th>14.7878</th><th>14.7736</th><th>14.7593</th><th>14.7450</th><th>14.7308</th></tr></thead><tbody><tr><td>60.9206</td><td>3.54 ± ?</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td></tr><tr><td>60.9137</td><td>3.32 ± 0.48</td><td>3.00 ± 0.22</td><td>3.00 ± 0.33</td><td>3.26 ± 0.36</td><td>4.75 ± 0.36</td><td>X</td></tr><tr><td>60.9067</td><td>3.92 ± 0.29</td><td>3.10 ± 0.28</td><td>3.03 ± 0.25</td><td>3.01 ± 0.16</td><td>5.12 ± 0.45</td><td>X</td></tr><tr><td>60.8998</td><td>3.21 ± 0.24</td><td>3.48 ± 0.24</td><td>3.60 ± 0.22</td><td>3.19 ± 0.19</td><td>2.94 ± 0.21</td><td>3.15 ± 0.84</td></tr><tr><td>60.8928</td><td>3.04 ± 0.30</td><td>3.28 ± 0.15</td><td>3.06 ± 0.19</td><td>3.92 ± 0.23</td><td>3.04 ± 0.23</td><td>X</td></tr><tr><td>60.8859</td><td>2.98 ± 0.22</td><td>3.04 ± 0.16</td><td>3.30 ± 0.16</td><td>3.44 ± 0.36</td><td>3.56 ± 0.36</td><td>X</td></tr><tr><td>60.8790</td><td>2.67 ± 0.12</td><td>2.78 ± 0.21</td><td>2.38 ± 0.18</td><td>2.29 ± 0.14</td><td>2.49 ± 0.15</td><td>X</td></tr></tbody></table>

<table><caption>**Table 6**
Rv Map and Uncertainty (Same Orientation as Figure 11)</caption><thead><tr><th>Decl./R.A. (deg)</th><th>14.8021</th><th>14.7878</th><th>14.7736</th><th>14.7593</th><th>14.7450</th><th>14.7308</th></tr></thead><tbody><tr><td>60.9206</td><td>2.76 ± ?</td><td>X</td><td>X</td><td>X</td><td>X</td><td>X</td></tr><tr><td>60.9137</td><td>3.44 ± 0.11</td><td>3.73 ± 0.14</td><td>3.48 ± 0.10</td><td>3.47 ± 0.15</td><td>3.75 ± 0.36</td><td>X</td></tr><tr><td>60.9067</td><td>3.51 ± 0.16</td><td>3.34 ± 0.10</td><td>3.34 ± 0.10</td><td>3.68 ± 0.11</td><td>3.68 ± 0.13</td><td>X</td></tr><tr><td>60.8998</td><td>3.43 ± 0.14</td><td>3.52 ± 0.11</td><td>3.45 ± 0.14</td><td>3.55 ± 0.07</td><td>3.59 ± 0.11</td><td>3.60 ± 0.25</td></tr><tr><td>60.8928</td><td>3.50 ± 0.31</td><td>3.58 ± 0.10</td><td>3.40 ± 0.09</td><td>3.58 ± 0.10</td><td>3.65 ± 0.15</td><td>X</td></tr><tr><td>60.8859</td><td>3.31 ± 0.09</td><td>3.53 ± 0.16</td><td>3.72 ± 0.14</td><td>3.63 ± 0.15</td><td>3.66 ± 0.26</td><td>X</td></tr><tr><td>60.8790</td><td>3.90 ± 0.06</td><td>3.61 ± 0.17</td><td>3.49 ± 0.15</td><td>3.61 ± 0.09</td><td>3.57 ± 0.11</td><td>X</td></tr></tbody></table>

with the tables oriented in the same way as the color maps in Figure 11.

## ORCID iDs

Dries Van De Putte https://orcid.org/0000-0002-5895-8268

Karl D. Gordon https://orcid.org/0000-0001-5340-6774

Julia Roman-Duval https://orcid.org/0000-0001-6326-7069

Benjamin F. Williams https://orcid.org/0000-0002-7502-0597

Maarten Baes https://orcid.org/0000-0002-3930-2757

Kirill Tchernyshyov https://orcid.org/0000-0003-0789-9939

## References

Abergel, A., Arab, H., Compiègne, M., et al. 2010, A&A, 518, L96

Anderson, J., & King, I. R. 2000, PASP, 112, 1360

Andersson, B.-G., Pirola, V., De Buizer, J., et al. 2013, ApJ, 775, 84

Andrews, H., Peeters, E., Tielens, A. G. G. M., & Okada, Y. 2018, A&A, 619, A170

Arab, H. 2015, HST Proposal, 23, 14186

Arab, H., Abergel, A., Habart, E., et al. 2012, A&A, 541, A19

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58

Bakes, E. L. O., & Tielens, A. G. G. M. 1994, ApJ, 427, 822

Barentsen, G., Farnhill, H. J., Drew, J. E., et al. 2014, MNRAS, 444, 3230

Berné, O., Joblin, C., Deville, Y., et al. 2007, A&A, 469, 575

Boersma, C., Bregman, J., & Allamandola, L. J. 2014, ApJ, 795, 110

Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978, ApJ, 224, 132

Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127

Bressan, A., Marigo, P., Girardi, L., Nanni, A., & Rubele, S. 2013, European Physical Journal Web of Conferences, 43, 03001

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245

Castelli, F., & Kurucz, R. L. 2003, arXiv:astro-ph/0405087

Compiègne, M., Abergel, A., Verstraete, L., & Habart, E. 2008, A&A, 491, 797

Compiègne, M., Verstraete, L., Jones, A., et al. 2011, A&A, 525, A103

Dalcanton, J. J., Williams, B. F., Lang, D., et al. 2012, ApJS, 200, 18

Desert, F.-X., Boulanger, F., & Puget, J. L. 1990, A&A, 237, 215

Diplas, A., & Savage, B. D. 1994, ApJ, 427, 274

Dolphin, A. 2016, DOLPHOT: Stellar photometry, Astrophysics Source Code Library, ascl:1608.013

Dolphin, A. E. 2000, PASP, 112, 1383

Draine, B. T. 1978, ApJS, 36, 595

Draine, B. T. 2003, ARA&A, 41, 241

Draine, B. T., & Li, A. 2001, ApJ, 551, 807

Draine, B. T., & Li, A. 2007, ApJ, 657, 810

Drew, J. E., Greimel, R., Irwin, M. J., et al. 2005, MNRAS, 362, 753

Dupac, X., Bernard, J.-P., Boudet, N., et al. 2003, A&A, 404, L11

Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10

Fitzpatrick, E. L. 1999, PASP, 111, 63

Fleming, B., France, K., Lupu, R. E., & McCandliss, S. R. 2010, ApJ, 725, 159

Fruchter, A. S., & Hook, R. N. 2002, PASP, 114, 144

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1

Galliano, F. 2018, MNRAS, 476, 1445

Girardi, L., Williams, B. F., Gilbert, K. M., et al. 2010, ApJ, 724, 1030

Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, ApJ, 594, 279

Gordon, K. D., Fouesneau, M., Arab, H., et al. 2016, ApJ, 826, 104

Gordon, K. D., Roman-Duval, J., Bot, C., et al. 2014, ApJ, 797, 85

Gorenstein, P. 1975, ApJ, 198, 95

Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3

Güver, T., & Özel, F. 2009, MNRAS, 400, 2050

Habart, E., Abergel, A., Boulanger, F., et al. 2011, A&A, 527, A122

Habing, H. J. 1968, BAN, 19, 421

H14PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116

Hollenbach, D., & Salpeter, E. E. 1971, ApJ, 163, 155

Hollenbach, D. J., & Tielens, A. G. G. M. 1997, ARA&A, 35, 179

Hollenbach, D. J., & Tielens, A. G. G. M. 1999, RvMP, 71, 173

Jansen, D. J., van Dishoeck, E. F., & Black, J. H. 1994, A&A, 282, 605

Jansen, D. J., van Dishoeck, E. F., Black, J. H., Spaans, M., & Sosin, C. 1995, A&A, 302, 223

Jansen, D. J., van Dishoeck, E. F., Keene, J., Boreiko, R. T., & Betz, A. L. 1996, A&A, 309, 899

Jones, A. P., Köhler, M., Ysard, N., Bocchio, M., & Verstraete, L. 2017, A&A, 602, A46

Juvela, M., & Ysard, N. 2012, A&A, 541, A33

Köhler, M., Habart, E., Arab, H., et al. 2014, A&A, 569, A109

Köhler, M., Ysard, N., & Jones, A. P. 2015, A&A, 579, A15

Krist, J. E., Hook, R. N., & Stoehr, F. 2011, Proc. SPIE, 8127, 81270J

Lanz, T., & Hubeny, I. 2003, ApJS, 146, 417
Lanz, T., & Hubeny, I. 2007, ApJS, 169, 83

Luri, X., Brown, A. G. A., Sarro, L. M., et al. 2018, A&A, 616, A9

Marigo, P., Girardi, L., Bressan, A., et al. 2008, A&A, 482, 883

Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 500, 501

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1

Pilleri, P., Montillaud, J., Berné, O., & Joblin, C. 2012, A&A, 542, A69

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

Predehl, P., & Schmitt, J. H. M. M. 1995, A&A, 293, 889

Reina, C., & Tarenghi, M. 1973, A&A, 26, 257

Rieke, G. H., Young, E. T., Engelbracht, C. W., et al. 2004, ApJS, 154, 25

Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523

Roman-Duval, J., Jenkins, E. B., Williams, B., et al. 2019, ApJ, 871, 151

Salgado, F., Berné, O., Adams, J. D., et al. 2016, ApJ, 830, 118

Schlafly, E. F., Peek, J. E. G., Finkbeiner, D. P., & Green, G. M. 2017, ApJ, 838, 36

Shetty, R., Kauffmann, J., Schnee, S., & Goodman, A. A. 2009a, ApJ, 696, 676

Shetty, R., Kauffmann, J., Schnee, S., Goodman, A. A., & Ercolano, B. 2009b, ApJ, 696, 2234

Steinacker, J., Baes, M., & Gordon, K. D. 2013, ARA&A, 51, 63

Stepnik, B., Abergel, A., Bernard, J.-P., et al. 2003, A&A, 398, 551

Tange, O. 2011, GNU Parallel 20150322, Zenodo, doi:10.5281/zenodo.16303

Thi, W. F., van Dishoeck, E. F., Bell, T., Viti, S., & Black, J. 2009, MNRAS, 400, 622

Valencic, L. A., Clayton, G. C., Gordon, K. D., & Smith, T. L. 2003, ApJ, 598, 369

Wakelam, V., Bron, E., Cazaux, S., et al. 2017, MolAs, 9, 1

Weingartner, J. C., & Draine, B. T. 2001, ApJS, 134, 263

Williams, B. F., Lang, D., Dalcanton, J. J., et al. 2014, ApJS, 215, 9

Yanchulova Merica-Jones, P., Sandstrom, K. M., Johnson, L. C., et al. 2017, ApJ, 847, 102

Ysard, N., Köhler, M., Jones, A., et al. 2015, A&A, 577, A110

Zhu, H., Tian, W., Li, A., & Zhang, M. 2017, MNRAS, 471, 3494

Zubko, V., Dwek, E., & Arendt, R. G. 2004, ApJS, 152, 211