Measuring dust in core-collapse supernovae with a Bayesian approach to line profile modelling
Antonia Bevan




firstpage–lastpage 2018

Analysis of trends in experimental observables and reconstruction of the implosion dynamics for direct-drive cryogenic targets on OMEGA
    E. M. Campbell
    Accepted on 31 July 2018
=======================================================================================================================================








Optical and near-IR (NIR) line profiles of many ageing core-collapse supernovae (CCSNe) exhibit an apparently asymmetric bluewards shift often attributed to greater extinction by internal dust of redshifted radiation emitted from the receding regions of the SN ejecta. The damocles Monte Carlo line radiative transfer code models the extent and shape of these dust-affected line profiles to determine the dust mass that has condensed, in addition to other properties of the dusty ejecta. I present here the application of an affine invariant Markov Chain Monte Carlo (MCMC) ensemble sampler (emcee) to the damocles code in order to investigate the multi-dimensional parameter space rigorously and characterise the posterior probability distribution. A likelihood function is formulated that handles both Monte Carlo and observational uncertainties.  This Bayesian approach is applied to four simulated line profiles in order to test the method and investigate its efficacy.  The majority of parameters can be tightly constrained using this method, and a strong (predictable) dependence between the grain size and the dust mass is quantified. The new approach is also applied to the Hα line and [O i] 6300, 6363 Å doublet of SN 1987A at 714 d post-outburst, re-examining a previous 5-dimensional smooth model and also investigating a new, more complex, 10-dimensional model that treats both features simultaneously.  The dust mass, dust grain size and a range of other parameters can be well constrained using this technique, representing a significant improvement over the previous manual approach.





	line: profiles –
	radiative transfer –
	methods: statistical –
	supernovae: general –
	supernovae: individual: SN 1987A  –  
	ISM: supernova remnants.






§ INTRODUCTION


There remain numerous questions surrounding the formation of dust in the universe. Significant challenges are still faced in the determination of dust formation rates, mechanisms and environments.  Motivated by seemingly inexplicably large masses of dust observed at high redshifts <cit.>, there is a widespread 
desire to understand the nature of the primary sources of dust in the universe.  

Core-collapse supernovae (CCSNe) are known to produce dust in their ageing ejecta. Theoretical models  predict that CCSNe are capable of producing >0.1M_⊙ of ejecta-condensed dust <cit.>. To date, however, dust masses in the majority of these objects have been inferred via fits to their near-infrared (NIR) and mid-infrared (MIR) spectral energy distributions (SEDs), which trace only warm dust.  Warm dust masses up to ∼10^-3 M_⊙ have been detected at late times (>1 yr) in several CCSNe <cit.>.  However, a few objects have also been observed in the far-IR, allowing their full SEDs to be fitted and therefore tracing the presence of cold dust as well as warm and hot dust.  Using this technique, dust masses ≳0.1 M_⊙ have been estimated to have formed in SN 1987A, Cassiopeia A and the Crab Nebula <cit.>. More recently, it has been suggested that a significant mass of dust (0.08 – 0.9 M_⊙) has also formed in the Galactic supernova remnant G54.1+0.3 <cit.> as well as very large masses of dust (>1 M_⊙) in a number of other Galactic supernova remnants (Chawner et al. in prep.). An average net dust production rate of 0.1 – 1.0 M_⊙ per CCSN is required to account for the dust masses observed in the early universe <cit.>. 
The dust budget problem would therefore be resolved if 
these few objects were, in fact, representative of the wider CCSN population and the dust was able to survive the passage of the reverse shock <cit.>. A larger sample of CCSN dust mass estimates is therefore required.


Following the end of the Herschel mission in 2013, there will be a long wait for instruments that are capable of detecting cold dust emission at far-IR wavelengths and so an alternative approach is needed.  The damocles Monte Carlo line radiative transfer code predicts dust masses in the ejecta of CCSNe by modelling the red-blue asymmetry frequently observed in their optical and NIR emission lines <cit.>.  This asymmetry is due to the condensation of dust in the ejecta causing redshifted radiation from the receding regions of the supernova to experience greater extinction than blueshifted radiation emitted from the approaching regions <cit.>. In addition to providing an alternative method for tracing both warm and cold dust in the ejecta, line profile modelling has the added advantage of tracing only newly-condensed dust within the ejecta. Pre-existing circumstellar dust may  contribute to the observed flux in the IR but the red and blue components of optical or NIR lines emitted from within the ejecta will be similarly attenuated by the surrounding circumstellar dust, i.e. any dust-induced red-blue asymmetry must be solely a result of internal, ejecta-condensed dust. 

The approach also allows other properties of the dust to be determined. Of particular interest is the dust grain size distribution.   Regardless of the masses of dust that form in the ejecta of CCSNe, the grains will eventually be subject to a reverse shock that will pass back through the ejecta, potentially destroying these newly-formed grains and significantly diminishing the dust mass that has formed.  The size of dust grains that condense within the supernova ejecta  determine their likelihood of survival. An understanding of dust grain sizes in CCSNe is therefore critical to determining the relative contribution of CCSNe to dust production in the universe.

B16 applied the damocles Monte Carlo code to the Hα and [Oi] 6300,6363 Å lines of SN 1987A between 714 d and 3500 d post-explosion. A steady increase in the ejecta dust mass over this period was inferred with a predicted current dust mass of 0.8 M_⊙, consistent with dust mass estimates derived from SED fitting and modelling <cit.>.  The damocles code was also applied to the late-time optical line profiles of SN 1980K, SN 1993J and Cas A by <cit.>.  Dust masses of 0.12 -– 0.3 M_⊙ at 30 yr, 0.08 –- 0.18 M_⊙ at 16 yr and ∼1.1 M_⊙ at ∼330 yr were predicted respectively <cit.>. Clearly, further examples are needed in order to establish whether this apparent trend towards larger (≳0.1 M_⊙) dust masses is an accurate representation of dust formation in CCSNe more generally. 

In working towards the overall goal of understanding the masses and properties of dust in the ejecta of CCSNe, I have explored the implementation of a Bayesian methodology for line profile fitting. The fundamental power of Bayesian statistics in providing a framework to understand the probability of a model when the the data is known has been increasingly exploited in astronomy over the last twenty years <cit.>. In particular, Markov Chain Monte Carlo (MCMC) methods have provided efficient, robust and rigorous procedures with which to explore highly multi-dimensional parameter spaces and quantify posterior probability distributions.  Increasingly available computing power has allowed these methods to be employed in a wide variety of fields with impressive results that yield significantly more insight than can be gained from a single best-fitting set of parameters. <cit.> presents a comprehensive review of MCMC methods for Bayesian data analysis in astronomy.  

 I have applied an affine invariant ensemble sampler <cit.> to the damocles code in order to map the multi-dimensional posterior probability distribution of a range of models and parameter spaces. I initially employed the sampler to model four simulated, or `theoretical', line profiles (generated by damocles) that were deemed representative of observed line profiles of CCSNe at late times but for which the `true solution' was known (models A1 - A4).  
I also revisited the models by B16 of the Hα line and the [O i] 6300,6363 Å doublet of SN 1987A at 714 d.  The new approach was applied to a smooth, 5-dimensional model based on their work (model B). A new model (model C) was also explored which treated both emission features simultaneously and used the newly-applied Bayesian methodology to characterise a significantly more complex 10-dimensional  parameter space.  
The ultimate goal is to assess the validity of this approach with regard to its application to both archival and future datasets, with a view to significantly expanding the current range of dust mass estimates for CCSNe years after outburst.  

In Section <ref>, the formulation of the problem is presented along with a discussion of how Monte Carlo and observational uncertainties are handled and a brief description of the affine invariant ensemble sampler.  The adopted priors, variable parameters and posterior distributions for all models are presented in Section <ref>.  I discuss the implications of these results in Section <ref> and compare them to results obtained from previous manual line profile fitting using damocles, as well as results obtained from SED fitting.  The constraints placed on other parameters such as dust grain size and dust density distribution are also discussed. I emphasise the potential for future application to other objects before summarising and concluding in Section <ref>.




	






§ FORMULATION OF THE BAYESIAN APPROACH






 §.§ damocles

damocles is a Monte Carlo radiative transfer code that models the effects of dust, composed of any combination of species and grain size distributions, on optical and NIR emission lines emitted from the expanding ejecta of a late-time (>1 yr) supernova.  For full details of the code and its testing, please see B16. By default, both the emissivity distribution and the dust distribution follow smooth  radial power-law distributions although any arbitrary distribution may be specified by providing the appropriate grid.  damocles will also treat a variety of clumping structures as specified by a clumped dust mass fraction, volume filling factor, clump size and clump power-law distribution. The emissivity distribution may also initially be clumped.  The code has a large number of variable parameters ranging from 5 dimensions in the simplest, smooth models to >20 in the most complex cases.  




 §.§ The Bayesian approach


The aim is to map the posterior probability distribution based on the observations and our prior understanding of the physical situation.  The posterior is defined by Bayes' Theorem as


    P(θ |  D) = P(θ )  P( D | θ)/P(  D)


where D represents the data that we wish to analyse (in our case, the observed or simulated line profile) and θ represents the parameters of our model. P(θ ) therefore represents our prior understanding of the probability of the model parameters (the prior), P( D | θ) is the probability of obtaining the data for a given set of model parameters (the likelihood) and P(D) is the probability of the data for all models (the evidence). Since P( D) is independent of θ, we will only be interested in the scaled posterior as defined by 


    P(θ |  D) ∝ P(θ )   P( D | θ)



The posterior distribution will allow us to understand relationships between the parameters and to visualise which are the most likely regions of parameter space.  The aim is not to identify the single `best-fitting' model but to map the variation of likelihood across the entire space.   The prior is the probability before looking at any data.  It can be driven by theoretical models, previous observations or physical intuition. The likelihood is, practically, a mechanism for forward modelling i.e. simulating the data given a model and its parameters. It is proportional to exp(-χ^2/2), where χ^2 is the standard metric typically used to compare data and models in frequentist techniques.

In order to characterise the target posterior distribution, we may draw samples from across the parameter space. A single sample in parameter space is translated to a point in the target posterior distribution via Equation <ref>.  A likelihood function that describes the relationship between the model and the data must therefore be defined, and a prior probability distribution for each parameter must also be specified based on our current knowledge (e.g. physical constraints). Once defined, the ensemble sampler can be employed with this likelihood function in order to map the complete posterior distribution.  





 §.§ Affine Invariant MCMC Ensemble Sampler


There are numerous MCMC algorithms that work out how to sample points in parameter space efficiently in order to converge on a stable solution as quickly as possible.  In this work, I used the Python package `emcee' <cit.>.  This package uses an affine invariant ensemble sampler as described by <cit.>, to which publication the reader is referred for full details of the algorithm. I present a summary below.

The ensemble sampler acts on a collection of points (n-dimensional position vectors in parameter space) termed `walkers'.  An initial position for each walker is sampled according to a distribution specified by the user.  The likelihood P( D | θ)  of the model corresponding to the current set of parameters is calculated. A new point in parameter space is sampled based on the current positions of the other  walkers and the likelihood of this new point is also calculated. The ratio of these likelihoods determines whether the position of the walker is updated or not (i.e. the new point is either accepted or rejected). 
As such, the walkers `walk' around the entire parameter space exploring the posterior distribution in such a manner that the value of the posterior distribution in a given region of parameter space is characterised by the density of walkers in that region.    Faster convergence is attained when the walkers are initialised near regions of high likelihood but they will explore the entire space regardless of their initial positions. The extent of the space to be explored is determined by the bounds of the prior distributions (if applicable).




It is, of course, possible simply to grid parameter space, evaluate the likelihood and prior at each point and multiply them to get the posterior. However, whilst this would be exact, it would also be incredibly intensive and likely impossible for >4 dimensions. MCMC methods approximate the posterior by exploring the parameter space intelligently and are therefore a popular alternative. damocles has between 5 and 20 variable parameters, strong degeneracies between certain of these parameters but no multimodality.  MCMC methods are extremely well-suited to this regime and are therefore an ideal choice.

The choice to use this particular MCMC methodology was made for a number of reasons.  Affine transformations are those which preserve the relative positions of points, lines and planes, for example reflection, rotation and scaling are all affine transformations.  This algorithm is designed to be affine invariant such that the parameter space can be `stretched' in order to sample points from a more isotropic distribution.  This ensures that it requires very little tuning in order to obtain good performance and is in contrast to a number of other MCMC algorithms such as Metropolis-Hastings <cit.>. This property makes the algorithm particularly useful for models with parameters that range over significantly different scales as here. In addition to this, the ease of use and implementation (via the emcee package) and its speed and efficiency for problems with dimensionality of this order led to the choice of this algorithm over other available options.



	



 §.§ Formulating the likelihood function


In order to quantify the likelihood of a particular set of model parameters given the observational or simulated data, we must define a function that relates the model to the data.  We here base the likelihood function on the typical χ^2 comparison which is defined as 


    χ^2 = ∑_ i=1^n( f_ mod,i-f_ obs,i) ^2/σ_ i^2


where f_ mod,i and f_ obs,i are the model flux and observed flux in frequency bin i respectively.  σ_i represents the overall uncertainty in frequency bin i and n is the number of frequency bins in the observed line profile.  

There are two primary contributions to the uncertainty σ_i in each bin: there is an inherent uncertainty on the observational data and there is also uncertainty arising from the statistical nature of the Monte Carlo radiative transfer simulation.  



	
   


The observational uncertainty in each frequency bin (σ_ obs,i) is usually determined when data are reduced and is often included in addition to fluxes in flux-calibrated spectral data files.  However, in a number of cases, particularly in cases of older, archival data, accurate uncertainties are not available.  In these cases, a region of flat continuum may be selected and the observational uncertainty estimated from the variance of fluxes in that region.  A number of different `flat' regions of the spectrum should be sampled and the mean variance calculated.  This value may be used as an approximation to σ_ obs^2 which is assumed to be constant over the whole line profile.  Whilst this is an approximation, over the small wavelength ranges of interest for a single line profile, it is generally reasonable to assume that there is little variation in the uncertainty although care should be taken if there was significant contamination to the profile by, for example, sky lines.  Where accurate errors are available, or a full set of raw observations is available such that accurate uncertainties can be calculated, these should be adopted.   The observational error should ideally include accurately calculated uncertainties from as many sources of observational uncertainty as possible (instrumental noise, calibration errors etc.) but particular care should be paid when handling  continuum subtraction.  This can be a significant factor that influences the results of line profile fitting and ideally should be included as a free parameter in any modelling.  This is discussed further in Section <ref>.
 
Each modelled line profile is also inherently uncertain due to the stochastic nature of Monte Carlo simulations.  The Monte Carlo uncertainty can be quantified analytically.  damocles propagates weighted energy packets through a dusty medium.  Once it has escaped, the weighted packet is added to the appropriate frequency bin.  Each frequency bin therefore receives weighted packets at a rate that is determined by the properties of the model.  Statistically, this is described by a compound Poisson distribution (i.e. identically, independently distributed weights arrive at a rate described by a Poisson distribution). In the limit of a large number of packets (as here), the compound Poisson distribution can be approximated by a normal distribution with associated Monte Carlo uncertainty in each frequency bin σ_ mod,i described by 


    σ_ mod,i = f_ obs√(∑_ j=1^n_ i  w_ ij^2)/∑_ i,j  w_ ij


where f_ obs is the total integrated flux of the observed line profile, n_i is the number of packets in bin i and w_ij is the weight of the j^th packet to arrive in bin i. The model flux in the i^th frequency bin is given by f_ mod,i = f_ obs∑_ j=1^n_ i  w_ ij/  ∑_ i,j  w_ ij. The fluxes are therefore scaled such that the total integrated flux of the model profile is equal to that of the observed profile. 

Since both the observational and Monte Carlo uncertainties can be assumed to follow normal distributions, the total error in the likelihood function (see Equation <ref>) therefore also follows a normal distribution and is described by 


    σ_i^2 =  σ_ obs,i^2 + σ_ mod,i^2


thus fully defining the likelihood function.




 §.§ Computational implementation

The Python package `emcee' <cit.> was coupled to the Fortran 95 damocles code using the F2PY Fortran-to-Python interface generator <cit.>.  Samples in parameter space are generated in Python, passed to damocles where the full Monte Carlo radiative transfer calculation is performed, before the model line profile is passed back to Python. The likelihood and prior are calculated and the algorithm progresses accordingly.  damocles is parallelised using OpenMP <cit.> and models were run on an 88-core machine with Intel Xeon CPU E5-4669 2.20GHz processors using half its capacity. The most complex, 10-dimensional model took approximately 2 weeks to converge (∼20,000 steps). 













§ RESULTS





 §.§ MCMC models of simulated line profiles (A)



I initially considered a number of simulated line profiles for which the true parameters were known.  Four simulated line profiles were produced that are similar to the types of asymmetric dust-affected optical and NIR line profiles observed in the spectra of late-time CCSNe with regard to the extent of their asymmetries, their shape and notable features.  The parameters used to generate these line profiles are described in Table <ref> with graphical representations of the geometrical structures presented in Figure <ref> and the profiles presented in Figure <ref>. All four profiles exhibit a blueshifted peak flux due to increased absorption by dust of redshifted radiation and an extended red scattering wing caused by repeated dust scattering events. Three of the profiles also display a `shoulder' or second peak at the position of the minimum radial velocity on the red side. This has been previously noted by B16 and occurs in scenarios with steeper dust and gas density distributions as a result of significant absorption in the central regions of the profile. 

A monochromatic line at 6563 Å (Hα) was modelled in each case assuming a post-explosion date of 1000 d and a symmetric shell ejecta in homologous expansion (v ∝ r) with maximum velocity at R_ out of v_ max = 4000 km s^-1.  The models adopted an intrinsic smooth power-law emissivity distribution which was coupled to the square of the density distribution of the dust (ρ_ d) such that for ρ_ d∝ r^-β the emissivity distribution followed i ∝ r^-2β, as appropriate for recombination lines assuming a constant dust-to-gas mass ratio.  Two schematics that illustrate the structure of the shell geometries and the smooth, radial power-law dust density distributions are presented in Figure <ref> using models A1 and A2 as examples. 

 100% amorphous carbon grains were used and the optical constants presented by <cit.> were adopted. The physical extent of the ejecta in each model was determined within the code based on the post-explosion time and the specified maximum velocity (R_ out = v_ maxt). The four simulated line profiles that were selected for investigation are presented in Figure <ref> with the parameters used to generate them detailed in Table <ref>. 
 
The ensemble sampler was applied to each of these four simulated profiles.  Five variable parameters were investigated in each case, namely: the maximum velocity (v_ max), the minimum velocity (v_ min), the index of the power-law dust density distribution (β), the grain radius (a) and the total dust mass (M_ d).  These parameters were selected on the basis of previous models (B16, <cit.>) which suggested that they are the  parameters to which a simulated line profile is most sensitive in a spherically symmetric scenario. Whilst dust optical depth and albedo could be substituted for dust mass and grain radius, the latter parameters were chosen in order to allow for more straightforward comparison to other works which present results in these terms. Prior distributions were adopted on all parameters and are described in detail in Table <ref>. Uniform priors were adopted for the maximum and minimum velocities and for the index of the density distribution.  Uniform priors were appropriate for these parameters since I sought to assume minimal prior knowledge and the range of feasible values that these parameters could take was easily encompassed within an order of magnitude.  This was not the case for the dust mass and the grain radius however, both of which could take values within a range covering more than three orders of magnitude.  As a result, these parameters were investigated in log space and uniform priors were adopted for log M_ d and  log a. The range of the prior for each parameter was either physically motivated (e.g. the minimum velocity of the expanding ejecta cannot be negative) or was based on realistic values given the observed line profile (e.g. there is no flux detected redwards of 8000 km s^-1). The adopted priors were the same for each of the four simulated lines.

In each case, 100 walkers were used and the code was run to convergence, which was determined based on the autocorrelation time. In general, several thousand steps were required to approach convergence.  In all cases, the runs were allowed to continue for several autocorrelation times past this point. An acceptance fraction in the range [0.2,0.5] was required in all cases and in most cases the acceptance fraction was ∼ 0.3.  

Figures <ref> to <ref> illustrate the results of these models.  For each simulated line profile, a two-dimensional contour plot of the posterior probability distribution for each pairing of the variable parameters is presented. The contours on these plots represent 0.5σ, 1.0 σ, 1.5 σ and 2.0 σ. Additionally, one-dimensional histograms of the probability density distribution for a single parameter (marginalised over the other parameters) are also presented with the 16^ th, 50^ th and 84^ th quantiles indicated, encompassing the central two-thirds of the data. The known, `true' values that were used to generate the simulated line profiles are marked on these plots.  For the sake of comparison, the single best-fitting model was tracked throughout the sampling process and this is also marked on the plots. As expected, the best-fitting model line profiles were virtually identical to the simulated line profiles input into the simulation at the start and so are not presented here.

In most instances, the parameters can be  tightly constrained.  However, there are certain parameters which exhibit a broad posterior probability distribution indicating that the line profile is largely insensitive to variations in this parameter. Dependencies and correlations between the parameters can be observed for some of the parameters, for example the maximum velocity and the density profile (see Fig. <ref>) or the grain radius and the dust mass (see Fig. <ref>).  For the majority of cases, the true values lie very close to or inside the most likely (1σ) regions of the contour plots.  Where there are exceptions to this, these can be understood as an insensitivity to  a specific parameter on which another parameter is dependent.  In particular, where the dust grain radius cannot be determined from the line profile, the dust mass is likely also to be ill-constrained.  I discuss the reasons for, and implications of, these results in more detail in Section <ref>. 




 §.§ MCMC model of SN 1987A




SN 1987A is an extremely well-studied, nearby CCSN that remains critical to our understanding of the formation and evolution of dust in CCSNe. Spectra of SN 1987A have been taken regularly since its outburst on 23 February 1987 and it is an ideal candidate for line profile modelling with asymmetric optical line profiles exhibited from ∼650 d <cit.>.





I applied the ensemble sampler to the Hα and [O i] 6300, 6363 Å lines of SN 1987A at 714 d post-outburst.  A region of the optical spectrum that was obtained with the CTIO-1.5m telescope on 6^ th February 1989 and includes these lines is presented in Figure <ref> <cit.>.  The spectrum is available on the CTIO archives. This epoch was selected to revisit due to the high signal-to-noise ratio of the spectrum and the distinct separation of the two features, which is not as clear at later epochs.  Additionally, the relative lack of  contamination of these broad lines by narrow nebular emission makes this epoch particularly attractive.

Both the Hα and [O i] 6300, 6363 Å lines have been previously investigated by B16 using damocles, who used a systematic, manual approach to determine a best-fitting set of parameters for both clumped and smooth dust density distributions.  I sought to compare the best-fitting parameters that they inferred with the results generated by the automated ensemble sampler.  Additionally, I was interested to understand whether a more sophisticated model that investigates a significantly higher-dimensional variable parameter space could be explored by employing the ensemble sampler. A grid-based or manual approach would not be feasible for higher-dimensional models and this was a primary consideration in the implementation of a MCMC procedure.  I have therefore investigated two models for SN 1987A at 714 d post-outburst.  



The first is a 5-dimensional, smooth model that allows for direct comparison with previous results.  The Hα line is modelled with a spherically-symmetric, smooth shell distribution.  The power-law emissivity distribution is coupled to the dust density distribution as described in Section <ref> and the same dust properties were used (i.e. 100% amorphous carbon dust with optical constants from <cit.>). This scenario is the same as that adopted by B16 for their smooth model of the Hα line at 714 d.  A five-dimensional parameter space is explored. Uniform priors were adopted for the maximum velocity (v_ max), the minimum velocity (v_ min) and the index of the power-law dust density distribution (β), and log-uniform priors were adopted for the grain radius (a) and the total dust mass (M_ d).  Full details of the priors can be found in Table <ref>. The range of the priors  was kept as wide as possible in an effort to identify any additional maxima in the posterior distribution and therefore to obtain all possible solutions.  In all cases, the best-fitting parameter set identified by B16 lies within the prior range adopted here. 

The modelled profiles were convolved to the resolution of the spectrum (16 Å) before the likelihood was calculate and the region between 440 km s^-1 and 1400 km s^-1 was excluded from this calculation since it is contaminated by the unresolved, narrow, nebular [N ii] 6583Å emission.   The high signal-to-noise ratio of the spectrum resulted in a negligible observational error as determined by assessing the variance in a flat region of the spectrum. The height of the continuum is a potentially important factor in determining a number of the model properties.  However, a preliminary investigation that included the continuum height as a free parameter revealed an insensitivity to the continuum height and so it was fixed at 2.1 × 10^-14 erg cm^-2 s^-1 Å^-1. 


The results of this model are presented in Figure <ref> with the best-fitting parameters as identified by B16 marked on the probability distributions for comparison.  In all cases, the previous results lie within 1σ of the marginalised 1D probability distribution and within the 1.5σ contour of the 2D joint-probability distributions.  This suggests good agreement between the two approaches but, as can be seen, significantly more information is yielded from the full investigation.  For example, the results indicate that the steepness of the density distribution does not significantly affect the likelihood.  They also highlight the relative insensitivity of the dust mass to all parameters except the grain radius.  The predictably strong correlation between grain size and dust mass, as noted by B16, is clear. However, whilst the grain radius has a fairly well-constrained minimum at around 0.05 μm, it is not tightly constrained at larger grain sizes and, as such, constraining the dust mass is difficult without further information.

The second model for SN 1987A treats both the Hα and [O i] 6300, 6363 Å lines simultaneously for the first time.  A spherically-symmetric shell-based geometry is once again adopted but, in this more complex scenario, the dust is located entirely in clumps that are stochastically distributed throughout the shell according to a power-law distribution.  The clumps all have equal volume equivalent to a single, cubical grid cell in the simulation of width R_ out/25, roughly consistent with what might be expected from Rayleigh-Taylor instabilities in the ejecta. The total volume of the ejecta occupied by dust clumps is described by the filling factor which is varied between 0.1 and 0.7.  All species (dust, Hα and [O i]) extend to the same maximum velocity, but the minimum velocities for Hα and [O i] 6300, 6363 Å are separate, variable parameters.  The minimum dust velocity is coupled to that of [O i] since it seems likely that most dust formation is occurring in regions of high metallicity where the constituent ingredients of dust grains are available to condense. All three species follow separate power-law density distributions (with the emissivity following the square of the density as described in Section <ref>). Finally, the flux ratio between the 6300Å and 6363Å components of the [O i] doublet is also left as a free parameter.  Intrinsically, the flux ratio is fixed. However, whilst it is assumed that the gas is optically thin, it is possible that there remain some gas optical depth effects that could influence this ratio and it is therefore included for clarity. This yields a total of 10 variable parameters which are summarised, along with the adopted priors for each parameter, in Table <ref>.  

The ranges of the adopted priors were motivated by the previous results of B16 and the results from the previous smooth 5D model.  For certain parameters (the dust mass and grain radius), the range was restricted slightly relative to the 5D simulation, without significant loss of information, in order to speed up the calculation. Physical factors also dictated the adopted ranges, e.g. the flux ratio F_ 6300/F_ 6363, which was capped at 3.3 since the theoretical value for an optically thin medium is 3.1 <cit.>, and the filling factor f, which must clearly be in the range [0,1]. The likelihood was calculated as per Equation <ref> but, for these purposes, the [O i] 6300, 6363 Å was scaled to the same peak flux as the Hα line in order to ensure that both features were weighted equally.  This aside, the adopted procedure for this model was identical to that of model B.  

The results of this 10D simulation are presented in Figure <ref>. A significant quantity of information is contained in this figure but it is of particular interest to note that the majority of parameters have been constrained and follow a distribution with a single peak.  The probability distribution peaks at an extreme of the range in the cases of the flux ratio F_ 6300/F_ 6363, the filling factor f and the clump number density distribution which is specified by β_ clump.  The line profile is not highly sensitive to the density distribution of any species but the minimum and maximum velocities can be restricted to a relatively narrow range, regardless of the values of the other parameters.  Of most interest however, is the strongly-peaked marginalised 1D probability distribution for the grain radius suggesting a large grain radius of the order of ∼0.2 μm. This has allowed the dust mass to be similarly constrained with the marginalised probability distribution yielding a 1σ range spanning only one order of magnitude. The best-fitting parameter set is marked in Figure <ref> for comparison and the corresponding line profile is presented in Figure <ref> for the purposes of illustrating the goodness-of-fit.

I discuss these results further in the context of dust formation in SN 1987A and other CCSNe in Section <ref>.




§ DISCUSSION





 §.§ Theoretical models



Of primary interest in investigating this approach to modelling asymmetric line profiles is whether the Bayesian methodology offers additional insight or rigour in comparison to manual or grid-based frequentist fitting.







The approach was initially tested against four simulated example line profiles that had been generated using damocles and that exhibited different shapes and features (models A1 – A4, see Figure <ref>). The posterior probability distribution for each of these models is presented in Figures <ref> to <ref> with the true parameters used to generate the profiles and the best-fitting parameter set from the chain indicated on the plots. A good test that the Bayesian calculation is being performed correctly is that the best-fitting model from the MCMC chain is, in nearly all cases, in broad agreement with the true values.  The resulting posterior successfully characterises the likelihood of a given parameter over a range of values and exhibits a single-peaked marginalised 1D probability distribution for most parameters, as well as revealing dependencies and correlations between parameters. The known, true values generally lie within the 1σ contour as would be expected.  Below, I discuss a number of interesting results from the simulated line profile modelling.



  §.§.§ Maximum velocity, v_ max


The maximum velocity of an emitting species is generally inferred from a line profile as the velocity at which the flux on the blue side of the profile goes to zero.  In the case of models A2 and A4, this approach would yield a reasonable estimate of the maximum velocity in agreement with the posterior distribution, which tightly constrains the value of the maximum velocity in both cases. However, models A1 and A3 (see Figure <ref> and <ref>) illustrate that it is not necessarily straightforward to determine the maximum velocity by eye. The transition between the blue wing of the emission line and the continuum (at zero flux in this theoretical scenario) is smooth for this steep emissivity distribution which makes it difficult to determine the exact velocity at which the transition occurs.  This issue is exacerbated for real observations where noise further obscures the inflection point.  

The full Bayesian calculation illuminates the relative likelihoods of different maximum velocities, considered both in isolation from the other parameters (via the 1D marginalised posterior) and jointly with the other parameters.  It highlights, for example, the postive correlation between the maximum velocity and the emissivity distribution (i ∝ r^-2β). Similarly, for model A1 (Figure <ref>), the relationship between grain radius and the maximum velocity is clear.  This can be interpreted as due to the fact that, for amorphous carbon grains, the single-scattering albedo is a monotonically increasing function of grain radius for fixed λ = 6563Å. This results in a red scattering wing extending to higher velocities than those observed on the blue side.  This feature of the line profile could be approximated, for small grains, by adopting larger maximum velocities.  This relationship is clearly revealed by the Bayesian calculation whilst the results still prefer the correct grain radius of  ∼ 0.1μm and maximum velocity of approximately  ∼3250 - 3750 km s^-1.

Determining the maximum velocity accurately is particularly important since it determines the  size of the ejecta and therefore has a significant effect on the overall dust optical depth to which radiation is exposed for a given dust mass. 
The co-dependence of the maximum velocity with several of the other parameters is handled rigorously by the ensemble sampler and can be easily quantified and communicated via the posterior distribution.




  §.§.§ Minimum velocity, v_ min


The range of viable values for the minimum velocity can be very narrowly constrained in all optically thin cases (i.e. A1 – A3).  It is only when the dust becomes significantly optically thick that the minimum velocity becomes harder to determine. In practice, it is normally the case that the blue-shifted peak flux is coincident with the minimum velocity of the emitting ion.  In this case, an asymmetry is observed as a result of absorption in the central regions causing an intrinsically flat-topped, boxy profile (as is produced by an expanding shell) to peak sharply at the blue `corner' of the flat top.  A secondary peak coincident with minimum velocity on the red side is also a possibility <cit.>.  These peaks are important for determining the minimum velocity.  Where dust optical depths are high enough that the peak flux  shifts beyond the minimum velocity (see simulated profile A4, lower right panel of Figure <ref>), there is less information available in the profile to constrain the minimum velocity. The results for A4 (Figure <ref>) suggest that the minimum velocity cannot exceed ∼2000 km s^-1, presumably because the profile would become too wide, but yield similar likelihoods for all v_ min<2000 km s^-1 with a broad peak centered around ∼ 1600 km s^-1.



  §.§.§ Density distribution index, β


The steepness of the density distribution (and hence also the emissivity distribution) is not tightly bound for any of the simulated lines.  However, this is most noticeable in run A2, where there is only a little variation in the likelihood of β across the full range explored.  The width of an intrinsic flat-top profile at its peak is determined by the minimum velocity but the shape of the wings is determined (for homologous expansion) by the steepness of the power-law emissivity distribution.  Where this is steeper, the profile appears more concave in its wings.  In A3, only a small fraction of the profile is in the wings, with the majority of the width of the profile arising from the intrinsically flat-topped region. As a result, there is limited information in the profile to allow β to be determined. However, the other parameters can still be reasonably estimated by marginalising over β since they are not strongly dependent on it.

This is a good illustration of the fact that, under certain conditions, observed line profiles will not contain sufficient information to determine some of parameters of interest.  Even in this case, however, it may be possible to constrain the other parameters by marginalising over these less sensitive parameters. The full posterior distribution clarifies the sensitivity of the line profile to the variable parameters.



  §.§.§ Grain radius, a, and dust mass, M_ d


The strong correlation between grain radius and dust mass is recovered by the Bayesian approach. The absorption and scattering efficiencies of dust grains of any species depend strongly on the grain radius, as does the cross-sectional area available for interaction.  As a result, there is a strong relationship between the opacity and the grain radius, and hence also the required dust mass and the grain radius. 

By making a number of assumptions,  the relationship between dust mass and grain radius can be determined analytically for a simplified version of the scenarios modelled in A1 – A4. bf We can then compare the 2D marginalised likelihood distributions for dust mass and grain radius to this analytic relationship as a test of the Bayesian approach.
We consider the dust number density to be independent of radius ( which is not the case for models A1 – A4). The dust optical depth at a given wavelength is then τ_ d = Q_ extπ a^2  n_ d R, where Q_ ext is the extinction efficiency, n_ d is the dust number density and R is the distance to be traversed by the photon.  
We also have the total dust mass described by M_ d = Vn_ d4π a^3/3ρ_ g, where V is the volume of the ejecta, ρ_ g is the mass density of a dust grain and other parameters are as previously defined. 
We require a specific dust optical depth in order to reproduce a line profile.  If we additionally assume that the physical extent of the ejecta is fixed and that the dust is entirely absorbing (such that Q_ ext = Q_ abs, the absorption efficiency), then we can conclude that M_ d∝ a/Q_ abs.  
At λ = 6563Å, Hα, we can determine the exact correlation between grain radius a and a/Q_ abs using Mie theory.  The resulting relationship for amorphous carbon grains (presented in Figure <ref>) is echoed in the joint 2D posterior distribution of log a and log M_ d in all theoretical runs as expected (A1 – A4; see Figures <ref> to <ref>).

However, the 2D likelihood distributions marginalise over the other parameters and therefore deviate from the analytic relationship derived above. Deviations from this relationship are due to, for example, the polychromatic nature of the transported packets, a dust number density that is non-constant with radius, a significant scattering component to the extinction etc. Silicate grains would be expected to follow a different, more complicated relationship. 

The dependency between grain radius and dust mass has significant implications for determining the ejecta-condensed dust mass via line profile fitting and has been discussed in detail by B16 and <cit.>. By quantifying the posterior distribution, the exact relationship between these two parameters can be understood by marginalising over the other parameters of interest.  If further information can be obtained that would allow the grain radius to be constrained, then the joint probability distribution described by the posterior dictates the required dust mass for a given model.  This grain radius could be estimated from dust emission features or from other techniques such as SED fitting.  The approach used here could also be expanded to include multiple optical or NIR emission lines at a given epoch in order to exploit the wavelength dependence of dust extinction and hence constrain the dust grain radius simultaneously with the other parameters (see Section <ref>).

One further implication of the dependency of dust mass on grain radius is that, unless the grain radius can be constrained reasonably tightly, the 1D marginalised dust mass probability distribution will tend towards a specific peak.  This is because there is a wide range of dust grain radii that all yield similar values of a/Q_ abs.  Since the marginalised probability distribution is integrated over the whole prior range of all other parameters, a narrow band of dust masses will naturally be preferred. Care should be taken to ensure that the grain radius has been accurately constrained before inferring the dust mass from the posterior.

It is worth noting that, whilst emphasis is often placed on determining the mass of dust that has formed in the ejecta of CCSNe, and this is therefore naively the most interesting parameter, determining the dust grain radius is critically important in its  own right.  In order to determine how much dust CCSNe can eject into the ISM, we must not only understand how much is formed in the ejecta but also how much is destroyed by shocks, and in particular by the reverse shock that will inevitably pass back through the newly-formed dust <cit.>.  The rate of destruction of dust grains by sputtering in shocks is independent of the size of the grain  and, as such, the initial size of the dust grain is critical to understanding whether or not it will survive into the ISM or will eventually be destroyed <cit.>.




 §.§ Application to SN 1987A at 714 d



The need to isolate the grain radius motivated the production of two different models of SN 1987A, one significantly more complex that the other.  The initial smooth model in five dimensions (model B) couples the Hα emissivity distribution with the dust density distribution and is analogous to the models of Hα produced by B16.  Their results are indicated on the posterior distribution which is presented in Figure <ref>. They are generally in good agreement with the results produced by the ensemble sampler.  However, additional insight is gained into the range of viable values for the maximum and minimum velocities, with the most likely regions of parameter space leaning towards a slightly lower maximum velocity at 3000 km s^-1 (compared to the B16 estimate of ∼3250 km s^-1 ) and a slightly higher minimum velocity at ∼900 km s^-1 (compared to the B16 estimate of 813 km s^-1). The steepness of the emissivity distribution is not tightly established but does not affect the ability of the sampler to constrain the other parameters.  

Of most interest, however, is the estimation of the dust grain radius and the dust mass. As previously discussed, there is a strong correlation between these parameters.  Since there is a wide range of small dust grain radii that result in similar dust mass estimates, there is a peak in the marginalised dust mass probability distribution that suggests a dust mass of ∼10^-6 M_⊙.  We can infer that this is likely the case if only small dust grains are present in the ejecta.  However, the dust grain radius is not tightly constrained by the smooth fitting, with a wide range of values >0.15 μm yielding similar probabilities. 

Model C is significantly more detailed and includes additional variable parameters resulting in a 10-dimensional parameter space.  All of the dust is located in clumps. This is a more realistic dust distribution than a smooth radial power-law; dust has been observed to be located in clumpy or filamentary structures in a variety of different CCSNe and remnants <cit.>.  In addition to a higher-dimensional parameter space, the more complex model also treats both the Hα line and the [O i] 6300, 6363 Å doublet simultaneously.  By providing the sampler with more data, and in particular two lines separated in wavelength space, the grain radius can be reasonably constrained.  The median grain radius is ∼ 0.2 μm, with a maximum at 1σ of ∼ 1 μm.  This yields a dust mass that is constrained to within one order of magnitude, with a median dust mass of ∼ 4.5 × 10^-5 M_⊙ and a maximum dust mass at 1σ of ∼ 1.5 × 10^-4 M_⊙. These estimates are very similar to the separate Hα (5.5 × 10^-5 M_⊙) and [O i] 6300, 6363 Å (2.0 × 10^-4 M_⊙) estimates by B16 for a grain radius of 0.6μm.  However, they are somewhat lower than the dust mass estimates inferred from radiative transfer models of the SED of SN 1987A presented by <cit.> for this epoch, who deduce a dust mass of 1.0 × 10^-3 M_⊙ at 615 d post-outburst.  This  discrepancy may be a result of their adoption of an MRN dust grain radius distribution (n(a) ∝ a^-3.5 for 0.005<a<0.25, <cit.>) or the assumption here of a single grain size.  

This more complex model, which is clearly still a simplification of a highly complicated reality, yields considerable insight into the relative likelihoods of the velocity distributions of the different species and the mass and grain radius of the dust in the ejecta.  Additionally, however, we also gain insight into other properties of the geometry of the nebula at this epoch. The results suggest that the clumps are all likely concentrated towards the central regions (high β_ clump) and occupy only a small fraction of the total volume of the ejecta (low f).  Similarly, they indicate that the [O i] is also concentrated towards the central regions (median β_ [OI] of 2.59) with the hydrogen more diffusely distributed (median β_ Hα of 1.18). This suggests a geometry that would be consistent with observations of SN 1987A obtained by <cit.> using the Atacama Large Millimeter Array (ALMA).  These spatially-resolved observations of IR lines of CO and SiO reveal that both species are concentrated in the inner ejecta and occupy a clumpy distribution, suggesting that the heavier elements, and in particular oxygen, are likely located in these central regions.

The results would also be consistent with the structures and geometries predicted by hydrodynamic explosion models of CCSNe <cit.>.  These models predict that, at very early times, only a few seconds after the explosion, the heavier elements are  mostly located within the central regions of the ejecta with small clumps of fast-moving material escaping at higher velocities.  A more expansive, more diffuse hydrogen envelope is also present.  Once homologous expansion has set in, the geometry will remain self-similar for many hundreds of years, assuming that there is no encounter with significantly dense circumstellar material, and so it may not be unreasonable to compare these results.



  §.§.§ 


Supernovae and supernova remnants are highly complex objects.  I have not included in my models different dust species, nor dust grain size distributions, and I have also restricted my investigations to geometries that are, with the exception of a stochastically generated dust clump distribution in one case, spherically symmetric.  These are important factors that should be explored in future work.  Similarly, I have explored only a few particular models.  The results of these analyses do not make comment on the validity of the model itself, rather on the relative likelihoods of the parameters given that particular model.  Care should be taken in the future to assess the applicability of a given model and whether dust formation represents the most likely explanation for the observed properties of a given line profile.  This had already been established from previous work in this case (B16).  The application of a Bayesian procedure may prove useful in this regard also since it lends itself well to quantified model comparison.  However, the above results illustrate the overall power of this methodology to constrain parameters, identify the parameters to which the line profile is sensitive and characterise dependencies between the parameters.  Most importantly, I am able to investigate and analyse highly complex models for which manual parameter estimation would be extremely difficult. 




§ CONCLUSIONS



I have applied an affine invariant ensemble sampler to the Monte Carlo radiative transfer code damocles in order to explore the variable parameter space in a rigorous fashion and apply a Bayesian methodology to the inference of conclusions from the data.  I have utilised the algorithm presented by <cit.> and implemented by <cit.> in order to create a Fortran-Python hybrid code that is capable of fitting dust-affected optical and NIR line profiles from CCSNe at late stages in their evolution in order to construct a posterior probability distribution.

The code was applied to four different simulated line profiles that were generated by damocles in order to represent different sorts of dust-affected line profiles that are observed in the spectra of late-time dust-forming CCSNe.  A smoothly distributed, spherically symmetric geometry was adopted and five variable parameters investigated. The posterior distributions are in good agreement with the known, true parameters and suggest that the methodology is accurate and effective for parameter estimation.  The theoretical runs highlight a number of dependencies between specific parameters.  The power of the Bayesian inferential approach in revealing and quantifying these dependencies is beneficial for future research using this methodology, but also illustrates the need for care when using line profile fitting (or indeed any other method) to estimate model parameters from observations. 

I also revisited the Hα line and [O i] 6300, 6363 Å doublet of SN 1987A at 714 d. A simple model  with five variable parameters analogous to the smooth model of Hα investigated by B16 was initially adopted.  I also investigated a significantly more complex model in 10-dimensional parameter space that treated both Hα and [O i] 6300, 6363 Å simultaneously.  The dust mass and dust grain radius predictions  are in agreement with the previous manual approach but their relative likelihood is now quantified, as is their dependence on other parameters.  The affine invariant ensemble sampler has proved to be an efficient and effective method to investigate and analyse highly complex models for which manual parameter estimation would be extremely difficult. The Bayesian methodology allows for considerably more insight to be gained and communicated than the previous manual approach and there is significant potential for using this approach to determine accurate ejecta dust masses for a large number of CCSNe.





§ ACKNOWLEDGEMENTS


AB would like to thank Dr Boris Leistedt for his patient guidance and teaching on astrostatistics and Bayesian methodologies, as well as Dr Roger Wesson, Dr Ilse de Looze and Professor Mike Barlow for many discussions and their help in readying this paper for publication. Many thanks also to the anonymous referee for their helpful suggestions. 

This work was supported by European Research Council (ERC) Advanced Grant SNDUST 694520 and is based on publicly available observations from the archives of the CTIO.


mnras







