Planetesimal Drift in Eccentric Disks: Possible Outward Migration

Radial drift of solid particles in the protoplanetary disk is often invoked as a threat to planet formation, as it removes solid material from the disk before it can be assembled into planets. However, it may also concentrate solids at particular locations in the disk, thus accelerating the coagulation process. Planetesimals are thought to drift much faster in an eccentric disk, due to their higher velocities with respect to the gas, but their drift rate has only been calculated using approximate means. In this work, we show that in some cases previous estimates of the drift rate, based on a modification of the results for an axisymmetric disk, are highly inaccurate. In particular, we find that under some easily realized circumstances, planetesimals may drift outwards, rather than inwards. This results in the existence of radii in the disk that act as stable attractors of planetesimals. We show that this can lead to a local enhancement of more than an order of magnitude in the surface density of planetesimals, even when a wide dispersion of planetesimal size is considered.


INTRODUCTION
Rapid radial drift of solid particles in protoplanetary disks has been identified as a major barrier to planet formation because it leads to high collision velocities and depletion of solid material in the disk.Often referred to as the "meter-size barrier", radial drift affects larger objects as well, with typical disk models yielding drift timescales under one Myr for objects of hundreds of meters in size at one AU (Adachi et al. 1976;Weidenschilling 1977).
In the presence of a nearby binary companion, objects of tens or even hundreds of kilometers in size experience substantial radial drift because the binary companion excites eccentricity in both the disk and the planetesimal.Rafikov & Silsbee (2015b) estimated that the presence of a companion star on an eccentric orbit at 20 AU increases the inspiral rate at 1 AU (for bodies larger than a few kilometers) by factors of 10's to 100's compared with the case of a circular disk.However, this estimate was made in a rather approximate way, by using the formulae derived in Adachi et al. (1976) for an axisymmetric disk and replacing the planetesimal eccentricity with the relative eccentricity between planetesimal and gas.
The planet formation rate is known to be suppressed around stars that have a companion within ∼ 200 AU (Kraus et al. 2016;Moe & Kratter 2021).Multiple factors may lead to this suppression.
Gravitational perturbations from the companion stir up the planetesimal population, potentially resulting in destructive collisions (e.g.Heppenheimer 1978;Paardekooper et al. 2008;Rafikov & Silsbee 2015a).In addition, the masses and lifetimes of circumstellar disks are generally reduced in tight binary systems compared with single stars (Harris et al. 2012;Akeson et al. 2019;Barenfeld et al. 2019).
The role of enhanced planetesimal drift in tight binary systems has not been investigated in detail.The effect was included in the multi-zone coagulation simulations done by Silsbee & Rafikov (2021), using the same prescription for drift as in Rafikov & Silsbee (2015b).They found radial drift to have a non-zero, but not dominant effect on the planetesimal coagulation process.However the simplified expression for the drift rate, which always results in inwards-directed drift, likely artificially lowered the importance of the effect.
In this work we revisit the calculation of the radial drift of a planetesimal in an eccentric gaseous disk.In Section 2, we use numerical integration to calculate the drift rate exactly as a function of planetesimal eccentricity given a power law model for the radial variation of gas density and eccentricity.In Section 3, we calculate also the effect of gas drag on the planetesimal's eccentricity vector.Combining this with the secular equations for the evolution of planetesimal eccentricity, we identify two scenarios in which planetesimals undergo sustained outwards migration.In Section 4, we discuss the ramifications of these calculations, and some limitations and uncertainties.In Section 4.1, we calculate the evolution of the planetesimal surface density in one disk model, showing dramatic enhancement in specific locations.In Section 4.2, we quantitatively compare our results with the previously adopted approximation.In Sections 4.3 -4.5, we discuss some uncertainties and limitations of our model.In Section 5, we briefly summarize our conclusions.Appendix A develops a perturbative approximation that allows the main results to be understood analytically, but is not adequately accurate for quantitative calculations.

CALCULATION OF RADIAL DRIFT IN ECCENTRIC DISKS
In this section, we describe the equations used to numerically calculate the radial drift rate of planetesimals due to gas drag in an eccentric disk.We consider the case of planetesimals with Stokes number much greater than 1, so that they may be assumed to follow Keplerian orbits.We limit ourselves to a disk with no warp, and apsidally aligned streamlines, and we assume the planetesimal to orbit in the disk midplane.We take the eccentricity e d and surface density at pericenter Σ peri and temperature at pericenter T peri of the streamline with semi-major axis a to be power law functions of a: We take a power-law model for disk eccentricity to allow semi-analytic computations, but we note that the true eccentricity structure is likely more complicated.Simulations addressing the eccentricity of circumstellar disks in binary systems show a wide variety of behaviors.We further discuss the role played by the disk eccentricity profile in Section 4.3.We consider the motion of a planetesimal orbiting in the disk midplane, with semi-major axis a p , eccentricity e p , and longitude of pericenter ϖ p with respect to the apsidal line of the disk.Let θ denote the angle of the particle's radial vector with respect to the apsidal line of its orbit.A sketch of the disk with the relevant angles labeled is shown in Figure 1.To calculate the planetesimal's velocity, we begin with the expression for the shape of a Keplerian orbit: We use conservation of energy and angular momentum to calculate ṙp and θ.That is to say, the following equations in conjunction with Equation (2), result in these relations for the radial and azimuthal velocity: .
Let ϕ = θ + ϖ p be the angle of the planetesimal's radial vector with respect to the apsidal line of the disk.Then, given the assumption that the gas streamlines are Keplerian ellipses, we substitute Equation (1) for the gas eccentricity into Equation ( 2), to yield an equation for the semi-major axis a g of the gas streamline at position [r, ϕ]: where ã = a 0 /a g .Throughout this work, we use the subscript "g" to denote a property of the gas, and "p" to denote a property of the planetesimal.We then solve Equation ( 5) numerically for ã.
Once a g at the location of the planetesimal is known, the eccentricity e g is given by Equation (1).In the absence of a pressure gradient, the radial and azimuthal velocities of the gas are given in terms of a g and e g by Equations (4).We then assume that the velocity is reduced by a factor of 1 − 2η(a g ), independent of location along the streamline, where η is given by For an axisymmetric disk, this is equivalent to Equation 3.4 from Adachi et al. (1976), noting that their α is equivalent to our (p + 3/2 − s/2), their β is our s, and the mean thermal velocity c m is related to the isothermal sound speed c s by c m = c s 8/π.Here c s is the isothermal sound speed, calculated assuming a mean mass per particle of 2.36 Daltons and a gas temperature appropriate for r = a [see Equations ( 1), ( 10) and ( 11)].v K (a) = GM/a is the Keplerian speed of a circular orbit with radius a. G is Newton's constant, and M is the mass of the central star.This results in a gas velocity given by This treatment does not account for the azimuthal variation in velocity arising due to the pressure gradient, but as argued in Section 4.4, the effect of the azimuthal component of the pressure gradient on the gas velocity is small compared to that of the standard radial pressure gradient term present in an axisymmetric disk.As the main aim of this paper is to elucidate the impact of disk eccentricity on the radial drift, and as the effect of pressure on the gas streamlines is uncertain, we leave further discussion of pressure effects to future work.In order to calculate the drag force on the planetesimal, The disk is assumed to have the fiducial parameters listed in Table 1.
we need to know the gas density as a function of position.The surface density distribution of a disk described by Equation ( 1) is given by (see Statler 1999;Silsbee & Rafikov 2015) where E is the eccentric anomaly at azimuthal angle ϕ of the streamline with semi-major axis a g , and Σ(a g , ϕ) is the surface density at azimuthal angle ϕ along the streamline with semi-major axis a g .Σ peri (a g ) ≡ Σ(a g , 0) is the surface density at pericenter.Then, assuming assuming the disk to be in vertical hydrostatic equilibrium and vertically isothermal, the mid-plane density along the apsidal line of the disk is given by where H(r) = rc s /v K (r).
We assume an adiabatic equation of state on an orbital timescale, and that all parts of the disk are in vertical hydrostatic equilibrium.In Section 4.5, we explore the effect of other treatments of the disk thermodynamics and vertical relaxation timescale.
In order to calculate the variation of mid-plane density ρ g with E, we use the following equations.At eccentric anomaly E, we let C can be found easily from Equation ( 8), and D from the relation r = a(1 − e cos E).In order to solve for A, we need 4 other equations.These are: 1. Equation of State: Here γ is the adiabatic index.We take this to be 10/7, appropriate for a mixture with a 5 : 1 number density ratio of H 2 and He, assuming 5 accessible degrees of freedom per H 2 molecule.The adiabatic index of H 2 has some temperature dependence, due to the change in number of accessible degrees of freedom; see Figure 1 of Sharda et al. (2019).Solution of these equations leads to the relations At azimuthal angle ϕ, we take the mid-plane gas density ρ g given by Equations ( 9) and ( 11).
The instantaneous drag force on a planetesimal of radius R is given by where C d is a constant of order unity that we take to be 0.5 in our numerical calculations.
We then calculate the rate of change in semi-major axis using the Gauss planetary equation relating the force on a body to the change in semi-major axis; see Equation 1.200 from Tremaine (2023): where m p = 4πρ p R 3 /3 is the mass of the planetesimal.In our calculations, we assume the planetesimal density ρ p = 2 g cm −3 .Figure 2 shows the instantaneous value of ȧp as a function of θ for different values of the planetesimal eccentricity.We assume that the planetesimal is apsidally aligned with the disk (i.e.ϖ p = 0), and the disk parameters are those given in Table 1.In order to calculate the orbit-averaged inspiral rate, we average the result in Equation ( 13) over θ, weighted by the inverse angular velocity, to get ⟨ ȧp (a p , e p )⟩ = (1 − e 2 p ) 3/2 2π 2π 0 ȧp (a p , e p , θ) (1 + e p cos θ) 2 dθ.
(14)  The top left panel shows the results for a disk and planetesimal with the fiducial parameters in Table 1.In this case, the region of parameter space, corresponding to outwards migration is centered on the k p axis, at a value of k p roughly 1/2 the local disk eccentricity.
The cause of outwards migration can be understood qualitatively as follows.For an apsidally aligned particle with eccentricity less than the local disk eccentricity, at pericenter the particle is interacting with a streamline with larger semi-major axis.Therefore, ignoring the effect of pressure on gas velocity, i.e. setting η to zero in Equation ( 7), the gas velocity is higher than the particle velocity at pericenter, since larger semi-major axis corresponds to larger (less negative) specific energy.The opposite is true at apocenter.Therefore the particle gains energy at pericenter and loses it at apocenter.This behavior can be seen in the red curve in Figure 2. Note that this is the reverse of how it works in an axisymmetric disk, where the planetesimal gains energy at apocenter and loses it at pericenter.In an eccentric disk, the energy gain at pericenter can be larger than the corresponding loss as apocenter for two reasons.First the gas density is generally higher at pericenter than at apocenter.Second, even including a non-zero η, when the disk eccentricity gradient is sufficiently steep, the relative velocity at pericenter is larger than that at apocenter for moderate values of the eccentricity difference [see Equations ( A8) -( A10)].
Apsidally aligned planetesimals with eccentricity larger than the local disk eccentricity gain energy at apocenter and lose it at pericenter (see the blue curve of Figure 2).However, in this case the energy losses at pericenter are greater than the gains at apocenter, so the net drift is inwards.A particle with eccentricity exactly equal to the local disk eccentricity (Figure 2, green curve) experiences motion relative to the gas only because of the pressure gradient, and therefore loses energy over its entire orbit.
In the top right panel, the disk temperature is increased by a factor of two relative to the fiducial value, i.e.T 0 = 400 K.This has two effects.First, since surface density is held fixed, it decreases the mid-plane gas density, thus slightly slowing the radial drift.Second, it increases the value of η, [see Equation ( 6)], thus reducing the amount of area in the figure corresponding to outwards motion.
The middle left panel shows the case with p = 2, i.e. a steeper radial decrease of surface density.In the blue region of the plot, the local gas streamline at planetesimal pericenter has higher semi-major axis than the local gas streamline at planetesimal apocenter.Therefore, counterintuitively, a steeper radial decrease of surface density increases the ratio of the drag force at apocenter relative to that at pericenter.As discussed above, for orbits with eccentricity vector lying in the blue-shaded region, the planetesimal is gaining energy at pericenter, whereas at apocenter it is losing energy.Therefore, a steeper surface density profile decreases the area in which the planetesimals drift outwards.
The middle right panel shows the case in which q = 1, i.e. the magnitude of the disk eccentricity increases inwards, rather than outwards.In this case, along the k p axis, for values of k p less than the disk eccentricity, the inspiral rate is faster than the standard case with q = −1.This is because the headwind felt at apocenter increases faster with the quantity ϵ = e 0 − k p , than does the tailwind at pericenter [see Equations (A9) and (A10)].Along the k p axis, for values of k p greater than the disk eccentricity, the inspiral rate is slower than in the fiducial case, but still negative.Ignoring the effect of the pressure gradient in the disk, planetesimals with h p = 0 and k 0 > e 0 feel a headwind at pericenter, and a tailwind at apocenter for q = 1.While the tailwind at apocenter is stronger than the headwind at pericenter [see Equations (A9) and (A10)], the drift is directed inwards due to a combination of the higher gas density at pericenter and the effect of the pressure gradient term.
The bottom left panel shows a case in which disk streamlines with a given semi-major axis have a factor of two lower eccentricity.In this case, the effect of the pressure gradient on the radial velocity becomes more important than the disk eccentricity, and therefore the drift is always inwards.In the bottom right panel, we have increased the disk eccentricity, and the region of outward radial drift is correspondingly larger.

SUSTAINED OUTWARDS DRIFT
In the previous section we ignored the effect of gas drag on the eccentricity vector.When only gas drag is acting on the planetesimal, the planetesimal's eccentricity vector will become aligned with that of the local streamline, at which point drift is determined solely by the pressure gradient.In this case, in the absence of a local pressure maximum, outward migration will only be a transient phenomenon resulting from a non-equilibrium initial condition.However, as we show in this section, sustained outwards drift is possible in the case that the planetesimal eccentricity is perturbed by something in addition to gas drag.In this section, we calculate the evolution of the planetesimal eccentricity vector and semi-major axis under the combined influence of gas drag and secular gravitational perturbations from the disk and a binary companion star.
To calculate the evolution of k p and h p from the gas drag, we insert our expression for the drag force given by Equation ( 12) into the Gauss equations for the evolution of eccentricity and apsidal angle [see Equations 1.200 from Tremaine (2023) We then average these over θ as done in Equation ( 14).kp and ḣp are given by kp = ėp cos ϖ p − e p sin ϖ p πp ; ḣp = ėp sin ϖ p + e p cos ϖ p πp .
In addition to gas drag, we assume that the planetesimal is perturbed by an external companion star with mass ratio ν, which orbits with semi-major axis a b , and orbital eccentricity e b .The orbit of the companion star is assumed to be co-planar and apsidally aligned with the gas disk.We consider only the secular perturbations from the companion star on the planetesimal orbit.
where A is the precession frequency of the free eccentricity, given by the sum of terms corresponding to the disk and the binary and the eccentricity excitation terms B d and B b are given by n p is the mean motion of the particle.kdrag p and ḣdrag p are calculated numerically using Equations 16 -18 from the current paper, rather than the approximation in Rafikov & Silsbee (2015a).The coefficients ψ 1 and ψ 2 , which take into account the location of the disk edges, are given in Equations A.33 and A.34 of Silsbee & Rafikov (2015).
We identify two sets of parameters leading to a region of the disk in which sustained outwards drift is possible.One case is that of a low-mass disk apsidally aligned with the binary orbit, in which gravitational perturbations are dominated by the binary.The second set of parameters corresponds to a high-mass disk in which gravitational perturbations from the disk dominate over those from the binary companion.

Sustained outwards drift in a low-mass disk
Figure 4 shows the results of two simulations of the orbital evolution of planetesimals in a lowmass disk (Σ 0 = 300 g cm −2 ), subject to gravitational perturbations from a binary companion.The left panels show the evolution of the planetesimal's eccentricity vector.The right panels show the evolution of the semi-major axis.Time is given by the color scale (see color-bars on the right).The top panels correspond to a planetesimal with radius 10 km, and the bottom panels to one with radius 1 km.The time required to drift from 2 to 3 AU is about a factor of 10 higher for the 10 km planetesimal than for the 1 km one.
In the left panels, the black line shows the equilibrium eccentricity calculated with Equations 22 and 23 from Rafikov & Silsbee (2015a)1 .This is visibly, though not dramatically different from the colored curve in the case with R = 1 km, and lies very close to the colored line for R = 10 km.This shows that the approximate calculation in Rafikov & Silsbee (2015a,b) adequately describes the planetesimals' eccentricity evolution.We also calculated the true equilibrium eccentricity as a function of a p by evolving k p and h p to steady state while holding a p fixed.We find that for a p ≥ 2.02 AU, h p and k p given by the colored curve deviate from their equilibrium values by less than 10 −4 .In other words, after a brief transient behavior (shown by the spiral portion of the colored curve), at all times the eccentricity vector remains very near the equilibrium value at the current semi-major axis.This is a consequence of the fact that the timescale for evolution of the semi-major axis is much longer than the timescale for gas drag to equilibrate the eccentricity.

Sustained outwards drift in a high-mass disk
We consider also a case in which the disk is massive enough (for the chosen binary separation and mass ratio) to have a dominant effect on the gravitational perturbations acting on the planetesimals.In this case, for our fiducial values of p and q, the equilibrium planetesimal eccentricity is larger than the disk eccentricity, and therefore does not support sustained outwards drift.
However, there are values of p and q for which sustained outwards drift is possible.Figure 5 shows the drift rate at 1 AU for a system with the fiducial parameters, except Σ 0 = 10,000 g cm −2 , and p  and q vary as shown on the axes.In this case we kept the inner edge of the disk at a in = 0.1 AU, but placed the outer edge of the disk at a out = 2.5 AU.This avoids streamline crossing, which would occur for models with q < −1 if we allowed the disk to extend to 5 AU.The maximum value of the (outward-directed) radial drift rate is labeled on each panel.In this case we see that substantial drift within the disk lifetime is possible even for ∼ 100 km planetesimals.Since we assume that gravitational perturbations from the disk dominate over those from the binary, both the perturbing gravitational force and the force of gas drag are directly proportional to the disk mass.As can be shown rigorously from the equations in Rafikov & Silsbee (2015a), the equilibrium planetesimal eccentricity becomes independent of disk mass in this case.Therefore, the inspiral rate for planetesimals of a given size is directly proportional to the disk mass, since the drag force for a planetesimal on a particular orbit is directly proportional to the gas density.Above a certain size, the planetesimal eccentricity vectors become aligned (Rafikov & Silsbee 2015a), and therefore the drift rate is inversely proportional to planetesimal size.Thus, the peak drift rate for 300 km planetesimals is almost exactly a factor of 10 slower than that for 30 km planetesimals.Radial drift rate at 1 AU assuming the planetesimal's eccentricity vector to be at its equilibrium value.This is in a system with the fiducial parameters in Table 1, except that Σ 0 = 10,000 g cm −2 , a out = 2.5 AU, and p and q are labeled on the axes.The planetesimal size is labeled on each panel, and the maximum outwards drift rate is labeled on each panel.The color scale does not extend beyond 10 −4 , which leads to a saturation of the red color in the top two panels.

Impact on planet formation
In the previous section, we pointed out two scenarios in which outward migration of planetesimals could occur.Since the radial pressure gradient yields a bias for inwards drift, it is difficult to construct a scenario in which planetesimal drift is outwards throughout the entire disk.Thus outward migration, if present, will be limited to a specific part of the disk.For a given planetesimal size, we denote the boundaries of such a region a 1 and a 2 , with a 2 > a 1 .The point a 2 is a stable attractor, in the sense that all planetesimals in the neighborhood of a 2 drift towards a 2 .In contrast, a 1 is an unstable stationary point in the sense that in the neighborhood around a 1 , all planetesimals (except those with a exactly equal a 1 ) drift away from a 1 .Over time, planetesimals become concentrated around a 2 , and if the location of a 2 is almost independent of planetesimal size over a wide range of sizes, the overall surface density of planetesimals may become strongly enhanced at a 2 .The location of the region of outward migration (and therefore the location of the points a 1 and a 2 ) depends on the choice of disk parameters.We discuss some of these dependencies in Section 4.3, however because of the complicated dependence of the disk gravity on parameters, and because the sign of the drift is determined by multiple effects acting in different directions, we are not able to give, even in the context of the current model a concise analytic expression for the location of a 1 and a 2 .
Possible concentration of planetesimals due to radial drift was discussed briefly in Rafikov & Silsbee (2015b), however it was not found to be a dramatic effect.The authors of that work used a simplified model for the drift, taking Equation 4.21 from Adachi et al. (1976) (derived and intended to be used for an axisymmetric disk) and replacing the planetesimal eccentricity with the relative planetesimal gas eccentricity.Unlike the current treatment, this results in a drift rate that is always inwards, limited from below by the rate arising from the pressure gradient.Roughly speaking in such a model, the peak planetesimal density enhancement in a local region of slower drift is equal to the ratio of the drift rate in the surrounding region to that in the local region (assuming that the region of slow drift is small enough that the planetesimal flux is approximately constant across it).In contrast, in the current model which allows for convergent planetesimal migration, given sufficient time, the density contrast can grow to be arbitrarily large.
We simulated the evolution of the planetesimal surface density under the action of radial drift.Motivated by the results in Figure 4, we consider a system in which p = 0.3, q = −1.4,Σ 0 = 10,000 g cm −2 , a in = 0.1 AU, and a out = 2.5 AU, and the rest of the parameters take their fiducial values.We began with a population of planetesimals distributed smoothly throughout the disk with semi-major axes from 0.2 to 2 AU assuming a local planetesimal to gas mass ratio of 0.2%.We calculated the drift rate assuming planetesimals to orbit at their equilibrium eccentricity, and then used that drift rate to solve for the evolution of the surface density distribution.
Figure 6 shows the planetesimal surface density at pericenter as a function of semi-major axis for different times between 0 and 10 6 years.The top panel shows the case in which all planetesimals are one km in size.These planetesimals experience outwards drift between a 1 = 0.77 AU and a 2 = 1.19 AU (labeled on the panel).This results in a few noticeable features on the graph.First, a singlepixel spike develops at a 1 after 1000 years.This is related to the discretization of semi-major axis space, and the fact that the drift rate drops to zero at a 1 .Since this spike is just composed of the planetesimals that were initially at that location, it is no higher than the initial surface density at that location.
In the vicinity of a 2 , we see two peaks growing and drifting towards a 2 .One peak just outside of a 2 is composed of all the particles in the outer disk that have drifted inwards, and one just inside a 2 , is composed of planetesimals initially between a 1 and a 2 .By 1 Myr, they have merged together forming a single peak with width equal to the smallest resolvable width in the simulation, and height over 1000 times the original surface density at that point.In principle if we considered infinitely narrow bins in semi-major axis space and let planetesimals drift indefinitely in a static disk with no self-interaction, the two spikes would become arbitrarily high but never merge.However at high planetesimal density and small scales, other physics that we don't include (e.g.collisions and planetesimal-planetesimal scattering) play a role as well, thus altering the size distribution and smoothing the surface density distribution.We include this plot to demonstrate that a strong enhancement of the surface density may occur, not to claim an accurate calculation of the profile of the surface density peak.
The middle panel shows the case in which all planetesimals have a size of 100 km.These experience outward drift between a 1 = 0.45 AU and a 2 = 1.19 AU.The same features are visible as in the top panel, but because the larger planetesimals drift more slowly, the surface density evolves more slowly and the two peaks around a 2 are still clearly distinguishable even after 10 6 years.The complicated shape of the peaks reflects the behavior of the drift velocity in the vicinity of a 2 .
The bottom panel shows the surface density at pericenter when we consider a size distribution from 1 to 100 km, with the initial number density of planetesimals proportional to R −7/2 .Since a 1 varies significantly between different planetesimal sizes, we see multiple discrete spikes corresponding to individual planetesimal sizes.For the system parameters considered here, planetesimals of all sizes between 1 and 100 km have almost the same value of a 2 , so the peaks around at a 2 remain quite sharp.The surface density of planetesimals in this simulation reaches a peak of 180 times its initial value, but again, we caution that diffusive processes not included in this simulation likely smooth the density profile somewhat.The fact that a 2 is nearly the same for a wide range of planetesimal sizes is a consequence of the fact that planetesimals above a certain size orbit at nearly the forced eccentricity determined by the gravitational potential of the disk and binary (Rafikov & Silsbee 2015a).It is therefore expected that in many cases the pile-up location will be nearly size-independent.

Comparison with an axisymmetric disk
While this paper has focused primarily on the implications of outwards radial drift, it is worth asking in general how much the drift rate differs from a simple extrapolation from the case of an axisymmetric disk with an eccentric planetesimal.
We made this comparison for the disk models shown in Figure 3. Figure 7 shows the ratio of the inspiral rate calculated using Equations ( 13) and ( 14), to that calculated using Equation 4.21 of Adachi et al. (1976) (intended to be used for an axisymmetric disk), with the eccentricity replaced with the relative eccentricity between the planetesimal and the gas streamline with the same semimajor axis, as done in Rafikov & Silsbee (2015b).In performing this approximate calculation, for a planetesimal of semi-major axis a p , we approximated the gas density as equal to the gas density at pericenter for the streamline with semi-major axis a p .
Areas of major disagreement (i.e.areas in which the approximate treatment gets the sign wrong, or overestimates the rate by a factor of 10 or more) are present in all panels but one, but make up less than half of the area of the panel, except for the bottom right panel, in which such areas comprise almost exactly 50% of the panel. .Evolution of the planetesimal surface density under the effect of radial drift: each curve shows the surface density of planetesimals along the apsidal line of the disk as a function of semi-major axis at a different time since the beginning of the simulation, as labeled in the legend.The disk and binary properties are identical to that of the "high mass" case discussed in Section 3.2.The top panel corresponds to 1 km planetesimals, the middle panel to 100 km planetesimals, and the bottom panel to a size distribution between 1 and 100 km proportional to R −7/2 .In the top and middle panels we have labeled the values a 1 and a 2 corresponding to the endpoints of the range of outward migration for planetesimals of that size.
Looking at the bottom panels, we note that the amount of disagreement depends strongly on the disk eccentricity.This makes sense, because the two methods agree as to the effect of the pressure gradient.For low eccentricity, the inspiral rate is dominated by the pressure gradient, and therefore the difference between the approximate and exact treatments is not large (see bottom left panel).The most dramatic differences occur for the disk in which e 0 = 0.15 (bottom right panel), as in this case the eccentricity, rather than the pressure gradient, has the dominant effect on the particle-gas velocity.

Effect of disk parameters
The radial drift rate is affected most strongly by the eccentricity profile of the disk.For outward drift to be possible, depending on the temperature, the eccentricity must be increasing outwards and ≳ 0.1.Simulations of disks in binary systems often find a region of the disk with these properties (e.g Marzari et al. 2009;Martin et al. 2020;Jordan et al. 2021), however the details vary substantially from one work to another.Jordan et al. (2021) ran a series of 2D radiative hydrodynamic simulations using using the FARGO code.They found that the eccentricity profile was dramatically affected by choices regarding the simulation domain, the boundary condition at the inner edge of the disk, the treatment of thermodynamics, and numerical resolution.Marzari et al. (2009) found that inclusion of self-gravity substantially decreased the disk eccentricity, however even in their self-gravitating model with a perturber orbiting at 30 AU, the disk reached a maximum eccentricity of 0.25 at 1 AU; higher than in any of the models considered in this paper.
It is also of note that disks in many simulations precess.For our high disk-mass case, precession, provided it is not too rapid, should not affect our results.However, the case of a low mass-disk in which binary gravity is important (illustrated in Figure 4), outward radial drift is only possible in an aligned disk.Rough apsidal alignment (though with substantial libration and some systematic offset) has been found in the simulations of Müller & Kley (2012); Gyergyovits et al. (2014); Martin et al. (2020), however other simulations (Paardekooper et al. 2008;Marzari et al. 2009Marzari et al. , 2012) ) find anti-alignment or precession.
In general, we find that the radial drift of planetesimals in disks with eccentricity ≥ 0.1 may show strong deviations from the values calculated using the approximate method in Rafikov & Silsbee (2015b).If these moderate eccentricities are combined with an outward-increasing eccentricity gradients in part of the disk, then there may be a local region in which the radial drift is outwards.Although the true surface density and eccentricity profiles are unlikely to be well-represented by power laws, it seems plausible that for more realistic disk models, these conditions would be satisfied somewhere within the disk.
The disk temperature is also important, since increasing the disk temperature by a given factor increases the pressure gradient and therefore η by the same factor [see Equation ( 6)].We use a disk temperature of 200 Kelvin at 1 AU, as appropriate for a disk heated by the radiation of the central star (Chiang & Goldreich 1997).Models show that a massive disk undergoing viscous accretion will become several times hotter than this (D'Alessio et al. 1998;Dullemond et al. 2007).However, the relation between accretion rate and mid-plane temperature is not simple.For example, some models suggest that accretion occurs mainly near the disk surface, with a low-viscosity "dead" zone at the disk midplane.In this case, the midplane temperature may be substantially cooler at a given accretion rate than in a simple model with vertically uniform viscosity (Hirose & Turner 2011).More dramatically, Mori et al. (2019) found that a magnetic wind driven accretion model can achieve high   13) and ( 14) to that given by Equation 4.21 in Adachi et al. (1976) (derived and intended by the authors to be used for an axisymmetric disk), with planetesimal eccentricity replaced by relative planetesimal-gas eccentricity as done in Rafikov & Silsbee (2015b).The ratio is shown by the color bar for values between -2 and 2, and with labeled contours for higher values.We include also labeled contour lines at 0 and 0.5.accretion rates with very low rates of mid-plane heating.For example, they considered a disk with surface density 17,000 g cm −2 at 1 AU, with accretion rate of 8.2 × 10 −7 M ⊙ /yr.This resulted in a midplane temperature slightly under 200K, compared with ∼ 900K in the viscous heating model with the equivalent accretion rate.For our fiducial disk, the temperature must be increased to 600 K at 1 AU for there to be no value of the planetesimal eccentricity vector leading to outwards migration.If we increase e 0 to 0.15 (as in the bottom right panel of Figure 3), then this critical temperature increases to 1470 K.

Effect of azimuthal pressure gradients
Here, we justify the use of an azimuthally independent value for η in Equation ( 7) by showing that pressure gradients along the streamline have a small effect on the velocity.Starting from Equation (11), we use Bernoulli's law to calculate the change in velocity along a streamline: .
(22) To isolate the effect of the azimuthal pressure gradient, we Taylor-expand the velocity given by Equation ( 22 where ) is the velocity in the absence of a pressure perturbation along the streamline, and is the contribution due to the pressure gradient.Using Equations ( 8) and ( 11), and dropping all terms of order e 2 , we find that the change in velocity between pericenter and apocenter, attributable to the pressure gradient along the streamline is ∆v This is less by a factor of 4γ(3 − 2q)e g /[(γ + 1)(p + (s + 3)/2] with respect to the typical "headwind" velocity due to the radial pressure gradient in a circular disk [see Equation ( 6)].At location a 0 in our fiducial disk, this factor is 0.43 -not a terribly small number.However this analysis would suggest that by calculating η at a g as we have done in Equation ( 7), we are underestimating the gas velocity at apocenter, and overestimating it at pericenter, and thus these errors should partially cancel each other.At any rate, it was shown in Figure 3 that our results are robust to a factor 2 increase in disk temperature above the fiducial value (corresponding to a factor 2 change in η), suggesting that this additional uncertainty will not alter any of the conclusions of the paper.

Different treatments of disk thermodynamics
Throughout this work, we assume that the gas in the disk behaves adiabatically on an orbital timescale, and that the disk height at a given azimuthal angle is given by c s r/v k (r).Here we explore three other limiting cases, and show that, except for the rather unrealistic assumption of an azimuthally invariant disk height along a given streamline, the effects on our results are minimal.
First, we consider an isothermal model in which the gas temperature along a streamline does not vary with azimuthal angle.In this case, χ = 1 [see Equation ( 10)], and re-derivation of Equation ( 11) with this assumption yields Second, we consider a case in which the disk height is constant along a streamline.In this case, the mid-plane density is directly proportional to the surface density, so This scenario is not entirely implausible as the rough timescale for the disk to come into vertical equilibrium (H/c s ) is of the same order as the dynamical timescale r/v K .However, as the time between apocenter and pericenter is π times the dynamical timescale, it is likely a better approximation, as we have assumed in the rest of this paper, that the vertical structure of the disk is in local hydrostatic equilibrium.
Finally, we consider a case in which the disk material is assumed to rapidly equilibrate its temperature in less than an orbital time.In this case, the temperature anywhere along a streamline is given by Equation ( 1) with a replaced by r, resulting in the relation Figure 8 shows the results of these different limiting cases for the disk thermodynamics.The upper left panel shows our fiducial case of an adiabatic disk in local vertical hydrostatic equilibrium (and is therefore identical to the upper left panel of Figure 3).The upper right panel shows the case where we assume a constant temperature along a streamline [Equation ( 26)].Compared with the adiabatic case, this results in a disk in which the variation of density along a streamline between pericenter and apocenter is higher, as the disk is warmer at apocenter.Since in the disk models under consideration, an apsidally aligned planetesimal gains energy at pericenter and loses energy at apocenter, a disk with isothermal streamlines results in a slightly enlarged region of outward migration.
In contrast, assuming a constant disk height along the streamline [Equation ( 27)], significantly increases the ratio of the density at apocenter to that at pericenter compared with the fiducial case, and thus shrinks the region of outward migration, (see bottom-left panel of Figure 8).Finally assuming that material instantaneously equilibrates its temperature, and that the temperature is a sole function of r [Equation (28)], results in a somewhat cooler temperature at apocenter, similar to the adiabatic case.Except for the somewhat unrealistic case where we assume that disk height doesn't change along a given streamline the changes associated with assumptions about thermodynamics are minor.

CONCLUSIONS
In this paper, we calculated the rate at which planetesimals drift radially in eccentric disks under the influence of aerodynamic drag.We found in some cases that the drift rates differ by an order of magnitude from those calculated for planetesimals in circular disks with the same eccentricity relative to the gas.This is true even for relatively small (∼0.1) disk eccentricities.In contrast both with typical expectations in protoplanetary disks, as well as results of previous approximate calculations, we found that under some plausible assumptions about the disk eccentricity profile, it is possible that planetesimals may drift away from the central star.
In addition to calculating radial drift as a function of planetesimal eccentricity, we used secular perturbation theory to calculate the equilibrium eccentricity as a function of planetesimal size for planetesimals evolving under the combined influence of gas drag, and gravitational perturbations from the companion star and the disk.We verified that, as in the case of an axisymmetric disk, the eccentricity vector evolves much faster under the influence of gas drag than does the semi-major axis.As a result, as planetesimals drift in the disk, they remain very close to their equilibrium eccentricity for their present semi-major axis.We found two classes of model in which planetesimals orbiting with the equilibrium eccentricity experience outwards drift over a substantial range of orbital radii, suggesting that in circumstellar disks with a close companion star there will plausibly be regions in which the planetesimal migration is outwards.
The first such case is a low surface density disk apsidally aligned with the binary orbit.In such a scenario, we found outward migration times of hundreds of thousands of years for a 1 km planetesimal, at a location where the disk eccentricity is 0.1.Second, we studied planetesimals in a high-mass disk (∼6 times MMSN at 1 AU) at a location with the same value of the disk eccentricity.Here, for favorable surface density and eccentricity profiles, we found outward migration timescales under 100,000 years, even for 30 km planetesimals.In this case, the disk need not be aligned with the binary, and the disk mass itself is not critical, provided that it is large enough to dominate the gravitational perturbations of the binary at the location of the planetesimal.
The outer edge of the zone of outward migration is a stable attractor for planetesimals, and thus over time the surface density of planetesimals is enhanced near this location.Sufficiently large planetesimals have their eccentricity vectors driven to the so-called "forced eccentricity" determined by the gravitational perturbations.Since the eccentricity vector is independent of planetesimal size, then the region of the disk in which planetesimals migrate outwards is also independent of planetesimal size.Thus all planetesimals pile up at the same location.In our model, this resulted in an enhancement by more than a factor of 100 in the planetesimal surface density at the edge of the zone of outward migration, though diffusive processes not accounted for in our modeling likely lower this number. .This shows the inspiral rate, for planetesimals apsidally aligned with the gas streamlines (h p = 0, k p ≥ 0).The red curves are made using the exact solution given by the method described in Section 2. The blue curves are from the analytic approximation in Equation (A18).Each panel corresponds to the same set of disk parameters as in the corresponding panel of Figure 3.The green star corresponds to k p = e 0 , or ϵ = 0.

Figure 1 .
Figure 1.Sketch of a planetesimal embedded in an eccentric disk with the angles used in this work labeled.

Figure 2 .
Figure2.Dependence of the instantaneous value of ȧp on angle θ for a planetesimal apsidally aligned with the disk (i.e. with ϖ p = 0. Different color curves correspond to different values of the planetesimal eccentricity.The disk is assumed to have the fiducial parameters listed in Table1.

Figure 3 Figure 3 .
Figure3shows the inspiral rate of a 1-km planetesimal in a variety of model disks in the space of the two components k p and h p of the eccentricity vector of a planetesimal embedded in the disk.These are defined in terms of e p and ϖ p as

Figure 4 .
Figure 4. Evolution of planetesimal eccentricity vector (left panels) and semi-major axis (right panels) under the combined influence of gas drag, and secular gravitational perturbations from the disk and binary companion.The color corresponds to the amount of time since the beginning of the simulation, as shown in the color-bars on the right.The top panels show the trajectory of a 10 km planetesimal, and the bottom a 1 km planetesimal.This assumes the fiducial system parameters [see Table 1], except Σ 0 = 300 g cm −2 , and e 0 = 0.05.The black line (almost overlapping the colored line in the upper panel) shows the equilibrium eccentricity calculated using Equations 22 and 23 from Rafikov & Silsbee (2015a).
Figure5.Radial drift rate at 1 AU assuming the planetesimal's eccentricity vector to be at its equilibrium value.This is in a system with the fiducial parameters in Table1, except that Σ 0 = 10,000 g cm −2 , a out = 2.5 AU, and p and q are labeled on the axes.The planetesimal size is labeled on each panel, and the maximum outwards drift rate is labeled on each panel.The color scale does not extend beyond 10 −4 , which leads to a saturation of the red color in the top two panels.
Figure6.Evolution of the planetesimal surface density under the effect of radial drift: each curve shows the surface density of planetesimals along the apsidal line of the disk as a function of semi-major axis at a different time since the beginning of the simulation, as labeled in the legend.The disk and binary properties are identical to that of the "high mass" case discussed in Section 3.2.The top panel corresponds to 1 km planetesimals, the middle panel to 100 km planetesimals, and the bottom panel to a size distribution between 1 and 100 km proportional to R −7/2 .In the top and middle panels we have labeled the values a 1 and a 2 corresponding to the endpoints of the range of outward migration for planetesimals of that size.

Figure 7 .
Figure7.Ratio of the inspiral rate given by Equations (13) and (14) to that given by Equation4.21inAdachi et al. (1976) (derived and intended by the authors to be used for an axisymmetric disk), with planetesimal eccentricity replaced by relative planetesimal-gas eccentricity as done inRafikov & Silsbee (2015b).The ratio is shown by the color bar for values between -2 and 2, and with labeled contours for higher values.We include also labeled contour lines at 0 and 0.5.

Figure 8 .
Figure8.The inspiral rate for a one kilometer planetesimal at 1 AU in a disk with our fiducial parameters, in the space of the components k p and h p of the planetesimal eccentricity vector.Different panels correspond to different limiting cases of the disk thermodynamics as discussed in the text.As in Figure3, the color shows the rate of radial drift, with red indicating drift towards the star, and blue indicating drift away from the star.The peak (outwards) drift rate is labeled on each panel.The color bar saturates at ±10 −4 AU/yr.The green star in each panel represents the gas eccentricity.
Figure9.This shows the inspiral rate, for planetesimals apsidally aligned with the gas streamlines (h p = 0, k p ≥ 0).The red curves are made using the exact solution given by the method described in Section 2. The blue curves are from the analytic approximation in Equation (A18).Each panel corresponds to the same set of disk parameters as in the corresponding panel of Figure3.The green star corresponds to k p = e 0 , or ϵ = 0.