Dust evolution across the Horsehead Nebula

Micro-physical processes on interstellar dust surfaces are tightly connected to dust properties (i.e. dust composition, size and shape) and play a key role in numerous phenomena in the interstellar medium (ISM). The large disparity in physical conditions (i.e. density, gas temperature) in the ISM triggers an evolution of dust properties. The analysis of how dust evolves with the physical conditions is a stepping-stone towards a more thorough understanding of interstellar dust. The aim of this paper is to highlight dust evolution in the Horsehead Nebula PDR region. We use Spitzer/IRAC (3.6, 4.5, 5.8 and 8 {\mu}m), Spitzer/MIPS (24 {\mu}m) together with Herschel/PACS (70 and 160 {\mu}m) and Herschel/SPIRE (250, 350 and 500 {\mu}m) to map the spatial distribution of dust in the Horsehead over the entire emission spectral range. We model dust emission and scattering using the THEMIS interstellar dust model together with the 3D radiative transfer code SOC. We find that the nano-grains dust-to-gas ratio in the irradiated outer part of the Horsehead is 6 to 10 times lower than in the diffuse ISM. Their minimum size is 2 to 2.25 times larger than in the diffuse ISM and the power-law exponent of their size distribution, 1.1 to 1.4 times lower than in the diffuse ISM. Regarding the denser part of the Horsehead, it is necessary to use evolved grains (i.e. aggregates, with or without an ice mantle). It is not possible to explain the observations using grains from the diffuse medium. We therefore propose the following scenario to explain our results. In the outer part of the Horsehead, all the nano-grains have not yet had time to re-form completely through photo-fragmentation of aggregates and the smallest of the nano-grains that are sensitive to the radiation field are photo-destroyed. In the inner part of the Horsehead, grains most likely consist of multi-compositional, mantled aggregates.


Introduction
Interstellar dust plays an essential role within the interstellar medium (ISM) through different microphysical processes happening on dust surfaces that can heat the gas, such as the photoelectric effect (e.g. Bakes & Tielens 1994;Weingartner & Draine 2001b), or cool the gas through gas-grain collisions (Burke & Hollenbach 1983). By acting as a catalyst, allowing atoms and molecules to react on its surface, dust is strongly involved in the chemistry of the ISM (e.g. Hollenbach & Salpeter 1971;Bron et al. 2014;Jones & Habart 2015;Wakelam et al. 2017). Also, dust plays a role in the redistribution of UV-visible stellar radiation into IR-mm radiation, a process that depends on the dust mass and the volume of dust grains (e.g. Draine 2003; Compiègne et al. 2011). The efficiency of these processes strongly depends on the dust properties, such as the grain size, composition and shape. It is therefore important to constrain dust properties in order to understand the different phenomena that take place in the ISM. To this purpose, several dust models have been developed and are classified in three categories. Those composed of silicate and graphite (e.g Mathis et al. 1977;Draine & Lee 1984;Kim et al. 1994) with an extension of these models, using PAHs (polycyclic aromatic hydrocarbons) (e.g. Siebenmorgen & Kruegel 1992;Li & Draine 2001;Weingartner & Draine 2001a). As a result of fragmentation and coagulation processes in the ISM, dust models with grains that have a core-mantle structure (e.g. Desert et al. 1990;Jones et al. 1990;Li & Greenberg 1997) and dust composite models composed of silicate and carbon grain aggregates (e.g. Mathis & Whiffen 1989;Zubko et al. 2004) have been proposed. In this paper, we use the THEMIS dust model (see Jones et al. 2013;Köhler et al. 2014;Jones et al. 2014;Köhler et al. 2015;Ysard et al. 2015;Jones et al. 2017), developed in combination with the results of laboratory experiments and astronomical observations. The cornerstone of this model is its self-consistent view of the evolution of the dust constituents through the ISM. This view is required for understanding dust evolution in response to the local ISM conditions (i.e. density, radiation field).
Some of the first evidence of dust evolution was shown by Fitzpatrick & Massa (1986) through the variation in the 2175 Å interstellar bump from diffuse (R V = 3.1) to denser regions (up to R V ∼ 5.5). Similarly, other studies (e.g. Cardelli et al. 1989;Cardelli & Clayton 1991;Campeggio et al. 2007) found the same variations that have for the first time been explained by Kim et al. (1994) by stating that these observations are consistent with a decrease in the carbonaceous nano-grains abundance (relative to the gas) together with an increase in larger grain abundance. It is also possible to follow dust evolution from its emission in the mid-IR (due to stochastically heated nano-grains) and in the far-IR (where large grains at thermal equilibrium emit). This has led to a wealth of studies (e.g. Boulanger et al. 1990;Laureijs et al. 1991;Abergel et al. 1994;Bernard et al. 1999;Stepnik et al. 2003;Flagey et al. 2009) revealing that nano-grains disappear in dense regions as they coagulate onto larger grains. Dust evolu-Article number, page 1 of 17 arXiv:2003.05902v1 [astro-ph.GA] 12 Mar 2020 A&A proofs: manuscript no. aa tion is also highlighted by variation of its far-IR opacity with the local environment (e.g. Juvela et al. 2011;Planck Collaboration et al. 2011;Martin et al. 2012;Roy et al. 2013;Ysard et al. 2013;Köhler et al. 2015;Juvela et al. 2015) explained with dust coagulation and accretion of ice mantles, a scenario which is supported by numerical simulations of dust evolution in dense regions (e.g. Ossenkopf & Henning 1994;Ormel et al. 2011;Köhler et al. 2015).
Photon-dominated regions (PDRs) (Hollenbach & Tielens 1997, 1999 correspond to the interface between HII regions and molecular clouds that are irradiated by energetic stars close by. In these regions, the physical conditions are strongly constrated hence PDRs are a unique place to study how do dust, gas and local physical conditions evolve with depth. Based on dust emission variations in the mid-IR observed with Spitzer in several PDRs (Ced 201, NGC 7023, ρ Ophiuchi West filament), Berné et al. (2007) construed that such variations can be explained by the photo-processing of carbonaceous nano-grains, a scenario later reinforced in other PDRs (Abergel et al. 2010;Pilleri et al. 2012;Boersma et al. 2014;Pilleri et al. 2015). Using far-IR observations from Herschel, together with the near and mid-IR observations from Spitzer, Arab et al. (2012) found that the carbonaceous nano-grains abundance decreases together with an increase in the opacity of the large grains in the Orion bar. They claimed that these variations are likely due to coagulation processes in the denser part of this region. Evidence of dust evolution has also been shown in IC 63 based on extinction mapping (Van De Putte et al. 2019).
In this paper, we focus on a well-known PDR, the Horsehead, that has previously been studied from the perspective of dust observations (e.g. Abergel et al. 2003;Teyssier et al. 2004;Compiègne et al. 2007;Pety et al. 2005;Compiègne et al. 2008;Arab 2012), gas observations (e.g. Habart et al. 2005;Goicoechea et al. 2006;Gerin et al. 2009;Guzmán et al. 2011;Pety et al. 2012;Ohashi et al. 2013;Le Gal et al. 2017) and laboratory experiments (Alata et al. 2015). The most important question we try to answer is how the dust properties change with physical conditions. Thus, is it possible to understand these observations with grains from the diffuse ISM? Otherwise, is there a viable dust evolution scenario that can explain the observations and is consistent with the physical conditions in the Horsehead?
The paper is organised as follows. In Sect. 2, we describe the previous studies and the observations of the Horsehead. In Sect. 3, we detail the dust model we use, THEMIS, as well as the local dust emission tool, DustEM. We also present the effects of variations in dust properties on its emission in the optically thin case with DustEM in order to disentangle variations in the dust spectrum due to changes in dust properties and those due to radiative transfer effects. In Sect. 4, we present SOC, the 3D radiative transfer code we use, as well as the effect of variations in the dust parameters on dust emission in the optically thick case, in the case of the Horsehead. In Sect. 5, we compare our model with the observations and present the best parameters we obtain. In Sect. 6, we discuss our results and propose a scenario of dust evolution in the Horsehead. Finally, we present in Sect. 7 our conclusions.

A prototypical PDR: the Horsehead
As physical conditions are strongly constrasted and spatially resolved in nearby photodominated regions, they are the ideal place to study dust evolution as a function of physical conditions. First, we introduce the different studies that have been made of the Horsehead; second, we present the observations of the Horse-head obtained with Spitzer and Herschel; third, we describe the density profile that we use to perform radiative transfer across the Horsehead.
Regarding dust, Teyssier et al. (2004) found that although small hydrocarbons are supposed to be photo-destroyed by the intense UV field at the edge of the Horsehead, they are still existing. They suggest that the photo-erosion of carbonaceous nanograins into small hydrocarbons is more efficient than the photodestruction of small hydrocarbons at the Horsehead edge. This scenario is reinforced by observations in Pety et al. (2005) as they found hydrocarbons such as CCH, c-C 3 H 2 and C 4 H in the UV-irradiated outer part of the Horsehead. It is also supported by laboratory experiments on thermal processed and UV-irradiated dust grains analogues (see Smith 1984;Zubko et al. 2004;Alata et al. 2014Alata et al. , 2015Duley et al. 2015). Based on Spitzer observations, Compiègne et al. (2007) proposed a scenario where PAHs survive in HII regions and Compiègne et al. (2008) construed that spectral variations in the mid-IR cannot only be explained by radiative transfer effects and therefore are a consequence of dust evolution across the Horsehead.

Observations with Spitzer and Herschel
We use Spitzer and Herschel observations (see Appendix B) of the Horsehead in 10 photometric bands from 3.6 µm to 500 µm, which cover nearly the entire dust spectrum. The processing of the Spitzer maps is detailed in Bowler et al. (2009). Data were processed in the HIPE environment, with standard Herschel corrections for instrumental effects and glitches. PACS 70 µm and 160 µm maps were obtained after the superposition of two observations with a scan speed of 20 /s whose directions were perpendicular to one another. The overall duration of these observations is 4122 seconds and cover 8.8 × 4.5 of the Horsehead. Concerning SPIRE 250 µm, 350µm and 500 µm, they were obtained after the superposition of two observations with a scan speed of 30 /s whose directions were perpendicular to one other. The overall duration of these observations is 1341 seconds and they cover 8 × 8 of the Horsehead. Striping induced by offsets in the flux calibration from one detector to another was removed using the Scan Map Destriper module included in the HIPE environment.
We study the observed emission profiles through three different cuts across the Horsehead (see Fig. 1). The calibration uncertainty in the IRAC bands (IRAC 3.6 , IRAC 4.5 , IRAC 5.8 and IRAC 8.0 ) is 2 % (Reach et al. 2005), 4 % in MIPS 24 (Engelbracht et al. 2007), 5 % in PACS 70 , 12 % in PACS 160 (Stansberry et al. 2007) and 15 % in the 3 SPIRE bands (Swinyard et al. 2010). In this study, we considered all these errors to be independent of the wavelength to first order. Also, we consider that the emission in all of these 10 bands is coming from dust, which is not completely the case in IRAC 3.6 and IRAC 4.5 . We estimate with a model of atomic and molecular gas in PDRs, the Meudon PDR Code (Le Petit et al. 2006), that gas can contributes less than 10% of the flux. However, this contribution does not affect the bulk of our results hence we consider that the observed emission is dust emission. Nevertheless, one must be careful in the interpretation of the observations as gas emission   can be larger than dust emission in photometric bands covering shorter wavelengths (HST or NIRCAM onboard the JWST).

Density profile across the Horsehead
In this paper, radiative transfer calculations are performed, which require information on the density profile across the Horsehead. We use the profile described in Habart et al. (2005). As the H 2 1-0 S(1) fluorescent emission is very sensitive to both the radiation field and the gas density, they observed this line with the SOFI instrument at the NTT. This observation was combined with previous observations of H α and dust mid-IR emission in order to constrain the density profile at the edge of the Horsehead. Habart et al. (2005) also used CO mm observations from the IRAM 30-m telescope (Abergel et al. 2003;Teyssier et al. 2004) and the Plateau de Bure Interferometer ) as well as 1.2 mm dust continuum emission obtained with MAMBO at the IRAM 30-m telescope (Teyssier et al. 2004) to constrain the density profile in the inner part. All these observations were interpreted with the Meudon PDR Code. This density profile (see Fig. 2, upper panel) was also used in Compiègne et al. (2008) and Arab (2012) and is defined as follows : where : with z the position from the edge of the Horsehead, γ the powerlaw exponent of the gas density profile and z 0 the depth beyond which constant density n 0 is reached.
In this study they also estimated the length of the Horsehead along the line of sight, l PDR . They found that this parameter is constrained to be between 0.1 pc and 0.5 pc. We assume that the density profile is independent of the position along the line of sight (see Fig. 2, bottom panel).

Dust modelling
The interpretation of the multi-wavelength observations of the Horsehead depends on its structure, the incident radiation field and the dust model. We therefore need a dust model and modelling tools to compute dust emission based on the local physical conditions. First, we describe our adopted dust model THEMIS; second, we introduce DustEM, that is used to compute the local dust emission and we describe how dust emission evolves with its properties in the optically thin case using DustEM.

THEMIS
The Heterogeneous dust Evolution Model for Interstellar Solids, THEMIS 1 (e.g., Jones et al. 2013;Köhler et al. 2014;Jones et al. 2017), is based on observational constraints and laboratory measurements on interstellar dust analogues that are amorphous hydrocarbons, a-C(:H); (e.g., Jones 2012a,b,c) and amorphous silicates, a-Sil. This model includes dust evolution through processes such as photo-processing, fragmentation and coagulation resulting from wide variations in the ISM physical condition.
THEMIS for the diffuse ISM (Jones et al. 2013;Köhler et al. 2014;Ysard et al. 2015) is composed of amorphous silicates (a-Sil/a-C) surrounded by a mantle of aromatic-rich carbon, and amorphous hydrocarbon solids which encompasse a-C:H material that are H-rich hence aliphatic-rich and a-C material that are H-poor hence aromatic-rich. Assuming a typical penetration depth of a UV photon in an amorphous carbon grain is about 20 nm (see Fig. 15 in Jones 2012b), carbonaceous grains that are smaller than 20 nm are entirely photo-processed hence aromatic. Larger grains are composed of an aliphatic core and surrounded by an aromatic mantle that is 20 nm thick, which prevents photoprocessing of the core hence allows the core to remain aliphatic. This view provides us with a continuous description of carbonaceous grains from the smallest that mostly contain aromatic cycles and are stochatiscally heated to the largest that are at thermal equilibrium. Details about the size distribution can be found in table A.1. As these grains are composed of either an a-C:H core or a silicate core surrounded in both cases by an aromatic carbonaceous mantle, they are called Core-Mantle grains (CM). 1 THEMIS is available here : https://www.ias.u-psud.fr/themis/ In the dust evolution framework assumed by THEMIS , large grains can form a second mantle either through accretion of C and H atoms, available in the gas phase or through coagulation of a-C nano-grains on the larger grains surfaces. These grains are called Core-Mantle-Mantle (CMM). In denser regions, CMM grains coagulate together to form aggregates  called Aggregate-Mantle-Mantle (AMM) grains. Where the shielding from energetic photons is efficient enough, a mantle of water ice can form around AMM, leading to Aggregated-Mantle-Mantle-Ice (AMMI) grains.
In the following, we use several dust mixtures (see Fig. 1 in Jones et al. 2017). Parameters associated with the size distribution of these dust mixtures can be found in Table A.1 and the size distributions themselves in Fig. 3 (upper panel) with the associated spectra (see Fig. 3, bottom panel), computed with DustEM (see Sect. 3.2). In the near-IR (1 to 5 µm) and mid-IR (5 to 30 µm), dust emission comes mainly from the a-C grains. In the far-IR (from 50 to 500 µm), dust emission comes mainly from a-Sil/a-C and a-C:H/a-C grains.    Table. A.1) for CM-grains (grey line), AMM (blue line) and AMMI (orange line). Black line, dotted-line and dash-dot line correspond to a-C, a-C:H/a-C and a-Sil/a-C respectively. Bottom : associated spectra, computed with DustEM with a radiation field corresponding to a star at 34 600 K with G 0 = 100.

Influence of dust properties on its emission with DustEM
DustEM 2 (Compiègne et al. 2011) is a modelling tool that computes the extinction, the emission and the polarisation of interstellar dust grains heated by photons, in the optically thin case (i.e. no radiative transfer). In order to disentangle the effects of radiative transfer from variations in the dust properties on emission, we study the influence of such variations with DustEM. We modify the following parameters : 1. The a-C abundance, i.e. the a-C mass to gas ratio, M a−C /M H , varying from 0.01 × 10 −2 to 0.20 ×10 −2 with steps of 0.01 × 10 −2 . 2. The a-C minimum size, a min, a−C , varying from 0.4 nm to 0.9 nm with steps of 0.02 nm. 3. The slope of the a-C power law size distribution, α, varying from -6 to -4 with steps of 0.1.
The results are shown in Fig. 5 where the spectra in panels d, e and f are associated with the size distributions in panels a, b and c, respectively. All the spectra are obtained with a radiation field that corresponds to a blackbody at 34600 K scaled so that G 0 = 100 (i.e. the radiation field illuminating the Horsehead).
A decrease in M a−C /M H or an increase in a min, a−C or α leads to a decrease in the smallest a-C grains hence a decrease in the near-IR emission. As the total dust mass is fixed, an increase in a min, a−C or α leads to a redistribution of the dust mass from the smallest to the largest a-C grains hence an increase in the mid-IR emission. In the far-IR, dust emission is unaffected by variations in M a−C /M H , a min, a−C and α as a-C grains are barely responsible for any dust emission at these long wavelengths. However, the dust emission in the far-IR slightly increases with an increase in α as the mass of the largest a-C increases significantly, unlike an increase in a min, a−C .

Radiative transfer modelling within the Horsehead
The Horsehead is an optically thick region that requires a radiative transfer modelling to properly interpret our multiwavelength observations. We present the 3D radiative transfer code SOC we use in this study. Performing radiative transfer is time consuming, and so we here explore the influence of the Horsehead length along the line of sight l PDR , and dust properties (i.e. M a−C /M H , a min, a−C , α) on dust emission after radiative transfer calculations.

Radiative transfer code : SOC
SOC is a 3D Monte Carlo radiative transfer code, parallelised using OpenCL libraries (Juvela 2019), that computes dust emission and scattering. SOC has been benchmarked with other radiative transfer codes in Gordon et al. (2017) and used in Juvela et al. (2018aJuvela et al. ( ,b, 2019. The radiation field corresponds to that of a blackbody at 34 600 K produced by a star to which a dilution factor has been applied to obtain G 0 = 100 at the Horsehead edge. This radiation field is estimated on a logarithmic grid of 334 frequencies that extends from 3×10 9 Hz to 3×10 16 Hz. As the Horsehead edge is located outside the HII region, there are no photons above 13.6 eV hence we applied the Lyman cut to the radiation field that is 2 DustEM is available here : http://www.ias.u-psud.fr/DUSTEM  heating the Horsehead edge. Each frequency is simulated with 10 6 photons. In SOC, clouds can be defined on regular Cartesian grids or octree grids. In this paper, we model the Horsehead using a cartesian grid that contains N X × N Y × N Z cubes that measure 0.0025 pc per side. N X is equal to 77 and corresponds to the number of cubes along the Horsehead-star axis. N Y is equal to 7 and corresponds to the number of cubes along the axis perpendicular to the Horsehead-star axis and the line-of-sight axis (i.e. the observer-Horsehead axis). N Z correponds to the number of cubes in the Horsehead along the line-of-sight hence depends on the value of l PDR : N Z = l PDR / 0.0025 pc. For each cube, we associate a value of gas density as described in Sect. 2.3.
In our study, we compute only dust emission as regardless of the dust properties, dust scattering contributes up to less than 1 % to the total dust brightness in the near-IR photometric bands. After the integration along the line-of-sight, dust emission profiles across the Horsehead are integrated into the different photometric bands and convolved with the PSFs.

Influence of l PDR on dust emission
In the following, we study dust emission at two positions : the near-IR peak position (NIR PP) in the Horsehead and the far-IR peak position (FIR PP). These positions corresponds respectively to the peak of emission in IRAC 3.6 and SPIRE 500 , shown in the Fig. 2. Also, in order to lighten the reading of the results obtained, we introduce I mod, max (i) = max I mod, i (z) where I mod, i (z) is the dust modelled emission in the i-th band at the position z along the cut.
Whether it is at the NIR PP or at the FIR PP, dust emission increases in all bands with l PDR (see Fig. 4, top and middle panels) since the dust mass increases along the line of sight as the column density 3 increases with l PDR . One may also note that dust emission increases linearly with l PDR (see Fig. 4, bottom panel) revealing that dust self-absorption, which depends on both the column density and the wavelength, is negligible at these wavelengths in the l PDR range we are considering. Consequently, we can consider that the intensity increases linearly with l PDR in the near, mid and far-IR and does not affect the shape of the dust spectrum. In the following, we therefore consider l PDR as a multiplying factor on the dust spectrum.

Influence of dust properties on dust emission after radiative transfer
Conversely to Sect. 3.2 where we study the influence of dust properties on dust emission in the optically thin case, we study here the influence of these properties in the optically thick case by performing a radiative transfer calculation. These results are shown in Fig. 5 where the spectra in panels g/j , h/k and i/l are respectively associated with the size distributions in panel a, b and c. Spectra in panels g, h and i are located at the NIR PP and those in panels j, k and l are located at the FIR PP. Dust grains are warmer at the NIR PP (see Fig. 5, panels g, h and i) than at the FIR PP (panels j, k and l) as the maximum intensity shifts towards longer wavelengths. This effect is due to the damping of the radiation field with increasing depth into the Horsehead.
One may note that at the NIR PP, dust emission in the far-IR does not vary with M a−C /M H (see Fig. 5, panel g), conversely to that is seen in the inner part (panel j). As dust emission in the far-IR is unaffected by variations in M a−C /M H in the optically thin case (see Sect. 3.2), this is strictly a radiative transfer effect. As the a-C grains bear a large fraction of the total dust crosssection, an increase in M a−C /M H increases significantly the extinction. Therefore as M a−C /M H increases, the radiation field is increasingly damped at the NIR PP and there are less photons available at the FIR PP to heat the larger grains. Indeed, as we can see in panel j, the wavelength associates with the maximum of emission shifts towards longer wavelengths with an increase in M a−C /M H . Therefore, dust emission in the far-IR varies with M a−C /M H , due to radiative transfer effects.
Regarding the other changes in the spectra, they are due to variations in dust properties and are explained in Sect. 3.2.

With evolved grains
Previously, we used only CM-grains throughout the Horsehead.
To study the influence of dust evolution on the emission across the Horsehead, we use CM-grains with modified size distributions (i.e. CM-grains with values of M a−C /M H , a min, a−C and α 3 N H (z) = n H (z) l PDR that differ from the diffuse ISM) in the outer part of the Horsehead where the dust is likely to be more diffuse ISM-like, and aggregate-grains (AMM, AMMI) above a density threshold of 7 × 10 4 H.cm −3 , where dust grains are assumed to be coagulated. In order to simplify our study, we define 3 different cases depending on the dust we use : • Case a : CM-grains with modified size distributions across all the Horsehead. Dust modelled emission profiles for the three cases are shown in Fig. 6.
As the maximum of emission in the near and mid-infrared is located in the outer part of the Horsehead, there is no modification of dust emission at these wavelengths as we always use modified CM grains. Regarding dust emission in the far-infrared, this increases when coagulated (AMM, AMMI) dust grains are used because they are more emissive. One may note that AMMI are more emissive than AMM as the dust mass in AMMI is larger than in AMM because of the ice mantle.

Comparison with observations
In this section, we constrain our dust model with the observations. First, we present our results using diffuse ISM-like dust (i.e. CM-grains); second, we introduce the methodology we use in the following parts; third, we constrain the four parameters M a−C /M H , a min, a−C , α and l PDR for the 3 cases of evolved grains as defined in Sect. 4.4 and across the 3 cuts (see Fig. 1).

Diffuse case
The results are shown in Fig. 7. The 10 upper panels correspond to the modelled emission across the Horsehead using CM-grains, with l PDR varying from 0.1 pc to 0.5 pc, for the 10 photometric bands. The observed emission is shown for cut 2. The bottom panels show the corresponding ratios of maximum observed and modelled intensities.
Regardless of the cut considered, it is not possible to simultaneously fit the observations in all the photometric bands (see Fig. 7, upper panel), whatever the l PDR value. With l PDR = 0.1 pc, we are able to roughly reproduce the observations in the near and mid-infrared but in the far-infrared, the modelled dust emission is too low by a factor ∼ 10 (see Fig. 7, bottom panels). With l PDR = 0.5 pc, we are able to reproduce the observations in the far-infrared but in the near and mid-infrared, the modelled dust emission is too high by a factor of at least ∼ 10.
If l PDR is higher than 0.10 pc (see Sect. 2.3), near and midinfrared modelled dust emission will always be too high, which implies reducing the dust abundance that is responsible for the emission at these wavelengths hence decreasing the a-C dust-togas ratio, M a−C /M H (see Sect. 4.3). On the other hand, the ratio between the modelled dust emission and the observations is not the same in the five near and mid-IR bands. One is therefore required to change the shape of the spectrum by varying a min, a−C and α (see Sect. 4.3).
To summarise, it is not possible to reproduce the observations across any of the three cuts in the Horsehead for any value  . Panels d-l show dust spectra for the size distributions shown in the corresponding top-row panels. The spectra are computed with DustEM (panels d-f) and with radiative transfer for the NIR (panels g-i) and the FIR (panels j-l) peak positions. of l PDR , using only diffuse ISM-like dust. We must therefore consider evolved dust.

Methodology
For the sake of reducing computation time, instead of exploring the 4D-space defined by M a−C /M H , a min, a−C , α and l PDR , we ex-plore the 3D-space defined by M a−C /M H , a min, a−C and α as variation in l PDR does not affect the shape of the dust spectrum (see Sect. 4.2), conversely to variations in M a−C /M H , a min, a−C and α (see Sect. 4.3). Therefore, l PDR can be adjusted after the fact.
Adjusting the shape of the modelled dust spectra to the observed dust spectra means that the ratio between I obs, max (i) and I mod, max (i) has to be roughly the same in every band. Therefore, Article number, page 7 of 17 A&A proofs: manuscript no. aa we minimise the following parameter : with where r obs is the relative error for each filter and defined in Sect. 2.2 and I obs, max (i) = max I obs, i (z) with I obs, i (z), the dust observed in the i-th band at the position z along the cut.
The following procedure is thus applied : 1. We constrain M a−C /M H , a min, a−C and α with a fixed l PDR in order to adjust the shape of the modelled dust spectrum to the observed dust spectrum by minimising χ 2 . 2. We use the dust properties (M a−C /M H , a min, a−C and α) associated with χ 2 min (i.e. the minimum value of χ 2 in the 3D-space defined by M a−C /M H , a min, a−C and α) and we adjust the overall modelled dust spectrum to the observed dust spectrum by multiplying the flux in all bands by the same factor to get µ = 1, which constrains l PDR .
We choose to remove the IRAC 4.5 and PACS 70 bands because it is not possible to simultaneously fit the observations in the 10 bands with these 2 included. We discuss this decision further in Sect. 6.1.

Constrain l PDR and dust properties M a−C /M H , a min, a−C and α
First, we study the χ 2 distribution in the 3D-space (M a−C /M H , a min, a−C , α) for each of the 3 cuts and the 3 cases, in order to obtain the best set of parameters in these 9 cases. The 3D-space is defined as follows : 1. M a−C /M H varying from 0.001 × 10 −2 to 0.041 × 10 −2 with steps of 0.002 × 10 −2 .
2. a min, a−C varying from 5 nm to 10 nm with steps of 0.25 nm.
Second, we study the χ 2 distribution in the 2D-spaces (a min, a−C , α), (M a−C /M H , a min, a−C ) and (M a−C /M H , α). Finally, we conclude with the comparison between the observed and modelled dust emission profiles for each of the 3 cuts and the 3 cases with the best sets of parameters. For more clarity, we define χ 2 min, 2D (M a−C /M H ) that is the minimum value of χ 2 in the 2D-space (a min, a−C , α) for a given value of M a−C /M H . We also define χ 2 min , that is the minimum value of χ 2 in the 3D-space ( First and foremost, M a−C /M H is between 0.01 ×10 −2 and 0.03 ×10 −2 , i.e. between 6 to 10 times lower than in the diffuse ISM (0.17 ×10 −2 ) regardless of the cut or the case considered. Second, a min, a−C is between 0.8 and 0.925 nm, i.e. between 2 and 2.25 times larger than in the diffuse ISM (0.4 nm). Third, α is between -7 and -5.5, i.e. between 1.1 to 1.4 times lower than in the diffuse ISM (-5).
One may note that regardless of the cut, M a−C /M H increases from case a to case c (see Fig. 8). In case a, we only use modified CM grains (i.e. CM grains with values of M a−C /M H , a min, a−C and α that differ from the diffuse ISM) in the outer part and in the inner part of the Horsehead but we use modified CM-grains in the outer part of the Horsehead and AMMI in the inner part in case c. As AMMI are more emissive in the far-IR than CM grains, emission in this wavelength range must decrease to fit the observations, which is achievable by reducing l PDR (see Sect. 4.2), hence l PDR decreases from case a to case c. This decrease in l PDR implies a decrease in emission in the near and mid-IR hence M a−C /M H must increase to counterbalance this variation. We show in Fig. 9, the χ 2 distribution for cut 2 in the 2D-spaces (a min, a−C , α), (M a−C /M H , a min, a−C ) and (M a−C /M H , α). We choose to focus on only one cut as we are interested in the behaviour of the χ 2 distribution here, which is the same regardless of the cut.
The most important result is that, regardless of the case, there is a unique minimum in all 2D-spaces. Also, as explained in Sect. 3.2, a decrease in M a−C /M H is, to first order, similar to an increase in a min, a−C and α, regarding dust emission in the near and mid-IR. An increase in a min, a−C is therefore counterbalanced by a decrease in α to keep low values of χ 2 , and an increase in M a−C /M H is counterbalanced by an increase in a min, a−C and in α hence explaining the banana-shape of the low χ 2 values in each of the 2D-spaces.
It can also be seen that from case a to case c, the position of χ 2 minimum value moves. From case a to case c, dust emission in the far-IR increases (see Fig. 6) hence this effect is counterbalanced by a decrease in l PDR that also reduces dust emission in the near and mid-IR. To compensate for this decrease in dust emission in the near and mid-IR, the value of M a−C /M H associated with the χ 2 minimum value increases from case a to case c in the 2D-spaces (M a−C /M H , a min, a−C ) and (M a−C /M H , α). In the 2Dspace (a min, a−C , α), this effect is counterbalanced by a decrease in α and an increase in a min, a−C .

Comparison between dust modelled emission and dust observed emission profiles
Here, we use the best set of parameters (M a−C /M H , a min, a−C , α and l PDR ) that are listed in Table. 1 and compare the modelled emission profiles in the 10 photometric bands for the three cases with the observed emission profiles in the three cuts (see Fig. 10). We focus on three aspects: the maximum of intensity in each of the 10 bands, the position of these maxima and the width of these profiles.
In the near and mid-IR, except in IRAC 4.5 , the maximum emission is well reproduced, regardless the case or the cut. In PACS 70 , although the maximum of emission is never reproduced, the discrepancy between the maximum modelled emission and the maximum observed emission decreases from case a to case c. From SPIRE 250 to SPIRE 500 , the maximum emission is in the error bars, regardless of the case or the cut and the discrepancy between the maximum modelled emission and the maximum observed emission decreases from case a to case c. Regarding PACS 160 , the maximum emission is within the error bars only for case c for cuts 2 and 3 but never for cut 1, regardless of the case.
Concerning the position of the maximum emission, it is well reproduced from IRAC 3.6 to PACS 70 regardless of the cut and the case. Regarding the cut 1, there is a small discrepancy between the position of the maximum emission and the position of the observed emission from PACS 70 to SPIRE 500 . For cut 2, there is the same discrepancy in SPIRE 350 and SPIRE 500 , regardless of the case. For cut 3, all the positions are well reproduced.
Regarding the width of the profiles, they are well reproduced from IRAC 3.6 to PACS 160 but slightly different from SPIRE 250 to SPIRE 500 , which could be due to large structures in the Horsehead.
To summarise, the observed dust emission is well reproduced in the near and mid-IR, except in IRAC 4.5 , regardless of the case and the cut. In the far-IR, the discrepancy between observed dust emission and modelled dust emission decreases from case a to case c.

Discussion
First, we discuss the discrepancy between the dust modelled emission and the dust observed emission in IRAC 4.5 and in PACS 70 ; second, the results obtained are described; third, we propose a scenario of dust evolution in agreement with the results obtained. We end with a discussion of dust processing timescales in support to this scenario.

Discrepancy in IRAC 4.5 and PACS 70
In IRAC 4.5 , the modelled dust emission is always overestimated (see Fig. 10) by a factor 2 to 4. As this filter covers the dust continuum and the wings of the IR bands from a-C:H nano-grains, this suggests that the wings of the IR bands in this region are different (i.e., weaker and/or narrower, see for instance Boutéraon et al. (2019) for more details about the variability of the IR band widths) from those in the diffuse ISM. Indeed, we are here looking at dust that is evolving from dense cloud dust in response to interaction with UV photons.
Moreover, a-C:H nano-grains freshly produced may not yet have time to be entirely photo-processed, hence have a large band-to-continuum ratio because of their high fraction of aliphatic bonds, as opposed to aromatic bonds. As discussed in Jones et al. (2013), this requires a-C:H nano-grains with a band gap larger than 0.1 eV, the value adopted in the diffuse ISM.
However, as we do not have dust spectroscopic informations of the Horsehead in the near-IR, we are not able to answer to these previous questions. In the near future, JWST spectroscopic data should allow us to understand such changes in the structure of a-C:H nanograins.
In PACS 70 , models always overestimate the emission (see Fig. 10) by a factor 3 to 4. This suggests that large grains (a-Sil/a-C and a-C:H/a-C for case a, AMM for case b and AMMI for case c) are somewhat too warm and not emissive enough. This is supported by recent laboratory experiments in which the mass absorption coefficient of silicates in the far-IR is larger (up to an order of magnitude) than those currently used in THEMIS (see Fig. 5 in Demyk et al. 2017). As a consequence, the large grains we use in there are probably not emissive enough. The incorporation of these new laboratory results in THEMIS will most likely reduces the discrepancy in PACS 70 .

Main results
Using the 3D radiative transfer code SOC together with the dust model THEMIS, we can reproduce the Horsehead observations in 8 of the 10 photometric bands of Spitzer and Herschel. The main results for the outer part of the Horsehead are the following : 1. The nano-grains (i.e. a-C grains) dust-to-gas mass ratio, M a−C /M H , is 6 to 10 times lower than in the diffuse ISM. 2. The minimum size of the nano-grains, a min, a−C , is 2 to 2.25 times larger than in the diffuse ISM. 3. The power-law exponent of the nano-grains size distribution, α, is 1.1 to 1.4 times lower than in the diffuse ISM, i.e. the size distribution is steeper.
The best size distributions for the three cuts and case c are shown in Fig. 11. Concerning the inner part of the Horsehead, we tested 3 different kinds of dust, diffuse ISM-like dust (CM) with modified size distributions in case a, aggregates of grains (AMM) in case b, aggregates of grains with ice mantles (AMMI) in case c. At long wavelengths (from PACS 160 to SPIRE 500 ) The results are significantly better when using AMMI instead of CM grains. Regarding PACS 70 , even if we are not able to reproduce the observed emission with our model, using aggregates (AMM/AMMI) instead of diffuse ISM-like dust (CM) with modified size distributions, significantly improves the fit in this band.
Finally, the length of the Horsehead along the line-of-sight, l PDR , is found to be within the range of 0.26 and 0.30 pc which is in agreement with previous gas studies (Habart et al. 2005).

Dust evolution scenario
Our results show significant variations of the dust size distribution and in the following we outline a possible scenario of dust evolution across the Horsehead interface. Given the strong incident radiation field, we assume that the dominant process is the exposure of dust grains from the dense molecular cloud (the inner region) to the UV light of σ-Ori. This suggests two major photo-processing sequences: (i) the partial fragmentation of aggregate grains from the inner region and (ii) the destruction of the smallest a-C:H nano-grains. We discuss the significance of these sequences by comparing their timescales to the advection timescale τ a , i.e., the time that the incident UV light needs to heat up and dissociate the molecular gas at the cloud border.
The advection timescale is defined as τ a = L/v DF where L ∼ 0.05 pc is the width of the outer part of the Horsehead and v DF ∼ 0.5 km.s −1 is the velocity of the dissociation front (Hollenbach & Tielens 1999). With these values, we find τ a ∼ 10 5 years.
Due to the lack of studies, we take the photo-darkening timescale as a lower limit to photo-fragmentation of aggregate grains and photo-destruction of a-C nano-grains, described by τ ph . Indeed, photo-darkening involves the dissociation of CHbonds, a process that is more likely faster than the breaking of CC-bonds that must occur in photo-fragmentation (Jones & Habart 2015). We thus express τ ph at the cloud edge in terms of the photo-darkening rate Λ pd  : where F 0 UV 3.8 × 10 9 photons.s −1 .cm −2 is the unattenuated UV field, (a) = min(1, 2 a[nm] ) is a size-dependent photo-darkening efficiency, σ CH 10 −19 cm 2 is the CH bond photo-dissociation cross-section and Q abs (a), the dust absorption efficiency which depends almost solely on the radius in the UV range. In the case of AMMI, τ ph is larger in reality because the ice mantle needs to be vaporised first but we do not take into account this effect as the time τ ph we estimate is already a lower limit.
We show τ ph (a) in Fig. 12 for CM, AMM and AMMI. As discussed in Ysard et al. (2016), more than 50 % of the AMM(I) dust mass is contained in grains larger than 250 nm. From this figure, one can see that aggregate grains can be photofragmented because τ ph ∼ τ a . One can also see that a-C nanograins can be efficiently destroyed as τ ph < τ a for a-C nanograins. Similar results were found by Alata et al. (2014), from laboratory experiments on a-C:H grain analogues, later applied to the Horsehead (Alata et al. 2015).
From this analysis emerges the following scenario. Within an advection timescale, the a-C nano-grains formed by fragmentation of aggregate grains are also partially destroyed by UV photons. This naturally explains the depletion of a-C:H grains around a=10 nm seen in Fig. 11). We note that the size distribution of these freshly formed small grains is significanlty different from the diffuse ISM case (blue curve in Fig. 11)). This evolved size distribution could reflect the photo-evaporated layer described by Bron et al. (2018).

Conclusion
With Herschel and Spitzer data, we studied the Horsehead using 10 photometric bands, from 3.6 µm to 500 µm, covering the entire dust spectrum. We modelled the dust emission across the Horsehead using the THEMIS dust model together with the 3D radiative transfer code SOC.
We show that it is not possible to reproduce the observations in the Horsehead using dust grains from the diffuse ISM hence the necessity to modify their size distributions and compositions. Dust therefore evolves across the Horsehead.
In the outer part of the Horsehead, the a-C nano-grains dustto-gas ratio is 6 to 10 times lower and their minimum size, 2 to 2.25 times larger than in the diffuse ISM. The power law of the size distribution is steeper than in the diffuse ISM. In the inner part of the Horsehead, we show that using aggregate grains with or without ice mantles significantly reduces the discrepancy between our model and the observations. The discrepancy between the observations and our model at 4.5 µm could be due to the shape of the aromatic band wings whence the overestimation of the dust modelled emission. We also find that large grains are too warm because our modelled dust emission at 70 µm is overestimated. However, laboratory studies show that large silicate grains are more emissive than those used in dust models hence cooler. These new results will soon be implemented in THEMIS.
Based on a time-scale analysis, we propose a scenario where the a-C nano-grains form by the partial photo-fragmentation of aggregate grains and are processed by the UV photons, leading to a size distribution depleted in grains of size from 5 to 10 nm. In the denser regions of the Horsehead, the dust composition is typical of dense clouds.
Spectroscopic observations of the Horsehead are required to go further on the structure and size distribution of a-C nanograins. Indeed, the observations with the JWST will, for the first time, spatially resolve the individual IR dust signatures across the Horsehead, offering an unprecedented look at the evolution of the interstellar matter in photon-dominated regions.