Three-dimensional Interaction between a Planet and an Isothermal Gaseous Disk. III. Locally Isothermal Cases

We performed linear calculations to determine the Type I planetary migration rate for 3D locally isothermal disks with radial temperature gradients. For 3D disks with radial temperature gradients, the linear wave equation has a divergent term of the third pole, which makes corotation a nonremoval singularity. We suppressed the divergence with the Landau prescription to obtain the wave solutions. Despite the singularity at corotation, we derived a definite torque on the planet because the divergent term amplifies the waves only in the neighborhood of corotation and has little effect on the planetary torque. Consequently, we derived the formulas for the total, Lindblad, and corotation torques for locally isothermal disks. The resulting torque term due to the disk temperature gradient agrees well with the results of previous 3D hydrodynamical simulations for locally isothermal disks. Our linear calculation also provides the 3D horseshoe torque, which is close to the results of previous 3D hydrodynamical simulations.


INTRODUCTION
Type I planetary migration is a crucial process in the planet formation theory.It is caused by the gravitational interaction between a gapless gas disk and an embedded low-mass planet.Initially, the Type I migration had been studied with two-dimensional linear calculations (e.g., Goldreich & Tremaine 1979, 1980;Lin & Papaloizou 1979;Ward 1986;Artymowicz 1993;Korycansky & Pollack 1993), and three-dimensional calculations were later carried out by Takeuchi & Miyama (1998) and Tanaka et al. (2002, hereafter referred to as Paper I).The Type I migration torque arises from Lindblad and corotation resonances.In locally isothermal disks with moderate density and temperature gradients, the torque due to Lindblad resonances generally overcomes that due to corotation resonances and causes inward planetary migration.In Paper I, a torque formula was Corresponding author: Hidekazu Tanaka hidekazu@astr.tohoku.ac.jp provided for three-dimensional isothermal disks with a surface density of σ ∝ r −α as Γ = −(1.36+ 0.54α)q 2 h p r p −2 where q is the mass ratio of the planet to the host star, r p and Ω p are the orbital radius and angular velocity of the planet, and h p and σ p are the disk scale height and surface density at r = r p .This torque formula gives a short inward migration time of ∼ 10 5 yrs for a Jovian solid core with 10M ⊕ at 5AU in a typical protoplanetary disk with a surface density of 100g cm −2 and a temperature of 100K.Rapid Type I migration is a major obstacle to the formation of solid cores of giant planets due to the much longer core accretion time of ∼ 10 7 yrs (e.g., Kokubo & Ida 2002).Many attempts have been made to accelerate the core accretion with a significant supply of pebbles or tiny planetesimals (e.g., Ormel & Klahr 2010;Lambrechts & Johansen 2012;Kobayashi & Tanaka 2021, 2023).Several three-dimensional hydro-dynamical simulations have shown good agreements with the migration arXiv:2404.12521v1[astro-ph.EP] 18 Apr 2024 rates given by this torque formula (e.g., D'Angelo et al. 2003;Bate et al. 2003;D'Angelo & Lubow 2010).These three-dimensional hydrodynamical simulations have also been done for locally isothermal disks with non-zero temperature gradients.D'Angelo & Lubow (2010) have systematically investigated the dependence of the torque on the disk temperature gradient as well as the dependence on the surface density gradient.Their results have been confirmed by Jiménez & Masset (2017).However, three-dimensional linear calculations of the Type I migration torque have not yet been performed for radially non-isothermal disks.
While linear calculations are useful in verifying hydrodynamical simulations, the three-dimensional linear calculation of the Type I migration has only been conducted for completely isothermal disks.We demonstrated in Paper I that the linear wave equation for three-dimensional density waves has a third pole at the corotation point for disks with non-zero temperature gradients.This indicates that corotation is an irregular singular point, making it challenging to find a 3D wave solution for disks with radial temperature gradients.
In this study, we attempt to obtain the Type I migration torque for 3D locally isothermal disks with radial temperature gradients, by suppressing the divergence at corotation in the linear wave equation with the Landau prescription.Our linear calculation gives a definite torque on the planet despite the singularity at corotation.In the next section, we describe the formulation of our linear wave calculation with the Landau prescrip-tion for 3D locally isothermal disks.We also explain the origin of the strong divergence at corotation in 3D locally isothermal disks in the subsections 2.3 and 2.4.In Section 3, we present the wave solutions and the torques due to the waves.We will show that the obtained torques do not depend on the damping parameter of the Landau prescription.Then, we also present the linear torque formulas for 3D locally isothermal disks.In Section 4, we compare the linear torque formula with the results of the previous 3D hydro-dynamical simulations.As well as the linear corotation torque, our linear calculation can evaluate the 3D horseshoe torque, which is also compared with the previous 3D simulations.In the last section, we summarize our findings.

Model and Assumptions
We adopt the following assumptions for density waves in three-dimensional disks, as done by Paper I 1.The amplitudes of the density waves excited by a planet are small enough for the linear calculation to be valid.We also assume that the mass of the planet is so small that the disk gap does not form.The orbit of the planet is circular and lies in the midplane of the disk (i.e., e = i = 0).
2. We consider a locally isothermal disk that is vertically isothermal, but radially non-isothermal.We also assume an isothermal equation of state.The disk temperature and the surface density of the unperturbed disk are both power-law functions of the radial coordinate r.The self-gravity of the disk is negligibly small.
3. The boundary condition is the non-reflecting one for density waves propagating away from the planet.In this paper, the waves are damped near the corotation point with the Landau prescription.Thus, the density waves propagating towards corotation are not reflected at corotation (see Section 3).The Landau prescription is necessary to suppress the strong divergence at corotation.

Basic Equations and the Unperturbed Disk
The basic equations are the Euler equation and the equation of continuity given by where ρ, p, v are the density, the pressure, and the velocity of the disk gas, respectively, and M * is the mass of the host star.The gravitational potential of the planet with the mass M p is given by ϕ p = −GM p /|r − r p |, where r p is the position vector of the planet.To describe the disk, we use a cylindrical coordinate system, (r, θ, z).The z-axis is chosen to coincide with the rotation axis of the disk.The unperturbed surface density σ and the temperature T of the disk depend on r as where r p = |r p | is equal to the radius of the circular orbit of the planet.We adopt the isothermal equation of state p = c(r)2 ρ.Since the isothermal sound speed c depends on r as c = c p (r/r p ) −β/2 , the disk is baroclinic.The unperturbed disk is steady and axisymmetric.The vertical density distribution of the unperturbed disk is where the disk scale height h is given by h = c/Ω K and Ω K ≡ GM * /r 3 is the Keplerian angular velocity on the midplane.The disk scale height depends on r as h = h p (r/r p ) µ and the power-law index is given by µ = 3/2 − β/2.The small disk aspect ratio, h/r, is assumed in Equation ( 5).On the disk midplane (z = 0), the pressure also has the power-law radial distribution of p ∝ r −δ , where the power-law index is δ = 3/2+α+β/2.The angular velocity Ω of the unperturbed disk rotation is where the terms of O((h/r) 4 ) and the higher are neglected.Note that Ω depends also on z in radially nonisothermal disks.

Fourier-Hermite Components of the Perturbation Equations
We describe the equations for the velocity and pressure perturbations u and p 1 .We also introduce a perturbation η defined by Although η defined by Equation ( 7) is the enthalpy perturbation in completely isothermal disks, it is not for the radially non-isothermal case.We expand the perturbations with Fourier series and Hermite polynomials, as in Paper I.
where Z is the non-dimensional vertical coordinate defined by Z = z/h(r) and ℜ represents the real part.The coefficients of the summation (e.g., η m,n (r)) are called the Fourier-Hermite components.In this paper, we consider slowly growing perturbations, by adding the factor e ϵt in Equation ( 8) with a positive constant ϵ, for the Landau prescription.Since ϵ is set to be much smaller than Ω p , the Landau prescription is effective only near corotation.
The linearized Euler equation and equation of continuity for the "Fourier components" are given by where, η ′ m = η m + ϕ p,m .The terms of ϵ represent the Landau prescription.These terms suppress the divergence of waves at corotation.
We further expand these equations with Hermite polynomials (Takeuchi & Miyama 1998;Paper I).For the expansion, each equation is multiplied by H n (Z) e −Z 2 /2 /( √ 2πn!) and integrated with respect to z.For example, the term, im(Ω where and we used the normal orthogonal relation between polynomials H n given by In Equation ( 13), noting that |r − r p |/h p , mh p /r p ∼ O(1), we neglected the terms of O((h p /r p ) 3 ) and the higher.For radially non-isothermal disks, ∆Ω has a nonzero value and Ω n depends on n.Transforming Equations ( 9) and ( 12) in a similar way, we obtain Equations ( 10) and ( 11) are transformed as In Equations ( 17)-( 20), the terms proportional to β (or ∆Ω) come from the radial temperature gradient.Due to the non-zero ∆Ω, each equation has coupling terms with different n modes.By eliminating the velocity perturbations u i,m,n in these equations, we can obtain a second-order differential equation for η ′ m,n (i.e., a wave equation).Due to these eliminations, strong divergences appear at corotation of radially non-isothermal disks.
When β = 0, ∆Ω vanishes and Ω n is independent of n.From Equations ( 19) and (20), we see that the perturbations u θ,m,n and u z,m,n are proportional to 1/(Ω n − Ω p ) and diverge at corotation, where Ω n = Ω p .Then, in the vicinity of corotation, the wave equation has a form when ϵ = 0. Thus, the wave equation has the secondorder pole at corotation where Ω n = Ω p , but it is a regular singular point.
For radially non-isothermal disks with a non-zero ∆Ω, the corotation points of each n mode are located at slightly different radial positions.Furthermore, Equations (19) and ( 20) have the coupling terms with adjacent n ′ modes, and the adjacent n ′ mode waves are also coupled with other modes.Because of this coupling chain, each n mode wave is influenced by the corotation resonances of all the other modes.This coupling effect makes the divergence at corotation stronger than that in a completely isothermal disk.In the next subsection, we will find that the wave equation has the third-order pole at corotation under the modified local approximation, as also shown in Paper I. Note that the terms proportional to im(Ω n − Ω p ) in Equations ( 17) and ( 18) do not cause any divergence.Thus, we will neglect the terms with ϵ in Equations ( 17) and (18) below.

Modified local approximation
Furthermore, the perturbations and their equations are expanded with the small disk aspect ratio, h p /r p .We define the local radial coordinate x as2 where the corotation radius at the midplane, r c,mid , is defined by Ω(r = r c,mid , z = 0) = Ω p .By using x, the augular velociy Ω n defined by Equation ( 15) is given by The unperturbed surface density, sound velocity, and scale height are expanded as σ where λ = (M p /M * )(r p Ω p /c p ) 2 and the subscript i indicates r, θ, or z.The expressions of the potentials Φ m,n,j (j = 0, 1) have been given by the Equation ( 27) of Paper I. We also define W ′ m,n,j ≡ W m,n,j +Φ m,n,j , which is related to η ′ m,n .The zeroth order perturbation W m,n,0 corresponds to the shearing sheet model and the first order perturbation W m,n,1 represents the deviation from the shearing sheet model, which gives the net torque on the planet.
We write down the equations for zeroth and first-order perturbations.From the zeroth and first-order terms with respect to h p /r p in Equations ( 17) to ( 20), we obtain where k θ = mh p /r p , ε = ϵ/Ω p and j = 0 or 1 in each equation.The equation with j = 0 is that for the shearing sheet model and j = 1 corresponds to the equation for the first-order perturbation.In Equations ( 25) and ( 26), terms with ε are omitted.Moreover, the source terms Q 1 's are given by3 .
In the above definitions, Q (X) 1 consist of terms proportional to the temperature gradient (β), surface density gradient (α), or scale height variation (µ), and other terms (that are called "curvature terms").The terms proportional to the pressure gradient (δ) appear in the potential Φ m,n,1 (see Paper I).Since Equations ( 25) to (28) are linear, the first-order solution is written by the superposition of terms related to each effect.
The main target of this paper is the term related to the temperature gradient, W ′(β) m,n,1 since other terms have obtained in Paper I.
We check the divergence at corotation again for Equations ( 25)-( 28) with ε = 0.In the zeroth order equations, U θ,m,n,0 and U z,m,n,0 have the first-order pole at corotation and the wave equations with j = 0 have the following form in the vicinity of corotation.
It has the second-order pole but it is a regular singular point.
For the first-order equations, we have the same wave equation but j = 1.In this case, U z,m,n+1,1 has the source term , which includes the terms U z,m,n±1,0 and U z,m,n−3,0 for a non-zero β (see Equation [32]).Since these terms are proportional to 1/x, U z,m,n+1,1 is proportional to 1/x 2 and the wave equation has the third pole at corotation, which causes a non-removal singular point in the wave solution.In this paper, we solve the wave equation, by suppressing the strong divergence at corotation with the Landau prescription.
Note that the term with the third pole at x = 0 is caused by the modified local approximation.In Equations ( 17) to ( 20), the corotation points of each n mode split into slightly different radial positions due to the ndependence of Ω n .In Equations ( 25)-( 28), on the other hand, the corotation points of all n modes degenerate to x = 0.This degeneracy of corotation leads to the third pole.In the next section, we will examine how this strong divergence at corotation affects the density waves and the planetary torques.

Boundary Conditions
The radiation boundary condition is applied to the solutions at large |x|, as done in Paper I. The radiation condition is written by In the modified local approximation, s(r) is expanded as and S j (x)'s are given by Equation ( 35) gives the boundary conditions for W ′ m,n,j using S 0 and S 1 .To seek the solution related to the temperature gradient W ′(β) m,n,1 , we leave only terms proportional to β in S 1 and put β = 1.The other terms in S 1 are omitted.
As in Paper I, the general solutions of W ′ m,n,j are written as W ′ m,n,j = A 1,j y 1 + A 2,j y 2 + W ′ m,n,j,p , where y 1 and y 2 are two independent solutions to the homogeneous differential equations ( 25) and ( 26), and W ′ m,n,j,p is a particular solution to the inhomogeneous equations.The complex coefficients A i,j are determined by the boundary conditions at x = ±∞.We adopt the following boundary conditions at x = 0.
Noting that the solutions W ′ m,n,j have symmetries as (see Paper I), we find that the coefficients A 1,0 and A 2,1 are real and that A 2,0 and A 1,1 are purely imaginary numbers.Therefore, we only have to integrate the equations ( 25) and ( 26) from x = 0 to +∞.

Total and Corotation Torques
A torque exerted by the planet on a ring region with a radius r and a width of the unit length of the disk (i.e., the torque density) is given by where the relation ρ 1 = ηρ 0 /c 2 is used.By integrating Equation ( 41), we have the torque on the entire disk, and the opposite sign is the torque on the planet.The torque density can be written as the sum of the contributions from each Fourier-Hermite component.
Since η and ϕ p are expanded as Equation ( 8), the contributions from each Fourier-Hermite component are given by where ℑ represents the imaginary part.Moreover, the expansions in the modified local approximation of Equations ( 22) and ( 24) yield where the zeroth-and first-order integrated torques T m,n,j are given by and The definition of Γ 0 is different from that in Paper I. Note that T m,n,0 (∞) = 0 due to the symmetry of the sheering sheet model.Then, integrating Equation ( 44), we find that Γ m,n = Γ 0 (h p /R p )T m,n,1 (∞).The radial integration of the torque density gives the total torque on the entire disk Γ total .From Eq. ( 44), we obtain where we used m = (r p /h p ) dk θ .
In Paper I, the torques are obtained from the angular momentum fluxes of each wave component.However, the angular momentum fluxes vary due to the radial temperature gradient as well as due to the torque exerted by the planet.For radially non-isothermal disks, therefore, the torques should be calculated using the above equations.
Similar to the first-order solution W ′ m,n,1 in Equation (33), the torque contributions are written as the sum of the terms due to each effect.
In this paper, we will obtain the terms Γ T (n).The other terms have been obtained in Paper I.
Furthermore, taking the summation for n and noting δ = 3 2 + α + β/2 and µ = 3 2 − β/2, we obtain the total torque on the whole disk, Γ total , as where the coefficients are given by The corotation torque is also expressed in a similar form.We can also obtain the Lindblad torque by substructing the corotation torque from the total torque.We also derive the formula of the corotation torques due to two-dimensional waves (n=0) for the case of β ̸ = 0. Near corotation in radially non-isothermal disks, the wave equation of η ′ m,0 is given by with where we have used η ′ m,2 ≃ 0 at r = r c (e.g., Appendix of Paper I).Then, evaluating the jump in the angular momentum flux, we obtain the corotation torque by twodimensional waves, Γ C,m,0 , as When β = 0, the above equation reduces to the expression in Paper I (Eq.[55]), in which isothermal disks are considered.The corotation torque is also written as the sum of terms due to each effect as in Equaition (50).
The term due to the temperature gradient is given by In Paper I, the corotation torque also has contributions from 3D waves.However, as shown in the next section and Appendix A, these torque contributions come from the 3D waves excited by the inner and outer nearside Lindblad resonances and should be included in the Lindblad torque.Therefore, in the present paper, the corotation torque is considered to result only from the 2D waves and is given by equation (57).We obtain the Fourier-Hermite components of the enthalpy perturbation, η m,n , by integrating Equations ( 25) and ( 26) from x = 0 to +∞ with boundary conditions (39) and (35), and using Equation (24). Figure 1 shows typical solutions of η m,n .Focusing on the effect of the disk temperature gradient, we set the first-order solution to be W m,n,1 = βW (β) m,n,1 .We can the solution W (β) m,n,1 , by leaving only the source terms proportional to β in the first-order equations ( 25) to (28).Using this firstorder solution, we plot η m,n for n = 0, 2, 4 in Figure 1.To emphasize the effect of the temperature gradient, we set β = 2.The left panel shows the full view and the right panel is the magnified view near corotation.Both the real and imaginary parts are plotted for solutions with the damping parameter ε = 10 −4 .Slight asymmetries of the waves for r = r p are due to the temperature gradient.The imaginary parts of the solutions with ε = 10 −3 , 10 −5 , and 10 −6 are also plotted.The solu-tions with different ε overlap almost completely, except near corotation.

Wave Solutions
The wave solutions with the smaller ε have stronger oscillations near corotation because a small ε has a weak damping effect.These oscillations result from the strongly divergent term in the wave equation for threedimensional waves.Although the wave with n = 0 also has some oscillations near corotation, the amplitude near corotation is much smaller than that of the waves propagating away from the planet (in the right panel).The oscillations of the wave with n = 0 are caused by the interaction with the three-dimensional wave with n = 2.
Note that the strongly divergent term affects the wave only near corotation.The limited effect of the divergent term can be explained by the wave propagation.The excitation and propagation of density waves are described in Appendix A. The three-dimensional waves with n ≥ 2 are excited at four Lindblad resonances for each n mode (see Figure 6).The waves around the planet originate from the near-side Lindblad resonances and propagate towards corotation4 .These waves are expected to be dissipated at corotation due to their extremely short wavelength.Because of the direction of wave propagation, the divergent term at corotation has no effect outside the vicinity of corotation.
In Paper I, we considered the torque associated with the waves propagating away from the planet to be the Lindblad torque and classified the torque associated with the three-dimensional waves propagating towards corotation as the corotation torque.However, the latter torque should be classified as the near-side Lindblad torque, since these waves are excited at the near-side Lindblad resonances.Thus we adopt this new classification and consider the corotation torque to result only from the two-dimensional waves.

Torques
Next, we examine the torques exerted on the disk by the wave excitation.Figure 2 shows the cumulative radial distributions of the torques due to the density waves in Figure 1.The cumulative distribution of the torque due to η m,n is given by We plot the cases where ε = 10 −3 , 10 −4 , 10 −5 , and 10 −6 , The lines for each ε almost overlap and show only small differences between them even near corotation (see the right panel).Although there are some oscillations near corotation, the amplitudes of the torque oscillations are small compared to those of η m,n in Figure 1.Furthermore, the ε dependence of the integrated torques is similar between before and after corotation, indicating that the contributions near the corotation are quite small for the integrated torques.By integrating the torque densities to a sufficiently large r p , we obtain the convergence values of the torques due to each wave.Since the integrated torque converges slowly as it oscillates, we obtained the convergence value by calculating the average value of the oscillations.The convergence values of each integrated torque are shown in the left panel.From the convergence value T m,n,1 (∞), the torque, Γ m,n , is obtained.Since we only consider the effect of the disk temperature gradient here, we obtain Γ  2, there is little dependence on ε.For n = 0, the torque with ε = 10 −3 is slightly different from those with ε = 10 −4 and 10 −5 .Thus, the value of the torque is almost independent of ε as long as ε ≤ 10 −3 .For small k θ , it is difficult to solve the wave equation numerically due to a small ε of the Landau prescription.For the two-dimensional waves with n = 0, fortunately, we can solve the wave equations without the Landau prescription as done in Paper I, because the wave equation with n = 0 does not have the strongly divergent term.Thus, for 0.03 < k θ < 0.05, the torque with n = 0 is obtained by solving the wave equations with the method of Paper I instead of the Landau prescription, which is shown with a green line in Figure 3.The torques obtained by the two methods are almost continuous at k θ = 0.05, confirming that both methods with and without ε are consistent.
The bottom panel of Figure 3 shows the results with ε = 10 −4 for n=0 to 14 (solid lines are positive and dotted lines are negative).In addition, for n = 0, the corotation torque given by Equation ( 57) is also plotted.For k θ < 0.1, the corotation torque and the total torque are almost identical.This means that the Lindblad resonance is much weaker in this region.As n increases, the peak value of each curve decreases but the decrease is not so rapid.As will be shown below, the contribution of larger n is not negligible.
For k θ ≲ 0.05, we could not accurately obtain the wave solutions with n ≥ 2. However, for the threedimensional waves with such small k θ , the torques are negligible compared to their peaks at k θ ∼ 1. There- fore, we neglect the contribution of such small k θ for the three-dimensional waves.
Note-The unit of the torques is Γ0.power-law function with the exponent of -2.18.With the fitting formula, we can extrapolate Γ (β) T (n) for n > 14.The torque contribution of n > 14 estimated by the fitting formula is also shown in Table 1.The contribution of n > 14 is ≃ 20% 5 .Using this extrapolation, the total torque is finally obtained as The summed corotation torque, Γ CT (0), is obtained as -0.8575, which is much larger than the total torque directly due to the temperature gradient, n Γ (β) T (n).This is because the Lindblad and corotation torques almost exactly cancel each other out for the effect of the temperature gradient.A similar cancellation between 5 For other effects, Γ (X) T (n) also decrease with power-laws and similar extrapolations are done to obtain the summation of Γ the Lindblad and corotation torques occur for the effect of the surface density gradient, which also results in the small total torques.
Note that the effect of the temperature gradient also appears indirectly through the effects of the pressure gradient and the scale height variation.Since h = Ω/c s and p ∝ T Σ/h, the power-law indices depend on β as µ = 3/2−β/2 and δ = α+3/2+β/2, respectively.Thus, the torque term proportional to β is obtained as The torque terms due to the other effects were obtained in Paper I. However, we have found some errors in the corotation torques of three-dimensional waves.These errors in the corotation torques also affect the total torques.We list the revised coefficients of the corotation torque in Table 2 as well as the torques due to the temperature gradient obtained in this paper.The revised coefficients of the total torque are given by These revised values are used in Equation ( 61).Note that the torque contributions from the 3D waves in Table 2 are included in the Lindblad torque according to the new definition of the corotation and Lindblad torques although they belonged to the corotation torque in Paper I. Using these results of the linear calculations for each effect, we obtain the formula of the total torque for lo-cally isothermal disks as Γ total /Γ 0 = 1.436 + 0.537α + 0.439β. (63) The last term due to the temperature gradient is newly added in this study.The terms due to the curvature and the surface density gradient are slightly changed from the values of Paper I because of the errors in the 3D corotation torques in Paper I. Note that this is the torque on the disk.The torque on the planet has the opposite sign.We can divide it into the Lindblad torque and the corotation torque.The corotation torque is given by the terms due to the two-dimensional waves in Table 2 and the Lindblad torque is given by Γ total − Γ C .Thus we obtain Their value of dΓ total /dβ is close to ours, 0.439Γ 0 .The difference is only 2%.The other coefficients Γ α=β=0 total and dΓ total /dα have deviations of 5-15% from our values.This difference could be due to the effect of the nonlinear horseshoe torque.We will discuss the horseshoe torque in Section 4.3.Jiménez & Masset (2017) also examined the temperature gradient dependence of the total torque with their 3D hydrodynamical simulations for locally isothermal disks.Their result is dΓ total /dβ = (0.5 ± 0.1)Γ 0 , which agrees with our value within their error bar.Paardekooper et al. (2010) suggested that the coefficient dΓ total /dβ (or dΓ C /dβ) is not affected by the horseshoe torque in locally isothermal disks.The agreement in dΓ total /dβ with hydrodynamical simulations supports Paardekooper et al.'s suggestion.
In the present study, we obtained the temperaturegradient term of the linear corotation torque as β(dΓ C /dβ) = −0.858βΓ in Equation ( 65).In Jiménez & Masset (2017), the corresponding term of the linear corotation torque is set to be −Γ lin T = −1.0βΓ0 , which should be replaced by our value, −0.858βΓ 0 .2010) constructed the torque formulas based on their 2D linear calculations.Figure 5 shows a comparison of our 3D total torque and corotation torque with Paardekooper et al.'s torque formulas for locally isothermal disks.Since their 2D torque formulas depend on the smoothing length included in the gravitational potential, Figure 5 shows each torque coefficient as a function of the smoothing length.Furthermore, we also performed 2D linear calculations and plotted the results as well.Although the torque formulas of Paardekooper et al. (2010) are expressed by simple power-law fits, their formulas are almost in agreement with our 2D linear calculations.The torque coefficients of the 3D linear calculations agree with Paardekooper et al's formulas with a smoothing length of ≃ 0.4h for both the curvature term and the surface density gradient term of the total and corotation torques.On the other hand, the temperature gradient term obtained in the present 3D calculation agrees with the formula with s ≃ 0.6h for the total torque, while a good agreement between the 3D result and their formula is seen in the wide range of s=0.4-0.7h for the corotation torque.These comparisons indicate that it is difficult to reproduce all the dependencies in the 3D torques with the 2D torque formulas with a single smoothing length.

Horseshoe Torque in Three-Dimensional Disks
In hydrodynamical simulations, the torque contribution around corotation is often closer to the horseshoe torque described below than to the linear corotation torque because of nonlinear effects, though it depends on the saturation (e.g., Masset & Casoli 2010;Paardekooper et al. 2011).For 2D isothermal disks, the unsaturated horseshoe torque on the planet, Γ C,hs , is given by (e.g., Ward 1991;Masset 2001) where x hs is the half-width of the ring-like horseshoe region.Masset et al. (2006) showed that the half-width x hs is given by for 2D disks.In this equation, the suffix S indicates that the perturbation, η ′ (= η + ϕ p ), is evaluated at the stagnation point of the horseshoe flow.Since η ′ S ∼ O(λr p /h p ), the horseshoe torque is obtained as A(3/2 − α)Γ 0 , where A is a coefficient of the order of unity.The half-width of the horseshoe region, x hs , can also be obtained from the streamlines of the horseshoe flow in hydrodynamical simulations.In 2D disks, x hs and the horseshoe torque depend on the smoothing length (e.g., Masset et al. 2006;Masset & Casoli 2010;Paardekooper & Papaloizou 2009b;Paardekooper et al. 2010).Note that in completely isothermal disks, the horseshoe torque is proportional to the vortensity gradient (i.e., the factor, 3/2 − α) as well as the linear corotation torque of Equation ( 65).Masset & Benítez-Llambay (2016) performed 3D isothermal simulations to obtain the horseshoe torque.Their 3D simulations show that the horseshoe flow is almost two-dimensional and nearly independent of z.Thus, Equation (67) for the horseshoe torque in the two-dimensional case can apply to the three-dimensional case.From their 3D simulations, they obtained the width of the horseshoe region as where q = M p /M * , and the horseshoe torque is This value of the horseshoe torque for 3D isothermal disks is consistent with the torque formula by Paardekooper et al. (2010) if the smoothing length of ≃ 0.5h is used.Moreover, it is about 1.4 times as large as the vortensity term of our linear corotation torque of Equation ( 65).Note that the simple horseshoe torque formula of Equation ( 70) is valid only for low-mass planets with q ≤ 0.1(h p /r p ) 3 (Jiménez & Masset 2017).
Here we also obtain the horseshoe torque by estimating η ′ S with the 3D linear calculation.Since the horseshoe flow is shown to be two-dimensional, η ′ is estimated using the two-dimensional waves with n = 0, and the three-dimensional waves with n > 0 are negligible.Thus, η ′ is obtained by summing only the terms for n = 0 over m in Equation ( 8) for the horseshoe region.Since the stagnation point is approximately at the corotation radius, we take the minimum value of η ′ at x = 0 as η ′ S .In the summation, the term for m = 0 is omitted.In addition, the term for m = 1 originally contains an indirect term, which is also ignored in the modified local approximation.
For comparison with Masset & Benítez-Llambay (2016), we set h p /r p = 0.05.We also include the firstorder terms of the modified local approximation.In the first-order terms, we include the curvature terms and set α = 3/2, β = 0, and δ = 3/2 + α = 3 as in the nominal case of Masset & Benítez-Llambay (2016).As a result, we obtain This is close to the result of Masset & Benítez-Llambay (2016), with a deviation of several % in the horseshoe torque.Thus, the linear calculation gives a good approximation of the horseshoe torque.
The results of Equation ( 71) are derived using the parameter set of Masset & Benítez-Llambay (2016).Note that the numerical coefficients in Equation ( 71) depend on the parameters.To obtain η ′ S above, we included the first-order terms of the modified local approximation, which depend on the values of α and δ.
If we ignore the first-order terms and use the shearing sheet model, the horseshoe torque is 3% smaller than Equation (71).The normalized horseshoe torque Γ C,hs /Γ 0 also depends on the disk aspect ratio h p /r p ; for h p /r p = 0.03, Γ C,hs /Γ 0 = 0.952(3/2−α), which is about 10% larger than Equation (71).The similar dependence of the horseshoe torque (or the width) on h p /r p is seen in the simulation results of Jiménez & Masset (2017) (in their Fig. 6).
In locally isothermal disks, the horseshoe torque also has the term related to the radial temperature gradient.In Jiménez & Masset (2017), such a term of the unsaturated horseshoe torque on the planet is given by for a low-mass planet.This is also close to our temperature-gradient term of the linear corotation torque of Equation ( 65).This agreement between the horseshoe and linear corotation torques in locally isothermal disks has already been suggested by Paardekooper et al. (2010).

SUMAMRY
In the present study, we have performed linear calculations to determine the torques of Type I planetary migration for 3D locally isothermal disks with radial temperature gradients, whereas Paper I gave the torques for completely isothermal disks.In a disk with a radial temperature gradient, the wave equation has a term that strongly diverges at the corotation, which causes a non-removal singularity in the wave solution.As a result, obtaining the wave solutions was challenging.We obtained the wave solutions, suppressing the divergence with the Landau prescription.Our results are summarized as follows.
1.The waves excited by a planet have a large amplitude near the corotation due to the divergent term, and the degree of amplification depends on the parameter of the Landau prescription (Fig. 1).Our wave solutions show that this wave amplification at corotation does not affect the region away from corotation.This limited influence of the divergent term is explained by the fact that the waves propagate towards corotation and dissipate at corotation.We also found that this wave amplification does not affect the torque on the planet due to the very narrow range of the wave amplification (Fig. 2).Thus, we can obtain a definite torque even though the wave equation has a divergent term of the third pole (Fig. 3).
2. The total torque was obtained by summing the torque contributions of each wave component.The obtained torque directly due to the temperature gradient is smaller than the torque due to other effects obtained in Paper I. Thus, the torque due to the temperature gradient mainly results from the pressure gradient and the scale height variation indirectly.The total torque considering all effects, including the temperature gradient, is given by Equation (63).
3. Equation ( 57) shows the revised formula of the corotation torque for locally isothermal disks with temperature gradients.Due to the new definition of corotation torque, three-dimensional waves do not contribute to the corotation torque.The revised formula gives the corotation torque as Equation (65), and the Lindblad torque is given by Equation ( 64).2010) based on 2D linear calculations requires a large smoothing length of 0.6h to make the torque due to the temperature gradient consistent with the 3D result.Therefore, it is difficult to reproduce all the dependence of the 3D torque on each effect using a 2D calculation with a single smoothing length.
5. Equation (71) shows the horseshoe torque obtained through our 3D linear calculation.This result is also close to that of 3D hydrodynamical simulations by Masset & Benítez-Llambay (2016).
The temperature gradient term of the horseshoe torque is very close to that of our linear corotation torque for locally isothermal disks.

APPENDIX
A. EXCITATION AND PROPAGATION OF DENSITY WAVES The excitation and propagation of density waves in two-and three-dimensional disks can be explained in terms of dispersion relations, and of course, such an explanation has already been given by previous studies (e.g., Lin & Shu 1968;Goldreich & Tremaine 1979;Lubow & Pringle 1993;Takeuchi & Miyama 1998).Here we briefly review their results on non-self-gravitating gaseous disks.In the WKB approximation, the radial dependence of the Fourier-Hermite components is written in the form exp(i r k r dr) and substituted into Equations ( 17)-( 20) to obtain the following dispersion relation (Takeuchi & Miyama 1998) 6 . where In the above derivation, we have neglected the terms proportional to 1/r in the WKB approximation and also omitted the terms proportional to the small ∆Ω.As a result, Equation (A1) agrees with the dispersion relation in the shearing sheet model.The approximated dispersion relation is valid for studying the qualitative properties of wave excitation and propagation.
The r-dependence of k 2 r is illustrated in Figure 6.In the case of two-dimensional density waves (n = 0), the wavenumber k r vanishes at the (original) inner and outer Lindblad resonances, r in L,0 , r out L,0 , where D 1 = 0 or Ω K = m m±1 Ω p .Density waves are excited at these Lindblad resonances and propagate away from the planet.The propagation regions have a positive k 2 r .As the waves propagate far away from the resonances, both the wavenumbers |k r | and the amplitudes increase.As a result, the density wave decays with shock dissipation (Goodman & Rafikov 2001).In the region between the two Lindblad resonances, k 2 r is negative and thus the waves are evanescent.The corotation point, r c , and the nearby radial position of the planet, r p , are almost in the center of this region.
For three-dimensional waves with n ≥ 2, k 2 r is positive even near corotation.The wavenumber also vanishes at the new inner and outer Lindblad resonances, r in The former waves correspond to the p-mode with high frequencies as well as the two-dimensional waves due to the original Lindblad resonance, and the latter waves correspond to the g-mode (or the r-mode) with low frequencies (Lubow & Pringle 1993;Lubow & Ogilvie 1998).Since the wavenumber is divergent at corotation, the latter waves are dissipated in the neighborhood of corotation.
The dispersion relation (A1) also gives the radial group velocity of the waves, v g,r .Since the angular frequency of the waves is given by mΩ p in Euqation (8), we obtain The sign of v g,r determines the direction of radial propagation.We note that the directions of propagation in Figure 6 are consistent with the signs of v g,r in each wave region when the wavenumber k r is positive (i.e., the waves are trailing).In fact, trailing density waves are obtained in our linear calculations (see Figure 1 as well as in the previous three-dimensional linear calculations (Takeuchi & Miyama 1998;Paper I).This is obtained by applying the radiation boundary condition.

Figure 1 .
Figure 1.Waves in a radially non-isothermal disk.Fourier-Hermite components of the enthalpy perturbation are shown for n = 0 (top), 2 (middle), and 4 (bottom).The left panel shows the full view and the right panel shows the magnified view near corotation.We set α = µ = δ = 0 and omit curvature source terms in the wave equations to focus on the effect of the radial temperature gradient.Other parameters are β = 2, m = 10, and hp/rp = 0.05.The black solid lines and dotted lines show the imaginary and real parts of the wave with ε = 10 −4 .In the left middle and left bottom panels, the real parts of η ′ m,n (= ηm,n + ϕm,n) with ε = 10 −4 are also plotted with the gray dotted lines.The imaginary parts of the solutions with ε = 10 −3 , 10 −5 , and 10 −6 are plotted with the blue, red, green lines, respectively.
. For the waves with different m and n, similar integrations of the torque density give the torques, Γ (β) m,n , and the results are shown in Figure 3.In the top panel of Figure 3, we compare the torques with ε = 10 −3 , 10 −4 , and 10 −5 , for n = 0, 2, and 4. As expected from Figure

T
(n)  is well fitted by a

Figure 5 .
Figure 5. Coefficients of the total torque and the corotation torque in two-and three-dimensional linear calculations.The torque coefficients are plotted as functions of the smoothing length s.The light blue and light red lines show the coefficients in the torque formulas of Paardekooper et al. (2010) for the locally isothermal case.The horizontal broken lines represent the 3D linear results of the present paper.Our 2D linear results are plotted with the blue and red lines.

4.
The torque due to disk temperature gradient obtained by our 3D linear calculation agrees well with the result of the previous 3D hydrodynamical simulations for locally isothermal disks (D'Angelo& Lubow 2010;Jiménez & Masset 2017).The torque formula ofPaardekooper et al. ( L,n , r out L,n , where D n = 0 (or Ω K (r) = m m± √ n Ω p ).The new Lind-6 Lubow & Pringle (1993) obtained a similar dispersion relation for 3D axisymetric waves though they also include the buoyancy using the polytropic equation of state.blad resonances are farther away from the planet (or r c ) than the original Lindblad resonances.The threedimensional waves with each n mode are excited at these four Lindblad resonances.The waves excited at the far-side Lindblad resonances propagate in the far region where D n < 0, while the waves excited at the near-side Lindblad resonances propagate towards corotation in the near region where D 1 > 0.

Figure 6 .
Figure 6.The r-dependence of the square of the wavenumber k 2 r for two-dimensional waves with n = 0 (upper panel) and three-dimensional waves with n ≥ 2 (lower panel).Density waves are excited at the Lindblad resonances, where kr = 0, and propagate into the regions where k 2 r > 0.

Table 1 .
Torques due to the radial temperature gradient, Γ

Table 2 .
Coefficients of the Corotation Torques