Dust in AGB wind-ISM interaction regions
M. Maercker
T. Khouri 1
M. Mecina 2
E. De Beck 1

Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, 43992 Onsala, Sweden

Department of Astrophysics, University of Vienna, Türkenschanzstr. 17, 1180 Vienna, Austria

In this paper we aim to constrain the dust mass and grain sizes in the interaction regions between the stellar winds and the interstellar medium (ISM) around asymptotic giant branch stars. By describing the dust in these regions, we aim to shed light on the role of low mass evolved stars in the origin of dust in galaxies.

We use images in the far-infrared at  and  to derive dust temperatures and dust masses in the wind-ISM interaction regions around a sample of carbon-rich and oxygen-rich asymptotic giant branch (AGB) stars. The dust temperature and mass are determined in two ways. First directly from the data using the ratio of the measured fluxes and assuming opacities for dust with a constant grain size of 0.1 . We then perform 3-dimensional dust-radiative transfer models spatially constrained by the observations to consistently calculate the temperature and mass. For the radiative transfer models each model contains one constant grain size, which is varied between 0.01  to 5.0 .

We find that the observed dust mass in the wind-ISM interaction regions is consistent with mass accumulated from the stellar winds. For the carbon-rich sources adding the spatial constraints in the radiative transfer models results in preferentially larger grain sizes ($\approx$2 ). For the oxygen-rich sources the spatial constraints result in too high temperatures in the models, making it impossible to fit the observed far-infrared ratio irrespective of the grain size used, indicating a more complex interplay of grain properties and the stellar radiation field.

The results have implications for how likely it is for the grains to survive the transition into the ISM, and the properties of dust particles that later act as seeds for grain growth in the ISM. However, the results for the oxygen-rich sources show that the derivation of dust properties is not straight forward, requiring more complex modelling

§ INTRODUCTION

Dust grains are solid state particles with sizes ranging from a few nm to a few $\mu$m, a large variety of chemical compositions, and with complex surface morphologies. The formation and evolution of dust in astrophysical sources is strongly affected by the density and temperature of the surrounding gas, shocks and gas compression, the radiation fields, and the chemical abundances.
The bulk of the dust observed in galaxies is likely a consequence of grain growth in the interstellar medium (ISM) on seed particles released into the ISM <cit.>. The efficiency of the ISM grain growth is still heavily debated and depends on the chemical evolution of the galaxies and the properties of the seed particles <cit.>. There is little doubt that evolved stars provide the seeds for grain growth in the universe. The seeds are the result of reprocessing of grains that were originally formed in the winds and ejecta of evolved stars. High-mass stars (M$_{\rm{star}}\gtrsim8-10$ ) end their evolution in supernova (SN) explosions. As the ejecta from the SN cools behind the shock moving through the pre-SN wind, dust particles form which are released into the ISM where they can act as seed particles. However, the reverse shock from the SN likely destroys a large fraction of the formed dust <cit.>. In contrast, the majority of all stars that have died in the Milky Way evolved along the asymptotic giant branch (AGB) and ended their lives in a less violent manner. Through their mass-loss, AGB stars are major contributors of dust grains <cit.>. The dust grains are released to the ISM where the stellar wind collides with the surrounding interstellar matter, providing seed particles for ISM grain growth. In models of the Large Magellanic Cloud (LMC), total dust production rates fit the observed dust mass if destruction through a SN reverse shock is included <cit.>. In this case the dust from AGB stars may account for up to 70% of the total dust budget observed at the present time. At earlier times the contribution from SN explosions dominates the dust production, owing to the longer evolutionary timescales for low-mass stars. At which point dust production is dominated by AGB stars is very model dependent, and may happen already for galaxies at redshift 4-5 <cit.>.

It is possible that up to 90% of the dust produced in stars comes from low-mass stars <cit.>, and the detection of pre-solar grains from low-mass stars in meteorites shows that (at least) some of this dust survives and can act as seeds for grain growth in the ISM <cit.>. In fact, models of the origin of dust in the Milky Way can not explain the observed dust masses from stellar sources alone, strongly indicating that ISM grain growth must take place <cit.>.

The dust grains from AGB stars form in the extended atmosphere of the stars, where radiation pressure on the grains initiates the stellar wind <cit.>. The wind and resulting circumstellar envelope (CSE) may be asymmetric, are likely inhomogeneous, and experience continued grain growth and processing in denser regions out to a few stellar radii. Eventually the dust will be released into the ISM, when the stellar wind collides with the surrounding medium in the wind-ISM interaction region <cit.>. The dust grains may contain complex molecules, and shield them from destruction due to the interstellar radiation field and passage into the ISM. The mass, effectiveness of shielding, and the properties of seeds released into the ISM, depend on the grain sizes, the porosity, and the structure of these stellar grains <cit.>. These properties in turn depend on the physical conditions in the dust formation and destruction regions. Attempts to determine the properties of dust grains have been made by looking at distinct features <cit.> and continuum emission <cit.>.

The properties of these original dust grains formed in the winds of evolved stars are directly connected to the properties of the seeds in the ISM, and the efficiency of the ISM-growth depends directly on the properties of the stellar grains that are released into the ISM. Specifically, properties such as the size and structure of the seed grains are critical, as they determine the likelihood for the grains to grow <cit.>. Theoretical models of grain growth in the ISM combine critical dust properties in a generic "sticking coefficient", constituting one of the largest uncertainties in grain growth models <cit.>, emphasising the importance of determining these properties directly.
In this paper we study a sample of wind-ISM interaction regions observed at far-infrared wavelengths as a first step to study the properties of the dust grains released into the ISM by AGB stars and their significance for the origin of dust in the ISM. We focus our analysis on the interaction regions between the wind and the ISM to empirically study the properties of the dust in these regions, instead of the formation of the dust in the inner wind. A similar study using the same data was done by <cit.>, albeit with a different approach and a focus on the overall geometry of the interaction regions rather than the dust properties. Section <ref> describes the archival data we use, Section <ref> we use to determine dust masses and temperatures. The results are presented in Sect. <ref>, including a comparison to the analysis in <cit.>.. Finally, Sect. <ref> contains a discussion putting the results into the context of dust return to the ISM, an investigation of what constraints can be put on the grain sizes, and difficulties in reproducing the observations and determining the grain properties. Our conclusions on the dust in AGB wind-ISM interaction regions are then summarised in Sect. <ref>.

§ DATA

The wind-ISM interaction regions around a sample of AGB stars and supergiants were observed as part of the MESS (Mass-loss of Evolved Stars) programme <cit.> with PACS <cit.> aboard Herschel Space Observatory <cit.>. The interaction regions were observed in scan maps in the blue (at ) and red (at ) filters. The emission in the infrared and the spatial resolution make it possible to determine the dust masses and dust temperatures in the wind-ISM interaction regions (Sects. <ref> and <ref>). The data were fully reduced by the MESS team and presented in <cit.>. We use the highest-level of the JScanam maps of the reduced data from the archive (data products are provided at varying quality levels). At  and  the pixel scales are 16 and 32, respectively. The absolute flux calibration uncertainty of the observations is 15%. For more details on the data properties, see C2012 and the Herschel/PACS documentation[https://irsa.ipac.caltech.edu/data/Herschel/docs/nhsc/pacs/pacs.html].

C2012 classified the observed structures based on their geometrical appearance into "fermata", "rings", "eyes", and "irregular". Here we investigate the subsample of sources classified as "fermata" (a central source surrounded by a bow-shape, e.g.,  Fig. <ref> or <ref>) and "rings" (e.g., Fig. <ref>). All the data included in this study is presented in Appendices <ref> and <ref>. The origin of the structures for these sources is likely dominated by the interaction of a spherical wind with the ISM, while the "eyes" and "irregular" shapes may have additional shaping mechanisms (e.g. interactions with companions). The sources included in this study, with some basic parameters, are listed in Table <ref>.

§ METHODS

We use two methods to determine the temperatures and dust masses in the wind-ISM interaction regions: direct estimates from the data, and radiative transfer models. In the first the temperature is directly obtained from the ratio of the observed fluxes and assumed opacity values for the grains (Sect. <ref>), while in the second method we explicitly model the wind-ISM interaction regions using 3-dimensional dust-radiative transfer to consistently derive the temperatures (Sect. <ref>).

§.§ Direct estimates of temperatures and masses from the FIR images

In the far-infrared (FIR), dust grains emit thermal radiation as a black body with temperature  , modified by the dust opacity.

\begin{equation}
\label{e:flux}
F_{\lambda}= {{M_d B(\lambda,T_{\rm{d}})\kappa_{\lambda}} \over {D^2}},
\end{equation}

where $M_d$ is the dust mass, $B(\lambda,T_d)$ is the blackbody flux at the dust temperature $T_d$ at wavelength $\lambda$, $\kappa_{\lambda}$ is the absorption coefficient at wavelength $\lambda$, and $D$ is the distance. The ratio of observed fluxes ($=$/) in the far-infrared hence effectively constrain the (average) dust temperature and total dust mass, assuming optically thin dust emission and specific grain properties (that is, $\kappa_{\lambda}$ for a particular grain). The PACS observations can therefore be used to directly determine the dust temperature and mass in the observed wind-ISM interaction regions.

We measure the total flux towards the wind-ISM interaction regions in the  and  filters after subtracting a local background outside the wind-ISM interaction regions, and summing the flux above the rms value (determined in the background regions). The star and present-day mass-loss dominate the flux along the line-of-sight towards the central star, and these regions are masked out. The error in the individual measurements of the total flux is $\sim$1%. In some cases the rms cut-off is not sufficient to exclude emission outside of the wind-ISM interaction structure, and does not exclude the contribution from background sources. In these cases the flux is determined by hand within a polygon surrounding the main structure. By comparing the change in the measured flux using different polygon areas, we find that while the manual selection generally does not affect the fluxes measured in the  images (by less than 10%), this can have a much stronger effect on the measured  fluxes. Adding the absolute calibration uncertainty of 15%, we estimate the total uncertainties in the measured flux values to be $\approx$18% in the  images and $\approx$25% in the  images. The measured  and  fluxes, and the resulting  are shown in Table <ref>.

In all cases we assume optically thin dust emission from spherical grains with a constant radius of =0.1. This is comparable to what is typically assumed in AGB CSEs <cit.>, and allows to directly compare to C2012. However, see Sect. <ref> for s discussion on grain sizes. For oxygen-rich sources (that is, M-type and S-type stars) we use the absorption and scattering coefficients for astronomical silicates <cit.> with a grain density of 3.3 g/cm$^3$, giving a $\kappa_{\lambda}$ of 71.08 cm$^2$/g and 12.71 cm$^2$/g at  and , respectively, calculated using Mie theory. For carbon-rich sources we use the absorption and scattering coefficients for amorphous carbon grains <cit.> with a grain density of 1.8 g/cm$^3$, giving  a $\kappa_{\lambda}$ of 86 cm$^2$/g and 26 cm$^2$/g at  and , respectively. The measured fluxes and derived masses and temperatures ( and , respectively) are presented in Table <ref>. The uncertainties in the flux measurement lead to uncertainties of $\approx$60% in  and $\approx$15% in . Note that this method does not take any spatial constraints into consideration. The results merely provide the temperature for which the input dust grains will produce the observed ratios in the FIR. Different dust grains will produce the same FIR ratio, albeit at different temperatures. The actual temperature of the grains in the observed wind-ISM interaction regions depends on the stellar radiation field (the spectral energy distribution (SED) of the star), the effect of the present-day wind on the SED, the distance between the wind-ISM interaction region and the star, and on the grain properties. Using the spatial information from the PACS observations, and modelling the observed structures and radiation fields explicitly, hence provides additional constraints on the dust properties in the observed regions.

Estimated masses and temperatures in the wind-ISM interaction regions based on the measured fluxes in the PACS  and  images and assuming  = 0.1 . The observed far infrared-ratio is given by $F_{\rm{70\,\mu m}}^{im}$/$F_{\rm{160\,\mu m}}^{im}$. All corresponding PACS images and radiative transfer models are shown in Figs. <ref> to <ref>.

Source     Mask      $F_{\rm{70\,\mu m}}^{im}$     $F_{\rm{160\,\mu m}}^{im}$         $\Delta$    $^{im}$

[]     [Jy]      [Jy]             [$10^{-5}$ ]     [K]

2lCarbon dust

AQ And     10     3.2        1.8       1.8     0.5     23.3      38

U Ant      20     22.0       6.2       3.6     1.1     3.9       52

UU Aur     40     6.0        3.3       1.8     0.6     6.5       39

U Cam*      30     4.3        4.5       1.0     0.3     23.7      32

RT Cap     40     1.0        1.4       0.8     0.2     4.0       30

S Cep*      20     2.5        1.0       2.5     0.8     2.1       44

Y CVn*      30     2.2        11.8      0.2     0.1     126.2     22

TT Cyg     10     2.6        1.3       1.9     0.6     6.6       40

U Hya      30     28.9       14.7      2.0     0.6     10.2      40

CW Leo*     0      3459.6     598.6     5.8     1.8     42.8      70

W Ori      45     1.2        0.3       2.0     0.6     0.3       55

W Pic      20     1.2        0.3       1.4     0.4     0.6       55

TX Psc*     10     5.8        1.9       3.1     1.0     1.5       48

R Scl      10     28.1       5.8       4.8     1.5     4.6       62

S Sct      20     17.3       14.9      1.2     0.4     62.8      33

X TrA*      20     7.6        7.4       1.0     0.3     13.3      32

2lSilicate dust

$\theta$ Aps     35      6.6      1.4      4.8     1.5     0.4      46

W Aql            30      9.6      3.3      2.9     0.9     7.6      43

EP Aqr           20      8.2      2.9      2.9     0.9     1.9      38

R Cas            45      19.9     6.3      3.2     1.0     3.6      39

S Cas*            20      2.6      2.2      1.2     0.4     35.4     29

$\mu$ Cep        45      46.8     6.2      7.5     2.3     13.0     59

o Cet*            30      72.4     18.4     3.9     1.2     4.2      43

$\chi$ Cyg*       20      7.6      2.8      2.7     0.8     2.5      37

X Her*            20      10.4     3.4      3.1     1.0     2.2      39

R Hya            55      11.3     5.5      2.1     0.6     3.8      34

W Hya            50      32.1     8.6      3.7     1.2     2.7      42

R Leo            50      20.7     3.3      6.4     2.0     0.3      53

T Mic*            20      7.4      2.4      3.1     0.9     3.7      39

$\alpha$ Ori*     200     77.9     31.3     2.5     0.8     23.3     36

X Pav            30      9.3      4.5      2.1     0.6     16.5     34

V1943 Sgr        30      6.4      3.3      2.0     0.6     7.3      33

NML Tau          55      10.3     4.8      2.1     0.7     15.0     34

RT Vir           50      4.9      3.5      1.4     0.4     4.9      30

5l*Sources not modelled with radiative transfer (Sect. <ref>).

§.§ Radiative transfer models with RADMC-3D

In addition to deriving a first estimate of the masses and temperatures using the total PACS fluxes, we therefore also model the main structures in the wind-ISM interaction regions. The wind-ISM interaction regions do not have a simple spherical symmetry making the use of 3-dimensional models necessary. We calculate models using v2.0 of the 3D radiative transfer code RADMC-3D <cit.> RADMC-3D calculates the full frequency-dependent dust-radiative transfer using the Monte Carlo method, and produces images and spectra. Models can be calculated for any arbitrary 3D density distribution, allowing us to model the wind-ISM interaction regions explicitly, and to constrain the sizes and densities in the regions.

We focus here on the main observable structure of the bow-shock. The bow-shock is formed when the stellar wind collides with the surrounding ISM, and is affected by the velocity of the wind relative to the ISM and the density of the ISM <cit.>. It is assumed that the stellar wind expands freely within the bubble created by the bow-shock, and the position of the star and present-day wind with relation to the bow-shock depend on the angle of the star's space motion with respect to the local ISM.

In principle, the 3D surface of the bow-shock can be described by a hollow paraboloid <cit.>. This surface is derived assuming a ram pressure balance between the stellar wind and the flow of the ISM and the conservation of momentum flux across the shell. The geometry is described by the stand-off distance between the star and the wind-ISM interaction region. Using this description for the bow-shock surface has the advantage that it depends on physical parameters (the present-day mass-loss rate, the stellar wind velocity, the density of the surrounding ISM, and the stellar space velocity), which in principle can be constrained by fitting the observed structures. However, we find that the observed structures are generally better described by ellipsoids rather than wilkinoids, in particular at large angles relative to the direction of the apex of the bow shock (Sect. <ref> and Fig. <ref>).

In order to estimate the masses and densities in the wind-ISM interaction regions, we therefore approximate the 3D structure of the interaction regions with ellipsoidal shells that follow a Gaussian density distribution in the radial direction,

\begin{equation}
\rho_{\rm{IA}}(r)\propto exp({{-(r-r_{\rm{IA}})^2} \over {2\sigma_{\rm{IA}}^2}})\times [-cos({\theta \over 2})^A],
\end{equation}

where $\rho_{\rm{IA}}$ is the density in the wind-ISM interaction region,  is the distance from the star where the wind-ISM interaction region has its peak density, $\sigma_{\rm{IA}}$ is the width of the Gaussian distribution, $\theta$ is the angle to the apex of the bow shock, and $A$ is an exponent that allows the density to taper off behind the star. The density is scaled so that the integral over the entire volume of the wind-ISM interaction region gives the total mass . The distance from the star is given by the ellipsoid

\begin{equation}
r_{\rm{IA}}= {{R_0 \times R_{\rm{p}}} \over \sqrt{{R_0^2 cos(\theta) + R_{\rm{p}}^2 sin(\theta)}}},
\end{equation}

The parameters of the ellipse shape describing the wind-ISM interaction regions: $R_{\rm{0}}$ (stand-off distance), $R_{\rm{p}}$ (the distance perpendicular to $R_{\rm{0}}$), and $\theta$ (the angle to the apex of the bow shock). The star is assumed to move towards the right relative to the ISM. The black dashed line shows the shape of the wilkinoid for the same stand-off distance (see text for details).

where $R_{\rm{0}}$ is the distance between the apex of the bow-shock and the star (the stand-off distance), and $R_{\rm{p}}$ is the distance perpendicular to $R_{\rm{0}}$ (Fig. <ref>). Depending on the geometry, $R_{\rm{0}}$ can either be the major or minor axis of the ellipsoid (with $R_{\rm{p}}$ then being the respective other axis). We construct "open" ellipsoids by setting the density to 0 at angles $\theta>\theta_{\rm{lim}}$. For most sources $=100$°, with the exception of R Cas ($\theta_{\rm{lim}}=140^{o}$) and W Hya (a closed ellipsoid). We use a mild tapering off with $A=1 - 2$. Only for $\theta$ Aps the shell is more strongly tapered with $A=4$. Note that the parameters , , and  are physical parameters. Their observed projection will change depending on the inclination angle of the wind-ISM interaction region.

In addition to the wind-ISM interaction region, we also model the stellar radiation field and the present-day (dusty) mass-loss .  As the focus of this work is the modelling of the wind-ISM interaction regions, we use a simplified approach in modelling the stellar radiation field and present-day dust mass-loss to include a reasonable incident radiation field. We investigate the effect of changes in the incident radiation field in Sect. <ref>.
The central star is included assuming black body radiation with temperature  and luminosity . The density of the present-day mass-loss follows the standard 1/$r^2$ density profile, assuming the same dust expansion velocity as the gas. This is a simplified approach, as in reality the dust will have a drift velocity relative to the gas. Observationally this drift velocity can not be measured, as it is generally not possibly to measure the velocity of the dust. Theoretical models however predict significant drift-velocities, especially for low mass-loss-rate objects <cit.>.

The modelled dust mass-loss rates  are hence lower limits. However, the resulting density profile will still produce the required radiation field. The resulting SEDs are constrained by data downloaded from the Vizier online photometry viewer[http://vizier.u-strasbg.fr] and by varying the   until a satisfactory fit to the SED is found. The present-day dusty mass-loss is not strongly constrained  <cit.>, and for our purposes a reasonable fit to the data by eye is sufficient. The present-day wind expansion velocities  and gas-mass-loss-rates  are from C2012 (see references therein). The final , , and  are presented in Table <ref>.

The model grids are 400$\times$400$\times$400 with a step size of 2, and are converted to cm using the distances provided in Table <ref>. We tested models with a finer grid, which did not affect the results significantly for the wind-ISM interaction regions. In the central pixel the grid is refined to a 6$\times$6$\times$6 sub-grid to properly sample the increasing density of the present-day wind. The number of photons in the Monte Carlo simulation is set to 500000. We assume isotropic scattering by the dust grains. The details in the ray-tracing are described in the RADMC-3D manual.

Although RADMC-3D allows to model an arbitrary 3-dimensional structure, there are some limitations to what is feasible in this study. The description of the wind-ISM interaction region with an ellipsoid and a Gaussian density profile can not reproduce deviations from this smooth geometry, such as arcs and/or clumpy structures, or extended tails behind the star. Hence, we generally expect the models to not be able to reproduce all the emission that is observed in the images. Instead we focus on the main, overall structures of the wind-ISM region, and rather than fitting the total model flux to the observed flux, we try to reproduce the surface brightness and overall spatial distribution of the emission in the images. This fitting is done by eye by varying the radii, widths, angles, and dust masses in the wind-ISM interaction regions (see the Figs. <ref> to <ref> to see the observed and modelled structures). In general, the radii ( and ) are constrained to within $\sim$10% of the width of , with being constrained to within $\approx$15%. Although we are only considering sources classified as "fermata" or "rings", the structures in some sources are either too weak, and/or too complex to be feasibly modelled. These sources are indicated with a * in Table <ref>. We initially model all sources assuming the same carbon and silicate grains as in Sect. <ref> with $=0.1$ . Generally we find it difficult to fit the emission in both the  and  filters (i.e. to reproduce the ratio between the two filters, , see Sect <ref>), and we therefore focus on reproducing the flux at 70 . Based on our experience while manually fitting the interaction regions, we estimate the uncertainty in the dust mass to be on the order of $\approx\pm$20% (similar to the uncertainty in the flux measurement). The uncertainty in the derived  is estimated to be $\approx\pm$10%. This may however be affected by the degree of clumpiness in the interaction region (making it more difficult to fit the surface brightness) and/or projection effects.

The results of the radiative transfer modelling of the wind-ISM interaction regions are given in Table <ref>. The density  refers to the peak density of the dust in the model in the wind-ISM interaction region. The model images are compared to the observations in Appendix <ref>.

Basic stellar and present-day wind parameters used as input in the modelling of the wind-ISM interaction regions. All models assume  = 0.1 .

Source     D                        g/d

[pc]     [K]     []     [$10^{-10}$ ]     [$10^{-8}$ ]               []

AQ And     825       2200      5200      1         65        6500      17.0

U Ant     268       2500      5000      0.3       2         667       4.0

UU Aur     341       2500      8000      2.5       27        1080      11.0

RT Cap     291       2500      2200      0.32      3.2       1000      8.0

TT Cyg     562       2825      2700      0.2       3.2       1600      4.0

U Hya     208       3000      3750      0.63      4.9       778       8.5

W Ori     377       2200      7600      1.5       23        1533      11.0

W Pic     512       2200      4000      0.5       30        6000      7.0

R Scl     360       2325      5000      6         30        500       10.5

S Sct     386       2500      4000      0.08      2         2500      4.0

O-rich and S-type

$\theta$ Aps     113       2500      3500      0.5       11        2200      4.5

W Aql     340       1350      6000      100       1300      1300      20.0

EP Aqr     135       2500      3000      3.1       31        1000      11.5

R Cas     127       1900      3000      12        120       1000      13.5

$\mu$ Cep     390       2000      35000     100       200       200       35.0

R Hya     118       2500      6500      1.5       16        1067      12.5

W Hya     104       2200      9600      4         7.8       195       8.5

R Leo     71        2000      1800      0.92      9.2       1000      9.0

X Pav     270       2000      8000      10        52        520       11.0

V1943 Sgr     197       2300      6000      0.43      13        3023      5.4

NML Tau     245       1250      7000      175       320       183       18.5

RT Vir     136       2250      2000      3         50        1667      7.8

Results from the radiative transfer modelling with RADMC-3D assuming  = 0.1 . For R Scl and S Sct  is the total (detached shell $+$ wind-ISM interaction region) dust mass. The value in parentheses is the dust mass in the wind-ISM interaction region only. is the FWHM of the wind-ISM interaction region. The other parameters are explained in the text.

Source                             $F_{\rm{70\,\mu m}}^{mod}$     $F_{\rm{160\,\mu m}}^{mod}$

[$10^{-5}$ ]     []     []     []     [K]       [$10^{-24}$]     [Jy]      [Jy]      [10$^3$ yr]     [10$^3$ yr]

AQ And     25.0      54        54        7         38        1.0       3.1       1.8       12        2500

U Ant     3.0       43        43        4         62        9.4       28.9      5.9       14        1000

UU Aur     2.0       90        100       8         47        0.7       4.8       1.6       13        80

RT Cap     1.8       97        97        8         39        0.4       2.0       1.1       17        563

TT Cyg     3.6       35        35        2         47        3.7       2.9       1.0       23        1800

U Hya     6.0       120       120       10        47        2.1       30.3      10.3      14        952

W Ori     1.2       98        98        15        43        0.1       1.5       0.6       16        80

W Pic     1.2       52        46        5         48        0.8       1.1       0.4       18        240

R Scl     4.3 (1.3)     52        61        10        65        0.6       29.4      5.6       8         72

S Sct     18 (9.0)     132       132       10        42        2.6       13.4      5.9       61        22500

O-rich and S-type

$\theta$ Aps     0.2       75        67        7         52        5.6       5.3       0.9       9         40

W Aql     5.0       49        56        7         37        6.8       5.5       2.1       4         5

EP Aqr     0.4       37        50        7         56        14.4      11.7      1.7       2         14

R Cas     4.0       138       142       14        41        6.3       18.5      5.3       6         33

$\mu$ Cep     30.0      125       100       14        44        3.3       25.6      6.2       7         30

R Hya     0.8       94        125       15        48        1.7       11.9      2.3       4         50

W Hya     0.5       90        75        10        56        4.7       26.2      3.8       5         13

R Leo     0.3       110       118       13        45        5.2       9.5       2.1       4         33

X Pav     3.0       55        75        5         44        6.1       7.1       1.7       6         30

V1943 Sgr     1.0       65        80        10        47        2.1       4.8       1.0       11        233

NML Tau     8.0       95        95        10        35        5.5       12.0      5.1       6         5

RT Vir     0.6       85        120       10        39        1.9       3.2       1.0       7         20

§ RESULTS

§.§ Direct estimates vs. radiative transfer modelling

Figure <ref> shows the ratio between the dust masses and temperatures derived consistently in the radiative transfer modelling ( and , respectively; Sect. <ref>) to the dust masses and temperatures estimated directly from the flux-ratios in the images ( and , respectively; Sect. <ref>). Both the radiative transfer models and the direct estimates  assume 0.1  sized carbon or silicate grains. For the models the total dust mass is systematically less than the dust mass determined directly. On average, the ratio between  and  is 0.6, both for the C-rich and the O-rich sources. This difference is caused by the temperatures in the radiative transfer models being higher than those inferred directly from the flux density ratios. The consequence is that the dust mass has to be reduced in the model in order to still fit the emission at 70 . As mentioned in Sect. <ref>, determining the temperatures and dust masses from the ratios in the PACS images directly does not take spatial constraints into account. As a consequence, when including the spatial constraints in the models and calculating the temperatures solving the radiative transfer, the grains become too warm, and the models do not reproduce the observed . This discrepancy likely reflects issues in the assumed dust properties and/or radiation fields (Sect. <ref>).

Ratios between dust temperatures and dust masses obtained through two different methods using 0.1  grains: full radiative transfer ( and ) and derived directly from the PACS images ( and ) for C-rich sources (red dots) and O-rich and S-type AGB stars (blue squares).

§.§ Comparison to Cox et al. 2012

The data that is used here was first presented by <cit.>. They provided a classification for the different wind-ISM interaction regions (see Sect. <ref>) and determine the dust masses in the wind-ISM interaction similar to our direct method. <cit.> further discuss in detail the effect of the space motion and possible binarity on the observed structures, find an estimate of the density of the surrounding ISM, and run hydrodynamical models to investigate the morphology of the bow-shock and the formation of instabilities. While we make use of a number of the results in C2012 (e.g., the classification scheme, the reduced data, and the derived position angles and inclinations), there are some notable differences in the method and focus:

* C2012 measure the flux in apertures that follow the shape of the bow shock, within a certain radial range and angle, while we measure the entire flux in the image that is attributed to the wind-ISM interaction. As such, we generally measure more flux than C2012.

* We use the same formalism to derive the dust mass based on the flux in the images (Eq. <ref>). However, C2012 assume a fixed temperature of 30 K and determine separate dust masses based on the  and fluxes, respectively. We use the ratio of the measured fluxes to constrain the temperature and then determine one dust mass.

* C2012 use the same opacities for all sources <cit.>, while we separate between carbon dust and silicate dust.

* C2012 measure the stand-off distance (corrected for inclination) and use the description of wilkinoids to estimate the density of the surrounding .

* C2012 calculate hydrodynamical models of the wind-ISM interaction regions to study the morphology and instabilities. However, they do not perform radiative transfer models of the dust to compare with the observed images.

* C2012 derive estimates of the total dust mass in the wind-ISM interaction regions. However, they do not derive densities in the wind-ISM interaction regions.

In comparison to C2012, we focus on modelling the (typical) dust masses and dust temperatures, explicitly taking the shape and size of the wind-ISM interaction regions into account in 3D-radiative transfer models, and investigate the constraints that can be set on the dust properties.

Comparing our dust masses derived using Eq. <ref> to the dust masses from C2012 (note that C2012 presented a corrigendum to their original results), our masses are systematically lower at only $\approx$30% of their masses measured at 70, despite our generally larger fluxes. However, using our measured fluxes and fixing the temperature to T=30 K, our masses using Eq. <ref> increase by a factor 3.9 to an average of 17% larger than C2012. Additionally, using the opacities that C2012 use, the masses based on our measured fluxes increase by an additional factor of 1.25 to an average of 47% larger compared to C2012. From Table <ref> it is clear that the temperature in the dust emitting regions deviates significantly from 30 K (between 29 K to 62 K, with an average of 40 K) for several objects when assuming 0.1 -sized grains. Hence, although we in principle probe more of the mass in the wind-ISM interaction region by measuring the flux in the entire image (as opposed to using a limited aperture), by deriving the temperature instead of assuming one temperature, the total estimated mass decreases significantly compared to C2012. It is also worth noting that C2012 derive significantly different masses in the  and  images, indicating the same discrepancy in reproducing the observed  (Sects. <ref> and <ref>).

The other parameter we can directly compare is the stand-off distance between the star and the bow-shock. In C2012 this corresponds to the radius that defines the shape of the wilkinoid, while here it is the major or minor axis of the ellipsoid that is used to model the wind-ISM interaction region. C2012 give two radii, one that is predicted based on the stellar space-velocity, the stellar mass-loss rate, the gas velocity of the stellar wind, and the density in the surrounding ISM, and a radius measured in the observations (de-projected for the inclination angle). Their observed stand-off distances for the objects studied here are a factor of $\approx$2 smaller than what is predicted from the stellar space-velocity. The radii derived in our models are consistent with the observed distances in C2012 (and the distances measured in the images). In order to investigate the formation and evolution of the interaction regions, C2012 present hydrodynamical models of the wind-ISM interaction. They find that the predicted stand-off distance is consistent with the models. However, the shape of the interaction region deviates significantly from what is predicted by the wilkinoid in most models. Depending on the density contrast between the stellar wind and ISM, and the velocities of the wind, dust, star, and ISM, instabilities can form that affect the overall geometry of the interaction regions, and the structure within the regions (e.g., width and clumpiness; see C2012 for a detailed discussion). The discrepancy emphasises the importance of modelling the interaction regions explicitly using hydrodynamical models <cit.>.
In this context, it is worth noting that we generally can not reproduce the observed shapes in the wind-ISM interaction regions using wilkinoids in our 3-dimensional radiative transfer models – for any given stand-off distance, the wilkinoid becomes too wide at increasing angles $\theta$ (Fig. <ref>). While the ellipsoid manages to reproduce the overall structures, it naturally fails to produce the extended tails that develop behind the stars. However, up to an angle of $\approx$100° (relative to the direction of the stand-off distance) we find that ellipsoids follow the shape of the interaction regions better than wilkinoids. For the purpose of this paper, it is important to place the observable dust in our models at the correct distance from the star. As such, our ellipsoids are ad hoc descriptions of the observed structures without underlying physical processes, but allowing us to model the dust radiative transfer of the observed structures.

§ DISCUSSION

§.§ Structure and densities

The set of PACS observations of the wind-ISM interaction regions around AGB stars show a variety of morphologies (C2012). The different morphologies can generally be explained with the interaction of the stellar wind with the ISM, and/or the interaction between different periods of the stellar wind with varying mass-loss rates and expansion velocities. A full description of the structures requires hydrodynamical modelling of the winds and the surrounding material <cit.>. Keeping that in mind, the wind-ISM interaction regions modelled here can be overall well reproduced using an ellipsoid, and the derived stand-off distance is still connected to relevant physical parameters, such as the stellar- and wind-velocities, the mass-loss rate, and the density in the surrounding ISM.

The derived full-width-at-half-maximum (FWHM) of the interaction regions,   is 10%–15% of the stand-off distance. This is consistent with the detached shells around carbon stars created by wind-wind interaction <cit.>. Note, however, that for some objects a clear distinction can be made between the wind-wind interaction due to a thermal pulse and wind-ISM interaction. The carbon stars U Cam, R Scl, S Sct show clear wind-ISM interaction regions at further distances from the star than the spherical detached shells, indicating a predominantly radial wind-wind interaction due to a thermal pulse inside the wind-ISM interaction region. In the case of spherical shells, the origin of the observed structures may hence be due to wind-ISM interaction, and/or wind-wind interaction following a thermal pulse.

<cit.> present a simple model in order to describe the properties of the detached dust shell around the carbon star Y CVn. In their model the interaction with the ISM consists of a spherical shell with external material swept up from the ISM by the bow-shock, separated from circumstellar material by a contact discontinuity <cit.>. They estimate the total mass of the shell to be dominated by circumstellar material, the mass swept up from the ISM only constituting a few percent of the total mass.

In order to get a rough estimate of the timescales and densities for the sources discussed here, we assume a simplified approximation of the wind-ISM interaction regions along the lines of the work by <cit.>. For each source we assume a spherical half-shell with a radius  and width , and estimate the amount of mass that could have been swept up by the surrounding ISM in the corresponding volume

\begin{equation}
\label{e:volume}
\end{equation}

The total mass from the ISM swept up by this volume is estimated using the densities of the surrounding ISM, ,

\begin{equation}
$\Mism$=$\rhoism$ \,$\Via$=\mu_{\rm{H}}\,m_{\rm{H}}\,n_{\rm{H}}\,$\Via$,
\label{e:ismmass}
\end{equation}

where $\mu_{\rm{H}}=1.4$ is the mean nucleus number per hydrogen atom, $m_{\rm{H}}$ is the mass of the hydrogen atom, and $n_{\rm{H}}$ is the hydrogen density in the ISM. The surrounding $n_{\rm{H}}$ can be calculated assuming a relation between the particle density and the distance to the galactic plane $z$ <cit.>,

\begin{equation}
n_{\rm{H}}(z)=2.0\times\rm{e}^{-{{|z|} \over {100\rm{pc}}}},
\label{e:ismdens}
\end{equation}

where $z(\rm{pc}) = \emph{d}\, \rm{sin}(\emph{b}) + 15$, and $d$ is the distance and $b$ the Galactic latitude. Assuming a gas-to-dust ratio in the ISM of 1000, we find that the amount of swept-up mass in the wind-ISM interaction regions is less than 1% of the observed mass.

Finally, if we assume that the stars have been losing mass at constant rates  and constant velocities , we can calculate the time  it takes the stellar wind to build up the mass in the interaction region, and the crossing time  (the time it takes for the material in the stellar wind to reach the interaction region):

\begin{equation}
$\tia$=$\Mia$ / $\Mdotdust$,
\end{equation}

\begin{equation}
$\tcross$=$\Rso$ / $\vterminal$.
\end{equation}

Table <ref> shows the resulting build-up and crossing times. Note that both estimates are very uncertain because of uncertain present-day dust mass-loss rates  (which are not well constrained in the SED models), and the assumption that the dust has (and retains) the same velocity as the gas (Sect. <ref>). The estimated crossing times and build-up times are therefore only upper limits. However, for most sources  is significantly less that (changes in the dust velocity would affect both to the same degree), and the dust particles remain in the wind-ISM interaction regions on time-scales of a few 10000 –100000 years. The timescales are hence consistent with the build-up of the mass in the interaction region during AGB evolution. The average estimated gas density $n_{\rm{IA}}^{H}$ is on the order of a few 1000 $\rm{cm^{-3}}$ (assuming purely atomic hydrogen and a gas-to-dust ratio of 1000), corresponding to the densities in regular circumstellar envelopes at distances of a few 1 000 AU from the star. The densities are consistent with hydrodynamical models of interaction regions by <cit.>, and are consistent or slightly higher than the densities predicted by <cit.> and <cit.>. The derived densities are comparable to the number densities found in giant molecular clouds, dark clouds, and star forming clumps ($10^3 - 10^5$ cm$^{-3}$), and three to four orders of magnitude larger than the cold and warm ISM <cit.>.

It must be kept in mind that these estimates are based on over-simplifying assumptions. However, the estimates give an order-of-magnitude indication of the parameters determined, and indicate that the wind-ISM interaction regions contain significant amounts of dust and gas that is consistent with the build-up dominated by material from the AGB stars. The mass resides in the shell for a significant period of time (several 10 000 yr to a few 100 000 yr) at relatively high densities and cold temperatures.

§.§ Constraints on grain sizes

Although comparatively high, the derived densities when taken at face value are likely not high enough and the timescales are too short to allow for significant further grain growth <cit.>. However, observations in the submm of thermal emission from the dust in detached shell sources indicate the presence of relatively large grains (0.1  –2.0 ) compared to what generally is assumed in the inner AGB wind <cit.>, indicating possible continued growth and processing of the dust grains in the wind-wind interaction. The densities in the shells are an order of magnitude larger than what is derived here, but the timescales are significantly shorter (a few 1000 years) compared to the time the dust resides in the wind-ISM interaction. It is possible that increased densities in clumps caused by instabilities, and/or increased densities in connection with the shock front allow grains to grow to larger sizes.  Additionally, the models using =0.1  generally do not reproduce the observed  because the temperatures of the grains are too high (Sect. <ref>).

We therefore investigate the constraints observations in the FIR can provide on the grain sizes, and whether different grain sizes can reproduce the observed . For a given grain type, the ratio between the fluxes at  and   depends on the grain size, as the ratio between the absorption at optical and NIR wavelength and the emission in the FIR changes, effectively changing the grain temperature. We determine the FIR model flux ratios for the different grain sizes for each object, and compare the models to the observed ratios in the wind-ISM interaction regions, using opacities for = 0.01 , 0.1 , 0.25 , 0.5 , 0.75 , 1.0 , 2.0 , and 5.0 , for both amorphous carbon and silicate grains <cit.>.

Top panel: The observed values for  for C-rich sources (red dots) and O-rich sources (blue squares) and their uncertainty $\Delta$. The grey bars show the range of   that the models produce for grain sizes between 0.01 and 5.0. Bottom panels: Examples of the model  vs. grain size for the carbon star UU Aur (left) and the O-rich star V1943 Sgr (right). The horizontal solid lines show the observed , the dashed horizontal lines indicate the uncertainty $\Delta$.

Figure <ref> shows the range of modelled   for the different grain sizes for each source compared to the observed values, separated into C-rich and O-rich objects, as well as examples of the modelled   vs. grain size for a C-rich and O-rich object, respectively. For the sources with amorphous carbon grains, the observed  constrains the grain sizes for grains larger than $\approx$0.75–1.0 , while for smaller grains there is an ambiguity between sizes $<0.25$  on the one hand and 0.25  to 0.75  on the other (that is, grains smaller than 0.25   and grains between 0.25   and 0.75   produce the same ). For silicate grains the ambiguity between small and larger grains is even larger.

For the carbon-rich AGB stars, the observed ratios indicate relatively large grains ($\gtrsim$1 ). Note that we are only modelling grains of one constant size, while a distribution in grain sizes would be more realistic. However, given the observational information, such a distribution would be difficult to constrain. The effects of using different optical constants were investigated in the case of detached shells by <cit.>. They do not find a significant difference between different optical constants. Taken at face value, the results here indicate that a significant amount of dust must be present in relatively large grains in the wind-ISM interaction regions around the carbon AGB stars.

Although the uncertainties in the measured fluxes do not constrain the grain sizes for the C-rich objects very well, generally the models produce  values within the observed range. For the models with silicate grains for O-rich stars it is significantly more difficult to reproduce the observed . In half of the cases the observed  is a factor 1.5–2 lower than the smallest predicted ratio. In the remaining cases the uncertainty in the observed  is so large that the grain size cannot be constrained. This may partially be caused by bad measurements of the flux in the   images. However, comparing the observed values of  to the modelled values, we systematically over-predict   for models with silicate dust, while the models with amorphous carbon are comparable to the observed ratios. Unless this is a coincidence, it is more likely that the reason for the discrepancy in the O-rich sources lies in the details of the dust modelling.

§.§ Assumptions in the dust modelling and their effects on the modelled

The factors that determine   are the temperature of the grain, and the slope (and value) of the opacities in the FIR. In addition to the ratio between absorption in the optical and NIR and emission in the FIR, the temperature of the grain is also affected by the shape of the SED (that is, the incident light that the grain receives), in particular in the optical and NIR. This opens for (at least) two problems. The optical properties of the grains are not measured but are extrapolated at $\lambda>100$ . Failing to reproduce the observed ratios in the FIR may hence be due to unknown dust properties. A second issue is that we assume a stellar black body with a given luminosity and effective temperature in all our models. However, molecular and atomic absorption in the optical and NIR may significantly reduce the SED at short wavelengths compared to a plain black body.

In order to test to what extent these different mechanisms affect the observed , we construct a test model with representative parameters for the observed sources and wind-ISM interaction regions (Table <ref>). Figure <ref> shows model FIR ratios derived for our test model. The models are identical in their geometrical setup and density distribution, and are calculated using the standard silicate grains and amorphous carbon grains as used for all other sources. We additionally artificially changed the slope of the silicate grains to follow a $\lambda^{-1.5}$ law instead of $\lambda^{-2}$ at $\lambda$ > 100 , more similar to the amorphous carbon grains (corresponding to a direct change in dust properties). We also calculated models where we reduced the stellar black body to 25% at $\lambda$ < 2  (corresponding to the molecular and atomic absorption in the star's atmosphere). Both decreasing the stellar SED, and changing the grain properties at long wavelengths, would change the model  to be consistent with observations.

Parameters for the test model. For all test models the provided values were left fixed, and only grain properties and inclination angles were changes (see text for details). The values for , , and  are given for 0 inclination.

Parameter        Value

Luminosity        6000

Effective temperature        2300 K

Distance    D         200 pc

Present-day mass-loss

Dust mass-loss rate        $5\times10^{-9}$

Wind expansion velocity       10

Wind-ISM interaction

Dust mass       $10^{-5}$

Stand-off distance        65

Perpendicular size       80

Width         10

Opening angle     $\theta_{\rm{lim}}$     100$^{o}$

Density taper      A     1

Peak density         2.0$\times10^{-24}$

<cit.> recently summarised lab measurements of opacities in the FIR of silicate dust grains. They found that compared to the extrapolation from shorter wavelengths the measured opacities have values that are larger by a factor of $\approx$ 10. Significantly larger values of the opacity at FIR wavelengths would strongly affect the dust temperature, and hence . In order to test this, we increased the values of the opacities of our silicate dust at $\lambda$ > 50  by a factor of 10, and again changed the slope of the opacities to follow a $\lambda^{-1.5}$ law, making them similar to the grains in <cit.> (and references therein). The result on   is shown in Fig. <ref> and has by far the strongest effect on the observed ratios. However, the grains summarised in <cit.> are representative of dust observed in the cold ISM. They are formed under much colder conditions (20 – 100  K) than is typical for AGB stars (where the dust forms at $\sim$1000 K), and they contain a mix of silicate and carbon dust. Whether these opacities are representative of the dust that is observed in the wind-ISM interaction regions is hence unclear.

The silicate grains here include Fe, which affects the opacities at short wavelengths and hence the temperature of the grains. The composition of the dust around O-rich AGB stars is debated. Generally, grains that include Fe form too far from the star in order to accelerate the stellar wind, since they would be too warm to form in the acceleration region (however, they may still form further out in the wind). An alternative explanation for the driving of the stellar wind from O-rich AGB stars is the radiation pressure through scattering of photons off relatively large ($\approx$0.5 ) iron-free silicate grains <cit.>. Such grains would decrease the temperature, making them more consistent with the observed .

Finally, <cit.> investigated the effect of different grain properties on the opacities for silicate, amorphous carbon, and hydrogenated carbon grains. Effects that they studied were porosity, shape (spherical vs. oblate), size, and aggregates, and combinations of these parameters. Depending on the parameters changed, and the different combinations, opacities can change significantly from the optical to the FIR, hence also likely affecting the temperature of the grains and the observed emission at different wavelengths.

The detailed properties of the grains are hence very uncertain, and the derived values are strongly affected by a combination of our assumptions on the composition
(e.g. Fe-free silicates), the stellar SEDs (to account for the absorption in the atmosphere at short wavelengths), the structure (spherical vs. oblate, porous, aggregates), and the sizes. The FIR observations here alone are not able to constrain the dust properties in the wind-ISM interaction regions around AGB stars. Additional spatially resolved spectroscopic and imaging observations in the optical, near- and far-infrared of dust towards the wind-ISM interaction regions, combined with direct observations of the stellar SED, are necessary to constrain the radiation field, grain composition, and grain structures.

While it is not possible to derive tight constraints on the grain properties based on the observations here, taking our results at face value it appears that the grains in the wind-ISM regions may be different from the grains in the circumstellar envelopes around AGB stars. Generally, when explicitly modelling observations of the dust in winds around AGB stars, the same assumptions are often made as in this study: stellar black body, spherical solid grains, one grain size. This may be more easily possible since current observations do not provide high enough spatial resolution to constrain the dust distribution from the inner to the outer CSE, hence allowing more freedom in the dust temperature. Any effect described above that might affect the dust models in the wind-ISM interaction regions, would also affect the models of dust in circumstellar winds. However, these effects might influence the two regions differently, and whether the discrepancy between the CSE and wind-ISM interaction region remains needs to be determined through further radiative transfer modelling.

Comparison of the  for different grain sizes and different opacities: silicate grains (blue), amorphous carbon grains (red), silicate grains with $\lambda^{-1.5}$ slope at $\lambda$ > 100  (black-dashed), opacities that represent the opacities presented by <cit.>, and models with silicate grains but with a modified stellar black body (black-solid). The models are calculated for the test model. See text for details)

§ CONCLUSIONS

We investigated the far-infrared emission from wind-ISM interaction regions around AGB stars observed with the PACS instrument aboard Herschel at  and . We derive the temperatures and dust masses in the observed regions directly from the observed ratios between the two wavelengths (), as well as through explicit radiative transfer modelling of the observed structures, assuming amorphous carbon grains for C-rich AGB stars, and silicate grains for O-rich (M- and S-type) AGB stars.

When assuming grain-sizes of  = 0.1 , the models produce grains with relatively high temperatures and fail to reproduce the observed . For amorphous carbon grains models with  $\gtrsim$ 1  can explain the observed , indicating the presence of relatively large grains in the wind-ISM interaction regions. For O-rich and S-type AGB stars, the observed  can not be reproduced for most sources, irrespective of the grain size. In this case it is likely that a combination of grain properties (composition, structure, and size) and the detailed stellar SED affect the temperatures of the grains. We discuss several possible mechanisms that could explain the discrepancies between modelled and observed  for the O-rich sources:

* The slope in the opacities for silicate grains at $\lambda$ > 100 .

* A reduction in the incident SED at $\lambda$ < 2  as an effect of atomic and molecular absorption in the stellar atmosphere.

* The possibility of grains with the properties of cold interstellar grains at $\lambda$ > 50 .

* Possible contamination of the grains with Fe.

* Changes in the opacities owing to porosity, shape, size, and aggregates.

It is not clear which of these mechanisms dominate, and in particular the implications for the grain properties affect our understanding of the dust that is released into the ISM.

The derived dust masses in the wind-ISM interaction regions are consistent with the build-up of dust through the stellar mass-loss, with swept-up ISM dust only constituting $\approx$1% of the total observed dust mass. The indication of large grains is consistent with previous studies of the dust in detached shells, and indicates that dust particles from AGB stars continue to grow to significantly larger grains compared to what is derived for the inner AGB wind. However, the derived typical densities in the wind-ISM interaction regions do not seem to be high enough to allow for (obvious) grain growth, unless the interaction regions contain density enhancements with up to two orders of magnitude.

This has implications for the properties of  particles that are released into the ISM, where the AGB dust acts as seed particles for continued grain growth in the ISM. Generally it is either assumed that the dust grains released in to the ISM are unchanged from the original dust formed, or that they are sublimated by the interstellar radiation field and supernova shocks and do not contribute to the dust (growth) in the ISM at all. However, while we can not conclusively constrain the dust properties, our results indicate that assumptions that are often made to model the dust in the inner CSE do not successfully reproduce the observations of wind-ISM interaction regions. There is additionally a tendency towards larger grains compared to what is generally assumed. This increases the chance of grain survival, and affects the role dust grains from AGB stars play in forming the seeds for grain growth in the ISM. The detection of pre-solar grains from AGB stars is consistent with this picture.

While our models give rough estimates of the grain sizes and dust masses in the wind-interaction regions, and explicitly model the 3D structure of the regions, they are still very basic. In order to fully derive the conditions in the wind-ISM interaction regions, and the properties of the dust released into the ISM, it is necessary to obtain both spectroscopic and imaging observations of the dust, spatially resolved, in order to determine the grain compositions and temperatures, as well as direct observations of the stellar SEDs. Full hydrodynamical models of the interaction regions aimed at explicitly reproducing the observed cases are necessary to determine the physical parameters of the regions. Such investigations will be necessary to constrain models of grain growth and processing in the wind-ISM interaction regions to derive the detailed properties (size, composition, shape, porosity) of the grains released into the ISM.

M.Maercker acknowledges support from the Swedish Research Council under grant number 2016-03402. T.K. acknowledges support from the Swedish Research Council under grant number 2019-03777. E.D.B. acknowledges funding by the Swedish National Spaceboard.

§ MODELS AND IMAGES FOR SOURCES THAT WERE MODELLED

AQ And: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 0.4$\times10^{-3}$ (70 ) and 0.2$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

U Ant: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 7.6$\times10^{-3}$ (70 ) and 1.2$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

UU Aur: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 0.75$\times10^{-3}$ (70 ) and 0.2$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

RT Cap: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 0.12$\times10^{-3}$ (70 ) and 0.056$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

TT Cyg: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 2$\times10^{-3}$ (70 ) and 0.26$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

U Hya: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 1.3$\times10^{-3}$ (70 ) and 0.38$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

W Ori: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 0.066$\times10^{-3}$ (70 ) and 0.025$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

W Pic: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 0.68$\times10^{-3}$ (70 ) and 0.16$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

R Scl: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 0.7$\times10^{-3}$ (70 ) and 0.2$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

S Sct: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 1$\times10^{-3}$ (70 ) and 0.35$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

$\theta$ Aps: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 1.4$\times10^{-3}$ (70 ) and 0.19$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

W Aql: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 1.8$\times10^{-3}$ (70 ) and 0.44$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

EP Aqr: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 3.5$\times10^{-3}$ (70 ) and 0.44$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

R Cas: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 1.2$\times10^{-3}$ (70 ) and 0.3$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

$\mu$ Cep: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 3$\times10^{-3}$ (70 ) and 0.58$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

R Hya: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 1.1$\times10^{-3}$ (70 ) and 0.17$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

W Hya: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 3$\times10^{-3}$ (70 ) and 0.35$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

R Leo: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 1.1$\times10^{-3}$ (70 ) and 0.2$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

X Pav: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 2.2$\times10^{-3}$ (70 ) and 0.34$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

V1943 Sgr: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 1.2$\times10^{-3}$ (70 ) and 0.19$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

NML Tau: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 0.6$\times10^{-3}$ (70 ) and 0.2$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

RT Vir: Top to bottom: The Radmc3D model, the PACS image, and the PACS image with contours from the model. Images are for 70  (left) and 160  (right). Maximum contour levels are 0.4$\times10^{-3}$ (70 ) and 0.1$\times10^{-3}$ (160 ), respectively. Minimum contour levels are 10% of maximum. The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

§ PACS IMAGES OF SOURCES THAT WERE NOT MODELLED

U Cam: The PACS images at 70  (left) and 160  (right). The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

S Cep: The PACS images at 70  (left) and 160  (right). The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

Y CVn: The PACS images at 70  (left) and 160  (right). The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

CW Leo: The PACS images at 70  (left) and 160  (right). The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

TX Psc: The PACS images at 70  (left) and 160  (right). The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

X TrA: The PACS images at 70  (left) and 160  (right). The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

S Cas: The PACS images at 70  (left) and 160  (right). The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

$\omicron$ Cet: The PACS images at 70  (left) and 160  (right). The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

$\chi$ Cyg: The PACS images at 70  (left) and 160  (right). The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

X Her: The PACS images at 70  (left) and 160  (right). The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

T Mic: The PACS images at 70  (left) and 160  (right). The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.

$\alpha$ Ori: The PACS images at 70  (left) and 160  (right). The colour scale is in . The red dashed circle shows the mask used to measure the flux from the star and present-day mass-loss.


# Towards Playlist Generation Algorithms Using
RNNs Trained on Within-Track Transitions

Keunwoo Choi

György Fazekas

Mark Sandler

Centre for Digital Music  EECS, QMUL  London, UK<EMAIL_ADDRESS>Centre for Digital Music  EECS, QMUL  London, UK<EMAIL_ADDRESS>Centre
for Digital Music  EECS, QMUL  London, UK<EMAIL_ADDRESS>
###### Abstract

We introduce a novel playlist generation algorithm that focuses on the quality
of transitions using a recurrent neural network (RNN). The proposed model
assumes that optimal transitions between tracks can be modelled and predicted
by internal transitions within music tracks. We introduce modelling sequences
of high-level music descriptors using RNNs and discuss an experiment involving
different similarity functions, where the sequences are provided by a musical
structural analysis algorithm. Qualitative observations show that the proposed
approach can effectively model transitions of music tracks in playlists.

###### keywords:

recurrent neural networks; music playlists; music recommendation

††conference: Workshop on Surprise, Opposition, and Obstruction in Adaptive
and Personalized Systems (SOAP) June 13–17, 2016, Halifax, NS,
Canada††titlenote:  This work was part funded by the “Fusing Audio and
Semantic Technologies for Intelligent Music Production and Consumption” (FAST
IMPACt) EPSRC Grant EP/L019981/1 and the European Commission H2020 research
and innovation grant AudioCommons (688382). Sandler acknowledges the support
of the Royal Society as a recipient of a Wolfson Research Merit Award.

<ccs2012> <concept>
<concept_id>10010147.10010257.10010293.10010294</concept_id>
<concept_desc>Computing methodologies Neural networks</concept_desc>
<concept_significance>500</concept_significance> </concept> <concept>
<concept_id>10010405.10010469.10010475</concept_id> <concept_desc>Applied
computing Sound and music computing</concept_desc>
<concept_significance>500</concept_significance> </concept> </ccs2012>

[500]Computing methodologies Neural networks [500]Applied computing Sound and
music computing

##  1 Introduction

In music recommendation, the quality of transitions become important
particularly when the recommendation is provided in the form of a playlist.
This is due to a unique aspect of music consumption. Unlike other products,
music is consumed i) instantaneously, for instance, while listening using
streaming services, ii) repeatedly, i.e., listeners are willing to listen to
the same music multiple times, and iii) quickly, i.e., an item usually only
lasts a few minutes. Hence, recommended items are typically consumed or played
in a sequence. This behaviour introduces the need for good transitions between
items, that is, the relevance and subjective judgement a recommended track
depends on the previous track.

Recommendation approaches using collaborative filtering are prone to overlook
niche or new items, although the popularity bias of known items can be
compensated for. This is called the cold-start problem [20]. Content-based
approaches which are designed to solve the cold-start problem can suffer from
lack of diversity when recommended items are selected simply by similarity.
This is often called top-$N$ recommendation. It is well known that
unexpectedness, surprise or serendipity play an important role in the music
recommendation and discovery [4]. Compared to other strategies, focusing on
transitions can naturally provide these qualities.

There have been approaches that primarily focus on the transitions of tracks
[13], [3], [15]. They assumed the Markov property of hidden states or
embeddings of tracks. Using the Markov property, it is assumed that future
events only depend on the current one and does not depend on the past. This
has been successfully used for sequence modelling for instance in speech [19]
too. In music computing, playlist datasets [16], [14], [3] collaboratively
created for reference by DJs and listeners were used for training and
evaluation of sequence modelling approaches. Although these datasets consist
of a large number of tracks, e.g. 101k playlists in [16], the lack of audio
data fundamentally limits research based on audio content analysis.

Recently, recurrent neural networks (RNNs) have become widely used for
sequence modelling in tasks such as speech recognition, substantially
outperforming previous hidden Markov model-based approaches [22]. The success
of the application of RNNs largely relies on the introduction of Long Short-
Term Memory (LSTM) units [10]. The merit of LSTM comes from the gate cells of
LSTM units, that decide how much the units take input, release output, and
forget the previous events. Especially, the forget gate improves the training
efficiency by helping the gradients flow well. However, RNNs have not been
used for playlists generation and modelling, due in part to the lack of
sufficient training data. To solve this problem, we propose using an RNN
trained on within-track transitions to model playlists.

We assume that transitions between structural segments of music can be used as
a model for generating the desired high-quality transitions between tracks. In
general, segments in a track are different but coherent and their musical
features can be expected to match well in succession. This is due to the
careful and intentional design by the composer. Using this approach, the
number of transitions can easily outnumber that of existing playlist datasets,
and therefore it enables to train an RNN model.

The rest of the paper is organised as follows. The proposed method is first
described in Section 2. We then present experimental results and discussion in
Section 3 and conclude in Section 4.

##  2 The Proposed Model

Figure 1: A block diagram of the proposed algorithm, (a,b) training of RNN and
(c) prediction of a feature vector, $x_{pred}$.

Figure 1 illustrates the procedure of dataset construction, as well as the
training and prediction stages of the proposed algorithm. First, the training
tracks are segmented and $x_{i}$, the features for each segments are extracted
(Fig. 1a). Then an RNN of length $N$ ($N$=3 in the figure) is trained to learn
the transitions of the sequence of feature vectors (Fig. 1b). When a seed
track is provided, the features of the last $N$ segments are extracted and fed
into the trained RNN to predict the feature vector $x_{pred}$ (Fig. 1c). The
algorithm selects a track with a start segment that is most similar to
$x_{pred}$.

###  2.1 Structural Segmentation

Structural segmentation is a task aiming to find the boundaries of different
segments or parts in music, e.g. intro, verse, bridge, chorus. The most common
approach is to take advantage of self-similarity between frames of the track
[9]. In the experiment, we used a basic and efficient method that is proposed
in [9]. Although the results introduce some errors, the feature vector
sequences that are based on the imperfect segmentation still approximate the
information about how each feature changes along time in each track.

###  2.2 Feature Extraction

The proposed algorithm can use feature extraction methods that are relevant to
listeners’ musical preferences and able to represent a musical segment. This
includes estimated latent features from collaborative filtering [26], tags
such as genre and emotion [7] or implicit features such as the weights of the
last hidden layer of a neural network classifier [12]. Using explicit labels
such as genre can facilitate explaining the behaviour of the algorithm, which
is important for research and also to the listener.

In the experiment, an auto tagging algorithm in [5] is used to predict a
50-dimensional vector whose elements correspond to the probability of each
tag. The tagging algorithm is based on deep convolutional neural networks and
trained on Million Song Dataset [2]. It shows state-of-the-art performance
while the tags cover variety of categories such as genre, emotion, instrument,
and era. Although some of the tags such as genre typically characterise the
entire music track, they are not necessarily constant over the whole track.

###  2.3 RNN Model

The goal of RNN training is to predict the feature of the following track
($x_{pred}$) that maintains consistency and fluctuations, i.e., a certain
variation over the features. To this end, a 2-layer RNN with 512 hidden units
is employed. LSTM units [10] are used as they show state-of-the-art
performance among RNN variants for several sequence modelling tasks [11].

###  2.4 Similarity Measure

A similarity measure is necessary to find the subsequent track using the
feature vector predicted by the RNN. The similarity metric directly affects
the properties of the generated playlists and therefore it should be carefully
selected. Using the cosine distance may compensate for the popularity bias and
result in recommending more niche items [23]. The Discounted Cumulative Gain
(DCG) turned out to be effective in our experiment.

DCG is a weighted version of Cumulative Gain (CG). CG is designed to measure
ranking quality of a retrieved list and DCG weights on the top-$N$ relevant
items by discounting lower relevant items. Applying this measure was motivated
by the type of feature extraction algorithm we use. The extracted feature is a
vector of probabilities of each tag and tags with large probabilities should
be weighted more than the others to facilitate maintaining the consistency of
the generated playlists. Because DCG weights high-ranking elements more, it
can theoretically work better.

##  3 Results and Discussion

###  3.1 Configurations

For training, we used a private dataset with 28,430 commercial tracks of
modern popular music including Rock, Hip-Hop, and Jazz. We used 7,880 tracks
from the ILM10k dataset for testing [21] [1].

The segmentation is performed using [9], which is implemented in [18]. As
mentioned in Section 2, an automatic tagging algorithm was used as a feature
extractor [5]. An RNN with a length of 50 is trained. We compared DCG with the
cosine distance and the $l_{2}$ norm for computing similarity. Audio
processing and RNN are implemented using librosa [17], Keras [6], and Theano
[25].

Figure 2: Transitions of feature vectors by cosine distance, $l_{2}$, and DCG
on left, centre, and right, respectively. Y-axis is time (top to bottom) and
x-axis refers to the 50 feature dimensions. The first track (topmost feature
vectors) is the seed.

Figure 2 shows the transitions of feature vectors in three playlists given the
same seed track but different similarity metrics. The seed track is
represented by a 7-by-50 matrix, i.e., 7 segments of the track. White
horizontal lines indicate beginnings and ends of each track. The predicted
feature vectors (1-by-50) are illustrated in between.

The figure helps to explain several aspects of qualitative observations by the
authors while listening to the generated playlists.

###  3.2 Discussions

First, we found both consistency and fluctuations in the extracted features
within tracks. In general, several features show consistently large (blue) and
small (red) values, while the other features vary. It supports the selection
of feature extraction algorithm.

However, there are still rooms for further improvements. For example,
whitening of each feature dimension can be adopted to compensate the prior
distributions of each feature. Although RNNs are able to adapt to such
differences, the similarity measure may be affected from such pre-processing.

Second, the transitions usually successfully keep the coherence within
playlists as demonstrated by the figure. However, we noticed that the model is
prone to missing overall similarity in long playlists. It may be related to
the observation that the trained RNN occasionally predicts a vector that does
not have near neighbours. It means the selected track is not similar enough to
the predicted vector, which may result in an undesirable or suboptimal
transition.

This may first be due to the short lengths of segments in the training data.
The majority of training tracks have fewer segments (90$\%$ of the training
tracks has less than 17 segments), therefore the long-term dependency may not
be learnt and the prediction may be dominated by short-term features.

In future work, training with longer sequences such as concatenated features
of tracks from sequences in an albums, setlists or curated playlists may be
used to help learning better transitions. It can be also a typical behaviour
of RNNs. Although RNNs generally model the long-term dependency of sequences,
in many cases, RNNs have shown a behaviour puts more emphasis on recent inputs
rather than older ones.

However, the problem seems to be partly resolved when DCG is used as the
similarity measure as discussed alongside out last set of observations.

Third, DCG provides more coherent playlists compared to the cosine distance
and $l_{2}$-norm. This phenomenon is found not only in the sequences
corresponding to Figure 2 but also in other playlists. This can be explained
as follows. In each track, there are consistently strong, consistently weak,
and fluctuating features. This pattern, especially the consistencies, can be
easily learned by the RNN and the consistent features are maintained in the
predictions. Finally, DCG prioritises the features that are large in
predictions, resulting in successfully finding a track with those consistently
large features. This improves the coherence of the resulting playlists.

##  4 Conclusion

In this paper, we proposed a novel algorithm for playlist generation that
relies on learning desirable musical qualities from within-track transitions
between musical segments. The proposed combination of an RNN, within-track
structure and DCG showed an encouraging result. Within-track structures showed
the consistency and dynamics that are assumed in playlists. An RNN model
learnt the feature sequences and its predictions are successfully used for the
selections of following tracks. Different similarity measures resulted in
different playlists.

Future work will investigate advanced architectures such as bidirectional RNNs
[24] and more formal assessments. Using bidirectional RNNs can be used to
create playlists that have more constraints e.g. start and end tracks [8] and
steerability [14]. Formal assessments will include subjective measurement e.g.
satisfaction and objective measures with regards to consistency, fluctuations
and diversity.

## References

  * [1] Allik, A., Fazekas, G., Barthet, M., and Saldler, M. myMoodplay: an Interactive Mood-based Music Discovery App.  In Proc. of the 2nd Web Audio Conference (WAC), April 4-6, 2016, Atlanta, Georgia. (2016).
  * [2] Bertin-Mahieux, T., Ellis, D. P., Whitman, B., and Lamere, P. The million song dataset.  In Proceedings of the 12th International Society for Music Information Retrieval Conference, ISMIR 2011, Miami, Florida, USA, October 24-28, 2011 (2011), pp. 591–596.
  * [3] Chen, S., Moore, J. L., Turnbull, D., and Joachims, T. Playlist prediction via metric embedding.  In Proceedings of the 18th ACM SIGKDD international conference on Knowledge discovery and data mining (2012), ACM.
  * [4] Choi, K., Fazekas, G., and Sandler, M. Understanding music playlists.  In Proceedings of the ICML 2015: Machine Learning for Music Discovery Workshop (2015), ICML.
  * [5] Choi, K., Fazekas, G., and Sandler, M. Automatic tagging using deep convolutional neural networks.  In Proceedings of the Conference on International Society of Music Information Retrieval, New York, USA (2016), ISMIR.
  * [6] Chollet, F. Keras: Deep learning library for theano and tensorflow.  https://github.com/fchollet/keras, 2015.
  * [7] Dieleman, S., and Schrauwen, B. End-to-end learning for music audio.  In Acoustics, Speech and Signal Processing (ICASSP), 2014 IEEE International Conference on (2014), IEEE, pp. 6964–6968.
  * [8] Flexer, A., Schnitzer, D., Gasser, M., and Widmer, G. Playlist generation using start and end songs.  In ISMIR (2008), pp. 173–178.
  * [9] Foote, J. Automatic audio segmentation using a measure of audio novelty.  In IEEE International Conference on Multimedia and Expo, (ICME) 2000. (2000), IEEE.
  * [10] Gers, F. A., Schmidhuber, J., and Cummins, F. Learning to forget: Continual prediction with lstm.  Neural computation 12, 10 (2000), 2451–2471.
  * [11] Greff, K., Srivastava, R. K., Koutník, J., Steunebrink, B. R., and Schmidhuber, J. Lstm: A search space odyssey.  arXiv preprint arXiv:1503.04069 (2015).
  * [12] Liang, D., Zhan, M., and Ellis, E. D. P. Content-aware collaborative music recommendation using pre-trained neural networks.  In Proceedings of the Conference on International Society of Music Information Retrieval (2015), ISMIR.
  * [13] Liebman, E., Saar-Tsechansky, M., and Stone, P. Dj-mc: A reinforcement-learning agent for music playlist recommendation.  In Proceedings of the Conference on Autonomous Agents and Multiagent Systems (2015).
  * [14] Maillet, F., Eck, D., Desjardins, G., Lamere, P., et al. Steerable playlist generation by learning song similarity from radio station playlists.  In Proceedings of the Conference on International Society of Music Information Retrieval (2009), ISMIR, pp. 345–350.
  * [15] McFee, B., and Lanckriet, G. R. The natural language of playlists.  In Proceedings of the Conference on International Society of Music Information Retrieval (2011), ISMIR, pp. 537–542.
  * [16] McFee, B., and Lanckriet, G. R. Hypergraph models of playlist dialects.  In Proceedings of the Conference on International Society of Music Information Retrieval (2012), ISMIR.
  * [17] McFee, B., Raffel, C., Liang, D., Ellis, D. P., McVicar, M., Battenberg, E., and Nieto, O. librosa: Audio and music signal analysis in python.  In Proceedings of the 14th Python in Science Conference (2015).
  * [18] Nieto, O., and Bello, J. P. Msaf: Music structure analysis framework.  In Proceedings of the Conference on International Society of Music Information Retrieval (2015), ISMIR.
  * [19] Rabiner, L. R. A tutorial on hidden markov models and selected applications in speech recognition.  Proceedings of the IEEE 77, 2 (1989), 257–286.
  * [20] Ricci, F., Rokach, L., and Shapira, B. Introduction to recommender systems handbook.  Springer, 2011.
  * [21] Saari, P., Fazekas, G., Eerola, T., Barthet, M., Lartillot, O., and Sandler, M. Genre-adaptive semantic computing and audio-based modelling for music mood annotation.  IEEE Transactions on Affective Computing (TAC) (2015), 1–17.
  * [22] Sak, H., Senior, A., and Beaufays, F. Long short-term memory based recurrent neural network architectures for large vocabulary speech recognition.  arXiv preprint arXiv:1402.1128 (2014).
  * [23] Schedl, M., Knees, P., McFee, B., Bogdanov, D., and Kaminskas, M. Music recommender systems.  In Recommender Systems Handbook. Springer, 2015, pp. 453–492.
  * [24] Schuster, M., and Paliwal, K. K. Bidirectional recurrent neural networks.  Signal Processing, IEEE Transactions on 45, 11 (1997), 2673–2681.
  * [25] Team, T. T. D., Al-Rfou, R., Alain, G., Almahairi, A., Angermueller, C., Bahdanau, D., Ballas, N., Bastien, F., Bayer, J., Belikov, A., et al. Theano: A python framework for fast computation of mathematical expressions.  arXiv preprint arXiv:1605.02688 (2016).
  * [26] Van den Oord, A., Dieleman, S., and Schrauwen, B. Deep content-based music recommendation.  In Advances in Neural Information Processing Systems (2013), pp. 2643–2651.