{'fraction_non_alphanumeric': 0.0351975815158713, 'fraction_numerical': 0.026403093775151647, 'mean_word_length': 4.285536418344055, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analysis of the emission of very small dust particles from Spitzer spectro-imagery data using blind signal separation methods

This work was conducted as part of the SPECPDR program, dedicated to the study of very small particles and astrochemistry, in Photo-Dissociation Regions (PDRs). We present the analysis of the mid-IR spectro-imagery observations of Ced 201, NCG 7023 East and North-West and rho-Ophiuchi West filament. Using the data from all four modules of the InfraRed Spectrograph onboard the Spitzer Space Telescope, we produced a spectral cube ranging from 5 to 35 um for each one of the observed PDRs. The resulting cubes were analysed using Blind Signal Separation methods (NMF and FastICA). For Ced 201, rho-Ophiuchi West filament and NGC 7023 East, we find that two signals can be extracted from the original data cubes, which are 5 to 35 um spectra. The main features of the first spectrum are a strong continuum emission at long wavelengths, and a broad 7.8 um band. On the contrary, the second spectrum exhibits the classical Aromatic Infrared Bands (AIBs) and no continuum. The reconstructed spatial distribution maps show that the latter spectrum is mainly present at the cloud surface, close to the star whereas the first one is located slightly deeper inside the PDR. The study of the spectral energy distribution of Ced 201 up to 100 um suggests that, in cool PDRs, the 5-25 um continuum is carried by Very Small Grains (VSGs). The AIB spectra in the observed objects can be interpreted as the contribution of neutral and positively-charged Polycyclic Aromatic Hydrocarbons (PAHs). We extracted the 5 to 25 um emission spectrum of VSGs in cool PDRs, these grains being most likely carbonaceous. We show that the variations of the mid-IR (5-35 um) spectra of PDRs can be explained by the photo-chemical processing of these VSGs and PAHs, VSGs being the progenitors of free PAH.


Introduction
The Aromatic Infrared Bands (AIBs), observed in emission in UV irradiated interstellar matter, have been attributed to Polycyclic Aromatic Hydrocarbons (PAHs) by Léger & Puget (1984) and Allamandola et al. (1985). Later observations with the Infrared Space Observatory (ISO) (Kessler et al. 1996) have largely contributed to the study of these macromolecules and their evolution in the ISM, together with many fundamental studies involving both laboratory experiments and quantum chemistry calculations. Although these studies seem to confirm that the AIB carriers are PAH-like species, a precise identification of the mixture of the interstellar molecules has not been possible yet. In parallel to this quest for identification, Cesarsky et al. (2000) gave strong evidence of the existence of Very Small carbonaceous Grains (VSGs), from the ISO observations of the reflection nebula Ced 201. They proposed that these grains, emitting broad bands and a continuum, are destroyed by UV radiation and/or shock waves freeing the AIB carriers in the process. In addition, they suggested that these grains could be present everywhere in the ISM, but detectable in the mid-IR only under special conditions. Photo-Dissociation Regions (PDRs) constitute the transition regions between UV-rich regions and molecular clouds. In PDRs with limited excitation, the infrared emission is dominated by dust particles. To emit at mid-IR wavelengths, these particles must be heated to high temperatures implying that they are very small (from ∼ 1 to 10 nm) and not in thermal equilibrium with the incident radiation field. PDRs are extremely useful to probe the properties of very small dust particles, because the local UV radiation field varies as a function of the distance from the star. This change in the excitation conditions is expected to imply changes in the dust processing and therefore provoke significant variations in the observed infrared spectrum. The analysis of the observed spectral variations can provide valuable clues on the physical and chemical properties of the emitting very small particles. The study of such regions with a high spectral and spatial resolution was first made possible with ISOCAM onboard the Infrared Space Observatory (ISO) (Cesarsky et al. 1996). Most analyses of the ISOCAM spectral cubes have been achieved with simple mathematical decompositions of the observed spectra in bands (the AIBs) and continua (Abergel et al. 2003). However, more sophisticated methods were used elsewhere to describe the observed spectrum on each pixel as the combination of the emission of several chemical populations. This approach requires the use of methods such as Single Value Decomposition (SVD is applied to the covariance matrix of the data, which is equivalent to Principal Component Analysis). Applying SVD to the observations of different PDRs, Boissel et al. (2001) and then Rapacioli et al. (2005) (RJB hereafter) were able to extract the spectra of different types of very small dust particles: PAH-like species and VSGs. Although pioneering, this work was restricted to two objects of the ISOCAM database, because of the limited efficiency of the method.
In the field of signal processing, quite a large number of algorithms have been developed to extract the original signals from observations which are mixtures of these signals. Such algorithms are referred to as Blind Signal Separation (BSS) methods. The most widely used are Independent Component Analysis (ICA) approaches, with algorithms like FastICA (Hyvarinen 1999). The use of ICA is so widespread that it is sometimes presented as the only class of BSS methods, although some other algorithms like Non-negative Matrix Factorisation (NMF) (Lee & Seung 1999) have proven to be efficient.
This paper presents the results of a study conducted with the newly available spectro-imagery data from the Infrared Spectrograph (IRS) (Houck et al. 2004) onboard Spitzer. Our goal is to use both the spectral and spatial information available in IRS data cubes to reveal the origin of the observed spectral variations. The IRS spectrograph provides several improvements when compared to ISOCAM. These are a higher sensitivity (about 0.06 mJy from 6 to 15 µm, and 0.4 mJy from 15 to 38 µm), a better spectral resolution ( λ ∆ λ ∼ 60 − 127 in the low resolution mode) and a range of wavelengths extended to 35 µm. We use two classes of BSS algorithms, i.e. ICA and NMF, to extract the spectra of the different populations of very small particles. This work, carried out on four PDRs, extends previous studies on the following aspects: i) we have gained in sensibility and spectral resolution; ii) a wider range of wavelengths is available; iii) the BSS methods are more powerful than SVD-based methods; iv) as a consequence, our results involve fainter objects.
After describing the observations and data reduction in Sect. 2, we define the BSS methods we have used for our analysis in Sect. 3. We then present the results we have obtained applying these methods to the data (Sects. 4-5) and provide evidence of their efficiency (Sect. 6). We show that the extracted spectra can be attributed to different populations of dust (Sect. 7). Using the extracted in the studied PDRs in Sect. 8.

Selected PDRs
This study was conducted as part of the SPECPDR 1 program  in which a total of 11 PDRs were observed with all the instruments onboard Spitzer. The PDRs selected for this paper, Ced 201, the ρ-Ophiuchi filament, NGC 7023 East and North, have mild radiation fields (200 to 1300 in units of the Habing (1968), G 0 =1.6 × 10 −10 W cm −2 ) and relatively simple geometries. They were chosen among the 11 observations because their IR emissions show the strongest spectral changes across the PDR, thus providing the best information on the nature of the carriers.

Ced 201
Ced 201 is a compact reflection nebula lying in the Cepheus constellation about 420 pc from the Sun (Casey 1991). Unlike common reflection nebulae, which give birth to their own star, Ced 201 is the result of a chance encounter between a B9.5 V star (BD+69 • 1231) and a molecular cloud, at a relative velocity of about 12 km s −1 (Witt et al. 1987). Because of this particularity, Ced 201 is a very attractive object for the spectroscopic study of interstellar dust. Cesarsky et al. (2000) revealed the presence of very small carbonaceous grains, with a spectrum exhibiting bands and continuum emission, identified as the general mid-IR VSG spectrum by RJB.

NGC 7023 East and North-West
NGC 7023 is a well studied reflection nebula located at 430 pc from the sun (van den Ancker et al. 1997). Three PDRs are visible around the Herbig Be illuminating star, HD 200775. The brightest one, NGC 7023 North-West (NGC 7023-NW hereafter), lies about 40 " North-West of the star and another one lies around 70 " South. The last one, NGC 7023 East (NGC 7023-E hereafter), is located further East of the star (∼170").
The spectroscopic studies of the North-West and South PDRs were performed by Cesarsky et al. (1996) using ISOCAM, and by Werner et al. (2004) with IRS, leading to the discovery of new infrared features. In the present paper, we analyse the IRS spectro-imagery observations of NGC 7023-E and NGC 7023-NW.

The ρ-Ophiuchi filament
The ρ-Ophiuchi molecular cloud is a nearby star-forming region. The western part of this cloud is bounded by a filament structure (Oph-fil hereafter), forming an edge-on PDR, illuminated by the B2 HD147889 star (Abergel et al. 1999;Habart et al. 2003).

IRS observations and data reduction
The four PDRs were observed using IRS onboard Spitzer, in the low resolution ( λ ∆ λ = 60 − 127) mode. The data were obtained in the spectral mapping mode with a spatial resolution from 3.6 to 10.2 arcseconds), where the spacecraft moves the spectrograph's slit before each integration, in order to cover a given area of the extended PDR. Each full spectral cube was assembled using the CUBISM software (Smith et al. 2006). The original Basic Calibrated Data (BCD) files are from the S12 pipeline of the Spitzer Science Center. Calibration was achieved using the FLUXCON tables provided by the Spitzer Science Center together with the retrieved data. Additional corrections for extended sources were applied (Smith et al. 2004) and the remaining bad pixels were removed by hand. In the low spectral resolution mode, IRS is divided in two modules: the Short wavelengths Low resolution (SL) module, and Long wavelengths Low resolution (LL) module. Each module is subdivided into two orders: SL2 (5 to 8.7 µm) and SL1 (7.4 to 14.5 µm) for SL and LL2 (14 to 21.3 µm) and LL1 (19.5 to 38 µm) for LL. For Ced 201, both SL and LL modules were used, providing the data from 5 to 38 µm. For NGC 7023-NW we assembled the SL data from the PID 28 Spitzer program 2 with the LL2 order data obtained as "bonus" in the SPECPDR program, while the LL1 order of the IRS was observing NGC 7023-E. For NGC 7023-E and Oph-fil, the SPECPDR observations were only performed with the LL module, to complement the available and good-quality ISOCAM spectral cubes. The full spectral cubes, ranging from 5 to 38 µm, were constructed assembling the cubes from each module, and reprojecting the data on a grid with the lowest resolution (i.e. the LL grid for Ced 201 and NGC 7023-NW, and the ISOCAM grid for NGC 7023-E and Oph-fil). The spectral data from 35 to 38 µm is disregarded because of its poorer quality. The resulting set of data consists of four spectral cubes with wavelengths ranging from 5 to 35 µm. Each spectral cube is thus a 3-dimensional matrix C(p x , p y , λ), which defines the variations of the recorded data with respect to the wavelength λ, for each considered position with coordinates (p x , p y ) in the cube.

Zodiacal light
The contribution from zodiacal light to the observations can become significant when looking at relatively faint objects. This is the case of Ced 201. Fortunately, this nebula is situated in a region of the sky which is slightly contaminated by zodiacal light. We however achieved a background correction for this object. This consists in subtracting the spectrum from a nearby region of dark sky (including background emission and zodiacal light) to the observations. The other objects show stronger emission and are thus less subject to relative contamination by zodiacal light (less than 10 % at 20 µm). For these objects no correction was applied.

MIPS-SED and IRAS data
We have combined the Infrared-Red Astronomical Satellite (IRAS) photometry and Multi-Band Imaging Photometer (MIPS, onboard Spitzer) data with the IRS data to build the spectral energy distribution of Ced 201, between 5 and 100 µm. The IRS spectrum of Ced 201 was obtained by averaging the full IRS cube (see Sect. 2.2). The IRAS data was retrieved from the archive. Because the IRAS beam is larger than the whole region observed in the spectral mapping mode with IRS, it is in principle not possible to compare the fluxes observed with the two instruments. Thus we have made the assumption that the emission of Ced 201 is dominated by the region close to the illuminating star and covered by IRS observations. Therefore the flux calibration of the IRAS points at 12, 25, 60 and 100 µm points was achieved by scaling the 12 and 25 µm points to the IRS-LL data. The MIPS-SED cubes were built stacking the BCD files from the S12 pipeline. The MIPS-SED intensity spectrum was done by averaging the data over the common region with IRS.

Problem overview
Blind Signal Separation is commonly used to restore a set of unknown "source" signals from a set of observed signals which are mixtures of these source signals, with unknown mixture parameters (Hyvarinen et al. 2001). It has e.g. been applied in acoustics for unmixing recordings, or in the biomedical field for separating mixed electromagnetic signals produced by the brain (Sajda et al. 2004) . BSS is most often achieved using ICA methods. In particular, an algorithm called FastICA (Hyvarinen 1999) has proven to be efficient for recovering source signals. An alternative class of methods for achieving BSS is NMF, which was introduced in Lee & Seung (1999) and then extended by a few authors.
In the field of Astrophysics, ICA has been sucessfully used for spectra discrimination in infrared spectro-imagery of Mars ices (Forni et al. 2005), elimination of artifacts in astronomical images (Funaro et al. 2003) or extraction of cosmic microwave background signal in Planck simulated data (Maino et al. 2002). To our knowledge, NMF has not yet been applied to astrophysical problems. However, it has been used to separate spectra in other application fields, e.g. for magnetic resonance chemical shift imaging of the human brain (Sajda et al. 2004) or for analysing wheat grain spectra (Gobinet et al. 2005) The simplest version of the BSS problem concerns socalled "linear instantaneous" mixtures. It is modelled as follows: where X is an m × n matrix containing n samples of m observed signals, A is an m × r mixing matrix and S is an r × n matrix containing n samples of r source signals. The observed signal samples are considered to be linear combinations of the source signal samples (with the same sample index). It is assumed that r ≤ m in most investigations, including this paper. The objective of BSS algorithms is then to recover the source matrix S and/or the mixing matrix A from X.
The correspondence between the generic BSS data model (1) and the 3-dimensional spectral cube C(p x , p y , λ) to be analysed in the present paper may be defined as follows. Each observed signal, consisting of all samples available for this signal, corresponds to a row of X in Eq. (1). In this paper, the sample index is associated to the wavelength λ, and each observed signal consists of the overall spectrum recorded for a cube pixel (p x , p y ). Moreover, each observed spectrum is a linear combination of "source spectra", which are respectively associated to each of the (unknown) types of dust particles that contribute to the recorded spectral cube. Therefore, the recorded spectra may here be expressed according to (1), with unknown combination coefficients in A, unknown source spectra in S and an unknown number r of source spectra. Sections 3.2 and 3.3 introduce the two BSS methods we have tested on these data, i.e. FastICA and an NMF algorithm.

FastICA
FastICA is a statistical BSS method intended for stationary, non-Gaussian and mutually statistically independent random signals (Hyvarinen 1999). It is expressed for zero-mean signals hereafter. In practice, such signals are obtained by first subtracting the sample mean of each observed signal to all signal samples in the corresponding row of X in (1). For the sake of simplicity, the notation X refers to these zero-mean signals hereafter.
The next step of FastICA consists in applying SVD to the covariance matrix of the observed data, i.e. in deriving a matrix Z of transformed signals defined as where M is selected so that: i) all signals in Z are mutually uncorrelated, ii) each of these signals has unit power and iii) the number of signals in Z isr. The value ofr is selected as follows.
When applied to the m signals in X, SVD intrinsically yields m output components. Keeping all these components therefore corresponds to selectingr = m. Instead, if r < m, one may choose to only keep ther output components which have the highest powers, withr selected so that r ≤r < m (see details on p. 129 of Hyvarinen et al. (2001) ). This reduces the dimensionality of the processed data and allows one to combine the following two features: i) usingr ≥ r still makes it possible to recover all source signals from Z and ii) selectingr < m decreases in Z the influence of noise components which exist in real recordings X but were not taken into account in the above data model (2). The basic version of FastICA then extracts a first source signal from the matrix Z. The criterion used to this end consists in maximising the non-Gaussianity of an output signal defined as a linear combination of the signals in Z. Therefore, denoting z a column of Z corresponding to a given sample index, the corresponding sample of the output signal reads where the column vector d is constrained to have unit norm. Several versions of the FastICA method have been defined, depending on which parameter is used to measure the non-Gaussianity of y. The most standard parameter is the absolute value of the non-normalised kurtosis, defined for a zero-mean signal y as where E[.] stands for expectation. Various algorithms may then be used for adapting d so as to maximise that absolute kurtosis parameter. Before the FastICA algorithm was introduced, Delfosse & Loubaton (1995) optimised it by using a standard gradient ascent procedure. FastICA is an alternative, fixedpoint, optimisation algorithm described in Hyvarinen (1999). It has been shown to yield much faster and more reliable convergence than gradient procedures. Moreover, it does not require one to select any tunable parameter (such as the adaptation gain of gradient algorithms). Once a first source signal has thus been extracted as an output signal y, FastICA removes its contributions from all observed signals contained by X. This yields a matrix X ′ which only contain r − 1 sources. The same procedure as above is then applied to X ′ in order to extract another source. This "deflation" procedure is repeated until all sources are extracted from the observations.

Non-Negative Matrix Factorisation
Unlike ICA, NMF is based on the assumed non-negativeness of the source signals and mixing coefficients without requiring independence of the source signals. It aims at recovering the r source signals by approximating the supposedly non-negative matrix X with the following factorisation: where W is a m × r non-negative weight matrix and H is a r × n non-negative matrix of approximated "source" signals. The approximation quality in eq. (5) can be optimised by adapting the non-negative matrices W and H so as to minimise the squared and The adjustement of W and H can be achieved using classical gradient descent but can be problematic (see Lee & Seung 2001). Therefore, the authors have proposed a new algorithm showing that the Euclidian distance is non increasing under the iterative update rule and that the divergence is also non increasing under the rule Thus, the following iterative algorithm can be derived from this result in order to minimise either the euclidian distance or divergence: 1. fix r 2. initialise randomly matrices W and H 3. update these matrices with the update rule (8) or (9) 4. if convergence is reached, then end. Otherwise go back to Step 3.
When convergence is reached, H provides an approximation of "source" signals.

Suitability of BSS methods for analysis of Spitzer-IRS cubes
In order to apply the above BSS methods to the IRS data cubes, it is necessary to make sure that the "linear instantaneous" mixture condition is fulfilled (see Sect. 3.1). Here we consider that each observed spectrum is a linear combination of "source spectra", which are due to the emission of different populations of dust particles. The main effect that can disturb the linearity of the model is radiative transfer as shown by Nuzillard & Bijaoui (2000), because of the non-linearity of the equations. In our case however, this effect is not expected to be a major concern. The mid-IR emission of dust particles occurs in the external layers of the cloud i.e. at optical depth lower or comparable with H 2 emission (Joblin et al. 2007). Using the Meudon PDR code (Le Bourlot et al. 1993, Le Petit et al. 2006 we found that H 2 emission occurs at extinctions of A V ≤ 0.2 in NGC 7023-E and Oph-fil. For Ced 201 which is a particular case, as the star has penetrated the cloud, we use the visual extinction A V = 0.21 of the illuminating star BD+69 • 1231 derived by Witt et al. (1987). From this we conclude that very small particles emit in the infrared at a maximum visual extinction of A V ∼ 0.2 in the studied PDRs. We can then roughly estimate τ 9.7 , the optical depth of the 9.7 µm silicate band which is the main absorption in the mid IR range. Typically, the τ 9.7 /A V ratio is between ∼ 10 and ∼ 20 (Mathis 1990), which yieds τ 9.7 < 2 × 10 −2 . This value is very small, indicating that the effect of radiative transfer on the mid-IR spectra is very weak. Thus, the source spectra are not significantly affected by this reabsorption, and we conclude that the "linear instantaneous" mixture model is fulfilled. This is validated a posteriori as we show that the extraction has the same efficiency around 10 µm, where the silicates absorb, as at the other wavelengths (see Fig.10).

Preliminary tests
To test the BSS methods presented in Sect. 3, we generated two synthetic spectra: one representing a PAH-type emission spectrum, consisting of Lorentzian profiles at 6.2 µm, 7.7 µm, 8.6 µm and 11.3 µm, and one representing the emission of VSGs (RJB), consisting of a broad emission band at 7.8 µm, the 11.3 µm band and a second-order polynomial continuum. These spectra were mixed with a random 100 × 2 mixing matrix A, thus providing 100 artificial observed spectra, which are linear combinations of the initial spectra. These spectra were then analysed using FastICA and NMF. With both methods, we managed to recover the original spectra efficiently. In this case we have artificially mixed the sources, and therefore we can compare the extracted spectra with these sources to measure the efficiency of the methods. With both methods, the correlation between the original and unmixed signals is excellent (correlation coefficients higher than 0.995). However a realistic test must include the effect of noise. The noise in the IRS instrument is quite complex, and strongly depends on wavelength because of the different modules that are used. In order to simply quantify the effect of noise on the BSS we applied a white spatially homogeneous noise. Though this is not representative of the real IRS noise wich would require a detailed study to be estimated, this provides good knowledge on the response of the algorithms to noise. With both methods the results are still good for Signal to Noise Ratios (SNRs) over 2 (i.e. 3 dB). As a comparison, the estimated IRS instrumental noise using the SpecPET tool 3 developed by the Spitzer Science Center, is between 20 and 100 (13-20 dB) depending on the wavelength for Ced 201. We have also estimated the noise using the spectrum of Ced 201 presented on Fig. 9 (position 2) and found a SNR of about 30 (14.8 dB) in the LL region which is consistent with the SpecPET estimation. Below the limit SNR equal to 2, the efficiency of FastICA drops dramatically while NMF is still able to significantly recover the original signals (correlation coefficient greater than 0.85).

Application to Spitzer data
We have applied FastICA and NMF to the Spitzer data. The spectra from a given cube of observations C(p x , p y , λ) are placed in the rows of the X matrix defined in Eq. (1) or (5).
To maintain the efficiency of the methods the noisiest spectra were removed. Thus, X contains the m IRS spectra of a considered PDR, over n points in wavelength. For NMF, Euclidian distance or divergence were minimised using Lee & Seung's al-gorithm (Lee & Seung 2001) implemented with Matlab, to find W and H. The best results were found using the divergence criterion after about 1000 iterations (which takes less than one minute with a a 3.2 GHz processor). The value of r (number of "source" spectra) is not imposed by the NMF method. The strategy for extracting the sources was the following: -Apply the algorithm with a minimum number of assumed sourcesr = 2, to a given dataset, providing 2 sources. → If the found solutions are physically coherent and linearly independent we consider that at leastr = 2 sources can be extracted. → Else, we consider that the algorithm is not suited for analysis (never occurred in our case). -try the algorithm on the same dataset but withr = 3 sources. → If the found solutions are physically coherent and linearly independent we consider that at leastr = 3 sources can be extracted. → Else, we consider that only two sources can be extracted, extraction was over atr = 2 and thus r = 2. -same as previous step but withr = 4 sources.
→ If the found solutions are physically coherent and linearly independent we consider that at leastr = 4 sources can be extracted. → Else we consider that only three sources can be extracted, extraction was over atr = 3 and thus r = 3. . . .
Physically uncoherent spectra exhibit sparse peaks (spikes) which cannot be PDR gas lines. We found r = 3 for NGC 7023-NW and r = 2 for the other PDRs, implying that we could respectively extract 3 and 2 spectra from these data cubes. Similarly, with FastICA, all sources were extracted with r = 2 for Ced 201, Oph-fil and NGC 7023 E, and r = 3 for NGC 7023-NW. Thus we were able to extract the "source" signals from the original cubes using both methods. The choice of the method was therefore based on its suitability, considering the information we have on the sources. In the following section, we argue on why NMF is more consistent in our case.

Choice of the BSS method
NMF and FastICA have different constraints. The strength of FastICA is that the algorithm is guaranteed to converge towards a solution which yields separated spectra. However, this property only holds if the source signals are statistically independent, due to the non-gaussianity criterion optimised by FastICA. This criterion is not ideal in our application, where the "source" spectra are likely to be correlated due to chemical similarities between the different emitting populations (RJB). Moreover, the presence of noise can degrade the efficiency of the separation with FastICA (see Sect. 4.2). On the contrary, the non-negativity of the "sources" and mixing coefficients, used in NMF, are consistent with the non-negativity of emission spectra and of their linear combination. Thus NMF is likely to be more appropriate for the analysis of our data cubes. The only identified drawback of NMF is that its convergence point may depend on the random initialisation of the matrices (see To test this effect, we have ran NMF on the data of NGC 7023-E using 200 different initialisation conditions. We found that the outputs only slightly vary with the initialisation conditions. Since NMF is most likely to have the greatest efficiency for our application, we now detail the results obtained with this algorithm when applied to our data cubes.

Results
In this section we present the extracted spectra for Ced 201, Oph-fil, NGC 7023-E and 7023-NW (Figs. 1-4) and their distribution maps. These spectra are normalised to 1. This has to be done because they have a different intensity depending on the considered pixel. The molecular hydrogen lines have been substracted in order to provide the spectra of dust only and the hydrogen/dust interaction will be the subject of a subsequent paper (Joblin et al. 2007).

Extracted spectra
From each of the cubes of Ced 201, Oph-fil and NGC 7023-E, we were able to extract two spectra. These independently found pairs of spectra have the same characteristics from one object to another. Moreover, for each object, the two spectra exhibit different features. The Signal 1 spectra, presented in the upper part of Figs. 1-3 , are characterised by a rising continuum combined with broad emission bands located at 6.2, 7.8, and 11.4 µm (see Table 1). On the contrary, the Signal 2 spectra, presented on the lower parts of Figs. 1-3, are dominated by the AIB emission at 6.2, 7.6, 8.6, 11.3, and 12.7 µm. In NGC 7023-NW we could extract three signals from the original data. These spectra are presented in Fig. 4. Signal 1 is very similar to Signal  1 found in the other PDRs studied above: it exhibits a broad 7.8 µm band and a steep continuum. Signal 2 and Signal 3 spectra are dominated by AIBs but with different relative intensities. Signal 2 and Signal 3 are respectively dominated by the 11.3 and 7.6 µm bands. Fig. 4. BSS extracted spectra in NGC 7023-NW. In this case, the 7.7 µm band of signal 1 was artificially normalised to 0.4 to be comparable with the spectra extracted in the other objects. The dotted line spectra are those extracted from the ISOCAM observations of the over-all NGC 7023 PDR by RJB.

Spatial distributions of the extracted spectra
The observed data in matrix X. may be expressed as the product of two matrices, W and H eq. (5). Therefore, each observed spectrum from X is defined as a linear combination of the extracted spectra in H. In our case, each IRS spectrum at position (p x , p y ) in a given cube can be written as follows: Obs(p x , p y , λ) = n w(p x , p y ) n S n (λ) where Obs(p x , p y , λ) is the observed spectrum from the pixel indexed by (p x , p y ), S n (λ) is the n th extracted spectrum in H and w(p x , p y ) n are the weights relative to spectrum n. In our case n is restricted to a maximum of 2 extracted spectra, except in the case of NGC 7023-NW where n = 3. The NMF algorithm extracts the sources up to arbirary scaling factors, i.e. it provides y p = η p S p , where η p is an unknown scale factor and S p is the p th "source". Each output can be compared to the observations using the correlation between observations and NMF output. We centered the observations (and therefore the outputs) and considered that the extracted spectra are not correlated. Thus the correlation parameter c p we consider here can be written: where E[.] stands for expectation. We can then calculate γ(p x , p y ) = c p /c p ′ the ratio of the correlations of Signal p and Signal p ′ to the observations: where κ is a constant for a given object. γ traces the ratio of the weights w(p x , p y ) at the position (p x , p y ), and thus the emission ratio between sources p and p ′ , up to the κ scale factor, which does not depend on position (p x , p y ). The maps of Figs. 5-7 show the spatial distribution of the γ emission ratio between Signal 1 and Signal 2, and Fig. 8 shows the c p correlation map for each signal.
In Ced 201 the star is located in the central part of the nebula, and the Signal 2 emission is also concentrated in this region (Fig. 5), whereas Signal 1 dominates in the periphery of the nebula (Fig. 5). In Oph-fil, the star is West of the filament. Again, Signal 2 is dominant in regions closer to star (Fig. 6) whereas Signal 1 emits in deeper regions (Fig. 6). In NGC 7023-E, the illuminating star is located West of the filament. The distribution maps show that the maximum Signal 2 emission is on the edge of the filament (Fig. 7). Further East in the cloud, the emission is highly dominated by Signal 1 (Fig. 7). Finally, in the case of NGC 7023-NW Signal 1 is dominant behind the filament, Signal 2 is dominant on the edge of the filament facing the star, and Signal 3 in front of the filament, in the closest region to the star (Fig. 8).

Efficiency of the reconstruction of the observations using the extracted spectra
BSS methods are quite new to astrophysics, and for this reason they should be used carefully. However, several arguments confirm that in our particular case, NMF can be used efficiently. Contrary to most approaches, there is absolutely no a priori information included in the algorithm, and no subjective physical constraints on the source spectra, so that the extracted spectra are purely based on the information contained in the data. Thus, it is quite convincing that we are able to recover a spatial structure without any information on this included in the algorithm. In this paper, we show that for the studied PDRs, the mid-IR emission spectrum can be fitted everywhere using linear combinations of two or three simple spectra extracted from the considered PDR. Figure 9 shows the reconstruction of the observations taken on three points of the Ced 201 PDR. The first position is at the periphery of the PDR, the second one is closer to the star and the last one is very close to the star. Using linear combinations of the extracted signals (see Sect. 5) we
can reproduce the observations accurately, for the three positions. The residuals from this fit are shown in Fig. 10. The ratio of the power of the observed signal to the power the of residuals is about 30 for position 2 (see Fig. 9), which is consistent with the value of the SNR we have estimated (see Sect. 4.2). The residuals are at noise level, proving that the reconstruction from the extracted spectra of observations is efficient. It is remarkable that we are able to show that in some regions of Ced 201, nearly all the mid-IR emission is due to Signal 1. Also note that, even though we chose to use NMF for the reasons Fig. 6. Map of the emission ratios Signal 1/Signal 2 (upper) and Signal 2/Signal 1 (lower) in Oph-fil. The 5-35 µm integrated mid-IR emission of the filament is in contours (contours 1, 2, 3, 4 are respectively 3.5, 3.7, 3.9, 4.2 ×10 −4 W m −2 sr −1 ). The illuminating star lies South-West of the filament. exposed in Sect. 4.4, FastICA which is a completely different method provided similar results, confirming their relevance.

Continuum-dominated spectrum
In all the extractions, Signal 1 exhibits a very clear association of continuum, and wide AIB emission (Tables 1, 2). This particularity was attributed to the emission of carbonaceous VSGs by Cesarsky et al. (2000) in Ced 201, from ISOCAM observations covering a wavelength range from 5 to 17 µm. The authors also proposed that these grains should be found everywhere in the interstellar medium. RJB have shown that using mathematical analysis, it is possible to identify these grains in other PDRs than Ced 201, from ISOCAM observations. They also mention that in the emission of this population, the position of the "7.7" µm band is shifted towards higher wavelengths which is also what we observe (Tables 1, 2). When compared to spectra 2 and 3, the main characteristics of the Signal 1 bands are widths typically 50 % higher and a clear red-shift of the 7.7 and 11.3 µm bands. The 8.6 µm band is weak or absent and it is not possible to extract the 12.7 µm band from the 12-15 µm plateau. The   Fig. 7. Map of the emission ratios Signal 1/Signal 2 (upper) and Signal 2/Signal 1 (lower) in NGC 7023-E. The 5-35 µm integrated mid-IR emission of the filament is in contours (contours 1, 2, 3 are respectively 1, 1.2, 1.4 ×10 −4 W m −2 sr −1 ). The illuminating star lies West of the filament.
15-18 µm region is difficult to study in details because of the presence of a strong H 2 emission at 17 µm.
The IRS-LL module enabled us to build spectral cubes up to 35 µm (note that the residual fringes observed in the 20-30 µm region are instrumental effects and are not due to the algorithm or signatures of the dust). We have shown that a single spectrum, thus a single continuum, can reproduce the continuum emission all over the Ced 201 PDR (see Sect. 6). This indicates that the shape of the continuum does not vary significantly with distance to the star, though its intensity varies. From Table 1. Central wavelengths of the emission bands in the extracted spectra. The νI ν integrated intensity ratios of the bands compared to the 7.7µm band are written in parentheses. The bands were fitted using lorentzian profiles.

VSGs
PAHs (  this observation, we can conclude that these grains are most likely transiently heated. If this was not the case, their temperature, and therefore the continuum shape, should strongly change as a function of the distance from the star, as it happens for classical big grains (BGs). However we cannot completely rule out the presence of a fraction of grains at thermal equilibrium emitting in the 25-35 µm range. On this basis, and because of the strong similarity of our Signal 1 spectrum and the ISO VSG spectra in the 5 to 16 µm region, either observed (Cesarsky et al. 2000) or extracted mathematically (RJB), we propose to attribute Signal 1 in the 5-25 µm range to the emission of VSGs which appear to be carbonaceous, and probably corresponding to the VSGs used in the Désert et al. (1990) model. In the following section we fit the emission of the Ced 201 PDR using the model of Désert et al. (1990) in order to provide some clues on the contribution of the two populations of grains (VSGs and BGs) to the continuum.

Modelling the emission of VSGs versus BGs
The problem of calculating the infrared emission from VSGs in the ISM has been addressed by several authors (Désert et al. 1990;Siebenmorgen et al. 1992;Draine & Li 2001;Zubko et al. 2004). One question we would like to address is whether BGs contribute to the continuum below 35 µm in the studied PDRs. Although this point would deserve a full study by itself, we provide here first elements. To estimate the emission of BGs, we fitted the total IR emission from Ced 201 gathering the IRS, MIPS, MIPS-SED and IRAS data, using the model of Désert et al. (1990). The input parameters we have used for the model are listed in Table 3. The radiation field was calculated using the Meudon PDR code (Le Bourlot et al. 1993;Le Petit et al. 2006). The stellar spectrum is from the Kurucz library with a temperature of 10 500 K (Kurucz 1991) and a radius of 2.26 R ⊙ . The distance to the cloud, d, was adjusted to obtain the value of G 0 =300 at the interface as given by Young Owl et al. (2002). This gives d=0.0145pc. This PDR has the strongest UV field of the three PDRs for which we have the  Fig. 10. Residuals from the fit presented in Fig. 9 (i.e. observation minus reconstruction).
to emit in the mid-IR continuum. The resulting fit (Fig. 11) shows that the emission from VSGs is dominant up to 50 µm. However, because this model is not completely adapted to the grain populations we consider, and because of the difficulties to calibrate the IRAS data, it is not possible to rule out BGs emission over 25 µm. Moreover, as mentioned by Zubko et al. (2004), the optical properties of BGs and VSGs are not well constrained, implying that their size distribution cannot be derived precisely. Still, the model provides qualitative information on the emitting population, and we can conclude from the fit that the emission in the 25 µm IRAS band is dominated by stochastically heated VSGs. * α is the exponent in the power law of the size distribution: n(a) ∝ a −α where a is the size in nanometers and n(a) gives the numer density of grains with a radius between a and a + da 7.3. Pure band spectra: PAHs Signal 2 is dominated by the AIBs and is similar to the PAHlike emission of RJB. Therefore, we propose that Signal 2 is due to free PAH emission. The three extracted PAH spectra of Figs.1-3, are very similar. In this case, the 11-14 µm region  Désert et al. (1990) with the parameters from Table 3. is dominated by the 11.3 and 12.7 µm bands and no plateau such as the one found in the VSG spectra is detected. Here also the 15-18 µm region is difficult to study because of the high contamination by the H 2 line at 17 µm line. However, the spectrum extracted in the case of Ced 201 exhibits a lower 11.3 µm band than in NGC 7023-E and Oph-fil when compared to the intensity of the "7.7" µm band (see Table 2). This variation of the 11.3/7.7 µm ratio was interpreted by Joblin et al. (1996); Sloan et al. (1999);Bregman & Temi (2005) as an effect of variation in the relative abundance of neutral and positively ionised PAHs (respectively PAHs 0 and PAHs + hereafter), and Flagey et al. (2006) have shown that the effect of size on this ratio was weak. Thus, the PAH extracted spectra in Ced 201, Oph-fil and NGC 7023-E are mixtures of neutral and cationic PAHs. In NGC 7023-NW, RJB were able to identify the emission spectra of two populations respectively dominated by PAHs 0 and PAHs + . In this object we were also able to extract two pure band spectra similar to those determined by RJB which we also attribute to PAHs 0 and PAHs + (see Fig.4).This strengthens the validity of both results. Note however that the PAH + spectrum we extract here has a higher 11.3 to 7.7 µm integrated intensity ratio than the one of RJB (see Tables 1, 2). This is because the observations we analyse here do not cover the regions close to the star. Thus our PAH + population is "less" ionised than the one extracted by RJB, which we will consider as reference hereafter. On the other hand, the PAH 0 we have extracted is close to the one of RJB (see Table 2 and Fig. 4). This extracted spectrum is also similar to the one observed in the HII region of the Horsehead nebula by Compiègne et al. (2006), where PAHs are found to be purely neutral because of the extreme abundance of electrons. From this we can conclude the following: -the PAH + spectrum extracted by RJB is dominated by PAH cations. In the following, we consider this spectrum to be representative of cation emission. This was confirmed by Flagey et al. (2006); -the PAH 0 spectrum can be attributed to a population containing only neutral PAHs.
7.4. Ionisation of PAHs in Ced 201, Oph-fil and NGC 7023-E As mentioned in Sect. 7.3, the PAH spectra of Ced 201, Ophfil and NGC 7023-E are mixtures of neutral and cation PAH emission. The reason why we were not able to extract the individual PAH 0 and PAH + spectra in these objects is that the spectral variations are not strong enough in these cases to provide sufficient information to the algorithm to disentangle the two populations. In order to estimate the normalised fraction of positively ionised PAHs n + /(n 0 + n + ) in these PDRs, we used the PAH 0 and PAH + spectra extracted from the data of NGC 7023-NW by RJB. By comparing the νI ν integrated intensity ratios between the 7.7 and 11.3 µm bands of these extracted spectra (Table 2), to the pure band spectra of Ced 201, NGC 7023-E and Oph-fil (Table 1), we could estimate n + /(n 0 + n + ), using NGC 7023-NW as a reference. This ratio is equal to 0.97, 0.53, and 0.40 for Ced 201, Oph-fil and NGC 7023-E respectively (Table 5). This result shows that Ced 201 has a much higher proportion of ionised PAHs, probably because of its low density (see Table 4) implying a low recombination rate with electrons. RJB and Joblin et al. (1996) also show that a large proportion of ionised PAHs is found in NGC 7023-NW and NGC 1333-SVS3 respectively. This is likely due to the fact that these PDRs are close to the illuminating star. NGC 7023-E and Oph-fil, are dominated by neutral PAHs, probably because these PDRs are far from the star and have a higher density with respect to Ced 201 (see Table 4). Flagey et al. (2006) have defined the R 7.7/11.3 parameter as the ratio between the λI ν integrated intensity of the 7.7 and 11.3 µm band. We have calculated this parameter for Ced 201, NGC 7023-E and Oph-fil (Table 5). Flagey et al. (2006) have shown that the R 7.7/11.3 parameter is related to the ionisation parameter G 0 √ T /n e , where G 0 is the UV field in units of the Habing field, T is the gas temperature and n e the electronic density. Using this result, we calculated the ionisation parameter from the PAH spectra of our PDRs with the R 7.7/11.3 ratio from Table 5. We also calculated this parameter independently using the physical parameters from Table 4. For this calculation, we used a n e /n H ratio of 1.4 10 −4 (Snow & Witt 1995). The results of these two calculations are reported in Table 5. The estimates of the ionisation parameter found here are consistent with the ones calculated from the physical parameters. This shows that our method to quantify the ionisation state of PAHs is consistent with the ones of Flagey et al. (2006) and could be applied to any interstellar PAH spectrum.

Chemical evolution of the very small dust particles in the observed PDRs
The distribution maps presented in Sect. 5.2 clearly show that VSGs and PAHs 0 /PAHs + are not situated in the same regions  Using the results from the model described in Flagey et al. (2006) 2 Using the data from Table 4 3 Using the method described in this paper of the studied PDRs. Going towards the star in the PDRs, VSGs followed by PAHs and eventually PAHs + (case of NGC 7023-NW) are successively dominant. This evolution of the dust populations across the PDRs while the UV flux is becoming greater suggests that these three populations are chemically linked. The disappearance of VSG emission, while the PAH emission increases, suggests that there is transformation of VSGs into PAHs under the action of the UV field, as proposed by Cesarsky et al. (2000) for Ced 201 and by RJB for ρOph-SR3 and NGC 7023-NW PDRs. Here, we report new evidence of this transformation in three regions: Ced 201, NGC 7023-E and Oph-fil PDRs. This is consistent with a scenario in which VSGs could be PAH clusters as suggested by RJB. Interestingly, recent laboratory experiments have shown that coronene clusters can be photo-evaporated into free coronene units (Bréchignac et al. 2005). PAH aggregates would likely form in cold dense regions Bernard et al. 1993;Rapacioli et al. 2006), and their photochemistry would start when the UV field is strong enough. Thus, the observed variations of the mid-IR interstellar spectra (Peeters et al. 2002;Werner et al. 2004;Bregman & Temi 2005) have to be explained as the evolution of a mixture of PAHs and VSGs under the effect of UV flux. The processing of VSGs can explain the diminution of continuum emission while approaching the irradiating star reported by Werner et al. (2004). Bregman & Temi (2005) have found that the 7.7 µm emission band shifts towards shorter wavelengths as the UV field increases, which they attributed to the evolution from anionic to cationic species. According to our scenario, this shift can be explained as the chemical evolution of the mixture, from VSGs, with a band between 7.74 and 7.87 µm, into PAHs with a band between 7.63 and 7.70 µm. In that previous work, the authors show that the R 7.7/11.3 ratio is increasing when moving away from the exciting source, tracing the ionisation state.
However, after a certain distance, the ratio reaches its maximum and starts decreasing again. This was explained by the authors as a result of a greater proportion of anionic PAHs, but is more likely consistent with the presence of VSGs which exhibit a low R 7.7/11.3 band ratio (see Tables 1,2).

Conclusion
Using the data from the Infrared Spectrograph onboard Spitzer, combined with powerful Blind Signal Separation methods, we were able to extract the spectra of two types of very small interstellar dust particles: one carrying mainly the AIB features and the other a mid-IR continuum and broad AIBs. These two populations are identified as PAHs and VSGs respectively. Concerning the pure AIB spectra, we could identify in the case of NGC 7023-NW the emission from PAH 0 and PAH + dominated populations. Using the spectra of RJB, we were able to simply estimate the ionisation fraction of PAHs in various PDRs. Cesarsky et al. (2000) proposed that VSGs should be found everywhere in the interstellar medium. However their detection is far from being obvious. Here we have shown that BSS methods enable us to extract their spectrum taking advantage of both the spatial and spectral information available in the IRS spectral cubes. Thanks to the wide range of wavelengths covered, it was possible to confirm that these grains are responsible for the interstellar mid-IR continuum and could dominate the emission up to 50 µm in cool PDRs. The similarities of their spectral features with the AIBs show that they are carbonaceous. This is consistent with the predictions of Désert et al. (1986) that VSGs are mostly graphitic and not silicates. The distribution maps of PAHs and VSGs further support this conclusion. We indeed show that there is a transformation of VSGs into PAHs, probably under the effect of the UV flux. PAHs and VSGs clearly dominate the emission in the 12 and 25 µm IRAS bands respectively. This result is consistent with the work of Fuente et al. (1992) who proposed that the 25 µm emission in reflection nebulae was due to another type of grains than the ones emitting at 12 µm. It is also consistent with the results of Verter et al. (2000), who have shown that the cloud-to-cloud variations of the mid-IR emission of molecular cirrus can only be explained by changes in the relative abundance of PAHs and VSGs. The evolution from VSGs to PAH molecules can explain the commonly observed mid-IR spectral variations. The recently available data from the MIPS photometer in the Spectral Energy Distribution mode (MIPS-SED) is currently under analysis in order to learn more about the properties of VSGs at longer wavelengths. The PAH 0 , PAH + and VSGs spectra extracted here can be used as a basis to probe the composition of the very small particles in various environments, including external galaxies. To a certain extend, this could provide information on the local physical conditions.