A dense ring of the trans-Neptunian object Quaoar outside its Roche Limit

Planetary rings are observed not only

The 2019 June 5 occultation provided the ring profiles in four bands and our data do not show any significant dependence of E p and A τ with wavelength.This can be paralleled with multi-wavelength occultations by the core of Saturn's F ring, implying that this core contains a substantial population of particles larger than around 10 µm (refs. 12,13,14 )Moreover, our numerical integrations show that isolated particles smaller than 100 µm should escape the ring over a few years' time scale due to radiation pressure, see Methods.
Contrary to Chariklo's ring 11 , Quaoar's ring is strongly irregular in azimuth.As such, it is reminiscent of Saturn's F ring that contains azimuthal features (clumps) or even local opaque structures interpreted as km-sized moonlets 15 .This clumpy nature is thought to be caused by the presence of thousands of small parent bodies (1-to 0.1-km in size) that collide and produce dense strands of micrometer-to centimeter-sized particles that re-accrete over a few months onto the parent bodies in a steady state regime 16,17,18 .A similar process may explain the patchy structure of Quaoar's ring.Another comparison can be made with respect to Neptune's Adams ring which also shows substantial longitudinal variations in optical depth.
In summary, our observations are consistent with a dense, irregular Quaoar's ring.The term 'dense' means that collisions play a key role in its dynamics.However, in contrast to all other known dense rings, Quaoar's ring orbits well outside the classical Roche limit.This excludes very tenuous dusty rings with τ < 10 −4 seen for instance outside Saturn's Roche limit 19 .These tenuous rings are dominated by radiation effects that permanently remove them from their place of formation.
Objects orbiting at distance a from Quaoar with mass M Q should have a density greater than the Roche critical density in order to avoid being disrupted by tidal forces.An often quoted value of γ is 1.6 which would correspond to a lemon-shaped particle aggregate filling its Roche lobe 20 .Using equation (1) and the values listed in Table 1, we obtain ρ Roche ∼ 30 kg m −3 , corresponding to extremely porous or even fluffy material.
The dense rings of the giant planets 1 , as well as Chariklo 2 and Haumea 3 lie within the Roche limit of the central bodies, assuming ρ=400 kg m −3 , typical of the small inner saturnian satellites 21 .With this value, Quaoar's classical Roche limit is near 1,780 km, much smaller than the ring radius of 4,100 km (Table 1 and Extended Data Fig. 2).
Numerical simulations of impact-generated disks show that material outside the Roche limit is expected to accrete over time scales of a few decades only 22,23 .This rapid processes would imply a very recent ring and, thus, a very low probability of being observed at the current epoch.That leaves us with a very young or extremely low density ring particles, or more likely, with the need for revisiting the Roche limit notion.In fact, equation (1) applies to a fluid satellite that is disrupted close to a planet.However, the reverse process -the accretion of colliding particles into a satellite -implies mechanisms not accounted for in equation (1).
Gravitational accretion of ring particles at a given distance depends not only on their bulk density ρ, but also on the radial velocity dispersion c between the particles.If c ≫ v esc , where v esc is the two-body escape velocity of the two ring particles involved in a pairwise impact, these impacts may avoid accretion outside the classical Roche limit One way to increase the radial velocity dispersion is the ring material to be perturbed by external forces, such as resonances with Quaoar itself, Weywot, or undiscovered satellites.Another, is if collisions between pairs are sufficiently elastic so that their post-impact Hill-energies are positive 24 .In this scenario, the evolution of c is governed by a competition between collisional dissipation and viscous gain of energy from the orbital motion, leading to a steady-state velocity dispersion c st .Its value depends on the relation ϵ(v n ) between the collisional coefficient of restitution ϵ and the perpendicular impact velocity between particles, v n .In particular, steeper drops of ϵ vs. v n result in smaller values of c st (Fig. 3 and Extended Data Fig. 4).
Laboratory measurements show that for icy particles, ϵ depends sensitively on their surface properties.Figure 3 displays simulations of particle accretion in Quaoar's ring with τ =0.25, with ϵ(v n ) models covering the range of existing laboratory measurements.Initial measurements 25 with frost-covered ice at temperature T =210 K (referred to as Model 1) suggest that ϵ drops rapidly with v n , attaining values below 0.5 for v n ≈ 1 mm s −1 , relevant to Quaoar's ring (Fig. 3).In this case, c st is of the order of v esc and accretion sets in rapidly in the self-gravitating ring simulations as soon as ρ exceeds ρ Roche , see Methods.In terms of the non-dimensional r H parameter, defined in equation (8), this corresponds to r H ≈ 1.1-1.2, in agreement with simulations of accretion threshold in Saturn's rings 26 .
However, better-controlled experiments 8 (referred to as Model 4) for particles with compactedfrost surfaces and lower temperatures 123 K indicated a much shallower ϵ(v n ) relationship.In this case the onset of accretion in Quaoar's ring requires ρ cr larger than about 5,000 kg m −3 , corresponding to r H > 5.This leaves a safe margin to prevent accretion of icy ring particles with ρ=400 kg m −3 .The reason for this dramatically different behavior is the much higher c st maintained by less dissipative impacts.
The accretion limits found above apply to optical depth τ ≲ 0.25.At larger τ the steady state value c st decreases because the reduced mean free path between impacts leads to a less effective viscous gain (Extended Data Fig. 5, panel a).For example, for Model 4, the value of c st decreases by a factor of nearly three as τ increases from 0.25 to 1, and the minimum value of ρ cr required for accretion drops to about 300 kg m −3 .This is small, but still in line with the expected density of icy particles.Finally, the critical density ρ cr also depends on the particle size R, ρ cr increasing as R decreases, making a ring with smaller particles less prone to accretion, see Methods, (Extended Data Fig. 5 panel b).In summary, there are no objections for having a τ ≲ 0.25 non-accreting ring far outside Quaoar's classical Roche limit, provided that particles have compacted-frost surfaces described in ref. 8 (Model 4).
However, two problems remain concerning the presence of Quaoar's ring.One is that external perturbations may lead to local condensations.Enhancing τ reduces c st and makes the ring prone to accretion.For example, Model 4 with ice particle densities of 400-900 kg m −3 would lead to this type of transition, as accretion is taking place for τ ≳ 1 but not at smaller τ .The other problem is that impacts cause viscous spreading, which needs to be balanced by some mechanism to maintain the ring narrow, such as resonances.
Quaoar's ring is in fact close to both the outer Quaoar 1/3 Spin-Orbit Resonance (SOR) and the inner Weywot 6/1 Mean Motion Resonance (MMR), see Methods and Extended Data Fig. 2. As a small body, Quaoar is expected to have an irregular shape that creates significant SORs stemming from the non-axisymmetric terms of its potential 9 .A similar configuration is observed around Chariklo's and Haumea's rings (Methods), suggesting that the 1/3 SOR may play a key role in the ring confinement.
An assessment of the effect of the 1/3 SOR can be achieved through N-body collisional simulations with a ring that completely surrounds Quaoar.Such simulations applied to Chariklo actually show that the 1/3 SOR not only excites the orbital eccentricities of the particles as expected from equation (11), but also leads to a ring confinement near this resonance 27,28 .However, these simulations do not account yet for self-gravity.These results are encouraging but they require further simulations using more realistic parameters, as well as theoretical models to understand the confinement mechanism.If confirmed, the 1/3 SOR would play two roles.One is to confine the ring and the other one is to maintain a velocity dispersion between parent bodies or particle aggregates high enough to prevent further accretion.
The 6/1 Weywot MMR is a more complex resonance than the 1/3 Quaoar SOR because it involves one corotation-type resonance and five eccentricity-type resonances, see Methods.Their strengths and couplings depend on poorly constrained parameters, that is Weywot's mass, orbital eccentricity and apsidal precession rate.Our calculations show that the 6/1 corotation resonance creates one L 4 -type Lagrange point that might concentrate Quaoar's ring material in a finite interval of longitudes if Weywot's orbital eccentricity e ′ is larger than about 0.1.At this stage of knowledge, neither of the two resonances is found to dominate the other.More observations are now necessary to pin down Quaoar's shape and Weywot's orbital elements, so that to better constrain the locations and strengths of the resonances that may interact with Quaoar's ring.Also, new detections of the dense region of the ring can be used to track its motion over time as another way to test models of resonant confinement.
Quaoar's ring is the third example of a dense ring around a small body found in the Solar System, suggesting that more still await discovery.Meanwhile, the large distance of this ring from Quaoar's means that the classical notion that dense rings survive only inside the Roche limit of a planetary body must be revised.• The opportunity to review the results and comment on the manuscript were given to all authors.The physical and orbital parameters of Quaoar and Weywot are obtained from the method described in ref. 7 .Details in Methods.
b From Kepler's third law a = (GM Q /n 2 ) 1/3 , assuming a spherical Quaoar of mass M Q , where G is the gravitational constant and n is the mean motion.c For each ring's pole orientation, an opposite direction is possible depending on the particles' direction of motion.Here, solution 2 is preferred, because it is closer to Weywot's orbital pole orientation.
d The uncertainty of the inclination between Weywot orbital pole (RA = 266 ± 10 degrees and Dec. +50 ± 6 degrees) and the ring pole is mostly dominated by Weywot.e In each station, "(bef.)"stands for before the closest approach and "(aft.)"for after it.
f The light curve S/N does not allowed the detection of the ring.g These values are not available as the ring is unresolved, only the Equivalent Width (Radial Width × Normal opacity) can be computed.
h Due the small number of points within each detections, the 1-sigma uncertainty has a large asymmetry.1) showing the steady-state velocity dispersion c st versus the bulk density ρ of the ring particles with assumed radius of 1 m.We used four velocity-dependent ϵ(v n ) models for the coefficient of restitution in impacts, spanning the range of various laboratory experiments of collisions between icy particles.Model 1 stands for Frost-covered ice spheres at temperature 234 , with v c = 0.0077 cm s −1 (ref. 25); Model 2 for Frost-covered ice at T = 123 K, ϵ(v n ) = 0.48v n −0.20 (ref. 8); Model 3 for as Model 1 but with v c = 0.077 cm s −1 ; and Model 4 for particles of radius R = 20 cm with compacted frost at T = 123 K, ϵ(v n ) = 0.90 exp(−0.22vn ) + 0.01v n −0.6 (ref. 8), here v n is the normal component of impact velocity (in cm/s).Also shown are simulations performed with a constant ϵ = 0.5.The results are displayed both in terms of ρ and the dimensionless Hill parameter r H .All simulations have optical depth τ = 0.25, and follow a co-moving ring region with a size of 160 m × 160 m.Open circles along the lines denote simulations which did not lead to particle aggregation.Filled symbols stand for simulations where particles accreted to gravity-bound aggregates.The boxes show a snapshot of the co-moving ring region with and without the accretion of the particles for Model 1 and Model 4. The shaded region sketches the region where accretion takes place (r H ≳ 1.2, c/v esc < 1.65 + 0.16 ln (ρ/80 kg m −3 ) ∼ 2), while the dashed line indicates the escape velocity at the surface of 1-m particles for a given bulk density ρ.

Methods 1 Prediction and observations
The four stellar occultations presented here were predicted using the standard procedures of the ERC Lucky Star project (https://lesia.obspm.fr/lucky-star)and are publicly available.The stars positions were obtained using Gaia Early Data Release 3 form 32 .Quaoar's ephemerides were derived from the Numerical Integration of the Motion of an Asteroid (nima 33 ) integrator, taking advantage of Quaoar's accurate positions derived from previous occultations.For Weywot's ephemeris, we use the genoid algorithm (GENetic Orbit IDentification 7,34 ) to fit an orbital model on the ten known observations of Weywot between 2006 to 2011 from the Keck observatory 6 and from a stellar occultation by Weywot, observed on 04 August 2019.It is a genetics-based algorithm that relies on a meta-heuristic method to find the most appropriate (i.e.minimum χ 2 ) set of dynamic parameters.On these data, the best results are from a Keplerian model of Weywot's motion around Quaoar.
Extended Data Table 1 provides astrometric and photometric information on the occultations and Supplementary Table 1 lists the circumstances of the observations.Ground-based stations and one space-borne instrument (the ESA/CHEOPS space telescope 31 ) were mobilized.

Data Analysis
The images of the occulted stars (plus the occulting object) were analyzed applying standard aperture photometry procedures 35 , using nearby reference stars to correct for sky transparency fluctuations.Secondary events (i.e., not caused by Quaoar itself) were observed on 02 September 2018, 05 June 2019, 11 June 2020 and 27 August 2021, revealing the presence of semi-transparent ring around Quaoar.
Abrupt opaque edges models, including the effects of diffraction, finite band width, exposure time and the Gaia DR2 stellar diameter, were fitted to the star dis-and re-appearances behind Quaoar's main body using the procedures within the Stellar Occultation Reduction and Analysis (sora) package 36 .For the ring events, an additional parameter is the apparent ring opacity p ′ that measures the fractional drop of stellar flux.This parameter is related to the apparent optical depth τ ′ by where 'apparent' refers to the quantity measured in the sky plane (see ref 11 and references therein).Knowing the ring pole orientation, these values can be converted to their values normal to the ring plane, p and τ .For cases where the Airy scale (see details in ref 37 ) is larger than the width of the ring as seen in the sky plane and assuming a monolayer ring, p is given by and assuming a polylayer ring, τ is given by 38 where B is the ring opening angle, B = 0 • (resp.B = 90 • ) corresponding to an edge-on (resp.pole-on) viewing, while Quaoar's ring has currently B ∼ −20 deg (Extended Data Table 2).The integrals of p and τ over the radial width W r of the rings define the equivalent width E p = pW r and the equivalent depth A τ = τ W r of the profile.The values of E p and A τ is proportional to the amount of material present in the profiles in the monolayer and polylayer case, respectively 10 .
If the ring profile is not resolved (i.e. drop in flux occurs over less than three data points), p ′ is not known.However, the integral E p is still measurable 10 .This is the case of the CHEOPS detections that has only one or two data points within the flux drops (Supplementary Fig. 1).
A noteworthy observational result is the large variation of W r and τ observed among detections.For instance, the 2021 August 27 occultation shows a narrow and dense ring with W r ≈ 8 km and τ ≈ 0.5.This feature was detected from three different sites, covering about 5.1 degrees along the ring, corresponding to a local and dense ring-arc feature of at least 365 km in length.Conversely, the 2019 June 5 event reveals a wider and more transparent feature (W r ≈ 300 km, and τ ≈ 0.014).Meanwhile, the equivalent widths and equivalent depths vary in narrower ranges of about 0.3-1.7 km and 0.3-3.8km, respectively (Table 1 and Extended Data Table 2).
Finally, the timings for the secondary events obtained at the various stations provide their sky-plane positions, to which an elliptical model is fitted using the χ 2 -statistical test, with where r i is the observed distance of the i th data point to the ring center, r ′ i is the corresponding distance of the model (the ellipse) and σ i is the uncertainty on the distance associated with the timing uncertainty of the i th data point.Finally, the σ model parameter accounts for extra uncertainties associated with our model 11,36 , stemming in particular from the unknown shape of Quaoar itself, resulting in model bias on its center position.Considering we have a single chord event that we fit an area equivalent circle of radius 555 km, but Quaoar's shape may be a Maclaurin spheroid 4 with an equatorial radius of 569 km and an oblateness of 0.087, this would cause a bias in the central position typically of ∼27 km, the value we choose for the σ model .
Here, we assume that Quaoar's ring is circular, so that its apparent flattening f ′ projected in the sky plane is related to its opening angle B through Assuming that Quaoar's ring pole orientation has been fixed over the 2018-2021 time interval, we used equation (5) to test a range of pole orientations and ring radii.This provides the two mirror solutions given in Table 1.Both solutions yield a satisfactory fit to the data, with best-fit χ 2 value per degree of freedom of χ 2 pdf = 0.28 and χ 2 pdf = 1.02, respectively.Solution 2 corresponds to a ring pole orientation that is consistent with Weywot's orbital pole orientation to within 6 ± 8 degrees, while Solution 1 results in a difference of 43 ± 8 degrees.Because we expect the ring and satellite to be coplanar, Solution 2 is our preferred one.
In summary, Quaoar's ring has a normal optical depth τ ranging from 0.004 to 0.77 depending on the longitude (Extended Data Tables 2 and 3).Using an impact frequency of ∼ 20τ impacts per particle and per orbit 39 , this implies that Quaoar's ring particles suffer between 1 and 10 collisions per revolution, which qualifies this ring as a dense one, that is a dynamics dominated by collisions.As the viewing geometry of the ring has changed very little between 2018 and 2021, and due to the paucity of observations, it is not possible yet to discriminate between a monolayer or a polylayer ring.
As a first approximation, we can parallel Quaoar ring with Saturn's rings, in that case, τ ∼ 0.25 corresponds to typical surface densities of Σ ∼ 500 kg m −3 (ref. 40).Adopting a radial width of W r = 8 km for the densest part of the ring, and considering that the equivalent width of the ring remains roughly constant in longitude (Extended Data Table 2) -so that its mass per unit length is also roughly constant -we obtain a ring mass estimation of M r = 2πaΣ ∼ 10 14 kg, where a ∼ 4, 100 km is the ring radius.If accreted into a single satellite with bulk density of 400 kg m −3 , typical of the small inner saturnian satellites 21 , this would yield body with a radius of the order of 5 km.
3 Multi-band profiles of the ring from GTC.
The 2019 June 05 event was monitored at the 10.4-meters Gran Telescopio Canarias (GTC) with the HiPERCAM instrument, using a four band system with g s (0.40-0.55 µm), r s (0.55-0.69 µm), i s (0.69-0.82 µm), and z s (0.82-1.00 µm) filters 41 .HiPERCAM also recorded the event in the u sband (0.3-0.4 µm), but the occulted star could not be detected at these wavelengths.Extended Data Fig. 1 displays the light curves obtained in each band, normalised between the unocculted and the zero stellar fluxes, respectively.
As expected, the noise decreases as the wavelength increases because the occulted star is brighter in the red.We fitted the ring parameters for each filter individually, see Table 1 and Extended Data Table 2.Even though we are fitting a ring detection with uniform opacity, we note that substructures may exist in some GTC profiles and thus may affect the obtained parameters, although any difference should be within the uncertainty of the parameters.No significant differences between the bands are observed in the parameters at the 3σ confidence level.Also, χ 2 -tests were performed to compare pairs of observations around the ring detections.The largest difference showed a p-value of 0.021, meaning that there is a small 2.1% chance that the light curve obtained with the r s (or z s ) filter is different from the g s filter.All the remaining possibilities had a p-value smaller than 0.2%.
The absence of wavelength dependence indicates that Quaoar's ring contains a substantial population of particles larger than around 10 µm 12,13,14 , which is in line with the dense and narrow rings of Saturn, Uranus and Chariklo.
4 Search for material in other light curves.
No light curves besides those described above show significant drops at the times expected from the ring geometry.This essentially stems from the lack of photometric sensitivity and/or temporal resolution of these data.We derived the detection limits provided by these data sets by calculating the standard deviation of the equivalent width where ϕ(i) is the observed flux and ∆r(i) is the radial interval covered by each exposure in the ring plane 37 .
The limits of detection were evaluated using the data at their original spatial resolutions, and also after a resampling over 300-km windows, which corresponds to the widest ring profile in this project.When applied, resampling allows for increased sensitivity to small variations caused by large, diffuse structures.This procedure consists of applying a Savitsk-Golay digital filter on curves of equivalent width as a function of radial distance in the plane of the ring.
For example, W. Hanna's data set (27 August 2021) from Yazz had a closest distance approach to Quaoar's center of 1,690 km counted in the ring plane and covered a region about 20,000 km before the event, up to 16,000 km after it.No detection of secondary events were found up to a 3σ level of E p =166 m per data point, with either an opaque ring (p=1) with radial width of ∼165 m, or a semi-transparent ring with W r =2.7 km and p =0.06, or any other intermediate solutions.However, a wider and more transparent (p ⩽ 0.01) feature similar to what was observed with the GTC would be lost in the noise.All the obtained upper limits values can be found in Extended Data Table 3.
Besides the four events used here, we also revisited ten other Quaoar occultations observed since 2011 (ref. 42) to search for secondary events.None of them provided data with sufficient quality to detect the ring.

N-body simulations of accretion.
We use numerical simulations with the local method where a co-moving ring patch of size 160 × 160 particle radii with periodic boundary conditions is followed to assess the accretion of particles for the proposed Quaoar ring.The particle impacts are calculated with a soft-particle method 43 .The number particles is N =2,000 in simulations with τ = 0.25 and N =8,000 when τ = 1; here τ stands for the dynamical optical depth defined as the total cross section of particles divided by the area of the simulated region.Such values of N are sufficient, since we are only interested in the onset of accretion, not the subsequent growth or mutual evolution of the aggregates, The advantage of a small N is that we can use a direct particle-particle method for calculation of self-gravity, to ensure that short-range binary interactions are accurately modeled.Typically, simulations lasted for 50-500 orbital periods (corresponding to time scales of a few years).
Several simulations were performed, varying the bulk density of particles ρ, and for a range of elasticity models (see below).We used Quaoar's mass M Q and ring radius a given in Table 1.Identical particles with a radius R=1 m were assumed in most of our runs.The dimensionless r H parameter characterizes the strength of self-gravity against the tidal force of the planet.It is defined as the ratio of the mutual Hill-radius (R Hill ) of a pair of particles to the sum of their physical radii, where µ = M 1 /M 2 = (R 1 /R 2 ) 3 is the mass ratio of the particles and a.
If we consider a test particle at rest on the surface of a synchronously rotating spherical large particle (i.e., µ = 0), the case r H (0) = 1 corresponds to the limiting case of a zero net attraction (tidal + self-gravity + centrifugal).Thus, in general, accretion (resp.disruption) of the pair is expected for r H larger (resp.smaller) than unity.In the case of equal-mass particles (µ = 1), we have where r H (1) is denoted by r H for simplicity.In the initial state of the simulation the radial velocity dispersion c exceeds by a large factor its steady-state value and also the 2-body escape velocity, During the simulation, c gradually decreases from its initial value due to collisional dissipation, until a steady-state value c st is reached, when the dissipation and viscous gain balance each other, or until c drops to about 2v esc , in which case the particles start to accrete and rapidly coalesce into a single aggregate, see Extended Data Fig. 3.This transition from non-accreting to accreting behavior is very sharp when ρ is increased.In case of accretion, we record the value of c at the onset of aggregate growth, in practice at the instant when the impact frequency had sharply increased 2-fold compared to its average value.Usually, the aggregate becomes visually evident within one orbital period from this instant of time.
Elasticity models.Extended Data Fig. 4 illustrates the ϵ(v n ) models used in our simulations.It also shows the theoretical relation between optical depth τ and the critical ϵ cr required for a balance between dissipation and viscous gain of energy 44 .The c st depends on the ϵ(v n ) model via the ϵ cr (τ ) relation and is independent of R as long as systems with c st ≫ Rn are considered.
The specific values used in simulations, τ = 0.25 and τ = 1.0, are highlighted and the corresponding values of ϵ cr (τ ) are marked on the left-hand plot.The small downward arrows for Model 4 illustrate the expected drop in the average steady-state v n (proportional to c st ) when optical depth increases: basically this drop follows from less effective viscous gain at large τ due the reduction of mean free path between impacts.Similarly, for a fixed τ , the Models 3, 2 and 1 imply successively smaller average v n 's and c st (refs. 45,29 . Accretion in simulations.The classical Roche limit applies for tidal disruption of a fluid satellite when brought gradually closer to the central body.It can be written in the form of Roche critical density for a given distance (equation (1)), where γ = 0.849 as in the original analysis by Roche 46 .However, besides a and ρ, the accretion of ring particles also depends on c, and thereby on the elasticity of particles.
Previous local N-body simulations 26 explored the onset of particle accretion in Saturn's rings, using both constant ϵ and the ϵ(v n ) relationship of Model 1 (ref. 25).The dependence on distance and bulk density was parameterised with the r H parameter (equation ( 8)).These simulations showed that r H > 1.1−1.2leads to accretion: the upper value is for identical frictionless particles like in the current experiments, whereas the lower value includes simulations with friction and/or particle size distribution.Note that the above condition for accretion is more stringent than the tidal destruction limit γ = 1.6 usually adopted (ref. 20), which corresponds to r H > 0.87.However, it is quite close to the classical value γ = 0.849 which corresponds to r H > 1.07.
All simulations of ref. 26 were limited to elasticity models (constant ϵ ≤ 0.5 or the relationship ϵ(v n ) of Model 1) which correspond to fairly inelastic particles.In this case, the balance between dissipation and viscous gain leads to small velocity dispersion for all ring optical depths.In particular, the velocity dispersion maintained by impacts alone is between ∼ 2Rn and ∼ 3Rn (where n is the mean motion) and is of the same order as the 2-body escape velocity.However, these highly dissipative models were compared in ref 29 to impact models with more elastic particles, and corresponding to our Model 4 (based on Fig. 22 of ref. 8 ).In this case, no accretion was possible even with r H = 1.23, which was the largest value explored.This stems from the high velocity dispersion which now exceeds the 2-body escape velocity by a large margin.
The simulations presented here extend these studies to larger r H .The behavior in N-body simulations is also consistent with 3-body integrations 24 which indicate that the probability of sticking in binary impacts increases rapidly when r H ≳ 1.2.The 3-body integrations (see Fig. 14.18 of ref. 45 ) also indicate that for large r H , the accretion probability goes rapidly down when impact velocities exceed 2-body escape velocity, in agreement with our current N-body experiments.
6 The dynamical environment of Quaoar's ring.
Table 1 lists the orbital radii of the 1/3 spin-orbit resonance (SOR) with Quaoar and the 6/1 mean motion resonance (MMR) with Weywot, assuming Keplerian motion around a spherical Quaoar of mass M Q .The error bars on the resonance locations are dominated by the uncertainty on M Q .Accounting for a possible large oblateness of Quaoar would shift these locations by a few kilometers only, leaving our conclusions unchanged.Moreover, Quaoar's rotational light curves yield two possible rotation periods 30 .The short rotation period (8.8394 ± 0.0002 h) would correspond to a single-peaked lightcurve, which would be caused by an oblate body with albedo features on its surface.However, we know that Quaoar cannot be an oblate spheroid but a triaxial ellipsoid due to the varying projected shape of the main body that we have determined in several stellar occultations by Quaoar (whose analysis is beyond the scope of this paper).Because triaxial ellipsoids give rise to double peaked lightcurves, the preferred rotation period is 17.6788 ± 0.0004 h.Besides, this double-peaked solution provided a better fit to the photometric data 30 .
The first-order (in the particle orbital eccentricity) SOR resonances -also called Lindblad resonances -exert torques that clear over short time scales (a few Myr) the region straddling the corotation region near 2,020 km, where the orbital period of the particles matches the rotation period of the body 9 .This clearing proceeds up to the outermost 1/2 Lindblad resonance near 3,200 km.Moving outwards the next resonance is the second-order 1/3 SOR, that occurs at a semi-major axis of a 1/3 = 4, 197 ± 58 km.This matches within the error bars with the two possible solutions for the ring orbital radius a R1 = 4, 097.3± 9.5 km and a R2 = 4, 148.4 ± 7.4 km, at the 1.7 and 0.9σ levels, respectively.Similarly, the 6/1 Weywot mean motion resonance is at a 6/1 = 4, 021 ± 57 km, again matching the two solutions a R1 and a R2 , at the 1.2 and 2.1σ levels, respectively.
We now consider resonances between a massless ring particle and either a Quaoar's mass anomaly or Weywot.For simplicity, Quaoar's mass anomaly is described as a hemispheric mountain at the surface of the body, see main text, however, it may also be an internal feature within an otherwise oblate, smooth body.We denote a, e, λ, ϖ, n and κ = n − π′ the semi-major axis, orbital eccentricity, mean longitude, longitude of pericenter, mean motion and epicyclic frequency of the particle, respectively.Depending on the case, λ ′ is the orientation of Quaoar's mass anomaly or Weywot's mean longitude, while n ′ is either Quaoar's spin rate or Weywot's mean motion.In the case of Weywot, the orbital eccentricity and the longitude of pericenter of the satellite, e ′ and ϖ ′ , must also be accounted for (contrarily to the mass anomaly, since the latter moves along a circle).Then κ ′ = n ′ − π′ is Weywot's epicyclic frequency.The quantities a 0 and n 0 are the radius and the mean motion at exact resonance.Finally, Quaoar's and Weywot's mass anomaly are denoted by µ Q and µ W , respectively, both normalized to Quaoar's mass.
We consider the resonance condition n 0 /n ′ = m/(m − j), where m is a positive or negative integer and j is positive.This resonance splits into j + 1 resonances, each described by a Hamiltonian of the form where ψ k = mλ ′ −(m−j)λ−kϖ −(j −k)ϖ ′ , and where k = 0, ...j is the order of the resonance in the particle's eccentricity.The case k = 0 corresponds to a corotation-type resonance with critical angle ϕ 0 = ψ 0 = mλ ′ −(m−j)λ−jϖ ′ , while the cases k ̸ = 0 correspond to eccentricitytype resonances with critical angles Each value of k is in fact associated with a mean motion resonance n 0 /n ′ = (m−j +k)/(m−j) of order k between the particle and a potential moving at the pattern speed For k ̸ = 0, we define the eccentricity vector as (X, Y ) = [e cos(ϕ k ), e cos(ϕ k )], where e 2 = X 2 + Y 2 .The Hamiltonian H is parameterized by the Jacobi constant where ∆a = a − a 0 is the distance to the resonance.Particles with the same value of ∆J but different initial conditions for (X, Y ) follow level curves of H under the Hamiltonian flow Ẋ = −∂H/∂Y and Ẏ = +∂H/∂X, eventually forming the phase portrait of the resonance.
The resonant forcing is encapsulated in the perturbing term containing ϵ k in the right-hand side of equation (10), not to be confused with the coefficient of restitution ϵ used in the main text and in Section 5.It is proportional to µ Q or µ W , and its numerical values is derived from expansion of the disturbing potential of the mass anomaly 48 or the satellite 49 .
The 1/3 SOR resonance corresponds to m = −1 and j = 2 in equation (10), so that n/Ω = 1/3, where the pattern speed Ω is Quaoar's spin rate.Because the mass anomaly moves on a circle, only the resonance k = j (=2 here) is allowed, so that where the subscript 1/3 refers to the 1/3 SOR, and where we have omitted, for brevity, the index k = 2 that should appear for the coefficient ϵ.Similarly, the resonant angle is denoted ϕ 1/3 = (−λ ′ + 3λ − 2ϖ)/2, and ∆J = [∆a/a 1/3 − (3/2)e 2 ]/2.From the expression of H(X, Y ; ∆J), it results that for an interval of ∆J of width 8|ϵ 1/3 |/9 centered on the resonance, the origin (X, Y ) = (0, 0) of the phase portrait is a fixed hyperbolic (and thus unstable) point.Collisions tend to damp the eccentricities, i.e., force the particles to move towards the origin of the phase portrait (i.e., on circular orbits).Conversely, in the interval mentioned above, the resonance tends to force them to follow eight-shaped trajectories, and thus, to acquire non-zero eccentricities.
From the expression of ∆J, this happens for a interval of semi-major axes of width W = (16|ϵ 1/3 |/9)a 1/3 centered on the resonance value a 1/3 .
The forced eccentricity then reaches a peak value (Extended Data Fig. 6) of e peak,1/3 = 8 3 (11) where we have used the expression of ϵ 1/3 in Extended Data Table 4.
In the second equation, we have assumed that the mass anomaly takes the form of a hemispheric mountain of height h at Quaoar's surface, see main text.
The peak eccentricity e peak,1/3 tends to maintain a velocity dispersion ∆v among the particles or putative parent bodies of Quaoar's ring (main text).If the motions of these bodies are not coherent, we have ∆v = 2e peak,1/3 v orb , where v orb = GM Q /a 1/3 is the orbital velocity at the resonance.The escape velocity at the surface of a particle with radius R p and density ρ is v esc = 8πGρR p /3. Consequently, the velocity dispersion ∆v is comparable to v esc for where we have used equation (11), ρ = 400 kg m −3 and the numerical values of Table 1.
In the presence of collisions, the second-order nature of the 1/3 SOR leads to mathematical difficulties.In particular, the periodic resonant streamlines intersect at one point, even for vanishing eccentricities 48 .This yields to a multi-valued velocity field and infinite densities, and therefore singularities in the hydrodynamical equations.
Considering the simple case of a mass anomaly in the form of a hemispheric mountain of height h on Quaoar's surface, the mass anomaly normalized to Quaoar's mass is µ Q = (h/R Q ) 3 /2, where R Q is Quaoar's radius.The 1/3 SOR resonance forces an eccentricity e peak,1/3 (equation (11)) which can in turn maintain locally a velocity dispersion among the putative parent bodies in Quaoar's ring (equation (12)).This velocity dispersion is comparable to the escape velocity at the surface of the parent bodies for a certain value of µ Q (equation (13)), or equivalently, for a certain mountain height of Assuming parent bodies with size R p ∼ 0.1 km, typical of Saturn's F ring parent bodies, we obtain h ∼ 10 km.Thus, plausible topographic features on Quaoar are indeed able to maintain a velocity dispersion among the parent bodies that prevent them from accreting near the 1/3 SOR.This is the case a fortiori for ring particle aggregates that have much smaller sizes.This assumes, however, that the resonant responses of the bodies are not coherent.In the opposite case, the velocity dispersion in equation (12) could be largely overestimated.Meanwhile, more realistic numerical simulations using N-body collisional codes do show that confinement of a collisional disk (without self-gravity) is observed near the 1/3 SOR with Chariklo's 28,27 .This confinement is associated with angular momentum flux reversal at certain longitudes of the ring, but is not yet backed up by analytical calculations.
The 6/1 MMR Weywot resonance corresponds to m = 6 and j = 5 in equation (10).It splits into six resonances as k takes the values 0,...5, with the respective resonant angles listed in Extended Data Table 4.The associated resonant terms of Weywot's disturbing potential have amplitudes of the form ϵ 6/1,k e k e ′5−k , where the coefficients ϵ 6/1,k are given in Extended Data Table 4.
Corotation resonance.The case k = 0 corresponds to the corotation resonance condition n = n ′ + 5κ ′ ≈ 6n ′ .It creates one stable elliptic (L 4 or L 5 -type) Lagrange point and one unstable hyperbolic (L 3 -type) point.Here we consider the possibility that the accumulation of ring material around the L 4 point might explain the longitudinal variability of Quaoar's ring.Classical calculations provide the full width of the corotation zone: that corresponds to the spread in semi-major axes of the particles trapped in this corotation resonance.Using the value of ϵ 6/1,0 (Extended Data Table 4), a 0 ≈ 4020 km (Table 1), and assuming that Quaoar and Weywot have the same bulk density, we obtain where R W is Weywot's radius.Both the quantities R W and e ′ (Weywot's orbital eccentricity) are poorly constrained.Assuming the same albedo for Quaoar and Weywot, we obtain R W ≈ 40 km (ref. 5), and Weywot's orbital solution provides e ′ < 0.15 (using the method presented in ref. 34 ).
The equation above then yields W COR ≲ 1 km.This is comparable to the spread in semi-major axis of Neptune's arcs, which is much smaller than the physical width of Neptune's ring-like arcs, about 15 km.This is classically explained by the coherent radial motion of Neptune's arc particles forced by the nearby satellite Galatea 50 .This process might also apply to Quaoar's ring, but still remains to be established.Meanwhile, we note that W COR rapidly decreases with e ′ .For instance, for e ′ = 0.01, we obtain W COR ∼ 1 m, which appears to be too narrow for maintaining an arc structure around Quaoar, as W COR is then comparable to the expected typical particle sizes.
Lindblad resonance.The resonance with k = 1 corresponds to a Lindblad (first-order) resonance with associated disturbing potential ϵ 6/1,1 ee ′4 cos(ϕ 6/1,1 ).Technically, it is a 2/1 resonance of the particle with the component of Weywot's potential that has the pattern speed n ′ + 2κ ′ ≈ 3n ′ .The same exercise as for the Quaoar 1/3 SOR, but now considering a first-order resonance, shows that this resonance forces a peak eccentricity for particles that start with a circular orbit of The value of e peak,6/1 is again poorly constrained.However, we may compare e peak,6/1 and e peak,1/3 (equation 11) for R W = 40 km and h = 10 km.Then the two values of e peak are comparable for a Weywot's orbital eccentricity as small as 0.005.In other words, Weywot's is expected to force eccentricities of the ring at the 6/1 Lindblad resonance that are comparable to or larger than the eccentricities forced by Quaoar's 1/3 SOR.
Other resonances.The cases k = 2, ...5 correspond to four additional resonances of orders k in the particle eccentricities.At this stage, it is difficult to assess their interactions and effects on the corotation and Lindblad resonances, and thus, on the ring behavior.In particular, their mutual radial distances depend on the value of the apsidal precession rate π′ , which is unknown.Moreover, their couplings also depend on the value of e ′ , which is largely unconstrained.9 The effect of radiation pressure.
Even though Quaoar is located in a faraway region compared to the giant planets, the weak gravity field of the central body allows micrometric ring particles to be strongly disturbed by the radiation pressure (RP).To first order, this force causes no secular effect in the semi-major axis of the particles, but induces periodic oscillations in the eccentricity.In a first approximation, the maximum eccentricity e RP attained by a particle of radius R, bulk density ρ initially in a circular orbit of radius a is 51 where n and n Q are the mean motion of the particle and Quaoar, respectively, c ′ is the speed of light, F ⊙ is the solar flux at Quaoar's orbit, and Q pr is the radiation pressure coefficient, equal to unity for an ideal material.
This equation shows that ring particles smaller than about 25 µm in radius are ejected from the system, while grains with radius ≲ 45 µm reach eccentricities large enough to collide with Quaoar.In fact, particles smaller than 0.2 cm suffer radial excursions that surpass the ring width.
These results are backed up by numerical simulations performed with a modified version of the Mercury package 52 that includes the radiation pressure effects.These simulations show that 1 µm-particles are ejected from the system in less than a couple of years, while ≲ 45 µm do not survive over a few decades.These results show that Quaoar's ring is likely to be quickly depleted of sub-cm grains, in agreement with the GTC multi-filter observations.elasticity, the system adjusts its impact velocities (via velocity dispersion) so that the effective mean ϵ corresponds to the ϵ cr (τ ).In case of constant ϵ < ϵ cr the system flattens to a nearmonolayer state with a minimum c ≈ 2 to 3Rn ≈ 0.01 cm s −1 R/1 m (where n is the mean motion) determined by the non-local viscous gain associated with the finite size of the particles.If the constant ϵ > ϵ cr , no thermal balance is possible and the system disperses via exponentially increasing c.
Extended Data Figure 5 | The influence of optical depth and particle size.In the upper panel, we show Model 4 simulations with τ = 0.25 and τ = 1.0.In both series of simulations particle size is R = 1 m.The same conventions as in Fig. 3 are used, but now the scale for c is linear, not logarithmic.For larger τ the steady-state velocity dispersion is reduced, and thus the condition c st /v esc ≲ 2 is achieved for smaller ρ.The bottom panel compares the three assumed particle radii (0.33 m, 1 m, and 3 m) on accretion, using a common optical depth of τ = 0.25 and assuming the velocity-dependent elasticity law of Model 2. Since c st is nearly independent of R, the critical density corresponding to c st /v esc ∼ 2 scales roughly as ρ cr ∝ 1/R 2 since v esc ∝ √ R (eq. ( 9)).Note that if a constant ϵ ≲ 0.5 is assumed, the particle size has no effect on the accretion limit, since in this case c st scales linearly with R in a similar fashion to v esc .
Extended Data Figure 6 | Topology of the Quaoar 1/3 Spin-Orbit Resonance (SOR).The bottom graph shows the maximum eccentricity e max reached by a particle starting on an initially circular orbit of semi-major axis a perturbed by the Quaoar 1/3 SOR resonance.This resonance is driven by a mass anomaly whose amplitude is quantified by the dimensionless parameter ϵ 1/3 .The exact resonance radius a 1/3 is marked by the dashed vertical tick mark.The top plots show the phase portraits in the eccentricity vector space (X, Y ) corresponding to particular values of a, with X = e cos(ϕ 1/3 ), Y = e sin(ϕ 1/3 ), where e is the orbital eccentricity and ϕ 1/3 = (−λ ′ + 3λ − 2ϖ)/2 is the resonant critical angle, see Methods for details.In an interval of width W = (16|ϵ 1/3 |/3)a 1/3 in semi-major axis centerd on the resonance, the origin of the phase portrait is an unstable hyperbolic point.Particles are then forced to reach a maximum eccentricity e max = 4/3 (W/2 − ∆a)/a 1/3 , where ∆a = a − a 1/3 is the initial distance of the particle to the resonance.The value of e max peaks at e peak,1/3 = (8/3) |ϵ 1/3 |/3 for ∆a = −(8/9)|ϵ 1/3 |a 1/3 .Outside this interval, e max = 0. Units are arbitrary in all the plots.Extended Data Table 4 | Coefficients ϵ k (Eq. 10 and refs 49,48 ) associated with the Quaoar 1/3 Spin Orbit Resonance (SOR) and the Weywot 6/1 Mean motion Resonances (MMRs).Note: depending on the line, λ ′ denotes either the orientation of Quaoar's mass anomaly or Weywot's mean longitude.The corotation term ϵ 6/1,0 accounts for both the direct and indirect parts of Weywot's disturbing potential.See the Methods for the definition of the quantities.

Figure 1 |
Figure1| Example of Quaoar's rings detections.The observed flux (black points) and the models (red lines) are plotted against the time relative to the observer closest approach.The blue shaded regions are enlarged in corresponding underlying panels.Panels a, b and c: the light curve observed by HiPERCAM (i s band) on the 10.4m Gran Telescopio Canarias (GTC) with closest approach time at 03:00:31.858UTC (05 June 2019).Panels d, e and f: the light curve observed in Algester by R. Langersek (citizen astronomer) with closest approach time as 10:59:00.442UTC(27 August 2021).No definite detection is made after closest approach, only an upper limit is derived.More detections are displayed in Supplementary Fig.1.

Figure 2 |
Figure 2 | Sky-plane projection of all our Quaoar's ring detections.The occulting chords (green dashed line) and the ring detections with their 1-sigma error bars (red lines) are plotted for the 02 September 2018, 05 June 2019, 11 June 2020 and 27 August 2021 events.The black ellipses are the two simultaneous best-fitting solutions (Table1), which are indistinguishable at this scale.The J2000 celestial North and East directions are shown in the upper right corner of each panel.Please note, that one of the chords in the 11 June 2020 occultations was obtained with CHEOPS space telescope, details in ref31 .
initiative funded by HEFCE and by the Wolfson Foundation J.d.W. and MIT gratefully acknowledge financial support from the Heising-Simons Foundation, Dr. and Mrs. Colin Masson and Dr. Peter A. Gilman for Artemis, the first telescope of the SPECULOOS network situated in Tenerife, Spain.The ULiege's contribution to SPECULOOS has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) (grant Agreement n°336480/SPECULOOS), from the Balzan Prize and Francqui Foundations,

Table 1 -
Physical parameters of Quaoar, Weywot and the ring, with their 1-σ error bars.

Table 1 )
31which are indistinguishable at this scale.The J2000 celestial North and East directions are shown in the upper right corner of each panel.Please note, that one of the chords in the 11 June 2020 occultations was obtained with CHEOPS space telescope, details in ref31.
Figure 3 | Local simulations of rings outside the classical Quaoar's Roche limit.Results of self-gravitating local simulations performed with Quaoar's mass and ring radius (Table