INFERRING PLANET MASS FROM SPIRAL STRUCTURES IN PROTOPLANETARY DISKS

Jeffrey Fung (馮澤之); Ruobing Dong (董若冰)

doi:10.1088/2041-8205/815/2/L21

1. INTRODUCTION

Young, forming planets embedded in protoplanetary disks are expected to generate global spiral structures that are potentially observable. Recent observations show a variety of spiral structures in the scattered light images of protoplanetary disks, such as those in HD 100453 (Wagner et al. 2015), MWC 758 (Grady et al. 2013; Benisty et al. 2015), SAO 206462 (Muto et al. 2012; Garufi et al. 2013), HD 142527 (Casassus et al. 2012; Rameau et al. 2012; Canovas et al. 2013; Avenhaus et al. 2014b; Christiaens et al. 2014), and HD 100546 (Boccaletti et al. 2013; Avenhaus et al. 2014a; Currie et al. 2014). In all cases, the observed disk contains more than a single spiral arm, and while this may point to the existence of multiple planets, it has been shown that a massive planet alone can excite multi-armed spiral structures; e.g., Dong et al. (2015a) found that the two-armed spirals in both MWC 758 and SAO 206462 are consistent with the influence of a single massive planet.

In the foundational theory of disk–planet interaction (Goldreich & Tremaine 1979; Artymowicz 1993), planets excite density waves in disks over a range of $m{\rm{th}}$ order azimuthal modes, and these waves interfere to create an one-armed spiral structure (Ogilvie & Lubow 2002). This picture is based on the assumption that the planet's perturbation on the disk is sufficiently small, so that the perturbation equations can be linearized. Typically this refers to planets whose "thermal mass", defined as ${q}_{{\rm{th}}}\equiv q/{H}_{0}^{3},$ where q is the planet-to-star mass ratio, and H₀ the disk aspect ratio at the planet's location, is much less than unity. If we consider a more massive planet with ${q}_{{\rm{th}}}\sim 1,$ then it can also excite "ultraharmonic" density waves, generating a secondary, or even tertiary, set of waves that can interfere with the primary set to create multi-armed structures.

Ultraharmonic wave excitation can be understood as the result of high order coupling between perturbations of different $m{\rm{th}}$ azimuthal modes, and is a natural way to extend the conventional perturbation theory into the weakly nonlinear regime. Originally developed under the context of galactic dynamics, Artymowicz & Lubow (1992) showed that the quadratic perturbation terms, or the self-coupling terms, dropped in the linear analysis by Goldreich & Tremaine (1979) are in fact the driving terms for a second order perturbation that excites density waves at the 2:1 resonance with the primary arm. This can be generalized to even higher order perturbations, such as a third order perturbation at the 3:1 resonance, driven by the coupling between the first and second order terms. Evidence for these features may have already been found in previous work, such as Zhu et al. (2015), who reported and characterized a secondary and tertiary arm in their simulations of disk–planet interaction; and Juhász et al. (2015), who also reported the existence of a secondary arm when ${q}_{{\rm{th}}}\gtrsim 1.$

These spiral features are potential diagnostics for inferring disk and planet properties; for instance, Zhu et al. (2015) proposed the use of pitch angle, spiral separation, and spiral contrast to infer planet mass. Despite that, the properties of multi-armed structures remain poorly understood, due to the complexity of its nonlinear nature. In this paper, we investigate how the shape and location of the secondary arm vary with planet mass and disk temperature, with an emphasis on its applications to protoplanetary disk observations. Throughout this study, we focus on the part of the disk interior to the planet, because the outer disk is usually much fainter than the inner disk and less likely to be observed. We defer the thorough investigation of ultraharmonic resonances in disk–planet interaction to a future paper.

2. NUMERICAL METHOD

As the basis for our analysis, we perform global 2D disk–planet interaction simulations using the Graphics Processing Unit-accelerated hydrodynamics code PEnGUIn (Fung 2015). PEnGUIn solves the Navier–Stokes equations for a compressible fluid in the Lagrangian frame:

$\begin{eqnarray}&&\displaystyle \frac{D{\rm{\Sigma }}}{{Dt}}=-{\rm{\Sigma }}\left({\rm{\nabla }}\cdot {\boldsymbol{v}}\right),\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray}&&\displaystyle \frac{D{\boldsymbol{v}}}{{Dt}}=-\displaystyle \frac{1}{{\rm{\Sigma }}}{\rm{\nabla }}P+\displaystyle \frac{1}{{\rm{\Sigma }}}{\rm{\nabla }}\cdot {\mathbb{T}}-{\rm{\nabla }}{\rm{\Phi }},\end{eqnarray} \tag{ 2 }$

where Σ is the disk surface density, ${\boldsymbol{v}}$ the velocity field, P the vertically averaged gas pressure, ${\mathbb{T}}$ the Newtonian stress tensor, and Φ the combined gravitational potential of the star and the planet. We adopt a locally isothermal equation of state, such that $P={c}^{2}{\rm{\Sigma }},$ where c is the sound speed of the gas. ${\mathbb{T}}$ is proportional to the kinematic viscosity ν, which we parameterize using the α-prescription by Shakura & Sunyaev (1973), such that $\nu =\alpha {c}^{2}/{\rm{\Omega }},$ where Ω is the orbital frequency of the disk. In all our simulations, we choose $\alpha =0.001,$ and verify by experimenting with a 10 times higher α that the specific choice of α has negligible effects on the spiral structure of the disks.

Technically, the simulation is run in the corotating frame of the planet, but the Coriolis force is absorbed into the conservative form of Equation (2) as suggested by Kley (1998). We solve the equations in the rest frame of the star, such that the star is always stationary at r = 0. This allows Φ to be written as

$\begin{eqnarray}{\rm{\Phi }} & = & -{GM}\left[\displaystyle \frac{1-q}{r}\right.\\ & & \left.-\displaystyle \frac{q}{\sqrt{{r}^{2}+{r}_{{\rm{s}}}^{2}+{r}_{{\rm{p}}}^{2}-2{{rr}}_{{\rm{p}}}\mathrm{cos}{\phi }^{\prime }}}-\displaystyle \frac{{qr}\mathrm{cos}{\phi }^{\prime }}{{r}_{{\rm{p}}}^{2}}\right],\end{eqnarray} \tag{ 3 }$

where G is the gravitational constant, M the total mass of the star and the planet, ${r}_{{\rm{p}}}$ the semimajor axis of the planet's orbit, ${r}_{{\rm{s}}}$ the softening length of the planet's potential, and ${\phi }^{\prime }=\phi -{\phi }_{{\rm{p}}}$ denotes the azimuthal separation from the planet. The third term in the in bracket of Equation (3) is the indirect potential introduced due to the acceleration of the frame. ${r}_{{\rm{s}}}$ is chosen to be $0.5{H}_{0}{r}_{{\rm{p}}},$ appropriate for approximating the vertically averaged gravitational force (Müller et al. 2012). We fix the planet to a circular orbit (i.e., planetary migration is ignored), so it is always positioned at $({r}_{{\rm{p}}},{\phi }_{{\rm{p}}}).$ We set ${GM}=1,$ and ${r}_{{\rm{p}}}=1,$ such that the planet's orbital frequency ${{\rm{\Omega }}}_{{\rm{p}}}$ is 1, and orbital period ${P}_{{\rm{p}}}$ is $2\pi .$ For convenience, we also denote the Keplerian orbital speed ${v}_{{\rm{k}}}$ and frequency ${{\rm{\Omega }}}_{{\rm{k}}}.$

2.1. Disk Model

Our disk model uses a power law profile for the initial surface density:

$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0}={{\rm{\Sigma }}}_{{\rm{p}}}{\left(\displaystyle \frac{r}{{r}_{{\rm{p}}}}\right)}^{-1},\end{eqnarray} \tag{ 4 }$

where ${{\rm{\Sigma }}}_{{\rm{p}}}$ is set to 1.³ Our choice of a fixed sound speed profile is

$\begin{eqnarray}&&c={c}_{0}{\left(\displaystyle \frac{r}{{r}_{{\rm{p}}}}\right)}^{-1/4}.\end{eqnarray} \tag{ 5 }$

This sound speed profile creates a disk with a flaring shape, since the disk aspect ratio $H=c/{v}_{{\rm{k}}}$ is proportional to ${r}^{1/4}.$ The normalization c₀ is determined in terms of the aspect ratio at the planet's location ${H}_{0}={c}_{0}/{v}_{{\rm{k}}}({r}_{{\rm{p}}}).$ Finally, we complete our initial conditions with the velocity field ${\boldsymbol{v}}.$ We set the initial radial velocity to zero everywhere, and the angular velocity to

$\begin{eqnarray}&&{\rm{\Omega }}=\sqrt{{{\rm{\Omega }}}_{{\rm{k}}}^{2}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial P}{\partial r}}.\end{eqnarray} \tag{ 6 }$

To obtain a thorough picture of the planet-excited spiral arms, we survey over two parameters: q and H₀. We vary q from $6.25\times {10}^{-5},$ just under Neptune mass, to $6.4\times {10}^{-2},$ around brown dwarf mass; and H₀ from 0.04 to 0.16. For reference, we note that Jupiter has a q of ${10}^{-3},$ and a change in H₀ can be interpreted as a change in the planet's location— ${H}_{0}=0.1,$ for instance, is equivalent to the planet being placed at about 25 AU away from a 0.5 solar mass ( ${M}_{\odot }$ ) star in a passively heated disk (Chiang & Goldreich 1997; also see Table 4 of Andrews et al. 2011 for disk scale heights in observed systems).

2.2. Setup

Our computational grid spans $r=\{0.125{r}_{{\rm{p}}},2{r}_{{\rm{p}}}\},$ and covers the full $2\pi$ range in azimuth. The azimuthal boundary conditions are periodic, and the radial boundaries contain fixed values as given by Equations (4)–(6). For cases when ${H}_{0}\geqslant 0.1,$ we use a resolution of $400(r)\times 900(\phi ),$ with the cells distributed logarithmically in the radial direction, and uniformly in azimuth. This gives a cell size of ∼0.007r in both directions. For models with ${H}_{0}\lt 0.1,$ we increase the resolution to $600(r)\times 1400(\phi ),$ so that even for the smallest H₀ used, 0.04, we still retain a resolution of more than eight cells per scale height near the planet.

The spiral structure is typically fully established after about one sound-crossing time, which, for our models, is about 1–4 ${P}_{{\rm{p}}}.$ We therefore run all our simulations for a minimum of ${t}_{{\rm{end}}}=10{P}_{{\rm{p}}},$ and we verify that in most cases, the spiral structure is steady after 5 ${P}_{{\rm{p}}}$ or less. Since the more massive planets open a gap over a viscous timescale much longer than ${P}_{{\rm{p}}},$ we additionally run our q = 0.001 and H₀ = 0.1 model to 300 ${P}_{{\rm{p}}}$ until the gap has fully opened, and confirm this does not alter the morphology of the spirals significantly.⁴

For models with $q\gt 0.001,$ where the strong influence from the planet generates fluctuations that require a longer time to damp out. In these cases, we run the simulations for 20 ${P}_{{\rm{p}}}$ when $q=0.004,$ and extend it to 30 ${P}_{{\rm{p}}}$ when $q\geqslant 0.016.$ These very massive companions are also able to carve out a deep gap within the timescale of our simulations, so the reported spiral morphologies of these models do include the effects of gap opening.

We additionally perform a three-dimensional (3D) simulation to compare with our 2D results. It uses a similar setup as our 2D simulations, and the differences are detailed in Section 3.2.

3. RESULTS

We perform a total of 20 simulations with different sets of $\{q,\;{H}_{0}\}$ parameters. The results are recorded in Table 1, and Figure 1 shows a montage of snapshots illustrating the surface density profile of nine different models. In these maps, we plot the normalized fractional variation in the surface density, Δ, in reference to the initial profile (Equation (4)):

$\begin{eqnarray}&&{\rm{\Delta }}=\displaystyle \frac{1}{A}\displaystyle \frac{{\rm{\Sigma }}-{{\rm{\Sigma }}}_{0}}{{{\rm{\Sigma }}}_{0}},\end{eqnarray} \tag{ 7 }$

where the normalization A is

$\begin{eqnarray}&&A={\left(1+3\;{q}_{{\rm{th}}}\right)}^{\frac{1}{3}}-1.\end{eqnarray} \tag{ 8 }$

This choice of A is motivated by the notion that in the linear limit of ${q}_{{\rm{th}}}\ll 1,$ the amplitude of the density perturbation is well described by $A\approx {q}_{{\rm{th}}},$ while in the nonlinear, ${q}_{{\rm{th}}}\gt 1,$ regime, we speculate that the amplitude scales as ${r}_{{\rm{H}}}/({H}_{0}{r}_{{\rm{p}}})\propto {q}_{{\rm{th}}}^{1/3},$ the ratio between the planet's Hill radius ${r}_{{\rm{H}}}$ and the disk scale height at the planet's orbit. Our prescription for A connects these two regimes to provide a smooth scaling across our sample. As can be seen in Figure 1, this amplitude scaling matches our results well.

**Figure 1.** Snapshots at the end point of nine selected simulations. The color shows the normalized fractional variations in surface density Δ, as defined by Equation (7). H₀ increases from right to left, where the left column contains snapshots of models with ${H}_{0}=0.16,$ middle column has ${H}_{0}=0.1,$ and right column has ${H}_{0}=0.063$ . q increases from top to bottom, where the top row has $q=0.001,$ middle row has $q=0.004,$ and bottom row has q = 0.016. The green dashed circles have a radius of $0.5\;{r}_{{\rm{p}}},$ and are provided as a length scale. All snapshots has the same orientation where the planet is located to the right of the picture, outside of the frame; see Figure 2 for a zoomed-out version that includes the planet.
Download figure:
Standard image High-resolution image

**Figure 1.** Snapshots at the end point of nine selected simulations. The color shows the normalized fractional variations in surface density Δ, as defined by Equation (7). H₀ increases from right to left, where the left column contains snapshots of models with ${H}_{0}=0.16,$ middle column has ${H}_{0}=0.1,$ and right column has ${H}_{0}=0.063$ . q increases from top to bottom, where the top row has $q=0.001,$ middle row has $q=0.004,$ and bottom row has q = 0.016. The green dashed circles have a radius of $0.5\;{r}_{{\rm{p}}},$ and are provided as a length scale. All snapshots has the same orientation where the planet is located to the right of the picture, outside of the frame; see Figure 2 for a zoomed-out version that includes the planet.
Download figure:
Standard image High-resolution image

Table 1. Simulated Models

H₀	q	${\phi }_{{\rm{sep}}}$ (degree)	${q}_{{\rm{th}}}$	${t}_{{\rm{end}}}\;({P}_{{\rm{p}}})$	{r₁, r₂} (a)
0.16	$6.25\times {10}^{-5}$	⋯^a	1/64	10	⋯
0.16	$2.5\times {10}^{-4}$	⋯^a	1/16	10	⋯
0.16	$1.0\times {10}^{-3}$	110(+6)(−5)	1/4	10	{0.4, 0.6}
0.16	$4.0\times {10}^{-3}$	125(+3)(−4)	1	20	{0.4, 0.6}
0.16	$1.6\times {10}^{-2}$	179(+22)(−19)	4	30	{0.4, 0.55}
0.16	$6.4\times {10}^{-2}$	178(+6)(−7)	16	30	{0.4, 0.55}
0.1	$6.25\times {10}^{-5}$	⋯^a	1/16	10	⋯
0.1	$2.5\times {10}^{-4}$	81(+8)(−6)	1/4	10	{0.4, 0.7}
0.1	$1.0\times {10}^{-3}$	99(+10)(−9)	1	10	{0.4, 0.7}
0.1	$4.0\times {10}^{-3}$	125(+21)(−18)	4	20	{0.4, 0.65}
0.1	$1.6\times {10}^{-2}$	172(+6)(−11)	16	30	{0.4, 0.55}
0.1	$6.4\times {10}^{-2}$	195(+4)(−2)	64	30	{0.4, 0.5}
0.063	$6.25\times {10}^{-5}$	56(+6)(−7)	1/4	10	{0.4, 0.7}
0.063	$2.5\times {10}^{-4}$	75(+8)(−15)	1	10	{0.4, 0.7}
0.063	$1.0\times {10}^{-3}$	102(+11)(−15)	4	10	{0.4, 0.7}
0.063	$4.0\times {10}^{-3}$	143(+19)(−57)	16	20	{0.4, 0.6}
0.063	$1.6\times {10}^{-2}$	183(+20)(−22)	64	30	{0.4, 0.5}
0.04	$6.25\times {10}^{-5}$	60(+23)(−12)	1	10	{0.4, 0.7}
0.04	$1.0\times {10}^{-3}$	100(+19)(−22)	16	10	{0.4, 0.7}
0.04	$1.6\times {10}^{-2}$	173(+13)(−30)	256	30	{0.4, 0.6}

Note.

^aThe secondary arms in these models are too weak to be precisely located.

Download table as: ASCII Typeset image

Ultraharmonic theory predicts that the second order perturbation is excited at the 2:1 resonance with respect to the first order perturbation, or equivalently, the planet, since the primary arm corotates with the planet. Consequently, if the secondary arm is a product of the second order perturbation, we expect to find it being launched at $r/{r}_{{\rm{p}}}\sim 0.63.$ One can even predict a tertiary arm being launched at the 3:1 resonance, or $r/{r}_{{\rm{p}}}=0.48.$ On the left panel of Figure 2, as well as the snapshots in Figure 1, we do find the secondary and tertiary arms launching at locations similar to these predicted values. This is evidence that these arms are likely related to ultraharmonic resonances. In the weakly nonlinear regime ( ${q}_{{\rm{th}}}\lesssim 1$ ) where ultraharmonic theory applies, the amplitude of the $n{\rm{th}}$ order perturbation decreases with increasing n following ${q}_{{\rm{th}}}^{n}.$ Accordingly, we find the tertiary and higher order arms are typically too weak to be seen. In the fully nonlinear case ( ${q}_{{\rm{th}}}\gt 1$ ), we find that the tertiary arm also cannot be easily identified, since it appears to collide and merge with the primary. Hence, in this paper, we will focus only on the morphology of the primary and secondary arm.

**Figure 2.** The left panel is the 2D Δ (Equation (7)) profile of the q = 0.001, H₀ = 0.1 model, same as the top middle panel of Figure 1. The middle panel shows the Δ map for the same model, but from a 3D simulation, and the right panel shows a face-on scattered light image corresponding to the 3D simulation. See Section 3.2 for the details of this 3D model. Similar to the way surface density is rescaled into Δ, we also rescale the scattered light image by a factor of r² to level the brightness profile. The magenta dot marks the location of the planet, and the concentric green dashed circles have radii of 0.5 and $1\;{r}_{{\rm{p}}}.$
Download figure:
Standard image High-resolution image

**Figure 2.** The left panel is the 2D Δ (Equation (7)) profile of the q = 0.001, H₀ = 0.1 model, same as the top middle panel of Figure 1. The middle panel shows the Δ map for the same model, but from a 3D simulation, and the right panel shows a face-on scattered light image corresponding to the 3D simulation. See Section 3.2 for the details of this 3D model. Similar to the way surface density is rescaled into Δ, we also rescale the scattered light image by a factor of r² to level the brightness profile. The magenta dot marks the location of the planet, and the concentric green dashed circles have radii of 0.5 and $1\;{r}_{{\rm{p}}}.$
Download figure:
Standard image High-resolution image

Two clear trends can be identified from Figure 1: (1) increasing pitch angle as the disk aspect ratio increases (right to left in Figure 1), and (2) increasing the azimuthal separation between the two arm, which we denote ${\phi }_{{\rm{sep}}}.$ as planet mass increases (top to bottom in Figure 1). The trend in pitch angle is predicted by linear perturbation theory, where the pitch angle at large radial separations from the planet can be approximated as ${\mathrm{tan}}^{-1}(c/[r| {{\rm{\Omega }}}_{{\rm{k}}}-{{\rm{\Omega }}}_{{\rm{p}}}| ])$ (e.g., Rafikov 2002). We note that, on top of this trend, there is an additional, weaker dependence on the planet mass, which is the consequence of the propagation of nonlinear ultrasonic wave (Rafikov 2002). Therefore, it is possible to evaluate planet mass using measurements of the pitch angle, but there are two challenges with this method: (1) a priori knowledge of the H profile of the disk is needed to model the pitch angle; and (2) the pitch angle is typically a strong function of its distance to the planet, and so a priori knowledge of the location of the planet is also necessary. For these reasons, we consider the trend in ${\phi }_{{\rm{sep}}}$ a more promising method for measuring planet mass.

3.1. ${\phi }_{{\rm{sep}}}$ Dependence on Planet Mass

For a given annulus, we locate the azimuthal positions of the primary and secondary arm by finding the local maxima in Σ associated with the two arms, which then allows us to construct ${\phi }_{{\rm{sep}}}$ as a function of r. On the left panel of Figure 3, we give examples of ${\phi }_{{\rm{sep}}}(r)$ using our set of models with ${H}_{0}=0.1.$ Typically, ${\phi }_{{\rm{sep}}}$ fluctuates by about ±10% over the range of $r/{r}_{{\rm{p}}}=\{0.4,\;0.7\},$ which is sufficiently small for us to clearly identify an increasing trend in ${\phi }_{{\rm{sep}}}$ as q increases. In Table 1, we report ${\phi }_{{\rm{sep}}}$ as the average over the radial range $\{{r}_{1},\;{r}_{2}\},$ where r₁ is fixed at $r/{r}_{{\rm{p}}}=0.4,$ while r₂ is the largest r at which the secondary is still clearly defined, up to $r/{r}_{{\rm{p}}}=0.7.$ We also report the maximum and minimum values of ${\phi }_{{\rm{sep}}}$ within this range as error margins.

The right panel of Figure 3 plots all of our ${\phi }_{{\rm{sep}}}$ measurements as a function of q. We find a well-defined correlation between the two that is independent of H₀. By inspecting Figures 2 and 3 of Zhu et al. (2015), we find our ${\phi }_{{\rm{sep}}}$ measurements to be consistent with both their isothermal and adiabatic simulations, which lends confidence to the accuracy of our results, and also implies that ${\phi }_{{\rm{sep}}}$ is not only insensitive to H₀, but may also be insensitive to the equation of state of the gas.

To further test this apparent insensitivity to disk thermodynamics, we perform two additional $\{q,\;{H}_{0}\}=\{0.001,\;0.1\}$ simulations, one with a globally isothermal equation of state, $c={c}_{0},$ which results in a more flared disk, and another one with $c={c}_{0}\sqrt{{r}_{{\rm{p}}}/r},$ which results in no flaring. We find ${\phi }_{{\rm{sep}}}=103^\circ$ and $96^\circ$ respectively. This suggests the dependence on the radial disk temperature profile is also negligibly weak compared to that on q.

Overall, ${\phi }_{{\rm{sep}}}$ is proportional to ${q}^{0.2},$ up to $q=0.016,$ beyond which ${\phi }_{{\rm{sep}}}$ remains approximately constant at 180°. A fit using ${\chi }^{2}$ minimization gives

$\begin{eqnarray}q\leqslant \;1.6\times {10}^{-2} & : & \\ {\phi }_{{\rm{sep}}} & = & 102^\circ {\left(\displaystyle \frac{q}{0.001}\right)}^{0.2},\end{eqnarray} \tag{ 9 }$

which typically matches our results to within 5%. We note that ${\phi }_{{\rm{sep}}}$ is the separation measured face-on. In general, observed systems will have some inclination that must be taken into account before our scaling relation can be applied. Using the uncertainty in our fitted parameters, and assuming ${\phi }_{{\rm{sep}}}$ is only a function of q, we also derive ${\rm{\Delta }}q,$ the uncertainty in q. Given ${\rm{\Delta }}{\phi }_{{\rm{sep}}},$ the measurement error in ${\phi }_{{\rm{sep}}},$ error propagation dictates that

$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}q}{q}=5\sqrt{0.002+0.001{\mathrm{ln}}^{2}\left(\displaystyle \frac{{\phi }_{{\rm{sep}}}}{102^\circ }\right)+\displaystyle \frac{{\rm{\Delta }}{\phi }_{{\rm{sep}}}^{2}}{{\phi }_{{\rm{sep}}}^{2}}}.\end{eqnarray} \tag{ 10 }$

While this establishes a basis for measuring a planet's mass using the disk's spiral structure, we need to further investigate whether the secondary arm is observable, i.e., whether it is bright enough to be detected. For this purpose, the Σ maps of Figure 1 can be misleading, because, as discussed by Juhász et al. (2015), variations in Σ provides weak observability in comparison to variations in the vertical structure of the disk. Moreover, Zhu et al. (2015) has suggested that vertical linear perturbation modes are important at high altitudes, which can lead to enhancements in density perturbation by orders of magnitude compared to the midplane. These results point toward the importance of 3D effects, which we address in the following section.

3.2. 3D Effects and Observability of the Secondary Arm

To investigate the differences between 2D and 3D disks, we additionally perform a 3D simulation for our q = 0.001 and H₀ = 0.1 model. It is in spherical coordinates, with a polar angle θ extending from 0 to 0.675 radians above the midplane. We only simulate the upper half of the disk and assume midplane symmetry. The overall resolution is reduced by about 30% to $288(r)\times 648(\phi )\times 72(\theta ).$ In 3D, ${r}_{{\rm{s}}}$ is no longer required to be of order the disk scale height, and should be a small finite value only to avoid singularity. We therefore choose ${r}_{{\rm{s}}}=0.1\;{r}_{{\rm{H}}}.$ The rest of the setup is identical to our 2D simulations.

After the simulation has reached 10 ${P}_{{\rm{p}}},$ we feed the simulation result to the Whitney et al. (2013) Monte-Carlo radiative transfer (MCRT) code to generate an H-band scattered light image with an angular resolution achievable by Very Large Telescope (VLT)/SPHERE, Gemini/GPI, and Subaru/SCExAO. The MCRT calculations largely follow the procedures presented in Dong et al. (2015a, 2015b), and is only briefly repeated here. We scale the hydrodynamics model so that the planet is at 125 AU. The central source is a typical 1 ${M}_{\odot }$ star at 1 Myr with a temperature of 4350 K and a radius of 2.325 ${R}_{\odot }.$ Dust grains are assumed to be ISM grains and well mixed with the gas (i.e., volume density of the dust is linearly proportional to the gas). The total mass of the ISM dust is assumed to be $2\times {10}^{-5}\;{M}_{\odot }.$ We experimented with different choices of total dust mass between 10⁻⁵ and ${10}^{-4}\;{M}_{\odot }$ and different filters including Y, J, H and K-band. We find the measurement of ${\phi }_{{\rm{sep}}}$ does not depend on these choices.

Full resolution synthesized polarized intensity images at H-band are produced from the MCRT calculations, and they are convolved by a Gaussian point-spread function with a full width half maximum (FWHM) of 0 farcs 05 to achieve an angular resolution comparable with the current capability of Subaru, VLT, and Gemini. Other than the intrinsic Poisson noise of the Monte-Carlo method, we do not include any additional noise to the synthesized image. The target is assumed to be at 140 pc, a typical distance of nearby star formation regions.

Figure 2 gives a side-by-side comparison for the surface density maps generated from the 2D and 3D hydrodynamics simulations (left and middle panels) along with the scattered light image translated from our 3D simulation (right panel). We find near perfect agreement between the surface density profiles of our 2D and 3D simulations. This is expected since at a large distance from the planet, the disk is effectively thin. In the scattered light image, however, some differences can be identified. First, the tertiary arm is too faint to be seen. Second, the pitch angles of both arms appear to be larger than their counterparts in surface density maps. This is due to the fact that scattered light comes from the disk surface, while surface density is determined mostly by the distribution of material at the disk midplane, and also the phenomenon observed by Zhu et al. (2015)—a bending of the propagating shock wave that points toward the star at high altitudes. Since both arms are subjected to this effect, it does not noticeably affect ${\phi }_{{\rm{sep}}}.$ Third, and most importantly, the secondary arm has a comparable brightness to the primary in the scattered light image, even though the surface density enhancement is much higher in the primary arm. This strongly suggests there is significant vertical motion near the arms, providing kinetic, rather than pressure, support to the local vertical structure, allowing the scattering surface of the secondary arm to rise to a similar height as the primary. A future study is required to thoroughly understand this feature.

One aspect of 3D effects that we have not taken into account in this study is the presence of a vertical temperature gradient. In a thermally stratified disk, density waves can become either deflected toward the disk surface, or confined in the midplane, depending on the direction of the temperature gradient (Lubow & Ogilvie 1998; Lee & Gu 2015). Since we do not expect the primary and secondary arms to behave differently in this process of wave channeling, we do not expect ${\phi }_{{\rm{sep}}}$ to be dependent on this effect. However, we do expect it to have a significant impact on the brightness of the arms, making them more prominent if they are channeled toward the surface, or dimmer vice versa.

Finally, we measure ${\phi }_{{\rm{sep}}}$ in the scattered light image. As shown in Figure 4, the two arms have a near constant azimuthal separation, and we find the average value of ${\phi }_{{\rm{sep}}}$ between $r/{r}_{{\rm{p}}}=\{0.45,\;0.65\}$ to be ∼96°. Inserting this value into Equation (9), it predicts q to be about $7.4\times {10}^{-4},$ which successfully recovers the planet's mass to within 30%.

3.3. Example: Application to SAO 206462

Due to the simplicity of our scaling relation, given an image of the spirals, one can predict the mass of the planet responsible for the spiral structure with minimal effort. Here we use SAO 206462 as an example.

The circumstellar disk around SAO 206462 is a convenient choice to apply our results because it is nearly face-on with an inclination of only $11^\circ$ (Muto et al. 2012). Therefore we may neglect projection effects and directly measure the separation of the two arms using Figure 1 of Garufi et al. (2013). We find ${\phi }_{{\rm{sep}}}\sim 130^\circ .$ This translates to $q\sim 0.0034$ using Equation (9) to within 30% accuracy, or ∼6 Jupiter masses, given the star's mass is $\sim 1.7\;{M}_{\odot }$ (Müller et al. 2011). We can also infer that the arm on the east side is the primary arm; therefore, this planet should be located toward the trip of that arm.

4. CONCLUSIONS AND DISCUSSIONS

We performed a parameter study for the azimuthal separation between the primary and secondary spiral arm, ${\phi }_{{\rm{sep}}},$ in the multi-armed spiral structure generated by a single planet on a fixed circular orbit, over a range of planet-to-star mass ratio, q, from $6.25\times {10}^{-5}$ to $6.4\times {10}^{-2},$ and disk aspect ratio at the planet's orbit, H₀, from 0.04 to 0.16. We found an empirical scaling that allows one to infer q directly from ${\phi }_{{\rm{sep}}}$ (Equation (9)) for planet mass companions; and also found that ${\phi }_{{\rm{sep}}}\sim 180^\circ$ for brown dwarf mass companions. Using a near-infrared scattered light image synthesized by a combination of MCRT and 3D hydrodynamics simulation, we showed that this method is capable of measuring the mass of a Jupiter mass planet to within 30% when the disk is face-on. Finally, we applied our scaling relation to SAO 206462, the most face-on two-armed spiral so far, and inferred that it contains a planet with about 6 Jupiter masses. In addition to observational applications, the spiral morphology described in this paper can serve as both motivation and pointers for further theoretical studies. In the following, we briefly discuss the implications of our results on nonlinear disk–planet interaction theory.

In the beginning of this paper we suggested that ultraharmonic density waves may be responsible for the multi-armed spiral structures. Even if this is true, without an in-depth theoretical understanding, it is difficult to predict what ${\phi }_{{\rm{sep}}}$ should be, since the multi-armed structure is the result of the interference between both the linear and ultraharmonic density waves, and to evaluate it requires precise knowledge of how the amplitudes of different azimuthal modes compare to each other. Nonetheless, it is surprising that ${\phi }_{{\rm{sep}}}$ depends on q only, and not the disk temperature, because the linear density perturbation, which drives ultraharmonic waves, certainly depends on the disk temperature as well as q. Adding further complication to the picture, the location of the primary arm also changes with planet mass in the nonlinear regime, due to the spiral wave steepening into a shock (Rafikov 2002). It may well be that to explain the scaling we found in this paper, both ultraharmonic waves and nonlinear wave propagation need to be taken into account.

J.F. gratefully acknowledges support from the Center for Integrative Planetary Science at the University of California, Berkeley. We thank Pawel Artymowicz for computational resources and helpful feedback. We also thank Eugene Chiang for fruitful discussions during the process of writing this manuscript.

INFERRING PLANET MASS FROM SPIRAL STRUCTURES IN PROTOPLANETARY DISKS

Article metrics

Permissions

Author affiliations

Author notes

ORCID iDs

Dates

ABSTRACT

1. INTRODUCTION

2. NUMERICAL METHOD

2.1. Disk Model

2.2. Setup

3. RESULTS

3.1. ${\phi }_{{\rm{sep}}}$ Dependence on Planet Mass

3.2. 3D Effects and Observability of the Secondary Arm

3.3. Example: Application to SAO 206462

4. CONCLUSIONS AND DISCUSSIONS

Footnotes

INFERRING PLANET MASS FROM SPIRAL STRUCTURES IN PROTOPLANETARY DISKS

Article metrics

Permissions

Share this article

Author affiliations

Author notes

ORCID iDs

Dates

ABSTRACT

1. INTRODUCTION

2. NUMERICAL METHOD

2.1. Disk Model

2.2. Setup

3. RESULTS

3.1. {\phi }_{{\rm{sep}}} Dependence on Planet Mass

3.2. 3D Effects and Observability of the Secondary Arm

3.3. Example: Application to SAO 206462

4. CONCLUSIONS AND DISCUSSIONS

Footnotes

3.1. ${\phi }_{{\rm{sep}}}$ Dependence on Planet Mass