© ESO, 2016

1. Introduction

The object V838 Monocerotis (hereafter V838 Mon) is one of the most enigmatic observed in stellar astrophysics in recent decades. It came to attention when it underwent a powerful eruptive outburst in 2002 (Brown et al. 2002), increasing in luminosity by a factor of 100 over a period of three months. Immediately following this event a spectacular light echo was formed from the outburst light reflecting off the surrounding dust (Bond et al. 2003). Post outburst, V838 Mon varied through spectral types F, G, and K, to a very cool M-giant (Evans et al. 2003). Loebman et al. (2015) found the star to show features of mixed spectral type, from M3 dwarf to mid-L supergiant: the atmosphere is obviously not a conformist.

V838 Mon is situated in a molecular cloud and shares its neighbourhood with a sparse cluster of young B stars (e.g. Afşar & Bond 2007). It has a B-type companion, located 28−250 AU from V838 Mon (Munari et al. 2007; Tylenda et al. 2009) which was engulfed by the material expanding out from V838 Mon after its outburst (Bond 2006; Munari et al. 2007; Kolka et al. 2009). Photometry from before the outburst will be a mix of the companion star and the precursor; and using pre-outburst photometry to fix the status of the precursor requires a reliable measure of the distance. Analysis of high-resolution HST polarimetry images of the light-echo led Bond et al. (2003) to place a lower limit of 6 kpc to the distance of V838 Mon. Working on the same HST material, Tylenda (2004) revised the distance to 8 ± 2 kpc. Munari et al. (2005) favoured an even larger distance of closer to 10 kpc. However, the latest value is 6.2 ± 0.6 kpc, measured by Sparks et al. (2008) on the basis of a detailed study of the light echo. At this distance, the companion completely dominates the pre-outburst photometry (Afşar & Bond 2007), and hence the precursor must have been the fainter star. We will adopt a value of 6.2 kpc in this paper.

The type of star V838 Mon was before the outburst, and the nature of the outburst itself, is as yet unresolved. It is made difficult by the sparsity of data pre-outburst, and the uncertain photometry (i.e. the presence of the secondary star). It has been suggested that the star was massive and hot (e.g. Munari et al. 2005) and the outburst was a thermonuclear event; that the outburst was a born-again event in a solar-mass white dwarf in close orbit with a lower-mass companion (Lawlor 2005); that it was pre-/main-sequence star of 5−10 M_⊙ (Tylenda et al. 2005); and that it was a main-sequence star which suffered a merger with planets or a lower-mass star (e.g. Retter et. al 2007; Tylenda & Soker 2003). Tylenda & Soker (2006) discuss the theories to date to explain V838 Mon, and similar outbursts in two other stars, M31 RV (Rich et al. 1989) and V4332 Sgr (Martini et al. 1999), and they favour the stellar merger mode; in this paper we assume this model is correct, although it makes little difference to many of our calcuations.

There have been numerous studies of V838 Mon (the point source and the circumstellar region) in atomic lines and molecules from the optical to the mid-IR and a few in the far-IR. The spectrum has been evolving dynamically and in terms of composition from the first to the most recent observations. The response of the star to the outburst seems to have been to develop an expanding photosphere and a multi-component circumstellar shell (Lynch et al. 2004; Geballe et al. 2007; Tylenda et al. 2009), this material flowing out and some falling back. The velocities are generally from some tens to 200 km s^-1, and temperatures range from about 2000 K in the photosphere, to a few 100 K in the extended circumstellar region some 200−300 AU out.

The dust in the molecular cloud surrounding V838 Mon can be seen as the light echo sweeps through, and HST images show a very chaotic region of loops and whirls and cavities (Bond et al. 2003; Crause et al. 2005). The appearance of the surrounding dust has led to the speculation that the outburst was not an isolated event. van Loon et al. (2004) found that IRAS and MSX observations show evidence for multiple mass-loss events prior to the 2002 outburst: a 7′× 10′ dust shell is visible in the IRAS data, and the MSX data also show an extended object (1.5′ diameter) at the position of V838 Mon. However, estimates of the combined mass of this dust and gas are in the range 10s−100 solar masses (e.g. Banerjee et al. 2006; Kamiński et al. 2011; Tylenda & Kamiński 2012). This would rule out all of this material having an origin in V838 Mon itself: instead, we are looking at the interstellar cloud from which V838 Mon, and other members of its cluster, were formed. CO studies (Kamiński et al. 2011; Kaminski et al. 2007) also show this extended material, confirming it is massive and moreover cool (~20 K).

Spitzer observations from 2004/5 at 24 μm, 70 μm and 160 μm show extended emission distributed over an arcminute from the star, and moreover that this is coincident with the HST light echo (Banerjee et al. 2006). This infra-red “light echo” is believed to arise from reprocessed thermal emission from dust grains (~1 M_⊙) heated by outward propagating UV–visible photons. Unresolved emission at 24 μm and 70 μm (PSFs of 6′′ and 24′′) is also seen in the Spitzer data at the position of the central star; moreover the flux has increased between 2004/5 and 2007 (Wisniewski et al. 2008), which is suggested to indicate the presence of hot dust freshly condensed in the outburst ejecta. Gemini images taken in 2007 in the mid-IR, at 11.2 μm and 18.1 μm, do not show evidence of extended emission over radial distances of up to 15′′ from the central source (Wisniewski et al. 2008).

It is clear that V838 Mon is an intriguing star, which is varying with time and which is interesting to study across most of the spectrum. We decided to take advantage of the opportunity to observe V838 Mon in the further reaches of the IR with the good spectral and spatial resolution that the instruments onboard Herschel offered. We report on our imaging and spectral study of V838 Mon and its near-by environment.

We note that in this paper we will use the word “extended” to refer to the emission we detect arising from the molecular cloud in which V838 Mon is found (i.e. the region which is coincident with the HST light echo, with a distribution of an arcminute or so from the star). “Point source” refers to anything detected by Herschel as a point source: the full width at half maximum (FWHM) of the PACS 70 μm beam is ~9′′ (55 000 AU at a 6.2 kpc distance), and that of the SPIRE 500 μm beam is ~37.̋ These large sizes mean that the point source includes the star, its circumstellar environment, and probably also interstellar medium (ISM).

2. Herschel observations and data processing

Photometric imaging with the Herschel (Pilbratt et al. 2010) instruments PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 2010) was obtained on three epochs spread over 1.5 yr, and additionally PACS, SPIRE, and HIFI (de Graauw et al. 2010) spectroscopy was obtained on the first epoch. Details of the observations are given in Table 1. In this section we describe the pipeline data reduction and the production of the science-quality maps and spectra.

2.1. PACS photometry

For PACS we use the scan map AOT (astronomer observation template) with two orthogonal (45° and 135°) scanning angles, four scan legs with a length of 12.̋4 each and a cross-scan step of 155′′, and medium scan speed. We obtained four maps at each scan angle: 70 and 160 μm simultaneously, and 100 and 160 μm simultaneously. The final maps are then produced by concatenating the frames from both the scan and cross-scan observations, resulting in maps with a uniform coverage of a square area of just over 6′ on a side. The separate PACS observations were taken back-to-back and within a few months of the SPIRE observations.

The data were processed up to Level 1 with the HIPE (Ott 2010) track 12 pipeline for scan maps and with calibration tree 56. The data were then transferred to be processed with Scanamorphos (Roussel 2012, version 23). The blue and green maps were processed with the default options in the “minimap” mode; for the red frames we also used the option “flat” in Scanamorphos when combing the data from all four datasets to produce a clean background. Glitch detection was done in Scanamorphos. The pixel sizes were chosen for the best results in the subsequent data analysis: 1.0′′ at 70 μm, 1.5′′ at 100 μm and 3.0′′ at 160 μm. Finally, we subtracted the mean of the background flux to set all the maps to the same zero level.

There were slight shifts of the world coordinate system (WCS) of the maps between the epochs – the position of V838 Mon (and a faint background point source) differed by a few pixels. This could be a consequence of slight telescope astrometry errors. Since we later make a direct comparison of the morphology of the extended emission between epochs, we moved the WCS of epochs 2 and 3 to that of epoch 1, independently for each band. This was achieved by fitting the position of V838 Mon (using the HIPE task sourceFitting) and resetting the WCS of the maps to the appropriate sky and pixel coordinates.

The calibration uncertainties for PACS photometry are 5% from the models and 2% reproducibility (information provided on the PACS documentation page on the Herschel portal), and a recent report of a comparison of different map-makers shows that those commonly used on PACS data perform to within 10% of each other for extended emission (Paladini et al. 2013). The errors in the measured photometry for the point and extended sources are in fact dominated by the sky noise and the difficulties in disentangling the two sources.

The PACS maps are shown in Fig. 1. For ease of comparison between epochs and wavelengths we converted the units to Jy/beam. (For the conversion we used the beam areas taken from a PACS Photometer PSF report of 4 April 2012 (“bolopsf_20.pdf”) from the PACS documentation page on the HSC portal. We used beam area values of 40, 54, and 145 sq. arcsec for the blue, green, and red maps.)

2.2. SPIRE photometry

For SPIRE we used the standard small map AOT with a repetition factor three. The observations were done in three filters (PSW = 250 μm, PMW = 350 μm, PLW = 500 μm) simultaneously, and produced maps with a uniform coverage over a 5′ diameter circle.

The SPIRE data were reduced in HIPE (track 12, calibration tree spire_cal_12_2) using the pipeline recipe for small maps (with the pipeline script parameter applyExtendedEmissionGain set to True). The Planck offsets were applied via the task zeroPointCorrection. We opted for map pixel sizes slightly smaller than the default values as this worked better for our subsequent analysis: 4.5′′, 6.25′′ and 9.0′′ for PSW, PMW, and PLW respectively. We produced the maps calibrated for extended emission (units of MJy sr^-1, and used in the analysis of the extended emission) and those calibrated for point sources (units of Jy beam^-1, and used to measure the photometry of the point source). We also created maps with Scanamorphos, using the standard reduction path and created maps with pixel sizes of 3′′, 5′′ and 7′′ for PSW, PMW, and PLW respectively and calibrated for point sources. We subtracted the mean of the background flux to set all the maps to the same zero level (subsequent comparison of the maps was then easier). Background source subtraction was not necessary. The calibration uncertainties for SPIRE photometry are 4% systematic and 1.5% random for a point source, and an additional 4% for the extended emission (Bendo et al. 2013): we adopted a value of 6% overall. As with the PACS photometry, the errors are in fact dominated by the sky noise and the difficulties in disentangling the two sources. The SPIRE maps (those produced in HIPE and calibrated for point sources) are shown in Fig. 1.

2.3. PACS spectroscopy

We obtained ten PACS spectral line scans, split between three concatenated AOTs. We used the chop-nod (medium throw) Pointed mode, performing line-scan spectroscopy with line repetition factors between one and three and a chop-nod repetition factor one at the (default¹) high-density spectral sampling setting. The reference wavelengths of the scans are given in Table 1. The wavelength range of the scans is about 2 μm in the red and 1 μm in the blue.

Our PACS observations consist of 20 discrete segments: the ten requested line scans and ten accompanying parallel scans. We used the background normalisation pipeline script for line scans to reduce the data, in HIPE track 12 with calibration tree version 60. This pipeline uses the telescope background to calibrate the data, and is recommended for observations of long duration and of low continuum flux levels, this latter being the case for V838 Mon.

Diffuse emission from the extended source around V838 Mon pervades the entire spectroscopy field-of-view. The PACS integral field unit (IFU) extends by ~47′′×47′′ (covered by a 5 × 5 array of spaxels), while our PACS maps show that the extended emission covers a patchy area of 1′−2′. This means that flux from the extended emission may be present in the spaxels containing the point source: see Appendix A. To measure the level of contamination of this diffuse emission on the stellar spectrum of V838 Mon, we compared the flux levels of the spectra of the outermost spaxels from our on-source pointings (i.e. the spaxels sitting on the extended emission only, with no point-source contribution) to the spectra from the off-source/blank sky pointings (which contain only telescope background emission) to see how much flux from the extended emission could have been picked up in the on-source observations: there is an excess of 0.1−0.4 Jy in the outer spaxels of the on-source pointings, which we believe comes from the extended source. Under the assumption that this flux is the same over the whole field-of-view, we removed this contribution from the spectra in the final 20 cubes produced by the pipeline. For each cube, we averaged the spectra from the faintest ten outer spaxels; fit a low-order polynomial to the average spectrum; and subtracted the fit from the entire cube. The subsequently-extracted point source spectrum of V838 Mon should then be clear of contamination from the extended emission. It is important to do this before applying the next step – the point source calibration – because that step requires the point source to contain only emission from the itself.

We discovered that V838 Mon was off-centred during our observations: it is located almost exactly between the central spaxel and one of its neighbours. This has consequences for the extraction of the spectrum of the star. The standard task provided in the PACS pipeline scripts to extract a correctly-calibrated spectrum of a point source (extractCentralSpectrum) requires it to be located close to the centre of the central spaxel. However, we were able to modify the task to work on V838 Mon, and while the resulting spectrum will not be as well corrected as a properly-centred source, it is superior to applying no correction. For details see Appendix B.

Examples of the PACS spectra of V838 Mon are shown in Fig. 2 with some spectral line identifications. The absolute flux calibration uncertainty (RMS, combined absolute, reproduceability, and pointing scatter) is about 15% peak-to-peak². For a faint source such as V838 Mon, there is an additional uncertainty of the order ±0.7 Jy on the level of the continuum in a point-source corrected spectrum³; we note that this is not also an additional uncertainty on the integrated line fluxes. The higher uncertainty in this continuum level due to our modification of the point source correction task and the background subtraction will take this to about ±1 Jy.

2.4. SPIRE spectroscopy

We obtained a single SPIRE-FTS pointing with sparse image sampling in the high spectral resolution mode and with 23 scan repetitions of the FTS mechanism to build up the signal-to-noise in the interferogram. The SPIRE FTS simultaneously covers the (SSW) short wavelength band (190−313 μm; 31−52 cm^-1; 957−1577 GHz) and (SLW) long wavelength band (303−650 μm; 15−33 cm^-1; 461−989 GHz), which both have an unapodised spectral resolution of 0.048 cm^-1. The resulting data set consists of a spectrum from the central bolometer (i.e. V838 Mon) and spectra from the surrounding concentric circles of bolometers. The centre-to-centre distance of the bolometers is 33′′ in the SSW band and 51′′ in the SLW band.

We reduced the SPIRE spectra using the standard point-source pipeline described by Fulton et al. (2010) using the HIPE calibration scheme 12_1 described by Swinyard et al. (2014). The size of the SPIRE beam (17′′−40′′ depending on frequency − see Makiwa et al. 2013) is such that some background emission will be captured as well as the source itself (see Appendix A). A careful inspection of the spectra in the surrounding bolometers showed no evidence of strong emission lines and we only removed a polynomial fit to the continuum (to subtract all continuum sources present), resulting in the rectified spectrum illustrated in Fig. 3.

This continuum subtraction adds an extra uncertainty to the subsequent spectral fitting, since our source is faint and the continuum noisy. We estimate an uncertainty in the continuum placement of ± 0.1 Jy. The effect of this is partly accounted for in the fitting errors: a sinc profile placed too low below 0 will be fit worse than a sinc correctly set to the zero level. We have, however, also added this uncertainty, converted to an integrated flux (1.2 × 10^-18 W m^-2), to the errors taken into account in our modelling of the integrated line fluxes (Sect. 6.4).

As with PACS, V838 Mon is offset from the centre of the FoV. Comparing the location of the star on the PACS and SPIRE images and the location of the central SPIRE bolometers, we determined an offset of close to 6.̋ This has consequence for the point-source calibration, which does not account for the flux loss due to off-centring. We used a SPIRE script provided in HIPE (Track 14) to calculate the underestimate of the fluxes due to this off-centring: the loss factor is about 1.04 below 1000 GHz and from 1.25 to 1.40 above 1000 GHz. We note that this task and the correction factors were available only after we had reduced and measured our spectra, but all fluxes and calculations presented here do include this offset correction, unless otherwise stated.

The SPIRE FTS calibration is based on the spectrum of Uranus and detailed analysis of the flux uncertainties (Swinyard et al. 2014) which shows that the absolute flux on a source of this brightness is calibrated to within 6%.

2.5. HIFI spectroscopy

We acquired five double sideband pointing observations with HIFI, further referred to as settings A−E. Details are given in Table 1. Settings A (in band 1b), B (in band 2a), and C (in band 5a) were observed in the standard dual beam switch (DBS) mode, while the FastDBS mode was used for the high-frequency settings D (in band 7a) and E (in band 7b) in order to optimise baseline stability. Data were obtained at nominal resolution for both the high-resolution spectrometer (HRS) and the wideband spectrometer (WBS), in both the horizontal and vertical polarisations. The bandwidth of the HRS is narrow, at 230 MHz, that is 115 km s^-1 and 40 km s^-1 in settings A and E, respectively. The antenna’s half power beamwidth is 34.̋8, 31.̋2, 18.̋4, 12.̋0, and 11.̋9, for settings A through E (Roelfsema et al. 2012, see their Eq. (3)). The HIFI data obtained with the WBS were reduced up to level two using the pipeline available within HIPE. Baseline subtraction and standing-wave removal in the WBS spectra were severely hampered by the noise levels present in the data. Additionally, the lack of coinciding features between the data obtained in the vertical and horizontal polarisations led us to conclude that the HIFI spectra yielded no reliable emission line measurements. This is likely a consequence of the lines being too faint and/or too wide. The ortho-H₂O 7_3,4−7_2,5 transition at 166.8 μm (1797.159 GHz), detected in the PACS observations, is covered by setting E, but no emission feature is seen in either polarisation. Using the point source sensitivities for HIFI, listed by Roelfsema et al. (2012), we find that the flux measured by PACS can indeed not be retrieved at the sensitivities reached in our HIFI observations. Moreover, Tylenda et al. (2009) measured an outflow with a terminal velocity of 215 km s^-1 from numerous lines from V838 Mon, giving a total linewidth of ~2580 MHz at 1797.159 GHz. Considering the WBS bandwidth of 2.56 GHz in setting E, it is clear that if this lines were indeed this broad, we would not only struggle to see it but also to reduce the data in a reliable way.

3. V838 Mon: point source photometry

As is clear from Fig. 1, the emission on the maps consists of a bright point source (V838 Mon itself) and fairly bright extended emission surrounding the star; for SPIRE the two emission sources are particularly intimate. To measure the flux density from the star, aperture photometry using various tasks in HIPE was possible for the PACS data, but for SPIRE the relative strength and unevenness of the extended emission made such measurements highly uncertain. We therefore took an alternative approach for SPIRE, which we also tested for PACS: to scale the point spread function (PSF; beam) and subtract it from the maps until only residual emission remained at the position of the star. The photometry of the scaled beam is then the photometry of V838 Mon. As a refinement of this method, we also performed PSF-subtraction on deconvolved SPIRE maps. Although the results were not as useful as hoped, we believe a summary of this work will be useful for others (Appendix C).

Jumping ahead to the final results for the point source (Sect. 6.1): within the 2σ errors there is no change in the point-source photometry between the three epochs, which cover 16 months. The properties of the dust giving rise to this emission have not changed noticeably during this period.

3.1. PACS

For the aperture photometry we adopted aperture sizes narrow enough to avoid including too much extended emission, and adopted the appropriate aperture and colour corrections: see Appendix D for the details. The errors in the photometry arising from the scatter in the sky background on the maps are 0.02 Jy (70 μm), 0.04 Jy (100 μm), and 0.12 Jy (160 μm), and these values were folded into the errors arising from the photometric measurements. Another source of uncertainty comes from the contaminating flux from the extended emission at the position of the star, which while not strong – the ratio of the peak flux from the star to that in the immediately surrounding extended emission is of the order 40-30-5 for the blue-green-red maps – it is non-negligable. To account for this, we also measured the flux values from the extended emission close to the position of the star and consider this to be the sky flux to be subtracted from the aperture photometry. The photometry given in Table 2 is then mean ± range of these two sets of measurements.

We also performed PSF-subtraction on the PACS maps of V838 Mon, using maps of the PSF taken from the PACS calibration web-page: see below and Appendix E for more detail. The main reason for doing this was to obtain maps of the extended emission only, but the method can also be used to obtain photometry for V838 Mon. Details of the photometric corrections (aperture and colour corrections) made on the scaled PSF are given in Appendix D, and the method itself is explained in Appendix E. The photometry resulting from this are also given in Table 2. We find V838 Mon to have a FWHM consistent with that of the PSF source we used, but with a 2d profile that is somewhat different: residuals in the subtraction of the point source from V838 Mon indicate a difference in the core and the wings of the profile (see Appendix E). This is the main contributor to the photometry errors arising from this process.

The PSF-subtracted maps were used in the study of the extended emission. After subtracting out the point source, the region left behind (which contains the above-mentioned residuals) was replaced by a circle of constant flux (value taken from the local background).

3.2. SPIRE

To measure the photometry of V838 Mon from the SPIRE maps we used a PSF-subtraction method. Aperture photometry was not feasible because of the relative brightness of the extended source: indeed at 500 μm the extended source almost completely overwhelms the point source (see e.g. Fig. E.1). We used the maps produced in HIPE and calibrated for point sources, that is with units of Jy/beam, to obtain the point-source photometry, but performed the same procedure on the maps calibrated for extended sources to create point-source-free maps of the extended emission.

The PSF-subtraction method is essentially a scale-and-subtract, and the details of this work, which was carried out in HIPE, can be found in Appendix E. The photometry is then measured from the beam maps which are scaled sufficiently to remove the point source emission from the V838 Mon maps when image subtraction is done. We used a SPIRE “useful” photometry script in HIPE for the photometry measurements, using the source fitting photometry task (sourceExtractorSussextractor), which works directly on the point-source calibrated maps. See Appendix D for details of the photometric corrections. The errors in the photometry with this method are dominated by the range of values that result in acceptable PSF-subtraction residuals, and for the 500 μm maps are especially large as the beam is large compared to the size of the extended emission; here we add that the result from our PSF subtraction from the deconvolved maps (Appendix C) is that the fluxes at this wavelength may be slightly overestimated. The values and errors given in Table 2 are taken from the acceptable range of PSF scaling factor values.

After subtracting out the point source, the region left behind was replaced by a circle of constant flux (value taken from the local background). These cleaner PSF maps were then used to study the extended emission.

4. V838 Mon: extended source photometry

For PACS, the sum of the flux in all the pixels within a certain contour on the maps, including the point source, was measured. Subtracting the flux of the point source then gave the extended source flux density values. The contours were chosen for each map to encompass the flux close to the level of the sky, with a secondary consideration given to having the area similar for each epoch: the contours values taken lie between 2σ and 4σ above the mean background value on the residual maps, where σ is the standard deviation of the sky flux as measured from about 10 apertures.

The SPIRE maps used were those calibrated for extended emission (units MJy/sr). From these measurements we subtracted the point-source photometric measurements to obtain the photometry of the extended emission only. We compared these values to the values obtained from measuring the photometry of the extended source from the PSF-subtracted maps (Sect. 3.2: the extended emission with the point source subtracted): these agree to better than ± 5%. We then applied the extended-source colour corrections and beam-area corrections to the result: see Appendix D for details of the corrections and Table 2 for the measurements. The error for the fluxes of the extended source is dominated by the error in the point-source photometry. Calibration uncertainties are not included in the error values. The PSF-subtracted maps (SPIRE and PACS) are shown in Fig. 4.

The photometric results given in Table 2 show that there is no difference, within generally less than 2σ, in the point source fluxes from epoch to epoch. For the extended source, from looking at the maps of Fig. 4 and at the better detail in the point-source subtracted deconvolved maps (Appendix C), it does seem that the flux in the northern part of the source decreases with time compared to the flux in the southern part of the source, at all wavelengths (excepting perhaps 500 μm). However, these and any other differences between epochs are mostly within 3σ, and the derived dust properties do not change significantly with epoch (cf. Sect. 6.1).

5. V838 Mon: point source spectroscopy

5.1. PACS

The V838 Mon PACS spectra were measured in HIPE, fitting a low-order polynomial and a Gaussian to all features with the width and peak intensity to appear to be emission lines. We then compared the line list to that from the red supergiant VY CMa: this star has a well-studied PACS–SPIRE spectrum, and it has a cool and extended atmosphere (Royer et al. 2010; Matsuura et al. 2014), as V838 Mon is suggested to have. The comparison was used as a guide for a second pass over the spectra to look for additional features and eliminate those most likely to be noise. Most of the lines we identified this way are of water. The PACS spectra have units of Jy and μm; the line fluxes we report are the integrated intensity of the fitted Gaussian in units of W m^-2 (computed as peak × σ × (2π)^0.5 × 3 × 10^-12/λ², where σ is the width of the Gaussian and the last phrase is to convert the units). Measurement errors were propagated according to the standard rules. Table 3 gives the positions, fluxes and uncertainties of all the lines identified in the PACS (and SPIRE) spectra. The lines for which we could not assign a clear identification are given in Table F.1.

The PACS lines are redshifted by ~0.02−0.03 μm (blue to red wavelengths) compared to our identifications, corresponding to 60−80 km s^-1. The PACS pipeline produces wavelengths corrected to LSR: converting to heliocentric produced redshifts of 40−60 km s^-1 which is close to previous redshift reported values 58 km s^-1 (Tylenda et al. 2009) or 71 km s^-1 (Tylenda et al. 2011). Some of our measured shift, however, will be caused by V838 Mon being located slightly off from the centre of the central spaxel (PACS Observer’s Manual⁴). V838 Mon lies between the central spaxel (12) and its neighbour (17) by less than 1/2 spaxel, and this will cause the spectrum obtained from the central spaxel to have a redshift of ~0.1−0.2 μm (blue to red wavelengths) – although this effect is diluted somewhat since our final spectrum is a combination of spaxel 12 (slight redshift) and 17 (slight blueshift).

In Table 2 we have indicated the approximate values of the spectral continuum flux at the wavelengths of the photometric bands. The spectroscopy and photometry match within their respective uncertainties.

5.2. SPIRE

Fourier transform spectrometers produce spectra in which lines have sinc profile functions (Naylor et al. 2014) that are linearly sampled in frequency. Therefore, all our FTS spectral measurements have been made in frequency space and we have retained the natural sinc profile for all spectral measurements. The integrated line strengths were obtained using a built-in HIPE line fitting script that fits a sinc profile to a given line position using least square fitting of the instrumental line profile (Naylor et al. 2014). The uncertainty in the line is taken from the standard deviation of the residuals around the line after subtraction of the fitted profile, combined with the line fitting errors. Table 3 gives the positions, fluxes and uncertainties of all the lines identified in the SPIRE (and PACS) spectra.

Because our analysis focuses on CO, SiO, and H₂O, we fit only these lines in the SPIRE spectra. All of these species have clear lines present, in both the PACS and SPIRE spectral ranges. Our line-list fitting to the SPIRE spectra is shown in Fig. 3. We note that there are features clearly visible in the spectrum which we have not fit. Some of these could be emission lines that were not on our line list, while others will be noise. An advantage of fitting a sinc profile is that it is not just the main peak that needs to fit the data, but the minor peaks in the wings also. The fitting errors are based on the fit to the entire profile, and hence give a more reliable estimate than if fitting the SPIRE data apodised to a Gaussian profile. However, a disadvantage of spectra with a sinc profile is that the multiple-peak nature of this profile leads to some noise spikes being as bright as real lines. It is for this reason that we fit the SPIRE spectra with a line list: we know the exact position of the expected lines and fit only those. All lines from the line list are included in Table 3, and where the fitting errors exceed the integrated flux, the measured flux is clearly to be considered unreliable.

6. Analysis

6.1. Far infrared spectral energy distribution: (dust) continuum emission

The total integrated Herschel flux density for the central (point source) emission and the large-scale extended infrared emission, for each epoch, are plotted in Fig. 5. In addition, the full SPIRE spectrum of the central point source (smoothed) is shown in light grey. The near- to far-infrared spectral energy distribution of the central (point source) is shown in Fig. 6 (Herschel-only and with literature data added). The fitting results are given in Tables 4 and 5.

Cool extended source.

A modified black body was fitted to the total integrated extended far-infrared Herschel emission (including data from all epochs) given in Table 2. The emitted flux density is modelled with the functional form: S_ν ∝ κ_νB_ν, with κ_ν = 0.1 (ν/ 1000 GHz)^β = 0.1 (300 μm /λ)^β cm²/g (cf. Hildbrand 1983; Beckwith et al. 1990), assuming a gas-to-dust ratio of 100. The power-law index β is an indicator of dust grain composition. β ~ 2 for silicate dust and β ~ 1 for amorphous carbon dust (Mennella et al. 1998; Boudet et al. 2005). For interstellar dust β ranges from 1.7−2.6 (Paradis et al. 2010; Gordon et al. 2014). For this cool material the Herschel wavelength coverage provides a good constraint to the dust temperature, T_dust, dust mass, M_dust, and the emissivity index, β. The best-fit result is shown in Fig. 5 and values for the fitted parameters for the extended emission are labelled in the figure.

The dust temperature for the extended source as derived from Herschel data is T_dust = 19.3 ± 0.7 K, which is close to the value of 25 derived by Banerjee et al. (2006) from 2005 Spitzer-MIPS photometry. The fitted dust emissivity index, β, is 2.43 ± 0.14, which fits into the range of that for ISM dust. The required cold dust mass is 0.57 ± 0.05 M_⊙ (at 6.2 kpc).

Banerjee et al. (2006) derived a value of 0.9 M_⊙ (converted to a distance of 6.2 kpc). Their higher mass value is because of their higher measured flux densities for the extended emission taken from Spitzer data at 70 and 160 μm from 2004/5. The difference (~11 Jy at 70 μm compared to our ~3 Jy) is quite striking. However, our derived mass is sufficiently large that the conclusion of Banerjee et al. – that it is ISM material – still holds. We discuss this further in Sect. 6.3.

Although there is some scatter in the flux densities derived for each epoch, for both the extended source and the point source (see below), there is no evidence for significant variation in the total integrated spectral energy distribution, and the fits to the individual epochs always yield consistent result, as shown in Table 4.

Warm point source.

Using only the Herschel far-infrared observations of the central (unresolved) point source emission, we cannot extract values for the degenerate parameters of dust mass and dust temperature. Therefore we have extended the SED of the central point source to the near-infrared with Spitzer and Gemini mid-infrared results for the unresolved central source in 2007 (Wisniewski et al. 2008). We note that we exclude the Spitzer 2004/5 point source values which are several factors fainter for the emission from the central source (it is suggested this is because new dust was created in the point source in 2007; Wisniewski et al. 2008). The Spitzer 70 μm point-source flux value from 2007 is 7.34 Jy (Wisniewski et al. 2008) which is essentially the same as our measurements, giving some justification to combining the data of 2007 with our later epochs (2011/13). Also available are AKARI/IRC flux densities at 8 and 19 μm (Ishihara et al. 2010), and while these wavelengths are expected to arise primarily from the central source it is not clear from the data how much extended source contamination could be present.

We fit these data with a pure black body and with a “dust” emission model, that is using one or two modified black bodies similar to the procedure adopted for the extended emission. Since it was unclear which of the literature data we can safely use – given it is unknown whether there were any intrinsic variations between then and now – we chose to fit all possible combinations, and we report on the best of these here: the Gemini + Herschel + Spitzer data; the Gemini + Spitzer data; and Herschel + Gemini data. The results of these fittings are given in Table 5 and Fig. 6.

The ranges of values presented in Table 6 show that the pure black body fits give a temperature in the range 300–350 K and a radius of 150–200 AU. The uncertainties on fit parameters are those provided by the least-square-fitting routine. The modified black body fit yielded less satisfactory results, so those of the simple blackbody are more secure, but with this modelling it is possible to estimate a mass of the material. The temperature range is a slightly lower, 230–330 K, and the mass estimates lie around 9 × 10^-4 M_⊙. The values for β are all not far from zero: that is close to a black-body. We note that the reliability of our mass estimate depends on the appropriateness of the dust model.

The uncertainties we quote in Table 5 are the statistical ones from the least-square fitting. The systematic uncertainties for this are much more difficult to ascertain and can be easily factor two to few, especially since we have no constraints from, for example, dust features. A more elaborate model is beyond the scope of this work since we lack simultaneous multi-wavelength data at high spatial resolution and have in fact very limited information on the actual geometry and nature of this emission.

Loebman et al. (2015) studied V838 Mon with mid-IR data taken in 2008 and derive a blackbody temperature of the dust of 285 K (since they have the shorter wavelengths their result is more secure than ours).

6.2. Spatial information on dust properties

Dust temperature maps of the extended emission were created from pixel-by-pixel fitting of the measured flux densities with a modified black body (S_ν ∝ B_νν^β, as above). First the maps at 70, 100, and 160 μm were resampled (Swarp; Bertin et al. 2002) to a common grid with identical pixel sizes, and subsequently convolved to match the PSF of either the 250 μm or 350 μm SPIRE photometry map. For this we adopted the convolution kernels and algorithm from Aniano et al. (2011). The sky background was computed in an “empty” region of the map and subtracted. A fit was then only attempted if the observed flux densities are higher than 10σ_sky of the sky background at at least three wavelengths.

We let β run first as a free parameter and used the maps convolved to the 350 μm SPIRE PSF. The resulting temperature, dust mass, and emissivity index (β) maps (epoch 1 only), as well as the corresponding percentage error maps, are shown in Fig. 7. The noise-dominated low-emission edge regions suffer from degeneracy between dust temperature and dust emissivity, and so these regions have a higher uncertainty than the errors indicate. In the next iteration we omitted the 350 μm SPIRE data to improve the spatial fidelity of the dust images. To derive dust properties we use all maps convolved to 250 μm. Again we limit the fitting to pixels with significant flux above the sky noise level at at least three wavelengths. In addition we fix β to the value found in fitting also the lower resolution 350 μm data. This approach improves the spatial resolution of the final dust maps, while retaining as much as possible information on the dust opacity (through β). Fixing β to a single global value restricts the fits and limits the physical interpretation possible. This procedure was repeated for each of the three epochs.

The final dust temperature, dust mass and emissivity maps of V838 Mon for each epoch are shown in Fig. 8. Mild variations are seen between the three epochs, most noticeably in the temperature of epoch 2 (being slightly higher within the lobes) but the pixel value distribution of T_dust, M_dust, and β agree with each other between the three epochs. The mean results from 2d fitting are given together with the ODR fitting results in Table 4. The comparison in results for these two approaches is acceptable and we do not consider any of the differences to be significant.

6.3. Direct comparision with Spitzer maps

The point source and extended source of V838 Mon are very intimate at the FIR wavelengths. Even for the relatively small beam of Herschel it was necessary to use a point-source extraction method to separate the point-source from the extended source at the longer wavelengths; at the shorter wavelengths it was necessary to have the knowledge obtained from point-source extraction to find an appropriately small aperture with which to do aperture photometry.

For the point source: with the 2005 Spitzer data, Banerjee et al. (2006) used a point-source extraction method to obtain 70 μm flux of ~4 Jy and a non-detection at 160 μm; with the 2007 Spitzer data, Wisniewski et al. (2008) used aperture photometry to obtain a 70 μm flux of ~7 Jy; and we measure 70 μm and 160 μm fluxes of ~7 Jy and ~3 Jy.

For the extended source: Banerjee et al. report fluxes for the extended source at 70 μm and 160 μm of ~11 Jy and ~17 Jy, Wisniewski et al. make no specific statement (which implies a non-detection), and we measure 70 μm and 160 μm fluxes of ~3 Jy and ~11 Jy. The comparison of these epochs indicates that the point source flux increased with time and the extended source flux decreased.

In our opinion a comparison of the Spitzer and Herschel 160 μm results should only be taken as indicative: the point source is so much fainter than the extended source and the Spitzer beam so large, that the reliability of the results is not high. But a comparison of the results at 70 μm should be more reliable.

We decided to make a direct comparison of the total source (the point plus extended source) from the two epochs of Spitzer data with our Herschel data. The 2005 Spitzer maps that were used by Banerjee et al. are public. We downloaded these MIPS observations (AOR 10523648 and 10523904), and from the “maic” FITS files we created single images at each wavelength using the Montage⁵ code. We then convolved the Herschel images to the beam size of the Spitzer data (18′′ at 70 μm and 40′′ at 160 μm) in HIPE, and converted to the MJy/sr units of the Spitzer data. A single value was subtracted from all maps to do a basic background subtraction. We then measured the fluxes from the entire source in an aperture of 80′′ from each map (this being the aperture used by Banerjee et al.). We then did the same for the 2007 Spitzer data used by Wisniewski et al. (AORs 21476096, 21475840, and 21475072 from 2007, for 70 μm and 160 μm both). Applying the appropriate colour corrections to the measurements is made difficult by the different temperatures of the point and extended sources (we measure here a total flux), but in any case are only of the order a few percent. At 70 μm we took a straight average of the respective corrections at 300 K and 20 K and at 160 μm we used the cool value.

We obtain: Spitzer (2005): 15.1 Jy at 70 μm, 16.1 at 160 μm; Spitzer (2007): 12.4 Jy at 70 μm, 12.1 at 160 μm; PACS (2011/12): 11.4 Jy at 70 μm, 15.4 Jy at 160 μm. (Banerjee et al. report fluxes for the 2005 data of 14.7  and 17.5 Jy; the slight difference with our measurements can easily be account for by our very rough background subtraction.) We also compared all images directly to each other. The Spitzer images are slightly larger in the blue and slightly smaller in the red, but we conclude that there is no significant difference in the spatial extent of the emission between 2005 and 2011/12.

Previous investigations by the Herschel Science Centre showed that PACS and MIPS total flux and surface brightness agree within 5−20% and no systematics were found⁶. Taking into account also a 10% calibration uncertainty quoted for both instruments, and our rough background subtraction, the differences in the total source flux that we measure (a drop with time at 70 μm and drop then rise at 160 μm) are barely significant.

In conclusion: looking only at the 70 μm results (where the point source and extended source are the least blended) the flux of the point source rose between 2005 and 2007 and then stayed the same, while the flux of the total source stayed the same (or dropped slightly) and its extent on the sky remained the same. Therefore, the brightness of the extended source must have dropped. Without longer wavelength multi-epoch data we cannot unfortunately say more.

6.4. CO and SiO

To interpret the line emission detected from CO, we plot the energy diagram (Fig. 9). The data are taken from Table 3 and to the flux errors there we have added the calibration uncertainties and SPIRE continuum rectification uncertainty (Sects. 2.3 and 2.4). The CO energy diagram shows two discrete components: a cold component with a peak of the line intensity at the upper state energy E_up< 55 K (lower than J = 4−3), and a warm component with a rising trend with increasing E_up. Our working hypothesis is that the cold component is associated with ISM gas – although it is also possible that it arises from a cold shell (à la Lynch et al. 2004) – while the warm component is associated with the circumstellar envelope/environment (i.e. arising out of the outflow triggered by the stellar impact). We stress that the spectra measured by PACS and SPIRE come from what they see as the point source, and at the spatial resolution at our wavelengths (~9′′ in at the shortest wavelength and ~37′′ at the longest), this area is quite large (2.7 to 11 pc at the distance of 6.2 kpc).

We first compare the distribution of the CO fluxes on the energy-level diagram for V838 Mon with those of what can be considered similar stars (i.e. having an extended atmosphere), and which are well-studied at these wavelengths: the red-supergiant, VY CMa and the AGB star, W Hya. Here we are treating VY CMa and W Hya as templates of circumstellar envelopes created by a constant mass loss, and we wish to see how/where their mass loss differs to that of V838 Mon. It is true that the mass loss from these stars has been ongoing for over a thousand years, a timescale much longer than that for V838 Mon, but it is general trends we are looking for. It is also true that the mass-loss rate of VY CMa is highly variable (Decin et al 2006; Matsuura et al. 2014). However, using the values in Decin et al.⁷ we find that only beyond ~1000  R_⋆ which takes a timescale of over 1000 yr to reach, are the mass-loss rate variations likely to have an effect on the measurements. Hence we should avoid only this part of the comparison.

The CO line intensities of these stars were taken from Khouri et al. (2014) and Matsuura et al. (2014), and these intensities have been scaled to the distance of V838 Mon (6.2 kpc); taking 1.14 kpc for VY CMa (Choi et al. 2008) and 78 pc for W Hya (Knappet al. 2003; Justtanont et al. 2005). The current mass-loss rates of VY CMa and W Hya are estimated to be 2 × 10^-4 M_⊙ yr^-1, and 2 × 10^-7 M_⊙ yr^-1. The overall CO line intensities of the energy diagram reflect the mass-loss rates of these two objects. The mass-loss rate of VY CMa is higher than W Hya, so VY CMa has brighter CO line intensities (after adjusting for its greater distance). The energy diagram shows that V838 Mon has fairly strong CO lines compared to these other two stars (scaled to the same distance).

We can also see that V838 Mon has a different CO curve to those of W Hya and VY CMa: it has a two-component curve, with a clear peak at E_up ~ 50 K, and a gradual rise with increasing E_up, while for the comparison stars there is an initial rise followed by a flattening. Khouri et al. (2014) explained that higher J (higher E_up) transitions trace the inner and warm part of the circumstellar envelope, while lower J transitions trace the outer and cooler part. The clear peak for V838 Mon therefore lies part of the curve that traces the outer and cooler part, which would presumably correspond either to the first epoch of mass loss, that is, that occurring during the outburst, or to pre-existing material. The difference in the slopes at the higher E_up trace differences in the inner part of the circumstellar environment. The slow rise in the CO curve for V838 Mon is because that the warm gas component is a single temperature zone, as our later fitting results will show: see also explanation Fig. 10. For such a gas, the general driver of the curve on the CO diagram is the temperature, and it is the peak of the overall curve that indicates the gas temperature. The two evolved stars we compare to have a more complex story: for example, the temperature of gas from AGB stars is know to vary thanks to their more complex and longer mass-loss history. Their curve on the CO diagram will therefore not be the same as that for a single-temperature zone.

The curve of V838 Mon shows that the molecular envelope consists of a cooler outer part and a warmer inner part, and the comparison to the evolved stars suggests that its envelope not have been made with a constant mass loss from the star. A good (and perhaps obvious) suggestion is that the material is associated with material erupted in 2002, which has remained close to the star and warm. The figure (specifically the comparison of the distance-scaled fluxes) also shows that the scale of gas eruption was large: one event of an eruption could have ejected material equivalent to thousands years of constant mass loss from our comparison stars.

We next derived physical conditions of the two CO-emitting components, by modelling the energy diagrams using the non-LTE radiative transfer code RADEX (van der Tak et al. 2007). RADEX calculates the level populations of molecules, we used CO–H₂ and SiO–H₂ cross-sections and Einstein A-coefficients from the LAMDA molecular and atomic data base of Schöier et al. (2005), which were adopted from calculations by Yang et al. (2010) and Dayou & Balança (2006). We adopted a distance to the source as 6.2 kpc.

The parameters derived and the fitting errors are summarised in Table 5, and in Fig. 9 we show the model fitted results. We now move our attention to a word about the fitting. The errors (1σ measurements+calibration error) are greatest for the cold component, and so simple χ² fitting was weighted to the warm component and ignored the cold component. To overcome this limitation, we (χ²) fit the cold and warm components independently, and summed them. Strictly-speaking, a full radiative transfer modelling of the cold component behind the warm component is the correct path to follow, but given the gaps in our knowledge of V838 Mon, we consider this an unnecessary step, and we ignore the small fraction of the ISM behind the warm component. The straight sum of the cold and warm components is therefore shown in Fig. 9. Inspection of the fits on the plots resulted in some adjustment so we did not underestimate the cold-warm cross-over point at E_up ~ 150 K. The resulting χ²-fitting errors in the temperature are given in Table 6. The errors for N_CO were calculated for a fixed temperature.

We derive temperatures of 400 K and 18 K for the two components. In this and the following discussion, we assume that these two components originate from different gas clouds: the cold component is ISM which is ubiquitous over the observed area (the point source diameter) and our results show it is optically thin, and the warm component comes from the stellar outflow and is smaller than the Herschel point source diameter.

Cold component.

The cold component could either be cold material associated with the outer parts of the molecular envelope, or it could be pre-existing ISM. Our data cannot distinguish between these two cases. The χ² fitting gave a temperature of 18 K, and a best value for the H₂ density of ~1 × 10⁸ cm^-3 for the cold component, although this value is very poorly constrained and could be a factor of 100 lower. In the following we assume a constant density (which is a good assumption for ISM gas). The line width of the cold component has been fixed to 5 km s^-1, as a typical CO line width of molecular gas towards V838 Mon was found to be 2−7  km s^-1 from low J CO lines (Kaminski et al. 2007). For the cold component, we assume that the SPIRE beams are filled by this cold CO gas, and SPIRE beam sizes, which have a frequency dependence, are taken from Makiwa et al. (2013). The estimated column density of the cold CO component is then 1 × 10¹⁵ cm^-2.

Warm component.

When fitting the warm component we find that CO line intensities are almost independent of the (assumed) H₂ density when it is >1 × 10⁶ cm^-3: anything between 10⁶ and 10⁷ cm^-3 is possible. The parameters resulting from the fitting are given in Table 6. We note here that in the RADEX calculations, the escape probability used is that for a slab/cylinder morphology. Lynch et al. (2004) suggests a model of a hollow sphere, and this model is one we wish to follow as much as possible in our calculations. The typical optical depth of the CO lines for V838 Mon as estimated by the RADEX modelling of the Herschel spectra is 3. With this optical depth, the escape probability for a slab is about 30% different to that for a sphere: this difference is very small, especially compared to our other uncertainties, and so we attempted no changes to the RADEX modelling process.

With the fitting results, can we estimate the mass of the gas? This is a multi-step process, since we need to estimate the emitting area of the gas, which is not an output of the model fitting.

First, can we estimate the FWHM of the CO emission lines? One of the model input parameters is the line width (FWHM; Δv), for which we did not have a direct measurement (but see below). The Δv of the warm component is linked with two key model parameters: the CO emitting area, and the saturation limit of CO lines. In the optically thin case, CO line intensities generally increase with higher column density. Once the line becomes optically thick, the line intensities become saturated, and will not increase as the column density increases. Among commonly-observable molecules, CO lines tend to saturate with relatively small column densities. The saturation limit is also coupled with the line width at a given column density. RADEX defines the optical depth at the line centre (τ₀) as τ₀ ∝ N/ Δv, where N is the column density (van Langevelde & van der Tak 2012). Narrowing Δv results in lowering the line intensities at the saturation limit.

A key point of Δv is its link to the CO emitting area. The size of the gas shell created by the eruptive event in 2002 is limited by the expansion velocity and the time since the eruption. Let us assume that the warm gas has expanded at a constant velocity (v_exp) since the eruption, and, as we will discuss later, v_exp = Δv/ 2. We will set the CO emitting area to be Ω_co = πθ²/ 4.25 × 10¹⁰ sr, where θ = Δv/ 2 × t/d in arcsec, t is the time since the eruption, d is distance to V838 Mon, and 4.25 × 10¹⁰ is the conversion factor between sterad and arcsec². As the RADEX output in erg cm^-2 s^-1 is given as the emitting area for the whole sky (4π), to compare with the measured CO line intensities, which come from a smaller area than the whole sky, the model fluxes were converted as f_CO (W m^-2) = f_CO (erg cm^-2 s^-1) Ω_co/ (4π). Hence, the modelled values of f_CO (W m^-2) are ∝ Δv². This analysis method follows the one used for the explosive event of supernova 1987A (Matsuura et al., in prep.; Kamenetzky et al. 2013). This equation also gives the maximum possible line intensities for a given Δv, because of the CO saturation limit. Using this, we narrowed down the range of acceptable Δv values as derived from our line fluxes. By changing Δv over a 50 km s^-1 grid, we established a limit to the Δv which translates to FWHM> 200 km s^-1: narrower than 200 km s^-1 under-predicts the CO intensities.

A value of 200 km s^-1 is the lower limit for the intrinsic FWHM obtained from modelling the fluxes, but what can we measure from the spectra? The lines in the SPIRE spectra, which has an instrumental resolution of 280−970 km s^-1, are unresolved. One CO line is detected in PACS at 162 μm where the instrumental resolution is ~240 km s^-1. The line is rather noisy, but can be fitted with FWHM values of 240−325 km s^-1 – at this upper limit the intrinsic FWHM (decoupling the measured and instrumental widths) is ~220 km s^-1. So the upper limit line width measured from one CO line is close to the lower limit taken from the modelling of the fluxes. To be more certain whether the emission line widths support the modelling, it would be necessary to obtain higher resolution spectra. For the calculations we do here, we will adopt a value of 200 km s^-1 for the intrinsic FWHM of the gas.

Next, what expansion velocity does this FWHM translate into? The FWHM of a line arising from a hollow sphere (thin shell) is 2v_exp (Robinsonet al. 1982), and from a filled sphere is (although the average velocity could be lower; McCray 1993). Taking a value of 200 km s^-1 for the FWHM gives an expansion velocity of the gas of 100 km s^-1 for the hollow sphere case, and over 10 yr this gives a radius of 210 AU; for the filled sphere case the radius is 297 AU. Expansion velocities of various values are reported in the literature as measured from different atomic and molecular lines, with values varying from 100−400 km s^-1 (Loebman et al. 2015); Tylenda et al. (2009) measure a terminal velocity of 215 km s^-1 in 2005. However, there are no contemporaneous measurements of the CO velocity for us to compare to.

Finally, can we estimate the warm CO gas mass? We will estimate the CO mass as r²N_COπ, which equation assumes a slab geometry. Lynch et al. (2004) suggests a model of a hollow sphere, however as long as the (averaged) density along the line of sight remain the same, it does not matter if the geometry is a slab or a hollow sphere. With a column density of 1 × 10²⁰ and a radius of 210 AU, the CO gas mass is estimated to be 7.4 × 10^-5 M_⊙ (see Table 6). We note however that the sources of uncertainty in this value are many – the geometry, the radius, the uncertainty in the column density – and this value should be taken as an estimate only.

SiO and H₂O.

A similar analysis is applied to SiO, with RADEX modelling (Fig. 11) to all the measured transitions (measurement + calibration error included). The collisional cross-sections and Einstein coefficients were taken from the LAMDA molecular and atomic data base (Schöier et al. 2005), which were adopted from calculations by Dayou & Balança (2006). The SiO line intensities have larger uncertainties than CO ones, and the transitions covered are E_up> 380 K, much higher than CO. Figure 11 shows that the kinetic temperatures of SiO could lie in the range 400−1200 K, and that the spectral shape cannot be fit well with a simple component, but needs to involve multiple temperature components. The range of the temperature suggests that SiO has a a similar or higher temperature than the CO warm component – we cannot constrain the options further, and can only state that the SiO-emitting region is associated with CO-emitting regions or lies closer to the heating source, which is presumably the central star.

We did not carry out a similar analysis for the H₂O lines as they even more sensitive to the model used and assumptions adopted.

7. Discussion and conclusion

7.1. Point source

We remind the reader that the point source includes everything inside the PSF of PACS and SPIRE: the PACS 70 μm beam is ~9′′, or 55 000 AU at the distance of V838 Mon, and that of the SPIRE 500 μm beam is ~37.̋

The SPIRE and PACS spectra show CO, SiO and H₂O molecules. The analysis of CO suggests presence of two component temperature, namely warm (~400 K) and cold (~18 K). The inferred SiO kinetic temperature is 400−1200 K. In the past, molecules have been detected from V838 Mon at optical, near- and mid-infrared and millimeter wavelengths (e.g. Lynch et al. 2004; Rushton et al. 2005). This includes the detection of CO rotational-vibration bands at near-infrared wavelengths (Lynch et al. 2004) and rotational transition at millimetre wavelengths (Kaminski et al. 2007), and rotational-vibration bands of SiO (Lynch et al. 2004), as well as SiO masers (Deguchi et al. 2005). We will compare some of these past observations with our molecular line analysis.

Lynch et al. (2004) used observations of V838 Mon’s IR spectra and proposed a model of a swollen star, with a cool and extended atmosphere and a circumstellar dust+gas shell; this model was expanded on using later observations by Lynch et al. (2007) and then Loebman et al. (2015). According to the latter, the radius of the outer shell expanded, and the temperature dropped between 2003 and 2009: from 28 AU to 263 AU and 750 K to 285 K. Chesneau et al. (2014) used interferometric measurements to determine that the star is surrounded by a compact, flattened dust structure, with the dust distributed from about 150 to 400 AU around the star (data from 2011−13). They suggest that this structure is the dusty remains from the merger event that caused the outburst in 2002.

We estimated the radius of warm CO to be 210 AU, constrained by PACS measurements of 162 μm CO line width and CO saturation limit, and assuming the model of an thin spherical shell. Our estimated warm dust radius lies in the range 150 to 200 AU, a value obtained by fitting our data with previously-taken literature data. These values are consistent with the findings quoted previously.

Lynch et al. (2004) reported the detection of CO rotational-vibration bands in the spectra up to 1.3 yr after the outburst. The estimated rotational-vibrational temperature of CO was 650 K. Our SPIRE and PACS spectra, which were obtained about nine years after the eruption, showed the warm CO gas has a kinetic temperature of 400 K. These two measurements show that the temperature has dropped over the last eight years. If one assumes that the eruptive shell is very thin, with a uniform density within, and that this shell has been expanding constantly since 2002, with about ten years between 2002 and our Herschel observations, the corresponding increase in volume and adiabatic cooling would normally lead to the temperature dropping by a factor of 100. That we do not see this suggests that a nearby heating source is present.

V838 Mon has a B3 V companion, located 28−250 AU distant (Munari et al. 2007; Tylenda et al. 2009). Could the companion be heating the gas? The obvious answer is “yes”, given the luminosity of this type of star, however modelling its specific effect on the gas outflow from V838 Mon is beyond the scope of this paper. The gas temperature we measure is what we find, whatever the source of the heating. For the dust it is possible to make a quick and simple calculation of the possible heating by the companion. Using equations from Tielens (2005: 5.43 and 5.44): (1)where we can adopt L_∗ = 1000 L_⊙ and d = 200 AU: the radiation field, G₀ ~ 2.23 × 10⁷. For amorphous carbon dust (β = 1) and taking a representative grain radius, a, of 0.1 μm): (2)or, for silicate dust (β = 2): (3)Naturally the devil lies in the details. The grain radius, for example, depends on the environment in which they form – but increasing the radius to 1 μm only decreases the temperature by 1.6 and 1.15 for the two types of dust. These temperatures are close to the value we measure for the warm dust (~300 K), and so we conclude that at least part of the source of the dust heating could be the companion star.

One puzzle is why dust continued to cool from ~750 K to 280 K in 10 yr as found by Loebman et al. (2015) and this work (~300 K), while the gas temperature is slightly higher (400 K). One possibility is that gas and dust do not have identical temperatures, although this is more commonly found in low density regions with less frequent collisions, and is very unlikely to happen in reasonably dense (>1 × 10⁶ cm^-3) gas, such as V838 Mon appears to have.

In our RADEX analysis we have adopted a spherical geometry. In fact this is probably not really the case. As well as the finding of Chesneau et al. (2014), evidence for a non-spherical morphology comes from Wisniewksi et al. (2003b). They found that spectroscopic and spectropolarimetric variations took place during the outburst, this being indicative of asymmetric geometry; this result is backed up the study of Desidera et al. (2004), who also found polarimetric changes during 2002. Further polarisation variations were observed a few months after the outburst by Wisniewski et al. (2003a), again suggesting the presence of an asymmetrical geometry of scattering material close to the star, and moreover that the distribution of this material had experience significant changes. It is therefore most likely that the circumstellar environment of V838 Mon is not spherical. However, to conduct a more detailed modelling of the gas and dust in the region, its temperature and chemical structure, would require more information about this morphology.

We estimated a dust mass of 9 × 10^-4 M_⊙ for the point source. This is based on a very simple modified blackbody fitting, and we note that the uncertainty in this value is high since it is obtained by fitting to previous epochs of data.

We estimated a CO mass of 7.4 × 10^-5 M_⊙, based on an assumed geometry and shell radius, and with factors of few uncertainty from the underlying CO column density. To turn this into a total mass we need a CO/H₂ ratio: we adopt here a ratio of 4 × 10^-4, which is the O-rich ratio for late-type stars (Willacy & Millar 1997, and taken because of the physical rather than chemical similarity to V838 Mon’s extended environment). This then gives a total gas mass of ~0.2 M_⊙. Taken together, these two values give a gas-to-dust mass ratio of 200. However, this value should only be taken as indicative, as the sources of uncertainty are large. If the outflowing gas from which the CO lines arise is molecular, the conversion from CO to total mass has only a factor of a few uncertainty. However, if the gas is neutral or ionised, then the total gas mass could be much higher.

If the outflowing gas from which the CO lines arise is molecular, our adopted value has only a factor of a few uncertainty. However, if the gas is neutral or ionised, then this value is completely wrong. Hence, a measurement of the gas-to-dust mass ratio requires more and better data. We also point out that it is not clear that the gas and dust we have detected arise from the same component: bear in mind the large size of the point source that Herschel sees (the beam is 55 000 AU at the shortest wavelengths). One should also consider that the dust created in the outflow from V838 Mon – created most likely as a consequence of an infalling low-mass star – may not conform to the usual ISM. Little is known about the type and properties of dust created in these eruptive systems.

7.2. The extended source

The extended source detected in our PACS and SPIRE photometry is shown in the maps in Fig. 4. This region extends over ≈2.7 pc around V838 Mon. The surface-integrated infrared flux (signifying the thermal light echo), and derived dust properties do not vary significantly between the different Herschel epochs (1.5 yr). Our multi-epoch photometry observations show that there has been no clear evolution in the morphology or the fluxes of either the point or the extended source over the 1.5 yr. Our fluxes also do not differ strongly from what was found prior to our Herschel observations. Loebman et al. (2015) conclude that no new dust has formed around the point source since a dust event taking place between 2004/5 and 2006/7 (Wisniewski et al. 2008), as they find no increase in the IR fluxes since then. Our 70μm fluxes are the same as the Spitzer 2007 value (Wisniewski et al. 2008), suggesting no further evolution has occurred.

We compared the fluxes of the emission of the total source (point+extended) between the PACS 2011/12 maps and the Spitzer maps of 2005 and 2007 (as used by Banerjee et al. 2006; and Wisniewski et al. 2008), and also compared the size of the emission from the extended source. There is slight evidence for a decrease in the 70 μm flux from the total source between 2005 and then 2007/2011, but no difference in the extent of the emission when comparing the three epochs of maps at the same flux levels. However, Wisniewski et al. (2008) find that the point source flux increased between 2005 and 2007 (from about 3.8 Jy at 70 μm to 7.3 Jy) and our PACS 2011/12 value is the same as the Spitzer 2007 value. If the flux of the total source has not changed (or slightly dropped), and the extent of the emission from the extended source has not changed, this implies that the brightness of the extended source has dropped.

From the Herschel maps we measured a dust mass of 0.57 M_⊙, which while less than the 0.9 M_⊙ measured by Banerjee et al. (2006) (converted to a distance of 6.1 kpc), is still high enough to result in a large gas + dust mass (if taking the standard ISM value of 100), and hence we agree with the conclusion that this material is ISM rather than arising from the star. We measured a temperature of ~19 K dust for this source, which is only slightly warmer than the T_kin ~ 12 K measured by Kamiński et al. (2011) from CO lines from the ISM around V838 Mon.

As mentioned in the Introduction, van Loon et al. (2004) found evidence for multiple mass-loss events prior to the 2002 outburst: a 7′×10′ dust shell is visible in the IRAS data, and the MSX data also show an extended object (1.5′ diameter) at the position of V838 Mon. The latter exactly matches with the size of the infrared emitting region detected with Herschel. The SPIRE images of V838 Mon reveal patches of faint emission at ~3.̋5 south of V838 Mon, but the maps are not large enough to verify if these structures are part of a larger structure corresponding to the IRAS dust shell.

With our PACS and SPIRE spectra we also detected a cold (~18 K) component of CO gas, which has a temperature close to that measured for the dust. This low temperature suggests that the cold component is associated with ambient ISM gas. The ISM is also responsible for millimetre CO lines which were detected by Kaminski et al. (2007), and are most likely associated with the light echo detected by Bond et al. (2003). Kamiński et al. (2011) found the kinetic temperature of CO to be 12−15 ±5 K, and our kinetic temperature of CO gas is slightly higher but consistent within the uncertainty. There is a large discrepancy in the column density of the cold CO component: in the range 2−6 × 10¹⁶ cm^-2 from Kamiński et al. (2011) while we obtained 1 × 10¹⁵ cm^-2. Some of the discrepancy could be associated with the uncertainty of H₂ density, as a higher density of the collisional partners would require a lower column density to reproduce the same line intensity within that range (this does not apply to warm component). Another possibility is that the cold CO is not from the ISM, but rather from a cold shell that is associated with the warm shell (e.g. the Lynch et al. 2004 model), and hence is part of the Herschel point source.

Our Herschel study of V838 Mon has left us with many questions. What is the relationship between the gas and the dust – are they co-spatial, and how were they each created by the outburst or within the outflow? How many temperature zones are present in the outflow, and does the cold CO gas arise from the ISM or is it directly related to the warm gas? How much dust has been created in the outburst outflow? When will the star return to its pre-outburst state? All very interesting questions for the future.

---

Acknowledgments

We thank John Wisniewski for a thoughtful refereeing process that resulted in an improved paper. N.L.J.C. and K.M.E. acknowledge support from the Belgian Federal Science Policy Office via the PRODEX Programme of ESA. PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC and UKSA (UK); and NASA (USA). This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA.

References

Appendix A: Footprint of the PACS and SPIRE spectrometers

As mentioned in Sects. 2.3 and 2.4, the footprint of the spectrometers extends over the extended emission around V838 Mon, and some of this extended emission “contaminates” the point source. In Fig. A.1 we show the footprint of the PACS and SPIRE spectrometers plotted on two of our epoch 1 maps.

Appendix B: Producing a calibrated point-source PACS spectrum

The HIPE task extractCentralSpectrum works by adding together the spectra of the spaxels that contain the flux of the point source, and then applying a correction for the shape of the beam (which is much wider than a single spaxel): effectively it performs a wavelength-dependent aperture correction to the flux levels in the spectrum. This correction, contained in the PACS calibration tree, is to be used on sources located close to the centre of the central spaxel of the 5 × 5 IFU, in particular because only for this central spaxel has the PACS beam been well calibrated (priv. comm.). The task can extract only the central spectrum or the sum of the central 3 × 3 spaxel box, and the point-source flux density loss correction applied takes this into account. For our data of V838 Mon, the star was located almost exactly between the central spaxel and a neighbour, and hence the task provided could not be used. However, the same procedure could be used, and the resulting spectrum would then be more correct than a simple combination of the brightest three spaxels which contain the flux of V838 Mon (priv. comm. with the PACS instrument team).

Therefore we modified the procedure of extractCentralSpectrum to take into account our non-central position for V838 Mon. For each of our 20 spectral segments we did the following:

• 1.

  Extracted the summed spectrum of the 3 × 3 spaxel box surrounding our brightest spaxels (those with most of the flux from V838 Mon).

• 2.

  Multiplied that by the values in the point-source flux-loss correction tables from the calibration tree (calTree.spectrometer.pointSourceLoss using the central-3 × 3-to-total values).

• 3.

  Then extracted the spectrum only of the brightest two spaxels within the previous 3 × 3 spaxel box: this spectrum has a measurably superior signal-to-noise ratio than that of step 1.

• 4.

  Fit the continuum of the spectra from steps 2 and 3, and divided the first by the second.

• 5.

  Then multiplied that ratio into the spectrum from step 3.

The result is a spectrum with the best possible signal-to-noise ratio but corrected for the point-source flux density losses.

Appendix C: Using deconvolution to assist with the point source photometry

The main difficulty with measuring the point source and extended source fluxes separately from the PACS and SPIRE maps is that the point source is located within the extended emission. For the PACS maps it is possible to still see the two separately, but for the SPIRE maps – with a fainter point source and a larger beam – the two are particularly intimate.

To try to improve the results of the PSF subtraction, we tested our PSF extraction method (Appendix E) on maps which had first been deconvolved: the separation between the point and extended source will then be greater and PSF subtraction should be easier. As a bonus, we would also have deconvolved images of the extended IR bright region to study.

The deconvolution was based on a maximum entropy method (MEM) and follows in general the scheme described by Hollis et al. (1992) for HST images and Ottensamer et al. (2011) for Herschel/PACS. In the deconvolution task, the multiplier is convolved with the PSF model. We tested with the PACS and SPIRE beams (as used also in our PSF-subtraction work) but chose instead to use observations of the point source AFGL 3068, which shows no circumstellar emission and gave very clean results. In the beginning the multiplier is the PSF itself, which is then scaled by the flux ratio of the image and the PSF, and finally re-convolved. At each iteration step the fluxes are compared and the residuals become part of the new multipliers. The iteration is stopped when artefacts become visible, which are mainly manifest as negative fluxes around the central source. The deconvolution was applied to the PACS and SPIRE maps of V838 Mon created by Scanamorphos, and we also deconvolved the maps of AFGL 3068 to then use them for the subsequent PSF-subtraction photometry, using the same technique as outlined in Appendix E.

The photometry resulting from the PSF subtraction method on the deconvolved maps is reported in Table C.1. Unfortunately, this process did not offer as much of an improvement as we had hoped. For the PACS 70 μm and 100 μm maps deconvolution did not make a much difference, and the problem of residuals resulting from a mis-match between the adopted PSF shape and the real PSF for our V838 Mon observations remained (see Appendix E). For the PACS 160 μm maps the deconvolution resulted in a slight improvement: the scatter in the resulting photometry is lower, and the results agree better with the aperture photometry (which are the preferred results for PACS) reported in Table 2. For the SPIRE maps, doing PSF subtraction on the deconvolved maps was a challenge, partly because the deconvolution appeared to alter the FWHM of the PSF differently to that of the V838 Mon maps – an extra smoothing step was necessary for the subsequent PSF subtraction. The photometry obtained from the PSF subtraction on these maps is in the same range as reported in Table 2, however, the 500 μm values are lower than those found previously, and this probably does indicate that our previous results are slightly too high.

In summary, using deconvolution to aid in the PSF subtraction did not help as much as hoped in our work because it was clear that a very good knowledge of the beam for your particular observation (i.e. subjected to the same observing plan, the same data reduction, and the same type of map-making) is crucial to achieving both a good deconvolution and a good PSF subtraction. Given the aims of this paper, we did not feel it was justified to spend more time on the deconvolution. However, in our opinion this method would work if more time could be invested in it.

The images resulting from the deconvolution are presented in Fig. C.1. For these we have subtracted the point source. It is clear that the overall features at each wavelength do not change with epoch, with the possible exception for the PACS maps, where it appears that the northern lobe becomes fainter with time.

Appendix D: Photometry corrections

For aperture photometry of the point source for the PACS and SPIRE data we used the aperture sizes and corrections given in Table D.1. The PACS aperture corrections were taken from photometer.apertureCorrection in calibration tree 65 (the fm7 values), except for the aperture photometry for the scaled PSF maps, for which we used apertures of 57.6′′, and the appropriate aperture corrections provided in the ASCII files that accompanied the tarball of the beam files (see below for details of the beams). The colour corrections for PACS were taken from report on the PACS calibration wiki page (dated April 12, 2011: cc_report_v1.pdf), and the colour and beam-area corrections for SPIRE were taken from the SPIRE calibration tree and a SPIRE photometry script provided in HIPE.

Appendix E: Details on the PSF-subtraction method

For performing PSF subtraction using the scaled PACS/SPIRE beams, we took the beams from:

• PACS: a tarball of FITS files (“PACSPSF_PICC-ME-TN-033_v2.0”) from the PACS documentation web-site from the Herschel portal. We used the examples created from Vesta, with scan speed 20, pixfrac 0.1, from OD 345, and with array-to-map angle +42.

• SPIRE: the normalised beams created in 2012 and provided on the SPIRE documentation web-page from the Herschel portal⁸.

The steps were:

• 1.

  Resample the PSF maps to the pixel sizes of the V838 Mon maps.

• 2.

  Rotate the PSF maps to match the position angle of the V838 Mon maps: the shape of the PSF depends on the scanning speed, angle, and the position angle (this is more true for PACS than for SPIRE). We note that rotation was done on the WCS, rather than on the map as an image.

• 3.

  Measure the sky position of the PSF and of V838 Mon to within a pixel, using sourceExtractor. Change the WCS values in the PSF maps so the PSF star and the point source of V838 Mon are in the same position.

• 4.

  Scale and subtract the fluxes of the PSF maps from the V838 Mon maps (in WCS space, using the HIPE task imageSubtract).

• 5.

  Determine the best scale-and-subtract value via visual inspection of the residual maps and plots of cuts taken at various position angles through the residual maps. The acceptable scaling factors are those that result in residuals that are the flattest.

• 6.

  Finally, replace the fluxes in a small circle around the subtracted star with the value of the flux surrounding that region, to create cleaner PSF-subtracted maps. These were used to create residual-free images for this paper, and to use as a check on the extended-source photometric measurements (Sect. 4).

Figure E.1 shows the SPIRE maps after having scaled-and-subtracted the PSF (Sect. 2.2). Also shown are cuts taken through the original and the residual maps at various position angles, and a cut through the scaled beam. It is from these plots and figures that the decision about the best scaling factors was made.

The PSF-subtraction process did not work for PACS as well as it did for SPIRE. We always had negative residuals, arising from a difference in the shape of the beam from the PSF maps and the effective shape of beam for our V838 Mon maps. It is noted in the report accompanying the PACS beams that the beam-shape depends also on the data reduction and map-making so that a re-reduction of the PSF data, following the same reduction and map-making methods used for the astronomical observations, is recommended. However, as we use the PSF-subtraction method for the PACS data only as a check, we did not do this. An additional consideration is that it was necessary to rotate the beams (the supports of the Herschel secondary lead to a tri-lobal shape for the beams, which orientation depended on the position angle). Rotation inevitably involves an interpolation of the images, and it is possible that this could also account for the mismatch between our point source and the beam profiles. The tri-lobal shape (and hence the resampling of that pattern) is more prominent in the PACS than the SPIRE beams, and this could also explain why the mismatch is greater for the PACS beams.

Appendix F: PACS spectral lines without a clear identification

Lines from the PACS spectra for which we could not find a clear identification, but for which the detection is formally significant, are given in Table F.1. We checked these all against a basic model spectrum based on the parameters of our point source, to look for more CO and H₂O lines, but no matches could be found.

Other potential lines include SiO, CO+, and SO₂, but a proper modelling of these lines is necessary before this can be taken further.

All Tables

All Figures