Torque wiggles – a robust feature of the global disc–planet interaction Open Access

ABSTRACT

Gravitational coupling between planets and protoplanetary discs is responsible for many important phenomena such as planet migration and gap formation. The key quantitative characteristic of this coupling is the excitation torque density – the torque (per unit radius) imparted on the disc by planetary gravity. Recent global simulations and linear calculations found an intricate pattern of low-amplitude, quasi-periodic oscillations in the global radial distribution of torque density in the outer disc, which we call torque wiggles. Here, we show that torque wiggles are a robust outcome of global disc–planet interaction and exist despite the variation of disc parameters and thermodynamic assumptions (including β-cooling). They result from coupling of the planetary potential to the planet-driven density wave freely propagating in the disc. We developed analytical theory of this phenomenon based on approximate self-similarity of the planet-driven density waves in the outer disc. We used it, together with linear calculations and simulations, to show that (a) the radial periodicity of the wiggles is determined by the global shape of the planet-driven density wave (its wrapping in the disc) and (b) the sharp features in the torque density distribution result from constructive interference of different azimuthal (Fourier) torque contributions at radii where the planetary wake crosses the star–planet line. In the linear regime, the torque wiggles represent a weak effect, affecting the total (integrated) torque by only a few per cent. However, their significance increases in the non-linear regime, when a gap (or a cavity) forms around the perturber’s orbit.

1 INTRODUCTION

Gravitational coupling between a gaseous disc and a perturbing mass (a planet, a satellite, or a binary companion) has been actively explored starting with the seminal studies of Lin & Papaloizou (1979) and Goldreich & Tremaine (1980, hereafter GT80). This tidal interaction is recognized as leading to a number of important effects such as gap formation (Lin & Papaloizou 1986), planet migration, planetary/binary eccentricity evolution (GT80), disc eccentricity excitation (Lubow 1991), etc.

A key ingredient of the disc–planet¹ interaction is the excitation of the density waves in the disc–pressure (sound) waves modified by the differential rotation of the disc. These waves gain energy and angular momentum due to the gravitational potential of the perturber at their launching sites (Lindblad resonances, GT80), propagate through the disc, and eventually deposit their angular momentum and energy to the disc material once they are damped. This makes disc–planet coupling an inherently non-local process (Lunine & Stevenson 1982; Greenberg 1983; Goldreich & Nicholson 1989; Goodman & Rafikov 2001; Petrovich & Rafikov 2012; Rafikov & Petrovich 2012). Moreover, the gravity of the perturber acts on the density perturbation due to the density wave resulting in the angular momentum exchange between the perturber and the wave (and eventually the disc).

One of the main quantitative characteristics of disc–planet coupling is the so-called excitation torque density dT_ex/dR, which is defined in polar cylindrical coordinates r = (R, ϕ, z) as the z-component of the torque exerted by the (Newtonian) potential of the planet per unit radial distance in the disc. One can also interpret dT_ex/dR as the amount of angular momentum added to the density wave per unit radial distance and per unit time by the planetary potential. The integral of dT_ex/dR over the full extent of the disc determines, via Newton’s third law, the evolution of the planetary angular momentum, which manifests itself as planetary migration (GT80). The radial profile of dT_ex/dR is one of the key inputs (together with the wave damping mechanism, see Goodman & Rafikov 2001; Rafikov 2002a) determining the amplitude evolution of the density waves as they travel away from the planetary orbit. Its knowledge is also important for determining the structure of planetary gaps (Rafikov 2002b). Thus, a detailed understanding of the behaviour of dT_ex/dR is of paramount importance for obtaining a complete picture of disc–planet interaction.

Some of the key features of dT_ex/dR behaviour in a disc which is roughly uniform near the planet (i.e. with no gap formed) have been predicted already in GT80, with further refinements provided by, e.g. Meyer-Vernet & Sicardy (1987), Artymowicz (1993), and Ward (1997). In particular, they have shown using linear theory that dT_ex/dR should rapidly decrease as |R − R_p| ≲ H_p (were R_p is the semimajor of the planet on circular orbit and H_p is the disc scale height at R_p), an effect that is known as ‘torque cut-off’. The overall shape of dT_ex/dR near the planet is set by the tidal disc–planet coupling at numerous Lindblad resonances, with the largest contribution coming from resonances located at |R − R_p| ∼ H_p. This is indeed what one sees in Fig. 1, in which we show a typical run of the excitation torque density obtained using our linear calculations (see Section 2.2 for details). One can see prominent peak and trough² of dT_ex/dR located just outside and inside of the planetary orbit, at |R − R_p| ∼ H_p. The asymmetry in their magnitudes is the reason behind planet migration (GT80). Moreover, both Lin & Papaloizou (1979) and GT80 predicted that the excitation torque density should monotonically decay as | dT_ext/dR | ∝ |R − R_p|⁻⁴ for |R − R_p| ≳ H_p.

Some of these statements have been refined recently. In particular, the two-dimensional (2D) hydrodynamic simulations of Dong et al. (2011) carried out in the local (sharing sheet) approximation found that dT_ex/dR does not decrease monotonically but actually reverses its sign around |R − R_p| ≈ 3.2H_p. This ‘negative torque density’ phenomenon is illustrated in the inset in Fig. 1(b), where dT_ex/dR first turns negative outside the main positive peak in the outer disc, R > R_p (in this case it is less pronounced in the inner disc), marked with an arrow (annotated ‘NT’). This phenomenon was explained by Rafikov & Petrovich (2012) as resulting from interference of the density wave contributions launched at different Lindblad resonances. They have also demonstrated analytically that the decay of dT_ex/dR beyond this point follows | dT_ext/dR | ∝ |R − R_p|⁻⁴ scaling but with a coefficient different in both the sign and magnitude from the prediction of GT80. This torque reversal was subsequently observed in the global fully non-linear simulations of Duffell & MacFadyen (2012) and Kley et al. (2012). It has been missed in the semi-analytic studies (GT80; Ward 1997; Dempsey, Lee & Lithwick 2020), which utilized only the amplitude information for the Fourier components of the torque in constructing the torque density profiles, whereas Rafikov & Petrovich (2012) showed that the information about the phase of the perturbation must also be used for recovering this behaviour.

Recent global simulations with radially extended domains revealed even more complexity. In particular, 3D simulations by Arzamasskiy, Zhu & Stone (2018) have shown that dT_ex/dR exhibits multiple sign reversals at |R − R_p| ∼ R_p (see their fig. 3). Evidence for this behaviour can also be found in the linear calculations of Miranda & Rafikov (2019a, see bottom rows of their figs 7 and 8) and simulations of Miranda & Rafikov (2019b, see their fig. 1), visible as the oscillations of the density wave angular momentum flux (AMF) F_J (its radial derivative is equal to dT_ex/dR in the absence of wave damping), which are especially pronounced in the outer disc. Sign reversals of dT_ex/dR in the outer disc can also be seen in the results of Dempsey et al. (2020, see the top row of their fig. 5) obtained for low perturber masses, although their simulations also included viscosity and formation of a gap around the planet. Our Fig. 1(b) also shows such dT_ex/dR features very clearly both as multiple, rather regularly spaced wiggles in the outer disc for R ≳ 2R_p and the lower amplitude, irregular features in the inner disc.

The goal of our present work is to shed light on the origin of the sign reversals and oscillations of dT_ex/dR. In particular, we demonstrate that these features are indeed real and not a numerical artefact. It is true that in all aforementioned cases, which fall in the linear regime of tidal coupling, the amplitude of these features is rather small, not exceeding several per cent of the main peak of the torque density (and decreasing with the distance from the planet). However, the magnitude of dT_ex/dR oscillations grows as the mass of the perturber increases and the disc–perturber coupling may become non-linear (see Dempsey et al. 2020). Thus, clarifying the origin of the features of dT_ex/dR behaviour in the linear regime also helps us understand the torque density behaviour in the non-linear regime, relevant for circumbinary discs around stellar and supermassive black hole binaries (Cimerman & Rafikov 2024).

Our work is organized as follows. After describing our set-up and methods in Section 2, we provide a heuristic explanation for the torque wiggles in Section 3. We then present a detailed theoretical model for the origin and properties of the torque wiggles in Section 4, which may be skipped at first reading, and explore their sensitivity to the disc parameters in Section 5 and to thermodynamic assumptions in Section 6. We further discuss our results in Section 7 and briefly summarize them in Section 8.

2 SET-UP AND METHODS

2.1 Model set-up

We consider a 2D model of a (initially axisymmetric) razor-thin disc that orbits a central star of mass M_⋆. We adopt polar (R, ϕ) coordinates centred on the central star. The disc has aspect ratio h ≡ H/R ≪ 1, where H = c_s/Ω is the vertical scale height, c_s is the sound speed, and Ω is the angular orbital frequency. We assume that the effective viscosity of the disc gas is low, implying that the flow is laminar (not turbulent), and do not include any explicit viscosity in our model. Nor do we include the disc self-gravity.

This disc is perturbed by the gravitational field of a coplanar planet with mass M_p, moving on a circular orbit with radius R_p, giving rise to non-axisymmetric perturbations to the basic state. Although our coordinate frame is centred on a (moving) central star, for simplicity we do not include the indirect potential when computing the response of the disc due to the presence of the planet. In this regard, we follow a large number of existing studies of disc–planet coupling that also neglected the indirect potential (Bate et al. 2003; D’Angelo & Lubow 2008, 2010; Duffell & MacFadyen 2012; Rafikov & Petrovich 2012; Miranda & Rafikov 2019a, b, 2020a; Fairbairn & Rafikov 2022). We note that the studies fully accounting for the indirect potential of the planet (Kley et al. 2012; Arzamasskiy et al. 2018) did not find noticeable differences from the calculations neglecting it. We also do not allow for the orbit of the planet to evolve with time.

In this work, we limit ourselves to planet masses that are well below the thermal mass (Goodman & Rafikov 2001)

where c_p = c_s(R_p) is the sound speed at R_p, Ω_p is the orbital angular frequency at R_p, H_p = H(R_p), and h_p = h(R_p). The assumption M_p ≲ M_th implies that the disc response to the planetary gravity is linear, i.e. δΣ/Σ ∼ M_p/M_th ≪ 1, where δΣ = Σ − Σ₀ is the planet-induced perturbation of the surface density Σ relative to the background surface density Σ₀(R).

Regarding the structure of the unperturbed (by the planet) disc, we assume that surface density and temperature obey a power-law ansatz:

such that the isothermal sound speed in the disc is

The disc is initially in radial centrifugal balance, taking into account the radial pressure gradient:

where |$\Omega _\mathrm{K}^2 = \sqrt{GM_\star /R^3}$| is the Keplerian orbital frequency, and P is the gas pressure (see below). The unperturbed velocity of gas is given by u_R,0 = 0, u_ϕ,0(R) = RΩ₀(R), with Ω₀(R) defined by equation (5) using Σ₀, P₀.

Regarding the equation of state (EoS), we explore several options. For most of our calculations, we adopt the adiabatic relation

where we fix the adiabatic exponent γ = 7/5 and K is the adiabatic constant. In our simulations, we explicitly solve the energy equation (see Appendix  A), thus properly accounting for the evolution of K (which would be conserved for each fluid element in a Lagrangian sense only in the absence of shocks or viscosity). It is important to keep in mind that the equation (3) does not imply the locally isothermal EoS: it is used, together with the equation (2), to set the initial radial profile of K in our adiabatic calculations. The adiabatic sound speed is c_s,adi = γ^1/2c_s,iso.

Additionally, in Section 6 we consider discs that allow for thermal relaxation of perturbations towards the unperturbed background following the so-called β-cooling prescription. We describe its detailed implementation in Appendix  B.

The treatment of thermodynamics has important implications on the global propagation of density waves excited by a planet. In particular, thermal physics directly affects the evolution of the AMF associated with these waves, defined as

with velocity perturbations u_R and δu_ϕ(R, ϕ) = u_ϕ(R, ϕ) − u_{ϕ, 0}(R). Miranda & Rafikov (2020a, b) have shown that for adiabatic discs obeying (6) the AMF carried by the density wave is conserved in the absence of wave damping (see also Section 4.2). On the contrary, discs with thermal relaxation, including the locally isothermal discs (Miranda & Rafikov 2019b), do not conserve the wave AMF, which may decay or get amplified even in the absence of non-linear damping. To distinguish these possibilities, we will call our models AMF-preserving when the disc thermodynamics is such (i.e. adiabatic) that the wave AMF is conserved in the absence of explicit damping, linear, or non-linear.

Linear calculations of the integrated one-sided Lindblad torque (GT80) find the characteristic magnitude of the AMF of the planet-driven density waves to be

which we will use as a reference value for F_J. To allow for meaningful comparison, radial torque densities are given in units F_J,0/R_p, as in Fig. 1.

2.2 Methods and torque calculation

We employ two methods to calculate the structure of a disc perturbed by a planet. First, we determine the disc structure in the linear approximation valid in the limit M_p/M_th ≪ 1. The method (involving a Fourier decomposition) for solving the linear problem has been developed in Miranda & Rafikov (2019a, 2020a) and we employ it here as well. Some details of its implementation and the parameters used can be found in Appendix A1.

Secondly, we also perform fully non-linear simulations (see Appendix A2 for details and numerical parameters) of the problem using Athena++³ (Stone et al. 2020), even though we still use M_p = 0.01M_th in these runs, which is well in the linear regime. This allows us to cross-check our linear calculations.

The perturbed pattern of the disc surface density Σ(R, ϕ) obtained by either of these two methods is then used to verify the results of our semi-analytical calculations of the torque wiggles in Section 4. More specifically, we compute the excitation torque density (torque per unit radius dR) by integrating R × dF/dR (with dF being the direct gravitational force exerted by a planet on to a disc element dS = RdRdϕ) over ϕ at fixed R as

which points in the z-direction. Note that this calculation accounts only for the direct planetary potential. In other words, we do not include the indirect potential in the torque calculation, to be consistent with all existing studies. Equation (9) uses a purely Newtonian potential for torque calculation, which is somewhat different from the softened potential (A1) used in our simulations. However, the two are essentially the same at the separations of interest for us.

We perform a parameter study of the problem by varying disc properties: surface density and temperature slopes p and q and aspect ratio at the planet radius h_p. In Table 1, we give an overview of the disc parameter sets that are explored in this work. We indicate the number of modes m_max obtained for the solution of the linear problem (see Appendix A1) and mark the parameter sets for which the athena++ simulations were performed. The adiabatic model with h_p = 0.1 and p = q = 1 highlighted in bold in Table 1 is the fiducial disc model that we will use to illustrate many of our results in Sections 3 and 4.

3 HEURISTIC EXPLANATION OF THE TORQUE WIGGLES

In this section, we provide a simple heuristic explanation for the origin of the torque wiggles seen in Fig. 1, and for the difference in their appearance in the outer and inner parts of the disc. This argument is then buttressed with quantitative calculations in Section 4.

In Fig. 2, we display a polar 2D map of the global surface density perturbation obtained via a linear calculation for the fiducial disc model: adiabatic with h_p = 0.1 (higher than h_p = 0.05 used in Fig. 1) and p = q = 1. One can easily see a one-armed spiral pattern that the planet excites inside and outside its orbit. To zeroth order, the shape of this spiral is given by the curve (R, ϕ_lin(R)), where⁴ (Ogilvie & Lubow 2002; Rafikov 2002a)

For a Keplerian disc (Ω → Ω_K) and sound speed profile in the form (4), one finds

This approximation works very well in the outer disc, as shown in Miranda & Rafikov (2019a). We see that over the entire domain R > R_p the spiral maintains a sharp crest (red), neighboured by a lower amplitude trough (blue). But in the inner disc this approximation starts failing for R ≲ 0.6R_p, as the spiral arm splits into multiple arms as it propagates. This evolution of a single-armed pattern into multiple spiral arms was found in simulations (e.g. Dong et al. 2015; Fung & Dong 2015) and explained by Bae & Zhu (2018) and Miranda & Rafikov (2019a). This difference⁵ has important implications for the shape of torque wiggles in different regions of the disc.

In Fig. 2, we also show the dashed grey line that passes through both the star and planet, i.e. ϕ = ϕ_p = const. As it turns out, this line has a special significance for explaining the nature of the torque wiggles. The crosses (not shown in the inner disc to avoid confusion) indicate the locations where the planet wake crosses this line, i.e. ϕ_peak = ϕ_p modulo 2|$\pi$|⁠. We will refer to this situation as wave- or wake-crossings. In the outer disc, where the approximation (10)–(14) works well, the radius R_n of nth such crossing (n = 0, 1, 2, …) is given by the condition φ(R_n) = 2|$\pi$|n, with R₀ = R_p.

To further support our statements, in Fig. 3 we show the same 2D map of the relative surface density perturbation δΣ/Σ₀ as in Fig. 2 but now in Cartesian geometry, see panels (a) and (b) for inner and outer disc, respectively. Also, in panels (c) and (d) we show the azimuthal profiles of δΣ in the outer and inner disc, respectively, at several fixed radii [shown in panels (a) and (b) using dashed lines of the same colour]. They are normalized by the characteristic wave amplitude δΣ_lin(R) expected in the linear regime (see Section 4.2) to highlight the evolution of the shape of the density wake. Finally, in panels (e) and (f) we show the radial profile of the excitation torque density dT_ex/dR on the same radial interval as in panels (a) and (b) to allow cross-matching of different features. We now examine this figure separately for the outer and inner discs.

3.1 Outer disc

In the outer disc (R > R_p), the planet-driven spiral retains a narrow one-armed shape as it wraps around the full azimuthal range of the disc multiple times. Moreover, panel (c) shows that the azimuthal profile of δΣ maintains its shape (upon normalization by δΣ_lin) to a good accuracy as the wake propagates, suggestive of a self-similar evolution. At the same time, the azimuthal location of the wake with respect to the ϕ = ϕ_p line (vertical dotted line) steadily changes, affecting the sign and value of dT_ex/dR.

At R = 2.17R_p (blue lines), the wake approaches its first passing of ϕ = ϕ_p and the torque density displays a sharp minimum. Panel (c) reveals that at this radius δΣ (blue curve) has a (positive) peak at ϕ > ϕ_p and (negative) trough at ϕ < ϕ_p, with the transition happening very close to ϕ = ϕ_p, when the wake is closest to the planet. In this optimal situation, both the wake overdensity (relative to unperturbed disc) at ϕ > ϕ_p and the wake underdensity at ϕ < ϕ_p pull the planet forward by their gravity, increasing its angular momentum. Newton’s third law then implies that the wake must be losing its angular momentum at this radius, resulting in strongly negative⁶ dT_ex/dR, which is what we see in panel (e).

As the wake crosses the line ϕ = ϕ_p (marked with a pink cross), the wake in panel (c) shifts to ϕ < ϕ_p, pulling back on the planet, which gives rise to a positive peak of dT_ex/dR at R = 2.28R_p (orange curve). This maximum is roughly only half of the magnitude of the preceding minimum, which is caused by the fact that the wake shape is not perfectly symmetric in ϕ (otherwise, we would expect similar magnitude of neighbouring extrema for tightly wound waves). Moving further out to R = 2.5R_p (green curves), the peak of the wake is close to ϕ = ϕ_p + π (i.e. behind the star relative to the planet), which lowers (still positive) dT_ex/dR.

Beyond that point dT_ex/dR becomes negative again and at R ≃ 2.8R_p (red lines) the wake profile in panel (c) almost coincides azimuthally with that at the previous minimum of dT_ex/dR (blue curve). As a result, dT_ex/dR reaches another (negative) minimum, although not as deep as the first one since this time the wake is further from the planet.

This pattern of alternating maxima and minima of dT_ex/dR (with steadily decaying amplitude) repeats in a very regular way. Comparing panels (a) and (e), it is clear that the radial periodicity of this pattern in dictated by the repeated wake crossings of the ϕ = ϕ_p line as the spiral propagates, with sharp minima occurring just before azimuthal alignments of the surface density peak with the planet. The overall shape of dT_ex/dR is only quasi-periodic: closer to the planet dT_ex/dR displays two local maxima (as wakes fully wraps around the star), while further out in the disc there is only one maximum. Nevertheless, the radial periodicity of this slowly evolving pattern of dT_ex/dR is very accurately set by the wake crossings, thanks to the essentially self-similar shape of δΣ in the outer disc.

3.2 Inner disc

We now turn to the inner disc (R < R_p) [see panels (b), (d), and (f)]. Close to the planet, at |R − R_p| ≃ (2 − 3)H_p = (0.2 − 0.3)R_p, the azimuthal profile of δΣ is similar to that in the outer disc, except that now it propagates in the opposite direction relative to the local mean flow. This can be seen in the grey curve in panel (d) drawn for R = 0.8R_p (sampling the main trough of dT_ex/dR), which is similar in shape (upon ϕ-reflection and rescaling) to the curves in panel (c).

However, closer to the star, and different from the wake in the outer disc, we notice an additional peak (red) in the surface density perturbation in panel (b) appearing next to the trough (blue) at R ≲ 0.6R_p. In the innermost regions, this secondary spiral arm is well-developed and there is a hint of a tertiary arm (peak) forming. This picture of multiple arm formation is consistent with Bae, Zhu & Hartmann (2017) and Miranda & Rafikov (2019a). These modifications are also reflected in the δΣ/δΣ_lin profiles for R = 0.265R_p, 0.18R_p, 0.118R_p in panel (d), which develop multiple peaks and troughs and certainly do not evolve in a self-similar fashion like in the outer disc (see also Figs 8 and 9).

The more complicated and steadily evolving (in R) azimuthal profiles of δΣ lead to the loss of coherence of dT_ex/dR at consecutive wake crossings in the inner disc, e.g. at R = 0.265R_p (pink) and R = 0.118R_p (magenta). As a result, one no longer sees the radial quasi-periodicity of the dT_ex/dR features, which are so obvious in the outer disc. Moreover, the amplitude of these features is also significantly reduced, partly because of the lower AMF carried by the wake in the inner disc (compare the magnitudes of the main peak and trough near the planet in Fig. 1) but also because of faster radial decay of the torque density in the inner disc (see Sections 4 and 4.5).

3.3 Toy model of the wiggles

The observations made in Sections 3.1 and 3.2 strongly suggest that the approximate self-similarity of the wake shape (along the azimuthal direction) as R varies is the critical factor for developing a regular pattern of the torque wiggles. This self-similarity exists in the outer but not in the inner disc, resulting in a rather irregular pattern of the wiggles for R < R_p.

To emphasize the role of self-similarity of δΣ even further, in Fig. 4 we show dT_ex/dR computed in the outer disc for a hypothetical ‘Gaussian’ wake in the form

where A is an (arbitrary) amplitude and we set w = h_p as the characteristic azimuthal width of the wake. This wake has a correct radial scaling of its amplitude to ensure the conservation of the wave AMF (see Section 4.2) and follows the curve ϕ = ϕ_lin(R) in polar coordinates.

One can see that this artificial wake gives rise to a regular pattern of torque wiggles (blue) with a clear radial periodicity and decaying amplitude. Because of the symmetric shape of the Gaussian wake, the consecutive peaks and troughs of the wiggles have similar amplitude as compared to the actual planet-driven dT_ex/dR (orange) also shown in that figure (see Section 3.1). But the pattern of radial variability is the same for both curves, highlighting that it is the ϕ_lin(R) behaviour (setting the radial frequency of wake crossings) that determines it.

4 THEORETICAL MODEL OF THE TORQUE WIGGLES

We now provide a more refined mathematical description of the torque wiggles. Our goal is to provide a semi-analytical description of these features in the outer disc, which would allow one to reproduce not only their radial (quasi-)periodicity but also their amplitude. We will demonstrate that this is possible in some cases of physical importance.

We start by introducing the planet-induced perturbation of the surface density δΣ(R, ϕ) = Σ(R, ϕ) − Σ₀(R). Since Σ₀ is axisymmetric and gives no net contribution to the torque, we can replace Σ(R, ϕ) with δΣ(R, ϕ) in equation (9). Always working in the frame co-rotating with the planet such that ϕ_p = 0, we can then write the magnitude of the torque density in the outer disc as

where we introduced α ≡ R_p/R < 1.

Integrating equation (16) by parts, we find

Using the series expansion

where

are the Laplace coefficients (e.g. Murray & Dermott 1999), one can re-write equation (16) as⁷

Up to a normalization, we recognize the remaining integral as the cosine Fourier coefficient of |$\partial \delta \Sigma /\partial \phi$|⁠:

such that we can write equation (20) as

where all information about the surface density perturbation is contained in the |$C_j^c(R)$|⁠.

In the inner disc, defining α ≡ R/R_p < 1, an analogous calculation gives

with an extra factor of α compared to the equation (22). This has important consequences for the amplitude of the torque wiggles in the inner disc (see Section 4.5).

So far we have not assumed any particular form of δΣ (and |$\partial _\phi \delta \Sigma$|⁠) such that our results (20) and (22) are fully general.

4.1 Torque wiggles for self-similar spiral arms

In the rest of the paper, unless mentioned otherwise, we will focus on exploring the properties of the torque wiggles in the outer disc. Inspired by our observations in Section 3.1 that, once formed, the outer spiral arm maintains its shape without significant distortion (see Fig. 3c), we provide a more in-depth characterization of the torque wiggles in the case of a fully self-similar wake in the outer disc. More specifically, we now assume the surface density perturbation due to the density wave to have the form

where function ψ describes the azimuthal structure of the wake, while the pre-factor |$\delta \tilde{\Sigma }(R)$| describes the evolution of the wave amplitude (equation (15) has such a form). In this self-similar picture, the wake travels with varying amplitude along the curves |$(R,\tilde{\phi }(R))$|⁠, and the only change of ψ with R is a translation in the ϕ-direction given by |$\tilde{\phi }(R)$|⁠.

With the ansatz (24), we have

where ψ′(x) = dψ(x)/dx.

Let us now introduce the complex Fourier coefficients of ψ′ as follows:

Writing |$\Psi _m = A_m e^{\mathrm{i}\theta _m}$|⁠, with A_m = |Ψ_m| and |$\theta _m = \arg \Psi _m$| being the amplitude and the phase of Ψ_m, respectively, we obtain

Substituting this expression into equation (20) gives, after straightforward manipulation,

Note that here |$\tilde{\phi }(R)$| is understood to lie in the interval (0, 2|$\pi$|⁠), i.e. is unique mod 2|$\pi$|⁠, since in equation (25) ϕ ∈ (0, 2|$\pi$|⁠).

This is the final expression for the torque density given the self-similar ansatz (24). Once the forms of δΣ(R), |$\tilde{\phi }(R)$|⁠, and Ψ_m (i.e. A_m and θ_m) are specified, the torque density can be computed by simple summation. Knowledge of Ψ_m is crucial for this calculation, and in Section 4.3, we discuss the properties of these coefficients for a certain class of disc models described next.

4.2 Self-similarity in AMF-preserving discs

Miranda & Rafikov (2019a) have shown that the self-similar ansatz (24) works well in the outer parts of the AMF-preserving (adiabatic) discs: the outer density wave maintains a single-armed, close to self-similar shape (see Fig. 4a of that paper that reveals only a weak evolution of the wake shape for R > R_p). This is also obvious from our Fig. 3(c) drawn for the fiducial disc model: azimuthal profiles of the wake at different radii look essentially identical (once one shifts each profile horizontally).

At the same time, the self-similar ansatz is clearly not applicable in the inner disc as a result of secondary (and tertiary) arm development there (see figs 4b, 5, 6 of Miranda & Rafikov 2019a). Our Fig. 3(d) also makes this clear since the wake profiles for different R have noticeably different shape and no horizontal shift would make them match.

This picture remains largely unchanged as one varies disc parameters, as we demonstrate in Section 5.1. Given all that, in the rest of this section we will focus on planet-driven density waves in the outer regions of AMF-preserving discs (with not too low h_p, see Section 5.1.2), for which the ansatz (24) works well.

By definition, the angular momentum of the freely propagating (i.e. not subject to further forcing) density waves is conserved in the AMF-preserving (e.g. adiabatic) disc models, i.e. |$\partial F_J/\partial R=0$|⁠. It has been shown (Rafikov 2002a) that conservation of F_J leads to the wave amplitude scaling as δΣ(R) ∝ δΣ_lin(R), where

In a power-law background disc (see equations 3 and 2), we find⁸

It is important to emphasize that the result (30) would not be valid in the non-AMF-preserving discs. For example, in discs with locally isothermal EoS F_J varies as |$F_J \propto c_\mathrm{s}^2(R)$| (Lin 2015; Miranda & Rafikov 2019a), so that the wave amplitude gets amplified (damped) relative to δΣ_lin in the inner (outer) disc. This can be easily accounted for by changing the power of c_p/c_s(R) in equation (30) from 3/2 to 1/2. Similarly, thermal relaxation (e.g. β-cooling) always damps F_J (Miranda & Rafikov 2020a, b), reducing δΣ(R) relative to δΣ_lin; in this case, no simple analytical correction to equation (30) is possible. We will discuss the impact of these alternative (to adiabatic) thermodynamic assumptions in Section 6.

The discussion above implies that the density waves in the outer parts of the AMF-preserving discs should be well approximated by the ansatz (24) with

with ϕ_lin(R) given by equations (10)–(14).

4.3 Properties of Ψ_m in AMF-preserving discs

Next, we explore the properties of the Fourier coefficients Ψ_m for the self-similar outer density waves in AMF-preserving discs. We employ the ansatz (24, 32) and consider an adiabatic disc with h_p = 0.1 and p = q = 1.

The Fourier coefficients of ψ can be obtained via the Fourier coefficients of δΣ, which are available to us from the solution of the linear problem. They are given by

Using equation (28), one can easily show that

In particular, the absolute magnitude of the Fourier coefficients of ψ and δΣ are related by |Ψ_m| = m|δΣ_m|/δΣ_lin.

In Fig. 5, we show the behaviour of the amplitude and phase of |$\delta \hat{\Sigma }_m$| for several values of R, which is helpful for understanding the properties Ψ_m later on. One can see in panel (a) the almost universal shape of |$\vert \delta \hat{\Sigma }_m (R)\vert$|⁠, weakly dependent on R and consistent with the self-similar description of δΣ. On the other hand, panel (b) looks like a scatter plot and shows that the phases of |$\delta \hat{\Sigma }_m (R)$| (computed in a frame with azimuthal axis aligned with ϕ_p) do not follow a clear pattern while varying rapidly. However, once we shift these phases by ϕ_lin(R) to account for the overall wrapping of the spiral wake (panel (c)), a clear pattern becomes obvious, showing only a slow evolution with R. This again supports the overall picture of a self-similar wake and shows that it closely follows the predicted path.

Next, we turn to examining the amplitude A_m = |Ψ_m| and phase θ_m = arg|$\, \Psi _m$| of the coefficients Ψ_m computed using |$\delta \hat{\Sigma }_m$| via equation (35). In Fig. 6, we show them as a function of mh_p, in the left and middle column, respectively, both in the outer disc (bottom) and, for completeness, in the inner disc (top).

Looking at the left column [panels (a) and (b)], we notice that A_m peak at around mh_p ≃ 1.45 (mh_p ≃ 1.2R_p) in the outer (inner) disc and do not show much variation with R. This is expected for free linear waves in flux-conserving (adiabatic and barotropic) disc models, since the radial scaling of δΣ with δΣ_lin(R) is already taken into account in ansatz (24) and (32). As we show in Section 6, the situation is very different for non-flux-conserving discs. Note that both (1) A_m peaking at mh_p ∼ 1 and (2) the decay of A_m at large mh_p ∼ 1 are directly related to the analogous behaviour of the modal contributions to the full planetary torque, which is ultimately caused by the torque cut-off phenomenon (GT80).

On the other hand, the phases θ_m (panels c and d), which include the phase shift by ϕ_lin(R) through equation (35), as in Fig. 5(c), show a substantial evolution with R. In the inner disc, the phases spread out with m as the distance from the planet increases. Around R = 0.1R_p, θ_m changes by 2π every Δ(mh_p) ≃ 2. On the other hand, in the outer disc, θ_m change a lot less with m and, in particular, show very weak variation with R at low m (similar to Fig. 5c).

Panel (f) shows the amplitudes A_m multiplied with Laplace coefficients |$b_{1/2}^{(m)}(\alpha)$| to illustrate the contribution of different azimuthal harmonics of δΣ to the torque density in the outer disc (see equation 29). Note that we show these over a smaller range of mh_p in order to zoom into the relevant region. For R closest to the planet, |$A_m b_{1/2}^{(m)}$| peak at around mh_p ≃ 0.4−0.5 and fall off rapidly at large m. As the distance from the planet increases, the location of the peak shifts towards lower m, with the overall amplitude of |$A_m b_{1/2}^{(m)}$| decaying. This is expected from the asymptotic behaviour of the Laplace coefficients (discussed in Section 4.5) and the weak dependence of A_m on R. In other words, the behaviour of the Laplace coefficients is the key factor determining which m provides the largest contribution to dT_ex/dR. Vertical dashed lines in all panels of Fig. 6 correspond to m_c defined as the critical m for which the product |$A_m b_{1/2}^{(m)}$| falls off below 1 per cent of its maximum value, such that only m < m_c need to be considered when computing dT_ex/dR (the values are m_c = 14, 9, 6, 4 for R = 1.5, 2, 3, 5R_p, respectively). These vertical dashed lines move to lower m as |R − R_p| increases, implying that closer to the planet, more modes contribute to dT_ex/dR, while for larger distances, only the lowest few m matter (see Section 4.5 for more details).

Similarly, in the inner disc we multiply A_m with |$\alpha b_{1/2}^{(m)}(\alpha)$| [panel (e), see equation 23]. Because of the extra factor of α (compared to the outer disc), the resultant curves precipitously diminish in magnitude as R decreases, explaining the weakness of wiggles in the inner disc. Dashed lines again mark m_c, for which |$A_m\alpha b_{1/2}^{(m)}(\alpha)$| fall to 1  per cent of the maximum value.

The behaviour of the phases θ_m in panel (c) makes it obvious that self-similarity is not a good approximation in the inner disc, consistent with the left columns of Figs (8) and (9). But the variation of θ_m at high m in panel (d) suggests that self-similarity may not work in the outer disc either. However, note that in panel (d) the variation of θ_m with R is not very dramatic for m < m_c (to the left of the corresponding vertical lines), i.e. for all azimuthal harmonics that matter for the calculation of dT_ex/dR. This allows us to consider δΣ as approximately self-similar in the outer disc, which will be used next.

4.4 Semi-analytic reconstruction of torque wiggles in AMF-preserving discs

Now, we demonstrate how the approximate self-similarity of the density perturbation in the outer disc allows us to reconstruct the main features of the torque wiggles behaviour. We compute dT_ex/dR using equation (29) with A_m, θ_m obtained under the assumption of a self-similar wake (see Section 4.3). Since Fig. 6(d) shows some evolution of θ_m with R, we pick a particular calibration radius R_cal = 3, at which we measure A_m, θ_m and assume these values to not depend on R (an approximation which should work reasonably well as we argued earlier) when computing the radial profile of dT_ex/dR. The resultant torque density is then compared with the dT_ex/dR obtained by solving the linear problem exactly.

In order to highlight the contribution of different individual modes to dT_ex/dR, we also use a truncated (at m = m_tr) version of the series expansion (29):

where here we have made explicit that under our self-similar approximation the values of (A_m, θ_m) are taken at a calibration radius R_cal.

In Fig. 7, we show the exact numerical solution for dT_ex/dR (blue solid) and the truncated self-similar series |$\mathrm{d}T^\mathrm{tr}_\mathrm{ex}/\mathrm{d}R$| (36) for several values of m_tr in panel (a); the deviation between them is shown in panel (b). One can see that close to the planet (R ≤ 2R_p) the self-similar approximation deviates significantly from the global linear solution for all m_tr. This is to be expected since θ_m still shows some dependence on R (compared to their values at R_cal) for m ≲ m_c close to the planet (see Fig. 6d). Further away, the first trough of dT_ex/dR is matched well both in shape and amplitude for all m_tr > 2. In the radial interval between this trough and the next one, we see rather small deviations for all m_tr. Beyond the second trough, the self-similar curves closely track the full solution for m_tr ≥ 2. For R ≥ 4R_p, the curve for m_tr = 2 merges with those for higher m_tr. The outermost torque oscillation period is matched reasonably well even for a single mode when m_tr = 1. Panel (b) reveals that for R ≥ 4R_p, the self-similar approximation shows a periodic deviation from the exact solution, with the amplitude that steadily decays as m_tr increases.

These results indicate that for sufficiently large distances from the planet, m_tr = 2−3 modes are sufficient to describe the behaviour of dT_ex/dR, if the Ψ_m are sampled at an intermediate distance R_cal within the radial region of interest. At the quantitative level, these results still depend on the value of R_cal, but this dependence should not curtail the good performance of the self-similar approximation described in Sections 4.1 and 4.2 in reproducing dT_ex/dR.

4.5 Understanding the key features of torque wiggles

As demonstrated in the previous subsection, equation (29) is able to successfully reproduce the correct behaviour of dT_ex/dR. Armed with this knowledge, we now use it to understand some key properties of the torque wiggles.

First, equation (29) allows us to predict reasonably well the shape of the outer envelope of the oscillating torque wiggles. This can be done by setting the cosine factors in equation (29) to unity, thus allowing the expression in the right-hand side to its maximum possible value,

This very simple constraint is plotted as the grey dot–dashed curve in Fig. 7(a) and it works remarkably well, with the troughs of peak wiggles almost touching this curve (see more on this below).

Secondly, equation (29) lets us analyse the torque density behaviour in the asymptotic limit R ≫ R_p (or α ≪ 1). Recall that in thin discs, H/R ≪ 1, the Fourier amplitudes of planet-driven δΣ_m peak at m_* ∼ R_p/H_p (i.e. m_*h_p ∼ 1) and exponentially decay for m ≳ m_* (GT80). The coefficients Ψ_m have an extra factor of m compared to δΣ_m, leading to A_m peaking at m slightly higher than m_*, which can be seen by comparing Figs 5 (see black and grey dotted lines) and 6. When computing dT_ex/dR via formula (29), A_m are multiplied by the Laplace coefficients |$b_{1/2}^{(m)}$|⁠, which asymptotically behave far from the planet as (Murray & Dermott 1999)

This suggests that as |R − R_p| increases, the main contributions to dT_ex/dR should be provided by the low-m components of δΣ. Eventually, far from the planet, only the m = 1 component would matter. In other words, we expect a transition from the interference of many modes closer to the planet to a single-mode dominated behaviour far away, just as we observed in Section 4.4.

Retaining only the first term in the sum in equation (29) and using the fact that |$b_{1/2}^{(1)} (\alpha) \rightarrow \alpha$| as α ≪ 1, we obtain the following asymptotic prediction for the outer envelope of dT_ex/dR:

where the dependence of δΣ_lin on R is given by the equation (30). This prediction is illustrated in Fig. 7 via the dotted magenta curves and it should be working well at large R ≫ R_p; at intermediate R it underestimates the envelope (37) of dT_ex/dR because of still significant contribution of other m harmonics.

Note that in the inner disc close to the centre, in the asymptotic limit α = R/R_p → 0, torque density is suppressed by an additional factor of α: dT_ex/dR ∼ δΣ₁(R/R_p)², where δΣ₁ is the m = 1 Fourier component of δΣ (see the equation 23). As a result, torque contributions of different harmonics fall off very rapidly as R decreases (see Fig. 6e). This, together with the decoherence of the planetary wake due to the formation of multiple spirals, acts to suppress the torque wiggles for R < R_p.

Thirdly, equation (29) also allows us to explain why in the outer disc dT_ex/dR features sharp (negative) peaks at radii, where the density wake crosses the line ϕ = 0 between the star and the planet. Indeed, note that θ_m in Fig. 6(d) exhibit rather small spread Δθ_m ≲ |$\pi$|/2 for all R and m ≤ m_c. This allows us to assume that θ_m is approximately the same for all m that contribute significantly to the sum in equation (29). In that case, these significant terms interfere constructively whenever mϕ_lin(R) ≲ 1 for all m ≲ m_c (this also explains why the troughs of dT_ex/dR overshoot the envelope (39) at intermediate R). As a result, whenever ϕ_lin(R) → 0 mod 2|$\pi$|⁠, i.e. the wake (centred on ϕ_lin(R)) crosses the line ϕ = 0 from the star to the planet, constructive interference of many terms in equation (29) leads to a sharp increase of the amplitude of dT_ex/dR, producing the characteristic appearance of the torque wiggles. In fact, the interference at these locations is so effective that the trough wiggles (negative values) almost reach the maximum possible absolute value of |dT_ex/dR| given by equation (37).

Fourthly, we can also use equation (29) to estimate the radial width of these sharp features of dT_ex/dR (although the outcome of this exercise may be of limited utility). Indeed, constructive interference of harmonics with m ≲ m_c leading to torque wiggles requires |$\delta \phi _\mathrm{lin}\lesssim m_\mathrm{c}^{-1}$|⁠, constraining the radial width of the sharp troughs of dT_ex/dR as

(see equation 11). Because of the rapid evolution of m_c with R (see Fig. 6), it may be somewhat challenging to extract a more explicit dependence of δR on R and disc parameters. Nevertheless, it is clear that δR should be smaller in thinner discs with lower c_s and h_p.

5 RESULTS: VARIATION OF DISC PARAMETERS IN FLUX-PRESERVING DISCS

Next, we examine how the characteristics of the torque wiggles in the AMF-preserving discs are affected by variation of the disc parameters. Recall that δΣ_lin depends on p and q – the slopes of the background surface density and temperature (which set the initial profile of entropy, see equations 2, 3, and 30), while ϕ_lin also depends on h_p (see equations 4  and 10–11). We will first explore in Section 5.1 how the self-similarity of the wake in the AMF-preserving discs is affected by the variation of p, q, and h_p, and then move on to study the effect of these parameters on the torque wiggles in Section 5.2.

5.1 Self-similarity of the wake structure

We illustrate the variation of the wake profile with disc parameters by examining the azimuthal dependence of the dimensionless profile shape ψ(x) defined by the equation (24) as a function of the (naturally scaled) azimuthal variable |$(\phi -\phi _\mathrm{lin})h_\mathrm{p}^{-1}$|⁠. This allows us to eliminate the obvious dependencies of δΣ on the disc parameters via δΣ_lin and ϕ_lin and to focus on more subtle effects. These ψ-profiles are obtained via the linear calculation for several adiabatic disc models described below.

5.1.1 Variation of the slopes p and q

In Fig. 8, we show ψ in the inner (left) and outer (right) disc with h_p = 0.1 for different values of p and q (colours) at several values of R (columns). These plots reveal a pattern which is reminiscent of the discussion of the fiducial disc model in Section 4.2.

Namely, in the outer disc the wake profile remains roughly independent of R suggestive of its approximate self-similarity. The variation of p and q introduces some deviations of ψ(x) from its shape in the fiducial disc model (blue curve), but these differences are not significant. At the same time, in the inner disc we again observe the formation of additional maxima of ψ which appear due to wake splitting into multiple spiral arms (Bae et al. 2017; Miranda & Rafikov 2019a). The structure of the inner wake does change noticeably as the disc parameters are varied, with the differences getting amplified towards the disc centre, e.g. see panel (c). It is also clear that the shape of the wake in the inner disc is affected mainly by the temperature slope q, which directly affects the emergence of additional spiral arms.

To summarize, the variation of p and q does not greatly modify the overall picture of wake propagation (for h_p = 0.1): the outer wake maintains a self-similar shape with good accuracy, while the inner wake does not follow this universal pattern.

5.1.2 Variation of h_p

We now investigate the effect of varying h_p on the wake shape, keeping p = q = 1 fixed. In Fig. 9, we show ψ for several values of h_p, similar to Fig. 8. While the rescaling with δΣ_lin and stretching in ϕ by a factor |$\propto h_\mathrm{p}^{-1}$| allow us to compare the wake shape on similar scales in this plot, we see clear deviations between the ψ profiles for different h_p, even close to the planet (panels a  and d).

As the distance from the planet increases, the wake suffers from the dispersive effects, which are especially pronounced for small h_p, consistent with the predictions of Miranda & Rafikov (2019a). This process is particularly noticeable in the inner disc (left), where the tertiary spiral becomes clearly visible at R = 0.1R_p for h_p = 0.025, while the h_p = 0.1 disc shows only two clear spiral arms.

What is more remarkable, for very low h_p a secondary arm starts to develop even in the outer disc. This can be noticed by comparing the h_p = 0.025 curves at R = 1.1R_p and 3R_p in the right columns of Fig. 9, revealing the development of a weak secondary arm far from the planet consistent with the analysis of Miranda & Rafikov (2019a). This indicates a modest violation of self-similarity within the radial region of interest occurring even in the outer disc for small enough h_p.

To summarize, h_p is a very important parameter affecting the density wave propagation more significantly than the changes in p and q.

5.2 Effect of disc parameters on torque wiggles

Next, we investigate the effect of changing disc parameters on the characteristics of torque wiggles. In Fig. 10, we show dT_ex/dR computed both using full linear theory and using direct athena++ simulations for a variety of disc parameters; a parameter that changes compared to the fiducial model (panel a) is labelled in each panel. Vertical magenta (dashed) lines correspond to the radial locations where the peak of the wake crosses the line ϕ = 0 from the star to the planet, while the green (dotted) lines correspond to the analytical prediction ϕ_lin = 2|$\pi$|n, n = 1, 2,...; one can see that the two agree very well. Note that panels (a)–(d) all have h_p = 0.1.

Comparing panels (a) and (b), we see that changing the surface density slope p does not affect the radial periodicity of dT_ex/dR, i.e. the positions of prominent troughs (see also dashed vertical lines). This is expected since ϕ_lin does not depend on Σ₀. However, the amplitude of the torque wiggles changes and in a non-uniform fashion: p = 0 disc shows lower (higher) dT_ex/dR in the inner (outer) disc compared to p = 1. This can be understood from the dependence of the excitation torque density on Σ₀(R): |$\mathrm{d}T_\mathrm{ex}/\mathrm{d}R\propto \delta \Sigma \propto \delta \Sigma _\mathrm{lin}\propto \Sigma _0^{1/2}(R)$| (see equation 16). Thus, given their approximately identical ψ(x) for discs with different p (see Figs 8d and f), we expect that curves of |$\Sigma _0^{-1/2}\mathrm{d}T_\mathrm{ex}/\mathrm{d}R$| would have a universal shape. This is indeed the case, as illustrated by Fig. 11 for p = 1 (solid) and p = 0 (dashed): with this normalization the two curves track each other closely over the entire radial range.

Panels (c) and (d) illustrate torque wiggles in discs with different initial entropy profiles by considering a reduced temperature slope q = 0.5 and a constant temperature (q = 0). It is clear that the variation of q leads to changes of not only the amplitude of the wiggles but also their radial periodicity – the distance between the consecutive troughs increases (decreases) in the outer (inner) disc as q goes down. This is explained via the behaviour of ϕ_lin, which determines the global shape of the wake. Regarding the amplitude, equation (30) predicts dT_ex/dR ∝ δΣ_lin ∝ (R/R_p)^3q/4 in AMF-preserving discs, such that lower q should increase (decrease) dT_ex/dR in the inner (outer) disc. This is consistent with the trends seen in panels (c) and (d).

Finally, in panels (e) and (f) we consider disc models with lower h_p, while keeping p = q = 1 as in panel (a). Comparing the locations of wake crossings with the radial structure of dT_ex/dR it is evident that when h_p is halved, the number of wake-crossings doubles, as expected from |$\phi _\mathrm{lin}\propto h_\mathrm{p}^{-1}$|⁠. Closer examination of panel (e) reveals that in the outer disc with h_p = 0.05, wake-crossings are again aligned with minima of dT_ex/dR for the first few radial periods. However, these minima are preceded by a maximum of similar magnitude, contrary to the fiducial case, where the troughs clearly dominate. This indicates that the variation of the shape of ψ as h_p is decreased (see Fig. 9) has a direct effect on dT_ex/dR. For even lower h_p = 0.025 (panel f), the wake-crossings in the outer disc align with the dominant peak or trough of dT_ex/dR even less; one also notices variation of the wiggle pattern, even between neighbouring (e.g. first and second) wake-crossings. Again, this is clearly caused by the erosion of the wake self-similarity for low h_p (see Section 5.1.2).

The amplitude of variations in dT_ex/dR (normalized by F_{J, 0} at fixed R) in panels (e) and (f) also clearly decreases with h_p, compatible with the expectation |$\mathrm{d}T_\mathrm{ex}/\mathrm{d}R\times F_{J,0}^{-1} \propto \delta \Sigma _\mathrm{lin}h_\mathrm{p}^{3} \propto h_\mathrm{p}^{3/2}$| (see equations 8 and 30).

6 TORQUE WIGGLES IN NON-AMF-CONSERVING DISCS

In this section, we study the effect of relaxing the wave AMF conservation on the existence and properties of torque wiggles. As an example of a non-AMF-conserving disc, we will consider a disc with thermal relaxation in the form of β-cooling (see Section 2.1 and Appendix  B). The propagation of planet-driven waves in such discs has been previously explored by Miranda & Rafikov (2020a). We illustrate our results using several values of β ranging from 10⁻³ to 10², which allows us to explore the transition from an isothermal to adiabatic thermodynamics. We also consider the pure locally isothermal limit β → 0 (Miranda & Rafikov 2019b).

We show dT_ex/dR results for different values of the dimensionless cooling time β (with the fiducial set of disc parameters) in Fig. 12, again using both the linear calculation and simulations (which are in excellent agreement). A first glance reveals that the variation of β clearly has an effect on the amplitude and radial periodicity of the torque wiggles, both in the inner and outer disc.

As β decreases, the radial pattern of repeating troughs shifts closer to the planet. This can be understood as follows: for β ≫ 1 the thermodynamic response is approximately adiabatic and c_s ≈ c_s,adi = γ^1/2c_s,iso – i.e. the adiabatic sound speed, which because of the factor γ > 1 is higher than the isothermal sound speed c_s,iso, to which c_s reduces in the opposite limit β ≪ 1. As a result, equation (10) predicts slower evolution of ϕ_lin with R, i.e. a less tightly wound density wave for β ≳ 1 as compared to β ≲ 1 case (see also Miranda & Rafikov 2020a). This point is illustrated using vertical blue (dotted) lines corresponding to ϕ_lin = 2|$\pi$|n, n = 1, 2,..., which are computed using c_s ≈ c_s,adi in (10) for β ≥ 1 and using c_s ≈ c_s,iso for β < 1; these lines are well aligned with the troughs of dT_ex/dR.

The amplitudes of torque wiggles are quite similar in the adiabatic (β = 10²) and locally isothermal limits, compare panels (a) and (e). At the same time, the wiggles are strongly suppressed for intermediate values of β, especially in β = 0.1 disc (panel c) and at large R. These trends can be easily understood by examining the behaviour of the amplitudes A_m and phases θ_m of the Fourier coefficients Ψ_m (see Section 4.3) in the outer disc, which are shown in Fig. 13. In doing this, it is helpful to keep track of the vertical coloured lines at m_c – the maximum m significantly contributing to dT_ex/dR – at every R (see Fig. 6 and Section 4.3).

For β = 10² (panel a), close to adiabatic limit, A_m and θ_m behave as they do in Fig. 6 – the former are independent of R, while the latter slowly evolve with R. As a result, the wiggle pattern is essentially indistinguishable from that of Fig. 10(a). For β = 1, we see A_m being uniformly (in m) suppressed as R increases, which is caused by the decay of the wave AMF due to thermal relaxation for intermediate values of β, as described in Miranda & Rafikov (2020a). In particular, note the obvious reduction of A_m to the left of the corresponding m = m_c lines. This reduces the amplitude of the wiggles seen in Fig. 12(b).

Wave decay is most dramatic for β = 0.1 (Fig. 13c), leading to the precipitous drop of A_m for m < m_c and the corresponding rapid decay of the torque wiggle amplitude with R in Fig. 12(c). Moreover, for this value of β the phases θ_m also show rather erratic behaviour, indicative of the loss of constructive interference of different Fourier harmonics, again consistent with Miranda & Rafikov (2020a) findings.

The recovery of the wiggle amplitude for β = 10⁻³ (see Fig. 12d) is due to the disc starting to transition to the locally isothermal regime. Fig. 13(d) shows that A_m decays with R slower than in panel (c), in particular, the most significant A_m(m) (for m < m_c) for every R are not that different even from panel (a). This agreement is improved even further in the locally isothermal limit shown in Fig. 13(e) (although here the suppression of A_m with R is more uniform⁹ in m, more akin to panel b), explaining the qualitative similarity of dT_ex/dR in panels (a) and (e) of Fig. 12 (modulo the differences in the radial periodicity).

To summarize, the evolution of the torque wiggles with the cooling time in discs with β-cooling is entirely consistent with the results of Miranda & Rafikov (2020a) on the damping of planet-driven density waves in such discs.

7 DISCUSSION

In this work, we provided a detailed exploration of the properties and origin of the torque wiggles – a recently recognized feature of the global disc–planet interaction which has been seen in only a handful of studies so far (Arzamasskiy et al. 2018; Miranda & Rafikov 2019a, b; Dempsey et al. 2020). Together with the ‘negative torque density phenomenon’ discovered in the local limit (Rafikov & Petrovich 2012), this feature represents an interesting departure from the simple dT_ex/dR ∝ |R − R_p|⁻⁴ picture that has been prevalent¹⁰ in the field for decades (e.g. Armitage & Natarajan 2002; Chang et al. 2010; Zagaria, Rosotti & Lodato 2021), and paves a way to understanding some key aspects of the tidal disc–perturber coupling in other settings.

Our study clearly shows that torque wiggles owe their existence to the global structure of the planet-driven density wave and its direct coupling to the planetary gravity. Periodicity of the wiggles arises because of the wave wrapping due to the differential rotation, making possible azimuthal alignments of the wake with the planet, which result in sharp features of dT_ex/dR (Section 4.5). This situation would not be possible in the (infinitely extended) shearing sheet geometry, meaning that torque wiggles do not arise in the local limit.¹¹ This is consistent with the calculations of Rafikov & Petrovich (2012) who showed that the negative torque density phenomenon (i.e. no additional dT_ex/dR reversals) is the only possibility far from the planet in the shearing sheet. It is also important to emphasize that the torque wiggles are not related to the fluid behaviour in the vicinity of the low-order Lindblad resonances (GT80) which are not present far from the planet. This issue is explored in more detail in Cimerman & Rafikov (2024).

In a fully global setting, we have shown the torque wiggles to be a very robust and ubiquitous phenomenon. They are present for different disc parameters (Section 5), different assumptions about the disc thermodynamics (Section 6), with and without accounting for the non-linear effects (more on this in Section 7.1). Torque wiggles were originally seen in 3D simulations of Arzamasskiy et al. (2018), but we clearly see them also in the 2D setting. Also, the work of Arzamasskiy et al. (2018) accounts for the indirect potential¹² when computing the response of the disc to the planetary perturbation, whereas our study (as well as Miranda & Rafikov 2019a, b and many other works) neglects the indirect potential. This implies that the dimensionality of the problem and indirect potential are not the key factors for the origin of the torque wiggles (the latter will be explicitly demonstrated in Rafikov & Cimerman, in preparation).

Self-similarity of the wake (even if approximate) is another important ingredient for the appearance of a regular pattern of torque wiggles. Rapid dispersive decoherence of the wake in the inner disc resulting in the wake splitting into multiple spiral arms is one of the reason why the wiggles are aperiodic and largely suppressed for R < R_p (another one is the faster radial decay of dT_ex/dR there, see Section 4.5). In the outer disc this decoherence is much weaker as shown by Miranda & Rafikov (2019a) and the wake maintains a self-similar form with good accuracy (except for very low h_p). We note that in the presence of a well-developed turbulence, the wake coherence may be destroyed even in the outer disc (Zhu, Stone & Rafikov 2013), resulting in the suppression of the torque wiggles (although they may still be evident after time-averaging).

Shear viscosity in a (weakly-turbulent) hydrodynamic disc is unlikely to challenge the overall picture of the torque wiggles. First, Miranda & Rafikov (2020b) have shown that viscous dissipation is ineffective at damping the density waves as compared to other processes – radiative and non-linear damping. Thus, the wake is likely to maintain its global structure unless the viscosity is very large, α ≳ 0.1. Second, viscous dissipation of the wake would mainly damp the high-m Fourier modes contributing to the wake. However, we have seen in Section 4.4 that it is mainly the low-m modes that dominate in the genesis of the torque wiggles, and increasingly so at large R. Our own simulations (which we do not show here) demonstrate that even for effective viscosity as high as α ≤ 10⁻² the appearance of torque wiggles does not change much, in agreement with the viscous simulations of Dempsey et al. (2020) which also exhibit torque oscillations.

7.1 Effects of wave non-linearity

Non-linear effects are important for propagation of the planet-driven density wave: the wake steepens into a shock, which leads to wave dissipation. This process is inevitable even for low-mass planets; see Fig. 14(a) in which we illustrate the non-linear evolution of the azimuthal profile of the wave in a fiducial disc model for M_p = 0.01M_th. While the profiles computed using linear theory (dashed curves) show almost no change with R, δΣ derived from simulations fully accounting for the wave non-linearity evolve significantly: as R increases, the wave stretches azimuthally and its amplitude goes down. This is reflected also in the behaviour of the AMF carried by the wave as illustrated in Fig. 15, where in panel (b) the orange curve shows a clear decay¹³ of F_J (extracted from our simulations) in the outer disc after the wave has shocked. One would expect this wave damping to have a direct impact on the appearance of torque wiggles.

However, comparison of dT_ex/dR profiles computed using linear theory and simulations in Figs 10 and 12 shows that this is not the case and, at least for the low M_p used in simulations, non-linear damping of the wave has rather little impact on the torque wiggles. This unexpected outcome can be understood based on the fact that far from the planet dT_ex/dR behaviour is determined only by a handful of low-m Fourier harmonics of δΣ (see Section 4.4 and Fig. 7). Moreover, Fig. 14(b) shows that while the amplitudes (represented by A_m) of high-m harmonics rapidly die out, the amplitudes of the low m = 1, 2 harmonics tend to grow with R, which happens because of the non-linear azimuthal stretching of the profile, see panel (a). These two effects conspire to keep the amplitudes of the torque wiggles essentially the same for M_p = 0.01M_th, despite the non-linear damping of the wake.

We expect that for higher M_p the non-linear wave damping will have a more dramatic effect on the torque wiggles, reducing their amplitude; this is supported by our preliminary runs in high M_p regime which we do not show here. Nevertheless, fig. 1 of Miranda & Rafikov (2019b) still shows some torque wiggles up to M_p = M_th (see also Dempsey et al. 2020). When wave non-linearity is strong, the shape of the torque wiggles may become sensitive to the slope of the surface density p (unlike in the linear case, see Fig. 11), since the rate of non-linear evolution of the wave depends on it (Rafikov 2002a; Cimerman & Rafikov 2021).

7.2 Effect of torque wiggles on the total torque

Fig. 15 also illustrates the impact of the torque wiggles on the full (integrated) one-sided torque exerted by the planet on the disc, by showing the behaviour of the cumulative torque (dashed curves in panel b)

One can see that close to the planet T_ex(R) follows F_J(R) quite closely, except for the small vertical offset related to the torque in the corotation region. The two start to diverge after the wave shocks and starts transferring its angular momentum to the disc fluid, with AMF decaying below T_ex(R) in the nonlinear case. Despite that, both T_ex(R) and F_J(R) clearly show the low-amplitude oscillations caused by the torque wiggles. It is interesting that the torque wiggles have an impact (albeit small) on the wave AMF even far from the planet, where the wave is usually considered to be freely propagating and no longer affected by the planetary potential.

These oscillations are more easily seen in panel (c), where we plot ΔT = T_ex(R) − T_ex(R_r) – the difference between the cumulative torque at R and the cumulative torque at some reference R_r. We pick somewhat arbitrarily R _r,i = 0.7R_p and R_{r, o} = 1.3R_p in the inner and outer disc, respectively (dotted black vertical lines), just outside the main peaks of dT_ex/dR. In the inner disc, the total torque hardly changes and ΔT remains close to zero for R < R_{r, i}. In the outer disc, however, we clearly see the variations due to the integrated effect of the torque wiggles. The greatest deviation from the maximum T_ex is ΔT/T ≃ −0.03 and its location coincides with the first wake crossing of ϕ_p around R ≃ 2.2R_p. These oscillations eventually die out as R → ∞; however, their integrated effect is non-zero and ΔT(∞)/T ≃ −0.015 (see panels b and c). Thus, the integrated effect of the torque wiggles in the outer disc is to slightly lower T_ex from the maximum value near R ≈ 1.3R_p.

In our models, the highest relative variation of T_ex in the outer disc due to torque wiggles is ΔT/T_ex ≃ 0.07 (after the first wake-crossing) for p = 0, q = 0, h_p = 0.1 disc. Clearly, this effect is small, as is the amplitude of torque wiggles in general, as the main peak of dT_ex/dR within a couple of scale heights from the planet provides the dominant contribution to T_ex (see Fig. 1). Moreover, T_ex is the integral of the radially oscillating dT_ex/dR (see equation 29), which additionally lowers ΔT. However, things would be different if the planet were more massive¹⁴ and opened a gap around its orbit. The reduced surface density would suppress the near-planet peak of dT_ex/dR, accentuating the effects of the torque wiggles (see fig. 5 of Dempsey et al. 2020) for an illustration of this trend. A study of the torque behaviour in this regime has been presented in Cimerman & Rafikov (2024).

8 SUMMARY

In this work, we explored the physics of torque wiggles – the conspicuous low-amplitude features in the radial profile of the torque density (due to the direct force from the planet) at large distances (|R − R_p| ≳ R_p) from a planet that launches a density wave in a protoplanetary disc. These features are characterized by a remarkable pattern of radial variability (Section 3) and have been previously seen in a number of simulations. We explored their properties using linear theory and direct hydrodynamic simulations (finding good agreement between them in appropriate limit of low planet masses), and developed a theory explaining their origin. We summarize our main findings as follows:

• The wiggles arise due to the global gravitational coupling of the planetary potential to the planet-driven density wave propagating through the disc. They appear predominantly in the outer disc, where their amplitude is substantial and their radial periodicity is obvious (Section 3.1). They are suppressed and are irregular in the inner disc (Section 3.2). The following statements are for the outer disc.

• The wiggles appear for all disc parameters that we explored in this work (Section 5) and for various assumptions about the disc thermodynamics, although they are strongly suppressed in discs with (the dimensionless cooling time) β ∼ 1 (Section 6). They appear in both 2D and 3D discs (Arzamasskiy et al. 2018).

• While we mainly explore the wiggle properties using linear theory, they are also remarkably insensitive to the non-linear effects (Section 7.1). Their amplitude scales as |$M_\mathrm{p}^2$|⁠, just as that of the full planetary torque.

• We developed analytical theory (Section 4) showing that quasi-periodic wiggles arise because of the global nature of the planet-driven density wave (its multiple wrappings in the disc) and the approximately self-similar behaviour (Section 4.2) of the wave in the outer disc (which starts getting violated as h_p → 0).

• This theory demonstrates that the behaviour of dT_ex/dR at |R − R_p| ≳ R_p is dominated with good accuracy by a small number of low-m Fourier components of the perturbation δΣ (Section 4.4).

• This theory also explains the key features of dT_ex/dR behaviour (Section 4.5). In particular, it accounts for the sharp features in dT_ex/dR that appear at the radii where the planetary wake crosses the star–planet line as a result of constructive interference of the most significant low-m modes of the density perturbation (in the inner disc violation of the constructive interference related to the formation of multiple spiral arms suppresses the wiggles).

• While the amplitude of the wiggles is low in the linear regime and they affect the total (one-sided) torque imparted by the planet on the disc only at the level of several per cent (Section 7.2), their significance will grow for higher planet masses when the gap opening becomes important.

Cimerman & Rafikov (2024) have also explored torque behaviour in the non-linear regime of circumbinary discs, when the perturber-to-star mass ratio q ∼ 1. They demonstrated that the torque wiggles become much more important in this regime.

ACKNOWLEDGEMENTS

We are grateful to Lev Arzamasskiy for useful discussions and to an anonymous referee for useful suggestions. NPC is funded by an Isaac Newton Studentship and a Science and Technology Facilities Council (STFC) studentship. RRR acknowledges financial support through the Ambrose Monell Foundation, and STFC grant ST/T00049X/1. A large part of the long term simulations were performed on the HPC cluster FAWCETT at DAMTP, University of Cambridge. Part of this work was performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (www.csd3.cam.ac.uk), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP/P020259/1), and DiRAC funding from the Science and Technology Facilities Council (www.dirac.ac.uk).

Software:numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), ipython (Perez & Granger 2007), matplotlib (Hunter 2007), athena++ (Stone et al. 2020).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

References

APPENDIX A: NUMERICAL METHODS

A1 Set-up of the linear solver

To obtain linear solution for the planet-induced perturbation we employ the exact same machinery as in Miranda & Rafikov (2019a, 2020a). The only difference is that now we use a more accurate approximation for the planetary potential with a softening prescription that converges to a Newtonian potential at higher (fourth) order close to the planet (Dong et al. 2011):

where d = |r − r_p| is the distance from the planet and r_s = ϵH(R = R_p) is the smoothing length for which we adopt ϵ = 0.6. This choice is consistent with our previous work (Cimerman & Rafikov 2021, see Appendix A therein).

The radial mode solutions are obtained on a grid ranging from R_in = 0.05R_p to R_out = 5R_p consisting of |$N^\mathrm{lin}_R = 10^5$| logarithmically spaced cells for all values of h_p > 0.025. For this lowest value, we consider a smaller radial domain ranging from R_in = 0.1R_p to R_out = 3R_p and keep |$N^\mathrm{lin}_R$| fixed to increase the radial resolution, which is needed as the wave pattern becomes very tightly wound.

The number of modes m_max that are typically needed for a converged solution depends on the disc scale-height h_p. In particular, we used m_max = 100, 220, 320 for h_p = 0.1, 0.05, 0.025 respectively. These values were chosen on a case by case basis until convergence was found.¹⁵

A2 Set-up of non-linear simulations

Our set-up is similar to the one used in Cimerman & Rafikov (2021); however, here we additionally solve an equation for the total energy of the fluid, to allow for initial temperature gradients (setting the initial entropy) in the disc. The total set of equations solved is

where E = ϵ + ρ|u|²/2 is the total energy density, ϵ = eρ is the internal energy density, P = (γ − 1)ϵ is the gas pressure with the ratio of specific heats γ, I is the identity tensor, Φ = Φ_⋆ + Φ_p is the total gravitational potential, where Φ_⋆ = −GM_⋆/R, and Φ_p (given by equation A1), are due to the central star and the planet, respectively (i.e. no indirect potential). The source term Λ allows for thermal relaxation to a fixed background and is described in Appendix  B.

The simulation domain ranges over 0.1 ≤ R/R_p ≤ 5.0 in radius and is 0 ≤ ϕ ≤ 2π in azimuth. We use N_R × N_ϕ = 2240 × 3600 cells in the fiducial set-up (when h_p = 0.1) that are logarithmically spaced in radius and uniformly spaced in azimuth. For cases with h_p = 0.05, this resolution is doubled. We do not perform an athena++ simulation for h_p = 0.025 due to the excessive resolution requirements. We have checked that our main findings are not influenced by resolution by performing a double resolution simulation for the fiducial model.

Wave-damping zones (de Val-Borro et al. 2006) are implemented close to the radial boundaries in the zones 0.1 ≤ R/R_p ≤ 0.15 and 4.5 ≤ R/R_p ≤ 5 in order to avoid reflections (see Cimerman & Rafikov 2021).

APPENDIX B: β-cooling prescription

Following Miranda & Rafikov (2020a), we use a prescription where thermal relaxation towards the unperturbed specific internal energy density |$e = \epsilon /\Sigma = c_\mathrm{s,adi}^2/[\gamma (\gamma -1)]$| occurs on a time-scale t_rel(R) = β/Ω(R) as

where δe is the thermal energy perturbation. Here, β is a spatially constant parameter (dimensionless cooling time), which we vary between models. We implement this term in athena++ by retrieving the internal energy density from the total energy density, which is then changed each time-step using the exact solution to equation (B1):

and added back to the kinetic term to update the total energy density. The implementation of this additional effect is verified by comparison with the linear solver.

The prescription (B1) includes two important limits: for β → ∞ the equation of state becomes adiabatic (i.e. individual fluid elements conserve their entropy in the absence of shocks), e.g. as in our fiducial model; for β → 0 one recovers the locally isothermal limit such that c_s = c_s,iso(R). Depending on the thermodynamic model, the effective sound speed varies, as a function of β (Miranda & Rafikov 2020a).

In practice, we set the cooling parameter β = 10² to achieve a situation close to an adiabatic limit. We tested and confirmed that these runs give virtually identical results to even greater β = 10³, indicating convergence towards the adiabatic limit, in line with the findings of Miranda & Rafikov (2020a).