SPIRAL WAVES TRIGGERED BY SHADOWS IN TRANSITION DISKS

The simulations were run in code units, assuming the central star mass, M₀ = M_⋆, and r₀ = 1 au, as the units of mass and length. The code unit of time, t₀, corresponds to the orbital period at r₀ divided by 2π, that is ${t}_{0}={({{GM}}_{\star }/{r}_{0}^{3})}^{-1/2}$. The gravitational constant is G = 1 in code units. The temperature unit is GM_⋆μm_p/(k_Br₀), with μ being the mean molecular weight (μ = 2.35 in all our simulations), m_p as the proton mass, and k_B as the Boltzmann constant. We adopt two values for the disk-to-star mass ratio, q = 0.25 and 0.05.

The computational domain in physical units extends from r = 10–150 au over n_r = 400 logarithmically spaced radial cells. The grid samples 2π in azimuth with n_ϕ = 800 equally spaced sectors. Gas material is allowed to outflow at the disk edges.

The initial density profile scales with r⁻¹:

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

where the value of Σ₀ is set by fixing the disk mass to 0.05 and 0.25 M_⋆ $\left({M}_{{\rm{disk}}}={\int }_{{r}_{{\rm{min}}}}^{{r}_{{\rm{max}}}}2\pi {\rm{\Sigma }}(r){rdr}\right)$. The disk initial aspect ratio, H/r, is fixed to 0.05. All models include self-gravity. Our choice of parameters is similar to other numerical models of transition disks (e.g., Dipierro et al. 2015).

The stellar irradiation heating per unit area is given by (Fröhlich 2003):

$$\begin{eqnarray}&&{Q}_{\star }^{+}=(1-\beta )\displaystyle \frac{{L}_{\star }}{4\pi {r}^{2}}\mathrm{cos}\phi ,\end{eqnarray} \tag{ 2 }$$

where β = 0.1 is a reflection factor (albedo) and ϕ is the angle formed between the incident radiation and the normal to the surface, given by $\mathrm{cos}\phi \simeq \tfrac{{dH}}{{dr}}-\tfrac{H}{r}$. The disk scale height is assumed to be in hydrostatic equilibrium, i.e., H = r c_s(T)/v_k, with c_s(T) as the sound speed and v_k as the Keplerian velocity. L_⋆ is the stellar luminosity.

In our simulations, shadows are cast by an inner region (inside the computational domain), which blocks a fraction of the stellar irradiation ${Q}_{\star }^{+}$ within an angle δ. The irradiation heating per unit area, including shadows projected onto the disk, ${Q}_{d}^{+}(r,\phi )$, reads;

$$\begin{eqnarray}&&{Q}_{d}^{+}(r,\phi )=\left\{\begin{array}{l}F(t)f(\phi ){Q}_{\star }^{+}\\ \quad {\rm{}}\;{\rm{if}}\;{\rm{}}| \phi | \gt (\pi -\delta /2)\wedge | \phi | \lt \delta /2\\ {Q}_{\star }^{+}\quad \quad {\rm{otherwise}},\end{array}\right.\end{eqnarray} \tag{ 3 }$$

where illuminated regions are connected by a smooth azimuthal function f(ϕ) = 1 − A exp(−ϕ⁴/σ²) − A exp(−(ϕ − π)⁴/σ²), with A = 0.999 and σ = δ/10. The time-dependent function F(t) is a shadow-tapering factor which gradually enables the shadows over a timescale of 30 orbits.⁵ Figure 1 shows the stellar irradiation prescription described by Equation (3).

Zoom In Zoom Out Reset image size

The equation for the thermal energy density, e, reads as (e.g., D'Angelo et al. 2003)

$$\begin{eqnarray}&&\displaystyle \frac{\partial e}{\partial t}+{\boldsymbol{\nabla }}\cdot (e{\boldsymbol{v}})=-P{\boldsymbol{\nabla }}\cdot {\boldsymbol{v}}+{Q}_{d}^{+}-{Q}^{-},\end{eqnarray} \tag{ 4 }$$

where ${\boldsymbol{v}}$ is the gas velocity, P is the pressure, and ${Q}_{d}^{+}$ is the stellar heating rate per unit area described by Equation (3). We implemented a radiative cooling per unit area function, Q⁻ = 2σT⁴/τ, where σ is the Stefan–Boltzmann constant and τ is the optical depth given by $\tau =\tfrac{1}{2}{\rm{\Sigma }}\kappa$, a constant opacity given by an effective κ = 130 cm² g⁻¹, obtained by taking into account an average opacity for a mix of dust species, informed by spectral energy distribution fitting of HD 142527 (Casassus et al. 2015). Dependency on opacity is further discussed in Section 3. For more details about this energy prescription, see Montesinos et al. (2015).

To close the system of equations, an ideal equation of state is used:

$$\begin{eqnarray}&&P={\rm{\Sigma }}T\bar{R},\end{eqnarray} \tag{ 5 }$$

where T is the mid-plane gas temperature and $\bar{R}={k}_{{\rm{B}}}/\mu {m}_{{\rm{p}}}$ is the gas constant. The thermal energy density is related to the temperature through

$$\begin{eqnarray}&&e={\rm{\Sigma }}T\left(\displaystyle \frac{\bar{R}}{\gamma -1}\right),\end{eqnarray} \tag{ 6 }$$

where the adiabatic index is fixed to the diatomic gas value γ = 1.4.

Control simulations without shadows were performed in order to test for azimuthal structures unrelated to illumination effects. Additionally, random noise at the ∼0.1% level was injected into the initial surface density of these simulations in order to test for stability and fragmentation.

Azimuthal features due solely to the disk self-gravity develop in control runs with M_d = 0.25 M_⋆ and ${L}_{\star }=1\;{L}_{\odot }$, on timescales of ≳10³ orbits. This is expected for such gravitationally unstable disks (Dipierro et al. 2015). For control models with M_d = 0.05 M_⋆ (gravitationally stable), no spiral structures emerge, independent of the strength of the stellar irradiation field.

When shadows are enabled, spiral-like structures emerge in the surface density of a simulated 0.05 M_⋆ disk irradiated by a $100\;{L}_{\odot }$ star. The onset of these structures occurs approximately after 130 orbital periods. On the other hand, no spirals appear if the 0.05 M_⋆ disk is illuminated by only 1 L_⊙.

We consider the onset of the shadow-induced spirals to be the time at which their scale heights suffice δH/H ≃ 0.3. This criterion is evaluated within the first 10–20 au of radius from the star. The perturbation on the scale height is calculated as $\delta H\equiv H-{H}_{0}$, where H₀ represents the background scale height without shadows. We note that a relative change of ∼0.2 in scale height is sufficient to produce detectable azimuthal signatures in scattered light predictions (Dong et al. 2015b; Juhász et al. 2015).

For the more massive 0.25 M_⋆ disk, with ${L}_{\star }=1\;{L}_{\odot }$, the first non-axisymetric structures appear approximately after 2500 orbits. For the same disk mass, but increasing stellar irradiation to 100 L_⊙, the first spirals arise shortly after only 150 orbits. These different timescales will be discussed in the next subsection.

Figure 2 shows the impact of different stellar irradiation values on the evolution of the density field, for a 0.25 M_⋆ disk. The top row shows the 100 L_⊙ model after 150, 250, and 400 orbits (from left to right, respectively), while the bottom panels show the 1 L_⊙ model at 2500, 3500, and 4000 orbits (from left to right). Each first snapshot corresponds to the moment spirals appear for the first time, as defined above. The rightmost top and bottom panels correspond to control runs, in which no shadows are present.

Zoom In Zoom Out Reset image size

In order to test the dependency of the disk evolution with the opacity, we also explore the case with κ = 1 cm² g⁻¹, finding the same results. This insensibility to the opacity, in order to produce shadow-induced spiral arms, is not surprising. The radiative cooling is computed according to T⁴/τ, therefore changes in κ imply (roughly speaking) changes in the temperature field of about ∼κ^1/4. Moreover, the amplitude of a thermal perturbation responsible for spiral development is given by δT/T (see the next section); therefore, in an equilibrium situation in which Q⁺ ∼ Q⁻, the assumption of a constant opacity results in a thermal amplitude perturbation independent of the opacity choice.

Our 2D simulations are not radiation hydrodynamics but just an ideal approximation. The role of the opacity will be considered in a future work, where we will study the illumination-induced spirals using full 3D radiation hydrodynamic simulations.

3.1.1. Pitch Angle

The spiral arms' pitch angles were computed by fitting the spirals with an Archimedean equation $r(\phi )={A}_{0}+{A}_{1}{\phi }^{n}$. The pitch angle Π is then obtained through $\mathrm{tan}{\rm{\Pi }}=(1/r){dr}/d\phi$. From this parameterized curve, a global pitch angle is calculated as the mean value along the curve.

Π varies over time for different model parameters. The shadow-induced spiral arms extend from A₀ (∼60–90 au, see Table 1) to the outer region of the disk ∼150 au (see Figure 2), as opposed to spirals resulting from GI, which persist over scales ≲100 au (e.g., Dong et al. 2015b). Table 1 summarizes the most important features of the shadow-induced spirals.

Table 1.  Main Spiral and Disk Characteristics

Model Parameters	Orbit Number	Spiral Parameters	Toomrea	Birth Timescale of
		A0, A1, n, Π°	Qmin	the Structuresb (in Orbits)
Md = 0.05 M⋆; L⋆ = 1 L⊙	3500	n/a	1.1	n/a
	4000	n/a	1.1
Md = 0.05 M⋆; L⋆ = 100 L⊙	250	68.7; 14.6; 2.1; 2259 ± 3.5	3.4	∼150
	500	56.9; 4.2; 1.64; 1382 ± 3.5	2.6
Md = 0.25 M⋆; L⋆ = 1 L⊙	3500	70.3; 4.6; 1.68; 1382 ± 3.5	0.6	∼2500
	4000	60.3; 5.4; 1.45; 1135 ± 3.5	0.5
Md = 0.25 M⋆; L⋆ = 100 L⊙	250	60.3; 28.2; 1.35; 1504 ± 2.5	1.0	∼150
	500	87.2; 17.8; 1.43; 1309 ± 2.0	1.2

Notes. The spirals can be described with the Archimedean function $r(\phi )={A}_{0}+{A}_{1}{\phi }^{n}$, where the pitch angle is given by $\mathrm{tan}{\rm{\Pi }}=(1/r){dr}/d\phi$. The value A₀ (in au) gives the inner radius of the spirals, which, in our models extends to the outskirts of the disk (∼150 au).

^aMinimum local value of the Toomre parameter at the given orbit number. ^bOrbit number when spirals are distinguishable for the first time in the density field (see Section 3.1).

Download table as:  ASCIITypeset image

Perturbations in a differentially rotating gaseous disk tend to wind up into spiral patterns (e.g., see the control model in the bottom rightmost panel in Figure 2). In our models, projected shadows act as forcing perturbations, which, efficiently and independently of the presence of random sources of symmetry breaking (e.g., gravitational instabilities), produce spirals. Because of the point-symmetric nature of this perturbation, density waves are excited with m = 2 morphology. Each shadow appears to trigger one spiral arm.

Shortly after the shadows have been cast onto the disk, and before the disk thermalizes, the gas pressure (Equation (5)) plummets at the shadowed regions. Figure 3 (left panel) shows the relative change in pressure $(P-{P}_{0})/{P}_{0}$, where P₀ is the background unshadowed gas pressure. This snapshot corresponds to a 100 L_⊙ and 0.25 M_⋆ model disk immediately after shadows are enabled. Figure 3 shows the azimuthal pressure acceleration, ${{\boldsymbol{a}}}_{\phi }$, after 30 orbits for two 0.25 M_⋆ disk models with L_⋆ = 1 L_⊙ and 100 L_⊙, middle and right panels, respectively. The pressure acceleration magnitude increases with stellar irradiation: the maximum value for the pressure acceleration increases from 0.2 to 24 cm s⁻² when the stellar irradiation passes from L_⋆ = 1 L_⊙ (middle) to L_⋆ = 100 L_⊙ (right panel).

Zoom In Zoom Out Reset image size

As the gas in the disk flows around the star (in a counterclockwise direction for all figures), it periodically enters and exits shadowed regions. At the interfaces between shadowed and fully illuminated sections, the gas is subjected to an azimuthal acceleration due to pressure gradients (i.e., ${{\boldsymbol{a}}}_{\phi }=-(1/r{\rm{\Sigma }}(r)){{\rm{\nabla }}}_{\phi }P$). As the gas enters a shadowed region (coming from high to low pressure), it experiences a positive azimuthal acceleration (anticlockwise). On the other hand, as it exits the shadow (moving from low to high pressure), the gas experiences a negative azimuthal acceleration (clockwise). The difference between the gain and loss of azimuthal acceleration is uneven, resulting in a net azimuthal acceleration in the counterclockwise direction (gas rotation direction), which varies with radius within the first ∼50 au. This creates a strong source of axisymmetry breaking, locally pushing the gas to move faster around the central star at the inner regions of the disk, piling up material at the shadow's frontier close to this inner region, which then leads to the formation of a trailing mode of spiral-like structures outspread by differential rotation.

Gas pressure acceleration depends on stellar luminosity as more luminous stars feed more energy into their disks, raising the temperature of the illuminated sections while leaving the shadowed region unaffected. This produces higher contrasts in the pressure field. This is the reason why spirals appear early (∼130 orbits) in models with L_⋆ = 100 L_⊙, when compared to 1 L_⊙ models in which the first spirals emerge only after ∼2500 orbits. The key factor is the contrast in pressure between dark and illuminated regions.

Depending on the gas cooling time, thermal perturbations can wear off over time. As the disk thermalizes, the illumination pattern becomes less effective at forcing a perturbation on the gas. In our simulations, this tendency to homogenize is seen after 10⁵ orbits. Ultimately, the induced spiral arms tend to blend and lose contrast with the background gas unless GI kick in. In these large viscous timescales, GI-induced spiral structures are expected to remain quasi-steady (e.g., Dipierro et al. 2015).

GI can be locally characterized by the Toomre parameter⁶ given by $Q=\tfrac{{c}_{{\rm{s}}}{\rm{\Omega }}}{\pi G{\rm{\Sigma }}}$, where c_s, Ω, and Σ correspond to sound speed, angular velocity, and surface density, respectively (Toomre 1964). The disk becomes unstable when Q < Q_crit, where Q_crit defines the range for which a disk is marginally unstable: 1 ≲ Q_crit ≲ 2.

According to the local Toomre parameter, denser and colder regions of the disk tend to have lower Q values. The ${M}_{{\rm{d}}}=0.05\;{M}_{\star }$ and L_⋆ = 100 L_⊙ model is stable everywhere in the disk, with a minimum Q value of 3.35. The most unstable model (M_d = 0.25 M_⋆ and L_⋆ = 1 L_⊙ attains a minimum local value of Q = 0.5. It is important to note that spiral structures emerge even if we disable self-gravity from our simulations. The shadow-induced spirals do not require the presence of GI for their development.

We input our simulations into a 3D radiative transfer code to produce scattered light predictions in the H-band (1.6 μm). We used the RADMC3D ⁷ code (version 0.39), assuming that the dust follows the same density field as the gas in the simulation, with a gas-to-dust ratio of 100. The dust distribution model consists of a mix of two common species: amorphous carbon and astronomical silicates. We used the Mie model (homogeneous spheres) to compute dust opacities for anisotropic scattering (Bohren & Huffman 1983). The optical constants for amorphous carbon (intrinsic grain densities 2 g cm⁻³) were taken from Li & Greenberg (1997), and those for silicates (intrinsic grain densities 4 g cm⁻³) were taken from Draine & Lee (1984).

To produce a three-dimensional (3D) volume, the vertical density structure was solved assuming hydrostatic equilibrium. The disk scale height is obtained from the temperature field (Section 2.2).

The temperature in the vertical direction is assumed to be constant and equal to the mid-plane temperature with T_gas = T_dust. The final image prediction is rendered using a second order volume ray-tracing.

Figure 4 shows a model prediction in the H-band based on a 0.25 M_⋆ disk with stellar irradiation of 1 L_⊙, after it had evolved for 3500 orbits. We recover the spiral arms in scattered light. We also obtain observable spiral structures in the L-band. The scattered light spirals follow the same equation as the one obtained from the density field, i.e., Π = 1382 (see Figure 2, bottom row, orbit 3500, and/or Table 1).

Zoom In Zoom Out Reset image size

It is worth mentioning that the assumption of hydrostatic equilibrium in 2D models likely underestimates the contrast of spirals images when compared with with full 3D hydrodynamical models (Zhu et al. 2015). In that case, the illumination effects impacting the disk should produce even more prominent scattered light images of spirals than those reported here.