THE STRUCTURE OF SPIRAL SHOCKS EXCITED BY PLANETARY-MASS COMPANIONS

Recent high-resolution direct imaging observations have revealed spiral structure in three protoplanetary disks around Herbig Ae/Be stars: SAO 206462 (Muto et al. 2012; Garufi et al. 2013), MWC 758 (Grady et al. 2013; Benisty et al. 2015), and HD 100546 (Currie et al. 2014). The polarized intensity has been measured in these observations to gain higher contrast between the disk and the central star. While the thermal emission from the central star is unpolarized, the scattered light from the disk is polarized. In these near-infrared (near-IR) polarized intensity images, two spiral arms with roughly 180°rotational symmetry are present in both SAO 206462 and MWC 758, similar to the grand design in a spiral galaxy (e.g., the Whirlpool Galaxy M51). The spiral arms also exhibit a high contrast against the background disk. The polarized intensity of the spiral arm is several times higher than that of the region outside the spiral arm. It should also be noted that, since the dust scattering opacity is quite large, these observations only probe structure high up at the disk atmosphere (e.g., several disk scale heights) where the last dust scattering surface is.

In addition to spiral patterns, these three disks also have gaps or holes (SAO 206462 with a submillimeter cavity of 46 AU, MWC 758 with a cavity of 73 AU from Andrews et al. 2011, and HD 100546 with a cavity of 10 AU revealed by SED fitting from Bouwman et al. 2003), which indicates that they are members of the protoplanetary disk class called transitional disks (Espaillat et al. 2014). One scenario to explain both spiral patterns and gaps/holes is that these disks harbor low-mass companions (e.g., young planets) which can open gaps and excite spiral waves at the same time (e.g., Baruteau et al. 2014).

However, there are two difficulties in explaining the observed spiral patterns using planet-induced spiral wakes. First, the large pitch angle of all the observed spiral arms suggests that the disk has a relative high temperature (e.g., ∼250 K at 70 AU for MWC 758, Benisty et al. 2015). In linear theory, spiral waves are basically sound waves in disks, and the pitch angle of the spiral arms is directly related to the sound speed in the disk. Using the linear theory, the best fit models for both SAO 206462 (Muto et al. 2012) and MWC 758 (Grady et al. 2013; Benisty et al. 2015) suggest that the disk aspect ratio (H/R with H ≡ c_s/Ω) at R ∼ 100 AU is around 0.2, which is too large for any realistic disk structure. For example, even if the stellar irradiation is perpendicular to the disk surface,⁵ the maximum disk temperature due to the stellar irradiation is $\sigma T{(R)}^{4}={L}_{*}/4\pi {R}^{2}$ so that $H/R\propto {T}^{1/2}\propto {L}_{*}^{1/8}.$ Assuming a 2 M_⊙ central star with 10 L_⊙ luminosity, the maximum temperature is ∼70 K at 100 AU and H/R is only ∼0.1. Since $H/R\propto {L}_{*}^{1/8},$ it is very difficult to make H/R ∼ 0.2.

Second, the observed spiral arms exhibit much higher brightness contrasts than suggested by the synthetic observations based on two-dimensional (2D) planet–disk simulations. Juhász et al. (2015) calculated the polarized scattered light images by combining 2D hydrodynamical simulations with three-dimensional (3D) Monte-Carlo radiative transfer (MCRT) simulations. Vertical hydrostatic equilibrium has been assumed to extend the 2D simulation to the third-dimension (the vertical direction). They find that a relative change of about 3.5 on the spiral arms in the surface density is required for the spirals to be detectable. This value is a factor of eight higher than what is seen in their hydrodynamical simulations.

In this paper, we first point out that the pitch angle formula derived from the linear theory, which has been used in almost all previous spiral arm modeling efforts, does not apply to the high planet mass cases. Spiral wakes that are excited by high mass planets (e.g., 1 M_J) become spiral shocks that propagate at speeds faster than the local sound speed (Goodman & Rafikov 2001; Rafikov 2002). Thus, the pitch angle difficulty above can be alleviated by considering the nonlinear extension of the spiral shock theory. We also show that spiral arms (especially the inner arms) have complicated non-hydrostatic 3D structure. Such structure can lead to strong density perturbation at the disk surface resulting in a corrugated shape of its atmosphere. Since near-IR observations only probe the shape of disk surface, this effect alleviates the second difficulty mentioned above. In Dong et al. (2015a), we have combined MCRT simulations with hydrodynamical simulations from this paper and demonstrated that planet-induced inner spiral arms can explain recent near-IR direct imaging observations of SAO 206462 and MWC 758. We note that, since the planets that we have proposed are outside the spiral arms, we cannot explain the gaps discovered at small radii in these transitional disks. Other mechanisms, e.g., another planet or photoevaporation, are needed to explain these gaps.

Before we introduce our numerical method in Section 3, we provide the theoretical background in Section 2. The shape of the spiral wakes will be studied in Section 4, and their 3D structure will be presented in Section 5. After a short discussion in Section 6, we summarize our results in Section 7.

As a result of planet–disk interactions, a spiral arm forms due to the constructive interference of density wakes with different azimuthal wavenumbers m excited by the planet at Lindblad resonances. In the linear density wave theory, the mth Fourier component of the planet potential excites the density wave having m spiral arms

$$\begin{eqnarray}&&\delta (R,\phi ,t)={\delta }_{0}(R){e}^{i[\int {k}_{R}(R){dR}+m(\phi -{{\rm{\Omega }}}_{p}t)]}\end{eqnarray} \tag{ 1 }$$

where δ is any perturbed quantity associated with the wave, δ₀(R) is its complex amplitude, k_R(R) is the radial wave vector, and Ω_p is the planet orbital frequency. Thus, the wave has the same phase along the curve satisfying dR/dϕ = −m/k_R(R). The pitch angle (β) of the equal phase curve satisfies $\mathrm{tan}\beta$ = $| {dR}/({Rd}\phi )| ,$ so $\beta$ = ${\mathrm{tan}}^{-1}| m/[{k}_{R}(R)R]| .$ Using the dispersion relationship for density waves in the large m limit and far from the launching point, ${m}^{2}{({\rm{\Omega }}(R)-{{\rm{\Omega }}}_{p})}^{2}$ ≈ ${c}_{s}^{2}{k}_{R}{(R)}^{2},$ we have $\beta ={\mathrm{tan}}^{-1}[{c}_{s}/(R| {\rm{\Omega }}(R)-{{\rm{\Omega }}}_{p}| )].$ Because β is independent of m, different m modes can constructively interfere to form the one armed spiral wake (Ogilvie & Lubow 2002). If the equal phase curve is integrated from the planet's position (R_p, ϕ_p), the shape of the wake far from R_p is given by Rafikov (2002) and Muto et al. (2012) as

$$\begin{eqnarray}\phi (R) & = & {\phi }_{p}-\displaystyle \frac{{\rm{sgn}}(R-{R}_{p})}{{h}_{p}}\\ & & \times \left[{\left(\displaystyle \frac{R}{{R}_{p}}\right)}^{1+\eta }\left\{\displaystyle \frac{1}{1+\eta }-\displaystyle \frac{1}{1-\alpha +\eta }{\left(\displaystyle \frac{R}{{R}_{p}}\right)}^{-\alpha }\right\}\right.\\ & & -\left.\left(\displaystyle \frac{1}{1+\eta }-\displaystyle \frac{1}{1-\alpha +\eta }\right)\right]\end{eqnarray} \tag{ 2 }$$

where h_p = H/R is the disk aspect ratio at R_p, Ω(R) ∝ R^−α, and the sound speed c_s(R) ∝ R^−η.

However, when the planet is massive enough, the above linear density wave theory breaks down. Linear waves excited by planets will steepen to shocks (Goodman & Rafikov 2001; Rafikov 2002; Dong et al. 2011; Duffell & MacFadyen 2012; Zhu et al. 2013) after they propagate over a distance

$$\begin{eqnarray}&&| {x}_{\mathrm{sh}}| \approx 0.93{\left(\displaystyle \frac{\gamma +1}{12/5}\displaystyle \frac{{M}_{p}}{{M}_{\mathrm{th}}}\right)}^{-2/5}H\end{eqnarray} \tag{ 3 }$$

where γ is the adiabatic index, and M_th is the disk thermal mass

$$\begin{eqnarray}&&{M}_{\mathrm{th}}\equiv \displaystyle \frac{{c}_{s}^{3}}{G{{\rm{\Omega }}}_{p}}\approx 1{M}_{J}{\left(\displaystyle \frac{{h}_{p}}{0.1}\right)}^{3}\left(\displaystyle \frac{{M}_{*}}{{M}_{\odot }}\right).\end{eqnarray} \tag{ 4 }$$

When M_p > M_th, the spiral waves will immediately become spiral shocks after they are excited around the planet. Unlike the linear wake that follows Equation (2), the spiral shock will expand away from the trajectory predicted by Equation (2). Thus, if there is a massive planet in the disk, using Equation (2) to fit the shape of the spiral shocks will predict an incorrect disk aspect ratio and temperature.

3.1. Method

To study density wakes/shocks excited by planets, we have carried out both 2D and 3D hydrodynamical simulations using Athena++. Athena++ is a newly developed grid based code using a higher-order Godunov scheme for MHD and the constrained transport (CT) to conserve the divergence-free property for magnetic fields (J. M. Stone et al. 2015, in preparation). But in this paper, we do not include magnetic fields and only solve hydrodynamical equations using Athena++. Compared with its predecessor Athena (Gardiner & Stone 2005, 2008; Stone et al. 2008), Athena++ is highly optimized and uses flexible grid structures, allowing global numerical simulations spanning a large radial range. Furthermore, the geometric source terms in curvilinear coordinates (e.g., in cylindrical and spherical-polar coordinates) are carefully implemented so that angular momentum is conserved exactly (to machine precision), which makes the code ideal for global disk simulations.

Our simulations use the adiabatic equation of state (EOS) with the adiabatic index γ = 1.4. A simple orbital cooling scheme has been applied to mimic the radiative cooling process in disks. In 3D simulations, we have adopted

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dt}}=-\displaystyle \frac{E-{c}_{v}\rho {T}_{\mathrm{irr}}}{{t}_{\mathrm{cool}}},\end{eqnarray} \tag{ 5 }$$

where ρ and E are the density and the internal energy per unit volume, while in 2D simulations, we have adopted

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dt}}=-\displaystyle \frac{E-{c}_{v}{\rm{\Sigma }}{T}_{\mathrm{irr}}}{{t}_{\mathrm{cool}}}\end{eqnarray} \tag{ 6 }$$

where Σ and E are the disk surface density and the internal energy per unit area. ${c}_{v}\equiv k/(\mu {m}_{u}(\gamma -1))$ is the heat capacity per unit mass, k is the Boltzmann constant, μ is the mean molecular weight, and m_u is the atomic mass unit. The cooling time t_cool can be written in the dimensionless form as T_cool = t_coolΩ(R). We fix T_cool to be a constant in each simulation. With this scheme, the disk temperature is relaxed to the background disk temperature (T_irr) determined by stellar irradiation. In our simulations, T_irr is set to be the initial disk temperature. To estimate T_cool in a realistic disk, we use the gray atmosphere approximation (Hubeny 1990) for the radiative cooling,

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dt}}=-\displaystyle \frac{16}{3}\sigma ({T}_{c}^{4}-{T}_{\mathrm{irr}}^{4})\displaystyle \frac{\tau }{1+{\tau }^{2}},\end{eqnarray} \tag{ 7 }$$

where σ is the Stefan–Boltzmann constant, τ = (Σ/2)κ_R is the optical depth in the vertical direction, κ_R is the Rosseland mean opacity, and T_c is the midplane temperature. Assuming E = c_vΣT_c and using Equations (6) and (7), we can derive

$$\begin{eqnarray}&&{t}_{\mathrm{cool}}=\displaystyle \frac{3{\rm{\Sigma }}{c}_{v}}{16\sigma ({T}_{c}^{2}+{T}_{\mathrm{irr}}^{2})({T}_{c}+{T}_{\mathrm{irr}})}\displaystyle \frac{1+{\tau }^{2}}{\tau }.\end{eqnarray} \tag{ 8 }$$

Approximating the polynomial of T_c and T_irr in the denominator of Equation (8) with [max(T_c, T_irr)]³ and assuming the central star is 1 M_⊙, we have

$$\begin{eqnarray}{T}_{\mathrm{cool}} & = & 0.002\left(\displaystyle \frac{{\rm{\Sigma }}}{10\;{\text{g cm}}^{-2}}\right){\left(\displaystyle \frac{100\;{\rm{AU}}}{R}\right)}^{1.5}\displaystyle \frac{{(60\;{\rm{K}})}^{3}}{{[{\rm{max}}({T}_{c},{T}_{\mathrm{irr}})]}^{3}}\\ & & \times \displaystyle \frac{1+{\tau }^{2}}{\tau }.\end{eqnarray} \tag{ 9 }$$

Thus, T_cool can vary dramatically at different radii in disks. Using the minimum mass solar nebular model, T_cool is ∼10⁵ at 1 AU and 10⁻² at 100 AU. We have carried out three sets of simulations with T_cool = 10⁻⁵, 1, and 100. They are respectively labeled as ISO, T1, T2 at the end of their names in Table 1. Simulations with fast cooling (T_cool = 10⁻⁵) are equivalent to locally isothermal simulations. Simulations with T_cool = 100 are basically adiabatic simulations considering the timescale of the simulations is several tens of orbits, and simulations with T_c = 1 should be between isothermal and adiabatic simulations. Somewhat surprisingly, we find that simulations with T_cool = 1 are qualitatively similar to those with T_cool = 100. This similarity suggests that the spiral waves in disks with T_cool = 1 behave adiabatically. We think that this is due to the short timescale for the flow to move across the spiral wake. The spiral wake has a typical width smaller than the disk scale height, and the background flow moves across the wake at nearly Keplerian speed. Thus, when the fluid travels in and out of the wave/shock, its response time is much smaller than the orbital time and should behave adiabatically even in disks with T_cool = 1. Considering this similarity, in most part of the paper, we only show results with T_cool = 10⁻⁵ and 1.

Table 1.  Models

2D
Run	Mp	Tcool	Domain	Resolution
	MJ		R	R × ϕ
CM1ISO	0.01	10−5	[0.2, 10]	1280 × 2048
CM1T1	0.01	1	[0.2, 10]	1280 × 2048
CM1T100	0.01	100	[0.2, 10]	1280 × 2048
CM2ISO	1	10−5	[0.2, 10]	1280 × 2048
CM2T1	1	1	[0.2, 10]	1280 × 2048
CM2T100	1	100	[0.2, 10]	1280 × 2048
CM3ISO	6	10−5	[0.2, 10]	1280 × 2048
CM3T1	6	1	[0.2, 10]	1280 × 2048
CM3T100	6	100	[0.2, 10]	1280 × 2048
3D
Run	Mp	Tcool	Domain	Resolution
	MJ		r × θ	r × θ × ϕ
STHIN	0.000316	10−5	[0.5, 2] × $[\frac{\pi }{2}-0.2,\frac{\pi }{2}+0.2]$	456 × 128 × 2048
SM1ISO	0.01	10−5	[0.3, 3] × $[\frac{\pi }{2}-0.6,\frac{\pi }{2}+0.6]$	256 × 128 × 688
SM1T1	0.01	1	[0.3, 3] × $[\frac{\pi }{2}-0.6,\frac{\pi }{2}+0.6]$	256 × 128 × 688
SM1T100	0.01	100	[0.3, 3] × $[\frac{\pi }{2}-0.6,\frac{\pi }{2}+0.6]$	256 × 128 × 688
SM2ISO	1	10−5	[0.3, 3] × $[\frac{\pi }{2}-0.6,\frac{\pi }{2}+0.6]$	256 × 128 × 688
SM2T1	1	1	[0.3, 3] × $[\frac{\pi }{2}-0.6,\frac{\pi }{2}+0.6]$	256 × 128 × 688
SM2T100	1	100	[0.3, 3] × $[\frac{\pi }{2}-0.6,\frac{\pi }{2}+0.6]$	256 × 128 × 688
SM3ISO	6	10−5	[0.3, 3] × $[\frac{\pi }{2}-0.6,\frac{\pi }{2}+0.6]$	256 × 128 × 688
SM3T1	6	1	[0.3, 3] × $[\frac{\pi }{2}-0.6,\frac{\pi }{2}+0.6]$	256 × 128 × 688
SM3T100	6	100	[0.3, 3] × $[\frac{\pi }{2}-0.6,\frac{\pi }{2}+0.6]$	256 × 128 × 688
SM2ISOlong	1	10−5	[0.3, 3] × $[\frac{\pi }{2}-0.6,\frac{\pi }{2}+0.6]$	256 × 128 × 688
SM2ISOcyl	1	10−5	[0.3, 3] × [−0.8, 0.8]a	296 × 176 × 690b
SM2ISOconstT	1	10−5	[0.3, 3] × $[\frac{\pi }{2}-0.6,\frac{\pi }{2}+0.6]$	256 × 128 × 688

Notes.

^aThis is the domain size in the R × Z direction since cylindrical coordinates have been used. ^bThe resolution is for R × Z × ϕ with the cylindrical coordinate system.

Download table as:  ASCIITypeset image

We have also varied the planet mass to be 0.01, 1, and 6 M_J in our main set of simulations, which are labeled as M1, M2, and M3 in their names respectively. Here, we have defined M_J as 0.001 of the central star's mass. The thermal mass (Equation (4)) for the h = 0.1 disk is ∼M_J. Thus, waves excited by a 0.01 M_J planet are in the linear regime and waves from a 6 M_J planet are in the highly nonlinear regime. To compare with Figure 2 in Tanaka et al. (2002), we have also carried out a thin disk simulation with H_p/R_p = 10^−1.5 (STHIN in Table 1). The thermal mass for such a thin disk is only M_th = 0.0316 M_J. Thus, in order to ensure that the waves are in the linear regime, we choose the planet mass of 0.01 M_th = 3.16 × 10⁻⁴M_J in this thin disk simulation. To avoid the divergence of planet potential, a smoothing length of 0.1 R_p has been applied for M2 and M3 cases (0.1 R_p is close to the disk scale height and the Hill radius of the planet). For the thin disk case which has a very small mass planet, we choose a smoothing length of 6 × 10⁻³R_p, roughly the length of two grid cells. For the low mass planet cases (M1), a smoothing length of 0.02 R_p, which is also roughly the length of two grid cells in these simulations, has been adopted. Planets are fixed in circular orbits at R = 1, and the indirect potential, which is due to the center of the coordinate system is at the star instead of the center of the mass, has been included. We have run the simulations for 10 planetary orbits. We choose this timescale because it is longer than the sound crossing time throughout the whole disk so that density waves/shocks have established, while it is shorter than the gap opening timescale to avoid complicated gap structures (e.g., vortices at the gap edges) and other longterm effects (e.g., radial buoyancy waves, Richert et al. 2015). To verify that the revealed wave mechanics still hold in long terms we have run one simulation for 120 orbits (SM2ISOlong in Table 1), which will be discussed in Section 6.1. A constant α viscosity with α = 10⁻⁴ has been applied in our main sets of simulations.

At inner and outer boundaries, all quantities are fixed at the initial states. For a numerical code using Godunov scheme, which calculates the flux by decomposing wave characteristics at the left and right grid cells, such boundary can absorb wave characteristics coming from the active zones and limit waves traveling from the ghost zones to active zones. Such boundary condition shows little wave reflection and is similar to the non-reflecting boundary condition (Godon 1996) used in the FARGO code (Masset 2000). The detailed code comparison is given in Appendix C of Zhu et al. (2014).

3.2. 2D Simulations

Compared with the 3D simulations in the next subsection, 2D simulations allow us to study density wakes in a bigger domain using a higher numerical resolution. The initial radial profile of the disk is

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

$$\begin{eqnarray}&&{T}_{0}(R,\phi )={T}_{0}({R}_{0}){\left(\displaystyle \frac{R}{{R}_{0}}\right)}^{-1/2}.\end{eqnarray} \tag{ 11 }$$

We choose R₀ = 1, Σ₀(R₀) = 1, and T₀(R₀) = 0.01 to make ${(H/R)}_{R={R}_{0}}=0.1.$

Cylindrical coordinates have been adopted. To make every grid cell have equal length in the radial and azimuthal direction throughout the whole domain, the grids are uniformly spaced in log(R) from R = 0.2 to 10, and uniformly spaced from 0 to 2π in the ϕ direction. Our standard resolution is 1280 in the R direction and 2048 in the ϕ direction, which is equivalent to 32 grids per H at R = 1 in both directions. In Table 1, 2D runs are denoted with a "C" in front of the model names, while 3D runs are denoted with a "S" in front of the model names.

3.3. 3D Simulations

To study the 3D structure of density wakes/shocks, we have run 3D hydrodynamical simulations in spherical polar coordinates except for one case in cylindrical coordinates. The initial density profile of the disk at the disk midplane is

$$\begin{eqnarray}&&{\rho }_{0}(R,z=0)={\rho }_{0}({R}_{0},z=0){\left(\displaystyle \frac{R}{{R}_{0}}\right)}^{p},\end{eqnarray} \tag{ 12 }$$

and the temperature is constant on cylinders

$$\begin{eqnarray}&&{T}_{0}(R,z)={T}_{0}({R}_{0}){\left(\displaystyle \frac{R}{{R}_{0}}\right)}^{q}.\end{eqnarray} \tag{ 13 }$$

We want to emphasize that R should not be confused with r. In this paper, we use (R, ϕ, z) to represent positions in cylindrical coordinates while using (r, θ, ϕ) for spherical polar coordinates. ϕ represents the azimuthal direction (the direction of disk rotation) in both coordinate systems. Considering the disk structure is more natural to be described in cylindrical coordinates, we have transformed 3D simulation results from spherical polar coordinates to cylindrical coordinates. Most results presented below are plotted in cylindrical coordinates with R representing the distance to the axis of the disk, even though most simulations are carried out in spherical polar coordinates.

Hydrostatic equilibrium in the $r-\theta$ plane requires that (e.g., Nelson et al. 2013)

$$\begin{eqnarray}{\rho }_{0}(R,z) & = & {\rho }_{0}({R}_{0},z=0){\left(\displaystyle \frac{R}{{R}_{0}}\right)}^{p}\\ & & \times {\rm{exp}}\left[\displaystyle \frac{{GM}}{{c}_{s}^{2}}\left(\displaystyle \frac{1}{\sqrt{{R}^{2}+{z}^{2}}}-\displaystyle \frac{1}{R}\right)\right],\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}&&{\rm{\Omega }}(R,z)={{\rm{\Omega }}}_{K}{\left[(p+q){\left(\displaystyle \frac{H}{R}\right)}^{2}+(1+q)-\displaystyle \frac{{qR}}{\sqrt{{R}^{2}+{z}^{2}}}\right]}^{1/2},\end{eqnarray} \tag{ 15 }$$

where ${c}_{s}=\sqrt{p/\rho }$ is the isothermal sound speed at R, ${{\rm{\Omega }}}_{K}=\sqrt{{{GM}}_{*}/{R}^{3}},$ and H = c_s/Ω_K as defined before.

We choose p = −2.25 and q = −1/2 in our main sets of simulations so that Σ₀ ∝ R⁻¹, similar to 2D simulations above. ${\rho }_{0}({R}_{0},z=0)$ is 1 and H/R is 0.1 at R = R₀. The grids are uniformly spaced in log(r), θ, ϕ with 256 × 128 × 688 grid cells in the domain of [log(0.3), log(3)] × [π/2 − 0.6, π/2 + 0.6 ] × [0, 2π] for the main sets of simulations. In runs with T_cool = 1 and 100, the cooling time decreases exponentially beyond z = 3H with ${T}_{\mathrm{cool}}(z)={T}_{\mathrm{cool}}{\rm{exp}}(-({z}^{2}/{H}^{2}-{3}^{2}))$ to mimic fast cooling at the disk surface (D'Alessio et al. 1998). Numerically, this treatment also maintains better hydrostatic equilibrium at the disk surface.

The boundary condition in the θ direction is chosen that v_r = v_θ = 0 in the ghost zones. We set v_ϕ and T in the ghost zones having the same values as the last active zones. Density in the ghost zones is set to be

$$\begin{eqnarray}&&\rho ({\theta }_{g})=\rho ({\theta }_{a}){\left|\displaystyle \frac{{\rm{sin}}({\theta }_{g})}{{\rm{sin}}({\theta }_{a})}\right|}^{{v}_{\phi }^{2}/T}\end{eqnarray} \tag{ 16 }$$

to maintain hydrostatic equilibrium in the θ direction, where θ_g and θ_a are the θ coordinates of the ghost and last active zones. We have also tried the boundary condition which sets the quantities in the ghost zones as the initial values, and found that the results are not affected by the choice of boundary conditions.

To further study the numerical effect from the boundary, we have applied a wave damping zone (de Val-Borro et al. 2006) operating at both r and θ boundaries in the run SM2ISOlong. The wave damping zone in the r direction is from R_in to 1.25R_in and from 0.84R_out to R_out. The damping zone in the θ direction is from θ_in to ${\theta }_{\mathrm{in}}+0.1$ and from ${\theta }_{\mathrm{out}}-0.1$ to θ_out. And R_in, R_out, θ_in, and θ_out are the boundary of the simulation domain. In these damping zones, the physical quantities are relaxed to the initial states on a timescale varying from infinity at the damping zone edge 1.25R_in, 0.84R_out, ${\theta }_{\mathrm{in}}+0.1,$ and ${\theta }_{\mathrm{out}}-0.1$ to a timescale of 0.1 orbit at the boundary of the simulation domain. The damping zone gradually damps waves traveling to the boundary of the simulation domain. Besides the wave damping zone and the long timescale, the run SM2ISOlong is different from our main set of simulations in other ways. We ramp the planet mass linearly for 10 orbits in run SM2ISOlong to test how our results will be affected if we insert the planet gradually in the disk. We also choose α = 0 in run SM2ISOlong to confirm that the small viscosity (α = 10⁻⁴) in the main set simulations will not affect the results.

To verify our results, especially for the 3D structure of spiral shocks, we have used Athena with cylindrical coordinates (run SM2ISOcyl) to carry out a 3D simulation having the same disk set-up as in SM2ISO. This Athena simulation (SM2ISOcyl) is different from the Athena++ simulation (SM2ISO) in several ways. First, SM2ISOcyl uses the Corner Transport Upwind Integrator which is different from the Van-Leer Integrator used in Athena++. Second, cylindrical coordinates with uniform radial grids have been used in SM2ISOcyl while spherical-polar coordinates with logarithmic radial grids have been used in SM2ISO. The detailed disk set-up and the boundary conditions can be found in Zhu et al. (2014). The results will be presented in Section 6.1. Overall, SM2ISOcyl confirms the results in SM2ISO, which greatly limits the chance that our results are due to numerical artifacts.

Volume rendering of δρ/ρ₀ in simulation SM1ISO is shown in Figure 1. δρ is the density difference between 10 orbits and the initial condition, and ρ₀ is the initial density at that position. Thus, δρ/ρ₀ highlights the density perturbation (e.g., spiral shocks) in the disk. In this paper, we use "spiral shocks" to refer to peaks of the density wakes and are associated with spiral arms seen in observations. It is apparent in Figure 1 that the spiral shocks are not perpendicular to the disk midplane and they have complicated 3D structure. In the figure, both the inner arms inside the planet and the outer arms outside the planet curl toward the central star at higher altitudes. This curled 3D shock structure will be studied in more detail in Section 5, while in this section we focus on the shape of the spiral wakes in the horizontal plane.

Zoom In Zoom Out Reset image size

Figure 2 shows the shape of the spiral wakes in both 2D and 3D simulations. The x-axis is plotted in log R, so that the pitch angle of the spiral wake can be easily estimated by using its slope in the figure (d log R/d ϕ = tan β).

Zoom In Zoom Out Reset image size

When a very low mass planet is present in a disk, it excites density waves that are in the linear regime. The linear theory for density waves in the 2D $R-\phi$ plane (Equation (2)) can accurately describe the shape of the excited spiral wakes in 2D simulations. This is demonstrated in the upper left panel (CM1ISO) of Figure 2 where Equation (2) with α = 3/2 and η = 1/4 fits the peak of the density wakes very well.⁶ Even at the midplane of 3D simulations, Equation (2) still provides a good fit to the density wakes (the upper middle panel).

However, the shape of the spiral arms at the disk surface is affected by the 3D structure of density wakes. At the disk surface in 3D simulations (even in the linear regime, shown in the upper right panel of Figure 2), both inner and outer arms are at smaller R than Equation (2) due to the tilted shock shape in Figure 1. When these shocks are far away from the planet, they are more tilted toward the central star at higher altitudes. This leads to the inner spiral arms becoming slightly more open (with a larger pitch angle) and the outer spiral arms becoming slightly less open (with a smaller pitch angle) than Equation (2) would predict.

When the planet has a mass larger than M_th (middle and bottom panels of Figure 2), it can launch spiral shocks immediately around the planet, and the shape of spiral shocks can deviate from the trajectory predicted by Equation (2) significantly. Shocks excited by a more massive planet deviate from linear theory more and they have larger pitch angles. As shown in Figure 2, spiral shocks in the 6 M_J cases (bottom panels) are more open and deviate from the prediction of linear theory (dotted curves) more than shocks in the 1 M_J cases (middle panels). The deviation from the linear theory has also been seen in previous simulations, e.g., Figures 2 and 10 of de Val-Borro et al. (2006), but it has not been explored and the physical reason for the deviation is left to be unexplained.

This deviation from linear theory shown in Figure 2 is consistent with the predictions from the weakly nonlinear density wave theory by Goodman & Rafikov (2001) and Rafikov (2002). In weakly nonlinear theory, the spiral shock can expand in both azimuthal directions away from Equation (2), and, at each radius, the shock density profile along the azimuthal direction is N-shaped (Figure 2 in Goodman & Rafikov 2001). The N-shaped shock profile expands in the azimuthal direction at a speed which is proportional to the normalized amplitude of the shock ($({{\rm{\Sigma }}}_{\mathrm{shock}}-{{\rm{\Sigma }}}_{0})/{{\rm{\Sigma }}}_{0}$). Thus, the higher is the planet mass, the stronger are the shocks and these shocks expand faster away from Equation (2). Then, the spiral shock has a larger pitch angle as a result.

Similar to the 0.01 M_J case, the inner spiral shocks in 1 and 6 M_J cases are even more open at z = 3H than at the midplane, while the outer arms become less open at the disk surface. For outer spiral arms, this 3D effect compensates the increased pitch angel due to the shock expansion, and coincidently the outer arms almost overlap with the prediction from linear theory.

Another important feature shown in Figure 2 is that, besides the primary inner arm which originates from the planet, a secondary inner spiral arm appears with some azimuthal shift from the primary arm. For some cases (e.g., CM2ISO), we even see a tertiary arm at the very inner disk. The secondary spiral arm has also been seen in previous simulations having massive planets, e.g., Figure 2 in Kley (1999) and Figure 10 of de Val-Borro et al. (2006). However, it has hardly been explored in earlier simulations. Figure 2 shows that, even with a very low mass planet (0.01 M_J, the upper panels of Figure 2), another density peak (the secondary arm) emerges close to the primary inner arm with low density region (the rarefaction wave of the primary arm) in between. After the secondary inner arm appears close to the primary arm, it can become shock during the propagation and later it will become N-shaped which is similar to the primary arm. Then in some cases, a tertiary inner arm appears at the rarefaction wave part of the secondary arm. For 0.01 M_J cases, the primary and secondary inner arms are separated by δϕ ∼ 1 at R = 0.3. When the planet gets more massive, the secondary inner arm is excited earlier and the separation between the primary and secondary arm increases. In 1 M_J cases, the two inner arms are roughly separated by δϕ ∼ 2, and in 6 M_J cases, the two arms are roughly separated by δϕ ∼ 3. This has an important application that we can use the separation between two arms to estimate the mass of the embedded planet. Similar to the primary inner arm, the secondary inner arm is also stronger at the disk surface than at the disk midplane. For outer arms, the secondary arm appears in disks that have massive planets embedded, but the secondary outer arm is less apparent than the primary outer arm. More discussions on the 3D structure of secondary arms will be presented in Section 5.

Spiral wakes/shocks are slightly more open in a disk whose EOS is not isothermal (Figure 3). This is because density waves propagate slightly faster in a fluid with a non-isothermal EOS than in a fluid at the same temperature with the isothermal EOS (e.g., Goodman & Rafikov 2001). In Figure 3, even with a moderate cooling rate (T_cool = 1), the spiral wakes excited by a low mass planet (0.01 M_J) can only be fitted by Equation (2) using a larger disk scale height that is calculated with the adiabatic sound speed instead of the isothermal sound speed (h_p in Equation (2) is thus ${c}_{s,\mathrm{adi}}/R{\rm{\Omega }}=\sqrt{\gamma }{c}_{s,\mathrm{iso}}/R{\rm{\Omega }}$). It is a little bit surprising that the spiral shape in disks with T_cool = 1 follows the spiral shape in adiabatic disks. We think this is due to the short timescale for the flow to move across the spiral wake. The spiral wake has a typical width smaller than the disk scale height, while the background flow moves across the wake at nearly Keplerian speed. Thus, when the fluid travels in and out of the wave/shock, its response time is much smaller than the orbital time and thus behave adiabatically.

Zoom In Zoom Out Reset image size

Other aspects of the spiral shocks in non-isothermal cases are similar to the isothermal cases, e.g., a higher mass planet excites a more open spiral shock, and the inner spiral shocks become slightly more open at higher altitudes.

To illustrate the shape of the spiral shocks in the physical space, we plot the relative density perturbation in Cartesian coordinates in Figure 4. Clearly, the more massive the planet is, the more the spiral shocks deviate from linear theory. Two well separated inner arms are also apparent when the planet mass is large, and the separation between these two arms is larger when the planet is more massive (comparing the middle and bottom panels in Figure 4).

Zoom In Zoom Out Reset image size

To see how successful the weakly nonlinear density wave theory of Goodman & Rafikov (2001) and Rafikov (2002) fits the shape of the shocks in simulations, we plot in Figure 5 the density contour and density profiles along the azimuthal direction at R = 0.3, 0.5, 1.5, and 2 for run CM2ISO. In the density contour panel, we can clearly see that the secondary arm appears at the edge of the low density rarefaction wave region of the primary arm. After the secondary arm propagates for a short distance, a tertiary arm appears close to the low density rarefaction wave region of the secondary arm. In the panels showing the density profiles, we follow Goodman & Rafikov (2001) and Rafikov (2002) to shift the density profiles so that ϕ' = 0 corresponds to the wake position from linear theory (black curve in Figure 5 or Equation (2)). The density profiles clearly show that the shocks are N-shaped. A rarefaction wave follows the shock front, and δρ at the rarefaction wave region can be negative before it merges to the background flow. The shock fronts deviate from ϕ' = 0, and due to the shock expansion, the deviation is larger when the shock is further away from the planet. Under the shearing-sheet approximation, the amplitude and width of the N-shaped shock scale as $| R-{R}_{p}{| }^{-5/4}$ and $| R-{R}_{p}{| }^{5/4}$ at $| R-{R}_{p}| \gg 0$ based on the weakly nonlinear density wave theory of Goodman & Rafikov (2001). In a global disk spanning a large range of radii, the amplitude and width of the N-shaped shock scale as t^−1/2 and t^1/2, where t is given in Equation (43) of Rafikov (2002). For our disk parameters, we have

$$\begin{eqnarray}&&t\propto |{\displaystyle \int }_{1}^{R/{R}_{p}}{| {s}^{3/2}-1| }^{3/2}{s}^{-13/8}{ds}|.\end{eqnarray} \tag{ 17 }$$

Since the middle point of the N-shaped shock is at ϕ' ∼ 0 (Rafikov 2002), we expect that in our global simulations the azimuthal distance between the shock front and the path predicted by Equation (2) (ϕ' = 0) should also scales as t^1/2.

Zoom In Zoom Out Reset image size

To test this prediction, we have measured the shock positions at R = 0.5 and 1.5 in Figure 5, which are ϕ' = −0.9 and 0.43 respectively. These two positions are labeled as the dashed lines in R = 0.5 and 1.5 panels. Then we calculate the shock positions at R = 0.3 and 2 to be −1.7 and 0.9, using their positions at R = 0.5 and 1.5 together with the scaling relationship t^1/2 where t at R = 0.3, 0.5, 1.5, and 2 are calculated from Equation (17). These predicted shock positions are labeled as the dashed lines in R = 0.3 and 2 panels. We can see that they agree with the actual shock positions in the simulation very well. This confirms that the shock positions are determined by the nonlinear expansion of spiral shocks. Using the same approach, we have calculated the nonlinear shock position at every radius for R < 1 and R > 1 with the normalization based on shock positions at R = 0.5 and 1.5. This new predicted shock shape is plotted as the dotted curve in the left panel of Figure 5. Despite some offset at R close to the planet, which is expected since the simple relationship t^−1/2 holds only when $| R-{R}_{p}| \gg 0,$ a good agreement has been achieved between the nonlinear density wave theory and the simulations. Note that in this comparison, we did not calculate the shock strength directly from nonlinear theory and compare its amplitude with simulations, instead we use the scaling relationship to verify the propagation of the shock. In future, direct comparison is desired when we have a more complete nonlinear theory which can calculate the shock excitation directly.

Although the primary arm can be fitted by the weakly nonlinear density wave theory, the excitation of the secondary (or even tertiary) arm still lacks a good theoretical explanation. It may be related to the low m mode (e.g., m = 2, 3, similar to disks in binary systems) or some nonlinear wave coupling. Figure 5 suggests that a secondary spiral arm is excited at the other end of the N-shaped primary shock (the R = 0.5 panel). After it is excited, it steepens to shocks and becomes another N-shaped shock later (the R = 0.3 panel). Its shock front can even travel into the rarefaction wave of the primary arm (e.g., in the R = 0.3 panel, the secondary shock is almost at ϕ' = 0 where the rarefaction wave of the primary arm should reside.). Unlike the primary arm which already dissipates when it travels inward from R = 0.5 to 0.3, this secondary arm is excited later and becomes stronger from R = 0.5 to 0.3. At R = 0.3, the secondary arm is even stronger than the primary arm. By comparing R = 0.5 and R = 0.3 panels, we also notice that the secondary arm almost keeps the same azimuthal separation with the primary arm (Δϕ' ∼ 1.7) during its propagation. Furthermore, at R = 0.3, a tertiary arm starts to appear at the other end of the N-shaped secondary arm.

To summarize the results in this section, Figure 6 shows the pitch angle from the linear theory and those measured in numerical simulations. When the planet mass is low (e.g., CM1ISO, the green dots), the measured pitch angle agrees with the linear theory. When the planet mass increases, the pitch angle also increases. For the 6 M_J case, the measured pitch angle of the spiral wake in the h_p = 0.1 disk is close to the pitch angle predicted in a much thicker disk (h_p = 0.2) using the linear theory. The secondary and even the tertiary arms have similar pitch angles as the primary arms. Thus, if we know the disk thermal structure very well, we can use the deviation of the measured pitch angle from the linear theory to estimate the embedded planet mass.

Zoom In Zoom Out Reset image size

Since near-IR scattered light observations only probe the shape of the disk surface, 3D structure of spiral shocks can affect the observational signatures of these spiral shocks. Intuitively, we would expect that the spiral shocks have complicated 3D structure. First, the wave excitation must have 3D structure since, at the same R in the disk, the distance between the planet and the disk surface is larger than the distance at the midplane, and the force is thus weaker at the disk surface. Second, the wave propagation may have 3D structure considering the disk becomes thinner at smaller R. Waves/shocks are more converged when they propagate inwards. They can also channel to the disk surface (Lubow & Ogilvie 1998; Bate et al. 2002; Lee & Gu 2015), and, during their propagation from the high density region (e.g., the midplane) to the low density region (e.g., the disk surface), the amplitudes of perturbations have to increase to conserve the wave action. Since the amplitudes of perturbations determine when the waves will break into shocks, the dissipation can also be quite different between the surface and the midplane. All these effects can contribute to the 3D structure of spiral waves/shocks.

Due to these complicated effects, it is difficult to develop an analytic theory to study the planet-induced 3D shock structure, and we rely on numerical simulations to study such structure. However, before delving into the highly nonlinear shock regime, we can use the linear theory developed in Tanaka et al. (2002) to estimate the 3D effect of the density waves. Following their theory for locally isothermal disks, the structure of the waves in the z direction can be studied with Hermite polynomials (H_n(Z)). We first expand perturbed quantities (η) from our simulations into Fourier series

$$\begin{eqnarray}&&\eta =\displaystyle \sum _{m}{Re}\left[{\eta }_{m}{e}^{{im}(\phi -{{\rm{\Omega }}}_{p}t)}\right],\end{eqnarray} \tag{ 18 }$$

where the Fourier components η_m are complex functions of R and z. Then, η_m can be further expanded with Hermite polynomials in the z direction,

$$\begin{eqnarray}&&\eta =\displaystyle \sum _{m=0}^{\infty }\;\displaystyle \sum _{n=0}^{\infty }{Re}\left[{\eta }_{m,n}{H}_{n}(Z){e}^{{im}(\phi -{{\rm{\Omega }}}_{p}t)}\right],\end{eqnarray} \tag{ 19 }$$

where Z is the normalized height as Z = z/H(R), and the first three Hermite polynomials are

$$\begin{eqnarray}&&{H}_{0}(Z)=1,\quad {H}_{1}(Z)=Z,\quad {H}_{2}(Z)={Z}^{2}-1.\end{eqnarray} \tag{ 20 }$$

By using the normal orthogonal relation between H_n, we have

$$\begin{eqnarray}&&{\eta }_{m,n}=\displaystyle \frac{1}{\sqrt{2\pi }n!}{\displaystyle \int }_{-\infty }^{\infty }{e}^{-{Z}^{2}/2}{H}_{n}(Z){\eta }_{m}{dZ}.\end{eqnarray} \tag{ 21 }$$

We can use η_m,n at different n to estimate the relative importance of different Hermite components. To compare with Figure 2 in Tanaka et al. (2002), we show m = 10, n = 0, 2, 4 Fourier–Hermite components for the perturbed density (δρ/ρ₀)⁷ in Figure 7. With the same parameters, the right panel of Figure 7 is very similar to Figure 2 in Tanaka et al. (2002). By comparing the left and right panels in Figure 7, we can see that although higher-order Hermite components have similar strength between thin and thick disks, n = 0 Hermite component is much weaker in a thick disk, suggesting that the wakes in a thick disk have more significant 3D structures than wakes in a thin disk.

Zoom In Zoom Out Reset image size

Figure 7 suggests that higher-order vertical components can dominate the disk structure at the atmosphere. Although it shows that the n = 4 component (η_10,4) is 10 times weaker than the n = 2 component (η_10,2), and the n = 2 component is 10 times weaker than the n = 0 component (which led Tanaka et al. 2002 to conclude that most of the angular momentum excited by the planet will be carried by 2D free waves), the base function (Hermite polynomials) at z = 3H gets ∼10 times larger sequentially from H₀ to H₂ and to H₄. Thus, ${\eta }_{\mathrm{10,2}}{H}_{2}$ and η_10,4H₄ are still comparable with η_10,0H₀. The density structure at the disk atmosphere can be significantly affected by higher-order vertical modes.

Although the modal analysis is useful to verify numerical simulations and can be suggestive on the relative amplitudes of various modes, it is the 3D structure in the real space that determines the observational signatures of waves/shocks.

By studying the shock structure in real space, we first find that the 3D shock structure is dramatically different between inner and outer arms. For the inner arms, the density perturbation is much larger at the disk surface than at the disk midplane. The 3D structure of the inner spiral arms at R = 0.5 is shown in Figure 8. At R = 0.5 and z = 3H (solid curves), the differences between the maximum and minimum density in the logarithmic scale are 0.015, 0.4, and 1.3 for M1, M2, and M3 cases respectively, in comparison with 0.004, 0.2, and 0.4 at the disk midplane (dotted curves). At the same radius (R = 0.5), the position of the wakes in non-isothermal disks (lower panels) are at smaller $\phi -{\phi }_{p}$ compared with those in isothermal disks (upper panels). This is because the wakes are more open in non-isothermal disks, as discussed in Section 4.

Zoom In Zoom Out Reset image size

The secondary inner spiral arms/shocks are also more prominent at the disk surface than at the disk midplane. At the disk midplane, the secondary arms have lower amplitudes compared with the primary arms (dotted curves in Figure 8), while at z = 3H, the secondary arms have almost the same amplitudes as the primary arms (solid curves). The large amplitude of primary and secondary inner arms at the disk surface is due to the corrugated motion in the v_θ direction. In Figure 8 which is in the corotating frame with the planet, the disk material flows in the direction from the left side to the right side of the figure. Before meeting with the shock, the disk is in vertical hydrostatic equilibrium with the background density and v_θ = 0 (Figure 9). After the shock, the disk material loses angular momentum and moves inwards with v_r < 0 (Figure 9). At the same time, v_θ also becomes negative, compressing the disk material at the midplane. This downward motion decreases the density of the rarefaction wave at z = 3H. Before meeting the secondary shock, v_θ starts to increase and becomes positive, leading to a higher density at the disk surface. At the secondary shock, v_θ reaches the maximum positive velocity and leads to the highest density at the disk surface for the secondary shock in Figure 8. This corrugated motion, first negative and then positive v_θ, leads to an enhanced contrast between the spiral shock and the rarefaction wave after the shock.

Zoom In Zoom Out Reset image size

On the other hand, the density perturbation of outer spiral arms is similar between the disk surface and the disk midplane, especially for isothermal disks, as shown in Figure 10. At R = 2, regardless of height, the differences between the maximum and minimum density in the logarithmic scale are both 0.002 for SM1ISO, 0.1 for SM2ISO, and 0.3 for SM3ISO. This lack of vertical variation is also reflected in Figure 11 where v_θ is very small compared with v_r (v_θ is almost two orders of magnitude smaller than v_r). Thus, the density structure of the outer spiral is mainly determined by v_r and v_ϕ in the horizontal plane in stead of the corrugated motion in the θ direction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For non-isothermal runs (bottom panels in Figure 10), the density perturbation of the outer spiral arms at the disk surface is slightly higher than the perturbation at the disk midplane. Disk material flows from the right-hand side of the figure to the left-hand side in Figures 10 and 11. When it meets the shock, it develops a v_θ toward the disk surface, which enhances the density at the disk surface. Although it is tempting to contribute such difference between isothermal and non-isothermal runs to the nonlinear hydraulic jumps (shock bores) (Boley & Durisen 2006), such disk structure also appears even in the linear regime for the 0.01 M_J case, implying that it may be a linear effect and related to the eigenfunctions of the 3D waves excited by the planet.

Finally, to illustrate the increase of the density perturbation with height in disks and the qualitative difference between inner and outer arms, we plot the relative density perturbation along the radius at different heights (z = 0, 1, 2, 3, 4 H) in Figures 12 and 13. The relative density perturbation is defined as ${\rho }_{\mathrm{max}}(R,z)/{\rho }_{\mathrm{min}}(R,z)-1,$ where ρ_max(R, z) and ρ_min(R, z) are the maximum and minimum density along the azimuthal direction (ϕ = [0, 2π]) at the fixed R and z. We can see that the inner arms and outer arms are qualitatively different. For inner arms, the relative density perturbation is getting larger at higher altitudes, and the perturbation can increase by more than a factor of 10 from the midplane to the disk surface. For the outer arms, the relative density perturbation is almost unchanged between the midplane and the disk surface in isothermal disks (Figure 12) and only increases slightly from the midplane to the disk surface in non-isothermal disks (Figure 13). Overall, for the inner arms, the large density perturbation at the disk surface has a significant effect on the near-IR observations as shown below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.1. Numerics, Longterm Evolution, and Different Disk Structures

As shown in Figure 1, the 3D spiral shocks are not perpendicular to the disk midplane. They curl in a way that they are almost along the θ direction in the spherical-polar grid. To quantify the curl, the upper panel of Figure 14 shows the relative density perturbation at the disk midplane (black curves) and at $\theta =\pi /2+0.35$ (cyan curves) for run SM2ISO (solid curves). Please note that the x-axis in Figure 14 is r (the radial position in spherical-polar coordinates) instead of R (the radial position in cylindrical coordinates). The fact that the shock density peaks at almost the same r at either the disk midplane or $\theta =\pi /2+0.35$ demonstrates that the shocks almost curl along the θ direction. Such coincidence makes us suspect that they could be numerical artifacts due to the adopted spherical-polar grid structure. Thus, we use Athena to carry out a similar simulation but using cylindrical coordinates (run SM2ISOcyl introduced in Section 3.3). After transforming simulation outputs from cylindrical coordinates to spherical-polar coordinates, the relative density perturbation for SM2ISOcyl is shown in the upper panel of Figure 14 as the dotted curves. The shocks are at similar positions as those in SM2ISO, which confirms that the curled shock structure is real instead of numerical artifacts.

Zoom In Zoom Out Reset image size

The physical reason for the curled shock structure is unclear. One would suspect that it may be related to the radial temperature gradient in the local isothermal disk. In such disks, the disk rotates at different angular velocities at different disk heights. Such vertical shear can change disk dynamics, such as leading to the vertical shear instability (Nelson et al. 2013) and it may also lead to curled shock structures. To test this idea, we have run one simulation with a constant temperature in the whole disk (SM2ISOconstT). In such a disk, the disk rotates at a constant angular velocity at a given R independent on the disk height. The relative density perturbation is shown in the bottom panel of Figure 14. However, the shock density still peaks at almost the same r at different θ, implying that the curl of the shock is not caused by the vertical shear in disks. Thus, the physical mechanism for the curled shock remains unclear and still deserves further explore.

To demonstrate that our derived shock structure is independent on specific numerical choices in our model, we show the density perturbation for run SM2ISOlong in Figure 15. Run SM2ISOlong is different from SM2ISO in several aspects (as introduced in Section 3.3): (1) the simulation runs for 120 orbits instead of 10 orbits; (2) it has a wave damping region, (3) the planet mass is ramped up slowly over 10 orbits, and (4) it uses zero viscosity. Despite these differences, its density perturbation at 15 orbits (bottom panel of Figure 15) is quite similar to the density perturbation in SM2ISO (middle panel of Figure 12). When a gap is induced at 120 orbits, the spiral shocks become weaker at the disk region close to the gap, while it still maintains the same strength at the region far away from the gap. We notice that vortices start to develop at the outer gap edge (R ∼ 1.4) at 120 orbits, which slightly enhances the density perturbation at the outer gap edge.

Zoom In Zoom Out Reset image size

We have also studied how the shock structure can be affected by different disk structures. Figure 16 shows the relative density perturbation of the wakes in a thinner disk (STHIN). Compared with Figure 12, we can see that the density perturbation increases by a factor of 2 from the midplane to 3 H, and a factor of 4 from the midplane to 4 H in the thin disk, compared with a factor of 4 and 10 respectively in the thick disk. This suggests that the wakes/shocks have more significant 3D structure in a thicker disk. Although this does not favor detecting spiral arms in thinner disks in future, we need to keep in mind that the thermal mass (Equation (4)) is smaller in a thinner disk so that the same mass planet corresponds to a more massive planet in the scale of thermal mass and it will excite stronger density waves. In this case the inner arms may still be observable in a thin disk.

Zoom In Zoom Out Reset image size

6.2. Near-IR Images

To understand how the 3D structure of density waves/shocks affects observations, we post-process our hydrodynamical simulations with MCRT calculations to generate near-IR scattered light images (Figures 17–19). The details on the MCRT calculations are presented in Dong et al. (2015a, 2015b). To assign physical scales to our simulations, we assume that the planet is at 50 AU and the central source is a typical Herbig Ae/Be star (2 M_⊙) with a temperature of 10⁴ K and a radius of 2R_⊙. ISM dust grains have been used and their distribution is assumed to follow the gas distribution. Following Dong et al. (2015a), we assume that the gas disk is 0.02 M_⊙. With the gas to dust mass ratio of 100:1, the total dust mass is 2 × 10⁻⁴M_⊙. Further assuming 10% of dust is in the form of ISM dust, the total mass of the ISM dust is 2 × 10⁻⁵M_⊙. In this model, Toomre Q parameter at 50 AU is ∼30, so neglecting disk self-gravity is justified. In MCRT simulations, photons from the central star are absorbed/reemitted or scattered by the dust in the surrounding disk. The convolved images in Figures 17–19 are derived by convolving full resolution images with a Gaussian point-spread function having a full width half maximum (FWHM) of 006. This resolution is comparable with NIR direct imaging observations using Subaru, VLT, and Gemini. In the right two panels of Figures 18 and 19, we assume that the object is 140 pc away, while in Figure 17 and the left two panels of Figures 18 and 19 we assume that the distance is 70 pc so that the inner arms are shown more clearly.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Both full intensity and polarized intensity images are calculated and shown in Figure 17. The spiral arms are evident in these images. When the disk is viewed face-on, both full intensity and polarized intensity images are almost identical except that the full intensity image is almost a factor of 2 brighter than the polarized intensity image. The polarization fraction, which is defined as the ratio between the polarized intensity and full intensity, is almost a constant (∼0.45). When the disk is viewed at some inclination angle, the polarization fraction is not a constant due to the dust forward scattering. The dependence of the polarization fraction on disk inclination can be useful to determine the disk inclination angle.

To highlight the importance of the 3D wave structure, we have also computed models by only using the disk midplane density from simulations, which is labeled as $2{\rm{D}}\to 3{\rm{D}}$ in Figures 18 and 19. In these models, we assume that the disk is in vertical hydrostatic equilibrium and puff up the midplane density to higher altitudes as $\rho ={\rho }_{\mathrm{mid}}{\rm{exp}}({-z}^{2}/2{h}^{2})$ where $h/R=0.1{R}^{0.25}$ (the same as the scale height used in 3D simulations).

Figures 18 and 19 show that, in 3D models, the inner arms are considerably more prominent than the outer arms, and normally a secondary arm can be as bright as the primary arm. The shape of the inner arms clearly deviate from the prediction of linear theory. As expected from the nonlinear expansion of spiral shocks, the pitch angle of the inner spiral arms in the more massive planet case (SM3ISO, Figure 19) is larger than those in the less massive planet case (SM2ISO, Figure 18). The inner arms are also quite sharp, while the outer arms are quite broad and sometimes indistinguishable from the background disk. This difference is partly because the sharp shock fronts are facing the star for the inner arms, while they are facing away from the star for the outer arms. The different geometry at the disk surface can greatly affect the intensity of the scattered light images (Takami et al. 2014). We calculate an approximate scattering surface, defined as the disk surface where the column density is 0.01 (in code units), for the SM3ISO model at $\phi -{\phi }_{p}=78^\circ ,$ shown in Figure 20. Clearly, for the inner arms, the shock fronts are facing the star, while, for the outer arms, the smooth rarefaction waves are facing the star. Since the rarefaction waves change gradually with radius, they are illuminated by the star more uniformly than the shock. Thus, the outer arms appear quite broad. However, when the planet mass is not very high (1 M_J case), the width of rarefaction waves in the radial direction can be smaller than the size of the observational beam, and we will not be able to distinguish the inner and outer arms based on the sharpness of the arms.

Zoom In Zoom Out Reset image size

In 3D models of Figures 18 and 19, the secondary inner arms are as apparent as the primary arms, even though the primary arms have higher surface density than secondary arms. This is due to the corrugated motion discussed above in Figures 8 and 9, which increases the density of the secondary arms at the disk surface. The secondary arm is offset from the primary arm with some azimuthal angle, as also shown in the surface density plot (Figure 2). This offset is smaller in SM1ISO (Figure 18) than that in SM6ISO (Figure 19). The two spiral arms are ∼100° apart in SM1ISO (Figure 18), while almost 180° apart in the more massive planet case (SM6ISO, Figure 19).

By comparing 3D models with $2{\rm{D}}\to 3{\rm{D}}$ models in Figures 18 and 19, we find that the inner spiral arms are more prominent in 3D models, as expected due to inner arms' higher density perturbation at the disk surface in 3D models (Figure 12). Even in convolved images, the polarized intensity of inner arms is at least twice stronger in 3D models than in $2{\rm{D}}\to 3{\rm{D}}$ models.

The outer shocks are not very apparent, and they are similar between 3D and $2{\rm{D}}\to 3{\rm{D}}$ models since the density perturbation of outer arms is almost height independent (Section 5). As discussed in Section 4, the outer arms coincidently follow linear theory (Equation (2)). When the planet mass is very large as in Figure 19, the secondary outer arm starts to become visible.

Another noticeable difference between 3D and $2{\rm{D}}\to 3{\rm{D}}$ models is that the planet (or the circumplanetary region) is bright in $2{\rm{D}}\to 3{\rm{D}}$ models while it is dim in 3D models. This is because, when we puff the disk from 2D to 3D, we have ignored the planet's gravity so that the higher density in the circumplanetary region leads to the higher density at the disk surface. The circumplanetary region even casts a shadow to the outer disk in $2{\rm{D}}\to 3{\rm{D}}$ models. In realistic 3D models, the gravity of the planet has been self-consistently included, which pulls the circumplanetary material toward the disk midplane and leads to a lower density at the disk surface. Thus, the circumplanetary region receives less irradiation by the central star, and becomes dark. On the other hand, the outer disk beyond the planet is better illuminated and thus becomes bright instead of being shadowed.

However, we need to keep in mind that we have ignored the luminosity from the planet and the circumplanetary disk in our models. Zhu (2015) point out that accreting circumplanetary disks can be very bright (∼0.001 L_⊙ if the circumplanetary disk accepts onto Jupiter at a rate ∼10⁻⁸M_⊙ yr⁻¹). Such high luminosity may be able to illuminate the circumplanetary region significantly. We may also be able to directly detect such accreting circumplanetary disks in direct imaging observations operating at mid-IR wavelengths.

6.3. Observational Implications

We have carried out 2D and 3D hydrodynamical simulations to study spiral wakes/shocks excited by young planets. Simulations with different planet masses (0.01, 1, and 6 M_J) and different equations of state (isothermal and adiabatic) have been carried out.

• 1.  
  We find that the linear density wave theory can only explain the shape of the spiral wakes excited by a very low mass planet (e.g., 0.01 M_J). Spiral shocks excited by high mass planets clearly deviate from the prediction of linear theory. For a more massive planet, the deviation is more significant and the pitch angle of the spiral arms becomes larger. This phenomenon can be nicely explained by the wake broadening from the nonlinear density wave theory (Goodman & Rafikov 2001; Rafikov 2002). A more massive planet excites a stronger shock which expands more quickly, leading to a larger pitch angle.
• 2.  
  A secondary (or even tertiary) inner spiral arm is also excited by the planet. It seems to be excited at the edge of the N-shaped primary arm. The more massive is the planet, the larger is the separation between the primary and secondary arm. At the disk surface, the secondary inner arm can be as strong as the primary arm. The secondary inner arm almost keeps the same azimuthal separation with the primary arm at every radius in the disk. The excitation mechanism for the secondary (tertiary) arm deserves further study.
• 3.  
  The spiral shocks have significant 3D structure. They are not perpendicular to the disk midplane. They are curled toward the star at the disk surface. This further increases the pitch angle of the inner arms at the disk surface, but reduces the pitch angle of the outer arms at the disk surface. For outer arms, this effect compensates the increased pitch angle due to wake broadening. Eventually at the disk surface, the shape of outer spiral arms still roughly follows the prediction of linear theory, while the inner arms are considerably more opened than predicted by the linear theory.
• 4.  
  The inner spiral shocks also have significant vertical motion. The corrugated motion increases the density perturbation of the inner spiral arms by more than a factor of 10 at z ∼ 3–4H compared with the perturbation at the disk midplane. This can dramatically increase the contrast of the spiral patterns in near-IR scattered light images. The outer spiral shocks have little vertical motion in isothermal disks. With a non-isothermal EOS, there are some vertical motions for the outer arms, which can make the outer arms more apparent.
• 5.  
  We have combined our hydrodynamical simulations with MCRT calculations to generate near-IR scattered light images. We find that the inner spiral arms are prominent features that are observable by current near-IR imaging facilities. Besides the apparent spiral patterns in some transitional disks, planet-induced spiral arms may also be responsible to some marginally detected near-IR features within cavity of some transitional disks. We further demonstrate that inner spiral arms in synthetic near-IR images using full 3D hydrodynamical models are much more prominent than those based on 2D models assuming hydrostatic equilibrium, consistent with the 3D structure of the inner arms. On the other hand, the outer shocks are not very apparent and they are similar in synthetic images using 3D and 2D models. This indicates the need to model observations (especially for inner arms) with full 3D hydrodynamics.
• 6.  
  The different geometry between the inner and outer arms also affects their appearance in near-IR images. The sharp shock fronts of the inner arms face the central star directly, producing sharp narrow spiral features in observations. On the other hand, for the outer arms, the smooth rarefaction waves face the central star, producing broad and dimmer spiral features.
• 7.  
  In near-IR images, the circumplanetary region is very dim since the planetary gravity reduces the density at the disk atmosphere. However, the disk region behind the planet can be better illuminated and becomes bright.

In the Appendix, we have shown that buoyancy resonances are confirmed in global adiabatic simulations even if the disk has a moderate cooling rate. They can lead to sharp density ridges around the planet, which may have observational signatures.

Overall, spiral arms (especially inner arms) excited by low mass companions are prominent features in near-IR scattered light images. Most importantly, we have three independent ways to infer the companion's mass: (1) the pitch angle of the spiral patterns, (2) the separation between the primary and secondary arms, and (3) the contrast of the spiral patterns.

In a companion paper (Dong et al. 2015a), we have combined MCRT and hydrodynamical simulations from this paper, and shown that planet-induced inner arms can explain spiral patterns revealed by recent near-IR direct imaging observations for SAO 206462 and MWC 758. We want to caution that our proposed model can not explain the gaps at the inner disks. Other mechanisms (e.g., another planet or disk photoevaporation) have to be invoked to explain the inner gaps.

We thank the anonymous referee for a thorough report that improved the initial manuscript. All hydrodynamical simulations are carried out using computer supported by the Princeton Institute of Computational Science and Engineering, and the Texas Advanced Computing Center (TACC) at The University of Texas at Austin through XSEDE grant TG- AST130002. This project is supported by NASA through Hubble Fellowship grants HST-HF-51333.01-A (Z.Z.) and HST-HF-51320.01-A (R.D.) awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555.

The inner and outer spiral shocks are due to the steepening of spiral density waves which are excited by the planet at Lindblad resonances. At Lindblad resonances, the Doppler shifted frequency matches the disk epicyclic frequency ($m({{\rm{\Omega }}}_{p}-{\rm{\Omega }})=\pm \kappa$) and density waves are excited.

However, besides the epicyclic frequency, the disk also has other natural frequencies. When the disk is not strictly isothermal, it has a non-zero Brunt–Väisälä frequency

$$\begin{eqnarray}&&N(z)=\sqrt{\displaystyle \frac{\gamma -1}{\gamma }}\displaystyle \frac{g(z)}{{c}_{s,\mathrm{iso}}}\end{eqnarray} \tag{ 22 }$$

where ${c}_{s,\mathrm{iso}}^{2}=p/\rho .$ Matching the Brunt–Väisälä frequency with the Doppler shifted frequency, we have

$$\begin{eqnarray}&&\sqrt{\displaystyle \frac{\gamma -1}{\gamma }}\displaystyle \frac{{{\rm{\Omega }}}_{K}(R)z}{H}{\left(1+\displaystyle \frac{{z}^{2}}{{R}^{2}}\right)}^{-3/2}=\pm m({{\rm{\Omega }}}_{p}-{\rm{\Omega }}).\end{eqnarray} \tag{ 23 }$$

Given a m, Equation (23) gives the position of the resonances. These buoyancy resonances were discovered in shearing box simulations (Zhu et al. 2012) and studied analytically in Lubow & Zhu (2014). They have significant contributions to the planetary torque, especially around the planet, which may affect planet migration. These resonances are infinitely thin, and no waves are excited to carry the deposited angular momentum and energy away. Their dissipation relies on microscopic viscosity or radiative cooling. Thin density ridges with large temperature and velocity variations appear at these resonances.

When various m modes overlap with each other, we can roughly estimate the position of the final density ridges caused by buoyancy resonances following Equations (10) and (11) in Zhu et al. (2012). First, given an m, we can calculate the corresponding resonance position at R and z. Then using the azimuthal wavelength for this mode λ = 2π/m, the geometric location of the constant phase 2nπ (n is integer) is given by ϕ = nλ (assuming the phase of buoyancy waves is 0 at the planet position). Thus,

$$\begin{eqnarray}&&\phi =\pm 2n\pi ({{\rm{\Omega }}}_{p}-{\rm{\Omega }})\sqrt{\displaystyle \frac{\gamma }{\gamma -1}}\displaystyle \frac{H}{{{\rm{\Omega }}}_{K}(R)z}{\left(1+\displaystyle \frac{{z}^{2}}{{R}^{2}}\right)}^{3/2}.\end{eqnarray} \tag{ 24 }$$

We plot temperature fluctuations and v_θ at z = 2H for SM1T1 in Figure 21. The positions of buoyancy resonances given by Equation (24) are plotted with dotted lines/curves. Figure 21 shows that both temperature fluctuations and v_θ are nicely tracked by Equation (24), suggesting that buoyancy resonances exist in disks even with T_cool = 1.

Zoom In Zoom Out Reset image size

Even though buoyancy resonances can affect planet migration, they may not be observed through direct imaging technique since the density fluctuations caused by these resonances are much weaker than the spiral density waves excited by Lindblad resonances. For example, in the z = 3H panels of Figure 3, we can see some density fluctuations close to the corotation region (especially for SM2T1, the 1 M_J case), but they are much weaker than the spiral shocks.

These buoyancy resonances also have vertical structure, as shown in Figure 22. We sliced through the disk at a fixed ϕ, and the buoyancy resonance curves are plotted as dotted curves in Figure 22. At the disk midplane, N = 0 and there are no buoyancy resonances. Since N increases with disk height, $| {{\rm{\Omega }}}_{p}-{\rm{\Omega }}|$ also needs to increase with height to match the Brunt–Väisälä frequency considering λ and m are the same at the same ϕ. Thus, the curves move away from the planet position toward the disk atmosphere.

Zoom In Zoom Out Reset image size

Figure 22 also shows the inner and outer spiral shocks which are hotter at the shock position.