Dust in the Eye of Andromeda
K. A. Marsh A. P. Whitworth M. W. L. Smith O. Lomax S. A. Eales

Mesoscale model of dislocation motion and crystal plasticity
    Jorge Viñals
    December 30, 2023
============================================================




 
We present new Herschel-derived images of warm dust in the Andromeda 
Galaxy, M31, with unprecedented spatial resolution (∼30 pc), 
column density accuracy, and constraints on the three-dimensional 
distributions of dust temperature and dust opacity index (hence
grain size and composition), based on the new ppmap Bayesian analysis 
procedure. We confirm the overall radial variation of dust opacity index 
reported by other recent studies, including the central decrease 
within ∼3 kpc of the nucleus. We also investigate the detailed
distribution of dust in the nuclear region,
a prominent feature of which is a ∼500 pc bar-like structure
seen previously in Hα. The nature of this feature has been the subject 
of some debate. Our maps show it to be the site of the warmest dust, with a
mean line-of-sight temperature ∼30 K.
A comparison with the stellar distribution, based on 2MASS data, provides
strong evidence that it is a gravitationally induced bar. A comparison
with radial velocity maps
suggests the presence of an inflow towards the nucleus from opposite directions
along this bar, fed presumably by the nuclear spiral
with which it appears to connect. Such behaviour is common in large-scale
bars in spiral galaxies, as is the phenomenon of nested bars
whereby a subkiloparsec nuclear bar
exists within a large-scale primary bar. We suggest that M31 represents
an example of such nesting.




galaxies: individual (M31) – galaxies: spiral – galaxies: bulges – galaxies: kinematics and dynamics – (ISM:) dust, extinction – techniques: high angular resolution 




§ INTRODUCTION


High resolution mapping of dust in spiral galaxies can provide information
regarding star formation rates, the interstellar radiation field,
and the global dynamics of the galaxies themselves. Such studies have 
been greatly aided by data from the Herschel Space Observatory whose
wavelength coverage (70–500 μm) includes the
peak of the spectral range of cool dust emission. In terms of spatial
resolution, the most favourable candidate for study is
the Andromeda galaxy, M31, 
being our nearest spiral neighbor, at an estimated distance of 780 kpc 
<cit.>. Knowledge of the distribution of cool dust provides
constraints on a number of physical phenomena in this galaxy and hence the
Herschel observations represent an important component in
the suite of multiwavelength data which extends from 
radio wavelengths to X-rays. 

The spatial structure of M31, as viewed in the far infrared, is dominated by 
a set of bright concentric rings whose approximate radii, in order of decreasing
prominence, are 10, 5, and 15 kpc <cit.>.  The disk of the galaxy is
inclined at 77^∘ to the plane of the sky, with a tilt axis at
position angle 38^∘ <cit.>.  There is also an inner 
ring, of approximate radius 0.7 kpc, which is less inclined than the main 
structure <cit.> and which may represent the aftermath of a 
collision with the M32 dwarf galaxy <cit.>. 
Studies of the inner isophotes at
optical and near-infrared wavelengths have provided evidence of a large-scale
bar of length ∼10 kpc (<cit.>, although the
existence of such a bar has been questioned by <cit.>.

Evidence for significant radial variations in dust properties has been obtained
by <cit.> using fits of model spectral energy distributions (SEDs)
to Herschel continuum data at six bands in the wavelength range 
70–500 μm. The angular spatial resolution was 36”, corresponding to a
linear resolution of ∼140 pc at the distance of M31.
They found that the dust opacity index, β, increases
inwards, reaching a maximum value at a radial distance ∼ 3 kpc,
inside of which it then decreases towards the nucleus. A similar conclusion
has been reached by <cit.>.

In the present paper we present the results of our investigation of 
the dust distribution in Andromeda using the same
Herschel data as used by <cit.>, but processed
using a new Bayesian analysis procedure, ppmap 
<cit.> which
yields significantly higher spatial resolution (30 pc).
This procedure, outlined in the next section, overcomes a number 
of limitations in standard analysis techniques. These include the requirement 
that all input images be degraded 
to the lowest resolution (normally 36” in the case of Herschel
data) and the assumption that the dust has uniform temperature and composition 
along each line of sight. Such techniques return typically a pair of 
two-dimensional maps, one giving a representative 
dust temperature, T̂__ D(x,y), and the other an approximate 
line-of-sight column density, N̂(x,y).
Here N is the total column density of material (dust plus gas), expressed 
in units of hydrogen molecule masses
per unit area at the source, and the hat symbols (e.g. N̂) 
denote values obtained with the standard analysis procedure.

Expressing the dust distribution in this way, i.e., in terms of column density, 
involves some grain model assumptions. Specifically,
we need to assume a value for the opacity per unit mass at 
some reference wavelength which we take as 300 μm, as discussed in the 
next section. However, ppmap can alternatively
produce maps of optical depth at that wavelength, 
free from any grain model assumptions. The results
could then be compared with optical and/or near-infrared extinction measurements
to extend the extinction curve to 300 μm in a model-independent way.
Such an analysis forms the basis of a companion paper (Whitworth & Marsh,
in preparation).





§ METHODOLOGY


In contrast with previous approaches to the dust mapping problem,
the ppmap procedure drops the assumptions of uniform 
dust temperature and opacity index along the line of sight
<cit.>[The original version treated the opacity index,
β__ D, as constant. We have extended the algorithm by
treating β__ D as an additional state variable. Further details
are given in Appendix A.].
After taking proper account of the point-spread function and spectral response of the telescope, it returns a pair of four-dimensional hypercubes containing, respectively, the expectation values of differential column density,

    N__Tβ   ≡   ∂^2N/∂ T__ D ∂β__ D ,

and their associated uncertainties, as a function of (x,y,T__ D,β__ D).
Here β__ D is the opacity index of the dust in the far-infrared 
i.e. 

    β__ D   =   - dlogκ_λ/dlogλ ;

κ_λ is the mass opacity coefficient of the dust at wavelength 
λ, and the derivative is evaluated in the wavelength interval 
covered by Herschel. Thus N__Tβ
(x,y,T__ D,β__ D) is the best estimate of the line-of-sight
column density at (x,y), in unit 
dust temperature interval about T__ D and unit 
opacity index interval about β__ D. 
The procedure assumes only that the emitting dust is in radiative equilibrium, and that, at the observed wavelengths,  the emission is optically thin. 
It works by defining a grid of discrete values of x, y, T__ D 
and β__ D, which represent sampling locations in the continuous
state space of (x,y,T__ D,β__ D). 
It then iterates towards the set of differential column densities, in the
vicinities of those grid points, that best reproduces the observed 
monochromatic intensities.
The resolution on the plane of the sky (x,y) is increased by a factor 
of about 4.5 relative to the standard analysis procedure, yielding maps of
angular resolution of ∼8”, sampled by 4” square pixels 

The dust temperature sampling interval is arbitrary. For this work, we have used twelve representative dust temperatures that are equally spaced in logarithm between 10.0  K and 50.0  K, i.e. T_1=10.0  K, T_2=11.6  K, T_3=13.4  K, T_4=15.5  K, T_5=18.0  K, T_6=20.8  K, T_7=24.1  K, T_8=27.8  K, T_9=32.2  K, T_10=37.3  K, T_11=43.2  K and T_12=50.0  K. 

The opacity index sampling interval is also arbitrary, and we have used 
four representative values equally spaced linearly 
between β__ D=1.5 and β__ D=3.0, 
i.e. β_1=1.5, β_2=2.0, β_3=2.5 
and β_4=3.0. If we were to invoke more representative 
dust temperatures and/or opacity indices, the computational cost of the 
analysis would increase proportionately but further accuracy would not 
necessarily be gained.

Our column density
scale is based on an assumed opacity of 0.1 cm^2 g^-1 at a wavelength
of 300 μm. The reference opacity is
defined with respect to total mass (dust plus gas). Although observationally
determined it is consistent with a gas to dust ratio of 100 <cit.>,
and we use it here on the assumption that interstellar dust grains in 
M31 have similar properties to those in our own Galaxy.

Given the N__Tβ hypercube, we can integrate out β to obtain two three-dimensional data-cubes containing, respectively, representative values of 

    N__T   ≡   dN/dT__ D  =  ∫_ all β__ D ∂^2N/∂ T__ D ∂β__ D dβ__ D

and their associated uncertainties. N__T(x,y,T__ D) is the best estimate of column density along the line of sight at (x,y), in unit dust temperature interval about ∼T__ D. This data cube can be interpreted in two interesting ways. First, we can take slices at the different representative dust temperatures, T_i, and plot separately the column density of dust in each dust temperature interval, ((T_i-1T_i)^1/2≤ T__ D≤ (T_iT_i+1)^1/2). Second, we can consider individual pixels separately, and analyse the mass-weighted distribution of dust temperature along the line of sight through the pixel. Since the dust temperature depends on the ambient radiant energy density, U__ RAD (in relatively unshielded regions, T__ D∝ U__ RAD^0.20±0.05), this can help to resolve regions along the line of sight that have different ambient radiation fields.

Similarly, we can integrate out T__ D to obtain two three-dimensional data-cubes containing, respectively, representative values of 

    N__β   ≡   dN/dβ__ D  =  ∫_ all T__ D ∂^2N/∂ T__ D ∂β__ D dT__ D

and their associated uncertainties. N__β(x,y,β__ D) is the best estimate of the line-of-sight column density at (x,y), in unit opacity index interval about ∼β__ D. Again, this data cube can be interpreted in two interesting ways. First, we can take slices at the different representative opacity indices, β_j, and plot the column density of dust in each opacity index interval, (β_j-0.25≤β__ D≤β_j+0.25). Second, we can consider individual pixels separately, and analyse the mass-weighted distribution of opacity index along the line of sight through the pixel. Since the opacity index depends on the constitution of the dust (size, composition, shape), this can help to distinguish regions along the line of sight that have different chemical histories. 

Next, we can integrate out both β__ D and T__ D, to obtain two two-dimensional maps of

    N   =   ∫_ all T__ D dN/dT__ D dT__ D

and its associated uncertainty. Here N(x,y) gives the best estimate for the 
line-of-sight column density at (x,y).

We can obtain an estimate of the mass-weighted mean dust temperature, T̅__ D(x,y), along the line of sight at (x,y),

    T̅__ D   =   1/N ∫_ all T__ D T__ D dN/dT__ D dT__ D ,

and higher-order moments of the dust temperature distribution along the line of sight (standard deviation, skewness and kurtosis), as discussed below.

Finally, we can obtain an estimate of the mass-weighted mean opacity index, β̅__ D(x,y), along the line of sight at (x,y)
using an integral analogous to Eq. (<ref>).



§ OBSERVATIONS


M31 was observed using the Herschel PACS and SPIRE instruments 
which provided continuum images in bands centred on
wavelengths of 70, 100, 160, 250, 350, and 500 μm. 
Full details of the observing strategy, calibration, and map-making 
procedures are given by <cit.> and <cit.>. As an example,
Fig.<ref> shows a SPIRE image of M31 in the 250 μm band.

The spatial resolution values of the scan maps, i.e., the beam sizes at full 
width half maximum (FWHM), are approximately
8.5, 12.5, 13.3, 18.2, 24.5, and 36.0 arcsec at the six wavelengths, 
respectively.  For our analysis we use PSFs
based on the measured Herschel beam profiles <cit.>. 



§ RESULTS


The above data have been used to generate a 4D hypercube of differential
column density as a function of sky location, dust temperature, and
dust opacity index, from which we can derive various line-of-sight
integrated quantities such as the
integrated column density of dust plus gas, N, the mean 
dust temperature, T̅__ D, and the mean dust 
opacity index β̅__ D. 
Here, and in the following 
discussion, an overscore (e.g. T̅__ D 
and β̅__ D) denotes mass-weighted averages along the 
line of sight.

Fig. <ref> shows the radial variation of azimuthally averaged values of
these quantities, <N>, <T̅__ D>, 
and <β̅__ D>, at each
galactocentric radius. Galactocentric radii are estimated 
assuming that the midplane of Andromeda makes an angle of θ =77^ o 
with the plane of the sky, and that the position angle of the tilt axis 
is 38^∘.  We note that the radial plot of 
<β̅__ D> in this figure indicates
a central dip in the dust opacity index, consistent with previous reports by
<cit.> and <cit.>. In particular there is detailed
correspondence with Fig. 13 of <cit.>, including the local
maximum at r∼3 kpc.





 §.§ The central region: “ZoomZone" 


In the remainder of the paper, we focus on the small region of Andromeda that 
is marked with a square of width 700 arcsec (2.65 kpc) on Fig. 1, and refer to it as the ZoomZone.
It is chosen because it includes the circumnuclear environment, the site
of some interesting physical processes which are not well understood
<cit.>.
Fig. <ref> shows maps of temperature-differential column density, N_T, 
in this region, at all 12 temperatures. These maps, and the associated
uncertainty maps, are
available in the on-line material as described in Appendix B.



On the first  (top) row of Fig. <ref> are maps of (a) N̂, (b) 
T̂__ D and (c) β̂__ D obtained by 
<cit.> using the standard analysis procedure, cropped at
the boundaries of the ZoomZone. The second (middle) row shows the 
corresponding maps obtained with ppmap, viz. (d) the integrated 
column density, N, (e) the mass-weighted mean dust temperature, 
T̅__ D, and (f) the mass-weighted mean opacity index, 
β̅__ D. These maps illustrate the greatly increased 
resolution obtained with ppmap. In addition, they demonstrate the 
ability of ppmap to dig out large quantities of cold dust that is 
underestimated by the standard analysis procedure. This is illustrated by 
the presence of lower mean dust temperatures (by up to 2 K) in panel (e)
than panel (b), while Fig. <ref> shows the presence of significantly
cooler dust (∼10 K as opposed to ∼15 K) than indicated by the
standard analysis in the corresponding region.  As a consequence, in 
many areas ppmap obtains a lower mean dust temperature 
(T̅__ D<T̂__ D), and a higher net column density 
(N>N̂).  The third (bottom) row shows (g) the relative standard 
deviation, σ_T/T̅_ D, (h) 
the skewness, Skew(T), and (i) the kurtosis, Kurt(T)  of the 
dust temperature distribution along each line of sight. 
These three quantities are defined by:

    σ_T   =   (
    1/N∫_ all T__ DdN/dT__ D (T__ D-T̅__ D)^2 
      dT__ D)^1/2
     Skew(T)   =   (
    1/N∫_ all T__ DdN/dT__ D (T__ D-T̅__ D)^3 
      dT__ D)/ σ_T^3 
     Kurt(T)   =   (
    1/N∫_ all T__ DdN/dT__ D (T__ D-T̅__ D)^4
      dT__ D)/ σ_T^4






They provide information on the form of the temperature 
distribution along each line of sight. In particular, for a pure Gaussian 
distribution we would expect a skewness of
zero and a kurtosis of 3. Examination of Fig. <ref> shows that
these conditions are met over much of the ZoomZone, although there are
localised regions in which the kurtosis is large (>_∼100
in many cases). Large values of the kurtosis indicate the presence of
a significant number of
high-temperature outliers in the tail of the distribution. It is not clear 
what high-kurtosis implies physically;
comparison with published positions of planetary nebulae <cit.>,
X-ray sources <cit.>, and supernova remnants <cit.>
indicate little or no correspondence with high-kurtosis features in the 
ZoomZone, while none of the HII regions listed by <cit.> falls within 
that zone. However, comparison over the full field of M31 shows a 
high degree of correlation with the locations of HII regions from the 
<cit.> survey, in the sense that HII regions are preferentially
located in regions of high temperature kurtosis and, interestingly, low
standard deviation of temperature. This suggests that high kurtosis might be 
a signpost of star formation, and this is confirmed by comparison with
a map of star formation rate <cit.>. It would clearly be interesting
to investigate what kurtosis and the other temperature moments tell us
about star formation, but since the current focus is on the nuclear
region, we defer a discussion of more global phenomena to a future paper.

The uncertainty maps corresponding to the middle row of Fig. <ref>
(column density, temperature, and opacity index) are presented in
Fig. <ref>.

Fig. <ref> shows a set of scatter plots, for the ZoomZone, between all 
pairs of the six quantities, i.e., column density, temperature, opacity 
index, and temperature standard deviation, skewness and kurtosis.
Various correlations are apparent, including an
anticorrelation between kurtosis and standard deviation which, as mentioned
above, appears to be a global characteristic. Also, particularly relevant here 
are the correlations involving the nuclear
bar-like feature, plotted in red. This feature clearly shows some
systematic behaviour, one example being that the dust opacity index,
β̅__ D(x,y),
is positively correlated with dust temperature in the interior of the
nuclear spiral (within 270 pc of the nucleus). This goes against the 
larger-scale trend in which β̅__ D(x,y),
is anticorrelated with T̅__ D within the
central ∼3 kpc. As discussed by <cit.>, various explanations are
possible for that anticorrelation, such as the destruction of grain
ice-mantles in the increasingly intense interstellar radiation field
towards the centre of the galaxy. However, the fact that it switches to
a positive correlation near the very centre would seem to argue against
that explanation, although it might instead indicate 
that a separate population of grains is present.



Fig. <ref> represents a three-temperature composite for the ZoomZone.
It shows that the innermost portion is dominated by a bar-like feature
which appears to consist of a pair of hot, roughly parallel, strands with a 
central crossover. It is distinctly visible at other wavelengths,
such as Hα <cit.> and 70 μm (in our Herschel
image), although less so at 24 μm <cit.>
due to confusion by a dense concentration of bulge stars. 
Its nature has been the subject of debate, but
with a length of only ∼500 pc, it clearly does not correspond to 
the much larger (∼10 kpc) bar inferred from near-infrared morphology 
and dynamical simulations <cit.>.







 §.§ The nature of the bar-like feature


Some information on the 3D geometry of the bar-like feature can be obtained
from a comparison of the column density distribution derived here
on the basis of dust emission, and a
previously published extinction map <cit.>, as shown in panels
(a) and (b) of Fig. <ref>.
These maps show some distinct differences. For example
the filamentary structure which extends between reference points “A" and “B"
in the column density map is absent in the extinction map, suggesting that
it lies behind the bulge and hence behind the source of illumination which
enables foreground features to be detected in extinction. This behaviour was
also noted by <cit.> who identify the feature with a dusty clump
detected in CO and conclude, similarly, that it is located on the far side 
of the bulge. We will thus refer to the feature in question 
as the “rear filament". However, an adjacent filament which extends 
from the nucleus (“+") to point “C" is seen prominently in 
both emission and absorption, indicating that it is closer to us. 
This segment is part of a longer filament which passes
directly over the nucleus in line-of-sight projection. We refer to this
filament as the “central filament".
Based on its elevated dust temperature 
(see Fig. <ref>e), it likely passes directly 
through the centre of the bulge where the interstellar radiation field 
due to evolved giant stars is strongest. 
Its overall position angle is ∼50^∘, although it contains a
central “jog" whereby it runs approximately east-west over the nucleus.

If the central filament is the manifestation of a bar, representing 
orbital motion under the influence of a bar-like gravitational potential
<cit.>,
we would expect to see some evidence of it in the distribution 
of bulge stars. Panel (c) of Fig. <ref> shows that this is indeed
the case—the bar structure is clearly visible in a map of stellar density.
The latter represents the number density of 
2MASS[Two Micron All Sky Survey] sources brighter than
14^ th magnitude at J band, estimated using a circular Gaussian kernel
whose width (FWHM) corresponds to 2 pixels (8 arcsec).  The position
angle of the stellar bar is approximately 45^∘, reasonably consistent
with that of the central filament.
We thus interpret the central dust structure as a bar, and explain its 
apparently double appearance in the column density map as a 
line-of-sight juxtaposition
with the “rear filament" which is actually located further back.

Further elucidation of the nature of the central structure can be obtained
from a comparison with radial velocity information. Since the morphology of the
dust distribution closely resembles that of the ionized gas <cit.>,
we make use of radial velocities of [NII], published previously by 
<cit.>. Such a comparison is facilitated by panel (d) of Fig. 
<ref>, which represents a superposition of the [NII] radial velocities 
on the dust column density distribution.
The general trend of the radial velocities (red-shifted at the upper left
and blue-shifted at the lower right) is consistent with the general
rotational motion in the kinematic model presented by <cit.>,
involving an inner disk tilted at a position angle of 70^∘ and an 
inclination of 43^∘ to the plane of the sky. For the reasons discussed
above, we now believe that this disk is actually a compact spiral containing
a classical bar.

Fig. <ref>(d) also shows that the profile of radial velocity 
along the central filament is
quite distinct from that of the general rotational motion. It changes
abruptly from a red-shift to a blue-shift as the nucleus position is
traversed, at which point the velocity gradient is approximately east-west,
consistent with the local orientation of the jog in the filament.
This velocity structure, together with the apparent positional 
coincidence between the central filament and nucleus, would indicate motion 
either towards or away from the nucleus. However, we can eliminate the latter 
possibility because of the orientation—the lower right portion of the
filament must be pointed away from us since it experiences more optical
extinction then the upper left portion. Such an orientation is, in fact,
consistent with the inferred tilt of the nuclear spiral as a whole,
based on the <cit.> kinematic model.
We thus conclude that the red- and blue-shifts indicate motion towards
the nucleus, i.e., that the filament
constitutes an accretion channel. Since the filament appears to
connect with the nuclear spiral at diametrically opposite points, it
seems likely that these spiral arms provide the source of the inflow.

We can estimate the inflow rate given
the column density of the central filament (∼2×10^20 H_2 masses
per cm^2), the radial velocity ∼100 km s^-1,
the inclination (∼47^∘ to the line of sight) and an estimated width
of 8 arcsec (60 pc). We thereby obtain a total flow rate (from both sides) of
∼0.05 M_⊙ yr^-1. This value greatly
exceeds, by 8 orders of magnitude, the minimum accretion rate
(based on mass-energy equivalence) required to power the central X-ray source, 
the estimated luminosity of which is 
≃1.2×10^36 erg s^-1 <cit.>. It is also exceeds, by
3 orders of magnitude, the Bondi accretion rate estimated by Li et al.
It is, however, close
to the inflow rate required for continuous replenishment of the hot corona
according to the evaporation model of Li et al.

The role of the bar as a conduit for gas flow inwards from the spiral arms 
to the nucleus in barred galaxies has been established observationally
<cit.> and has been the subject of a number of theoretical 
and modelling studies in which the inflow is attributed to
bar-induced torques <cit.>. 
It is not uncommon for
nesting of bars to occur, whereby a subkiloparsec (secondary) nuclear bar
exists within a large-scale primary bar <cit.>. 
We suggest
that M31 represents an example of such nesting. Future dynamical modelling
will be needed to determine whether our inferred nuclear bar can
supply sufficient torque to drive the suggested inflow,
and to establish the dynamical relationship between the primary and 
secondary bars.




§ ACKNOWLEDGEMENTS



We thank Anne-Laure Melchior for interesting discussions on this work.
We gratefully acknowledge the support of a consolidated grant (ST/K00926/1), 
from the UK Science and Technology Funding Council. This work was performed 
using the computational facilities of the Advanced Research Computing at 
Cardiff (ARCCA) Division, Cardiff University. This publication makes use of
data products from the Two Micron All Sky Survey, which is a joint project of
the University of Massachussetts and the Infrared Processing and Analysis
Center/California Institute of Technology, funded by the National Aeronautics
and Space Administration and the National Science Foundation.

99
[Athanassoula & Beaton2006]ath06 
  Athanassoula, E. & Beaton, R. L. 2006, MNRAS, 370, 1499

[Azimlu et al.2011]azi11 
  Azimlu, M., Marciniak, R., & Barmby, P. 2011, ApJ, 142, 139

[Beaton et al.2007]beat07 
  Beaton, R. L., Majewski, S. R., Guhathakurta, P. et al. 2007, ApJ, 658, L91

[Block et al.2006]block06 
  Block, D. L., Bournaud, F., Combes, F. et al. 2006, Nature, 443, 832

[Boulesteix et al.1987]boule87 
  Boulesteix, J., Georgelin, Y. P., Lecoarer, E., Marcelin, M. &
  Monnet, G. 1987, A&A 178, 91

[Ciardullo et al.1988]ciar88 
  Ciardullo, R., Rubin, V. C., Jacoby, G. H. et al. 1988, AJ, 95, 438

[Diaz et al.2017]diaz17 
  Diaz, M. B., Wegg, C., Gerhard, O. 2017, MNRAS, 466, 4279

[Dong et al.2016]dong16 
  Dong, H., Li, Z., Wang, Q. D. et al. 2016, MNRAS, 459, 2262

[Draine et al.2014]draine2014 
  Draine, B. T., Aniano, G., Krause, O. et al. 2014, ApJ, 780, 172

[Ford et al.2014]ford13 
  Ford, G. P., Gear, W. K., Smith, M. W. L. et al. 2013, ApJ, 769, 55

[Fritz et al.2012]fritz12 
  Fritz, J., Gentile, G., Smith, M. W. L. et al. 2012, A&A, 546, 34

[Galliano2018]gal18 Galliano, F.
  2018, MNRAS, 476, 1445

[García-Burillo et al.2009]gar09 
  García-Burillo, S., Fernández-García, S., Combes, F. et al. 2009
  A&A, 496, 85

[Gordon et al.2006]gordon06 
  Gordon, K. D., Bailin, J., Engelbracht, C. W. et al. 2006, ApJ, 638, L87

[Griffin et al.2013]griffin13 
  Griffin, M. J., North, C. E., Amaral-Rogers, A. et al. 2013, MNRAS, 434, 992

[Groves et al.2012]groves13 
  Groves, B., Krause, O., Sandstrom, K. et al. 2012, MNRAS, 426, 892

[Habing et al.1984]habing84 
  Habing, H. J., Miley, G., Young, E. et al. 1984, ApJ, 278, L59

[Hildebrand1983]hil83 
  Hildebrand, R. H. 1983, QJRAS, 24, 267

[Hopkins & Quataert2011]hop11 
  Hopkins, P. F. & Quataert, E. 2011, MNRAS, 415, 1027

[Jacoby, Ford & Ciardullo1985]jac85 
  Jacoby, G. H., Ford, H. & Ciardullo, R. 1985, ApJ, 290, 136

[Juvela et al.2013]juv13 
  Juvela, M., Montillaud, J., Ysard, N., Lunttila, T. 2013, A&A, 556, A63

[Kelly et al.2012]kelly12 
  Kelly, B. C., Shetty, R., Stutz, A. M. et al. 2012, ApJ, 752, 55

[Kirk et al.2015]kirk2015 
  Kirk, J. M., Gear, W. K., Fritz, J. et al. 2015, ApJ, 798, 58

[Knapen, Pérez-Ramírez & 
  Laine2002]knapen02 Knapen, J. H., Pérez-Ramírez &
  Laine, S. 2002, MNRAS, 337, 808

[Laine et al.2002]laine02 
  Laine, S., Shlosman, I., Knapen, J. H., & Peletier, R. F. 2002, ApJ, 567, 97

[Lee & Lee2014]lee14 
  Lee, J. H. & Lee, M. G. 2014, ApJ, 786, 130

[Lewis et al.2015]lewis2015 
  Lewis, A. R., Dolphin, A. E., Dalcanton, J. J. et al. 2015, ApJ, 805, 183

[Li, Wang & Wakker2009]li09 
  Li, Z., Wang, Q. D. & Wakker, B. P. 2009, 397, 148

[Marsh et al.2015]mar15 
  Marsh, K. A., Whitworth, A. P. & Lomax, O., 2015, MNRAS, 454, 4282

[Melchior & Combes2011]melch11 
  Melchior, A.-L. & Combes, F. 2011, A&A, 536, A52

[Melchior & Combes2013]melch13 
  Melchior, A.-L. & Combes, F. 2013, A&A, 549, A27

[Merrett et al.2006]merr06 
  Merrett, H. R., Merrifield, M. R., Douglas, N. G. et al. 2006, MNRAS, 369, 120

[Opitsch2016]op16 Opitsch, M. 2016,
  The bar of the Andromeda galaxy revealed by integral field spectroscopy.
  Ludwig-Maximilians-Universität München, Germany 

[Pilbratt et al.2010]pilb10 Pilbratt,
G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1

[Poglitsch et al.2010]pog10 Poglitsch, A.,
  Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

[Querejeta et al.2015]quer16 
  Querejeta, M., Meidt, S. E., Schinnerer, E. et al. 2016, A&A, 588, 33

[Sellwood & Wilkinson1993]sell93 
  Sellwood, J. A. & Wilkinson, A. 1993, Rep. Prog. Phys., 56, 173

[Shlosman & Heller2002]shlos02 
  Shlosman, I. & Heller, C. H. 2002, ApJ, 565, 921

[Smith et al.2012]smith2012 
  Smith, M. W. L., Eales, S. A., Gomez, H. L. et al. 2012, ApJ, 756, 40

[Stanek & Garnavich1998]stan98 
  Stanek, K. Z. & Garnavich, P. M. 1998, ApJ, 503, L131

[Stiele et al.2011]stiele2011 
  Stiele, H., Pietsch, W., Haberl, F. et al. 2011, A&A, 534, A55

[Vilardell et al.2006]vil06 
  Villardell, F., Ribas, I., & Jordi, C. 2006, A&A, 459, 321

[Wozniak2015]woz15 
  Wozniak, H. 2015, A&A, 575, A7







§ EXTENSION OF PPMAP TO ENABLE Β ESTIMATION


The original version of ppmap <cit.> treated the opacity index,
β__ D, as constant. We have extended the algorithm by
treating β__ D as an additional state variable, thus
enabling the generation of 4D hypercubes of differential column density
as a function of RA, Dec, T__ D, and β__ D.
In doing so it is necessary to guard against the well-known degeneracy
between T__ D and β__ D, whereby the presence of
measurement noise induces an anticorrelation between the estimates of those
quantities. Various approaches have been proposed for dealing with this
problem, most of which are oriented towards fitting the observed spectral 
energy distributions (SEDs) of models of multiple uniform clumps in a
given field. In the Hierarchical Bayesian approach <cit.> the
solution is constrained by the use of a prior distribution designed to
represent the statistics of the source population as a whole and whose
parameters are estimated adaptively from the observations. 
The algorithm then works towards finding solutions which
minimize the dispersion of the estimated parameter values, and appears to
work well provided one can characterise the source population 
by a single mode. However, extending the hierarchical Bayesian scheme to the 
more general imaging problem would present difficulties since it would
belong in the category of single-pixel SED fitting techniques requiring all
input images to be smoothed to the same resolution, thus losing an important
advantage of ppmap. In addition, one would need to assume a single 
temperature and opacity index along each line of sight and this, in itself,
leads to difficulties in deducing dust properties <cit.>.  The
ppmap procedure is not subject to that restriction.

We can, however, achieve the main goal of suppressing
the T__ D-β__ D degeneracy using a
non-hierarchical Bayesian scheme in which we specify
a priori the population statistics regarding β__ D.
For this purpose we assume that the a priori probability distribution 
of possible β__ D values is described by a Gaussian function 
characterised by mean and standard deviation 
β̅ and σ_β, respectively. The estimation
procedure is then essentially the same as described in <cit.>,
except as noted below. As before,
the desired differential column density, represented by 
vector ρ, is obtained 
numerically by a stepwise integration of the differential equation:

    ∂ρ/∂ t + ϕ_1ρ = 0 ,

where ϕ_1 is the conditioning factor which contains information on
both the measurement model and the data, and t is a progress variable which
increases from 0 to 1 during the solution procedure. The only modification
necessary for incorporating the prior information regarding β__ D 
is that the initial condition, ρ_(t=0)=η is now replaced with:

    ρ_n(t=0)= η×αexp[-(β_n-β̅)^2/(2σ_β^2)]

in which the index n refers to the nth cell in the state space
for which the local value of β__ D is equal to β_n, and 
α is a normalisation constant designed to ensure that the average
value of ρ_n(t=0) over all n is equal to η. The latter quantity
represents the a priori dilution, as defined by <cit.>,
for which we have assumed a value of 0.3 in our M31 analysis.

For the Gaussian prior on β__ D 
we have assumed β̅ = 2.0
and σ_β = 0.35. In regions where the signal to noise ratio is high, 
these parameters have essentially no effect on the solution, so that no
bias results <cit.>, and this is true of both
the hierarchical and non-hierarchical Bayesian approaches. However, in regions
of low signal to noise, the algorithm will default to the prior itself so that,
in essence, we are taking β__ D to be 2.0 except where
the data compel us otherwise.




§ SUPPLEMENTARY MATERIAL


We have made available online[http://www.astro.cf.ac.uk/research/PPMAP_M31/] the set of ppmap output files for M31. They consist of a
set of files in FITS format which include:

 
  * 4D hypercube of differential column density, 
	with axes RA, Dec, dust temperature (T__ D), and dust 
	opacity index (β__ D). The sets of grid values 
	of T__ D and β__ D are specified in the FITS 
	header. 
 
  * 4D hypercube of 1σ uncertainty values of
				differential column density
 
  * 2D map of integrated line-of-sight column
	density in units of 10^20 H_2 molecules cm^-2,
				truncated at the 5σ level.
 
  * 2D map of mean line-of-sight dust opacity 
				index, truncated at the 5σ level.
 
  * 2D map of mean line-of-sight dust temperature,
				truncated at the 5σ column density level.
 
  * 2D map of the standard deviation of dust
				temperature along the line of sight, truncated
				at the 5σ column density level.
 
  * 2D map of the skewness of the dust temperature
                                distribution along the line of sight, truncated
				at the 5σ column density level.
 
  * 2D map of the kurtosis of the dust temperature
                                distribution along the line of sight, truncated
				at the 5σ column density level.




