SHOCK-DRIVEN ACCRETION IN CIRCUMPLANETARY DISKS: OBSERVABLES AND SATELLITE FORMATION

Planets form and grow in circumstellar disks. Initally when a protoplanet forms, the solid core is surrounded by a hydrostatic gaseous envelope that is in contact with the planet's Hill sphere and the rest of the circumstellar disk. As the core's mass increases, the planet envelope contracts, due to both the stronger gravity and cooling through radiation. Eventually, when the core mass reaches $\sim 10{M}_{\oplus }$, the planet undergoes runaway accretion and contracts significantly. At this stage, the envelope shrinks significantly and detaches from the planet's Hill sphere (Papaloizou & Nelson 2005). Material that resides beyond the Hill sphere can still flow into the Hill sphere, but it forms a circumplanetary disk (CPD) around the protoplanet to conserve angular momentum (Lubow et al. 1999; Ayliffe & Bate 2009). The accretion of the CPD onto the planet allows the continuous growth of the planet even after the planet's runaway growth stage.

CPDs may provide observational signatures of young forming planets in disks. Giant planets can be too faint to be detected, but accreting CPDs can be bright and detectable. For example, to form a giant planet of 1–10 M_J mass within the circumstellar disk's lifetime (about a few million years, Hernández et al. 2007), the CPD needs to accrete at a rate of $\dot{M}\gtrsim {10}^{-9}-{10}^{-8}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. Such an accretion disk will have a luminosity (e.g., Owen 2014; Zhu 2015) of

$$\begin{eqnarray}&&{L}_{\mathrm{disk}}=\displaystyle \frac{{{GM}}_{p}\dot{M}}{2{R}_{J}}=1.5\times {10}^{-3}{L}_{\odot }\displaystyle \frac{{M}_{p}}{1{M}_{J}}\displaystyle \frac{\dot{M}}{{10}^{-8}{M}_{\odot }\,{\mathrm{yr}}^{-1}},\end{eqnarray} \tag{ 1 }$$

which is as bright as a late M-type or early L-type brown dwarf and can be detected by current direct imaging techniques (Zhu 2015). Unlike a planet or brown dwarf that has an almost constant surface temperature, an accreting CPD has a lower temperature at larger disk radii, and those outer disk regions will emit significant infrared flux. Thus, the emission from an accreting CPD is redder than the emission from a planet or a brown dwarf. Direct imaging observations have found several red sources within circumstellar disks (Kraus & Ireland 2012; Quanz et al. 2013; Biller et al. 2014; Reggiani et al. 2014; Currie et al. 2015), and their photometry at near-IR bands is consistent with accreting CPDs (Zhu 2015). Emission lines for ${H}_{\alpha }$, which are another observational signature of accretion disks, are also found in some low-mass substellar objects (Zhou et al. 2014; Bowler et al. 2015). The most direct evidence for accreting CPDs is LkCa 15b, where two accretion tracers, the ${{\rm{H}}}_{\alpha }$ line and near-IR thermal emission, have both been detected (Sallum et al. 2015).

CPDs are also essential for satellite formation. In our solar system, most satellites around giant planets are in prograde, nearly circular, and coplanar orbits, implying that they formed in a shared CPD orbiting within the planet's equatorial plane. The ratio between the total mass of four major Galilean satellites and Jupiter's mass is $\sim 2\times {10}^{-4}$. The Saturnian satellite system also has a similar mass fraction with respect to Saturn's mass (Canup & Ward 2006). Assuming the gas-to-dust mass ratio is 100, a minimum gas mass of ∼0.02 planet mass is required to produce these satellites (Canup & Ward 2002, 2006). There are two main scenarios for Galilean satellite formations. One assumes that the satellites form in situ in a CPD containing this amount of material (e.g., Lunine & Stevenson 1982; Mosqueira & Estrada 2003). This disk is referred to as the "minimum-mass subnebula," which is the analog of "minimum-mass nebula" for the circumstellar disk. Spreading 0.02 Jupiter mass within 30 Jupiter radii (where the four major satellites reside), the minimum-mass subnebula has a very high gas surface density of $\sim {10}^{5}\,{\rm{g}}\,{\mathrm{cm}}^{-2}$. In the other scenario, the satellites form in a CPD with a much lower surface density, but the disk is dynamically evolving and being supplied by the circumstellar disk (Canup & Ward 2002, 2006). This "gas-starved" scenario only requires that the total supplied mass during the whole satellite formation timescale is larger than 0.02 planet mass.

The key question for understanding the structure of CPDs is how CPDs accrete. They may accrete in ways similar to how circumstellar disks accrete. Turner et al. (2014) suggest that the surface of CPDs can be ionized by X-rays so that the magnetorotational instability (MRI) can operate at the disk surface, leading to accretion, a process similar to the layered accretion proposed in protoplanetary disks (Gammie 1996). However, Fujii et al. (2011, 2014) find that the active layer is so thin (${\rm{\Sigma }}\sim {10}^{-3}\mbox{--}{10}^{-2}\,{\rm{g}}\,{\mathrm{cm}}^{-2}$) in CPDs that the accretion through MRI is negligible. On the other hand, other nonideal MHD effects (e.g., Hall effects; Kunz & Lesur 2013; Bai 2014; Lesur et al. 2014) that are important in protoplanetary disks can also be important in CPDs, and CPDs may accrete through magnetic braking (Keith & Wardle 2014). A magnetocentrifugal disk wind can also be launched in CPDs, carrying away angular momentum and leading to disk accretion (Quillen & Trilling 1998; Gressel et al. 2013). Almost all these proposed CPD accretion mechanisms (e.g., layered accretion, nonideal MHD effects, disk wind) can find their roots in circumstellar disk models. Thus, they are facing the same uncertainties as the accretion mechanisms in circumstellar disks: they sensitively depend on the net magnetic fields assumed and the detailed microphysics in the disk (e.g., dust size distribution).

On the other hand, CPDs are different from circumstellar disks in that they are subject to the tidal torque from the central star and are truncated within the Hill sphere of the planet (Martin & Lubow 2011a). In addition, circumstellar disk material flows through the Hill sphere, continuously replenishing CPDs. These properties make CPDs similar to disks with inflows from companion stars in close binary systems or cataclysmic variables (CV). In these binary systems, the tidal torque from the companion star will excite spiral density waves in disks, and when these waves shock in disks, they can transport angular momentum to the disk, leading to disk accretion. Previous inviscid isothermal simulations by Rivier et al. (2012) have suggested that accretion due to spiral shocks in CPDs is inefficient with Jupiter's mass doubling time ∼ 5 Myrs. Szulágyi et al. (2014) have measured a much higher accretion rate in their 3-D isothermal inviscid simulations (10⁻⁴ M_J yr⁻¹), but this measured accretion rate is not numerically converged with their higher resolution simulations. Instead, by measuring the torque exerted by the star in the simulations, they estimate that the real accretion rate is ∼ 2.5 × 10⁻⁶ M_J yr⁻¹. Recent studies on spiral shocks (Ju et al. 2016) have shown that the accretion due to spiral shock dissipation is sensitive to the disk thermodynamics assumed. When the disk is hot and the Mach number, defined as the ratio between the Keplerian speed and the sound speed, is small (<10), the equivalent α of shock-driven accretion can reach ∼ 0.01–0.02. The Mach number of CPDs depends on the disk accretion rate and is sensitive to the equation of state applied in simulations (D'Angelo et al. 2003; Ayliffe & Bate 2009; Machida 2009; Szulágyi et al. 2016).

Considering the potential importance of spiral shocks driving accretion in CPDs, in this paper we construct two-dimensional inviscid hydrodynamical simulations to study CPDs. These inviscid simulations differ from most previous works that use artificial viscosity (e.g., α viscosity) to sustain the disk accretion. Furthermore, to reduce the numerical viscosity in the CPD region, we adopt a grid structure centered on the planet. Since the thermodynamics that controls the disk Mach number is crucial for the shock-driven accretion, a simple radiative cooling scheme has been included in the simulations, which differs in this work from inviscid simulations by Rivier et al. (2012) and Szulágyi et al. (2014). We have also measured the mass accretion rate directly from simulations, unlike these two works, where they use the torque exerted by the star onto the CPD to estimate the disk accretion rate. As Ju et al. (2016) have pointed out, shock-driven accretion is determined by the shock dissipation (the difference between the tidal torque exerted on the disk and the angular momentum flux carried away by the wave), instead of the total torque alone. We will show that, by allowing the disk to be heated up by the shock, the inviscid simulation is numerically converged, and the disk can reach a steady state, transferring inflow material from the Hill radius all the way to the central planet quasi-steadily. This shock-driven accretion is very efficient, with $\alpha \sim 0.001\mbox{--}0.01$.

One caveat in our simulations is that our simulations are limited to two dimensions, and numerous previous three-dimensional (3D) simulations have shown that the infall from circumstellar to CPDs occurs at high altitudes (Bate et al. 2003; Machida et al. 2008; Tanigawa et al. 2012; Morbidelli et al. 2014; Szulágyi et al. 2014, 2016). Furthermore, since the CPD is significantly heated and puffed up, our simple cooling treatment based on the thin-disk approximation is inaccurate. Three-dimensional simulations with realistic radiative transfer (Szulágyi et al. 2016), thermodynamics, planet evolution (Ward & Canup 2010), and even planetesimal accretion (D'Angelo & Podolak 2015) are needed in the future to confirm if spiral shocks can lead to efficient accretion in CPDs.

In Section 2, our numerical method is introduced. The results are presented in Section 3, including both the disk structure and the accretion rate. After a short discussion in Section 4, the paper will be concluded in Section 5. The detailed energy budget of the accretion process is given in the Appendix.

We solve the Euler equations to study dynamics in CPDs using Athena++ (J. Stone et al. 2016, in preparation). Athena++ is a newly developed grid-based code using a higher-order Godunov scheme for MHD and constrained transport (CT) to conserve the divergence-free property for magnetic fields. Compared with its predecessor Athena (Gardiner & Stone 2005, 2008; Stone et al. 2008), Athena++ is highly optimized for speed 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 is crucial for the angular momentum budget analysis in Section 3.2.

In this paper, magnetic fields are ignored. We solve hydrodynamical equations in the cylindrical coordinate system (R and ϕ) that is centered on the planet. This grid choice is different from most previous works, which have adopted cylindrical or spherical coordinate systems that are centered on the star. Having the cylindrical coordinate system centered on the planet allows the material in CPDs to flow along the azimuthal grid direction, thus significantly reducing numerical truncation errors. Adopting the cylindrical coordinate system with uniform grids in ln(R) also allows each grid to have the same length along the azimuthal and radial directions throughout the whole simulation domain.

Furthermore, we adopt the rotating frame so that the star is stationary at R = 1, ϕ = 0. The equations being solved are

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial t}+{\rm{\nabla }}\cdot ({\rm{\Sigma }}{\boldsymbol{v}})=0,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Sigma }}{\boldsymbol{v}}}{\partial t}+{\rm{\nabla }}\cdot ({\rm{\Sigma }}{\boldsymbol{vv}}+P{\boldsymbol{I}})=-{\rm{\Sigma }}{\rm{\nabla }}{{\rm{\Phi }}}_{p}+{\boldsymbol{F}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial E}{\partial t}+{\rm{\nabla }}\cdot [(E+P){\boldsymbol{v}}]=-{\rm{\Sigma }}{\rm{\nabla }}{{\rm{\Phi }}}_{p}\cdot {\boldsymbol{v}}+{\boldsymbol{F}}\cdot {\boldsymbol{v}}-{Q}_{c}\end{eqnarray} \tag{ 4 }$$

, where Σ is the disk surface density, ${\rm{\Sigma }}{\boldsymbol{v}}$ is the linear momentum, ${{\rm{\Phi }}}_{p}$ = −GM_p/R is the potential of the planet having the mass of M_p, and $E=\epsilon +\tfrac{1}{2}{\rm{\Sigma }}{v}^{2}$ is the total energy density per unit area with being the internal energy. Note that since the planet is at the coordinate center, coordinates R and ϕ are the positions relative to the planet instead of the star. Here, ${\boldsymbol{F}}$ is the total external force excluding the gravitational force from the planet, and thus ${\boldsymbol{F}}$ includes the centrifugal force and Coriolis force due to the adoption of the rotating frame, the direct gravitational force from the star, and the indirect force (since the grid is centered at the planet instead of the center of mass). The detailed forms of these terms are given in Ju et al. (2016). Solving these equations in cylindrical coordinates also introduces geometrical source terms, which are written in forms conserving the angular momentum with machine-error precision (as in Ju et al. 2016). As will be shown in Section 3.2, maintaining angular momentum conservation is crucial for the angular momentum diagnostics.

No physical viscosity has been applied in our simulations. Thus, angular momentum and energy transport can only come from shock dissipation and numerical viscosity. Since simulations having twice the resolution produce almost identical results, which will be shown below, we conclude that numerical viscosity is negligible and that angular momentum and energy transport are dominated by the shock dissipation.

The equation of state for an ideal gas has been assumed so that $P=\epsilon (\gamma -1)$. We further assume γ = 7/5. The density-weighted disk temperature is $T=\mu P/({\mathfrak{R}}{\rm{\Sigma }})$, where ${\mathfrak{R}}$ is the gas constant and $\mu =2.4$. Since most disk mass is concentrated at the disk midplane, we use T to approximate the disk midplane temperature T_c.

The cooling rate per unit area, Q_c, is approximated by

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

where $\tau =({\rm{\Sigma }}/2){\kappa }_{R}({\rho }_{c},{T}_{c})$ is the optical depth to the disk midplane at radius R (Hubeny 1990). The midplane density ${\rho }_{c}$ is ${\rm{\Sigma }}/(\sqrt{2\pi }H)$, where H is the disk scale height ${c}_{s}({T}_{c})/{\rm{\Omega }}$. The Rosseland mean opacity ${\kappa }_{R}({\rho }_{c},{T}_{c})$ is assumed to be 10 ${\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$, which is the typical value (D'Alessio et al. 2001) at temperatures below the dust sublimation temperature of ∼1500 K.³ The particular form $\tau /(1+{\tau }^{2})$ is chosen so that the cooling term has the correct form in both optically thick and thin limits.⁴ The term $\sigma {T}_{\mathrm{ext}}^{4}$ represents the energy flux into the disk from the stellar irradiation, leading to a background disk temperature of ${T}_{\mathrm{ext}}$.

Due to the inclusion of the thermodynamics in the problem, we have to specify physical units in the simulations. We assume that the distance between the star and the planet is 5 AU, the star is a solar-mass star, and the planet's mass is 0.001 M_⊙ (M_J). And ${T}_{\mathrm{ext}}$ is 128 K so that the aspect ratio of the circumstellar disk at the planet position (${c}_{s}/{v}_{\phi }$ where ${v}_{\phi }$ is Jupiter's orbital velocity around the Sun) is 0.05. Besides the simulations with radiative cooling, we have also carried out simulations with isothermal and adiabatic equations of state for comparison.

Our cylindrical grids are uniformly spaced from 0 to $2\pi$ in the ϕ direction, and uniformly spaced in ln(R) from R = 2 or 4 Jupiter radii to half the distance between the planet and the star (which is 2.5 AU). Given an adopted resolution in the radial direction, we adjust the resolution in the azimuthal direction to make sure every grid cell is square in physical size. Since a Jupiter-mass planet should induce a gaseous gap in the circumstellar disk (Lin & Papaloizou 1986; Kley & Nelson 2012), we have manually carved out a gap in the simulations, as shown in Figure 1. Motivated by detailed gap-opening studies (e.g., de Val-Borro et al. 2006; Fung et al. 2014), we parameterize the gap density profile as

$$\begin{eqnarray}&&{\rm{\Sigma }}({\boldsymbol{r}})=\left\{\begin{array}{c}{{\rm{\Sigma }}}_{l}+\displaystyle \frac{{{\rm{\Sigma }}}_{h}-{{\rm{\Sigma }}}_{l}}{2}\times \left(2-\exp \left\{\displaystyle \frac{-{\rm{\Delta }}{R}_{{dp}}+{w}_{g}}{{w}_{t}}\right\}\right)\\ \quad \mathrm{if}\ {\rm{\Delta }}{R}_{{dp}}\gt {w}_{g}\\ {{\rm{\Sigma }}}_{l}+\displaystyle \frac{{{\rm{\Sigma }}}_{h}-{{\rm{\Sigma }}}_{l}}{2}\times \left(\exp \left\{\displaystyle \frac{{\rm{\Delta }}{R}_{{dp}}-{w}_{g}}{{w}_{t}}\right\}\right)\\ \quad \mathrm{if}\ {\rm{\Delta }}{R}_{{dp}}\lt {w}_{g}\end{array}\right.\end{eqnarray} \tag{ 6 }$$

where ${{\rm{\Sigma }}}_{l}$ and ${{\rm{\Sigma }}}_{h}$ are the disk surface density within and outside the gap, w_g is the gap half-width, w_t represents the sharpness of the gap, ${\rm{\Delta }}{R}_{{dp}}\equiv | {R}_{* }-| {\boldsymbol{R}}-{{\boldsymbol{R}}}_{* }\parallel$ represents how far away each point is from the gap center, and ${{\boldsymbol{R}}}_{* }$ is the star's position in the planet-centered frame. In all our simulations, we choose ${w}_{t}=0.25\,\mathrm{au}$, which is roughly 1/6–1/4 of the gap half-width. The velocity field is initialized by maintaining hydrostatic equilibrium in the radial direction by balancing the pressure gradient, the central stellar gravity, and centrifugal and Coriolis forces. We have tested that the gap shape is nicely maintained in the simulation where there is no planet at the center. We also initialize a Keplerian CPD with a very low surface density ($\sim {10}^{-4}$ g cm⁻²) within 0.4 Jupiter Hill radii.

Zoom In Zoom Out Reset image size

To maintain the gap structure during the simulation, we fix physical quantities at the outer boundary to be the initial values. At the inner boundary, we adopt the "outflow, no inflow" open boundary condition. The density and radial velocity are copied from the last active zones to the ghost zones. When the radial velocity is toward the computational domain, it is set to be zero. The azimuthal velocities in the ghost zones are set to their local Keplerian values. Internal energy in the ghost zones is the same as that in the last active zones.

Our initial and boundary conditions are quite different from those used in previous CPD simulations. Instead of capturing the planet gap-opening process in the simulation, we prescribe the gap density profile based on detailed gap-opening studies (e.g., de Val-Borro et al. 2006; Fung et al. 2014). By this approach, the gap shape is fixed, and the inflow from the circumstellar disk to the CPD can reach a steady state. In reality, the detailed gap width and depth depend on the disk turbulence level, the disk thermal structure, and even the magnetic field geometry (Zhu et al. 2013). In inviscid disks, the gap will even become deeper with time (Zhu et al. 2013). To simplify the gap-opening process, we use different gap widths and depths in different simulations to represent the gaps in disks with different properties or during different stages of the gap-opening process. We also use these different gap widths and depths to control the inflow rate to the CPD. In our simulations, we allow the circumplanetary region to settle down by itself and have no direct control of the inflow rate to the CPD. Thus, we take advantage of the fact that the simulation with a wider and deeper gap has a lower inflow rate to the CPD, and, through varying gap shapes, we indirectly control the inflow rate to cover four orders of magnitude.

All simulations are summarized in Table 1. The number after the capital letter "G" (where p2 means 0.2) represents the gap half-width (w_g) in the code unit (1 is 5 AU). The number after the capital letter S represents the circumstellar disk surface density (${{\rm{\Sigma }}}_{h}$) relative to the surface density in the fiducial case (Gp2). Our fiducial case (Gp2) has a circumstellar disk with 10 g cm⁻² surface density, which is one order of magnitude smaller than the surface density at 5 AU based on the minimum-mass solar nebula model (Hayashi 1981). The simulations having "iso" and "adi" in their names are isothermal and adiabatic simulations. All other simulations have considered simple radiative cooling (Equation (5)). For simulations with low circumstellar surface density and thus low inflow rates (Gp3Sp01 and Gp3Sp001), the CPDs take much longer to reach a steady state. Thus, we have to adopt a bigger inner boundary (4 Jupiter radii) to speed up the calculations and a higher initial CPD surface density to reduce the time to reach a steady state.

Table 1.  Models with Prescribed Gap Structure

Run	EoS	Resolution	${{\rm{\Sigma }}}_{h}$ a	${{\rm{\Sigma }}}_{l}$	wg	Time
		$R\times \phi$	g cm−2	g cm−2	au	year
Gp2	Cool	432 × 352	10	0.1	1	210
Gp2H	Cool	864 × 704	10	0.1	1	85
Gp2adi	Adi.	432 × 352	10	0.1	1	250
Gp3	Cool	432 × 352	10	0.1	1.5	250
Gp3iso	Iso.	432 × 352	10	0.1	1.5	297
Gp3Sp1	Cool	432 × 352	1	0.01	1.5	515
Gp3Sp01b	Cool	400 × 352	0.1	10−3	1.5	1216
Gp3Sp001b	Cool	400 × 352	0.01	10−4	1.5	2485

Notes.

^a ${{\rm{\Sigma }}}_{h}$ and ${{\rm{\Sigma }}}_{l}$ are the disk surface density outside and inside the gap, and w_g is the gap half-width. ^bThe inner boundary of the simulation domain is at 4 Jupiter radii instead of 2 Jupiter radii. The initial surface density of the CPD is set to be 10 g cm⁻², and the mass quickly drains to the planet at the early stages.

Download table as:  ASCIITypeset image

The typical disk surface density during the simulation is shown in Figure 1. The morphology of the CPD is very similar to previous studies (e.g., Lubow et al. 1999). Gas material in the circumstellar disk enters the Roche sphere of the planet when the gas undergoes horseshoe orbits around the planet. The planet excites spiral arms in the circumstellar disk, while the star excites spiral arms in the CPD. These spiral arms in CPDs are crucial for angular momentum transport in CPDs, as shown in Section 3.2.

3.1. Disk Structure

We have run these simulations until they reach a quasi-steady state, except for the isothermal case, which shows no sign of reaching a steady state during the whole run.⁵ The total mass of the CPD within the Hill sphere ($R\lt 0.347$ AU) is shown in Figure 2. When a quasi-steady state has been reached, the accretion rate onto the planet equals the inflow rate, and the disk mass remains constant. Disks with narrower (smaller w_g) gaps and higher surface density (higher ${{\rm{\Sigma }}}_{h}$ and ${{\rm{\Sigma }}}_{l}$) have a higher inflow rate from the circumstellar disk to the CPD. For the case Gp2, which has the highest accretion rate, the CPD reaches a steady state within 50 years. For cases with low accretion rates (e.g., Gp3Sp001), the disk reaches a quasi-steady state on the timescale of 1000 years. For a comparison with the dynamical timescale, the orbital time at our inner boundary (2 R_J) and half the Hill radii is only 8 hr and 2.3 years, respectively. Thus, 100 years, which is the typical duration of the simulation, is $1.1\times {10}^{5}$ innermost orbits.

Zoom In Zoom Out Reset image size

The accretion onto the planet is not strictly steady with time. As shown in Figures 2–4 and also in Ju et al. (2016), there are short-timescale variabilities. The disk accretion rates at the inner boundary and at R = 10 R_J are shown in the upper left panels of Figures 3 and 4 for both Gp2 and Gp3Sp001, respectively. The accretion rates in these panels have been averaged over 1 ${{\rm{\Omega }}}_{5\,\mathrm{au}}^{-1}$. For case Gp2, there is significant variability even within the timescale of 1 ${{\rm{\Omega }}}_{5\,\mathrm{au}}^{-1}$. Thus, we plot the accretion rates at a much higher cadence (averaged over 0.01 ${{\rm{\Omega }}}_{5\,\mathrm{au}}^{-1}$) within a short period of time. To understand the short-timescale variability, we plot the disk surface density profile at different times during the outburst, as shown in the upper right panels of Figures 3 and 4. The 2D surface density contours at these two times are shown in the lower left two panels of these figures (the lower left panel corresponds to the black curve in the upper right panel, while the lower middle panel corresponds to the red curve). We can see that the disk develops significant asymmetric structure during these accretion peaks. Vortices are generated in these disks first. As clearly shown in Figure 4, banded structures are first induced by spiral shocks, and these structures are subject to Rossby wave instability (Lovelace et al. 1999) and become vortices. In the lower left panel of Figure 3, a vortex forms at the outer disk and migrates inward. Eventually it crosses the inner boundary (lower middle panel), leading to an outburst. In Figure 4, the vortex does not migrate significantly in the disk, but it excites spiral arms, leading to additional angular momentum transport.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Thus, vortices provide some other ways to transport angular momentum in the disk, which are different from the accretion process that is due to the star-induced spiral shocks mentioned above. First, vortices can excite spiral density waves, which can become spiral shocks and transport angular momentum. Second, vortices themselves will migrate in disks (Paardekooper et al. 2010), and this vortex migration directly leads to mass transport since the vortex region has a high disk surface density. Thus, it is important to know how much disk accretion is driven by star-induced spiral shocks and how much accretion is due to vortices generated. Unfortunately, using Reynolds stress alone, as in Section 3.2, we cannot distinguish these accretion mechanisms. On the other hand, we can use the change of the Reynolds stress when a vortex appears or disappears in the disk to qualitatively study the role played by the vortex. As shown in the bottom right panels of Figures 3 and 4, the Reynolds stress does not change significantly when there is a vortex in the disk. This suggests that the spiral shocks play a more important, if not dominant, role in transporting angular momentum in disks.

On the other hand, vortices regulate the disk accretion rate, causing episodic accretion. And, for the disks with higher accretion rates, vortices are generated at smaller radii (by comparing Figures 3 and 4), leading to a shorter timescale and more violent outbursts.

To remove the short-timescale variability and study the general properties of disk structure, we average the mass accretion rate and surface density at each time step over 10 ${{\rm{\Omega }}}_{5\,\mathrm{au}}^{-1}$ (18 years, for the high-resolution case, Gp2H), 50 ${{\rm{\Omega }}}_{5\,\mathrm{au}}^{-1}$ (90 years, for fiducial resolution cases with 2 R_J inner boundary), or 100 ${{\rm{\Omega }}}_{5\,\mathrm{au}}^{-1}$ (180 years, for fiducial resolution cases with 4 R_J inner boundary) before the end of the simulations, and we plot the azimuthally averaged values in Figure 5. For reference, the Hill radius is ∼740 R_J. The time-averaged disk accretion rates in our models span a large range from 10⁻⁹ to 10⁻⁵ M_J yr⁻¹ (the upper left panel). In all cases except the isothermal case, the accretion rate is almost a constant throughout the whole disk, implying that the disk reaches a quasi-steady state. In the isothermal case, the disk has a low temperature, and the spiral shocks are so tightly wound up at the inner disk and spiral shocks inefficiently carry away angular momentum there (Ju et al. 2016, also in Section 3.2). Among all the models with radiative cooling, the disk surface density only spans one order of magnitude despite the large range of disk accretion rates. The maximum surface density is between 50 and 500 g cm⁻² in these models (the upper right panel).

Zoom In Zoom Out Reset image size

The disk midplane temperature and aspect ratio at the end of the simulations are shown in the bottom panels of Figure 5. The aspect ratio of the CPD is defined similar to the aspect ratio of the circumstellar disk. It is the ratio between the disk scale height (H) and the radius (R) in the CPD, which is equal to ${c}_{s}/{v}_{\phi }$ with ${v}_{\phi }$ being the orbital velocity around the planet. The disk midplane temperature varies significantly among cases. In simulations that have realistic cooling, the maximum disk temperature ranges from 200 to 5000 K. The disk with a higher accretion rate has a higher temperature. Due to the large temperature range in different cases, the disk aspect ratio also varies significantly, from 0.05 in the cases with low accretion rates to 0.2 in the cases with high accretion rates.

Both the isothermal and adiabatic cases are significantly different from other cases with realistic radiative cooling. The isothermal case has not reached a steady state. The mass is continuously being piled up in the disk. The adiabatic case has a too-high temperature with H/R larger than 0.5 at $R\gt \,10{R}_{J}$, so the disk cooling treatment is inaccurate (Equation (5)). The disk is so hot that it also limits the inflow rate to the CPD. The infall rate in the adiabatic case is more than one order of magnitude smaller than the infall rate in the isothermal case. After normalizing the accretion rate with the disk surface density, α is, however, larger in the adiabatic case than the isothermal case since shocks can propagate to the inner disk more easily.

3.2. Accretion Mechanism

In our inviscid hydrodynamical simulations, dissipation of spiral shocks is the main angular momentum transport mechanism. To show how spiral shocks propagate in CPDs, we plot the temperature color contours from two simulations (Gp2 and Gp3Sp001) with dramatically different disk accretion rates in Figure 6. Two spiral arms are apparent in the figure. Spiral waves are excited by the tidal force at Lindblad resonances, ${R}_{\mathrm{Lind}}={R}_{0}{(1\pm 1/m)}^{2/3}{(1+q)}^{-1/3}$, where R₀ is the separation between the perturber and the central source, m is the order of the resonance, and q is the mass ratio between the perturber and the central source. In our case with the star being the perturber, q is 1000, so only the inner m = 2 Lindblad resonance at ${R}_{\mathrm{Lind}}={(1/4004)}^{1/3}{R}_{0}$ is within the Hill radius of the planet ${(1/3000)}^{1/3}{R}_{0}$, which is why two spiral arms are excited in Figure 6. The strength of the excited spiral waves is strongly correlated with the disk surface density at Lindblad resonances. After spiral waves are excited, they will propagate to the inner disk. During the propagation, linear spiral waves can quickly steepen into spiral shocks (Goodman & Rafikov 2001) and dissipate. On the other hand, these spiral shocks are rather weak, and the shock front only deviates slightly from the spiral wave front calculated by the linear dispersion relation (Zhu 2015). As in Ju et al. (2016), we plot the wave front by integrating the linear dispersion relation

$$\begin{eqnarray}&&\displaystyle \frac{d\phi }{{dR}}=-\displaystyle \frac{1}{{c}_{s,a}(R)}\sqrt{{({\rm{\Omega }}-{{\rm{\Omega }}}_{p})}^{2}-{\hat{\kappa }}^{2}/{m}^{2}}.\end{eqnarray} \tag{ 7 }$$

shown in the bottom panels of Figure 6 as the dotted curves, where $\hat{\kappa }$ is the epicyclic frequency and $\hat{\kappa }={\rm{\Omega }}$ in a Keplerian rotating disk. Here, m is assumed to be 2, and ${c}_{s,a}(R)$ is the sound speed calculated using the adiabatic equation of state and the azimuthally averaged temperature at each R.⁶ Figure 6 clearly shows that the linear dispersion relation reproduces the spiral shock front very well.

Zoom In Zoom Out Reset image size

The openness of the spiral arm reflects the disk temperature and is also related to the efficiency of angular momentum transport by spiral shocks (Ju et al. 2016). This is clearly demonstrated in Figure 6. The slope of the curves in the bottom panels of Figure 6 is directly related to the disk pitch angle:

$$\begin{eqnarray}&&\cot \theta =\left|\,\displaystyle \frac{{Rd}\phi }{{dR}}\,\right|=\left|\,\displaystyle \frac{d\phi }{d\mathrm{ln}R}\,\right|.\end{eqnarray} \tag{ 8 }$$

With Equations (7) and (8), we see that tan$\theta \propto {c}_{s,a}$. Thus, spiral arms in a hot disk (the left panels in Figure 6) have larger pitch angles and are more open than spiral arms in a cold disk (the right panels). When the spiral arm is more open, it is easier to propagate to the inner disk and should lead to stronger angular momentum transport there (Ju et al. 2016). Thus, we expect a more efficient angular momentum transport and higher α value in the left-hand panels than in the right-hand panels.

To understand the accretion efficiency, we follow Ju et al. (2016) to separate the angular momentum budget in the angular momentum equation. The derivations below are basically the same as that in Shakura & Sunyaev (1973) and Balbus & Hawley (1998), but we wrote them in quantities that can be directly calculated from numerical simulations. After being averaged over the azimuthal direction, the angular momentum equation is reduced to

$$\begin{eqnarray}&&{\partial }_{t}\langle {\rm{\Sigma }}{{Rv}}_{\phi }{\rangle }_{\phi }=-\displaystyle \frac{1}{R}{\partial }_{R}({R}^{2}\langle {\rm{\Sigma }}{v}_{R}{v}_{\phi }{\rangle }_{\phi })+\langle {\boldsymbol{R}}\times {\boldsymbol{F}}{\rangle }_{\phi }\end{eqnarray} \tag{ 9 }$$

where $\langle X{\rangle }_{\phi }$ represents ${\int }_{0}^{2\pi }{Xd}\phi$, and ${\boldsymbol{F}}$ is the external force as in Equation (3). After subtracting the Keplerian motion of the disk with $\delta {v}_{\phi }={v}_{\phi }-{v}_{K}$, we have

$$\begin{eqnarray}&&{\partial }_{t}\langle {\rm{\Sigma }}{\rangle }_{\phi }({{Rv}}_{K})+{\partial }_{t}\langle {\rm{\Sigma }}R\delta {v}_{\phi }{\rangle }_{\phi }\\ &&\quad =-\displaystyle \frac{1}{R}{\partial }_{R}\langle R{\rm{\Sigma }}{v}_{R}{\rangle }_{\phi }({{Rv}}_{K})-\displaystyle \frac{1}{R}\langle {\rm{\Sigma }}{{Rv}}_{R}{\rangle }_{\phi }{\partial }_{R}({{Rv}}_{K})\\ &&\qquad -\displaystyle \frac{1}{R}{\partial }_{R}({R}^{2}\langle {\rm{\Sigma }}{v}_{R}\delta {v}_{\phi }{\rangle }_{\phi })+\langle {\boldsymbol{R}}\times {\boldsymbol{F}}{\rangle }_{\phi }.\end{eqnarray} \tag{ 10 }$$

After canceling out the first terms on the left-hand and right-hand sides using the mass continuity equation, we have

$$\begin{eqnarray}&&{\partial }_{t}\langle {\rm{\Sigma }}R\delta {v}_{\phi }{\rangle }_{\phi }=-\displaystyle \frac{1}{R}\langle {\rm{\Sigma }}{{Rv}}_{R}{\rangle }_{\phi }{\partial }_{R}({{Rv}}_{K})\\ &&\qquad \qquad \qquad -\displaystyle \frac{1}{R}{\partial }_{R}({R}^{2}\langle {\rm{\Sigma }}{v}_{R}\delta {v}_{\phi }{\rangle }_{\phi })+\langle {\boldsymbol{R}}\times {\boldsymbol{F}}{\rangle }_{\phi }.\end{eqnarray} \tag{ 11 }$$

The left-hand side is the time derivative of the perturbed angular momentum. To be consistent with Ju et al. (2016), we multiply this quantity by R² and define it as

$$\begin{eqnarray}&&{\mathrm{AM}}_{t}(R)={R}^{2}{\partial }_{t}\langle {\rm{\Sigma }}R\delta {v}_{\phi }{\rangle }_{\phi }.\end{eqnarray} \tag{ 12 }$$

The first term on the right-hand side is the angular momentum change due to mass accretion. After being multiplied by R², it is defined as

$$\begin{eqnarray}&&{\mathrm{AM}}_{\dot{M}}(R)=\dot{M}R{\partial }_{R}({{Rv}}_{K}),\end{eqnarray} \tag{ 13 }$$

where $\dot{M}=-\langle {\rm{\Sigma }}{{Rv}}_{R}{\rangle }_{\phi }$. The second term on the right-hand side is the Reynolds stress gradient or the angular momentum flux gradient. With the additional factor of R², it is defined as

$$\begin{eqnarray}&&{\mathrm{AM}}_{\mathrm{FH}}=-R{\partial }_{R}({R}^{2}\langle {\rm{\Sigma }}{v}_{R}\delta {v}_{\phi }{\rangle }_{\phi }).\end{eqnarray} \tag{ 14 }$$

The last term is the torque by all forces. With the additional factor of R², it is defined as

$$\begin{eqnarray}&&T(R)={R}^{2}\langle {\boldsymbol{R}}\times {\boldsymbol{F}}{\rangle }_{\phi }.\end{eqnarray} \tag{ 15 }$$

For steady linear waves that are generated by the torque and propagate in disks, the ${\mathrm{AM}}_{\dot{M}}$ term is zero since linear perturbation cannot change the background state, and $\dot{M}$ is zero. In this case, all angular momentum generated by the torque (the T(R) term) is carried away by the wave (the AM_FH term). However, when the waves steepen into shocks, part of the angular momentum flux carried by the wave is lost into the background disk so that the disk accretes. Thus, shock dissipation can be quantified by the difference between the T(R) term and the AM_FH term. When the disk is in a quasi-steady state as in our cases, the terms on the right-hand side of Equation (11) need to balance each other, so

$$\begin{eqnarray}&&-{\mathrm{AM}}_{\dot{M}}(R)\approx {\mathrm{AM}}_{\mathrm{FH}}(R)+T(R).\end{eqnarray} \tag{ 16 }$$

This equation suggests that accretion is not led by the torque alone but by the dissipation $({\mathrm{AM}}_{\mathrm{FH}}(R)+T(R))$, which is defined as the difference between the total torque and the angular momentum flux carried away by the wave.

The detailed angular momentum budget for case Gp2 is shown in Figure 7, where all quantities are averaged over each time step for 50 ${{\rm{\Omega }}}_{5\,\mathrm{au}}^{-1}$. The dissipation term ${\mathrm{AM}}_{\mathrm{FH}}(R)+T(R)$ is balanced by the disk accretion ${\mathrm{AM}}_{\dot{M}}(R)$. Figure 7 also shows that most of the torque is exerted on the disk at $R\gt 50$ R_J, and part of the torque is used to launch spiral waves. These spiral arms can propagate all the way to Jupiter's surface and dissipate throughout the whole disk with a nonzero dissipation term.

Zoom In Zoom Out Reset image size

When the disk is in a quasi-steady state (as shown in Figure 5), we can directly integrate Equation (11) to derive the relationship between the mass accretion rate and the stress:

$$\begin{eqnarray}&&\dot{M}=\displaystyle \frac{{R}^{2}\langle {\rm{\Sigma }}{v}_{R}\delta {v}_{\phi }{\rangle }_{\phi }+C-\int R\langle {\boldsymbol{R}}\times {\boldsymbol{F}}{\rangle }_{\phi }{dR}}{{{Rv}}_{K}},\end{eqnarray} \tag{ 17 }$$

where C is determined by the boundary condition. Motivated by the standard α disk theory where $\dot{M}=\alpha 3\pi {\rm{\Sigma }}{c}_{s}H$, we divide $3\pi {\rm{\Sigma }}{c}_{s}H$ in the above equation to write

$$\begin{eqnarray}&&{\alpha }_{\mathrm{eff}}={\alpha }_{\mathrm{Rey}}+{\alpha }_{\mathrm{const}}+{\alpha }_{T},\end{eqnarray} \tag{ 18 }$$

where the effective α calculated by $\dot{M}$ is

$$\begin{eqnarray}&&{\alpha }_{\mathrm{eff}}=\displaystyle \frac{\dot{M}}{3\pi {\rm{\Sigma }}{c}_{s}H},\end{eqnarray} \tag{ 19 }$$

the pressure-normalized Reynolds stress is

$$\begin{eqnarray}&&{\alpha }_{\mathrm{Rey}}=\displaystyle \frac{\langle {\rm{\Sigma }}{v}_{R}\delta {v}_{\phi }{\rangle }_{\phi }}{3\pi {\rm{\Sigma }}{c}_{s}^{2}},\end{eqnarray} \tag{ 20 }$$

and the α values due to the boundary condition and the torque are

$$\begin{eqnarray}&&{\alpha }_{\mathrm{const}}=\displaystyle \frac{C}{3\pi {\rm{\Sigma }}{c}_{s}{{HRv}}_{K}},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{\alpha }_{T}=\displaystyle \frac{-\int R\langle {\boldsymbol{R}}\times {\boldsymbol{F}}{\rangle }_{\phi }{dR}}{3\pi {\rm{\Sigma }}{c}_{s}{{HRv}}_{K}}.\end{eqnarray} \tag{ 22 }$$

These α values with respect to the radius are shown in Figure 8 for all of the cases that have radiative cooling. The isothermal and adiabatic cases are shown in Figure 9. For the isothermal case, the black and red curves do not overlap, suggesting that the disk has not reached a steady state. The adiabatic case has profiles similar to Gp2. To derive these α values, various components in Equation (17) are averaged over each time step for 10 ${{\rm{\Omega }}}_{5\,\mathrm{au}}^{-1}$ (for the high-resolution case, Gp2H), 50 ${{\rm{\Omega }}}_{5\,\mathrm{au}}^{-1}$ (for cases with 2 R_J inner boundary), or 100 ${{\rm{\Omega }}}_{5\,\mathrm{au}}^{-1}$ (for cases with 4 R_J inner boundary) before the end of the simulations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

All these α values can be directly measured in the simulations except ${\alpha }_{C}$, which depends on C and the boundary condition used in the simulation. We notice that ${\alpha }_{\mathrm{eff}}$ is far larger than ${\alpha }_{\mathrm{Rey}}$ and ${\alpha }_{T}$ at the inner boundary for all the cases in Figure 8. Thus, we can apply the zero-stress boundary condition to our simulations. Specifically, we use ${\alpha }_{\mathrm{eff}}={\alpha }_{\mathrm{const}}$ (Equation (18)) at the inner boundary to derive $C=\dot{M}{R}_{\mathrm{in}}{v}_{K,\mathrm{in}}$ in Equation (17), where R_in and ${v}_{K,\mathrm{in}}$ are the radius and the Keplerian velocity at the inner radius. Then with this C we can calculate ${\alpha }_{C}$ shown in Figure 8.

We want to emphasize that the zero-stress boundary in our simulations is exactly the same boundary condition used in the thin-disk theory. The nonzero ${\alpha }_{C}$ is the generic property of the zero-stress boundary condition. The large ${\alpha }_{C}$ value at the boundary does not mean that the disk accretion there is induced by our adopted numerical boundary condition in the simulation. Even using the thin-disk theory, an analytical model with a large ${\alpha }_{C}$ at the boundary can be constructed as long as Σ goes to small values at the boundary. In order to understand the effects of the boundary condition, we need to choose different boundary conditions or simulate the boundary layer directly, which will be the focus of a future paper.

Figure 8 shows that ${\alpha }_{\mathrm{eff}}$ is quite close to ${\alpha }_{\mathrm{Rey}}$ in most parts of the disk. The effect of ${\alpha }_{T}$ only becomes noticeable starting from half of the Hill radius outward (the Hill radius is 740 R_J). Here, ${\alpha }_{\mathrm{const}}$ is important at the inner disk within 10 R_J, and ${\alpha }_{\mathrm{eff}}$ is at the 0.001–0.01 level despite disk accretion rates that span four orders of magnitude among different models.

The α value is slightly higher in a disk with a higher accretion rate. This is due to the feedback from shock dissipation to disk accretion. With a higher infall rate, more material is present within the Roche sphere (or Hill sphere). As shown in Figure 5, the disk is larger with a higher infall rate. Since spiral waves are more easily excited in disk regions closer to the Lindblad resonances, the larger disks have stronger spiral waves. After these stronger waves dissipate in disks, they heat up the disk more significantly, and the disk temperature becomes higher. Since spiral arms are more open in hotter disks (Figure 6), they propagate farther into the inner disk and lead to stronger accretion there.

Quantitatively, ${\alpha }_{\mathrm{eff}}$ is ∼0.01 in Gp2, which has $\dot{M}\sim {10}^{-5}$ ${M}_{J}\,{\mathrm{yr}}^{-1}$ and $H/R\sim 0.25$, while ${\alpha }_{\mathrm{eff}}$ is ∼0.001 in Gp3Sp01, which has $\dot{M}\sim 4\times {10}^{-8}$ ${M}_{J}\,{\mathrm{yr}}^{-1}$ and $H/R\sim 0.1$. Larson (1990) derived the α value from a steady, self-similar shock (see also Spruit 1987) as

$$\begin{eqnarray}&&\alpha \sim 0.013{[{({c}_{s}/{v}_{K})}^{3}+0.08{({c}_{s}/{v}_{K})}^{2}]}^{1/2}.\end{eqnarray} \tag{ 23 }$$

Using this equation, we can estimate α to be ∼0.002 when ${c}_{s}/{v}_{K}=H/R\sim 0.25$, and α to be ∼0.0006 when $H/R\sim 0.1$. These estimates are similar to the simulated values within the same order of magnitude, which suggest that angular momentum transport by spiral shocks can indeed be efficient in CPDs. On the other hand, these estimated values are a factor of five and two times smaller than those derived in simulations. Shocks in our simulations are neither steady nor self-similar. To understand the difference between simulated values and analytical estimates, a more detailed analytical model needs to be developed.

3.3. Energy Budget

4.1. Observational Signatures

With the CPD structure uniquely determined in our simulations, we can calculate the observables for CPDs and test if current telescopes can detect CPDs.

4.1.1. SEDs at Optical and Near-infrared

Detailed spectral energy distributions (SEDs) of accreting CPDs have been calculated by Zhu (2015) when the disk extends to the planet surface or when the disk is truncated by magnetic fields of the young planets. These SEDs are much redder than SEDs produced by brown dwarfs or planets (e.g., Spiegel & Burrows 2012) since the outer disk emits significant infrared flux. Thus, to use direct imaging to find the accretion disks around low-mass planets and distinguish them from brown dwarfs or hot high-mass planets, it is crucial to obtain photometry at mid-infrared bands ($L^{\prime}$, M, N bands).

Direct imaging observations have found several red sources within CPDs (Kraus & Ireland 2012; Quanz et al. 2013; Biller et al. 2014; Reggiani et al. 2014; Currie et al. 2015). Because of their redness, these sources could be potential candidates for CPDs. Thus, we have used the theoretical SEDs in Zhu (2015) to fit the observed photometry, trying to constrain the disk properties. Table 2 summarizes recent detections of point sources in protoplanetary disks and the predictions from theoretical models. We fit the L-band observations using the L-band magnitudes given in Zhu (2015) to constrain the $M\dot{M}$. Then we give predicted magnitudes at other wavelength bands. The distance and disk inclination are from the given references for each source. To simplify the fitting procedure, we only use the SED models with the disk inner radius of 1 Jupiter radius. Even with this simplification, the agreement between the model and observations at various bands is quite good, especially considering that L-band observation is the only one used in the model to make predictions at all of the other bands. If the assumption of the inner disk radius is relaxed, multiple bands of data can be fitted much better (e.g., as done in Currie et al. 2015 for HD 100546b; extinction also plays a role at the Hband for HD 100546b). Generally, for all these CPD candidates, $M\dot{M}$ is around 10⁻⁶–10⁻⁵ ${{\rm{M}}}_{J}^{2}\,{\mathrm{yr}}^{-1}$.

Table 2.  Magnitudes of Accreting Circumplanetary Disks Compared with Observations

	$M\dot{M}({M}_{J}^{2}\,{\mathrm{yr}}^{-1})$	J	H	K	L'	M	N	870 μm	1.3 mm
HD169142b
Obs.a		>13.8			12.2 ± 0.5
The.	10−5	14.8	14.66	13.82	12.2	11.62	10.12	∼300 μJy	∼100 μJy
HD100546b
Obs.b			19.4 ± 0.32	>15.43 ± 0.06	13.92 ± 0.1	13.33 ± 0.16
The.	2 × 10−6	20.66	18.41	16.50	13.9	13.05	11.37	∼800 μJy	∼300 μJy
LkCa 15b
Obs.c				14.2 ± 0.5	13.2±0.5
The.	7 × 10−6	16.09	15.92	14.96	13.2	12.54	11.04	∼300 μJy	∼100 μJy

Notes.

^aReggiani et al. (2014), distance of 145 pc, inclination of 0°. ^bQuanz et al. (2015), Currie et al. (2015), distance of 100 pc, inclination of 40°. ^cSallum et al. (2015), distance of 145 pc, inclination of 50°.

Download table as:  ASCIITypeset image

To derive SEDs for CPDs, Zhu (2015) has used the temperature profile from the thin-disk theory as in Equation (45), where R_in is the inner radius of the disk (it is equal to R_p if the disk extends to the planet surface or is equal to the magnetic truncation radius when the planet has strong magnetic fields and is undergoing magnetospheric accretion; Lovelace et al. 2011).⁷

However, as shown in the Appendix, the cooling rate from the thin-disk approximation (Equation (45)) is inadequate to reproduce the real cooling rate from the simulation due to the radial energy advection. To understand how this difference affects the disk's SED, we generate the SEDs directly from simulations. First, we use the cooling rate in simulations to calculate the disk effective temperature at each radius. Then we assume blackbody emission at that radius to calculate the SED generated by that annulus of the disk. Finally, we integrate all of the disk emission from the inner boundary to the planet's Hill radius to get the SED of the entire CPD. These SEDs are plotted in Figure 10 as solid curves. For comparison, two SEDs calculated from the thin-disk approximation (Equation (45)) are also plotted: the dotted curves use the same mass accretion rates as those measured from simulations, while the dash-dotted curves assume only half of the disk accretion rates. Clearly, the dotted curves roughly fit the simulated SEDs at the short wavelengths (below where the flux peaks). The dash-dotted curves fit the simulated SEDs at the long wavelengths (the Rayleigh–Jeans part of the SEDs). This is expected from Figure 15, where at the outer disk the emission is only half of that from the thin-disk theory. Thus, if we use the Rayleigh–Jeans part of the observed SED to estimate the disk accretion rate based on the thin-disk theory, we will underestimate the disk accretion rate by a factor of two. Thus, the mass accretion rates given in Table 2 need to be multiplied by a factor of two to derive the real disk accretion rates.

Zoom In Zoom Out Reset image size

The dashed curves in Figure 10 are SEDs calculated using our empirical temperature relationship (Equation (46)). This empirical relationship works well for the SED at long wavelengths but it overpredicts the flux at short wavelengths, due to the reason presented at the end of the Appendix.

Note that SEDs in Figure 10 are only from the accretion disks. The planet's SEDs have not been added. If the planet has a cold atmosphere, as in the "cold-start" model (Marley et al. 2007), its flux is significantly lower than the flux from the disk (Zhu 2015). Since the zero-stress boundary condition has been used to derive the temperature in the thin-disk theory and this boundary condition also stands in our simulations, it implicitly assumes that the planet rotates at the breakup speed and the disk joins the planet smoothly. In reality, the planet may rotate slower than the disk and a boundary layer is formed around the planet, or the planet has a strong magnetic field to truncate the CPD. The observational signatures of the boundary layer and magnetospheric accretion will be discussed in Section 4.1.3.

4.1.2. Submillimeter to Millimeter Flux at ALMA/EVLA Bands

With the Atacama Large Millimeter/submillimeter Array (ALMA) and the Very Large Array (VLA)'s great sensitivity, we may detect CPDs at the submillimeter to millmeter bands. Isella et al. (2014) have used parameterized models to calculate submillimeter to millimeter flux for CPDs and suggest that ALMA can probe CPDs with mass down to 5 × 10⁻⁴ M_J. Due to the use of the parameterized models, disk mass, size, and accretion rate are degenerate, leaving a large parameter space to explore. In our first-principle calculations, the disk density and temperature structure are uniquely determined at a given accretion rate, so we can give a unique prediction of the CPDs' submillimeter to millimeter flux.

To calculate the disk's flux at submillimeter to millimeter bands, we use the dust opacity of $0.034\,\times (0.87\,\mathrm{mm}/\lambda )$ cm² g⁻¹ (Andrews et al. 2012). At one annulus, if the disk is optically thin at the given wavelength, we approximate the brightness temperature (T_b) as the product of the midplane temperature and the optical depth. When the disk is optically thick at the given wavelength, the brightness temperature is chosen as the temperature where the disk becomes optically thick at that wavelength. In detail, the temperature at the optical depth of ${\tau }_{R}$ is

$$\begin{eqnarray}&&T{({\tau }_{R})}^{4}=\displaystyle \frac{3}{16}{T}_{\mathrm{eff}}^{4}{\tau }_{R}+{T}_{\mathrm{ext}}^{4},\end{eqnarray} \tag{ 24 }$$

where ${\tau }_{R}$ is the optical depth calculated using the Rosseland mean opacity. At the disk height where ${\tau }_{R}={\kappa }_{R}/{\kappa }_{\lambda }$, the disk becomes optically thick at λ. If we plug the effective temperature from Equation (5) into Equation (24), we can derive

$$\begin{eqnarray}&&{T}_{b}^{4}=T{({\tau }_{R}={\kappa }_{R}/{\kappa }_{\lambda })}^{4}=({T}_{c}^{4}-{T}_{\mathrm{ext}}^{4})/{\tau }_{\lambda }+{T}_{\mathrm{ext}}^{4}.\end{eqnarray} \tag{ 25 }$$

where ${\tau }_{\lambda }={\rm{\Sigma }}{\kappa }_{\lambda }/2$.⁸ Then knowing T_b at each radius, we can integrate the emission from the whole disk to derive the total flux. The flux at various wavelengths for disks with different accretion rates is given in Table 3. The source is assumed to be 100 pc away from us.

Table 3.  Disk Properties and Derived Submillimeter to Millimeter Flux

Run	$\dot{M}$	CPD Mass	Fluxa 870 μm	1.3 mm	7 mm
	${{\rm{M}}}_{J}\,{\mathrm{yr}}^{-1}$	MJ	μJy	μJy	μJy
Gp2	1.07 × 10−5	6.4 × 10−4	200	79	1.2
Gp3	2.26 × 10−6	3.1 × 10−4	84	34	0.51
Gp3Sp1	3.89 × 10−7	1.0 × 10−4	31	12	0.12
Gp3Sp01	3.93 × 10−8	3.6 × 10−5	12	4.0	0.028
Gp3Sp001	1.92 × 10−9	2.0 × 10−5	6	2.1	0.013
Gp3iso	⋯	1.0 × 10−3	44	16	0.28

Note.

^aThe total received flux is calculated by assuming the source is at 100 pc.

Download table as:  ASCIITypeset image

Although the CPDs in our simulations are 5 AU away from their central stars, we can roughly scale our results to CPDs at other distances. Based on Figure 5, the disk surface density is almost a constant along radii, and the midplane temperature roughly scales as ${R}^{-1}$. Thus the total intensity scales as R. Since the size of the CPD is one-third of the planet's Hill radius (Martin & Lubow 2011a), the total submillimeter to millimeter flux in Table 3 should roughly scale as the distance of the planet to the central star.

Thus, we scale the derived submillimeter to millimeter flux for Jupiter at 5 AU in Table 3 to several CPD candidates and provide their fluxes in Table 2. We want to caution that this scaling is highly simplified. Intrinsically, it assumes that α in the CPD around Jupiter is the same as CPDs around young planets at 20–50 AU. Considering the complicated interplay between the thermodynamics, spiral excitation, and shock dissipation, this is unlikely to be true. Furthermore, we implicitly assume that the CPD candidates all have a 1 Jupiter-mass planet at the center, the same as in our simulations. Thus, the submillimeter flux given in Table 2 can only be considered as an order of magnitude estimate. Nevertheless, with ALMA's great sensitivity (assuming 1.4 mm and all bandwidths set for continuum, we can detect 230 μJy point sources with the signal-to-noise ratio of 10 in 57 minutes; private communication with John Tobin), we should easily detect these sources if their origins are CPDs.

4.1.3. Other Observables and Variability

Our simulations suggest that CPDs can accrete quite efficiently due to the spiral shocks. Besides the SED, optical/near-infrared emission lines (e.g., Hartmann et al. 1994; Muzerolle et al. 1998, 2001), ultraviolet excess (e.g., Gullbring et al. 2000; Herczeg & Hillenbrand 2008; Ingleby et al. 2013), and line veiling (e.g., Calvet & Gullbring 1998) are all accretion signatures (Bouvier et al. 2007, p. 479; Rigliaco et al. 2012; Alcalá et al. 2014). For disks around young stellar objects, the profiles of emission lines suggest that these lines are produced by infalling material during magnetospheric accretion (Hartmann et al. 1994). Ultraviolet excess is also believed to be produced by the accretion hot spots on the surface of the star where the magnetosphere connects to the star (Calvet & Gullbring 1998). If CPDs around young planets can undergo magnetospheric accretion, they should also produce these accretion features. However, in order to truncate accreting CPDs by magnetic fields, young planets need to have magnetic fields of ∼10–100 Gauss (Fendt 2003; Lovelace et al. 2011; Zhu 2015), much stronger than Jupiter's current magnetic field. If young planets do not have such strong magnetic fields, the disk will accrete onto the planet through the boundary layer. Although accretion through the boundary layer has been better understood (Belyaev et al. 2013; Philippov et al. 2016), boundary layers in CPDs can be quite thick (H/R ∼ 0.2, Figure 5) and directly exchange energy with the planet (Owen & Menou 2016). Whether boundary layers can still produce strong emission lines and UV excess needs to be explored.

Our simulations also suggest that the disk accretion rate varies significantly at the timescale of years when the disk accretion rate is $\sim {10}^{-5}\,{M}_{J}\,{\mathrm{yr}}^{-1}$. The disk cooling timescale is also around one year for such a disk, estimated with Equation (8) in Zhu (2015). Thus the change of the disk accretion rate would lead to the variability of the observables over a timescale of years. Thus, for these disks, the accretion tracers may vary between observations at different epochs. This is shown in Figure 11, where the dotted and dashed curves are only separated by 1.8 years but the peak intensity is different by almost a factor of 10.

Zoom In Zoom Out Reset image size

Besides disk accretion, CPDs could also launch jets or outflows (Quillen & Trilling 1998; Gressel et al. 2013) and have shock fronts that are due to the inflow from the circumstellar disks (Tanigawa et al. 2012). They may also have observational signatures. Last but not least, the Keplerian rotation of the CPD may even be probed by molecular lines using ALMA (Perez et al. 2015).

4.2. Satellite Formation

We can compare the disk surface density from our simulations with the disk surface density required by various satellite formation models. Our derived disk surface density is orders of magnitude smaller than the in situ satellite formation models or the so-called "minimum-mass subnebula" model (Figure 12). Our simulated CPDs have surface densities between 10 and 1000 g cm⁻² depending on the disk accretion rates, while the "minimum-mass subnebula" (Lunine & Stevenson 1982; Coradini & Magni 1984; Mosqueira & Estrada 2003) has a surface density of a few $\times \,{10}^{5}$ g cm⁻². The model in Mosqueira & Estrada (2003) is plotted as the dotted curve in Figure 12, which is much higher than the surface density from our 2D simulations. On the other hand, the surface densities in our simulations are consistent with that in the "gas-starved" satellite formation model (Canup & Ward 2002, 2006; Ward & Canup 2010).

Zoom In Zoom Out Reset image size

Since Galilean satellites are ice-rich, they need to form in conditions where the disk midplane temperature is below 150 to 200 K, depending on the pressure (Prinn & Fegley 1989, p. 78; Canup & Ward 2002). In order for the disk to be cooler than 200 K at $R\gt 30\,{R}_{J}$, Figure 5 suggests that the disk accretion rate needs to be smaller than $4\times {10}^{-8}\,{M}_{J}\,{\mathrm{yr}}^{-1}$ or ${10}^{-5}\,{M}_{\oplus }\,{\mathrm{yr}}^{-1}$. Thus, Galilean satellites should form in a CPD with a very low accretion rate. This accretion rate is broadly consistent with that (2 × 10⁻⁷ M_J yr⁻¹) suggested by Canup & Ward (2002). Using α ∼ 0.001, Heller & Pudritz (2015a, 2015b) have reached the similar conclusion that Jovian moons have to form during the final stages of CPD accretion. Despite this low accretion rate, it will only take $\sim 5\times {10}^{5}$ years, well within the protoplanetary disk lifetime, for such a CPD to supply enough material for building Galilean satellites.

Although the effects of spiral shocks can be represented by the equivalent α parameter as shown above, spiral shocks are substantially different from viscosity in that they can travel a long distance in disks and have a nonhomogeneous azimuthal structure. The density and temperature are high at the shock front. Since the spiral shocks are stationary in the corotating frame, material in the CPD will move in and out of the spiral shocks during each orbit around the planet. Disk material can frequently be heated up, dissociated, and condensed again. We plot the 1D density and temperature profiles along the azimuthal direction at 10 and 30 R_J in Figure 13. Thus, disk material is heated episodically (twice every orbit due to the two spiral arms) and may have imprints in the material accreted onto the satellites or left in the CPD.

Zoom In Zoom Out Reset image size

4.3. Caveats

CPDs control the growth of planets, supply material for satellites to form, and provide observational signatures of young forming planets. In this paper, we provide a new mechanism to explain how CPDs accrete.

We have carried out two-dimensional hydrodynamical simulations to study CPDs using Athena++. Simple radiative cooling has been considered in the simulation. Different from most previous simulations, we choose a coordinate system centered on the planet, which significantly reduces the numerical error and enables simulations with small numerical viscosity. Our simulation domain extends from the circumstellar disk all the way to Jupiter's surface. A gap in the circumstellar disk is prescribed, allowing us to control the inflow rate from the circumstellar disk to the CPD.

Two spiral shocks are present in CPDs, induced by the tidal force from the central star. We find that these spiral shocks can lead to significant angular momentum transport and energy dissipation in CPDs. Meanwhile, dissipation and heating by spiral shocks have a positive feedback on shock-driven accretion itself. As the disk is heated up by spiral shocks, the shocks become more open and propagate farther into the inner disk, leading to more efficient angular momentum transport at the inner disk. On the other hand, shock-driven accretion cannot guarantee strict steady disk accretion since the angular momentum transport depends on the global wave propagation and local shock dissipation. Mass will pile up in some regions of the disk, and sometimes vortices are produced, which produce short-timescale variability. The disk is adjusting itself through these variabilities and tries to maintain a quasi-steady state. Eventually, a quasi-steady state of accretion flow is reached in our simulations from the planet's Hill radius all the way to the planet surface. After averaging quantities over a long timescale, the angular momentum budget is carefully analyzed. The effective α coefficient characterizing angular momentum transport due to spiral shocks is ∼0.001–0.02 even though the disk accretion rates span four orders of magnitude. The α value is higher in a disk with a higher accretion rate due to the shock-heating feedback.

With an energy budget analysis, we show that radial advection of energy becomes important, and the disk generates less infrared radiation than that from the thin-disk approximation by a factor of about two. Thus, if we use the infrared flux calculated from the thin-disk approximation to derive the accretion rate of CPDs, we may underestimate the disk accretion rate by a factor of two.

Finally, we calculate the flux from CPDs at ALMA and EVLA wavelength bands and predict the flux for several recent CPD candidates (e.g., HD169142b, HD100546b, LkCa 15b). These CPD candidates should be relatively bright at ALMA wavelength bands. In future, ALMA should be able to discover many accreting CPDs.

Although our simulations are limited to 2D, the possibility that we may have understood how CPDs accrete and its huge implications for observations and satellite formation make it worth being studied in detail in future.

We thank Lee Hartmann and Nurial Calvet for useful comments. We thank the referee for a helpful and constructive report. All hydrodynamical simulations are carried out using computers 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. Part of the work is done when ZZ is at Princeton, supported by NASA through Hubble Fellowship grant HST-HF-51333.01-A 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.

Besides studying angular momentum transport by spiral shocks, we also investigate how energy is dissipated and transported in CPDs. Understanding the energy budget during the accretion process is essential for studying the disk thermal structure and the observational signatures of disks (e.g., spectral energy distributions). CPDs are geometrically quite thick, with the aspect ratio reaching 0.3 (Figure 5), so local cooling may not perfectly balance local viscous/shock heating, and radial advection of energy can be important.

To calculate the energy budget, we take a similar approach as for the angular momentum equation by averaging the energy equation in the azimuthal direction:

$$\begin{eqnarray}&&{\partial }_{t}\langle E{\rangle }_{\phi }+\displaystyle \frac{1}{R}{\partial }_{R}[\langle R(E+P){v}_{R}{\rangle }_{\phi }]\\ &&\quad ={-\langle {\rm{\Sigma }}{\rm{\nabla }}{{\rm{\Phi }}}_{p}\cdot {\boldsymbol{v}}{\rangle }_{\phi }+\langle {\boldsymbol{F}}\cdot {\boldsymbol{v}}{\rangle }_{\phi }-\langle {Q}_{\mathrm{cool}}\rangle }_{\phi }.\end{eqnarray} \tag{ 26 }$$

The second term on the left-hand side is the radial energy flux gradient, which is designated as En${}_{E+P}$. The first term on the right-hand side, which is the gravitational potential energy released during the accretion process, is En_pot. The second term on the right-hand side is the work done by the torque, En_T, while the cooling rate $\langle {Q}_{\mathrm{cool}}{\rangle }_{\phi }$ is En_cool. Thus, the above equation is

$$\begin{eqnarray}&&{\partial }_{t}\langle E{\rangle }_{\phi }+{\mathrm{En}}_{E+P}={\mathrm{En}}_{\mathrm{pot}}+{\mathrm{En}}_{T}-{\mathrm{En}}_{\mathrm{cool}},\end{eqnarray} \tag{ 27 }$$

which basically states that the energy change is due to the flux into the domain, the work done by forces (the gravitational force from the central object and other forces), and radiative cooling. After being normalized with ${GM}\dot{M}/{R}^{3}$, all these terms for case Gp2 are plotted in the left-hand panel of Figure 14. Figure 14 suggests that the work done by the torque is negligible at the inner disk, and the cooling rate equals the difference between the released potential energy and the energy flux gradient.

Zoom In Zoom Out Reset image size

After removing second-order terms, En${}_{E+P}$ can be expanded as

$$\begin{eqnarray}{\mathrm{En}}_{E+P} & = & \displaystyle \frac{1}{R}{\partial }_{R}[\langle R(E+P){v}_{R}{\rangle }_{\phi }]\\ & \sim & \displaystyle \frac{1}{R}{\partial }_{R}\left[\langle {{Rv}}_{R}\left(\displaystyle \frac{1}{2}{\rm{\Sigma }}{v}_{K}^{2}+{\rm{\Sigma }}{v}_{K}\delta {v}_{\phi }+\displaystyle \frac{\gamma }{\gamma -1}P\right){\rangle }_{\phi }\right]\\ & \sim & \displaystyle \frac{1}{R}{\partial }_{R}\left[R\left(\displaystyle \frac{1}{2}{v}_{K}^{2}\langle {\rm{\Sigma }}{v}_{R}{\rangle }_{\phi }+{v}_{K}\langle {\rm{\Sigma }}{v}_{R}\delta {v}_{\phi }{\rangle }_{\phi }\right.\right.\\ & & +\left.\left.\displaystyle \frac{\gamma }{\gamma -1}\langle {v}_{R}P{\rangle }_{\phi }\right)\right].\end{eqnarray} \tag{ 28 }$$

With $\dot{M}=-R\langle {\rm{\Sigma }}{v}_{R}{\rangle }_{\phi }$, the above equation becomes

$$\begin{eqnarray}{\mathrm{En}}_{E+P} & \sim & \displaystyle \frac{1}{R}{\partial }_{R}\left[-\displaystyle \frac{1}{2}{v}_{K}^{2}\dot{M}+{{Rv}}_{K}\langle {\rm{\Sigma }}{v}_{R}\delta {v}_{\phi }{\rangle }_{\phi }\right.\\ & & +\left.\displaystyle \frac{\gamma }{\gamma -1}R\langle {v}_{R}P{\rangle }_{\phi }\right].\end{eqnarray} \tag{ 29 }$$

The three terms on the right-hand side are the gradient of kinetic energy, the Reynolds stress, and the pressure. Thus, we designate the three terms as En_kic, En_Rey, and En_pre, and we have

$$\begin{eqnarray}&&{\mathrm{En}}_{E+P}\sim {\mathrm{En}}_{\mathrm{kic}}+{\mathrm{En}}_{\mathrm{Rey}}+{\mathrm{En}}_{\mathrm{pre}}.\end{eqnarray} \tag{ 30 }$$

Note that ${\mathrm{En}}_{\mathrm{pre}}$ actually consists of two parts: one is the advection of the internal energy and the other is related to the work done by the pressure.

If we assume that the disk accretion rate is a constant, we can replace the Reynolds stress with the disk accretion rate by using Equation (17):

$$\begin{eqnarray}{\mathrm{En}}_{\mathrm{Rey}} & = & \displaystyle \frac{1}{R}{\partial }_{R}({{Rv}}_{K}\langle {\rm{\Sigma }}{v}_{R}\delta {v}_{\phi }{\rangle }_{\phi })\\ & = & \displaystyle \frac{1}{R}{\partial }_{R}[{v}_{K}^{2}\dot{M}-{v}_{K}\displaystyle \frac{C}{R}+\displaystyle \frac{{v}_{K}}{R}\displaystyle \int R\langle {\boldsymbol{R}}\times {{\boldsymbol{F}}}_{\mathrm{ext}}{\rangle }_{\phi }{dR}]\\ & = & {\mathrm{En}}_{\dot{M}}+{\mathrm{En}}_{\mathrm{Const}},\end{eqnarray} \tag{ 31 }$$

where we have further neglected the external force term and designate the first and second terms as En${}_{\dot{M}}$ and En_Const. Thus, finally we have

$$\begin{eqnarray}&&{\mathrm{En}}_{E+P}={\mathrm{En}}_{\mathrm{kic}}+{\mathrm{En}}_{\dot{M}}+{\mathrm{En}}_{\mathrm{Const}}+{\mathrm{En}}_{\mathrm{pre}}.\end{eqnarray} \tag{ 32 }$$

If we assume ${\phi }_{p}=-{v}_{K}^{2}$ and ${v}_{K}=\sqrt{{GM}/R}$ and normalize Equation (26) with ${GM}\dot{M}/{R}^{3}$ (designating the corresponding quantities with a hat symbol), we have

$$\begin{eqnarray*}&&\overline{{\mathrm{En}}_{\mathrm{cool}}}+\overline{{\mathrm{En}}_{\mathrm{kic}}}+\overline{{\mathrm{En}}_{\dot{M}}}+\overline{{\mathrm{En}}_{\mathrm{Const}}}+\overline{{\mathrm{En}}_{\mathrm{pre}}}=\overline{{\mathrm{En}}_{\mathrm{pot}}}+\overline{{\mathrm{En}}_{T}},\end{eqnarray*}$$

and thus

$$\begin{eqnarray}&&\overline{\langle {Q}_{\mathrm{cool}}{\rangle }_{\phi }}+\displaystyle \frac{1}{2}+\left(-1+\displaystyle \frac{3}{2}\displaystyle \frac{C}{\dot{M}{({GMR})}^{1/2}}\right)\\ &&\qquad +\displaystyle \frac{{R}^{2}}{{GM}\dot{M}}{\partial }_{R}\left(\displaystyle \frac{\gamma }{\gamma -1}R\langle {v}_{R}P{\rangle }_{\phi }\right)=1.\end{eqnarray} \tag{ 33 }$$

If we assume zero stress at R_in and neglect the pressure term, as in the thin-disk approximation, we have derived $C=\dot{M}{R}_{\mathrm{in}}{v}_{K,\mathrm{in}}$ in Section 3.2. Thus

$$\begin{eqnarray}&&\overline{\langle {Q}_{\mathrm{cool}}{\rangle }_{\phi }}=\displaystyle \frac{3}{2}-\displaystyle \frac{3}{2}{\left(\displaystyle \frac{{R}_{\mathrm{in}}}{R}\right)}^{1/2}).\end{eqnarray} \tag{ 34 }$$

In the total cooling rate, 1 comes from the release of the gravitational energy, and 1/2–3/2 ${({R}_{\mathrm{in}}/R)}^{1/2}$ comes from the radial energy transport. Except for the very inner disk, the radial energy transport normally heats up the disk. The term −3/2 ${({R}_{\mathrm{in}}/R)}^{1/2}$ in the radial transported energy is from the zero-torque boundary condition (Shakura & Sunyaev 1973).

To compare all the energy terms derived directly from our simulations with the corresponding terms in the thin-disk approximation, we summarize all the normalized energy terms in the thin-disk theory:

$$\begin{eqnarray}&&\overline{{\mathrm{En}}_{\mathrm{cool}}}=\overline{{\mathrm{En}}_{\mathrm{pot}}}+\overline{{\mathrm{En}}_{T}}-\overline{{\mathrm{En}}_{E+P}}\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&\overline{{\mathrm{En}}_{E+P}}=\overline{{\mathrm{En}}_{\mathrm{kic}}}+\overline{{\mathrm{En}}_{\dot{M}}}+\overline{{\mathrm{En}}_{\mathrm{Const}}}+\overline{{\mathrm{En}}_{\mathrm{pre}}}\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&\overline{{\mathrm{En}}_{\mathrm{pot}}}=1\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&\overline{{\mathrm{En}}_{T}}=0\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&\overline{{\mathrm{En}}_{\mathrm{kic}}}=\displaystyle \frac{1}{2}\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&\overline{{\mathrm{En}}_{\mathrm{pre}}}=0\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&\overline{{\mathrm{En}}_{\dot{M}}}=-1\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&\overline{{\mathrm{En}}_{\mathrm{Const}}}=\displaystyle \frac{3}{2}{\left(\displaystyle \frac{{R}_{\mathrm{in}}}{R}\right)}^{1/2}.\end{eqnarray} \tag{ 42 }$$

Thus,

$$\begin{eqnarray}&&\overline{{\mathrm{En}}_{E+P}}=-\displaystyle \frac{1}{2}+\displaystyle \frac{3}{2}{\left(\displaystyle \frac{{R}_{\mathrm{in}}}{R}\right)}^{1/2}\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&\overline{{\mathrm{En}}_{\mathrm{cool}}}=\displaystyle \frac{3}{2}-\displaystyle \frac{3}{2}{\left(\displaystyle \frac{{R}_{\mathrm{in}}}{R}\right)}^{1/2}.\end{eqnarray} \tag{ 44 }$$

and, eventually,

$$\begin{eqnarray}&&\sigma {T}_{\mathrm{eff}}^{4}=\displaystyle \frac{3{{GM}}_{p}\dot{M}}{8\pi {R}^{3}}\left(1-{\left(\displaystyle \frac{{R}_{\mathrm{in}}}{R}\right)}^{1/2}\right).\end{eqnarray} \tag{ 45 }$$

The normalized quantities in the thin-disk approximation are plotted in the right-hand panel of Figure 14.

By comparing the left- and right-hand panels in Figure 14 where all quantities have been normalized by 8GM_J $\dot{M}$/R³, we can see that the cooling rate in the simulation is smaller than the cooling rate derived from the thin-disk approximation, especially at the outer disk. The ratio between the cooling rate from simulations and that from the thin-disk approximation is shown in Figure 15. Clearly, the cooling rate is only half of that predicted in the thin-disk theory at $10{R}_{J}\lt R\lt 100{R}_{J}$. Based on the relation in Equation (35), this lower cooling rate in real simulations is caused by the higher value of $\overline{{\mathrm{En}}_{E+P}}$. A higher $\overline{{\mathrm{En}}_{E+P}}$ value means that the radial energy transport depletes more energy from the disk.

Zoom In Zoom Out Reset image size

To understand why $\overline{{\mathrm{En}}_{E+P}}$ is higher in simulations, we plot different components of $\overline{{\mathrm{En}}_{E+P}}$ in Figure 16. By comparing the real simulations (left panel) with the thin-disk approximation (right panel), we can see that the pressure term (En_pre), which represents the radial advection of both the internal energy and the pressure, is nonnegligible in the simulation. It is negative at the inner disk within 4 R_J and positive at the outer disk. Thus, the inner disk loses energy, and the outer disk gains energy due to the nonzero En_pre term. This nonnegligible radial transport of energy makes the CPD more like a "slim disk" (Abramowicz et al. 1988) than a thin disk.

Zoom In Zoom Out Reset image size

We also notice that $\overline{{\mathrm{En}}_{\mathrm{pre}}}$ almost balances the Reynolds stress term ($\overline{{\mathrm{En}}_{\mathrm{Rey}}}$), so the total energy flux gradient term ($\overline{{\mathrm{En}}_{E+P}}$) is almost equal to the kinetic energy term ($\overline{{\mathrm{En}}_{\mathrm{kic}}}$). Thus, empirically, we can approximate $\overline{{\mathrm{En}}_{E+P}}$ as 1/2 instead of $-1/2+3/2{({R}_{\mathrm{in}}/R)}^{1/2}$ in the thin-disk theory. The new cooling rate $\overline{{\mathrm{En}}_{\mathrm{cool}}}$ is thus ∼1/2. If we relate the cooling rate with the effective temperature (${Q}_{c}=\sigma {T}_{\mathrm{eff}}^{4}$), we can derive the effective temperature in the CPD as

$$\begin{eqnarray}&&\sigma {T}_{\mathrm{eff}}^{4}\sim \displaystyle \frac{{GM}\dot{M}}{8\pi {R}^{3}},\end{eqnarray} \tag{ 46 }$$

which is significantly different from the temperature in the thin-disk theory (Equation (45)). Thus, compared with the thin disk, the simulated disk is hotter at the inner disk and cooler at the outer disk. We want to emphasize that Equation (46), which assumes $\overline{{\mathrm{En}}_{\mathrm{cool}}}\,=$ 1/2 is an empirical approximation, slightly overestimates the inner disk temperature compared with real simulations, as shown in Figure 14, where $\overline{{\mathrm{En}}_{\mathrm{cool}}}$ is smaller than 1/2 at the inner disk.

Overall, our simulations show temperature profiles different from the thin-disk theory. This difference is not caused by the boundary condition since the zero-stress boundary condition also stands in our simulations, as discussed at the end of Section 3.2. Instead, the difference is caused by the radial advection of energy since H/R is large and the thin-disk approximation is not valid anymore in the CPD.