{'fraction_non_alphanumeric': 0.04031108898312359, 'fraction_numerical': 0.012343983381118341, 'mean_word_length': 4.343362994220714, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the corotation torque in a radiatively inefficient disk

We consider the angular momentum exchange at the corotation resonance between a two-dimensional gaseous disk and a uniformly rotating external potential, assuming that the disk flow is adiabatic. We first consider the linear case for an isolated resonance, for which we give an expression of the corotation torque that involves the pressure perturbation, and which reduces to the usual dependence on the vortensity gradient in the limit of a cold disk. Although this expression requires the solution of the hydrodynamic equations, it provides some insight into the dynamics of the corotation region. In the general case, we find an additional dependence on the entropy gradient at corotation. This dependence is associated to the advection of entropy perturbations. These are not associated to pressure perturbations. They remain confined to the corotation region, where they yield a singular contribution to the corotation torque. In a second part, we check our torque expression by means of customized two-dimensional hydrodynamical simulations. In a third part, we contemplate the case of a planet embedded in a Keplerian disk, assumed to be adiabatic. We find an excess of corotation torque that scales with the entropy gradient, and we check that the contribution of the entropy perturbation to the torque is in agreement with the expression obtained from the linear analysis. We finally discuss some implications of the corotation torque expression for the migration of low mass planets in the regions of protoplanetary disks where the flow is radiatively inefficient on the timescale of the horseshoe U-turns.


Introduction
It is known since the early eighties that low mass planetary objects (that is, up to a few Earth masses) embedded in protoplanetary gaseous disks should undergo a fast decay towards their central object, on timescales much shorter than the lifetime of the disk. This process, known as type I migration, has constituted for a long time a critical stage for the theory of giant planet formation. While it may account for the discovery of close-in extrasolar planets, with orbital periods of a few days, it renders problematic the build up of giant planet cores at distances of their central stars of several astronomical units. Most published studies of the tidal interaction of low mass objects with their parent disk have used either a barotropic assumption (such as a polytropic equation of state), or a locally isothermal equation of state. All these studies, whether analytical or numerical, confirmed the vigorous tidal interaction of the planet with the disk, leading to its inward migration on short timescales.
There has been some exceptions to such assumptions: Morohoshi & Tanaka (2003) considered the case of a planet interacting with an optically thin disk, in the shearing sheet approximation, and found that radiative effects can significantly alter the one-sided torque between the planet and the disk. More recently,   (hereafter PM06) have performed global, high resolution 3D calculations with nested grids that include radiative transfer. For the setup that they considered, they found that the total torque exerted by the disk on the planet increases with the disk opacity. For sufficiently large values of the opacity (and in the limit case of an adiabatic flow, corresponding to an infinite opacity), they find that the total torque on the planet is positive. This result is of great importance, as it potentially solves the lingering problem of type I migration. PM06 identified the existence of a hot, underdense part of the co-orbital region lagging the planet, which accounted for the torque excess that they measured. The present work corresponds to an attempt to further investigate this topic, so as to identify the physical mechanism responsible for these effects. For this purpose, we consider a more restricted situation, namely two-dimensional adiabatic flows.
This paper is organized as follows. In section 2 we set up the problem and define the notation. We then present an analysis of the corotation torque in an adiabatic disk in the linear regime, at an isolated resonance, at section 3. Our original motivation for the study of the linear regime was that PM06 found that the total torque reverses in a radiatively inefficient disk both for a 5 M ⊕ and a 0.5 M ⊕ planet, which pointed out that the effect is likely a linear one. In section 4, we check by means of customized two-dimensional hydrodynamical simulations involving an isolated resonance the torque expression found in section 3. In section 5, we turn to the case of a planet embedded in an adiabatic disk, for which we check that there is an excess of corotation torque that scales with the entropy gradient. We also check in this section that the torque excess corresponds to the sum of the linear contributions of all co-orbital corotation resonances, for a sufficiently small planet mass. We discuss the implications of the modified corotation torque expression for the issue of planet-disk tidal interactions, and we suggest further research on this topic in section 6. We sum up our results in section 7.

Setup and notation
We consider an inviscid, radiatively inefficient (that is to say, for our purposes, adiabatic) two-dimensional disk. In order to avoid corotation torque issues, we shall consider either a potential slowly turned on (sections 3 and 4), or early stages after the introduction of a planet (section 5). The unperturbed state of the disk corresponds to a rotational equilibrium between the gravitational force of the central object, the pressure gradient and the centrifugal force. The unperturbed state is axisymmetric. The disk rotates with the angular speed Ω(r), where r is the distance to the central object. We denote by p the pressure, Σ denotes the surface density, u and v respectively the radial and azimuthal velocities, ϕ the azimuthal angle. We denote by a "0" subscript the unperturbed quantities, and with a "1" subscript the perturbed ones. For instance, p(r, ϕ) = p 0 (r) + p 1 (r, ϕ). We shall essentially consider disks in which the unperturbed pressure and density are power laws of the radius, respectively with index λ and σ: We shall make use of the two Oort's constants: which scales with the local shear in the flow, and: which is half the vertical component of the flow vorticity, and which is also (2r) −1 times the radial derivative of the specific angular momentum. We will also use the epicyclic frequency κ = (4ΩB) 1/2 .

Basic equations
We study the linear response of the disk to a perturbing non-axisymmetric potential Φ(r, ϕ) = Φ m (r) cos[m(ϕ − Ω p t)]. The perturbing potential rotates at constant angular velocity Ω p . In the inertial frame, the linearized Euler equations of the disk are: and ∂v 1 ∂t The linearized continuity equation is: We refer to the quantity S = pΣ −γ as the gas entropy, where γ is the adiabatic index. The energy equation is equivalent in our case to the conservation of the gas entropy. The linearized conservation of the entropy along a fluid element path reads where S 1 = S 0 (p 1 /p 0 −γΣ 1 /Σ 0 ). We furthermore assume that the gas is described by an ideal equation of state so that p 0 and Σ 0 are connected by p 0 = Σ 0 c 2 s /γ, c s being the adiabatic sound speed.
We assume a perturbation of the form x 1,m (r) exp(im{ϕ−Ω p t}) where x 1 stands for any perturbed quantity of the flow 1 . We note ∆ω = m(Ω p − Ω) and we use the prime notation to denote ∂/∂r. Eq. (8) can be recast as: Combining Eqs. (5), (6) and (9) we are led to: and where S and P are given by and where Ψ is defined as and where F is defined by with D = κ 2 − ∆ω 2 − SPc 2 s /r 2 . Substituting Eqs. (9), (10) and (11) into Eq. (7) leads to where: and Eq. (16) reduces to the equation (15) of Li et al. (2000) if one considers the propagation of free waves (Φ = 0), while it reduces to the equation (13) of Goldreich & Tremaine (1979) in the case of a homentropic (S = 0) flow.

Corotation torque
We now estimate the rate of angular momentum exchanged between the perturber and the radiatively inefficient disk described in section 3.1. This rate therefore corresponds to the disk torque, which we denote by Γ, and which we define as the torque exerted by the disk on the perturber (unless otherwise stated). It reads: We limit ourselves to the torque exerted by the disk material lying in the vicinity of corotation, hence to the corotation torque, which we denote by Γ c . In a linear analysis, this torque can be expressed as a series of contributions at each azimuthal wavenumber: Γ c = m Γ c,m . Each individual torque can be expressed, assuming that Φ is real, as: where ℑ denotes the imaginary part, r c is the corotation radius, and x = (r − r c )/r c . We substitute Σ 0 u 1 in Eq. (9) by the expression given by Eq. (10), and we keep only the terms which are large in the vicinity of corotation. As in Goldreich & Tremaine (1979), we assume that the disk responds to a slowly increasing perturbation and take ∆ω to have a small, positive imaginary part α: where ǫ = −α/[mr c Ω ′ (r c )] > 0. In the vicinity of corotation, we can finally write: We are primarily interested in the imaginary part of Σ 1 . In the limit ǫ → 0, we can write the terms that yield a non-vanishing contribution to the torque as: where δ(x) is Dirac's delta function. The first two terms of the R.H.S. of Eq. (25) yield respectively the following contributions to the corotation torque: -7 - The third term of Eq. (25) yields a contribution that can be shown to be negligible, in the planetary context, compared to Γ c,m,2 . This is shown in appendix A.
The first term, Γ c,m,1 , is the contribution of the function Ψ, such as in the barotropic case. The second term, Γ c,m,2 , corresponds to a singularity at corotation, associated to a non-vanishing entropy gradient. It corresponds to the torque arising from the advection of entropy in the corotation region, which results in a surface density perturbation if the entropy is not uniform. The perturbation is singular for the surface density and the entropy, but not for the pressure (see section 3.2.2). It remains confined to corotation, where it yields a singular contribution to the torque. Some further insight into the dynamics of this perturbation will be given in section 4.3.
We provide in the next section an expression for the corotation torque in the limit of a cold disk, then we turn to the general case.

Limit of a cold disk
We contemplate here the case for which |Ψ| ≪ |Φ|, which we shall refer to as a cold case. This condition depends on the strength of the perturbing potential, its radial scale, and on the disk temperature. In particular, in the planetary context, some corotation resonances may correspond to a cold situation, while others have |Ψ| ∼ |Φ|. Nevertheless, a given resonance eventually satisfies the cold case condition as the disk temperature tends to zero.
The evaluation of Eq. (26) requires an explicit expression for Ψ, obtained by solving the differential equation (16) in the vicinity of corotation. This has been done by Goldreich & Tremaine (1979) for a cold barotropic disk. An explicit solution can also be obtained for a cold adiabatic disk within the same level of approximation. Note however that some additional difficulties arise, in particular the existence of a double pole (term proportional to ∆ω −2 ) in the coefficients C and D, defined respectively by Eqs. (18) and (19).
We discard the double pole for the following reasons: • Unlike the simple pole, it scales with c 2 s , which indicates that when the disk aspect ratio tends to zero, it becomes negligible; differently stated, there should be a critical disk thickness under which it is safe to neglect this term.
• This term is the only one that depends both on the entropy and on the pressure gradients. As we shall see in section 5.2.2, our results of numerical simulations for a planet embedded in a disk with aspect ratio h = 0.05 show that the torque excess with respect to an isothermal situation essentially depends on S, the gradient of entropy, which indicates that already for h = 0.05 the double pole term is negligible.
• The double pole is regularized with a very small amount of dissipation. Even the molecular viscosity suffices to render it negligible in the disks that we consider (S.-J. Paardekooper, private communication).
The second contribution to the corotation torque, given by Eq. (27), is specific to the adiabatic case and involves the singularity arising from the entropy advection. Using Eq. (27) and noting that |ℜ(Ψ)| ≪ |Φ|, this contribution to the corotation torque reads From Eqs. (30) and (31), we infer the corotation torque for a cold, adiabatic disk, which reads: This expression does not depend on S. We note from Eqs. (15) and (20) that V can be approximated as since the disk aspect ratio at corotation h(r c ) = c s (r c )/[r c Ω(r c )] satisfies h(r c ) ≪ 1. Eq. (32) therefore corresponds to the corotation torque expression 2 of Goldreich & Tremaine (1979). This argues that the corotation torque for a cold case does not depend on whether the disk can radiate energy efficiently (assuming a locally isothermal equation of state) or not (assuming an adiabatic energy equation). This can be expected on general grounds: in the cold disk limit, the internal energy of the fluid is negligible with respect to its mechanical energy, hence the power (and the torque) of the tidal force correspond to the case of noninteracting test particles, for which the expression of Goldreich & Tremaine (1979) prevails.

General case
We consider in this section the general case where we cannot neglect Ψ with respect to Φ in Eqs. (26) and (27), as we have done in the previous section. Instead of resorting to a solution of Eq. (16), we shall use a method similar to the method used by Tanaka et al. (2002), based on the jump of angular momentum flux at corotation. In the case of Tanaka et al. (2002), this eventually yields a torque expression similar to the expression of Goldreich & Tremaine (1979), except that Φ has to be substituted by Φ + η (where η is the enthalpy perturbation). The drawback of this method is that it provides a torque expression that depends on the (unknown) solution of the differential equation. Nevertheless, it gives some insight into the dynamics of the corotation region, and allows to draw the general trends of the corotation torque in an adiabatic disk. In our case, the torque expression features Ψ = p 1 /Σ 0 . We note that in the isothermal case, Zhang & Lai (2006) have provided an explicit solution for the perturbed enthalpy at corotation, that leads to a corotation torque expression that only depends on the forcing potential.
We note that the jump of angular momentum flux at corotation misses the singular contribution of the entropy perturbation at corotation and as such leads only to an evaluation of Γ c,m,1 . The contribution Γ c,m,2 of the entropy perturbation to the corotation torque needs to be calculated similarly as in Eq. (31). The angular momentum flux F A is given by: where ℜ stands for the real part and the star superscript denotes the complex conjugate. Eq. (34) can be written as F A = m F A,m with: Combining Eqs. (10), (11) and (35), we obtain In the homentropic (S = 0) case, Eq. (36) reduces to the expression used by Tanaka et al. (2002). The contribution Γ c,m,1 to the corotation torque is then given by: where r + c > r c and r − c < r c are the radii of locations respectively beyond and before corotation, and where we evaluate the flux of advected angular momentum. Tanaka et al. (2002) showed that Φ + η is continuous at corotation. Here, since Eq. (16) cannot be recast as an ordinary differential equation involving only Φ+Ψ, we have to consider more stringent albeit reasonable assumptions, namely that both Φ and Ψ are continuous at corotation. The fact that Φ is continuous at corotation can be realized with an arbitrarily small softening length of the potential, in the case of an embedded point-like mass (for which the potential components would diverge logarithmically at corotation, in the absence of any softening). Assuming that Φ is continuous at corotation, Eq. (16) imposes that Ψ is also continuous at corotation (we would otherwise have a null linear combination of δ(x) and δ ′ (x) functions with non-vanishing coefficients, which is impossible).
Our continuity assumption implies that the terms proportional to S in the R.H.S. of Eq. (36) does not contribute to the torque. The jump in the advected flux therefore comes from the jump in d(Φ + Ψ)/dr. We integrate Eq. (16) over an infinitesimal interval containing r = r c . All finite terms in this equation yield a vanishing contribution, hence we are left only with the jump of d(Φ + Ψ)/dr stemming from the second derivative term of Eq. (16) and the poles of the terms CΨ and DΦ. This reads: where Using Eqs. (36), (37), (38) and 2B = κ 2 /2Ω, we find that Eq. (39) reduces to Eq. (30) in the cold disk limit.
We now come to the contribution Γ c,m,2 of the entropy perturbation to the corotation torque. Eq. (27) yields: Eq. (40) reduces to Eq. (31) in the cold disk limit.
The general expression for the corotation torque is obtained by accounting for the contribution given by Eq. (39), and that of the entropy perturbation, given by Eq. (40): Eq. (41) reduces to the expression of Tanaka et al. (2002) when S = 0, while it reduces to that of Goldreich & Tremaine (1979) for a cold disk.
A case of interest is that of a disk perturbed by a peaked potential (that of an embedded protoplanet for instance), for which |Φ+ℜ(Ψ)| ≪ |Φ|, and |Φ+ℜ(Ψ)| ≪ |ℜ(Ψ)| at corotation. For such case, |Φℜ(Φ + Ψ)| rc ≫ |Φ + Ψ| 2 rc , hence the corotation torque may be approximated as Γ c,m ≈ −2Γ 0 [S Φ ℜ(Φ + Ψ)] rc . The corotation torque may therefore be much larger in the non-homentropic case (S = 0) than in the homentropic case (S = 0). Furthermore, its sign is given by that of S rather than that of V. This enhancement of the corotation torque in an adiabatic flow may have a dramatic impact on the type I migration of an embedded protoplanet, as was noted by PM06.

Numerical study of an isolated corotation resonance
We check in this section the analytical predictions of section 3 by means of numerical simulations involving an isolated corotation resonance (hereafter CR).

Numerical issues
Our setup offers a number of similarities with the setup of Masset & Ogilvie (2004) for the case of an isothermal disk. The hydrodynamics equations for the disk described in section 3.1 are solved using the code Fargo. A description of the properties of this code is deferred to section 5.1, in which the code is used to simulate an embedded planet. As in Masset & Ogilvie (2004), we deal with the m = 3 CR. The disk is therefore torqued by an m = 3 external potential Φ that reads where φ(r) denotes the radial profile of the potential, Ω p its pattern speed (note that we work in the corotating frame), t is the time and where is a temporal tapering that turns on the potential on the timescale τ .
The total torque Γ c exerted by the disk on the perturber, given by Eq. (21), is evaluated by where N r (N s ) is the radial (azimuthal) number of zones of the mesh, S i,j is the surface area of zone (i, j), Φ i,j and Σ i,j are the external potential and surface density at the center of this zone, and ∆ϕ = 2π/N s is the azimuthal resolution. Furthermore, the contribution Γ c,1 of the function Ψ to the torque is obtained by substituting Σ 1 by p 1 /c 2 s in Eq. (21). It is therefore evaluated by where p i,j and c si,j are the pressure and sound speed at the center of zone (i, j). The contribution Γ c,2 of the entropy perturbation to the torque is eventually estimated as follows: The radial computational domain is narrow enough to avoid the location of the m = 3 inner and outer Lindblad resonances (see Masset & Ogilvie 2004). Despite this precaution, wave killing zones next to the boundaries were implemented to minimize unphysical wave reflections (de Val-Borro et al. 2006). Furthermore, the torque evaluation is performed by summing only on a domain of the grid that does not contain the wave killing zones, and the summation includes a spatial tapering on the edges of that domain.
The disk surface density and temperature are initially axisymmetric with power-law profiles: and where Σ c and T c are the surface density and temperature at the corotation radius r c , and where f is the flaring index of the disk. The disk aspect ratio is given by h is the disk scale height at radius r. A vanishing value of the flaring index f therefore corresponds to a uniform disk aspect ratio. The functions V and S are constant and read: The main numerical parameters are those taken by Masset & Ogilvie (2004), namely a h(r c ) = 0.01 disk aspect ratio at corotation, and Σ c = 1. Our disk is inviscid. The libration islands are resolved by 30 zones azimuthally. As the potential increases, the radial width of the islands also increases. Their maximal radial width W spans approximately 20 zones.
The results presented in next section have the following units: the mass of the central object M * is the mass unit, the corotation radius r c of our m = 3 CR is the distance unit and the Keplerian orbital period T orb at r = r c is 2π times the time unit.
This case, that we call the "flat potential case", has the following parameters: σ = 2 and f = −0.3, which implies, from Eqs. (48) and (49), that V = −0.5 and S ≈ −0.57, 2. A potential profile that corresponds to the m = 3 Fourier component of the smoothed potential of a M = 3.1 × 10 −6 M * point-like object. The softening length is ε = H(r c ), which is approximately equal to W . The object rotates at speed Ω p , with orbital radius r c . This neglects the pressure gradient effects, as we do not resolve the distance from orbit to corotation, but this distance is much smaller than the potential softening length, so this is not a concern in the present case. By contrast to the previous case, we call this situation the "peaked potential case". The The results shown are obtained with an adiabatic calculation, except in the close-up, where we compare the isothermal and adiabatic corotation torques over the whole duration of the calculations. Numerical results are displayed with a symbol while the theoretical expectations are displayed with curves. We plot as a function of time the adiabatic corotation torque (diamonds and solid curve), the contribution of the function Ψ to the torque (stars and dashed curve), and the contribution of the entropy perturbation (triangles and dot-dashed curve). The long-dashed curve, which is nearly superimposed to the solid curve, displays the corotation torque expression of Goldreich & Tremaine (1979). The vertical solid line gives an estimate of the final libration time (see text).
3. A calculation similar to the previous one, except that σ = 0.5 and f = −0.1, so that V = 1 and S ≈ −0.71. The results are depicted in Fig. 2b.
For the three pairs (V, S) quoted above, the tapering timescale value is τ = 150 T orb , which corresponds to the duration of the calculations. This is about three times larger than the final libration time, estimated as In each case we evaluate: • the total corotation torque (diamonds) with Eq. (43), to be compared to the analytical expression (solid curve) given by Eq. (41). In our units, Γ 0 ≈ 39.5, • the contribution of the function Ψ to the corotation torque (stars) obtained with Eq. (44), the expected expression of which (dashed curve) is calculated using Eq. (39), • the contribution of the entropy perturbation to the torque (triangles) using Eq. (45), which is to compare to the prediction of Eq. (40), represented by the dot-dashed curve.
In these figures, the corotation torque first increases with time since the potential is progressively turned on until it reaches its final value at the end of the calculation. After some time it starts to oscillate. This oscillation corresponds to the saturation of the CR, as the ratio t/T lib tends to unity (Ogilvie & Lubow 2003). Figs. 1, 2a and 2b therefore argue that our numerical simulations succeed in reproducing the results of our analytical study as long as t T lib , that is when a linear analysis is grounded (which requires that the time of the calculation be much smaller than the libration time).
The examination of the results of these calculations leads to the following comments: • In the flat potential case, depicted in Fig. 1 the calculation, where ψ(r) denotes the radial profile of Ψ. This situation therefore corresponds to a cold case. As expected from Eq. (32), the analytical corotation torque and the expression of Goldreich & Tremaine (1979) almost coincide. The closeup shows the torque evolution over the whole extent of the calculation, up to t = τ . The torque obtained with a locally isothermal equation of state is also depicted. Our isothermal runs have same radial temperature dependence as the adiabatic runs (see Eq. (47). Although there is an entropy gradient in these isothermal calculations, it does not contribute to the corotation torque as it would in an adiabatic disk: the appearance of the singular contribution at corotation in the adiabatic case is linked (i) to the advection of entropy and (ii) to the appearance of a singularity in the perturbed density and temperature fields. In the isothermal situation, neither the entropy is conserved along a fluid element path, nor is a temperature singularity allowed to appear. The comparison of isothermal and adiabatic calculations shows that, as expected for a cold case, the adiabatic and isothermal torques coincide, as long as we are in the linear regime. We note that both torques do not oscillate about 0 since the potential reaches a stationary value only at the end of the calculation.
• For the two calculations of the peaked potential case, depicted in Figs. 2a and 2b, we find that ℜ[ψ(r c )] ≈ −0.2φ(r c ). Thus, the term | − Φℜ(Φ + Ψ)| slightly dominates the term |Φ + Ψ| 2 in Eq. (41). Because S < 0 for these calculations, the corotation torque in the adiabatic case (diamonds and solid curve) is larger than the corotation torque in an isothermal disk (long dashed curve) with the same parameters, as predicted by Tanaka et al. (2002). In particular, in the case for which V = 0, the isothermal corotation torque vanishes, while we find a net, positive corotation torque for an adiabatic flow, in correct agreement with the analytical expression.

Dynamics of the corotation region
We discuss in this section the dynamics of the corotation resonance of an adiabatic disk and give some comments about the corotation torque expression of Eq. (41).
In the isothermal case, the corotation torque expression involves the product of the gradient of vortensity and the square of the effective potential (Φ + η), see e.g. Tanaka et al. (2002). The torque is then given by the angular momentum budget between material flowing outwards and material flowing inwards at corotation, regardless of the sign of Φ+ η. Eq. (41) displays a term that has a similar behavior, except that it does not feature the vortensity gradient only, but rather V + 2S. This factor scales with the (logarithmic) gradient of (Σ 0 /B)S 2/γ , which is a key quantity considered by Li et al. (2000) and by Lovelace et al. (1999), who pointed out that vortensity is not conserved in a two-dimensional adiabatic flow.
In addition to this term, Eq. (41) contains a term that scales with Φ[Φ + ℜ(Ψ)]. The sign of this term therefore depends on the relative signs of Φ and Φ + ℜ(Ψ). In order to get some insight into the physical meaning of this term, we show at Fig. 3 the response of the disk in the corotation region, for the entropy and the surface density. These fields correspond to the calculation with the flat potential profile considered at the previous section. The disk has a negative radial entropy gradient. Therefore, libration brings the (larger) inner entropy to the outer part of the libration islands, yielding a positive perturbed entropy (brighter areas), while it brings the (smaller) outer entropy to the inner part of the libration islands, yielding a negative perturbed entropy (darker areas). An opposite behavior is observed for the perturbed density, since the relative pressure perturbation (not represented) is much smaller.
The sign of this torque component can be understood as follows. Fig. 4 depicts the situation in two cases: Φ and Φ + ℜ(Ψ) have same sign (left), and Φ and Φ + ℜ(Ψ) have opposite signs (right). In the left case, the negative perturbed surface density on the outside of corotation is located in the region where ∂ϕΦ < 0, hence the perturbation yields a positive torque on the perturber. A similar conclusion applies to the material flowing inwards which has positive perturbation of surface density. The torque on the perturber is therefore positive, in agreement with the sign of −SΦ[Φ + ℜ(Ψ)]. An opposite conclusion holds for the case where Φ[Φ + ℜ(Ψ)] < 0.
The order of magnitude and functional dependence of this torque component can be justified as follows. As the sign has been justified at the previous paragraph, we give here an estimate of the absolute value. The perturbed surface density on the outside of corotation is ∼ |SΣ 0 δ/r c |, where δ = [(Φ + ℜ(Ψ))/(−8AB)] 1/2 is an order of magnitude of the width of the libration islands. The specific torque in the region of surface density perturbation is ∼ |mΦ|, while the area covered by the perturbation of surface density scales with r 2 c δ. The torque arising from this region therefore scales with |mr c δ 2 SΦΣ 0 |, which is exactly the scaling of |Γ 0 SΦ[Φ + ℜ(Ψ)]|, within a numerical factor in O(1).
The singular behavior of this torque component, which stems from Eq. (40), and which appears as a Dirac's delta function at corotation, can be understood as follows: as the strength of the perturbation decreases, the width of the libration islands tends to zero, while the libration time tends to infinity (libration disappears), hence we are left, in the linear regime, with a torque contribution that comes strictly from the corotation radius and therefore appears as singular.
It is worth noting that only half of the second term of Eq. (41) comes from Eq. (40). Eq. (39), which is obtained from the momentum flux jump, and which as such captures effects occurring at a finite (albeit small) distance from corotation, also displays a term similar to that of Eq. (40). The advection of entropy perturbations is not a silent process: it triggers the emission of pressure waves (Foglizzo & Tagger 2000). Our torque expression indicates that half of the energy required to advect entropy in the libration islands is evacuated through pressure waves.

Application to the case of an embedded protoplanet
In section 3, we derived an expression for the corotation torque between a radiatively inefficient disk and an external rotating potential. This expression is successfully reproduced by local numerical simulations of an isolated corotation resonance, in the linear regime. We now contemplate the case of an embedded protoplanet in a radiatively inefficient twodimensional disk, for which all co-orbital corotation resonances are simultaneously active.

Numerical features and setup
Our numerical simulations are performed with the code Fargo. It is a staggered mesh hydrocode that solves the Navier-Stokes, continuity and energy equations on a polar grid. It uses an upwind transport scheme with a harmonic, second-order slope limiter (van Leer 1977). Its particularity is to use a change of rotating frame on each ring of the polar grid, which increases the timestep significantly (Masset 2000a,b), thereby lowering the computational cost of a given calculation. The energy equation that we implemented in Fargo is: where e is the thermal energy density, v = (u, rΩ) T denotes the flow velocity, p is the vertically integrated pressure and Q is a heating source term that accounts for the disk viscosity (see e.g. D' Angelo et al. 2003). The energy equation solver is implemented as in Stone & Norman (1992).
In this work, the disk is taken inviscid so Q = 0. There is no radiative transfer either, since the disk is assumed to be radiatively inefficient. Furthermore, p and e are connected by an ideal equation of state p = (γ − 1)e, where the adiabatic index γ is set to 1.4. This equation of state can be expressed in terms of the disk temperature T and surface density Σ as p = ΣT . The adiabatic sound speed reads c s = √ γT , hence c s = √ γ c s,iso , where c s,iso refers to the isothermal sound speed. We comment that the Lindblad torque, which scales as c −2 s (Ward 1997), is therefore weakened by a factor of γ in an adiabatic disk. The same is true of the corotation torque, when there is no entropy gradient. We checked both effects with appropriate calculations, not reproduced here. This plays in favor of a total torque reversal in adiabatic disks with a negative entropy gradient.
The disk is initially slightly sub-Keplerian (the pressure gradient is accounted for in the centrifugal balance), axisymmetric, with power-law profiles for the surface density and temperature given by Eqs. (46) and (47).
For a comparative purpose, calculations involving a locally isothermal equation of state are performed. In isothermal calculations, no energy equation is solved: p and Σ are simply connected by p = Σc 2 s,iso . These isothermal calculations have same initial surface density and temperature profiles as the adiabatic runs.
The disk is perturbed by the smoothed potential of a protoplanet. We adopt a Plummer potential, with a softening length ε = 0.6H(r p ) (unless otherwise stated), r p being the planet orbital radius. This fiducial value is quite substantial for our purposes, but investigating the disk response at much smaller softening lengths, where the adiabatic effects on the corotation torque are increasingly important, requires a very large resolution. A high resolution systematic study at small softening length will be presented in a forthcoming work.
The protoplanet is held on a fixed circular orbit, at r = r p . The disk parameters are summed up in Table 1, where they are expressed in the following unit system: r p is the length unit, the mass of the central object M * is the mass unit and (GM * /r p 3 ) −1/2 is the time unit, G being the gravitational constant (G = 1 in our unit system). We denote by T orb the planet orbital period, M p the planet mass and q = M p /M * the planet to primary mass ratio.

An illustrative example
We show the results of an illustrative calculation with a q = 2.2 × 10 −5 planet to primary mass ratio (corresponding to M p = 7.3 M ⊕ if the central object has a solar mass). The horseshoe libration time is where Ω p is the protoplanet angular velocity and x s denotes the half-width of the horseshoe region. Masset et al. (2006) have given an estimate of x s in the isothermal case, that reads x s ≈ 1.16r p q/h(r p ). A streamline analysis was performed and confirmed that this estimate holds for an adiabatic disk, if one substitutes h(r p ) with √ γh(r p ). We find therefore τ lib ≈ 60 T orb . Numerical diffusion eventually alters the conservation of entropy. Nevertheless, the horseshoe region spans 20 zones radially, which is sufficient to follow the horseshoe dynamics over several libration times. Since we are concerned here with a fraction of the libration time, the entropy is conserved with a good level of accuracy over the duration of our runs, and it can be regarded as a Lagrangian tracer of the flow.
Two calculations were performed: an adiabatic and an isothermal one. Both lasted thirty orbital periods, hence half the horseshoe libration time. This calculation has σ = 0.5 and f = 0, as in PM06. This gives S ≈ −0.57. Fig. 5 displays the gas entropy, surface density and pressure obtained in the adiabatic calculation, after 15 T orb . Each field represents the relative perturbation of the corresponding quantity with respect to the unperturbed state. For instance, the top right panel shows [Σ(r, ϕ)−Σ 0 (r)]/Σ 0 (r). While the azimuthal range spans the whole [0, 2π] interval, the radial range depicted is restricted to a band of width 2.5x s around the corotation radius r c . We overplot streamlines to the entropy panel to give an idea of the extent of the horseshoe region. The vertical dashed line represents the corotation radius. Whereas the pressure panel does not display any significant perturbation, the entropy and density panels show the propagation of a perturbation inside the horseshoe region, which slides along the separatrices. This is reminiscent of the behavior commented in the case of an isolated resonance at section 4.3.
The interpretation of this dynamics is as follows: the entropy of the fluid elements is conserved as they perform a horseshoe U-turn in the co-orbital region. When there is initially an entropy gradient at corotation, the co-orbital dynamics yields an entropy perturbation that has a sign opposite of that of the entropy gradient on the outwards U-turns, and the sign of the entropy gradient on the inwards U-turns. Since the pressure field is only weakly perturbed, the entropy perturbation is related to a density perturbation of opposite sign and, in relative value, of same order of magnitude. Therefore, if there is a negative entropy gradient at corotation (S < 0, as in the example shown here), the co-orbital dynamics yields a negative density perturbation at ϕ < ϕ p and a positive density perturbation at ϕ > ϕ p , with straightforward consequences for the corotation torque. Using an expression inherited from the terminology of Riemann solvers, we call this perturbation a contact discontinuity. A contact discontinuity is characterized by a discontinuity in the density and temperature fields, while the pressure and velocity fields are continuous. A contact discontinuity is simply advected by the flow. Here it follows the horseshoe dynamics, and it remains confined to the horseshoe region.
We give hereafter a simple estimate of the relative perturbation of the disk surface density due to the advection of entropy. We consider a fluid element that performs a horseshoe U-turn from the inner part of the horseshoe region (where we assume that there is no entropy perturbation, which is true as long as t < τ lib /2) to the outer part. All physical quantities at the inner (outer) leg of the horseshoe streamline are denoted by a minus (plus) subscript. A first-order expansion yields, assuming no pressure perturbation: where 0 < x < x s is the distance of the streamline to corotation, and: On the outer horseshoe leg, the disk surface density is perturbed according to the entropy perturbation and reads: where R is the relative perturbation of surface density at r = r c + x (we assume a symmetric horseshoe U-turn), due to the entropy advection. Entropy conservation along the fluid element path (S − = S + ) leads to: The horseshoe U-turn that we have considered lags the planet (ϕ < ϕ p ). A similar conclusion holds for a horseshoe U-turn that switches from the outer leg to the inner one (at ϕ > ϕ p ), hence we finally have: The bottom right panel of Fig. 5 displays the slices of the perturbed density field at t = 15 T orb , for ϕ − ϕ p = 1 (diamonds) and ϕ − ϕ p = −1 (stars). The two horizontal dashed lines display the values of R(−x s ) and R(x s ), where x s is estimated through a streamline analysis. Similarly, the long-dashed curve shows R(x) = 2xS/r c , which is in correct agreement with the calculation results. The surface density structure in the horseshoe region is therefore dictated by the sign of S. In particular, we do not expect any contact discontinuity in the homentropic case (S = 0). We have checked this prediction with a numerical simulation (not presented here).

Excess of corotation torque and entropy gradient
An order of magnitude of the excess of corotation torque arising from the perturbation of the surface density field can be given by estimating how the standard horseshoe drag expression (Ward 1991;Masset 2001) is modified by the perturbation of surface density R(x)Σ 0 (r c ). We consider the outwards horseshoe U-turns that occur at ϕ < ϕ p . Assuming, in this order of magnitude estimate, that the rotation profile of the disk is unperturbed, we evaluate the variation of angular momentum flux of the horseshoe disk material after the U-turn attributable to the change of the disk's surface density: where j c is the specific angular momentum of the material at corotation. The first factor of the integrand of Eq. (58) represents the material velocity in the corotating frame, due to the shear. The last factor is the material specific angular momentum obtained from a first order expansion at corotation. Similarly, we obtain the change of angular momentum flux due to the perturbation of surface density on inwards horseshoe U-turns: Adding Eqs. (58) and (59), we are left with: Fig. 6 shows the excess of corotation torque between an adiabatic and isothermal calculation with same parameters, as a function of the half-width of the horseshoe region. This excess is obtained by subtracting the total torque of an adiabatic and an isothermal calculation (the isothermal torque being rescaled by a factor γ −1 , since c s = √ γ c s,iso ). We call this difference the torque excess for further reference. Each data point corresponds to a calculation with a given planet mass, for which we determine x s through a streamline analysis. We find that the torque excess approximately scales as x 4 s , and that it is within a factor 2 of our order of magnitude estimate, given by −∆Γ HS .
The torque expression of Eq. (41) as well as the horseshoe drag expression of Eq. (60) suggest that the torque excess scales with S, hence with the entropy gradient. In order to check that, we have undertaken a number of calculations with different values of S. These calculations have q = 2.2 × 10 −5 , and the disk parameters are those of Table 1. Each entropy gradient is realized with different combinations of the indexes of the pressure and surface density power laws. Adiabatic effects on the torque are assessed in two different ways: 1. By calculating the torque excess, as in Fig. 6. 2. By evaluating the following integral: which provides an estimate of the torque due to the contact discontinuity (this contribution arises from perturbations of Σ which do not have a pressure counterpart). In the linear regime, Eq. (61) amounts to a summation over m of the last term of Eq. (24). We shall check this statement in the next section.
These two estimates of adiabatic effects on the torque value are shown respectively in Figs. 7a and 7b. Remarkably, they coincide within ∼ 25 %. We will comment further this coincidence in the next section.
The main conclusion that can be drawn from the results of Fig. 7 is that the torque excess (or the contact discontinuity contribution) essentially depends on the entropy gradient, as expected. The excess is positive for a negative entropy gradient, hence we may expect the total torque exerted on a planet embedded in a radiatively inefficient disk to be a positive quantity if the radial entropy gradient is sufficiently negative.

Connection to the analytical expression
We have given at Eq. (40) an estimate of the singular torque contribution from the contact discontinuity at an isolated resonance, while we have estimated the total contribution in the planetary case of the contact discontinuity using Eq. (61) at section 5.2.2. We check in the present section that this total contribution corresponds to the sum over m of the torque expression of Eq. (40). For this purpose, we have adopted a planet to primary mass ratio q = 5 × 10 −6 , as the one adopted in the previous sections (q = 2.2 × 10 −5 ) led to poor agreement, presumably because of the onset of non-linear effects. For each azimuthal wavenumber m, we measure ℜ(Ψ m ) from the calculation output (at t = 5 T orb ), and we evaluate the sum over m of the torque Γ c,m,2 : where: is the partial sum of Γ c,m,2 . We compare the torque contribution given by Eq. (61) to Γ ∞ . The results are presented in Fig. 8. The agreement between the direct torque measurement and the linear estimate is excellent.
This confirms what we anticipated in section 3.2, and what is shown in appendix A, that the contribution of the last term of Eq. (25) to the torque is negligible in the planetary context. Also of interest is the torque density associated respectively to p/c 2 s and Σ − p/c 2 s . The sum of these two torque densities is the total torque density. They are represented at Fig. 9. The total torque density displays a smooth profile and a narrow peak at corotation. This is reminiscent of the torque density found by PM06 (their Fig. 2) or by Morohoshi & Tanaka (2003) (their Fig. 3). The decomposition above splits this total torque density in a smooth component arising from p/c 2 s , which reminds the torque density in an isothermal disk, and a sharp, localized torque density arising from Σ − p/c 2 s . This corresponds to the torque density of the contact discontinuity contribution given by Eq. (61). Fig. 9 shows that this contribution (which is singular at corotation in the linear case for an isolated resonance) is here bounded by the extent of the horseshoe region.
We comment the surprising agreement found at the previous section between the torque excess and the contribution of the contact discontinuity. The linear analysis suggests that the former should be the sum of Γ 0,m [2S|Φ m + Ψ m | 2 − 2SΦ m (Φ m + ℜ(Ψ m ))], which, in the limit where |Φ m + ℜ(Ψ m )| ≪ |ℜ(Ψ m )| and |Φ m + ℜ(Ψ m )| ≪ |Φ m |, should reduce to −2Γ 0,m SΦ m (Φ m + ℜ(Ψ m )), that is twice the contribution of the contact discontinuity (see section 3.2.2). Nevertheless, for the calculations presented here, we can check that 2 m m|Φ m + Ψ m | 2 is almost exactly compensated by m mΦ m (Φ m + ℜ(Ψ m )). Namely the ratio of the former to the latter quantity is found to be 1.07, which explains why the full excess expression essentially amounts to the contact discontinuity contribution. Presumably this coincidence is fortuitous and linked to the relatively large softening length that we use. At smaller softening length, the term in Φ m (Φ m + ℜ(Ψ m )) should largely dominate, yielding a ratio of 2 between the torque excess and the contribution of the contact discontinuity. We note that PM06 also quote that the torque estimate given by their equation (1) accounts for the total torque within 25 % (this equation can also be seen as an estimate of the contact discontinuity contribution). This seems to suggest that the softening length of 0.6H(r p ) that we adopted is a correct choice to reproduce the magnitude of the corotational effects in adiabatic three-dimensional disks.

Softening length
In an isothermal disk, the corotation torque scales with |Φ + η| 2 (Tanaka et al. 2002). Even if Φ diverges at corotation, Φ + η remains finite, which makes the isothermal corotation torque a quantity relatively insensitive to the softening length. The situation is quite different for the effects linked to the entropy advection that we present here: they involve the product Φ(Φ + Ψ), which diverges when Φ does. Adiabatic effects on the corotation torque should acquire a very large magnitude at small softening length. In particular, if the softening length is smaller than the distance from orbit to corotation, the magnitude of these effects should strongly depend on this distance, which scales with the pressure gradient. If one regards the softening length as a proxy for the altitude in a three-dimensional disk, the extent of the disk vertical scale height concerned by these very small softening length issues should be small, however, since the distance from orbit to corotation is a fraction of r p h 2 . Nevertheless, it is of interest to investigate the behavior of the corotation torque in an adiabatic flow at very small softening length to assess the importance of such effects. Owing to the very large resolution required to investigate this problem, we defer this investigation to a forthcoming work.

Saturation
The origin of the effects presented here is the advection of entropy in the corotation region, that triggers an entropy perturbation (and therefore a density perturbation) whenever there is an entropy gradient in the equilibrium profile. Libration occurs on different timescales for the different streamlines of the corotation region, which tends to stir the entropy and to flatten out the entropy profile across the corotation region (be it the horseshoe region in the planetary case or a libration island in the isolated resonance case). This is quite similar to the behavior of the corotation torque in an isothermal disk, which tends to saturate because the vortensity profile is flattened out by libration. In this case, it is the viscous diffusion which can prevent the flattening out of the profile if it acts sufficiently rapidly to establish the large scale gradients before a libration time. This has been studied for an isolated resonance by Goldreich & Sari (2003) and Ogilvie & Lubow (2003), and by Balmforth & Korycansky (2001) and Masset (2001) for a planetary co-orbital region. In both cases, the degree of saturation of the corotation torque in steady state depends on the ratio of the libration time and of the viscous time across the libration region. The dissipative processes required to prevent the torque saturation in the situation presented here should be able to impose the large scale entropy gradient over the corotation region in less than a libration time. Radiative processes (cooling and heating) should therefore occur on a timescale longer than a horseshoe U-turn (otherwise the flow can rather be considered as locally isothermal), but they should act on a timescale shorter than the libration time.
We provide an estimate of the horseshoe U-turn time and of the libration time for a small mass object embedded in a gaseous disk. The horseshoe half-width x s is ∼ r p q/h(r p ). Neglecting pressure effects and writing a simplified Jacobi constant for a test particle near a horseshoe U-turn as: J = −GM p /(2Br p |ϕ − ϕ p |) + A(r − r p ) 2 , we can estimate the distance of closest approach between the planet and a test particle flowing along a horseshoe separatrix as r p |∆ϕ| s = Ω 2 p H(r p )/2|AB| = O(H(r p )). The time required to perform a horseshoe U-turn can be deduced using the radial drift velocity of the test particle when it crosses the orbit, at its closest approach from the planet:ẋ = GM p /(2Br 2 p ∆ϕ 2 s ). That yields: where R H = r p (q/3) 1/3 is the Hill radius of the planet, and where the last equality holds for a Keplerian disk. When the planet emerges from the disk (H(r p ) ∼ R H ), the horseshoe U-turn occurs on the dynamical timescale. When dealing with an embedded object however (R H < H(r p )), the horseshoe U-turn time can be substantially longer than the dynamical time (e.g. 10 times longer for an Earth mass object embedded in a disk with h(r p ) = 0.05).
Using Eq. (52), we are led to: There is at least an order of magnitude difference between the horseshoe U-turn time and the libration time in a thin disk, hence it should be possible to find a location in the disk where the cooling time is much longer than the U-turn time and yet shorter than the libration time, so as to maintain an unsaturated corotation torque.

Entropy gradient and baroclinic instability
The effect that we present in this two-dimensional analysis is of particular interest when there is a negative entropy gradient at corotation, since this may suffice to halt type I migration. It would be of interest to generalize the present analysis to the case of a three-dimensional baroclinic disk. We comment also that in such systems, a negative entropy gradient may render the disk unstable to a baroclinic instability (Klahr 2004;Klahr & Bodenheimer 2003). It is certainly important to examine the interplay of the baroclinic instability and of the corotational effects presented here. The turbulence generated by the baroclinic instability, in particular, could provide a mechanism to prevent the saturation of the corotation torque, much like the turbulence arising from the MRI can prevent the corotation torque saturation in an isothermal disk.

Conclusions
We evaluate the corotation torque between an adiabatic gaseous disk and a uniformly rotating external potential. In the linear case for an isolated resonance, we find a singular contribution at corotation which scales with the entropy gradient, and which arises from the advection of entropy within the libration region. This effect neither exists in isothermal or locally isothermal flows, nor does it exist for barotropic fluids (such as fluids described by a polytropic equation of state). We provide a torque expression at an isolated resonance which involves the pressure perturbation at corotation. We then check the torque expression by two-dimensional adiabatic calculations that involve an isolated resonance. In particular, we exhibit a case with a flat vortensity profile, for which the corotation torque does not cancel out and is in correct agreement with the analytical expression. We then turn to the case of an embedded planet, for which we find an excess of corotation torque in the adiabatic case, which scales with the entropy gradient. For a sufficiently small planet mass, we check that this excess can be accounted for by a summation over the resonances of the torque excess that we found in the first part. This confirms that this effect is essentially a linear effect. We finally discuss in section 6 some open questions linked to the softening length, to the saturation, to the case of a three-dimensional baroclinic disk, and to the interplay with the baroclinic instability, on to which theoretical efforts should focus in a nearby future.
We wish to thank Sijme-Jan Paardekooper and John C.B. Papaloizou for interesting discussions on the topics covered in this manuscript. We also thank Alessandro Morbidelli for a thorough reading of a first version of this manuscript, and an anonymous referee for comments that led to an improvement of the paper.

A. Additional contribution to the corotation torque
We check hereafter that the contribution of the last term of Eq. (25) is negligible for an embedded planet. For this purpose, we compare G = dx ℑ[Ψ(x)]Φ(x)/x to −π[ℜ(Ψ)Φ] rc . Fig. 10a shows the m = 8 component of Ψ = p/Σ 0 for the calculation presented in sec-tion 5.2.3. We clearly see that the behavior of ℑ(Ψ) in the vicinity of corotation comes from the overlap of the behavior arising at the inner and outer Lindblad resonances. Between its outermost inner minimum at r − ∼ 0.87 and its innermost outer maximum at r + ∼ 1.12, ℑ(Ψ) can be considered as having a linear dependence in x. Notwithstanding the decrease of Φ as one recedes from corotation, the main contribution to G will come from ℑ[Ψ(x)]/x between these two radii, as it exhibits a flat behavior over this range. This yields G ∼ 2Φ(r c )|ℑ[Ψ(r ± )]|, that it to say a result comparable in order of magnitude to −π[ℜ(Ψ)Φ] rc . Nevertheless, the final contribution of the extra term is much smaller than the singular one at corotation for the following reasons: 1. The potential decreases sharply as one recedes from corotation, which provides a cut-off to the extra term, that is not localized at corotation.
2. The extra term is partially compensated for by the first term of Eq. (10), which we have neglected in writing Eq. (24), and which yields another term in Eq. (25) that reads −(F S/r 2 Ω) rc dℑ[Ψ(x)]/dx. Adding this additional term and the last term of Eq. (25), we are left with: which shows that, in addition to the potential cutoff, the last term of Eq. (25) should be decreased by a factor of 4 (we neglect, at this level of accuracy, the jump in ℑ(Ψ ′ ) at corotation).
We have checked on the calculation presented in section 5.2.3 that the contribution of these extra terms is indeed small compared to the singular contribution at corotation at all m. This is shown in Fig. 10b, from which we can conclude that the total contribution of the extra terms is about an order of magnitude smaller than the singular contribution. The agreement that we found in section 5.2.3 between the numerical simulation of an embedded planet and the torque series, which was of the order of a percent, might then be fortuitous. Nevertheless, we expect an agreement of the order of 10 %, still very satisfactory. These findings are also compatible with the fact that we hardly see any diffuse torque density outside of the horseshoe region in Fig. 9b. [-] [+] [-] [+] [-] [+] [-] [+] [-] -- . The minima and maxima of Φ are indicated at the left, while the minima and maxima of Φ + ℜ(Ψ) are indicated at the right. In the corotation region, material librates about the maxima of the effective potential Φ + ℜ(Ψ). We assume a negative entropy gradient, hence material flowing outwards has a negative perturbed surface density, while material flowing inwards has a positive perturbed surface density, as indicated by the minus and plus signs. of the gas entropy, surface density and pressure, at t = 15 T orb ≈ τ lib /4. The protoplanet is located in r = r p , ϕ = ϕ p . In the top left panel, streamlines are overplotted and the vertical dashed line stands for the corotation radius r c . In the top right and bottom left panels, the color scale is adjusted to highlight the advection of the entropy perturbation (see text). The nearly horizontal overdensity structure at ϕ = ϕ p is the protoplanet's wake. Bottom right: slices of the relative perturbed density field at the same time, at ϕ−ϕ p = 1 (diamonds) and ϕ−ϕ p = −1 (stars). The two horizontal dashed lines refer to the values of R(−x s ) and R(x s ), while the long-dashed curve displays the quantity 2(r − r c )S/r c (see text and Eq. (57)).