Modeling of deep gaps created by giant planets in protoplanetary discs

A giant planet embedded in a protoplanetary disk creates a gap. This process is important for both theory and observations. Using results of a survey for a wide parameter range with two-dimensional hydrodynamic simulations, we constructed an empirical formula for the gap structure (i.e., the radial surface density distribution), which can reproduce the gap width and depth obtained by two-dimensional simulations. This formula enables us to judge whether an observed gap is likely to be caused by an embedded planet or not. The propagation of waves launched by the planet is closely connected to the gap structure. It makes the gap wider and shallower as compared with the case where an instantaneous wave damping is assumed. The hydrodynamic simulations shows that the waves do not decay immediately at the launching point of waves, even when the planet is as massive as Jupiter. Based on the results of hydrodynamic simulations, we also obtained an empirical model of wave propagation and damping for the cases of deep gaps. The one-dimensional gap model with our wave propagation model is able to well reproduce the gap structures in hydrodynamic simulations. In the case of a Jupiter-mass planet, we also found that the waves with smaller wavenumber (e.g., $m=2$) are excited and transport the angular momentum to the location far away from the planet. The wave with $m=2$ is closely related with a secondary wave launched by the site opposite from the planet.


Introduction
During the past decade, more than a thousand extrasolar planets have been discovered. A survey of extrasolar planets has revealed the diversity of giant planets outside the solar system (e.g., Burke et al. 2014). Giant planets are born in protoplanetary disks, due to core accretion (e.g., Mizuno 1980;Kanagawa * & Fujimoto 2013) or to the gravitational instability of the gaseous disk (e.g., Cameron 1978;Zhu et al. 2012a). Once formed, they undergo orbital migration, and gas accretion results in growth their mass. The diversity of extrasolar planets is closely connected to such processes (e.g., Mordasini et al. 2012;Ida et al. 2013).
Planet forming regions in protoplanetary disks can now be directly imaged, such as by the Atacama Long Millimeter/Submillimeter Array (ALMA) and by eight-metre class optical and/or near-infrared telescopes. High-resolution images have revealed the presence of complex morphological structures in disks; these structures include spirals (e.g., Muto et al. 2012;Grady et al. 2013;Christiaens et al. 2014;Benisty et al. 2015;Currie et al. 2015;Akiyama et al. 2016) and gaps (e.g., Osorio et al. 2014; ALMA Partnership et al. 2015;Momose et al. 2015;Akiyama et al. 2015;Nomura et al. 2016;Yen et al. 2016;Tsukagoshi et al. 2016). Direct images are now finding possible signatures of planets being formed in disks (e.g., Sallum et al. 2015). To understand the origins of the disk structures and their possible connection to the formation of planets, it is important to construct appropriate quantitative models.
A planet interacts gravitationally with the gas in the surrounding disk, and as a result, the planet excites density waves (spirals). If the planet is sufficiently massive, it also creates gap structures (e.g., Lin & Papaloizou 1979;Goldreich & Tremaine 1980;Artymowicz & Lubow 1994;Kley 1999;Crida et al. 2006). Recently there have been a number of studies on the quantitative relationship between the mass of a planet and the gap structure (i.e., depth and width) (Duffell & MacFadyen 2013;Fung et al. 2014;Kanagawa et al. 2015a;Duffell & Chiang 2015;Fung & Chiang 2016;Kanagawa et al. 2016), and the application to actual observations has been discussed (Kanagawa et al. 2015a;Momose et al. 2015;van der Marel et al. 2016;Nomura et al. 2016;Tsukagoshi et al. 2016). However, there is still room for improvement of the current models. Kanagawa et al. 2016 presented an empirical formula for the gap width: they defined it to be the location at which the gap surface density is depleted by 50% from its original value. The overall structure of the gap created by the planet should be further investigated.
The gap structure is closely associated with the damping of spiral density waves excited by a planet, because the planet exchanges the angular momentum to the disk via the waves (e.g., Takeuchi et al. 1996;Goodman & Rafikov 2001;Rafikov 2002;Dong et al. 2011). Duffell 2015 constructed analytical models of the gap structure by using the wave damping model of Goodman & Rafikov 2001, which is particularly useful in the case of shallow gaps (created by relatively low-mass planets). However, the wave propagation in the case where the deep gap is formed by a high-mass planet is still poorly understood. The propagation of waves induced by the high-mass planet would be qualitatively different from that in the case of the low-mass planet, as implied by the parametrized study by Kanagawa et al. 2015b (hereafter, K15).
In this paper, we extend our previous model of the gap width and depth, and present a complete model of the gap shape; this is based on a number of two-dimensional long-term hydrody-namic simulations. We also provide a model of the wave propagation indicated by the hydrodynamic simulations, when the deep gap is formed. In Section 2, we briefly summarize the setup of the numerical simulations we performed. In Section 3, we obtain an empirical formula for the gap structure, showing the results of hydrodynamic simulations. From this empirical formula of the gap structure, we can build an empirical model of propagation of the density waves. In Section 4, we obtain the empirical model of the wave propagation. Adopting this model of the wave propagation, we provide a semianalytical radial one-dimensional model of the gap structure which is able to reproduce the gap structures given by the two-dimensional simulations. In Section 5, we discuss an observational application of our formula. We also discuss how the wave excitation and propagation are changed along with an increase of a planet mass, showing the results of the two-dimensional hydrodynamic simulations. A discussion about hydrodynamical instabilities such as the Rayleigh instability and Rossby wave instability is included in this section. Section 6 contains a summary and discussion. Kanagawa et al. 2016 performed a number of the twodimensional simulations over a wide range of parameter space (planet mass, disk scale height and viscosity), and derived a relationship between the gap width and the planet mass. In this paper, we extend the model by further analyzing of the numerical simulations. In this section, we briefly review the setup of hydrodynamic simulations.

Summary of setup of numerical simulations
We numerically calculated the gap formation processes in a protoplanetary disk in the presence of a planet. We assumed a geometrically thin and non-self-gravitating disk. We adopted a two-dimensional cylindrical coordinate system (R, φ), and the origin was located at the position of the central star. We adopted a simple locally isothermal equation of state, and the distribution of the temperature is independent of time. The mass of the central star is also constant during the simulations. Using FARGO (Masset 2000), which is widely used in the study of disk-planet interaction (e.g., Crida & Morbidelli 2007;Baruteau et al. 2011;Zhu et al. 2011), we solved the equations of continuity and motion of disk gas. In this paper, we adopt α-prescription of Shakura & Sunyaev 1973, and then the kinetic viscosity is written as ν = αh 2 Ω −1 k , where h and Ω k are the disk scale height and the Keplerian angular velocity, respectively.
The gap is opened by the disk-planet interaction and is closed by the viscous diffusion. In steady state, the viscous angular momentum flux is balanced by the planetary torque. The timescale of the gap opening would be scaled by the viscous timescale (e.g., Lynden-Bell & Pringle 1974). As an empirical formula, the gap width (∆gap), which is defined as a width of a region where the surface density is smaller than 0.5Σ0 in steady state, where Σ0 is the surface density outside the gap, is obtained as follows (Kanagawa et al. 2016): where The timescale of the gap opening (tvis) is Using equation (1), this timescale is obtained as To obtain the steady state, we have to calculate the evolution until t ∼ tvis (see, Appendix 1). Note that for nominal parameters, this timescale is usually shorter than or comparable with the migration timescale of the type II (∼ R 2 p /ν) or disk lifetime (∼ 1Myr). Hence, a full-width gap can be observed.

Empirical formula for gap structures
In this section, we obtain an empirical formula for radial structure of the planet-induced gap in steady state. Kanagawa et al. 2016 have found the scaling relation of the gap width, which is defined by a radial width of a region where the surface density is smaller than a density threshold of Σ th = 0.5Σ0. First we here show that this scaling relation can be extended even when the density threshold is different from 0.5Σ0. As illustrated in the left panel of Figure 1, we measure the gap widths for Σ th = 0.1Σ0, 0.3Σ0, 0.5Σ0, and 0.8Σ0. The right panel of Figure 1 shows the gap width with each Σ th in terms of K ′ . As can be seen in the figure, the gap widths defined by Σ th = 0.1Σ0 and 0.3Σ0 are proportional to K ′1/4 , as was seen with Σ th = 0.5Σ0. Note that we confirmed that the widths shown in the figure reach to the steady values at t ≃ tvis given by equation (3) (see, Appendix 1). Although the gap widths defined by Σ th = 0.8Σ0 deviate slightly from the line defined by ∝ K ′1/4 for K ′ < ∼ 1, they are roughly proportional to K ′1/4 . We obtain an scaling relation including that shown by Kanagawa et al. 2016 as, This scaling relation agrees with the results given by previous hydrodynamic simulations (see, Appendix 2). Note that the above relationship corresponds to equation (1) when Σ th /Σ0 = 0.5. Equation (5) is in good agreement with the gap widths shown in Figure 1. We can derive a formula of a radial distribution of the gas density from equation (5), setting Σ th = Σ(R) and ∆gap = 2|R − Rp| as, with where ∆R1 and ∆R2 are given by The surface density of the gap bottom Σmin is obtained by (e.g., Duffell & MacFadyen 2013;Fung et al. 2014;Kanagawa et al. 2015b) where We illustrate the surface density given by the empirical formula of equation ( (6) is still able to fit the gap structure of the simulations, but for (α, hp/Rp) = (10 −3 , 1/20), the gap obtained by the simulation is wider than that given by the equation. Equation (6) does not give a good estimate when the gap is shallow, that is, Σmin > ∼ 0.5Σ0. Crida et al. 2006 provided a semianalytical model for deep gaps. They introduced "pressure torque" and considered the balance between the pressure torque, the planetary torque and the viscous angular momentum flux. In Figure 3, we compare gap radial structures given by equation (6) and Crida et al. 2006. We obtain Crida's model in the figure by solving equation (14) of Crida et al. 2006. In the calculation of Crida's model, we adopt 1.5 times larger Σ0 than ours, to set Σ outside the gap in Crida's model to be the same level as that of our model. Crida's models result in significantly shallower gaps compared to the gaps obtained in the two-dimensional simulations with our parameter set. Equation (6) nicely reproduces the gap depths in the case of the parameters presented in the bottom right panel of Figure 3. In the case with α = 10 −3 , hp/Rp = 1/20 and Mp/M * = 10 −3 (the bottom right case), the gap depth obtained by equation (6) is about 5 times shallower than that in the simulation. This discrepancy partially comes from fact that the gas is not rotating at Kepler velocity. We provide a semianalytical model that takes into account this effect in Section 4. For the gap width, both equation (6) and Crida's model reasonably agree with the two-dimensional simulations. However, when the planet mass is small, Crida's model gives narrower gaps than the simulation, as can be seen from the bottom left panel in Figure 3.

Model of wave propagation
A planet exchanges angular momentum with the surrounding disk gas through density waves. As pointed out by previous studies (e.g., Petrovich & Rafikov 2012;Kanagawa et al. 2015b), a propagation of the density waves is directly connected with the surface density distribution of the gap. The planet gives angular momentum to density waves. Then, the angular momentum from the planet is carried by the waves. Finally the angular momentum is deposited to the disk gas by dissipation of the waves (e.g., Takeuchi et al. 1996;Goodman & Rafikov 2001). In this section, considering the wave propagation which is suggested by our empirical formula of the surface density distributions of the gap obtained in the previous section, we present a one-dimensional model of the gap. From the conservation of the angular momentum, where v avg φ is the azimuthal averaged v φ and FJ (Rp) is the angular momentum flux at the radius of the orbit of the planet (Kanagawa et al. 2015b). The viscous angular momentum flux F vis J is given by equation (A6) in Appendix 3. Furthermore, the cumulative torque deposited to the disk, T d (R) is defined by equation (A11), which is equal to the difference between the cumulative torque excited by the planet, Tp(R), and the angular momentum flux due to the density waves, F wave J (R) which is defined by equation A5. In Figure 4, we show a schematic picture of the effect of the wave propagation. (see Kanagawa et al. 2015b and Appendix 3 for detail).
The surface density is given by equation (12). If assuming that the gap width is smaller than Rp, we neglect the terms which is proportional to a power of hp/Rp, as e.g., R = Rp and Ω = Ω k,p , where the suffix p indicates values at Rp in the following. The viscous angular momentum flux is given by F vis J ≃ 2πR 2 p νΣ(dΩ/dR). The derivative of the angular velocity is dΩ/dR ≃ −3Ω k,p /(2Rp) 1 − (h 2 p /3)(d 2 ln Σ/dR 2 ) , of which second term is not negligible since d 2 ln Σ/dR 2 ∼ 1/h 2 p in the gap region (Kanagawa et al. 2015b). Then, we rewrite equation (12) as, As shown in equation (13), the surface density distribution is determined if the function of T d is given. The function of T d implicitly depends on the surface density, because T d depends on the torque exerted by the planet which is determined by the surface density distribution. To obtain the distributions of the surface density and deposited torque, an iterative method is required (see, Kanagawa et al. 2015b).
The propagation of the waves launched by the planet has been investigated by e.g., Goodman & Rafikov 2001. Their model indicates that the waves launched by a sufficiently massive planet (Mp > M1 = (2/3)(hp/Rp) 3 M * ≃ 1 × 10 −4 M * ) decay due to the weak shock that develops close to the launching point of the density wave. Because of the shock, the angular momentum of the waves decreases as |R − Rp| 5/4 from the shock location. Using the wave propagation model of Goodman & Rafikov 2001, Duffell 2015 proposed an analytical model of the gap structure. His model is able to accurately reproduce the structures of the gap, if the gap is shallow. When the gap is deep, however, the model of Duffell 2015 predicts a much narrower gap than that obtained by numerical simulations, which implies that the angular momentum deposition from the waves is less effective than that predicted by the model of Goodman & Rafikov 2001. For the deep gap, Kanagawa et al. 2015b presented the simple model of wave propagation using two free parameters x d and w d , which represent the position at which waves are damped and the width of the deposition site, respectively, as where T ∞ p is the total excitation torque exerted by the planet (see Appendix 3 for detail). If we assume that the form of the torque deposit is given by equation (14), we can obtain the gap shape using equation (12). It should be noted, however, that equation (14) is a purely parametrized model for T d (R), and its validity should be checked with numerical simulations. In particular, the gap surface density distribution depends strongly on the choice of the two free parameter x d and w d and therefore, realistic values for these parameters should be determined using numerical simulations.
Motivated by the surface density distribution of equation (6), we empirically determine the parameters x d and w d as where ∆R1 and ∆R2 are defined by equations (8) and (9), respectively, and the sign of x d is positive when R > Rp and negative when R < Rp. Figure 5 shows the cumulative torques deposited onto the disk given by equation (12) and the twodimensional hydrodynamic simulations. The cumulative deposited torque given by equation (14) with equations (15) and (16) reasonably reproduces the results of the two-dimensional hydrodynamic simulations. We use the WKB formula (e.g., Ward 1986) to calculate the excitation torque. However, we should consider nonlinear effects since we deal with relatively massive planets creating deep gaps in this paper. According to Miyoshi et al. 1999, the nonlinear effect reduces the excitation torque as compared with that expected by the linear theory. In order to include the nonlinear effect, we introduce a reduction factor fNL(< 1) and multiply this factor by the torque density (see Appendix 3 for detail). As shown below, fNL = 0.4 better reproduces the gap depth obtained by the numerical hydrodynamic simulations. Since we consider cases where Mp > (hp/Rp) 3 , the factor fNL is always set to 0.4.

Radial gap structures
In Figure 6, we show the radial structure of a gap, calculated by solving equation (13) with x d and w d given by equations (15) and (16). For comparison, we plot gap structures obtained by the two-dimensional hydrodynamic simulations and equation (6) in the figure. The solutions of equation (12) are in good agreement on the width and depth with the hydrodynamic simulations for various disk scale heights and viscosities, as long as the planet mass is relatively large. Since the deviation from Keplerian rotation is taken into account in the semianalytical model, the structure of gap bottom is smoother and the gap is slightly deeper than those given by equation (6). When Mp/M * = 10 −3 , α = 10 −3 , and hp/Rp = 1/25 (bottom right panel of Figure 6), the semianalytical model (equation 12) gives approximately three time shallower gap than that given by the two-dimensional simulation. We consider that the semianalytical model provides a reasonable fit to deep gaps created by a planet.
The left panel of Figure 7 shows the minimum surface density of the gap Σmin obtained by the one-dimensional model with fNL = 0.4 and 1.0. We also plot the minimum value of the azimuthal averaged surface density given by the twodimensional simulations. When the gap is shallow as K < 10 2 , the one-dimensional model (equation 13) gives us a similar Σmin as that given by the two-dimensional simulations, which are also consistent with equation (10), in both the cases of fNL = 0.4 and 1 (but the model with fNL = 1 looks more appropriate). As the gap is deeper in K > 10 2 , the one-dimensional   (10), and the solid and chain lines indicate Σmin, as given by the solution of equation (12) with fNL = 0.4 and 1.0, respectively. (Right) Surface density averaged throughout the gap bottom Σ ′ min , as given by Kanagawa et al. 2016 (circles). The crosses, squares, and triangles are the gap depths given by Fung et al. 2014, Varnière et al. 2004, and Duffell & Chiang 2015 model with fNL = 1 gives us a much smaller Σmin as compared with equation (10) (and 2D simulations). On the other hand, Σmin obtained by the one-dimensional model are consistent with these given by two-dimensional simulations and equation (10). For very deep gap with K > 10 3 , Σmin given by two-dimensional simulations is smaller than that predicted by the empirical formula of equation (10). Even in this case, the one-dimensional model with fNL = 0.4 partially fits Σmin given by the simulations. The assumption of fNL = 0.4 is reasonably acceptable in the cases with K < 10 4 . When the gap is very deep, Σmin very sensitively depends on the value of T ∞ p . More sophisticated model of wave excitation by a giant planet would be required to reproduce the depth of the very deep gaps with K > 10 4 .
The empirical formula of equation (10) predicts a larger value of Σmin than that given by the two-dimensional simulations for K > 10 3 , as shown above. However, equation (10) is still useful for Σ ′ min which is the surface density averaged throughout the gap bottom region 1 , as Fung et al. 2014 showed. In the right panel of Figure 7, we illustrate Σ ′ min given by the two-dimensional simulations. When the gap is not very deep (K < ∼ 10 3 ), both Σ ′ min and Σmin are reproduced by equation (10). Even when the gap is very deep with K > ∼ 10 3 , Σ ′ min is also reasonably consistent with equation (10)  around R = Rp ± 2RH , rather than Σmin = Σ(Rp). According to Tanigawa & Watanabe 2002, the gas accretion rate onto the planet depends on the surface density around the location that is twice the Hill radius apart from the planet. Hence, equation (10) may be useful to estimate the gas accretion rate, whereas it overestimates the minimum surface density of the gap when K > ∼ 10 3 .

Constraint from the gap radial structure
Here we would like to discuss observational applications of our empirical formula of the gap structure of equation (6). Recently, many gap structures have been discovered by ALMA; these have been found not only in dust (e.g., ALMA Partnership et al. 2015) but also in gas (e.g., Yen et al. 2016). Such gap structures in protoplanetary disks can be created by dust growth (Zhang et al. 2015), sintering (Okuzumi et al. 2016), the effects of MRI (Flock et al. 2015), secular gravitational instability due to gasdust friction (Takahashi & Inutsuka 2016), and disk-planet interaction. It is difficult to distinguish the origin from observations. However, using the relationship of the gap radial structure that is shown above, we may be able to assess whether the observed gap is due to the disk-planet interaction. Using equation (5), we can determine the ratio of the gap width for any two arbitrary surface densities Σ th,a and Σ th,b as As can be seen from the above equation, the ratio of the gap widths depends only on the surface densities defined by the widths. In Figure 8, we show two ratios of gap widths, defined by Σ th = 0.1Σ0 and 0.5Σ0, and by Σ th = 0.3Σ0 and 0.5Σ0. Using equation (8), we obtain that the ratios are 0.51 and 0.76, respectively. As can be seen in the figure, the ratios obtained from the simulations are reasonably consistent with those obtained with equation (17).

Gap depth-width relation
As shown in Kanagawa et al. 2016, there is a relationship between the depth and width of a gap. Combining equations (5) and (10), we can derive this relationship in terms of the locations of Σ th , as follows: If the gap is created by a planet, its width, as measured by the above-surface density, should satisfy equation (18). When the gap structure is completely resolved and the disk aspect ratio is precisely estimated, we can use the gap widths measured by the different surface density at the gap edge and equation (18) to strictly judge whether the gap was created by the planet.

Applicability of the model
We now discuss some considerations when equations (17) and (18) (17) and (18) should be used for observations of disk gas. Alternatively, the dust-gas coupling depends on the gas surface density, as well as on the size of the dust particles. If the gas density is sufficiently large, the gas and dust particles will be well mixed. In this case, equations (17) and (18) would provide a good estimate.
Second, we assume that the gap structure is in steady state. As shown in Appendix 1, the timescale required for the gap to be in steady state can be roughly estimated as tvis ∼ 0.1 Myr. If an observed gap is younger than this timescale, the gap will be narrower than it would be in steady state. Such young gaps cannot be estimated by equations (17) and (18). We note that equations (17) and (18) should be used for relatively deep gaps. As shown in Figure 1, for a shallow gap (K ′ < 1), the actual gap width measured for a larger surface density is slightly wider than that estimated by equation (18). Because of this, in this case, equation (18) would estimate the disk scale height to be about 1.5 times the actual value. Moreover, as also discussed in Kanagawa et al. 2016, the observational uncertainties of the orbital radius of the planet and aspect ratio should be taken into account.
Finally, we briefly make comments about our assumption of two-dimensional disks and spatially constant kinematic viscosity ν. Fung & Chiang 2016 performed three-dimensional simulations of disk-planet interaction in the case when a deep gap is induced by a planet. The gap profile (i.e., the depth and the width) is not very different from the results of two-dimensional calculations. Hence, our models may be valid even when threedimensional effects are properly taken into account. The assumption of constant viscosity may not be always satisfied. Zhu et al. 2013 performed ideal MHD simulations of the interaction between a low-mass planet and a disk and showed that the effective viscosity α within the gap is approximately twice as large as that outside the gap region. Therefore, the gap may be shallower than our model. However, since the dependence of α on the gap depth is not very strong (see, equation 10), this effect may not significantly influence our results.

Planetary torque and angular momentum deposition in the disk
In Section 4, we present the empirical model of wave propagation and obtain the semianalytical one-dimentional model of gap structure using this wave propagation model. In

Excitation and damping of waves in low-m modes
To further examine the wave propagation when a deep gap is present, we investigate the wave resonances associated with small wavenumber using the results of two-dimensional simulations. We consider the Fourier components for Σ, vR, v φ , and Ψ which are given by where f is Σ, vR, v φ , or Ψ. For convenience, we will let the subscript m indicate the m-th Fourier component. The torque density exerted on the m-th resonance is given by (Goldreich & Tremaine 1980) and the integrated torque exerted on the m-th resonance is given by The angular momentum flux doe to the m-th mode of the wave is The cumulative torque deposited by the m-th mode T d,m can be written as In Figure 10, we illustrate the cumulative torque exerted by the planet and that deposited on the disk for low-m modes. The parameters for the disk and the planet are the same as those in Figure 9. We plot the contributions from m = 1 -4 individually and the sum of contributions from modes larger than 4. In the case of Mp/M * = 10 −4 (left panel of the figure), the planetary torque is mainly exerted by the resonances with m > 4 since the wave excitation is in the linear regime (Goldreich & Tremaine 1980). The cumulative deposited torques are also dominated by the contribution from m > 4 waves.
In the cases where the plant mass is large as in the middle and right panels of Figure 10 1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9  variation of Tp at large distances pointed out in previous subsection (and shown in Figure 6) is originated from the contributions from m = 1, 2 resonances. Note that the contribution of m = 1 resonance is related with the indirect term of the gravitational potential (which is connected with the gravitational interaction between the planet and the central star), which just oscillates around zero at larger distances from the planet, and the contribution of m = 1 does not significantly influence the angular momentum flux of the waves. As in the planetary torque, the contributions of the cumulative deposited torque from the small m modes become larger as the planet mass increases. When m = 2, in particular, even at radii very far away from the planet (R/Rp ∼ 2), the cumulative deposited torque slowly increases, in the case of Mp/M * = 10 −3 . These behaviors of Tp and T d indicate that the waves with small m modes (except m = 1) significantly excites and carry the angular momentum to the large radii from the planet, with the increase of the planet mass. As the gap becomes deeper and wider, the contributions on the wave excitation from large m resonances (near the planet) becomes smaller because the surface density at the resonance decreases. In consequence, the wave excitation on smaller m resonances becomes significant, as discussed by Juhász et al. 2015. Moreover, Lee 2016 has shown that the excitation of waves with small azimuthal wavenumber can be enhanced due to nonlinear effects. The significantly excitation of waves with small m would be explained by the gap opening and the nonlinear effects.
We should note that the wave with m = 2 deposit the angular momentum outside the gap with Σ > Σ0 (R/Rp > 1.8). This indicates that wider gap is formed by the deposition of the angular momentum of m = 2 wave, whereas opening such the wide gap requires very long time that is comparable with the disk life time (see equation 3).
The development of the m = 1, 2 contributions significantly changes the morphology of the waves in addition to the transport of angular momentum. In Figure 11, we show the surface density produced by the m = 1, 2 components and the sum of the m > 2 components, for the same case as shown in Figure 10. The surface density produced by the m-th component is where Σm is the Fourier component of the surface density. The waves that are far from the planet are created by the contributions from the m = 1 and m = 2 resonances. In particular, the secondary wave that is launched from the site opposite the planet is mainly composed of the m = 1 and m = 2 components. Because the contribution from the m = 2 component is strong, the secondary wave originates at a site opposite from the planet. Fung & Dong 2015 have shown that the azimuthal separation of the primary and secondary waves depends on the planet mass. For a large planet, as shown in Figure 11, this separation is close to 180 • , which is consistent with our results. For a smaller planet, the separation is smaller than 180 • . In this case, the contributions from other lower modes (e.g., m = 3, 4), are not much smaller than these from m = 1, 2 (see Figure 10). Because of these contributions, the launching point of the secondary wave would move toward the planet. It is worth noting that the secondary wave interacts with the planet. As seen in Figure 11, the secondary wave passes through the location at which φ = φp around R/Rp = 1.4. Near this location, the planetary torque due to the m = 2 resonance increases and then decreases, as shown in the top panel of Figure 10; this is a result of the interaction between the planet and the secondary wave. The cumulative deposited torque of m = 2 also increases and decreases when the secondary wave passes through the location at which φ = φp; this is also a consequence of the interaction between the planet and the secondary wave, and especially of the m = 2 component.

Effect of instability for gap structures
Previous studies (e.g., Li et al. 2000;Tanigawa & Ikoma 2007;Lin 2012;Lin 2014;Ono et al. 2014;Kanagawa et al. 2015b;Ono et al. 2016) pointed out a possibility that hydrodynamic instabilities (i.e., Rayleigh Instability, Rossby Wave Instability) occur at the edge of the gap induced by the planet. The condition for the onset of the RI is given by d(R 2 Ω)/dR < 0 (Chandrasekhar 1961). In our parameter range, however, this condition is not satisfied because the gap structure is less steep. For a very deep gap with a massive planet as Mp/M * ≃ 4 × 10 −3 , the RI may be violated as shown by Fung & Chiang 2016. According to Lovelace et al. 1999, the RWI can occur when the function of ωzS 2/Γ = ωzc 2 s (when locally isothermal), where ωz = (∇ × V )z/Σ is the potential vorticity and cs is sound speed, respectively, has an extreme value (a maximum or min-imum), whereas Ono et al. 2016 have shown that Lovelace's condition is not a sufficient one. The left panel of Figure 12 shows radial distributions of azimuthal averaged values of ωzc 2 s for Mp/M * = 10 −4 , 5 × 10 −4 , and 10 −3 when hp/Rp = 1/20 and α = 10 −3 . When Mp/M * = 10 −4 and 5 × 10 −4 , there is no extreme around the gap edge. When Mp/M * = 10 −3 , we can find a local minimum and maximum around R/Rp = 1.2. Even in this case, the two-dimensional distribution of ωzc 2 s is almost axis-symmetric (see the right panel of the figure) and there is no vortex-like structure, it does not seem that the RWI affects the gap structure.
The RWI can occur at the edge of the planet-induce gap in the case of the low viscosity (e.g., Yu et al. 2010;Fu et al. 2014;,see also Appendix 4). The RWI may constraint the gap structure in this case, as stated by Hallam & Paardekooper 2017. However, with viscosity larger than α ∼ 10 −3 -10 −4 , the RWI is stabilized by the viscosity during the long-term evolution (e.g., Lin 2014;Fu et al. 2014;. In fact, there is no clear evidence which the RWI affects the gap structure in Figure 12. In this case, the effect of the RWI on the gap structure may not be essential.

Summary
Using the results of the survey of the hydrodynamic simulations (34 runs) over the wide parameter space, we derived a quantitative relationship between a planet and the gap radial structure. We also obtain an empirical model of the wave propagation, from the results of the survey. Using this empirical model of wave propagation to a semianalytical model of the gap structure provided by Kanagawa et al. 2015b, we can obtain the gap structures which are in very good agreement with these obtained by the two-dimensional simulations. Our results can be summarized as follows: 1. We extended the scaling relation of the gap width given by Kanagawa et al. 2016. From this scaling relation, we derived an empirical formula for the surface density distributions of the gap (equation 6), which accurately reproduces the results of hydrodynamic simulations; see Figure 2. 2. Our model puts a constraint on the origin of the observed gap structure, as discussed in Section 5.1. Using equations (17) and (18), we can judge whether an observed gap is induced by the planet; see Figure 8. 3. We confirmed the validity of the empirical model for wave propagation adopted in Kanagawa et al. 2015b (equation 14) with the parameters given by equations (15) and (16). This model describes wave propagation that is consistent with the results of the hydrodynamic simulations; see Figure 5. Using this wave propagation model, we can reproduce the gap structure obtained by the two-dimensional hydrodynamic simulations by the one-dimensional model; see Figure 6. 4. Our model of wave propagation indicates that the waves excited by the larger planet carry the angular momentum at larger distances from the planet. As the planet mass is larger, waves with smaller azimuthal wavenumber (e.g.,m = 2) are excited strongly. The angular momentum is transported to the distant locations from the planet, even when the planet mass is larger as Jupiter; see Figure 10. 5. The development of modes with small azimuthal wavenumber changes the morphology of spiral waves. When the planet mass is sufficiently large, the m = 2 mode creates a secondary wave launched from the site opposite from the planet (see Figure 11), in addition to a primary wave originated from the location of the planet. The secondary wave would be closely related with the transport of angular momentum.
As the planet mass increases, the nonlinear effects become significant in excitation and propagation of waves, and thus the mechanism of the gap formation is different from that when the planet is small in linear theory. However, these theoretical mechanisms are not yet completely understood. To understand the mechanism of the gap formation and the morphology of the density waves induced by a giant planet, it will be necessary to further investigate the wave excitation and propagation when there are deep gaps.

Appendix 1 Time variations of the gap depth and width
Here we consider time-variation of the depth and width of gaps. Figure 13 shows the time variation of the gap width and depth (left panel) and a snapshot of the azimuthally averaged surface density (right panel), for Mp/M * = 10 −3 , α = 10 −3 , and hp/Rp = 1/25. We will first consider the time variation of the gap width. As can be seen in the left panel of Figure 13, a narrow gap, ∆gap = 0.4Rp, is formed at t = 100 planetary orbits. The gap width gradually widens with time. Finally, the gap width reaches 0.79Rp at t ≃ tvis = 1.5 × 10 4 planetary orbits. After that, the gap width is almost saturated. The distribution of the azimuthally averaged surface density becomes wider with time as shown in the right panel. Note that the evolution of the surface density is slightly asymmetric with respect to the planet; this may be due to the influence of the inner boundary. However, the gap structure is almost symmetric, as can be seen in the right panel.
In Figure 14, we show the time variations of the gap width and depth in the cases with different viscosity. The gap width becomes saturated slower with smaller viscosity, as expected by equation (4). In either cases in the figure, after t ≃ tvis, the gap widths are saturated. The saturated time for the gap depth also depends on the viscosity. If α = 10 −4 (the right panel of Figure 14), the gap depth is saturated after 10 4 orbits, which is about 10 times larger than that in α = 10 −3 case (the middle panel). Hence, if the viscosity is very small, a long-term simulation must be required to obtain gaps in steady state.

Appendix 2 Comparisons of the gap widths with previous works
We compared our result to those of hydrodynamic simulations in previous studies. Figure 15 shows the gap widths given by our simulations and those of Varnière et al. 2004 andDuffell &MacFadyen 2013. Because these studies used Σ th = 1/3Σ0 to define the gap width 2 , we plot the gap widths as measured by Σ th = 1/3Σ0. For K ′ < ∼ 10, the gap widths obtained by previous studies increases as K ′1/4 , which is consistent with the results of equation (5). As K ′ continues to increase, the gap widths given by previous studies have a larger scatter respect our results. In particular, the widths given by Duffell & MacFadyen 2013 are significantly narrower than our results. They mainly investigated situations with very low viscosity, that is, α < ∼ 10 −4 . We here remind the reader of the time variations of the gap width described in Appendix 1. The time that the gap reaches the saturated width is estimated by the viscous timescale, as in 2 Strictly speaking, Varnière et al. 2004  equation (3). When α = 10 −4 , the time scale is about 10 5 orbits, although previous studies estimated this to be several thousand planetary orbits. The variance in the gap width seen in the previous studies can be explained by their short computational times. Note that when the viscosity is very low, because the time scale of equation (3) is very long, the gap cannot reach the saturated width within the lifetime of the disk. In this case, equation (5) underestimates the planet mass.
Duffell & Chiang 2015 investigated the gap widths in the viscous case (α ∼ 10 −3 ). They have measured the gap width using Σ th = 0.1Σ0 and derived the following empirical formula: ∆gap/Rp = 0.5(hp/Rp) 0.22 K ′0.22 . Their empirical formula is quite similar to equation (5) if hp/Rp = 0.04; this is because they have primarily considered cases for which hp/Rp ≃ 0.04.

A.3.1 Wave propagation in one-dimensional model
In order to connect the two-dimensional simulation with the one-dimensional model of K15, we obtain the expression of wave propagation which are given in two-dimensional simulations, for the semianalytical model. We used the basic equations for the two-dimensional hydrodynamic simulations to derive the basic equation for the one-dimensional model.
Combining the continuity equation and the equation of motion in two-dimensional disk, the equation of angular momentum becomes where f φ is the viscous force per unit area acting in the azimuthal direction, and the gravitational potential is Ψ. We assumed that background structure is axisymmetric and rotates with a constant angular velocity, while the perturbed structure (wave) is nonaxisymmetric and rotates with an angular velocity that depends on R and φ. Hence, we can decompose the azimuthal velocity into two parts, the background and the perturbation (Muto et al. 2010): Note that equation (A2) is not a linear approximation. The azimuthally averaged version of equation (A1) can be written as ∂RΣv φ ∂t + 1 2πR where FJ is the azimuthally averaged angular momentum flux, in which where the first term is the angular momentum flux induced by an axisymmetric advection, and the azimuthally averaged mass flux FM is 2πRΣvR. The angular momentum fluxes due to nonaxisymmetric advection (waves) F wave J and viscous diffusion F vis J are respectively defined by In Figure 16, we show F vis J and the first term of equation (A6) in the cases of Mp/M * = 5×10 −4 and 1×10 −3 , hp/Rp = 1/20 and α = 10 −3 . As shown in the figure, F vis