Torque wiggles — a robust feature of the global disc-planet interaction

Gravitational coupling between planets and protoplanetary discs is responsible for many important phenomena such as planet migration and gap formation. The key quantitative characteristics 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 should increase in the non-linear regime, when a gap (or a cavity) forms around the perturber’s orbit.


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 1 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 (Rafikov & Petrovich 2012;Petrovich & Rafikov 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 ★ E-mail: rrr@damtp.cam.ac.uk (RRR) 1 In the rest of the paper we will usually call the perturber a 'planet', regardless of its exact nature.
is the so-called excitation torque density d ex /d, which is defined in polar cylindrical coordinates r = (, , ) as the -component of the torque exerted by the (Newtonian) potential of the planet per unit radial distance in the disc.One can also interpret d ex /d 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 d ex /d 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 d ex /d 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 d ex /d is of paramount importance for obtaining a complete picture of disc-planet interaction.
Some of the key features of d ex /d behaviour in a disc which is roughly uniform near the planet (i.e. with no gap formed) have been predicted already in GT80.In particular, they have shown using linear theory that d ex /d should rapidly decrease as | −  p | ≲  p (were  p is the semi-major of the planet on circular orbit and  p is the disc scaleheight at  p ), an effect that is known as 'torque cutoff'.The overall shape of d ex /d near the planet is set by the tidal discplanet coupling at numerous Lindblad resonances, with the largest contribution coming from resonances located at | −  p | ∼  p .This is indeed what one sees in Fig. 1, in which we show a typical run of  8)) using linear calculation (Section 2.2) for disc parameters ℎ p = 0.05 and  =  = 1 (se Section 2.1).Top and bottom panels show the same data on different vertical scales.The arrows mark the main peaks ('main') of d ex /d and indicate the first sign change of d ex /d outside the main peak, also known as the negative torque density phenomenon ('NT', see Dong et al. 2011a;Rafikov & Petrovich 2012).Far away from the planet the label 'torque wiggles' marks the region where d ex /d exhibits multiple sign changes (Arzamasskiy et al. 2018;Miranda & Rafikov 2019a,b;Dempsey et al. 2020).
the excitation torque density obtained using our linear calculations, see Section 2.2 for details.One can see prominent peak and trough2 of d ex /d located just outside and inside of the planetary orbit, at | −  p | ∼  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 d ex /d ∝ | −  p | −4 for | −  p | ≳  p .Some of these statements have been refined recently.In particular, the two-dimensional (2D) hydrodynamic simulations of Dong et al. (2011a) carried out in the local (sharing sheet) approximation found that d ex /d does not decrease monotonically but actually reverses its sign around | −  p | ≈ 3.2 p .This 'negative torque density' phenomenon is illustrated in the inset in Fig. 1(b), where d ex /d first turns negative outside the main positive peak in the outer disc,  >  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 d ex /d beyond this point follows d ex /d ∝ | −  p | −4 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 nonlinear simulations of Duffell & MacFadyen (2012) and Kley et al. (2012).
Recent global simulations with radially extended domains revealed even more complexity.In particular, 3D simulations by Arzamasskiy et al. (2018) have shown that d ex /d exhibits multiple sign reversals at | −  p | ∼  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 & 8) and simulations of Miranda & Rafikov (2019b, see their Fig.1), visible as the oscillations of the density wave angular momentum flux   (its radial derivative is equal to d ex /d in the absence of wave damping), which are especially pronounced in the outer disc.Sign reversals of d ex /d 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. 1b,d also shows such d ex /d features very clearly both as multiple, rather regularly spaced wiggles in the outer disc for  ≳ 2 p as well as 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 d ex /d.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 d ex /d oscillations grows as the mass of the perturber increases and the disc-perturber coupling becomes nonlinear, see Dempsey et al. (2020).Thus, clarifying the origin of the features of d ex /d behaviour in the linear regime will also help us understand the torque density behaviour in the nonlinear regime, relevant for circumbinary discs around stellar and supermassive black hole binaries (Cimerman & Rafikov, in prep.).
Our work is organized as follows.After describing our setup 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.

Model setup
We consider a 2D model of a (initially axisymmetric) razor-thin disc that orbits a central star of mass  ★ .We adopt polar (, ) coordinates centred on the central star.The disc has aspect ratio ℎ ≡ / ≪ 1, where  =  s /Ω is the vertical scale-height,  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  p , moving on a circular orbit with radius  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 discplanet coupling that also neglected the indirect potential (Bate et al. 2003;D'Angelo & Lubow 2008, 2010;Duffell & MacFadyen 2012;Dong et al. 2011a,b;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 where  p =  s ( p ) is the sound speed at  p , Ω p is the orbital angular frequency at  p ,  p =  ( p ) and ℎ p = ℎ( p ).The assumption  p ≲  th implies that the disc response to the planetary gravity is linear, i.e. Σ/Σ ∼  p / th ≪ 1, where Σ = Σ − Σ 0 is the planet-induced perturbation of the surface density Σ relative to the background surface density Σ 0 ().
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 Ω 2 K = √︁   ★ /3 is the Keplerian orbital frequency, and  is the gas pressure (see below).The unperturbed velocity of gas is given by  ,0 = 0,  ,0 () = Ω 0 (), with Ω 0 () defined by equation (5) using Σ 0 ,  0 .
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  is the adiabatic constant.In our simulations we explicitly solve the energy equation, see Appendix A, thus properly accounting for the evolution of  (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 locallyisothermal EoS: it is used, together with the equation ( 2), to set the initial radial profile of  in our adiabatic calculations.The adiabatic sound speed is  s,adi =  1/2  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 angular momentum flux (AMF) associated with these waves, defined as with velocity perturbations   and   (, ) =   (, )− ,0 ().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 nonlinear.Linear calculations of the integrated one-sided Lindblad torque (GT80) find the characteristic magnitude of the AMF of the planetdriven density waves to be which we will use as a reference value for   .To allow for meaningful comparison, radial torque densities are given in units   ,0 / p , as in Fig. 1.

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  p / th ≪ 1.The method 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.
Second, we also perform fully non-linear simulations (see Appendix A2 for details and numerical parameters) of the problem using Athena++ 3 (Stone et al. 2020)  p = 0.01 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 Σ(, ) 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 d) by integrating R × dF/d (with dF being the direct gravitational force exerted by a planet onto a disc element d = dd) over  at fixed  as which points in the -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  and  and aspect ratio at the planet radius ℎ 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  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 ℎ p = 0.1 and  =  = 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 & 4.

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 ℎ p = 0.1 (higher than ℎ p = 0.05 used in Fig. 1) and  =  = 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 (,  lin ()), where4 (Rafikov 2002a;Ogilvie & Lubow 2002) 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  >  p the spiral maintains a sharp crest (red), neighboured by a lower The colormap is the same as in Fig. 2. Azimuthal slices in panels (c) and (d) correspond to the radii indicated by the horizontal coloured dashed lines (of corresponding color) in panels (a),(b),(e) and (f).Pink tripods indicate locations where the peak of the wake crosses the line connecting the central star and the planet ( peak =  p ) as in Fig. 2; here we also show them in the inner disc.See text for details.amplitude trough (blue).But in the inner disc this approximation starts failing for  ≲ 0.6 p , as the spiral arm splits into multiple peaks 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 difference5 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.
To further support our statements, in Fig. 3 we show the same 2D map of the relative surface density perturbation Σ/Σ 0 as in Fig. 2 but now in Cartesian geometry, see panels (a) and (b) for inner and outer disc, respectively.Also, in panels (c) & (d) we show the azimuthal profiles of Σ in the outer and inner disc, respectively, at several fixed radii (shown in panels (a) & (b) using dashed lines of the same color).They are normalized by the characteristic wave amplitude Σ lin () expected in the linear regime (see Section 4.2) to highlight the evolution of the shape of the density wake.Finally, in panels (e) & (f) we show the radial profile of the excitation torque density d ex /d 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.

Outer disc
In the outer disc ( >  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 d ex /d.
At  = 2.17 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 6 d ex /d, which is what we see in panel (e).
As the wake crosses the line  =  p (marked with a tripod), the wake in panel (c) shifts to  <  p , pulling back on the planet, which gives rise to a positive peak of d ex /d at  = 2.28 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  = 2.5 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) d ex /d.
Beyond that point d ex /d becomes negative again and at  ≃ 2.8 p (red lines) the wake profile in panel (c) almost coincides azimuthally with that at the previous minimum of d ex /d (blue curve).As a result, d ex /d 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 d ex /d (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 d ex /d is only quasi-periodic: closer to the planet d ex /d 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 d ex /d is very accurately set by the wake crossings, thanks to the essentially self-similar shape of Σ in the outer disc.

Inner disc
We now turn to the inner disc ( <  p ), see panels (b), (d), (f).Close to the planet, at | −  p | ≃ (2 − 3) p = (0.2 − 0.3) 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  = 0.8 p (sampling the main trough of d ex /d), which is similar in shape (upon -reflection and rescaling) to the curves in panel (c). 6Even despite the smaller, leading trough at  −  p ≃ 0.2 , which provides a small but positive contribution to d ex /d.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  ≲ 0.6 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 et al. (2017) and Miranda & Rafikov (2019a).These modifications are also reflected in the Σ/Σ lin profiles for  = 0.265 p , 0.18 p , 0.118 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,9).
The more complicated and steadily evolving (in ) azimuthal profiles of Σ lead to the loss of coherence of d ex /d at consecutive wake crossings in the inner disc, e.g. at  = 0.265 p (pink) and  = 0.118 p (magenta).As a result, one no longer sees the radial quasi-periodicity of the d ex /d features, which are so obvious in the outer disc.Moreover, the amplitude of these features is also significantly reduced, partly because of the lower angular momentum flux 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 & 4.5.

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  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  <  p .
To emphasize the role of self-similarity of Σ even further, in Fig. 4 we show d ex /d computed in the outer disc for a hypothetical 'Gaussian' wake in the form where  is an (arbitrary) amplitude and we set  = ℎ p as the charac-teristic azimuthal width of the wake.This wake has a correct radial scaling of its amplitude to ensure the conservation of the wave angular momentum flux (see Section 4.2) and follows the curve  =  lin () 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 d ex /d (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 () behaviour (setting the radial frequency of wake crossings) that determines it.

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, that 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 Σ(, ) = Σ(, ) − Σ 0 ().Since Σ 0 is axisymmetric and gives no net contribution to the torque, we can replace Σ(, ) with Σ(, ) 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  ≡  p / < 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 7 Up to a normalization, we recognize the remaining integral as the cosine Fourier coefficient of Σ/: 7 The term with  (0) such that we can write equation ( 20) as where all information about the surface density perturbation is contained in the    ().In the inner disc, defining  ≡ / 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   Σ) such that our results ( 20) & ( 22) are fully general.

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. 3(c)), 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  Σ() 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 (, φ()), and the only change of  with  is a translation in the -direction given by φ().
With the ansatz (24) we have where  ′ () = d()/d.Let us now introduce the complex Fourier coefficients of  ′ as follows: Writing Ψ  =    i  , with   = |Ψ  | and   = arg Ψ  being the amplitude and the phase of Ψ  , respectively, we obtain Substituting this expression into equation ( 20) gives, after straightforward manipulation, Note that here φ() is understood to lie in the interval (0, 2), i.e. is unique mod 2, since in equation ( 25)  ∈ (0, 2).This is the final expression for the torque density given the selfsimilar ansatz (24).Once the forms of Σ(), φ() and Ψ  (i.e.  and   ) are specified, the torque density can be computed by simple summation.Knowledge of Ψ  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 selfsimilar shape, see Fig. 4a of that paper that reveals only a weak evolution of the wake shape for  >  p .This is also obvious from our Fig.3c 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. 3d also makes this clear since the wake profiles for different  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 ℎ 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.   / = 0.It has been shown (Rafikov 2002a) that conservation of   leads to the wave amplitude scaling as Σ() ∝ Σ lin (), where In a power-law background disc, see equations ( 3) & (2), we find 8 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   varies as   ∝  2 s () (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  p / s () in equation ( 30) from 3/2 to 1/2.Similarly, thermal relaxation (e.g.-cooling) always damps   (Miranda & Rafikov 2020a,b), reducing Σ() 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 () given by equations ( 10)-( 14). 8The | (/ p ) −3/2 −1| factor is missing the absolute magnitude in equations (37) of Rafikov (2002a) and ( 16) of Miranda & Rafikov (2019a).It is correctly present in the function  () in Rafikov (2002a).10), ( 11).The grey dotted vertical line corresponds to the same radius as in Fig. 6, where   peak and the black dotted line show where |  Σ | peak (ℎ p ≃ 0.4).This plot illustrates the radial universality of amplitudes and evolution of phases of the planet-driven density wave.See Section 4.3 for details.

Properties of Ψ 𝑚 in AMF-preserving discs
Next, we explore the properties of the Fourier coefficients Ψ  for the self-similar outer density waves in AMF-preserving discs.We employ the ansatz ( 24), (32) and consider an adiabatic disc with ℎ p = 0.1 and  =  = 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 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 () to account for the overall wrapping of the spiral wake (panel (c)), a clear pattern becomes obvious, showing only a slow evolution with .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   = |Ψ  | and phase   = arg Ψ  of the coefficients Ψ  computed using  Σ via equation (35).In Fig. 6 we show them as a function of ℎ 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 & b), we notice that   peak at around ℎ p ≃ 1.45 (ℎ p ≃ 1.2 p ) in the outer (inner) disc and do not show much variation with .This is expected for free linear waves in flux-conserving (adiabatic and barotropic) disc models, since the radial scaling of Σ with Σ lin () is already taken into account in ansatz ( 24), (32).As we show in Section 6, the situation is very different for non-flux-conserving discs.
On the other hand, the phases   (panels c & d), which include the phase shift by  lin () through equation ( 35), as in Fig. 5c, show a substantial evolution with .In the inner disc, the phases spread out with  as the distance from the planet increases.Around  = 0.1 p ,   changes by 2 every Δ(ℎ p ) ≃ 2. On the other hand, in the outer disc,   change a lot less with  and, in particular, show very weak variation with  at low  (similar to Fig. 5c).
Panel (f) shows the amplitudes   multiplied with Laplace coefficients  () 1/2 (), 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 ℎ p , in order to zoom into the relevant region.For  closest to the planet,    () 1/2 peak at around ℎ p ≃ 0.4 − 0.5 and fall off rapidly at large .As the distance from the planet increases, the location of the peak shifts towards lower , with the overall amplitude of    () 1/2 decaying.This is expected from the asymptotic behaviour of the Laplace coefficients (discussed in Section 4.5) and the weak dependence of   on .In other words, the behaviour of the Laplace coefficients is the key factor determining which  provides the largest contribution to d ex /d.Vertical dashed lines in all panels of Fig. 6 correspond to  c defined as the critical  for which the product    () 1/2 falls off below one percent of its maximum value, such that only  <  c need to be considered when computing d ex /d (the values are  c = 14, 9, 6, 4 for  = 1.5, 2, 3, 5 p , respectively.).These vertical dashed lines move to lower  as | −  p | increases, implying that closer to the planet, more modes contribute to d ex /d, while for larger distances, only the lowest few  matter; see Section 4.5 for more details.
Similarly, in the inner disc we multiply   with  () 1/2 () (panel (e)), see equation ( 23).Because of the extra factor of  (compared to the outer disc) the resultant curves precipitously diminish in magnitude as  decreases, explaining the weakness of wiggles in the inner disc.Dashed lines again mark   , for which    () 1/2 () fall to 1% of the maximum value.
The behaviour of the phases   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 & 9.But the variation of   at high  in panel (d) suggests that self-similarity may not work in the outer disc either.However, note that in panel (d) the variation of   with  is not very dramatic for  <  c (to the left of the corresponding vertical lines), i.e. for all azimuthal harmonics that matter for the calculation of d ex /d.This allows us to consider Σ as approximately self-similar in the outer disc, which will be used next.

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 d ex /d using eq.( 36) eq. ( 38) Figure 7. Illustration of the convergence of the series (36) to the exact expression (29).Top: The blue solid line shows the full linear solution using  max = 100 modes with (   ,   ) varying with radius (see Fig. 6), i.e. without using the self-similar approximation.Dashed lines show the approximation (36) truncated at different  tr <  max (colours) with (   ,   ) sampled at  cal = 3 p (see Fig. 6).Dot-dashed grey and dotted magenta curves illustrate two predictions for the envelope of the torque wiggles given by the equations ( 37) and (39), respectively.Bottom: deviation Δ of the truncated series (36) from the exact linear solution.Agreement is best around  cal and the truncated series rapidly converge with  tr to the exact solution, especially at larger .
equation ( 29) with   ,   obtained under the assumption of a selfsimilar wake, see Section 4.3.Since Fig. 6d shows some evolution of   with , we pick a particular calibration radius  cal = 3, at which we measure   ,   and assume these values to not depend on  (an approximation which should work reasonably well as we argued earlier) when computing the radial profile of d ex /d.The resultant torque density is then compared with the d ex /d obtained by solving the linear problem exactly.
In order to highlight the contribution of different individual modes to d ex /d, we also use a truncated (at  =  tr ) version of the series expansion (29): where here we have made explicit that under our self-similar approximation the values of (   ,   ) are taken at a calibration radius  cal .
In Fig. 7 we show the exact numerical solution for d ex /d (blue solid) and the truncated self-similar series d tr ex /d (36) for several values of  tr in panel (a); the deviation between them is shown in panel (b).One can see that close to the planet ( ≤ 2 p ) the selfsimilar approximation deviates significantly from the global linear solution for all  tr .This is to be expected since   still show some dependence on  (compared to their values at  cal ) for  ≲  c close to the planet, see Fig. 6d.Further away, the first trough of d ex /d is matched well both in shape and amplitude for all  tr > 2.
In the radial interval between this trough and the next one, we see rather small deviations for all  tr .Beyond the second trough, the self-similar curves closely track the full solution for  tr ≥ 2. For  ≥ 4 p , the curve for  tr = 2 merges with those for higher  tr .The outermost torque oscillation period is matched reasonably well even for a single mode when  tr = 1.Panel (b) reveals that for  ≥ 4 p , the self-similar approximation shows a periodic deviation from the exact solution, with the amplitude that steadily decays as  tr increases.
These results indicate that for sufficiently large distances from the planet,  tr = 5 − 10 modes are sufficient to describe the behaviour of d ex /d, if the Ψ  are sampled at an intermediate distance  cal within the radial region of interest.At the quantitative level these results still depend on the value of  cal , but this dependence should not curtail the good performance of the self-similar approximation described in Sections 4.1,4.2 in reproducing d ex /d.

Understanding the key features of torque wiggles
As demonstrated in the previous subsection, equation ( 29) is able to successfully reproduce the correct behaviour of d ex /d.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 its maximum possible value This very simple constraint is plotted as the grey dot-dashed curve in Fig. 7a and it works remarkably well, with the troughs of peak wiggles almost touching this curve (see more on this below).Second, equation ( 29) lets us analyze the torque density behaviour in the asymptotic limit  ≫  p (or  ≪ 1).Recall that in thin discs, / ≪ 1, the Fourier amplitudes of planet-driven Σ  peak at  * ∼  p / p (i.e. * ℎ p ∼ 1) and exponentially decay for  ≳  * (GT80).The coefficients Ψ  have an extra factor of  compared to Σ  , leading to   peaking at  slightly higher than  * , which can be seen by comparing Figs. 5 (see black and grey dotted lines) and 6.When computing d ex /d via formula (29),   are multiplied by the Laplace coefficients  () 1/2 , which asymptotically behave far from the planet as (Murray & Dermott 1999) This suggests that as | −  p | increases, the main contributions to d ex /d should be provided by the low- components of Σ.
Eventually, far from the planet, only the  = 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 (1) 1/2 () →  as  ≪ 1, we obtain the following asymptotic prediction for the outer envelope of d ex /d: where the dependence of Σ lin on  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  ≫  p ; at intermediate  it underestimates the envelope (37) of d ex /d because of still significant contribution of other  harmonics.
Note that in the inner disc close to the centre, in the asymptotic limit  = / p → 0, torque density is suppressed by an additional factor of : d ex /d ∼ Σ 1 (/ p ) 2 , where Σ 1 is the  = 1 Fourier component of Σ, see the equation ( 23).As a result, torque contributions of different harmonics fall off very rapidly as  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  <  p .
Third, equation ( 29) also allows us to explain why in the outer disc d ex /d features sharp (negative) peaks at radii, where the density wake crosses the line  = 0 between the star and the planet.Indeed, note that   in Fig. 6d exhibit rather small spread Δ  ≲ /2 for all  and  ≤  c .This allows us to assume that   is approximately the same for all  that contribute significantly to the sum in equation ( 29).In that case, these significant terms interfere constructively whenever  lin () ≲ 1 for all  ≲  c (this also explains why the troughs of d ex /d overshoot the envelope (39) at intermediate ).As a result, whenever  lin () → 0 mod 2, i.e. the wake (centred on  lin ()) crosses the line  = 0 from the star to the planet, constructive interference of many terms in (29) leads to a sharp increase of the amplitude of d ex /d, 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 |d ex /d| given by equation ( 37).
Fourth, we can also use equation ( 29) to estimate the radial width of these sharp features of d ex /d (although the outcome of this exercise may be of limited utility).Indeed, constructive interference of harmonics with  ≲  c leading to torque wiggles requires  lin ≲  −1 c , constraining the radial width of the sharp troughs of d ex /d as see equation ( 11).Because of the rapid evolution of  c with  (see Fig. 6) it may be somewhat challenging to extract a more explicit dependence of  on  and disc parameters.Nevertheless, it is clear that  should be smaller in thinner discs with lower  s and ℎ p .

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  and  -the slopes of the background surface density and temperature (which set the initial profile of entropy), see equations ( 2)-( 3) & (30), while  lin also depends on ℎ p , see equations ( 4) & ( 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 , , and ℎ p , and then move on to study the effect of these parameters on the torque wiggles in Section 5.2.

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 () defined by the equation ( 24) as a function of the (naturally scaled) azimuthal variable ( −  lin )ℎ −1 p .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.

Variation of the slopes 𝑝 and 𝑞
In Fig. 8 we show  in the inner (left) and outer (right) disc with ℎ p = 0.1 for different values of  and  (colours) at several values of  (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  suggestive of its approximate self-similarity.The variation of  and  introduces some deviations of () 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 , which directly affects the emergence of additional spiral arms.
To summarize, the variation of  and  does not greatly modify the overall picture of wake propagation (for ℎ p = 0.1): the outer wake maintains a self-similar shape with good accuracy, while the inner wake does not follow this universal pattern.
. Same as Fig. 8 but now illustrating the effect of the disc scaleheight ℎ p variation on the wake profile evolution.The coordinate on the -axis has been rescaled ∝ ℎ −1 p , such that the wakes have similar horizontal scale.Note that for lower ℎ p the dispersion and splitting of the wave accelerate and selfsimilarity gets violated, even in the outer disc.See text for details.

Variation of ℎ p
We now investigate the effect of varying ℎ p on the wake shape, keeping  =  = 1 fixed.In Fig. 9 we show  for several values of ℎ p , similar to Fig. 8.While the rescaling with Σ lin and stretching in  by a factor ∝ ℎ −1 p allow us to compare the wake shape on similar scales in this plot, we see clear deviations between the  profiles for different ℎ p , even close to the planet (panels a & d).
As the distance from the planet increases, the wake suffers from the dispersive effects, which are especially pronounced for small ℎ 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  = 0.1 p for ℎ p = 0.025, while the ℎ p = 0.1 disc shows only two clear spiral arms.
What is more remarkable, for very low ℎ p a secondary arm starts to develop even in the outer disc.This can be noticed by comparing the ℎ p = 0.025 curves at  = 1.1 p and 3 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 ℎ p .
To summarize, ℎ p is a very important parameter affecting the density wave propagation more significantly than the changes in  and .

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 d ex /d com-puted 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 labeled 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,  = 1, 2, ...; one can see that the two agree very well.Note that panels (a)-(d) all have ℎ p = 0.1.
Comparing panels (a) and (b) we see that changing the surface density slope  does not affect the radial periodicity of d ex /d, i.e. the positions of prominent troughs (see also dashed vertical lines).This is expected since  lin does not depend on Σ 0 .However, the amplitude of the torque wiggles changes and in a non-uniform fashion:  = 0 disc shows lower (higher) d ex /d in the inner (outer) disc compared to  = 1.This can be understood from the dependence of the excitation torque density on Σ 0 (): d ex /d ∝ Σ ∝ Σ lin ∝ Σ 1/2 0 (), see equation ( 16).Thus, given their approximately identical () for discs with different  (see Fig. 8d-f), we expect that curves of Σ −1/2 0 d ex /d would have a universal shape.This is indeed the case, as illustrated by Fig. 11 for  = 1 (solid) and  = 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  = 0.5 and a constant temperature ( = 0).It is clear that the variation of  leads to changes of not only the amplitude of the wiggles and but also their radial periodicity -the distance between the consecutive troughs increases (decreases) in the outer (inner) disc as  goes down.This is explained via the behaviour of  lin , which determines the global shape of the wake.Regarding the amplitude, equation ( 30) predicts d ex /d ∝ Σ lin ∝ (/ p ) 3/4 in AMF-preserving discs, such that lower  should increase (decrease) d ex /d in the inner (outer) disc.This is consistent with the trends seen in panels (c) and (d).
Finally, in panels (e), (f) we consider disc models with lower ℎ p , while keeping  =  = 1 as in panel (a).Comparing the locations of wake crossings with the radial structure of d ex /d it is evident that when ℎ p is halved, the number of wake-crossings doubles, as expected from  lin ∝ ℎ −1 p .Closer examination of panel (e) reveals that in the outer disc with ℎ p = 0.05, wake-crossings are again aligned with minima of d ex /d 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 ℎ p is decreased (see Fig. 9) has a direct effect on d ex /d.For even lower ℎ p = 0.025 (panel f), the wake-crossings in the outer disc align with the dominant peak or trough of d ex /d 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 ℎ p , see Section 5.1.2.

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 .Torque density profiles obtained using linear theory (blue), and direct Athena++ simulations with  p = 0.01 th (orange), for different disc parameters (the value of the parameter varied relative to the fiducial model is indicated in each panel).Vertical magenta (dashed) lines correspond to the radial locations where the peak of the wake crosses the line  =  p from the star to the planet, while the green (dotted) lines correspond to the analytical prediction  lin = 2 ,  = 1, 2, ....For the smallest value ℎ p = 0.025, no Athena++ simulation was performed.Note that axes are broken and log-scale is different on left and right hand side.Grey shaded areas indicate damping zones in the outer disc in Athena++ simulations (inner disc damping zones are outside of plot range).
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 −3 to 10 2 , 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 d ex /d 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  s ≈  s,adi =  1/2  s,iso -i.e. the adiabatic sound speed, which because of the factor  > 1 is higher than the isothermal sound speed  s,iso , to which  s reduces in the opposite limit  ≪ 1.As a result, equation (10) predicts slower evolution of  lin with , 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,  = 1, 2, ..., which are computed using  s ≈  s,adi in (10) for  ≥ 1 and using  s ≈  s,iso for  < 1; these lines are well aligned with the troughs of d ex /d.
The amplitudes of torque wiggles are quite similar in the adiabatic ( = 10 2 ) 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 .These trends can be easily understood by examining the behaviour of the amplitudes   and phases   of the Fourier coefficients Ψ  (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 colored lines at  cthe maximum  significantly contributing to d ex /d -at every , see Fig. 6 and Section 4.3.
For  = 10 2 (panel (a)), close to adiabatic limit,   and   1/2 0 , to eliminate the scaling of d ex /d with Σ lin for fiducial  = 1 (solid) and for  = 0 (dashed).This illustrates that in the linear regime the overall pattern of the torque density remains unchanged by varying the surface density slope, i.e. the interference of the modes is unchanged, only their individual amplitudes are rescaled with Σ lin , see equation (30).See text for details.
behave as they do in Fig. 6 -the former are independent of , while the latter slowly evolve with .As a result, the wiggle pattern is essentially indistinguishable from that of Fig. 10a.For  = 1 we see   being uniformly (in ) suppressed as  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   to the left of the corresponding  =  c lines.This reduces the amplitude of the wiggles seen in Fig. 12b.
Wave decay is most dramatic for  = 0.1 (Fig. 13c), leading to the precipitous drop of   for  <  c and the corresponding rapid decay of the torque wiggle amplitude with  in Fig. 12c.Moreover, for this value of  the phases   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 −3 , see Fig. 12d, is due to the disc starting to transition to the locally isothermal regime.Fig. 13d shows that   decays with  slower than in panel (c), in particular, the most significant   () (for  <  c ) for every  are not that different even from panel (a).This agreement is improved even further in the locally isothermal limit shown in Fig. 13e (although here the suppression of   with  is more uniform9 in , more akin to panel (b)), explaining the qualitative similarity of d ex /d 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.

DISCUSSION
In this work we provided a detailed exploration of the properties and origin of the torque wiggles -a new 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 d ex /d ∝ | −  p | −4 picture that has been prevalent in the field for decades (e.g.Armitage & Natarajan 2002;Chang et al. 2010;Zagaria et al. 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 d ex /d (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 limit10 .This is consistent with the calculations of Rafikov & Petrovich (2012) who showed that the negative torque density phenomenon (i.e.no additional d ex /d reversals) is the only possibility far from the planet in the shearing sheet.It is important to emphasize that the torque wiggles are not related to the low-order Lindblad resonances (Goldreich & Tremaine 1980) which are not present far from the planet.
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 potential11 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.
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  <  p (another one is the faster radial decay of d ex /d 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 ℎ p ).We note that in the presence of a well-developed turbulence, the wake coherence may be destroyed even in the outer disc (Zhu et al. 2013), resulting in the suppression of the torque wiggles (although they may still be evident after time-averaging).amplitude; this is supported by our preliminary runs in high  p regime which we do not show here.Nevertheless, Fig. 1 of Miranda & Rafikov (2019b) still shows some torque wiggles up to  p =  th , see also Dempsey et al. (2020).When wave nonlinearity is strong, the shape of the torque wiggles may become sensitive to the slope of the surface density  (unlike in the linear case, see Fig. 11), since the rate of nonlinear evolution of the wave depends on it (Rafikov 2002a;Cimerman & Rafikov 2021).Figure 14.Top: Azimuthal profiles of the relative surface density perturbation obtained using Athena++ simulations with fiducial parameters and  p = 0.01 th (solid) and linear calculations (dashed, falling almost on top of each other) at locations where  lin =  p .Even for a low  p , non-linearity leads to substantial deviations from linear theory at large distances from the planet (see also Cimerman & Rafikov (2021)).Bottom: Amplitude   of Fourier components of Ψ  for the non-linear Athena++ simulations as a function of wave-number , measured at the same radii as in the top panel; vertical dashed mark  c .As opposed to the linear solution (cf.Fig. 6b), there is a clear evolution of   () with  leading to efficient transfer of power to low  = 1, 2 modes, as well as the overall decay at high  due to shock-damping.

Effect of torque wiggles on the total torque
One can see that close to the planet  ex () follows   () 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  ex () in the nonlinear case.Despite that, both  ex () and   () 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 Δ =  ex () −  ex ( r ) -the difference between the cumulative torque at  and the cumulative torque at some reference  r .We pick somewhat arbitrarily  r,i = 0.7 p and  r,o = 1.3 p in the inner and outer disc, respectively (dotted black vertical lines), just outside the main peaks of d ex /d.In the inner disc, the total torque hardly changes and Δ remains close to zero for  <  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  ex is Δ/ ≃ −0.03 and its location coincides with the first wake crossing of  p around  ≃ 2.2 p .These oscillations eventually die out as  → ∞, however their integrated effect is non-zero and Δ (∞)/ ≃ −0.015, see panels (b) and (c).Thus, the integrated effect of the torque wiggles in the outer disc is to slightly lower  ex from the maximum value near  ≈ 1.3 p .
In our models, the highest relative variation of  ex in the outer disc due to torque wiggles is Δ/ ex ≃ 0.07 (after the first wakecrossing) for  = 0,  = 0, ℎ p = 0.1 disc.Clearly this effect is  b) Integrated (total) torques  ex () (dashed) and AMF   () (solid) obtained using linear theory (blue) and simulations (orange).(c) Global variations in the total torque relative to its value at some reference radius  r , see text for details.As before, the vertical lines mark radii where either  lin =  p (green) and at the position of the peak of Σ,  peak =  p (purple).See Section 7.2 for details.small, as is the amplitude of torque wiggles in general, as the main peak of d ex /d within a couple of scaleheights from the planet provides the dominant contribution to  ex , see Fig. 1.Moreover,  ex is the integral of the radially oscillating d ex /d see equation ( 29), which additionally lowers Δ.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 d ex /d, 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 will be presented in Cimerman & Rafikov (in prep.).

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 (| −  p | ≳  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 below.
• 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  2 p , just as that of the full planetary torque.
• We developed analytical theory (Section 4) showing that quasiperiodic wiggles arise because of the global nature of the planetdriven 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 ℎ p → 0).
• This theory demonstrates that the behaviour of d ex /d at | −  p | ≳  p is dominated with good accuracy by a small number of low- Fourier components of the perturbation Σ (Section 4.4).
• This theory also explains the key features of d ex /d behaviour (Section 4.5).In particular, it accounts for the sharp features in d ex /d 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- 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.
Future work (Cimerman & Rafikov, in prep.) will explore torque behaviour in the non-linear regime of circumbinary discs, when the perturber-to-star mass ratio  ∼ 1.We will demonstrate that the torque wiggles become much more important in this regime.
in the fiducial setup (when ℎ p = 0.1) that are logarithmically spaced in radius and uniformly spaced in azimuth.For cases with ℎ p = 0.05, this resolution is doubled.We do not perform an Athena++ simulation for ℎ 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 ≤ / p ≤ 0.15 and 4.5 ≤ / p ≤ 5 in order to avoid reflections (see Cimerman & Rafikov 2021).where  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) (, ,  + Δ) = exp (−Δ/ rel ) (, , ) (B2) 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  s =  s,iso ().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 2 to achieve a situation close to an adiabatic limit.We tested and confirmed that these runs give virtually identical results to even greater  = 10 3 , indicating convergence towards the adiabatic limit, in line with the findings of Miranda & Rafikov (2020a).
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure1.Overview of the typical radial distribution of the excitation torque density d ex /d (normalized according to the equation (8)) using linear calculation (Section 2.2) for disc parameters ℎ p = 0.05 and  =  = 1 (se Section 2.1).Top and bottom panels show the same data on different vertical scales.The arrows mark the main peaks ('main') of d ex /d and indicate the first sign change of d ex /d outside the main peak, also known as the negative torque density phenomenon ('NT', seeDong et al. 2011a; Rafikov & Petrovich 2012).Far away from the planet the label 'torque wiggles' marks the region where d ex /d exhibits multiple sign changes(Arzamasskiy et al. 2018; Miranda & Rafikov 2019a,b;Dempsey et al. 2020).

Figure 2 .
Figure 2. Polar plot of Σ/Σ 0 × (  p / th ) −1 for the fiducial disc model.The grey dashed line goes through the origin where the central star is located and the planet.Pink crosses mark the locations where the peak of the spiral density wave  =  peak () crosses the line  =  p ; we only show them in the outer disc for clarity.The planet and disc rotate in a clockwise direction.

Figure 3 .
Figure 3. (a),(b) 2D Map of the surface density perturbation obtained by linear solution in the fiducial disc model (ℎ p = 0.1,  =  = 1) together with (c),(d) azimuthal slices of the surface density perturbation at several radii, and (e),(f) the radial profile of torque density, in the inner and outer disc respectively.The colormap is the same as in Fig.2.Azimuthal slices in panels (c) and (d) correspond to the radii indicated by the horizontal coloured dashed lines (of corresponding color) in panels (a),(b),(e) and (f).Pink tripods indicate locations where the peak of the wake crosses the line connecting the central star and the planet ( peak =  p ) as in Fig.2; here we also show them in the inner disc.See text for details.

Figure 4 .
Figure 4. Excitation torque density in the outer disc for our toy Gaussian model (blue) and the full linear solution of the linear problem (orange).Note that they share the same periodicity, confirming that it is the behaviour of  lin () that determines the radial structure.The toy model produces a more symmetric pattern with neighbouring minima and maxima being of almost identical height.The torque density associated with the Gaussian profile vanishes exactly when  peak =  lin , as expected from symmetry.

Figure 5 .
Figure 5. Amplitudes (a) and phases (b) of  Σ -the Fourier components of the surface density perturbation in the outer disc; panel (c) shows phases with an additional shift by  lin (R), see equations (10), (11).The grey dotted vertical line corresponds to the same radius as in Fig.6, where   peak and the black dotted line show where |  Σ | peak (ℎ p ≃ 0.4).This plot illustrates the radial universality of amplitudes and evolution of phases of the planet-driven density wave.See Section 4.3 for details.

)Figure 6 .
Figure 6.Amplitude of Ψ  (  , left), phase of Ψ  (  , middle) and   multiplied by the corresponding Laplace coefficient (right) and  in panel (e) as a function of ℎ p for the fiducial disc for several radii (colours) in the inner (top) and outer (bottom) disc.The dashed vertical lines indicate  c ℎ p , such that only the modes with  <  c (to the left of those lines) contribute significantly to d ex /d.See text for details.

Figure 8 .
Figure 8.The dimensionless function  (  ) describing the azimutal shape of the wake in the self-similar anzatz (24) sampled at different  (panels) in the inner disc (left column) and outer disc (right column) for different values of the surface density and temperature slopes  and  (colours).The coordinate on the -axis has been rescaled ∝ ℎ −1 p .See text for discussion.
Figure10.Torque density profiles obtained using linear theory (blue), and direct Athena++ simulations with  p = 0.01 th (orange), for different disc parameters (the value of the parameter varied relative to the fiducial model is indicated in each panel).Vertical magenta (dashed) lines correspond to the radial locations where the peak of the wake crosses the line  =  p from the star to the planet, while the green (dotted) lines correspond to the analytical prediction  lin = 2 ,  = 1, 2, ....For the smallest value ℎ p = 0.025, no Athena++ simulation was performed.Note that axes are broken and log-scale is different on left and right hand side.Grey shaded areas indicate damping zones in the outer disc in Athena++ simulations (inner disc damping zones are outside of plot range).

Figure 13 .
Figure 13.Amplitudes   and phases   of Ψ  (similar to Fig. 6) sampled at different  >  p (shown at the top) for disc models with -cooling and the fiducial disc parameters  =  = 1, ℎ p = 0.1.The cooling time  (indicated in each row) takes the same values as in Fig. 12.The vertical dashed lines have the same meaning as in Fig. 6.

Fig. 15
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))  ex () =

Figure 15 .
Figure 15.(a) Radial excitation torque density d ex /d (same data as in Fig. 10a.(b) Integrated (total) torques  ex () (dashed) and AMF   () (solid) obtained using linear theory (blue) and simulations (orange).(c) Global variations in the total torque relative to its value at some reference radius  r , see text for details.As before, the vertical lines mark radii where either  lin =  p (green) and at the position of the peak of Σ,  peak =  p (purple).See Section 7.2 for details.