An automated approach for photometry and dust mass calculation of the Crab nebula

TL;DR

This paper presents an automated method using Herschel Space Observatory data to accurately estimate the dust mass in the Crab Nebula, highlighting the importance of precise photometric analysis.

Contribution

The study introduces an automated image processing approach for constructing the spectral energy distribution and estimating dust mass in supernova remnants.

Findings

Dust mass of 0.056±0.037 solar masses in the Crab Nebula.

Temperature of the dust is approximately 42 K.

Automated techniques improve object-background separation for accurate photometry.

Abstract

Ample evidence exists regarding supernovae being a major contributor to interstellar dust. In this work, the deepest far-infrared observations of the Crab Nebula are used to revisit the estimation of} the dust mass present in this supernova remnant. Images in filters between 70 and 500 $\mu$m taken by the PACS and SPIRE instruments on-board of the Herschel Space Observatory are used. With an automated approach we constructed the spectral energy distribution of the Crab nebula to recover the dust mass. This approach makes use of several image processing techniques (thresholding, morphological processes, contouring, etc..) to objectively separate the nebula from its surrounding background. After subtracting the non-thermal synchrotron component from the integrated fluxes, the spectral energy distribution is found to be best fitted using a single modified blackbody of temperature $T=42.06\pm1.14$ K and a dust mass of $M_{d}=0.056\pm0.037$ M$_{\odot}$. In this paper, we show the importance of the photometric analysis and spectral energy distribution construction in the inference of the dust mass of the Crab nebula.

Tables

Table 1. Table 1: Photometric observations of the Crab Nebula using the Herschel Space Observatory .

Obs. ID	Obs. Date	Obs. Duration (s)	Instrument	Filter set $\lambda$ ($\mu$m)
1342183905	2009-09-15	2221	PACS	(70, 160)
1342183906	2009-09-15	2221	PACS	(70, 160)
1342183907	2009-09-15	2221	PACS	(70, 160)
1342183908	2009-09-15	2221	PACS	(70, 160)
1342183909	2009-09-15	2221	PACS	(100, 160)
1342183910	2009-09-15	2221	PACS	(100, 160)
1342183911	2009-09-15	2221	PACS	(100, 160)
1342183912	2009-09-15	2221	PACS	(100, 160)
1342191181	2010-02-25	4555	SPIRE	(250, 350, 500)
1342204441	2010-09-13	1671	PACS	(70, 160)
1342204442	2010-09-13	1671	PACS	(100, 160)
1342204443	2010-09-13	1671	PACS	(70, 160)
1342204444	2010-09-13	1671	PACS	(100, 160)

Table 2. Table 2: Results of calculation for all the images. The first column lists the reference wavelength of each photometric band. It is clear that as the resolution of the image increases (and pixel size decreases to maintain adequate sampling), the number of number of pixels inside the ROI is increasing as well, however the area itself is quite constant, statistically it amounts to 116493.5 ± 2 650.5 ′′ ⁣ 2 plus-or-minus 116493.5 superscript 2650.5 ′′ 2 116493.5\pm 2\,650.5^{\prime\prime 2} , where the uncertainty is computed as 1/6 × 1 / 6 ​ ∑ i = 1 6 ( A R ​ O ​ I ​ i − A R ​ O ​ I ¯ ) 2 absent 1 6 superscript subscript 𝑖 1 6 superscript subscript 𝐴 𝑅 𝑂 𝐼 𝑖 ¯ subscript 𝐴 𝑅 𝑂 𝐼 2 \times\sqrt{1/6\sum_{i=1}^{6}(A_{ROIi}-\overline{A_{ROI}})^{2}} .

$\lambda$	$\alpha$	$N$	$A_{ROI}$	$S_{\nu }^{tot}$	$S_{\nu }^{back}$	$S_{\nu }$	$\Delta S_{\nu }^{ROI}$	$\Delta S_{\nu }^{back}$	$\Delta S_{\nu }^{str}$	$\Delta S_{\nu }$	${\theta }_{opt}$
($\mu$m)	(”)	(pix)	(”2)	(Jy)	(Jy)	(Jy)	(Jy)	(Jy)	(Jy)	(Jy)	(Jy)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)
70	1.6	53858	137876	219.509	4.491	215.017	1.380	0.055	0.336	1.771	-0.000302
100	1.6	39920	102195	217.895	8.805	209.090	1.165	0.049	1.467	2.681	0.000303
160	3.2	9660	98918	140.946	9.940	131.007	1.111	0.090	8.902	10.104	0.00275
250	6	2910	104760	162.376	69.860	92.516	1.246	0.106	1.257	3.347	0.0277
350	10	1186	118600	144.968	48.730	96.238	2.224	0.104	1.019	5.059	0.0476
500	14	697	136612	127.976	24.716	103.260	3.448	0.069	1.538	8.433	0.0450

Table 3. Table 3: Flux densities used to model the synchrotron SED using data from Gomez et al. ( 2012 ) . We have specifically avoided using bands located too close to the infrared excess.

$\lambda$($\mu$m)	$S_{\nu }^{syn}$(Jy)	$\Delta S_{\nu }^{syn}$(Jy)
3.4	12.9	0.6
4.5	14.4	0.26
4.6	14.7	0.75
5.8	16.8	0.1
8	18.3	0.13
850	128.6	3.1
1382	147.2	3.1
2098	187.1	2
3000	225.4	1.1
4286	253.6	2.5
6818	291.6	1.3
10000	348.2	1.2
10000	348.2	1.2

Table 4. Table 4: This table summarizes that flux density calculation and the relevant error bars.

$\lambda$($\mu$m)	$S_{\nu }$(Jy)	$\Delta S_{\nu }$(Jy)	$S_{\nu }^{syn}$(Jy)	$\Delta S_{\nu }^{syn}$(Jy)	$S_{\nu }^{IR}$ (Jy)	$\Delta S_{\nu }^{IR}$(Jy)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
70	215.02	1.77	45.86	1.26	169.16	3.34
100	209.09	2.68	53.08	1.49	156.01	4.08
160	131.01	10.10	64.36	1.86	66.65	11.28
250	92.52	2.61	77.27	2.29	15.24	4.39
350	96.24	3.35	88.70	2.68	7.54	6.09
500	103.26	5.06	102.66	3.17	0.60	9.01

Table 5. Table 5: The characteristics of the dust based on the single component modified black-body model.

Parameter	value $\pm$ $\Delta$ value
$T$	42.06 $\pm$ 1.14 (K)
$M_{d,\Delta D=+0.5Kpc}$	0.084 $\pm$ 0.010 (M⊙)
$M_{d,\Delta D=-0.5Kpc}$	0.030 $\pm$ 0.004 (M⊙)
$M_{d,\Delta D=0Kpc}$	0.054 $\pm$ 0.007 (M⊙)
$\overline{M_d}$	0.056 $\pm$ 0.034 (M⊙)
$\overline{M_{d,cal=5\%}}$	0.056 $\pm$ 0.037 (M⊙)

Equations

$$S_{\nu }^{tot}=f(\theta ).$$

$$S_{\nu }^{tot{}^{\prime }}\equiv {\delta }_h[S_{\nu }^{tot}](\theta )=S_{\nu }^{tot}(\theta +\frac{1}{2}h)-S_{\nu }^{tot}(\theta -\frac{1}{2}h),$$

$$S_{\nu }^{tot{}^{\prime }}\equiv {\Delta }_h[S_{\nu }^{tot}](\theta )=S_{\nu }^{tot}(\theta +h)-S_{\nu }^{tot}(\theta ),$$

$$A_{ROI}=N.{\alpha }^2.$$

$$S_{\nu }^{tot}=\sum _{i=1}^{N}I(x_i,y_i{)}_{ROI}.$$

$$S_{\nu }^{back}=\sum _{i=1}^{N}I(x_i,y_i{)}_{ROI}^{back}.$$

$$\sigma_{ROI}=\sqrt{1/N\sum_{i=1}^{N}\bigg{(}I(x_{i},y_{i})_{ROI}-\overline{I_{ROI}}\bigg{)}^{2}}$$

$$\sigma_{back}=\sqrt{1/M\sum_{i=1}^{M}\bigg{(}[I(x_{i},y_{i})-I(x_{i},y_{i})^{back}]-\overline{[I-I^{back}]}\bigg{)}^{2}}$$

$$\Delta S_{\nu }=({\sigma }_{ROI}+{\sigma }_{back}+{\sigma }_{str})\sqrt{N}$$

$$S_{\nu }^{syn}=A{\lambda }^{\gamma }.$$

$$\chi^{2}_{w}=\sum_{i=1}^{N}\Bigg{(}\frac{S^{syn}_{\lambda_{i}}-y_{i}}{\sigma_{i}}\Bigg{)}^{2},$$

$$S_{\nu }^{IR}=S_{\nu }^{tot}-S_{\nu }^{back}-S_{\nu }^{syn}.$$

$$B(\nu ,T)=\frac{2h{\nu }^3}{c^2}\frac{1}{e^{\frac{h\nu }{k_{B}T}}-1},$$

$$S_{\nu }^{IR}=\Omega N{\kappa }_0{\left(\frac{\nu }{{\nu }_0}\right)}^{\beta }B(\nu ,T),$$

$$S_{\nu }^{IR}\propto {\nu }^{\beta }B(\nu ,T).$$

$${\kappa }_{\nu }=K_{3000GHz}(\frac{\nu }{3000\mathrm{GHz}}{)}^{\beta },$$

$$S_{\nu }^{IR}=\frac{M_d{\kappa }_{\nu }}{D^2}B(\nu ,T),$$

$$S_{\nu}^{IR}=\frac{8\times 10^{26}}{D^{2}}\frac{h\nu^{3}}{c^{2}}\bigg{(}\frac{\nu}{3\times 10^{12}Hz}\bigg{)}^{1.5}\,\Bigg{[}\frac{M_{d}}{e^{\frac{h\nu}{k_{B}T}}-1}\Bigg{]}.$$

Full text

††journal: Experimental Astronomy

1 Introduction

Stars in their Asymptotic Giant Branch (AGB) phase have been recognized as important sources of dust (with typical yields of $0.01-0.05$ M*⊙*, Helling & Woitke, 2006; Andersen, 2007; Gail et al., 2009; Schneider et al., 2014; Dell’Agli et al., 2015). But the amount of dust produced and ejected into the interstellar medium by these stars is not enough to compensate for the known destruction rates (Zhukovska et al., 2008). Stellar dust injection by itself also cannot explain the amount of dust in high-redshift objects where star formation has not been going on for large periods of time for a significant AGB star population to exist. Hence the dust mass calculated from observations must have other sources that can compensate this discrepancy.

In the literature, evidence has been presented which favors the concept of supernovae being a major contributor for dust formation (Krafton et al., 2016). This is supported by nucleation models (Silvia et al., 2010; Marassi et al., 2015), where dust is formed in dense regions during the first phases of supernovae explosions and allows explaining the large quantities of dust found in galaxies (Dwek & Cherchneff, 2011). Theoretical models predict that supernovae can be efficient dust producers ($0.1-1.0$ M*⊙*), potentially accounting for most of the dust production in the early Universe (De Looze et al., 2017). Observational evidence for this dust production efficiency is however currently limited to only few supernova remnants (SNR).

To revisit this issue, the Crab nebula presents many advantages: as one of the most commonly observed celestial objects, it possesses a wealth of observational data and at a relatively young age of the remnant, the mass of swept-up material is small compared to the mass in the supernova ejecta, which makes it possible to separate supernova dust from any swept-up circumstellar material. It also provides the cleanest view along the line of sight compared to other observed SNRs. It is a pulsar wind nebula (PWN) that is sweeping up supernova ejecta material observed in the form of bright optical filaments (Gaensler & Slane, 2006). Unlike other remnants, there is almost no dusty galactic material in front of, or behind the Crab Nebula, so the images are less effected by scattering or extinction (though we will see that the environment is far from innocuous). Spitzer infrared observations at short wavelengths allowed the discovery of warm dust in the filaments but it missed a potentially more massive cool dust component. Herschel, observing at longer wavelengths, is better suited for that purpose, allowing astronomers to measure the total amount of dust for the first time.

Trimble (1977) reported the first evidence for the existence of dust in the Crab Nebula. This emerged from the existence of infrared (IR) excess emission in the synchrotron spectrum. In further studies, discrepancies were noted in the dust masses based on different theoretical modelling assumptions. For example, in Gomez et al. (2012), the dust component observed along the ionized filament structure in spatially coincidence with the location of ejected material was modelled with modified black-bodies (MBB). A fit with a single MBB model lead to the derivation of 0.08 M*⊙* to 0.14 M*⊙* depending on the choice for dust composition (respectively amorphous carbon or silicate). This model however provided a poor fit to the entire spectral energy distribution (SED). A more adequate fit was obtained with two temperature components, leading to masses of 0.11 M*⊙* for amorphous carbon grains to 0.24 M*⊙* for standard silicates (we note that, in that paper, the 24 $\mu$m flux was included in the fitting process, which is debatable for a model assuming the emission is coming from a black-body). In Temim & Dwek (2013), the total dust mass derived from the best-fit models to the observed IR spectrum ranged between 0.019 M*⊙* and 0.13 M*⊙. In their work, the models compared three different dust compositions, each characterized by a power-law distribution in grain radii whose index was allowed to vary. In Owen & Barlow (2015) several models were applied and compared with a variety of grain dust compositions and grain size power law spectral indexes ranging from 2.7 to 3.5, resulting in dust mass values ranging between 0.1 to 1.0 M⊙*.

It is thus clear that the dust content of the Crab nebula is still not well constrained in the literature, hence we revisit in this work the estimation of the dust content in the nebula using improved photometric images and image analysis techniques, and attempting at a proper quantification of all the uncertainties associated to the process. We have applied an automated image analysis technique to optimize the extraction of the integrated object flux density at each wavelength. This approach increases the accuracy of the flux densities by excluding as much as possible the background flux contamination by edge detection processes. This paper is arranged as follows: in Section 2, we will describe the data and processing used to produce the improved images, in Section 3 image analysis and flux calculation will be explained, Section 4 will present the model used for the dust mass calculation, Section 5 will describe the calculations of the mass and temperature of dust, and finally in Section 6 we discuss our result and conclude this study.

2 Observations and Data

Images used for the calculation of the flux densities of the Crab nebula were taken by the Herschel Space Observatory (Pilbratt et al., 2010). The Crab nebula was observed between September 2009 and September 2010. The PACS (Poglitsch et al., 2010) and SPIRE (Griffin et al., 2010) instruments performed photometry at 70, 100, 160, 250, 350 and 500 $\mu$m, as part both of a calibration program for the PACS observing modes, and of a Principal Investigator observing program (see Gomez et al., 2012). A summary of the observations is listed in Table 1.

The PACS photometry data were obtained in scan-map mode with speed 20 arcsec/s. In order to obtain images in all three filters of the PACS photometer, observations are done by pairs, one using the 70 and 160 $\mu$m filters, and one using the 100 and 160 $\mu$m filters. Therefore the 160 $\mu$m image has twice the exposure time as the other two channels. In the present work we have for the first time included data that was obtained on the Crab nebula during the testing and qualification of PACS observing modes. As these data were obtained in an instrumental set-up which is identical to the operational one, they can be directly combined with the already published data, thus nearly doubling the depth of the resulting maps (at 70, 100, and 160 $\mu$m). The resulting PACS maps have a size of $25^{\prime}\times 25^{\prime}$, although the full depth of the maps is only achieved on a smaller field of view as not all observations mapped the same area. The SPIRE maps were obtained using the Large Map mode with a scan length of 30 arcmin at a speed of 30 arcsec/s, to produce a map of $32^{\prime}\times 32^{\prime}$.

Figure 1 displays the resulting map at 160 $\mu$m as it is quite representative of what we observed in the field. In almost all our images, there exists two objects: the Crab nebula and an unknown celestial object to its eastern side, as well as an extended background/foreground emission to the NW that reaches to the nebula itself. This unknown object has a significant IR flux especially in the long wavelengths bands and will need to be masked out when we will try to fit the structure of the extended background.

The PACS data were first reduced using the Herschel Interactive Processing Environment (HIPE, see Ott, 2010). Then the Unimap software (Piazzo et al., 2015), which is a generalized least-squares map maker for Herschel data, was used to eliminate low frequency noise in the maps and combine all PACS observations in a single map per photometer band. For SPIRE, the standard pipeline HIPE version 14.0.0 was used (Griffin et al., 2010).

3 Photometric Image analysis

**The main challenge regarding the photometry of extended objects in the far infrared (FIR) is the segmentation and the separation of the object of interest in foreground from that of background emission (note that we use these terms, foreground and background, in a generic way as some of the structure unrelated to the Crab Nebula could be in front of it). A proper estimation and removal of non-intrinsic contributions to the object of interest is very important. **The Crab nebula observations are no exception: Fig. 2 shows significant emission in the immediate vicinity of the nebula (this is clearly visible from 160 to 350 $\mu$m but this emission is detected at nearly all wavelength).

To calculate the integrated flux densities of the Crab nebula at each wavelength, the standard approach is applying aperture photometry with an aperture based on the shape of the object under study (such as circle or an ellipse) but for extended objects this means including higher number of background pixels which can create a substantial error by accumulation. Therefore it is necessary to detect the edges of the Crab nebula morphologically and as objectively as possible. The iterative procedure that we present here aims at detecting the outer edges of the Crab nebula, through the adaptation of image analysis techniques for edge detection. It also aims at increasing the accuracy of the integrated flux density by reducing the inclusion of background pixels.

Independently for each photometric band, our adopted edge detection technique is based on computing a threshold flux density value that will classify the pixels of the photometric image as either belonging to the object of interest or background. In other words, the algorithm searches for the transition between pixels dominated by the foreground object (Pulsar synchrotron and dust IR emission) and that of the background (Galactic sources) by searching for the location in the image where variations in pixel brightness are tracing the wings of the PSF applied to the extended object. The successive steps of the method we have retained are detailed in the following subsections.

3.1 Initial choice of an image thresholding level

An initial test was done using a classical image segmentation technique called Otsu’s method (Otsu, 1979). It utilizes the statistical analysis of the pixel intensity histogram of images to classify the pixels into two categories or distributions(background and foreground). Otsu’s method is an algorithm that tries to find a separation criteria or threshold value for the flux densities of the image assuming that it has a bimodal nature. Using this assumption, it derives a flux threshold value that maximizes the intraclass variance and minimizes the interclass variance of the data to generate a clear separation of the image in two sets of areas. Once this optimal threshold is determined, a mask is built from the map where pixels above the threshold are set to 1 (foreground) and pixels below are set to 0 (background), see Fig. 3 (a) and (b). Most astronomical images do not have a bimodal distribution, however, and situations of uni or multi-modality can exist. **The histograms of Herschel maps have a common nature (see Fig. 3a): the low background flux density values dominate with the highest occurrence rates, whereas the higher flux densities associated to the filaments or the pulsar synchrotron provide a long tail of less populated bins at high intensity values.

**

For this reason, the region initially selected as foreground using Otsu’s threshold ($\theta_{otsu}$) is not optimal for deriving the photometry of the Crab Nebula: it is creating voids in the inner parts of the nebula and rejecting outside of the mask fainter emission regions that still obviously belong to the Nebula (see Fig. 3). Thus $\theta_{otsu}$ is rather an estimation for the initial guess needed for the adopted iterative optimal threshold determination method that we will describe in section 3.3.

3.2 Opening and isophotal contouring

To remove the effect of high spatial frequency noise which manifests itself in the form of detached or attached fine structures around the object of interest in the masked image, we applied a morphological image processing step called opening (Smith, 1997, see Fig. 3, b and c). After this step, we measured the length of all the contours that can be built on the binary masks, and we used the longest of them to define the region of interest (ROI) for photometry: all pixels with coordinates falling inside the contour (including the hollow regions excluded by Otsu’s threshold, see Fig. 3c) contribute to the integrated flux associated to this value of the threshold. To each threshold value an ROI and thus a flux are associated, thus creating a functional relation from which the optimum threshold will be computed. The opening step was necessary only for the PACS as the spatial noise is higher than in the SPIRE images. For wavelengths 70 and 100 $\mu$m, we used a circular morphological filters, also known as structural elements (strels), with a radius of 12 and 5 pixels respectively, and for 160 $\mu$m we used also circular strels but with a radius of 8 pixels. The selection of the strels is essentially empirical. Opening helped decreasing the number of contours that can be build on the binary mask at each step and, as a result, decreased the memory consumption and increased the speed of the process.

3.3 Finding the optimum threshold and integrated flux density

As the Otsu-derived threshold always clearly cuts into the object, iterations on the threshold consist in exploring decreasing values of it. During the descent of the threshold values, the generated ROIs will start to include a larger region of the Crab nebula and will include the outer regions where the flux density values tend to have a very slow and linear variation. The aim of this iterative approach is to objectively detect the edge of the nebula through characteristic features of the gradient of the integrated flux density. To this aim, we repeat the previously described process of ROI pixel coordinate generation for every threshold values, generating contours such as those displayed on Fig. 4. As a result, we can build a piecewise functional relationship between threshold values ($\theta$) and the integrated flux densities and can be defined as:

[TABLE]

The first degree numerical derivatives (gradient) is calculated because it is considered as a more reliable indicators to the optimum threshold: in some cases the $S_{\nu}^{tot}$ by itself would be sufficient to locate the optimum threshold, but in other cases the values were difficult to locate because of the nature of the variability near the edges of the object. For all the wavelengths however, the first derivatives did reveal sudden variations at the transition between foreground and background (see below). If $h$ is the step of $\theta$, then the numerical derivatives of $S_{\nu}^{tot}$ can be computed using the central differences in the interior:

[TABLE]

and the first differences at boundaries:

[TABLE]

While decreasing the threshold values, the generated isophotal contours will start reaching the boundaries of the Crab nebula. Near the edges, these contours will start including pixel values which are lower (even negative) than that of the inner region of the Crab nebula. As a result, and in our tests, we realized that the $S_{\nu}^{tot}$ of ROI near the boundary with respect to $\theta$ will either have a linear shape or a sudden variation. When calculating the first order gradient, a constant value will be expected and we will see a plateau with a low signal to noise ratio followed by a sudden variation. The sudden variation following the plateau in the first derivative acts as the main indicator (see Fig. 5) for our optimum threshold.

To select the optimum threshold, the standard deviation of the neighborhood of each point of the first derivative was calculated. In our tests, a neighborhood of 6 step points and a standard deviation limit of 0.08 was used as a criteria for a threshold value to be considered belonging to the plateau. The smallest $\theta$ which passed this criteria was selected as the optimum threshold ($\theta_{opt}$). A posteriori visual inspection of the edges determined by this process revealed no issue that would lead to photometric errors, we therefore adopted it as our objective method to measure the Crab nebula’s integrated flux density at all wavelengths (See Appendix for results and graphs).

3.4 Derivation of the Nebula’s integrated flux

With the object’s ROI defined, the last step to determine flux of the object is to remove the background contribution in the ROI. We used the polygon generated by $\theta_{opt}$ at each wavelength and we used its interior pixel coordinates to mask the ROI. Then, a global surface interpolation is applied using a degree 2 polynomial on the pixels outside the ROI, to represent the large spatial frequencies in the background structure, and the result of this interpolation is used to compute a model of the background at the same sampling and resolution as the original map (Fig. 6). **for clarity, this model image is referred to as $I_{back}$.

**

A number of quantities attached to the ROIs can then be derived. To establish a mathematical formalism, we assume that the $(x_{i},y_{i})$ are the coordinates of pixels, $I(x_{i},y_{i})$ is the observed flux density at each pixel, $\overline{I}$ is their average, $\alpha$ is the angular size of each pixel, N is the number of the pixels in the ROI. Hence we can define:

The ROI area, to compare the sizes of the ROIs between different wavelengths:

[TABLE]

The total flux density in the optimum ROI:

[TABLE]

The predicted background flux density in the optimum ROI:

[TABLE]

The dispersion of the pixel fluxes in the optimum ROI:

[TABLE]

The dispersion of the pixel fluxes outside the optimum ROI which traces the existence of small spatial structures superimposed on the large one fitted by the polynomial (the background area is composed of M pixels):

[TABLE]

Further, since a polynomial is a crude representation of the complex background we are facing, we estimated a supplementary source of uncertainty: the area over which the polynomial fit is performed will impact the predicted level of the background in the center of the map, simply because a larger area potentially means a stronger deviation from our idealized polynomial representation of the background. We explored the variation of the predicted background contribution in the ROI ($S_{\nu}^{back}$) using areas that sample the map size beyond $\approx$3 times the object’s size (we consider smaller sizes as being too close to the object to provide a fair representation of the complexity of inferring the background level ). From this we derived an estimated $\sigma_{str}$ which is an error component related to the large scale structures observed in the map. As can be seen in 2, the largest contribution is at 160 $\mu$m, where the edges of the map also suffer from data processing artefacts.

Finally, the combined error on the flux is derived as:

[TABLE]

We note that in principle, calibration errors enter the error budget. However calibration errors are systematic, and therefore apply simultaneously to all bands in an unknow similar way. They should therefore not be mixed with the statistical errors studied above. These error amount to 5% in absolute photometry111see Herschel explanatory documentation on https://www.cosmos.esa.int/web/herschel/pacs-overview and https://www.cosmos.esa.int/web/herschel/spire-overview and in effect contribute uncertainty at the mass determination level.

The results of ROI calculations of the 6 images of PACS and SPIRE for wavelengths 70 to 500 $\mu$m are shown in Table 2. Before moving on to the derivation of the infrared excess attributed to dust, we point out that our derived fluxes are compatible, within the respective error bars, with the previous measurements of Gomez et al. (2012), with the notable exception of the 500 $\mu$m flux. We also remark that our derivation is systematically lower. However we note that our flux derivation is (a) purely data driven, i.e. we let the data indicate in an objective manner how the object of interest should be extracted from the map, and (b) even though we proceed for each band independently we obtain an area for the object of interest which is quite constant (and slightly larger at 500 $\mu$m, thus indicating that our lower flux value is not due to a smaller aperture). This will obviously have an impact when computing the dust mass in the nebula. **To further increase our confidence in our photometric process, we have explored how the background-corrected object flux vary if we artificially increase the ROI size by an amount similar to the observed dispersion of sizes as a function of wavelength (column 4 in Table 2), or greater as, for the longer wavelengths, expanding the perimeter by just one pixel already leads to a larger area increase than the dispersion. The observed flux variations are significantly smaller that the flux uncertainties listed in the table. Therefore the observed area variations from one band to another have no impact on the derived fluxes, and we postulate that differences with previous measurements indicate that the latter have been contaminated by foreground/background emission at the longest wavelengths.

**

3.5 Constructing the SED

After calculating the integrated flux densities, we continue constructing the SED of the Crab nebula. Due to the pulsar, the SED is dominated by the synchrotron radiation, on which is superposed an IR flux excess corresponding to the dust emission that we wish to model, to quantify the mass of dust in the remnant. Therefore it is necessary analyze the synchrotron component and remove its contribution. The synchrotron flux density model can be defined as:

[TABLE]

Thus to achieve this step, we applied a power law fit by minimizing a merit function ($\chi^{2}$) with respect to the spectral index ($\gamma$) and the amplitude of flux density ($A$). To account for the measurement errors, the merit function $\chi^{2}$ was weighted to handle the uncertainties. The uncertainty of the flux density at each wavelength was considered as the error weight ($\sigma_{i}$). Thus the weighted merit function becomes :

[TABLE]

where $S^{syn}_{\lambda_{i}}$ is the synchrotron flux density model evaluated at the $i^{th}$ band, $y_{i}$ is the observed data point. The data used for modeling the synchrotron radiation was adopted from Gomez et al. (2012), excluding the Herschel bands to avoid contaminating the fit by dust emission (see Table 3 and Fig. 7).** We have also acknowledged the existence of the spectral break value derived by Gomez et al. (2012) at 10 000 $\mu$m, and applied the fitting up to that point. **

The power law fit resulted in a spectral index $\gamma=0.410\pm 0.003$ and an amplitude equal to $A=8.039\pm 0.195$ Jy.$\mu$m*-γ*.

After calculating the power law fit, we predict the synchrotron emission flux densities and their uncertainties for PACS and SPIRE wavelengths (see Table 4, columns 4 and 5). To calculate the excess IR, we subtract the background and synchrotron contributions from the observed total flux:

[TABLE]

The results of the flux density calculations are presented in Table 4 (columns 6 and 7). As hinted before, we measure no excess due to dust at 500 $\mu$m.

4 Dust grain model

The Crab Pulsar (B0531+21) is the brightest source of synchrotron radiation in our Galaxy. Its emission in the NIR and MIR is dominated by continuum synchrotron radiation. The dust emission appears as an excess superimposed over that synchrotron continuum, with dust grains emitting, to a first approximation, a modified back-body spectrum (hereafter MBB, the black-body spectrum modified by the spectral dependence of the dust absorption cross section). For that emission, the basic parameters necessary to describe the emission from dust grains are the temperature Td and a form of the emissivity function $\epsilon_{\nu}$.

At FIR wavelengths ($\lambda\geq 60$ $\mu$m), the shape of the SED due to thermal emission from dust is empirically found to be well represented by the Planck function:

[TABLE]

which is evaluated at the dust temperature T, modified by a power law frequency dependent $\nu^{\beta}$ (Casey, 2012). More precisely the dust emission flux density can be written as:

[TABLE]

where $\Omega$ is the solid angle of the observing beam, N is the column density of emitting material, $B(\nu,T)$ is the Planck function, and $\kappa_{0}\left(\frac{\nu}{\nu_{0}}\right)^{\beta}$ is the opacity of the emitting dust. Note that $\tau_{0}=N\kappa_{0}$ is the optical depth at frequency $\nu_{0}$. An essential assumption in Equation (14) is that dust is optically thin, so that $\tau(\nu)\ll 1$ (Kelly et al., 2012), and hence all the dust grains within the beam contribute emission. Dropping the constants from this equation gives:

[TABLE]

The spectral index $\beta$ determines the opacity $\kappa_{\nu}$ of the dust, and encodes information about grain composition. It is found to be in the range 1 to 2 for the silicate and/or carbonaceous grain composition common in the diffuse ISM (Draine & Lee, 1984). Drawing from Hildebrand (1983), the dust mass absorption coefficient is:

[TABLE]

where for $\beta$ = 1.5, $\kappa_{3000GHz}=40$ cm2 g-1.

**

The IR excess attributed to dust in the Crab nebula can adequately be modeled by a single MBB. From this component we deduce from the SED fit a temperature and the dust mass $M_{d}$. Starting from equation (14), and relating the solid angle and column density to the distance of the object and total mass of grains (see for instance Laor & Draine, 1993; Weingartner & Draine, 2001) we get:

[TABLE]

where, $S_{\nu}^{IR}$ is the flux density, $\kappa_{\nu}$ is the dust mass absorption coefficient, and $D$ is the distance between to the Crab Nebula and is equal to $2.0\pm 0.5$ Kpc (Trimble, 1973). Combining equations (13), (16), and (17) we will get:

[TABLE]

Where now the equation is rearranged to be in Jy. In the modeling the optimization process will be based on the main two parameters which are the mass and the temperature. Statistically speaking, while optimizing the $\chi^{2}$, the predictor will be the frequency ($\nu$ in Hz) and the response will be the integrated flux density ($S_{\nu}^{IR}$ in Jy)

5 IR excess SED fitting and dust mass calculation

Although the modeling is based on $\nu$, in Figs. 9 and  8 the plots of the SED of the Crab nebula and the results from the fitting process are presented in a more traditional way as a function of the wavelength ($\lambda$). In both figures, we can see the fitting of the IR excess, in addition to a best fit to which the synchrotron flux is added. The results for the dust mass computation are presented in Table  5. We used the model-fitting open-source scipy package curve_fit() which applied fitting by including the uncertainties of the values of $S_{\nu}$222The description of the function with the related computation and references are explained in the documentation of the package: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html ..

In our modeling, we assumed the distance to the Crab nebula to be known, however to provide estimates of the mass, we calculated it using an upper and lower boundary values for the distance, since for a given SED, the mass and the square of the distance are proportional to each other. We calculated mass average of the mass $\overline{M_{d}}=0.056$ M*⊙*  and the uncertainty as the distance from the mean to the boundary values. The total uncertainty on the mass is sum of fitting and calibration errors (5%), or 0.037 M*⊙*.

As can be seen on Figure 8, the fit is poorer at 250 $\mu$m. We see no indication in our photometric process that this particular flux derivation would be incorrect. We have investigated whether a different assumption on the dust properties, in particular the exponent $\beta$ of the emissivity could provide a better fit to the SED, and indeed a better fit is obtained with $\beta=2$. While this can be instructive, we still keep to our initial selection of a value of 1.5 as it allows direct comparison with previous work, and is in effect tied to the absolute value of the mass absorption coefficient that we have used.

6 Discussion and conclusion

We have presented an improved method for calculating the flux densities using several image processing techniques applied on PACS and SPIRE images. With an optimal object-threshold selection method we constructed the spectral energy distribution of the Crab nebula in the Herschel band. We modeled the synchrotron component and removed its contribution to extract the IR excess. The resulting SED can adequately be modeled with a single component MBB using weighted least square fitting. **This leads to a dust component of parameters $T=42.06\pm 1.14$ K and $\overline{M_{d,cal=5\%}}=0.056\pm 0.037$ M⊙. The combined (synchrotron+dust) SED is shown in Fig 9 and  8 where we see that it under-predicts the flux in the range around 20 $\mu$m, i.e. it does not violate the constraint set by the Spitzer observation at 24 $\mu$m ($S_{\nu}=59.8\pm 6.0$ Jy; adopted from Gomez et al., 2012). As mentioned earlier, we consider that using a thermal process in this wavelength range is not justified, and remark that thermally fluctuating dust particles (responsible for emission in that range) will anyway not affect much the total mass in the nebula.

Our objective in this paper is primarily to estimate how the photometric methods impact the dust mass determination. We thus now compare our results with previous determinations. Gomez et al. (2012) and Temim & Dwek (2013) have provided reference derivation of the dust content of the nebula, where the main differences lies in the dust model that was used to convert from emission to mass. Our work is most comparable to that of Gomez et al. (2012) as Temim & Dwek (2013) used a more complex approach involving a simplified radiation transfer approach that led to temperature distributions for the dust component. While our value of the dust mass is comparable to the carbon-based single component model of Gomez et al. (2012), where the temperature is 40 K and the mass is 0.08 M*⊙*, there is a noticeable difference in the fact that the single component model is not the favored result of Gomez et al. (2012). Rather they favor a two-component model with temperatures bracketing our best-fitting temperature, and as a result a total dust mass 2 to 4 times larger due to the presence of cold dust. Two facts explain this difference: (1) we chose not to include observations below 70 $\mu$m in the fit for the reasons explained above, and (2) our dust excess SED is significantly different from that of Gomez et al. (2012), and in particular, with no excess at 500 $\mu$m, leaves no room for a cold component. **This stresses how complex and uncertain the process of determining the dust content in the Nebula is: cold dust easily becomes the dominant component in mass but rests mostly on photometry at long wavelength. This is a range where the Galaxy is still quite bright and therefore where background estimation and subtraction becomes critical. As can be seen, the background in the vicinity of the Crab is complex, asymmetric with respect to the object, such that determining it in an aperture around the target can lead to errors or underestimation of the uncertainties. Our approach attempts at capturing all the uncertainties involved in the photometry in presence of complex structural perturbations.

**

We thus draw the attention to an often overlooked aspect of the difficulty in measuring dust masses in resolved objects, and that is the construction of the spectral energy distribution of the object. This is even more complex in regions where the environment (background or foreground) of the object contributes emission at the location of the object. In this paper we have taken special care to the photometric approach of separating morphologically the background from the Crab nebula. Our proposed technique is a fully objective one which utilizes the gradient variation for choosing the optimum ROI and enabled to construct an objective SED. This is likely safer, and at least more reproducible, than techniques that rely on apertures, however carefully chosen, that appear to match the object, and it allows the adoption of a data driven-approach for the mass and temperature inference. We consider further that, after having applied this method independently on 6 different photometric images resulting from 2 different instrument (PACS and SPIRE), i.e. a data set where the point spread function, spatial sampling, and noise properties change quite significantly, our method arrives at geometric areas that are very close to each other in each band, as we would expect if the method always extract the same object. Thus even though our dust excess SED is quite different at long wavelengths from previous derivation using similar data, we believe the method used to extract it from the data is more robust, and therefore the derivation of dust mass less dependent on measurement choices.

We also, remark that the dust mass range obtained by Temim & Dwek (2013), between 0.02 and 0.13 M*⊙*, using very different modeling assumption, also brackets our result. We must however point that the photometric uncertainties that we have established in this paper implies that determining which dust model best fits the observation remains essentially unconstrained. More than sensitivity, improving in this situation will require a better understanding of how to separate the object from its background, which can likely be done by a combination of (a) an better statistical study of the background properties, and (b) higher spatial resolution mapping in the longest wavelengths.

Appendix A: Results of calculation

In this appendix, we present a set of visual support material for the results of the image processing of the Herschel maps. In Fig. 10, we present the results of the PACS images: For each first,second, and third row corresponding to 70, 100, and 160 $\mu$m  respectively, the right column represents background subtracted region of FOV (row to a 1500”) and the left column is the zoomed-in region of the crab nebula with the optimum ROI. In Fig. 11, we present the results of the SPIRE images: For each first,second, and third row corresponding to 250, 350, and 500 $\mu$m  respectively, the right column represents background subtracted region of the full FOV (row to a 1500”)and the left column is the zoomed-in region of the crab nebula with the optimum ROI.