# Spiral arms in thermally stratified protoplanetary discs

Spiral arms in thermally stratified protoplanetary discs Abstract Spiral arms have been observed in nearly a dozen protoplanetary discs in near-infrared scattered light and recently also in the submillimetre continuum. While one of the most compelling explanations is that they are driven by planetary or stellar companions, in all but one cases such companions have not yet been detected and there is even ambiguity on whether the planet should be located inside or outside the spirals. Here, we use 3D hydrodynamic simulations to study the morphology of spiral density waves launched by embedded planets taking into account the vertical temperature gradient, a natural consequence of stellar irradiation. Our simulations show that the pitch angle of the spirals in thermally stratified discs is the lowest in the disc mid-plane and increases towards the disc surface. We combine the hydrodynamic simulations with 3D radiative transfer calculations to predict that the pitch angle of planetary spirals observed in the near-infrared is higher than in the submillimetre. We also find that in both cases the spirals converge towards the planet. This provides a new powerful observational method to determine if the perturbing planet is inside or outside the spirals, as well as map the thermal stratification of the disc. accretion, accretion discs, hydrodynamics, protoplanetary discs, circumstellar matter 1 INTRODUCTION With the dawn of extreme high resolution imaging capabilities of current telescopes (e.g. SPHERE/VLT, GPI/Gemini, ALMA) more and more discs show spiral arms in near-infrared (NIR) scattered light (e.g. Muto et al. 2012; Benisty et al. 2015; Wagner et al. 2015; Stolker et al. 2016; Benisty et al. 2017) and recently also in the sub-millimetre (Pérez et al. 2016; Tobin et al. 2016). The observed spirals have similar morphologies in most cases; two symmetric arms, shifted in azimuth by approximately 180°, and with a pitch angle of about 15°–30°. In most cases, there is also a gap/cavity just inwards of the spirals. It has been known for a long time that spiral density waves are launched by planetary or stellar companions at Lindblad resonances (e.g. Goldreich & Tremaine 1979; Ogilvie & Lubow 2002; Rafikov 2002) both inwards and outwards of companion’s orbit. Therefore, a straightforward explanation for the observed spirals is that they are launched by a so far undetected companion, possibly a forming embedded planet (e.g. Muto et al. 2012; Grady et al. 2013); this would open up the exciting possibility of getting a direct insight into planet formation. Using 2D hydrodynamic simulations and 3D radiative transfer computation, Juhász et al. (2015) showed, however, that the pitch angle and contrast of the spirals driven by planets orbiting inwards of the spirals are not compatible with the NIR observations. Dong et al. (2015a) combined 3D hydrodynamic simulation with 3D radiative transfer calculations to show that the observed contrast and pitch angle of spirals are in fact in agreement with the spirals being launched by massive (several MJup) planets exterior to the spirals. However, planets are not the only explanations for the observed spirals. Multi-armed spiral structure can also form due to gravitational instability (GI; see e.g. Rice et al. 2003) with the number of arms depending on the mass and temperature of the disc via the Toomre Q parameter (Cossins, Lodato & Clarke 2009). For very massive discs (Mdisc/M⋆ ≳ 0.25), GI can produce two armed spirals in the disc, similar to the NIR scattered light observations (Dong et al. 2015b). Meru et al. (2017) showed that GI could explain the spirals in the submillimetre ALMA images of Elias 2–27. Additionally, if a marginally gravitationally stable disc is perturbed by an embedded planet, the tidal interaction might tip the balance and trigger GI in the disc (Pohl et al. 2015). This combined mechanism may also produce spirals with pitch angles similar to the ones in observations. Numerous theoretical studies have recently increased the level of complexity in the modelling of spirals, showing the importance to understand their origin. For instance, while tidally induced spiral density waves have been studied for a long time in 2D, Zhu et al. (2015) demonstrated the spiral amplitude and pitch angle (the complement of the angle between the radial direction and the tangent of the spiral) in the non-linear regime are higher in 3D than in 2D simulations. However, an important aspect of a realistic disc model has not yet been investigated. Protoplanetary discs in which spirals have been detected are passively heated by the absorbed stellar radiation. As a consequence of this, they will have a positive vertical temperature gradient (Chiang & Goldreich 1997), resulting in a vertical thermal stratification in the disc. It is known that the pitch angle of the spirals depends on the local sound speed (Rafikov 2002), thus we expect vertical thermal stratification to affect the morphology of the spiral wake. (see e.g. Lee & Gu 2015). The goal of this study is to investigate the morphology of tidally induced spirals by embedded planets in 3D thermally stratified discs. We note, however, that our results would also be applicable to stellar companions as the tidal interaction between the disc and the companion is the same as for the planetary case. 2 MODELLING 2.1 Hydrodynamic simulations We run 3D hydrodynamical simulations with the code pluto (Mignone et al. 2007) to investigate the disc response to the planet. We include the FARGO algorithm (Mignone et al. 2012) to reduce the computational cost and use the HLLC Riemann solver, linear reconstruction and second-order Runge–Kutta time integration. We employ a spherical grid with 400, 50, and 1024 cells, distributed in the [0.4, 2.5], [1.22, π/2], and [0, 2π], intervals in the radial, polar, and azimuthal directions. The planet was placed at the coordinates [1.0, 0, π/4]. In the radial and polar directions, we define ‘wave-killing zones’ where we damp disturbances, following de Val-Borro et al. (2006), extending up to 0.44 and from 2.25 in radius and until 1.27 in θ. We employ a locally isothermal equation of state, in which the sound speed cs is a function of position only, and a physical viscosity using the Shakura & Sunyaev (1973) prescription with a value of α = 10−3. We assume the planet mass to be 3 MJup. We run two simulations, the first (subsequently called ‘isothermal’) in which cs depends only on radius: cs ∝ r−1/4, with the normalization set such that the disc aspect ratio at the location of the planet is 0.07, a typical value at tens of au. In the second simulation (stratified), we allow the sound speed to vary as a function of the vertical coordinate. To this end, we use the prescription commonly employed when fitting observations of Dartois, Dutrey & Guilloteau (2003) and subsequently updated by Rosenfeld et al. (2013). We assume that the temperature in the upper layer of the disc is 2.5 times the value in the mid-plane (the same as in the isothermal case), and that the transition between the cold mid-plane and the hot upper layer (the parameter zq of Dartois et al. 2003) happens at three scale heights. The initial surface density Σ of the disc follows Σ ∝ r−1. We solve for the hydrostatic equilibrium in the vertical direction to assign the initial density at every point and we initialize the velocity so that the gas is initially in Keplerian rotation. We evolve the simulations for 100 orbits to reach a steady state in the structure of the spirals. 2.2 Radiative transfer simulations We calculate images from the hydrodynamics simulations using the 3D radiative transfer code radmc-3d (http://www.ita.uni-heidelberg.de/dullemond/software/radmc-3d/), using the same 3D spherical mesh of the hydrodynamic simulations. We assume that the central star is a Herbig Ae star (Teff = 9500 K, M⋆ = 2.5 M⊙, R⋆ = 2.0 R⊙), a constant gas-to-dust ratio of 100 and we calculate the dust opacity assuming astronomical silicate (Weingartner & Draine 2001) for a grain size distribution of n(a) ∝ a−3.5 between a = 0.1 μm and a = 1 mm. With radmc-3d, we first calculate the dust temperature with a thermal Monte Carlo simulation and we then calculate images at wavelengths λ = 1.65 μm and λ = 1.3 mm and for an inclination of 0°. Finally, we generate synthetic observations by convolving the calculated images with a Gaussian kernel with an FWHM of 0.04 arcsec, typical of current state-of-the-art NIR cameras such as SPHERE/VLT or GPI/Gemini and of ALMA (ALMA Partnership et al. 2015). We assume a source distance of 100 pc. 3 RESULTS 3.1 Hydrodynamic simulations We show in Fig. 1 the pitch angle of the spiral in the stratified simulation as a function of the height above the mid-plane. We plot with the colours shown in the legend the pitch angle at different radii, normalized to the value in the mid-plane. The red line shows how the sound speed varies with the vertical coordinate. We find that the pitch angle increases with z at all radii; in addition, the dependence of the pitch angle closely resembles the one of the sound speed. Figure 1. View largeDownload slide Dependence of the spiral pitch angle on the latitude in the stratified simulation. The red line shows the sound speed profile used in the hydrodynamic simulations. The large spread in the points close to the upper boundary comes from the difficulty of determining the spiral pitch angle in that region. Figure 1. View largeDownload slide Dependence of the spiral pitch angle on the latitude in the stratified simulation. The red line shows the sound speed profile used in the hydrodynamic simulations. The large spread in the points close to the upper boundary comes from the difficulty of determining the spiral pitch angle in that region. This behaviour is the natural consequence of the temperature dependence of the pitch angle. The spirals are located where constructive interference occurs between waves emitted from Lindblad resonances in the disc (Ogilvie & Lubow 2002); this location depends on the wave propagation speed, and thus the disc temperature. Indeed Rafikov (2002) showed that the pitch angle of a tidally induced spiral in a 2D Keplerian disc, in which the sound speed varies with radius as cs = cp(r/rp)−q, depends on the sound speed:   $$\beta = \pi /2-\tan ^{-1}\left\lbrace \frac{r_{\rm p}}{h_{\rm p}}\left(\frac{r}{r_{\rm p}}\right)^{q+1}\left[\left(\frac{r_{\rm p}}{r}\right)^{3/2} -1\right]\right\rbrace\! ,$$ (1)where β is the pitch angle, rp is the orbital radius of the planet and hp is the pressure scale height (which depends on the sound speed via hp = cp/Ωp). The resemblance between the dependence with the vertical coordinate of the pitch angle and the sound speed profile may imply that this relation holds also in a 3D disc, provided that the local temperature is used in the relation. We leave a detailed analysis of the pitch angle in 3D to further studies and conclude that the vertical temperature gradient increases the pitch angle of the spirals in the upper layers of the disc. 3.2 Simulated observations We show the simulated NIR scattered light and submillimetre continuum images calculated from vertically isothermal (see Figs 2a and b) and thermally stratified (see Figs 2c and d) hydrodynamic simulations. In terms of the general structure, the isothermal and stratified simulations look largely similar, with the only significant difference that the inner part of the disc in the stratified simulation is brighter with respect to the isothermal simulation. This is because the inner edge of the disc in the stratified simulation is more ‘rounded off’, increasing the flaring index of the disc locally and consequently the brightness of the disc. Figure 2. View largeDownload slide Synthetic images in NIR scattered light (Panels a and c) and in submillimetre (Panels b and d) based on 3D hydrodynamic simulations assuming a vertically isothermal (Panels a and b) and stratified (Panels c and d) temperature structure. The scattered light images are scaled with the square of the radial distance to the star. In the NIR, the isothermal simulation shows higher contrast spirals with lower pitch angle compared to the stratified simulation. At submillimetre wavelengths, the pitch angle and the amplitude of the spirals are rather similar in the vertically stratified and isothermal simulations. Figure 2. View largeDownload slide Synthetic images in NIR scattered light (Panels a and c) and in submillimetre (Panels b and d) based on 3D hydrodynamic simulations assuming a vertically isothermal (Panels a and b) and stratified (Panels c and d) temperature structure. The scattered light images are scaled with the square of the radial distance to the star. In the NIR, the isothermal simulation shows higher contrast spirals with lower pitch angle compared to the stratified simulation. At submillimetre wavelengths, the pitch angle and the amplitude of the spirals are rather similar in the vertically stratified and isothermal simulations. 3.2.1 Amplitude of the spirals In Fig. 3, we compare the azimuthal surface brightness profiles of the two models at two radii, 1.7 rp and 0.58 rp. From this plot, we measure the peak-to-peak variation of the intensity in the azimuthal direction and the full width at half-maximum (FWHM) of the spiral peak (for the inner spiral we considered the stronger, primary arm). Since the two models employ the same temperature in the disc mid-plane, in the submillimetre case (Figs 3b and d) the peak-to-peak intensity variation and the width of the spiral peak agree within 10 per cent. In the NIR, scattered light images (Figs 3a and c) instead the spiral contrast is higher in the isothermal simulations by 49 per cent and 45 per cent in the inner and outer arms, respectively. This is caused by the higher temperature in the upper layers in the stratified simulations; the enhanced pressure forces tend to wash out more the spiral wake and reduce its amplitude. While the azimuthal width of the spirals in the two simulations agrees to 4 per cent inwards of the planet, for the outer spirals in the isothermal models the azimuthal width is a factor 3.5 times larger compared to the stratified simulations. Figure 3. View largeDownload slide Comparison of the azimuthal surface brightness profile of the disc in synthetic images calculated from locally isothermal and stratified simulations. The comparison is shown for inner (a, c) and outer spirals (b, d) and both in NIR scattered light (a, b) and in the submillimetre continuum (c, d). The surface brightness profiles are normalized to the mean value in the [200°, 300°] interval. In the NIR, the amplitude and the azimuthal width of the spirals are significantly larger in isothermal simulations compared to the stratified one, while the vertical structure has little effect on the submillimetre continuum observations. Figure 3. View largeDownload slide Comparison of the azimuthal surface brightness profile of the disc in synthetic images calculated from locally isothermal and stratified simulations. The comparison is shown for inner (a, c) and outer spirals (b, d) and both in NIR scattered light (a, b) and in the submillimetre continuum (c, d). The surface brightness profiles are normalized to the mean value in the [200°, 300°] interval. In the NIR, the amplitude and the azimuthal width of the spirals are significantly larger in isothermal simulations compared to the stratified one, while the vertical structure has little effect on the submillimetre continuum observations. 3.2.2 Spiral wake In Fig. 4, we present the spiral wakes we derived from the synthetic images by extracting the positions of the local maxima in the radial surface brightness distribution. In submillimetre images, probing the disc mid-plane, the spirals have similar pitch angles in the locally isothermal simulations and in the vertically stratified simulations (see Figs 4a and d) since the mid-plane temperatures in the two simulations (locally isothermal versus stratified) are identical. Figure 4. View largeDownload slide Panels a and b show the spiral wake in the vertically isothermal simulation inwards and outwards of the planetary orbit, respectively. Panel c shows the difference in the azimuthal position of the spirals between submillimetre and NIR scattered light images as a function of radius. Panels d–f are for the vertically stratified simulation. In a vertically isothermal disc, the shape of the spiral wake is largely the same at both wavelengths. For the stratified simulation, the spirals are more open in the NIR than in the submillimetre both inwards and outwards of the planetary orbit. In the stratified simulation, the difference in the azimuthal position of the spirals between NIR and submillimetre increases with increasing distance from the planet, i.e. the spirals converge towards the planet. Figure 4. View largeDownload slide Panels a and b show the spiral wake in the vertically isothermal simulation inwards and outwards of the planetary orbit, respectively. Panel c shows the difference in the azimuthal position of the spirals between submillimetre and NIR scattered light images as a function of radius. Panels d–f are for the vertically stratified simulation. In a vertically isothermal disc, the shape of the spiral wake is largely the same at both wavelengths. For the stratified simulation, the spirals are more open in the NIR than in the submillimetre both inwards and outwards of the planetary orbit. In the stratified simulation, the difference in the azimuthal position of the spirals between NIR and submillimetre increases with increasing distance from the planet, i.e. the spirals converge towards the planet. Spirals in scattered light images computed from stratified simulations (see Fig. 4e) instead are more open compared to the ones computed from locally isothermal simulations (see Fig. 4b), consistently with the analysis of the simulation in Section 3.1. For example, while we find that outer spirals in NIR scattered light images have a low pitch angle (7$${^{\circ}_{.}}$$5 at 1.5rp) in the isothermal simulations, the pitch angle is higher in the stratified simulation (13° at 1.5rp); the effect is smaller for the inner spiral (15° versus 13° at 0.7rp). Dong et al. (2015a) suggested that a low pitch angle could be used to identify the spiral as an outer one, but it is not clear whether this method still applies to the stratified case. To better illustrate the difference in opening angle, we show the difference in azimuth angles between the spirals derived from scattered light and submillimetre images as a function of stellocentric radius in Figs 4c and f for the locally isothermal simulations and for the stratified simulations, respectively. In the locally isothermal simulation, the difference in the azimuthal position of the spiral wake is not larger than 8° for the inner arms at all radii with a weak trend of increasing difference towards smaller radii. For the outer spirals, the difference fluctuates between −19° and −6° with no clear dependence on the radius. In case of the stratified simulations, there is a clear trend for the azimuthal position difference of the spirals in NIR and submillimetre. The azimuthal difference in the spiral position increases both inwards and outwards from the position of the planet. The difference in the spiral position in NIR and submillimetre images is nearly −20° for the inner arms at 0.55rpl in the stratified simulations compared to the −8° in the locally isothermal simulations. For the outer arms, the azimuthal position difference is about 30° at 2rpl in the stratified simulations while only −15° in the locally isothermal simulations. 4 DISCUSSION AND CONCLUSIONS Direct detection of young planets is notoriously difficult due to the contrast limitation of the observations. Due to their larger spatial extent, spiral density waves excited by young planets are easier to detect. Thus, spiral arms could, in principle, be used as signposts of planets (see e.g. Dong et al. 2015a). However, there is ambiguity in whether the planet is inside or outside the spirals. In the previous section, we have shown that the spiral arms in the submillimetre and in scattered light converge towards the location of the planet due to the different temperature in the mid-plane compared to the upper layers of the disc. Therefore, we suggest to use the direction of convergence of the spirals in NIR and in submillimetre continuum images to infer whether the planet is inwards or outwards of the spirals. Observationally, this requires having high-resolution observations of the same disc in NIR scattered light, probing the upper layers of the disc, and in submillimetre continuum, probing the disc mid-plane. Although not explored in this paper, we note that optically thin gas observations also trace the disc mid-plane and could therefore be used in place of the submillimetre continuum. We also note that, if we know that a planet is responsible for the perturbation, the dependence of the spiral pitch angle on the local temperature could be used as a diagnostic to map the vertical thermal stratification of the disc. This method requires using different tracers (e.g. NIR scattered light, various CO isotopologues, submillimetre continuum) coming from different heights above the mid-plane. Mapping the vertical thermal structure of the disc with the spiral pitch angle has the advantage over other alternative methods (e.g. modelling of the dust continuum emission, molecular rotational transitions), that the pitch angle depends only on the dynamics and is not affected by uncertainties in chemistry and opacities. While the absolute value of the pitch angle might be weakly affected by the planet mass or viscosity (Dong & Fung 2017), the differential variation of it as a function of height is likely independent of these factors. A detailed mapping will require a non-linear theory of spiral density waves in thermally stratified discs, which has yet to be developed. ACKNOWLEDGEMENTS This work has been supported by the DISCSIM project, grant agreement 341137 funded by the European Research Council under ERC-2013-ADG. We thank the referee, Ruobing Dong, for the insightful review that improved our paper. REFERENCES ALMA Partnership et al.  , 2015, ApJ , 808, L3 CrossRef Search ADS   Benisty M. et al.  , 2015, A&A , 578, L6 CrossRef Search ADS   Benisty M. et al.  , 2017, A&A , 597, A42 CrossRef Search ADS   Chiang E. I., Goldreich P., 1997, ApJ , 490, 368 CrossRef Search ADS   Cossins P., Lodato G., Clarke C. J., 2009, MNRAS , 393, 1157 CrossRef Search ADS   Dartois E., Dutrey A., Guilloteau S., 2003, A&A , 399, 773 CrossRef Search ADS   de Val-Borro M. et al.  , 2006, MNRAS , 370, 529 CrossRef Search ADS   Dong R., Fung J., 2017, ApJ , 835, 38 CrossRef Search ADS   Dong R., Zhu Z., Rafikov R. R., Stone J. M., 2015a, ApJ , 809, L5 CrossRef Search ADS   Dong R., Hall C., Rice K., Chiang E., 2015b, ApJ , 812, L32 CrossRef Search ADS   Goldreich P., Tremaine S., 1979, ApJ , 233, 857 CrossRef Search ADS   Grady C. A. et al.  , 2013, ApJ , 762, 48 CrossRef Search ADS   Juhász A., Benisty M., Pohl A., Dullemond C. P., Dominik C., Paardekooper S.-J., 2015, MNRAS , 451, 1147 CrossRef Search ADS   Lee W.-K., Gu P.-G., 2015, ApJ , 814, 72 CrossRef Search ADS   Meru F., Juhász A., Ilee J. D., Clarke C. J., Rosotti G. P., Booth R. A., 2017, ApJ , 839, L24 CrossRef Search ADS   Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS , 170, 228 CrossRef Search ADS   Mignone A., Flock M., Stute M., Kolb S. M., Muscianisi G., 2012, A&A , 545, A152 CrossRef Search ADS   Muto T. et al.  , 2012, ApJ , 748, L22 CrossRef Search ADS   Ogilvie G. I., Lubow S. H., 2002, MNRAS , 330, 950 CrossRef Search ADS   Pérez L. M. et al.  , 2016, Science , 353, 1519 CrossRef Search ADS PubMed  Pohl A., Pinilla P., Benisty M., Ataiee S., Juhász A., Dullemond C. P., Van Boekel R., Henning T., 2015, MNRAS , 453, 1768 CrossRef Search ADS   Rafikov R. R., 2002, ApJ , 569, 997 CrossRef Search ADS   Rice W. K. M., Armitage P. J., Bate M. R., Bonnell I. A., 2003, MNRAS , 338, 227 CrossRef Search ADS   Rosenfeld K. A., Andrews S. M., Hughes A. M., Wilner D. J., Qi C., 2013, ApJ , 774, 16 CrossRef Search ADS   Shakura N. I., Sunyaev R. A., 1973, A&A , 24, 337 Stolker T. et al.  , 2016, A&A , 595, A113 CrossRef Search ADS   Tobin J. J. et al.  , 2016, Nature , 538, 483 CrossRef Search ADS PubMed  Wagner K., Apai D., Kasper M., Robberto M., 2015, ApJ , 813, L2 CrossRef Search ADS   Weingartner J. C., Draine B. T., 2001, ApJ , 548, 296 CrossRef Search ADS   Zhu Z., Dong R., Stone J. M., Rafikov R. R., 2015, ApJ , 813, 88 CrossRef Search ADS   © 2017 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society: Letters Oxford University Press

# Spiral arms in thermally stratified protoplanetary discs

, Volume 474 (1) – Feb 1, 2018
5 pages

/lp/ou_press/spiral-arms-in-thermally-stratified-protoplanetary-discs-foNjXb5kNS
Publisher
journal_eissn:11745-3933
ISSN
1745-3925
eISSN
1745-3933
D.O.I.
10.1093/mnrasl/slx182
Publisher site
See Article on Publisher Site

### Abstract

Abstract Spiral arms have been observed in nearly a dozen protoplanetary discs in near-infrared scattered light and recently also in the submillimetre continuum. While one of the most compelling explanations is that they are driven by planetary or stellar companions, in all but one cases such companions have not yet been detected and there is even ambiguity on whether the planet should be located inside or outside the spirals. Here, we use 3D hydrodynamic simulations to study the morphology of spiral density waves launched by embedded planets taking into account the vertical temperature gradient, a natural consequence of stellar irradiation. Our simulations show that the pitch angle of the spirals in thermally stratified discs is the lowest in the disc mid-plane and increases towards the disc surface. We combine the hydrodynamic simulations with 3D radiative transfer calculations to predict that the pitch angle of planetary spirals observed in the near-infrared is higher than in the submillimetre. We also find that in both cases the spirals converge towards the planet. This provides a new powerful observational method to determine if the perturbing planet is inside or outside the spirals, as well as map the thermal stratification of the disc. accretion, accretion discs, hydrodynamics, protoplanetary discs, circumstellar matter 1 INTRODUCTION With the dawn of extreme high resolution imaging capabilities of current telescopes (e.g. SPHERE/VLT, GPI/Gemini, ALMA) more and more discs show spiral arms in near-infrared (NIR) scattered light (e.g. Muto et al. 2012; Benisty et al. 2015; Wagner et al. 2015; Stolker et al. 2016; Benisty et al. 2017) and recently also in the sub-millimetre (Pérez et al. 2016; Tobin et al. 2016). The observed spirals have similar morphologies in most cases; two symmetric arms, shifted in azimuth by approximately 180°, and with a pitch angle of about 15°–30°. In most cases, there is also a gap/cavity just inwards of the spirals. It has been known for a long time that spiral density waves are launched by planetary or stellar companions at Lindblad resonances (e.g. Goldreich & Tremaine 1979; Ogilvie & Lubow 2002; Rafikov 2002) both inwards and outwards of companion’s orbit. Therefore, a straightforward explanation for the observed spirals is that they are launched by a so far undetected companion, possibly a forming embedded planet (e.g. Muto et al. 2012; Grady et al. 2013); this would open up the exciting possibility of getting a direct insight into planet formation. Using 2D hydrodynamic simulations and 3D radiative transfer computation, Juhász et al. (2015) showed, however, that the pitch angle and contrast of the spirals driven by planets orbiting inwards of the spirals are not compatible with the NIR observations. Dong et al. (2015a) combined 3D hydrodynamic simulation with 3D radiative transfer calculations to show that the observed contrast and pitch angle of spirals are in fact in agreement with the spirals being launched by massive (several MJup) planets exterior to the spirals. However, planets are not the only explanations for the observed spirals. Multi-armed spiral structure can also form due to gravitational instability (GI; see e.g. Rice et al. 2003) with the number of arms depending on the mass and temperature of the disc via the Toomre Q parameter (Cossins, Lodato & Clarke 2009). For very massive discs (Mdisc/M⋆ ≳ 0.25), GI can produce two armed spirals in the disc, similar to the NIR scattered light observations (Dong et al. 2015b). Meru et al. (2017) showed that GI could explain the spirals in the submillimetre ALMA images of Elias 2–27. Additionally, if a marginally gravitationally stable disc is perturbed by an embedded planet, the tidal interaction might tip the balance and trigger GI in the disc (Pohl et al. 2015). This combined mechanism may also produce spirals with pitch angles similar to the ones in observations. Numerous theoretical studies have recently increased the level of complexity in the modelling of spirals, showing the importance to understand their origin. For instance, while tidally induced spiral density waves have been studied for a long time in 2D, Zhu et al. (2015) demonstrated the spiral amplitude and pitch angle (the complement of the angle between the radial direction and the tangent of the spiral) in the non-linear regime are higher in 3D than in 2D simulations. However, an important aspect of a realistic disc model has not yet been investigated. Protoplanetary discs in which spirals have been detected are passively heated by the absorbed stellar radiation. As a consequence of this, they will have a positive vertical temperature gradient (Chiang & Goldreich 1997), resulting in a vertical thermal stratification in the disc. It is known that the pitch angle of the spirals depends on the local sound speed (Rafikov 2002), thus we expect vertical thermal stratification to affect the morphology of the spiral wake. (see e.g. Lee & Gu 2015). The goal of this study is to investigate the morphology of tidally induced spirals by embedded planets in 3D thermally stratified discs. We note, however, that our results would also be applicable to stellar companions as the tidal interaction between the disc and the companion is the same as for the planetary case. 2 MODELLING 2.1 Hydrodynamic simulations We run 3D hydrodynamical simulations with the code pluto (Mignone et al. 2007) to investigate the disc response to the planet. We include the FARGO algorithm (Mignone et al. 2012) to reduce the computational cost and use the HLLC Riemann solver, linear reconstruction and second-order Runge–Kutta time integration. We employ a spherical grid with 400, 50, and 1024 cells, distributed in the [0.4, 2.5], [1.22, π/2], and [0, 2π], intervals in the radial, polar, and azimuthal directions. The planet was placed at the coordinates [1.0, 0, π/4]. In the radial and polar directions, we define ‘wave-killing zones’ where we damp disturbances, following de Val-Borro et al. (2006), extending up to 0.44 and from 2.25 in radius and until 1.27 in θ. We employ a locally isothermal equation of state, in which the sound speed cs is a function of position only, and a physical viscosity using the Shakura & Sunyaev (1973) prescription with a value of α = 10−3. We assume the planet mass to be 3 MJup. We run two simulations, the first (subsequently called ‘isothermal’) in which cs depends only on radius: cs ∝ r−1/4, with the normalization set such that the disc aspect ratio at the location of the planet is 0.07, a typical value at tens of au. In the second simulation (stratified), we allow the sound speed to vary as a function of the vertical coordinate. To this end, we use the prescription commonly employed when fitting observations of Dartois, Dutrey & Guilloteau (2003) and subsequently updated by Rosenfeld et al. (2013). We assume that the temperature in the upper layer of the disc is 2.5 times the value in the mid-plane (the same as in the isothermal case), and that the transition between the cold mid-plane and the hot upper layer (the parameter zq of Dartois et al. 2003) happens at three scale heights. The initial surface density Σ of the disc follows Σ ∝ r−1. We solve for the hydrostatic equilibrium in the vertical direction to assign the initial density at every point and we initialize the velocity so that the gas is initially in Keplerian rotation. We evolve the simulations for 100 orbits to reach a steady state in the structure of the spirals. 2.2 Radiative transfer simulations We calculate images from the hydrodynamics simulations using the 3D radiative transfer code radmc-3d (http://www.ita.uni-heidelberg.de/dullemond/software/radmc-3d/), using the same 3D spherical mesh of the hydrodynamic simulations. We assume that the central star is a Herbig Ae star (Teff = 9500 K, M⋆ = 2.5 M⊙, R⋆ = 2.0 R⊙), a constant gas-to-dust ratio of 100 and we calculate the dust opacity assuming astronomical silicate (Weingartner & Draine 2001) for a grain size distribution of n(a) ∝ a−3.5 between a = 0.1 μm and a = 1 mm. With radmc-3d, we first calculate the dust temperature with a thermal Monte Carlo simulation and we then calculate images at wavelengths λ = 1.65 μm and λ = 1.3 mm and for an inclination of 0°. Finally, we generate synthetic observations by convolving the calculated images with a Gaussian kernel with an FWHM of 0.04 arcsec, typical of current state-of-the-art NIR cameras such as SPHERE/VLT or GPI/Gemini and of ALMA (ALMA Partnership et al. 2015). We assume a source distance of 100 pc. 3 RESULTS 3.1 Hydrodynamic simulations We show in Fig. 1 the pitch angle of the spiral in the stratified simulation as a function of the height above the mid-plane. We plot with the colours shown in the legend the pitch angle at different radii, normalized to the value in the mid-plane. The red line shows how the sound speed varies with the vertical coordinate. We find that the pitch angle increases with z at all radii; in addition, the dependence of the pitch angle closely resembles the one of the sound speed. Figure 1. View largeDownload slide Dependence of the spiral pitch angle on the latitude in the stratified simulation. The red line shows the sound speed profile used in the hydrodynamic simulations. The large spread in the points close to the upper boundary comes from the difficulty of determining the spiral pitch angle in that region. Figure 1. View largeDownload slide Dependence of the spiral pitch angle on the latitude in the stratified simulation. The red line shows the sound speed profile used in the hydrodynamic simulations. The large spread in the points close to the upper boundary comes from the difficulty of determining the spiral pitch angle in that region. This behaviour is the natural consequence of the temperature dependence of the pitch angle. The spirals are located where constructive interference occurs between waves emitted from Lindblad resonances in the disc (Ogilvie & Lubow 2002); this location depends on the wave propagation speed, and thus the disc temperature. Indeed Rafikov (2002) showed that the pitch angle of a tidally induced spiral in a 2D Keplerian disc, in which the sound speed varies with radius as cs = cp(r/rp)−q, depends on the sound speed:   $$\beta = \pi /2-\tan ^{-1}\left\lbrace \frac{r_{\rm p}}{h_{\rm p}}\left(\frac{r}{r_{\rm p}}\right)^{q+1}\left[\left(\frac{r_{\rm p}}{r}\right)^{3/2} -1\right]\right\rbrace\! ,$$ (1)where β is the pitch angle, rp is the orbital radius of the planet and hp is the pressure scale height (which depends on the sound speed via hp = cp/Ωp). The resemblance between the dependence with the vertical coordinate of the pitch angle and the sound speed profile may imply that this relation holds also in a 3D disc, provided that the local temperature is used in the relation. We leave a detailed analysis of the pitch angle in 3D to further studies and conclude that the vertical temperature gradient increases the pitch angle of the spirals in the upper layers of the disc. 3.2 Simulated observations We show the simulated NIR scattered light and submillimetre continuum images calculated from vertically isothermal (see Figs 2a and b) and thermally stratified (see Figs 2c and d) hydrodynamic simulations. In terms of the general structure, the isothermal and stratified simulations look largely similar, with the only significant difference that the inner part of the disc in the stratified simulation is brighter with respect to the isothermal simulation. This is because the inner edge of the disc in the stratified simulation is more ‘rounded off’, increasing the flaring index of the disc locally and consequently the brightness of the disc. Figure 2. View largeDownload slide Synthetic images in NIR scattered light (Panels a and c) and in submillimetre (Panels b and d) based on 3D hydrodynamic simulations assuming a vertically isothermal (Panels a and b) and stratified (Panels c and d) temperature structure. The scattered light images are scaled with the square of the radial distance to the star. In the NIR, the isothermal simulation shows higher contrast spirals with lower pitch angle compared to the stratified simulation. At submillimetre wavelengths, the pitch angle and the amplitude of the spirals are rather similar in the vertically stratified and isothermal simulations. Figure 2. View largeDownload slide Synthetic images in NIR scattered light (Panels a and c) and in submillimetre (Panels b and d) based on 3D hydrodynamic simulations assuming a vertically isothermal (Panels a and b) and stratified (Panels c and d) temperature structure. The scattered light images are scaled with the square of the radial distance to the star. In the NIR, the isothermal simulation shows higher contrast spirals with lower pitch angle compared to the stratified simulation. At submillimetre wavelengths, the pitch angle and the amplitude of the spirals are rather similar in the vertically stratified and isothermal simulations. 3.2.1 Amplitude of the spirals In Fig. 3, we compare the azimuthal surface brightness profiles of the two models at two radii, 1.7 rp and 0.58 rp. From this plot, we measure the peak-to-peak variation of the intensity in the azimuthal direction and the full width at half-maximum (FWHM) of the spiral peak (for the inner spiral we considered the stronger, primary arm). Since the two models employ the same temperature in the disc mid-plane, in the submillimetre case (Figs 3b and d) the peak-to-peak intensity variation and the width of the spiral peak agree within 10 per cent. In the NIR, scattered light images (Figs 3a and c) instead the spiral contrast is higher in the isothermal simulations by 49 per cent and 45 per cent in the inner and outer arms, respectively. This is caused by the higher temperature in the upper layers in the stratified simulations; the enhanced pressure forces tend to wash out more the spiral wake and reduce its amplitude. While the azimuthal width of the spirals in the two simulations agrees to 4 per cent inwards of the planet, for the outer spirals in the isothermal models the azimuthal width is a factor 3.5 times larger compared to the stratified simulations. Figure 3. View largeDownload slide Comparison of the azimuthal surface brightness profile of the disc in synthetic images calculated from locally isothermal and stratified simulations. The comparison is shown for inner (a, c) and outer spirals (b, d) and both in NIR scattered light (a, b) and in the submillimetre continuum (c, d). The surface brightness profiles are normalized to the mean value in the [200°, 300°] interval. In the NIR, the amplitude and the azimuthal width of the spirals are significantly larger in isothermal simulations compared to the stratified one, while the vertical structure has little effect on the submillimetre continuum observations. Figure 3. View largeDownload slide Comparison of the azimuthal surface brightness profile of the disc in synthetic images calculated from locally isothermal and stratified simulations. The comparison is shown for inner (a, c) and outer spirals (b, d) and both in NIR scattered light (a, b) and in the submillimetre continuum (c, d). The surface brightness profiles are normalized to the mean value in the [200°, 300°] interval. In the NIR, the amplitude and the azimuthal width of the spirals are significantly larger in isothermal simulations compared to the stratified one, while the vertical structure has little effect on the submillimetre continuum observations. 3.2.2 Spiral wake In Fig. 4, we present the spiral wakes we derived from the synthetic images by extracting the positions of the local maxima in the radial surface brightness distribution. In submillimetre images, probing the disc mid-plane, the spirals have similar pitch angles in the locally isothermal simulations and in the vertically stratified simulations (see Figs 4a and d) since the mid-plane temperatures in the two simulations (locally isothermal versus stratified) are identical. Figure 4. View largeDownload slide Panels a and b show the spiral wake in the vertically isothermal simulation inwards and outwards of the planetary orbit, respectively. Panel c shows the difference in the azimuthal position of the spirals between submillimetre and NIR scattered light images as a function of radius. Panels d–f are for the vertically stratified simulation. In a vertically isothermal disc, the shape of the spiral wake is largely the same at both wavelengths. For the stratified simulation, the spirals are more open in the NIR than in the submillimetre both inwards and outwards of the planetary orbit. In the stratified simulation, the difference in the azimuthal position of the spirals between NIR and submillimetre increases with increasing distance from the planet, i.e. the spirals converge towards the planet. Figure 4. View largeDownload slide Panels a and b show the spiral wake in the vertically isothermal simulation inwards and outwards of the planetary orbit, respectively. Panel c shows the difference in the azimuthal position of the spirals between submillimetre and NIR scattered light images as a function of radius. Panels d–f are for the vertically stratified simulation. In a vertically isothermal disc, the shape of the spiral wake is largely the same at both wavelengths. For the stratified simulation, the spirals are more open in the NIR than in the submillimetre both inwards and outwards of the planetary orbit. In the stratified simulation, the difference in the azimuthal position of the spirals between NIR and submillimetre increases with increasing distance from the planet, i.e. the spirals converge towards the planet. Spirals in scattered light images computed from stratified simulations (see Fig. 4e) instead are more open compared to the ones computed from locally isothermal simulations (see Fig. 4b), consistently with the analysis of the simulation in Section 3.1. For example, while we find that outer spirals in NIR scattered light images have a low pitch angle (7$${^{\circ}_{.}}$$5 at 1.5rp) in the isothermal simulations, the pitch angle is higher in the stratified simulation (13° at 1.5rp); the effect is smaller for the inner spiral (15° versus 13° at 0.7rp). Dong et al. (2015a) suggested that a low pitch angle could be used to identify the spiral as an outer one, but it is not clear whether this method still applies to the stratified case. To better illustrate the difference in opening angle, we show the difference in azimuth angles between the spirals derived from scattered light and submillimetre images as a function of stellocentric radius in Figs 4c and f for the locally isothermal simulations and for the stratified simulations, respectively. In the locally isothermal simulation, the difference in the azimuthal position of the spiral wake is not larger than 8° for the inner arms at all radii with a weak trend of increasing difference towards smaller radii. For the outer spirals, the difference fluctuates between −19° and −6° with no clear dependence on the radius. In case of the stratified simulations, there is a clear trend for the azimuthal position difference of the spirals in NIR and submillimetre. The azimuthal difference in the spiral position increases both inwards and outwards from the position of the planet. The difference in the spiral position in NIR and submillimetre images is nearly −20° for the inner arms at 0.55rpl in the stratified simulations compared to the −8° in the locally isothermal simulations. For the outer arms, the azimuthal position difference is about 30° at 2rpl in the stratified simulations while only −15° in the locally isothermal simulations. 4 DISCUSSION AND CONCLUSIONS Direct detection of young planets is notoriously difficult due to the contrast limitation of the observations. Due to their larger spatial extent, spiral density waves excited by young planets are easier to detect. Thus, spiral arms could, in principle, be used as signposts of planets (see e.g. Dong et al. 2015a). However, there is ambiguity in whether the planet is inside or outside the spirals. In the previous section, we have shown that the spiral arms in the submillimetre and in scattered light converge towards the location of the planet due to the different temperature in the mid-plane compared to the upper layers of the disc. Therefore, we suggest to use the direction of convergence of the spirals in NIR and in submillimetre continuum images to infer whether the planet is inwards or outwards of the spirals. Observationally, this requires having high-resolution observations of the same disc in NIR scattered light, probing the upper layers of the disc, and in submillimetre continuum, probing the disc mid-plane. Although not explored in this paper, we note that optically thin gas observations also trace the disc mid-plane and could therefore be used in place of the submillimetre continuum. We also note that, if we know that a planet is responsible for the perturbation, the dependence of the spiral pitch angle on the local temperature could be used as a diagnostic to map the vertical thermal stratification of the disc. This method requires using different tracers (e.g. NIR scattered light, various CO isotopologues, submillimetre continuum) coming from different heights above the mid-plane. Mapping the vertical thermal structure of the disc with the spiral pitch angle has the advantage over other alternative methods (e.g. modelling of the dust continuum emission, molecular rotational transitions), that the pitch angle depends only on the dynamics and is not affected by uncertainties in chemistry and opacities. While the absolute value of the pitch angle might be weakly affected by the planet mass or viscosity (Dong & Fung 2017), the differential variation of it as a function of height is likely independent of these factors. A detailed mapping will require a non-linear theory of spiral density waves in thermally stratified discs, which has yet to be developed. ACKNOWLEDGEMENTS This work has been supported by the DISCSIM project, grant agreement 341137 funded by the European Research Council under ERC-2013-ADG. We thank the referee, Ruobing Dong, for the insightful review that improved our paper. REFERENCES ALMA Partnership et al.  , 2015, ApJ , 808, L3 CrossRef Search ADS   Benisty M. et al.  , 2015, A&A , 578, L6 CrossRef Search ADS   Benisty M. et al.  , 2017, A&A , 597, A42 CrossRef Search ADS   Chiang E. I., Goldreich P., 1997, ApJ , 490, 368 CrossRef Search ADS   Cossins P., Lodato G., Clarke C. J., 2009, MNRAS , 393, 1157 CrossRef Search ADS   Dartois E., Dutrey A., Guilloteau S., 2003, A&A , 399, 773 CrossRef Search ADS   de Val-Borro M. et al.  , 2006, MNRAS , 370, 529 CrossRef Search ADS   Dong R., Fung J., 2017, ApJ , 835, 38 CrossRef Search ADS   Dong R., Zhu Z., Rafikov R. R., Stone J. M., 2015a, ApJ , 809, L5 CrossRef Search ADS   Dong R., Hall C., Rice K., Chiang E., 2015b, ApJ , 812, L32 CrossRef Search ADS   Goldreich P., Tremaine S., 1979, ApJ , 233, 857 CrossRef Search ADS   Grady C. A. et al.  , 2013, ApJ , 762, 48 CrossRef Search ADS   Juhász A., Benisty M., Pohl A., Dullemond C. P., Dominik C., Paardekooper S.-J., 2015, MNRAS , 451, 1147 CrossRef Search ADS   Lee W.-K., Gu P.-G., 2015, ApJ , 814, 72 CrossRef Search ADS   Meru F., Juhász A., Ilee J. D., Clarke C. J., Rosotti G. P., Booth R. A., 2017, ApJ , 839, L24 CrossRef Search ADS   Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS , 170, 228 CrossRef Search ADS   Mignone A., Flock M., Stute M., Kolb S. M., Muscianisi G., 2012, A&A , 545, A152 CrossRef Search ADS   Muto T. et al.  , 2012, ApJ , 748, L22 CrossRef Search ADS   Ogilvie G. I., Lubow S. H., 2002, MNRAS , 330, 950 CrossRef Search ADS   Pérez L. M. et al.  , 2016, Science , 353, 1519 CrossRef Search ADS PubMed  Pohl A., Pinilla P., Benisty M., Ataiee S., Juhász A., Dullemond C. P., Van Boekel R., Henning T., 2015, MNRAS , 453, 1768 CrossRef Search ADS   Rafikov R. R., 2002, ApJ , 569, 997 CrossRef Search ADS   Rice W. K. M., Armitage P. J., Bate M. R., Bonnell I. A., 2003, MNRAS , 338, 227 CrossRef Search ADS   Rosenfeld K. A., Andrews S. M., Hughes A. M., Wilner D. J., Qi C., 2013, ApJ , 774, 16 CrossRef Search ADS   Shakura N. I., Sunyaev R. A., 1973, A&A , 24, 337 Stolker T. et al.  , 2016, A&A , 595, A113 CrossRef Search ADS   Tobin J. J. et al.  , 2016, Nature , 538, 483 CrossRef Search ADS PubMed  Wagner K., Apai D., Kasper M., Robberto M., 2015, ApJ , 813, L2 CrossRef Search ADS   Weingartner J. C., Draine B. T., 2001, ApJ , 548, 296 CrossRef Search ADS   Zhu Z., Dong R., Stone J. M., Rafikov R. R., 2015, ApJ , 813, 88 CrossRef Search ADS   © 2017 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical Society: LettersOxford University Press

Published: Feb 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations