ALMA Reveals an Eccentricity Gradient in the Fomalhaut Debris Disk Joshua B. Lovell1aa, Elliot M. Lynch2, Jay Chittidi3aa, Antranik A. Sefilian4aa, Sean M. Andrews1aa, Grant M. Kennedy5aa, Meredith MacGregor3aa, David J. Wilner1aa, and Mark C. Wyatt6aa 1 Center for Astrophysics, Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138-1516, USA; joshualovellastro@gmail.com 2 Univ Lyon, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230 Saint-Genis-Laval, France 3 Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA 4 Department of Astronomy and Steward Observatory, University of Arizona, Tucson, AZ 85721, USA 5 Department of Physics, University of Warwick, Coventry CV4 7AL, UK 6 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK Received 2024 September 30; revised 2025 June 10; accepted 2025 July 8; published 2025 September 4 Abstract We present evidence of a negative eccentricity gradient in the debris disk of the nearby A-type main-sequence star, Fomalhaut. Fitting to the high-resolution, archival ALMA 1.32 mm continuum data for Fomalhaut (with a synthesized angular resolution of 0.76 × 0.55; 4–6 au), we present a model that describes the bulk properties of the disk (semimajor axis, width, and geometry) and its asymmetric morphology. The best-fit model incorporates a forced eccentricity gradient that varies with semimajor axis, e af npow, a generalized form of the parametric models of E. M. Lynch & J. B. Lovell, with npow = −1.75 ± 0.16. We show that this model is statistically preferred to models with constant forced and free eccentricities. In comparison to disk models with constant forced eccentricities, negative eccentricity gradient models broaden disk widths at pericenter versus apocenter, and increase disk surface densities at apocenter versus pericenter, both of which are seen in the Fomalhaut disk, and which we collectively term “Eccentric Velocity Divergence.” We propose single-planet architectures consistent with the model and investigate the stability of the disk over 440Myr to planet–disk interactions via N-body modeling. We find that Fomalhaut’s ring eccentricity plausibly formed during the protoplanetary disk stage, with subsequent planet–disk interactions responsible for carving the disk morphology. Unified Astronomy Thesaurus concepts: Debris disks (363); Circumstellar disks (235); Eccentricity (441); Planetary-disk interactions (2204); Planetary system evolution (2292) 1. Introduction Debris disks are cold (tens of K) dust belts formed via mutual planetesimal collisions (M. C. Wyatt 2008; B. C. Matt- hews et al. 2014; A. M. Hughes et al. 2018; S. Marino 2022). While such collisions grind down disks over time, debris disks can survive well into (and beyond) the stellar main sequence, producing observable dust signatures over 10 s–1000 s of Myr. Planet–disk interactions significantly alter the structures of debris disks, e.g., via depleting, overpopulating, and warping planetesimal distributions. On secular timescales, eccentric planets force their own eccentricities into the orbits of planetesimals via three-body interactions (i.e., in a star– planet–planetesimal three-body system C. D. Murray & S. F. Dermott 1999; A. J. Mustill & M. C. Wyatt 2012). In resolved observations of disks, diagnostics such as the stellar offset from the disk center, relative brightness asymmetries, and relative width asymmetries have all been interpreted as metrics to derive disk eccentricities and eccentricity distribu- tions (M. C. Wyatt et al. 1999; M. Pan et al. 2016; E. M. Lynch & J. B. Lovell 2021; J. B. Lovell & E. M. Lynch 2023). The most famous, and best-studied, eccentric debris disk, is that around the 7.66 pc-distant, 440Myr old, A-type main-sequence star, Fomalhaut (derived from the Arabic name “fam al-hūt (al-janūbī)” and translates to the “mouth of the (southern) whale”; R. H. Allen 1963). Fomalhaut (HD 216956) is the brightest member of southern Piscis Austrinus constellation, though its debris has only been known of since 1985 (via analysis of its Infrared Astronomical Satellite (IRAS) observations; H. H. Aumann 1985). Since then, it has been resolved at higher-resolution with the Hubble Space Telescope (HST; in optical scattered light; P. Kalas et al. 2005, 2008; A. Gaspar & G. Rieke 2020), the Spitzer Space Telescope and James Webb Space Telescope (JWST; in near- and mid-infrared thermal emission; K. R. Stapelfeldt et al. 2004; A. Gáspár et al. 2023), the Herschel Space Telescope (in far-infrared thermal emission; B. Acke et al. 2012), and the James Clerk Maxwell Telescope and Atacama Large Millimeter/submillimeter Array (ALMA; in (sub) millimeter thermal emission; W. S. Holland et al. 2003; A. C. Boley et al. 2012; L. Matrà et al. 2015; M. A. MacGregor et al. 2017; L. Matrà et al. 2017). Most recently, both G. M. Kennedy (2020) and J. Chittidi et al. (2025) presented evidence of a width asymmetry in the disk, with a broader pericenter width versus apocenter. These observations defini- tively proved the eccentric nature of the debris disk: the star is offset from the geometric disk center (in both scattered light and thermal emission) and the disk presents wavelength- dependent disk ansae brightness asymmetries, specifically exhibiting “pericenter glow” in infrared emission (see e.g., M. C. Wyatt et al. 1999) and “apocenter glow” in millimeter emission (see e.g., M. Pan et al. 2016). Both of these latter features are expected based on the structure and thermal conditions in Fomalhaut’s disk at the observed resolution scales (E. M. Lynch & J. B. Lovell 2021). Here, we present a millimeter model of Fomalhaut’s main (outer) debris disk belt that incorporates a forced eccentricity gradient parameter, which we fit to the high-resolution ALMA The Astrophysical Journal, 990:145 (18pp), 2025 September 10 https://doi.org/10.3847/1538-4357/adfadc © 2025. The Author(s). Published by the American Astronomical Society. aaaaaaa Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 1 https://orcid.org/0000-0002-4248-5443 https://orcid.org/0000-0002-4985-028X https://orcid.org/0000-0003-4623-1165 https://orcid.org/0000-0003-2253-2270 https://orcid.org/0000-0001-6831-7547 https://orcid.org/0000-0001-7891-8143 https://orcid.org/0000-0003-1526-7587 https://orcid.org/0000-0001-9064-5598 mailto:joshualovellastro@gmail.com http://astrothesaurus.org/uat/363 http://astrothesaurus.org/uat/235 http://astrothesaurus.org/uat/441 http://astrothesaurus.org/uat/2204 http://astrothesaurus.org/uat/2292 https://doi.org/10.3847/1538-4357/adfadc https://crossmark.crossref.org/dialog/?doi=10.3847/1538-4357/adfadc&domain=pdf&date_stamp=2025-09-04 https://creativecommons.org/licenses/by/4.0/ images of J. Chittidi et al. (2025). This model provides a novel interpretation of the disk’s morphology; incorporating a steep, negative eccentricity gradient with semimajor axis. In Section 2, we outline a physical description of the “eccentric velocity divergence” (EVD) effect that describes this type of disk structure. In Section 3, we describe our model setup, fitting methodology, and results. In Section 4, we derive planet properties consistent with observations. In Section 5, we consider the stability of these planetary scenarios with N-body modeling. We summarize our conclusions in Section 6. 2. Eccentric Velocity Divergence Central to the investigation that we outline in this work is the dependency between disk surface density distributions and their underlying eccentricity distributions. This has been demonstrated recently in the literature by E. M. Lynch & J. B. Lovell (2021) where simple power-law forced eccen- tricity profiles ( )e a af npow with fixed npow values of 0, −1, or +1 were modeled (where disk eccentricity gradients were either flat (∂ef/∂a = 0), positive (∂ef/∂a > 0), or negative (∂ef/∂a < 0), respectively). These models were shown to alter the resulting disk surface densities and emission morphologies. In this work, we adopt the generalized form ( )e a af npow (allowing for noninteger values) to examine power-law eccentricity profiles in the context of the Fomalhaut debris disk. We illustrate the influence of this eccentricity distribution in Figure 1 showing normalized dust surface densities for a subset of power-law profiles (0.5, 0, −0.5, −1, and −2, defined by ( ) ( )/=e a e a af f n ,0 0 pow) adopting parameters con- sistent with the Fomalhaut disk (i.e., a semimajor axis of a0 = 139 au, a Gaussian disk width of σr = 15.6 au, and a forced eccentricity ef,0 = 0.127 at a0) and an arbitrarily chosen ωf = 0.0. These plots show the same behavior as described by E. M. Lynch & J. B. Lovell (2021), i.e., that surface densities in disks with positive eccentricity gradients are enhanced at pericenter, whereas disks with negative eccentricity gradients have surface densities raised at apocenter. In addition, these models show that disk widths decrease where these are most dense, and broaden where these are least dense. For example, negative eccentricity gradient profiles preferentially broaden disk widths at pericenter, and preferentially narrow disk widths at apocenter. To simplify discussion of the combined impact on disk surface densities and widths, we term this effect Eccentric Velocity Divergence (EVD), which simply arises due to mass continuity in eccentric disks. One can consider the effect on the disk width analytically by calculating the radial distance between two concentric orbits at pericenter (Δrperi) and apocenter (Δrapo). In the case of fixed positive eccentricities, the definitions of the pericenter and apocenter radii result in widths Δrperi = Δa(1 − e) and Δrapo = Δa(1 + e), i.e., Δrapo > Δrperi always for e > 0. However, for the general case of e = e(a) (retaining only linear- order terms in the expansion of e = e(a)), the same derivation finds Δrapo ± Δrperi ≈ Δa(1 ± (e + a∂e/∂a)). Consequently one can show that apocenter and pericenter widths are dependent on eccentricity gradients, with Δrapo > Δrperi requiring ∂e/∂a ≳ − e/a. For the case of ( )e a af npow, one finds thatΔrapo > Δrperi is only true if npow ≳ −1, which can be seen by comparing the ef(a) ∝ a−1 model (column four of Figure 1) with other presented models. Overall, EVD can induce similar features to those observed in the Fomalhaut debris disk, so we consider the possibility of its operation in this disk. Power-law eccentricity profiles are physically motivated. For example, in the idealized situation of an eccentric three- body system with either an inner-planetary perturber (i.e., a star-eccentric planet-planetesimal system), or an external- planetary perturber (i.e., a star–planetesimal–eccentric planet system) the disturbing function predicts analytical relation- ships between the semimajor axes and eccentricities of the perturbing planet and the planetesimal. To first-order (i.e., under the conditions of low eccentricity and well-separated star–planet–(massless) planetesimals), these relationships can be derived as either ef(a) ∝ a (for external-planetary perturbations) or ef(a) ∝ a−1 (for internal-planetary perturba- tions). These gradient slopes can steepen, for example, when e a0.5 r e a0 r e a 0.5 r e a 1 r e a 2 r Figure 1. Normalized surface density maps (face-on projections, in (top) r–f space, from 0 to 2π and 0 to 180 au, and (bottom) x–y space) for different power-law eccentricity profiles (with their e = ef(a) functional form shown in the lower right of each panel). The models all have apse-aligned arguments of pericenter, as well as radii, widths, and eccentricities consistent with those of Fomalhaut, and an argument of pericenter directed south or toward 0 or 2π in r–f space. Contours in increasing brightness scale represent regions of higher surface density, set at the levels 10%, 30%, 50%, 70%, and 90% (10% levels are shown with white-dotted contour lines). 2 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. the semimajor axis of a perturbing planet approaches those of its perturbers or if planetesimals are instead modeled with mass (see e.g., A. A. Sefilian et al. 2021; A. A. Sefilian 2024). Indeed, such eccentricity distributions are present in the Kuiper Belt, whereby the distribution of trans-Neptunian Objects in the cold (classical) belt host steep ef(a) profiles (e.g., at the level of a few percent per au; R. I. Dawson & R. Murray-Clay 2012). Moreover, theoretical work on eccentric protoplanetary (gas) disks has shown that eccentricity gradients are necessary to avoid violent differential precession in ( ) =e a const.f disks (J. Teyssandier & G. I. Ogilvie 2016). If debris disks form from primordial material with such eccentricity distributions, or are shaped by eccentric planet–disk interactions, it is plausible that their structures may be well described by the EVD effect. 3. Observations and Modeling 3.1. ALMA Observations For this study, we fit to the images of Fomalhaut from J. Chittidi et al. (2025). That work presents the calibration and data alignment procedure for the three ALMA band 6 data sets included (P.I.: Aaron Boley ADS/JAO.ALMA#2013.1.00486. S, P.I.: Paul Kalas ADS/JAO.ALMA#2015.1.00966.S, and P. I.: Meredith MacGregor ADS/JAO.ALMA#2017.1.01043.S) and utilizes standard CASA ALMA reduction pipelines. We use the same tclean parameters for imaging as J. Chittidi et al. (2025) with a Briggs-weighting scheme and robust parameter of 0.5, which produces a mosaic image centered on the 2015 coordinates of Fomalhaut, with a pixel size of 0.075 and a synthesized clean beam with 0.76 × 0.55, and beam position angle (PA) −87°.4. At the 7.66 pc distance to Fomalhaut, this image provides a physical resolution of 6 × 4 au. We present the image in Figure 2 (left). The three observations that were combined to produce this image were all conducted with different mosaics, sensitivities, and ALMA configurations (and thus effective angular resolutions). As a result, the final image has non-Gaussian noise properties and primary-beam attenuation that is dominated by the two-pointing mosaic of ADS/JAO.ALMA#2017.1.01043.S, which had phase centers directed toward the disk ansae. As can be seen in Figure 2, the data is therefore much noisier near to the minor axis of the disk in comparison to the two ansae. 3.2. Model Setup We produce parametric models of optically thin disks consistent with those presented in J. B. Lovell et al. (2021), E. M. Lynch & J. B. Lovell (2021), J. B. Lovell et al. (2022), and J. B. Lovell & E. M. Lynch (2023), i.e., by imaging three- dimensional dust distributions within the RADMC-3D framework (C. P. Dullemond et al. 2012). We fix the model wavelength to that of the band 6 observations (λ = 1.32mm). We fix the stellar parameters to values consistent with E. E. Mamajek (2012), i.e., the effective temperature (T� = 8500K), stellar radius (R� = 1.8R⊙), and stellar mass (M� = 1.9M⊙), which together fully define the Kurucz stellar template spectra (see R. L. Kurucz 1979) from which the temperature conditions through the disk are derived by RADMC-3D. The flux density at the location of the star (as adopted later) is a fixed value independent of these stellar parameters and the disk temperature conditions. The models assume the dust to have a fixed size distribution between 0.9 μm–1.0 cm with a −3.5 power-law index (as per J. S. Dohnanyi 1969), and with grain densities of ρ = 3.5 g cm−3. We adopt a single Gaussian semimajor axis distribution, based on the (mean) disk semimajor axis (a0) and Gaussian-width (wr, which is related to the disk full width at half-maximum (FWHM), by =FWHM 2 2 ln 2 wr). The sur- face density of the disk is modeled with a (mean) forced eccentricity (ef,0, measured with respect to a0), an argument of forced eccentricity (ωf, measured anticlockwise from south), and a power-law index of npow (for a power-law eccentricity distribution of the form ( ) ( )/=e a e a af f n ,0 0 pow), as per ( ) ( )= w j a a w 1 2 exp 2 , 1 r r 0 2 2 -10.00.010.0 -20.0 -10.0 0.0 10.0 20.0 R el at iv e D ec l. O ff se t [a rc se c] −20 0 20 40 60 80 100 120 -10.00.010.0 Relative RA Offset [arcsec] ef ∝ a−1.75 −20 0 20 40 60 80 100 120 -10.00.010.0 In ten sity [µ J y b eam − 1] −40 −30 −20 −10 0 10 20 30 40 Figure 2. Left: ALMA data as presented in J. Chittidi et al. (2025) which we fit to in this work. Center: Best-fit ef(a) ∝ a−1.75 model (on same image scale). Right: Residual emission after we subtract our best-fit model from the ALMA data (left). Contours are shown at the ±3σ level in the residuals and at the 5σ level for the data/model. Beams are shown in the lower-left of each plot. Emission remains present in the minor axis region of the residuals; however, this is dominated by imaging noise (enhanced due to the primary-beam response) and background sources (such as the “Great Dust Cloud,” see A. Gáspár et al. 2023 and G. M. Kennedy et al. 2023). North is up, east is left. 3 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. with j the coordinate system Jacobian determinant (as also presented as Equation (8) in E. M. Lynch & J. B. Lovell 2021) ( ) ( ) ( )/ = + j e e a e a e q E 1 1 1 cos 2 2 for the eccentric anomaly (E) ( ) ( ) ( )= + + E e e cos cos 1 cos 3 f f and the orbital intersection parameter (q, in the case of an untwisted disk, i.e., ∂ωf/∂a = 0) ( ) ( )/ / = + q a e a e e a e a1 . 4 To avoid orbital intersections, we require |q| < 1. Equation (1) fully describes the plots presented in Figure 1 (gridded on r–f and x–y space).7 We adopt a single Gaussian-distributed vertical density profile, with a vertical aspect ratio (h = H/r) for H the classical scale-height for a given disk radius r, such that the vertical linear density is ( )=l H z H 1 2 exp 2 . 5 2 2 This parameterization defines the total density ρ = MdustΣlρ, where we introduce a dust mass scaling parameter (Mdust) which uniformly scales model pixel intensities to mass units. By default, our models are face-on and oriented north–south (with the pericenter direction southwards, apocenter direction northwards, as per Figure 1). We use RADMC-3D to project and image models on the sky-plane. We provide parameters for the PA (which rotates our model anticlockwise referenced to north), inclination (i, which inclines our model with reference to angles between i = 0° for face-on projections, and i = 90° for edge-on projections), and offsets in the right ascension (ΔR.A.) and declination (Δdecl.), respectively (applied uniformly to both the stellar location and disk focus). We present in Figure 2 (middle) one such imaged model with Fomalhaut’s disk and stellar parameters (convolved with the same clean beam parameters as the imaged ALMA data, using the best-fit parameters inferred later). 3.3. Model Fitting Our fitting approach is consistent with the procedures of J. B. Lovell et al. (2021, 2022), except in this work we fit in the image domain instead of the visibilities.8 We fit our model to the ALMA data with emcee, an “Affine Invariant Markov Chain Monte Carlo Ensemble Sampler”9 (J. Goodman & J. Weare 2010; D. Foreman-Mackey et al. 2013). Samples of the posterior distribution (generated by emcee) are deter- mined from the assumed prior distributions and a Gaussian likelihood function ( ( )/exp 22 ), with ( )= = M I I I l , 6 i j 2 , 1 N data,i,j model,i,j 2 data 2 pix 2 beam i,j pix for an image with size Npix × Npix, intensity in each pixel (I, for either the model or data), the image noise (δIdata, see further below), the beam size (Ωbeam), image pixel size (lpix, where the latter term /lpix 2 beam reweights χ2 to account for pixel–pixel correlation), and the fitting mask (Mi,j). We fit to the primary- beam corrected image (Idata,i,j) and measure the χ2 in the residual map corrected by the primary-beam response, assuming the CASA form of ALMA’s primary beam (as described in e.g., ALMA Knowledgebase, Mar 202510). The ALMA image has a non-Gaussian noise distribution most prominent near to the disk minor axis given the image derives from different mosaic patterns from different ALMA config- urations (see Figure 2). By comparing fits to data with different masking conditions, we proceeded to fit with a mask defined by a cut-off in the ALMA imaged primary beam (only fitting emission where the primary-beam response exceeds 0.66), and a projected radius cut-off (only fitting to emission within a projected 200 au distance from the stellar center). We present in Appendix A (Figure 5) the location of the fitting mask and the primary-beam response for the data, as well as the discussion of and results of fitting the General model to the data without including the primary-beam cut-off in the map. The mask simply multiplies pixel values with either =M 1i,j or =M 0i,j and is applied to both the data and model identically. The mask choice focuses the fitting on high signal- to-noise ratio (SNR) disk emission, avoiding imaging artifacts and background sources near the disk semiminor axes (e.g., such as the “Great Dust Cloud” (GDC) as discussed by G. M. Kennedy et al. 2023). We note that by masking in the fitting stage rather than during model generation, our models still comprise azimuthally complete disks. We adopt an uncertainty of δIdata = 8.0 μJy beam−1 based on the rms of pixel intensities in disk-free regions of the ALMA image, which is scaled by the primary-beam response during fitting, down-weighting noisier regions of the data. We find from test iterations of fitting models to the data with emcee that the residual maps have standard deviation ≈8.0 μJy beam−1 (i.e., as expected given the image data) and mean <1.0 μJy beam−1 (i.e., very close to 0, as expected for well-fitted models). To verify the fitting methodology, we tested this procedure on the Cycle 3-only ALMA data of Fomalhaut as described in Appendix B. 3.4. Modeling Results 3.4.1. Fitting an EVD Disk Model We fit the stellar flux and location separately to the disk to reduce the number of free parameters in our fitting (since we simply fit the stellar flux, we leave fixed the physical stellar inputs that derive the disk temperature conditions, e.g., R�, M�, and T�). We adopt a small aperture mask over the full image to fit the stellar flux and location, masking all but a 10 × 10 pixel 7 The code for generating eccentric ring models and reproducing these plots is available via DOI: 10.5281/zenodo.16758671 and at https://github.com/ astroJLovell/eccentricDiskModels. 8 It is in principle possible to fit such models in the visibility domain; however, this demands a much greater computational time investment. We note that model convergence can be achieved with this image-based setup in a modest 1–2 days, whereas in the visibility domain a much more intensive 1–2 weeks per model is required (see J. Chittidi et al. 2025), and thus the approach here enables us to investigate a broad array of model setups. 9 I.e., see: https://emcee.readthedocs.io/en/stable/. 10 https://help.almascience.org/kb/articles/pdf/how-do-i-model-the-alma- primary-beam-and-how-can-i-use-that-model-to-obtain-the-sensitivity-pr 4 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. https://doi.org/10.5281/zenodo.16758671 https://github.com/astroJLovell/eccentricDiskModels https://github.com/astroJLovell/eccentricDiskModels https://emcee.readthedocs.io/en/stable/ http://help.almascience.org/kb/articles/pdf/how-do-i-model-the-alma-primary-beam-and-how-can-i-use-that-model-to-obtain-the-sensitivity-pr http://help.almascience.org/kb/articles/pdf/how-do-i-model-the-alma-primary-beam-and-how-can-i-use-that-model-to-obtain-the-sensitivity-pr box. We find that the star is best-fitted with offsets consistent with 0.0 in both axes (ΔR.A. = 0.01 ± 0.04 and Δdecl. = 0.01 ± 0.04, in accordance with the realignment and imaging procedure of J. Chittidi et al. (2025), and a flux of 0.737 mJy (i.e., F� = 0.74 ± 0.08 mJy accounting for a 10% ALMA flux calibration uncertainty in quadrature with our fitting error). We run an initial fitting routine with emcee (with 48 walkers, 3000 steps) to estimate the convergence timescale and to further verify the fitting results for the higher-resolution data. This initial fitting routine showed that several thousand more steps were required to achieve convergence, which we discuss in detail in Appendix D. Nevertheless, the physical model parameters we measured were all consistent with previous works. Of note, we fitted PA = 336°.19 ± 0°.10, a (Gaussian) vertical aspect ratio of h = 0.0157 ± 0.0013 and ωf = 15°.6 ± 0°.6. The Gaussian vertical aspect ratio is slightly lower than the 0.02 assumed in M. A. MacGregor et al. (2017), and consistent with that fitted by G. M. Kennedy (2020).11 Given the resolution of the ALMA data, the low value that we determine for h, and the low h constrained by previous work, we note the vertical structure is plausibly still unresolved. Furthermore, in G. M. Kennedy (2020), the disk rotation was modeled with the longitude of ascending node Ω = 156°.4 (measured anticlockwise from north) which is directed toward the southwest ansa, whereas we model PA = 336°.19 antic- lockwise from north to align with the northeast ansa (i.e., these are simply 180° offset, and thus fully consistent). In addition, the definition of f of G. M. Kennedy (2020) is measured with respect to the longitude of ascending node, whereas the definition we adopt is with respect to south, i.e., independent of the PA (with both again measured anticlockwise). As such, the value of = ±41 1f reported in G. M. Kennedy (2020) would have been measured as = +f f , i.e., 17°.1 ± 1°.0 via the model setup we describe here. In this work, we consistently fit values within a few degrees of this value with similar uncertainties, and thus we deem these to likewise agree. To reduce the number of parameters we fit to, we proceeded by fixing ΔRA and ΔDec, and fitted only for the disk mass (Mdisk), semimajor axis (a0), disk width (wr), (mean) forced eccentricity (ef,0), argument of forced eccentricity (ωf), the power-law index for the forced eccentricity gradient (npow), the inclination (i), the PA, and the (Gaussian) vertical aspect ratio (h). We term this the “General” model. We present in Appendix C the posterior distribution of the converged emcee chains, and discuss its convergence criteria in Appendix D after rerunning the fitting routine for 6000 steps. We find a significant determination of the forced eccentricity gradient power-law index (npow = −1.75 ± 0.16) and parameters for the disk mass, semimajor axis, width, eccentricity, argument of pericenter, inclination, PA, and vertical aspect ratio that are all consistent with M. A. MacGregor et al. (2017) and (where presented) G. M. Kennedy (2020). We assess the properties of our residual emission map to investigate if our model is a good fit. For the complete masked region, we measure a pixel variance equal to the square of our measured σimage, a pixel mean equal to 0.1× σimage, and a distribution that follows a Gaussian-distribution. These values imply that only minimal excess emission remains in the final residual image (which may be fully explained by the presence of bright point sources, which may be background sources) and that it varies as would be expected if it is populated only by noise. The residual map as shown in Figure 2 appears to demonstrate the General model is a very good fit to the data. For example, in the ansae the broader and fainter pericenter are fitted within the noise, and in the apocenter the peak is reproduced within the noise. In the residual maps at the two ansae, we see only small (sub-beam sized) residuals that exceed 3σ (none of which ≳4σ). Along the minor axes, the model appears to be a worse fit, given the model is fainter than the data in those locations. Some of these larger peaks are either background sources (i.e., GDC, as noted by G. M. Ken- nedy et al. 2023) or noisier regions of the data. The residual emission coincident with the apocenter outer edge appears to be slightly oversubtracted by our model, which suggests that our width may not be quite narrow enough in this location.12 These residual features may therefore suggest that a steeper eccentricity gradient model may be needed, but the existing data does not prefer such a model. Similar negative apocenter residuals are present in the residual maps of G. M. Kennedy (2020) and also the fit we performed to the lower-resolution data (discussed in Appendix B). However, we note that the residuals of G. M. Kennedy (2020) host in addition a strong positive residual (located in-between the negative residual regions), i.e., those models cannot account for the apocenter brightness enhancement in the disk. In contrast, we find no positive residuals in the disk ansae in our fits to both the high and low- resolution data. Our findings thus present evidence that the Fomalhaut debris disk can be understood in the paradigm of EVD, with this hosting the first evidence of an eccentricity gradient in its debris belt. We collate our best-fit parameter results in Table 1 for all our fitting setups. We investigate in Section 3.4.2 whether the data prefers a model with free eccentricity and either a constant or negative forced eccen- tricity profile. 3.4.2. Sensitivity to Free Eccentricity? To test the possibility that the power-law eccentricity distribution model fit is biased by not accounting for the “free” (or “proper”) eccentricity parameter (ep), we produce a disk model which includes, ep, ef, and a power- law distribution in the forced eccentricity. We solve analytically the orbit equation for the case of a disk with both free and forced eccentricities following e anpow with fixed values of npow = 0, −1 (and leave to future work 11 Although we note that the upper limits on inclination in G. M. Kennedy (2020) were not as stringent as they should have been given the data, we have traced this to a coding error that resulted in the inclination parameters only affecting radial structure (i.e., the disk models were flat). With the vertical structure parameterization correctly included, we reran those models and instead yield an upper limit on the Gaussian vertical aspect ratio of h < 0.03, with only very small differences for other parameters (e.g., ef = 0.126 ± 0.001 compared to 0.125 ± 0.001 previously). 12 To investigate this, we attempted a more complex radial profile fit to the data, allowing distinct wr values internal/external to the semimajor axis in case these residuals resulted from fitting a Gaussian surface density profile to data that may not be Gaussian in emission. While the fit yielded improved residuals (with inner-edge and outer-edge widths of wr,in = 11.2 ± 0.7 au and wr,in = 15.7 ± 0.7 au, respectively, a smaller mean semimajor axis of a = 137.6 ± 0.4 au, and otherwise statistically consistent best-fit model parameters), this model neither fully removed these negative residual features, nor was the improvement in χ2 statistically significant versus the simpler Gaussian model. 5 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. the general case for any value of npow, i.e., for =e ( ) [ ( )] [ ( ])/ +e a a i e iexp expf n f p p,0 0 pow , for a sys- tem with (mean) forced and free eccentricities (ef,0 and ep, respectively) and arguments of forced and free eccentricities (ωf and ωp, respectively). To model the free eccentricity, we adopt the same approach as J. B. Lovell & E. M. Lynch (2023), i.e., we linearly sum independent “families” of orbits (each comprising a single argument of free eccentricity ωp) uniformly over the range of 0� ωp < 2π, where each family comprises a disk withMdisk, a0, wr, ef,0, ωf, ep, and (fixed) npow. J. B. Lovell & E. M. Lynch (2023) demonstrated that this parametric modeling approach for the free eccentricity distribution yields consistent results to methods that instead model individual particles (e.g., M. A. MacGregor et al. 2017 and G. M. Kennedy 2020, where npow = 0 was implicitly adopted). Given the negative best-fit npow that we earlier inferred (and that previous modeling efforts implicitly had npow = 0), we consider here only the two cases of npow = 0 and npow = −1. First, by fitting a six-parameter model (i.e., with npow fixed to 0) and ep allowed to freely vary, we run the same model fitting procedure described earlier with emcee (for 3000 steps) and find that the best-fit solution results in ep = 3.85% ± 0.17%, i.e., this requires a significant free eccentricity. We term this model “Proper 1.” Second, rerunning this process instead with npow fixed to −1 (for 5000 steps to ensure convergence), we find that the best-fit solution results in ep = 1.9% ± 0.9% (for which we place an upper limit on ep < 3.59% based on the 99.7th-percentile of all chains after burn-in). We term this model “Proper 2.” We present the model and residuals of these fitting processes in Figure 3, all best-fit values in Table 1, and the corner plots for the posteriors in Figures 9 and 10. One can see that while the npow = −1 provides a reasonable fit (albeit poorer than the model with a steeper eccentricity gradient and no free eccentricity), the model with npow = 0 and a significant free eccentricity is poorer at interpreting the broader width at pericenter and the brightness enhancement at apocenter consistent with the findings in J. Chittidi et al. (2025), and also those of M. A. MacGregor et al. (2017) and G. M. Kennedy (2020). We note that by definition there is a degeneracy between ep and wr, which results in the wide range of wr posterior distributions. 3.4.3. Which Model is Preferred? In Table 1, we define Δχ2 as the relative difference in the best-fit value χ2 of the “General” model versus the models with free eccentricities, where the General model returns a lower χ2 than the Proper 1 and 2 models by 70 and 22, respectively. We quantify these differences statistically by estimating the associated probability that (given each of the models has nine free parameters) the relative difference in χ2 is significantly improved by calculating the Incomplete Gamma function for each value of Δχ2, i.e., ( ) ( ) / =P N t tdt 1 2 exp . 7 N par 0 2 1 2 par In terms of confidence intervals, these respective Δχ2 values therefore imply that the General model is favored over the Proper 1 and 2 models by 6.7σ and 2.6σ, respectively, and that the Proper 2 model is favored over Proper 1 by 5.2σ. We note that in a third version of the Proper model which instead had a fixed eccentricity power-law gradient of −2, we found the General model still had a lower χ2 by 10.6, and thus performed statistically as well, with a preference of only 1.0σ. Overall, these differences imply that while the disk may host some small amount of free eccentricity, the eccentric disk is dominated by the presence of a forced eccentricity distribution which falls with semimajor axis, given the strong preference of the General and Proper 2 models over the Proper 1 model, and further the preference of the General model over the Proper 2 model. Although fitting models to the data in the visibility domain requires an order of magnitude longer to converge, we note that the overall result we present is robust to the methodology applied, and a visibility-domain analysis would conclude similarly. We verified this by producing from each of our three best-fit General, Proper 1 and Proper 2 models, visibility tables (gridded to the same UV-spatial frequencies as the data) from the Cycle 5, Apocentre-only pointing which contains the very highest SNR data (see J. Chittidi et al. 2025). By subtracting these visibilities from the data in the visibility domain (after applying the CASA statwt function to the data), we measured their relative χ2 values (directly from the visibilities) and imaged the residuals. Via this method, the residuals strongly favor the General model, which leaves comparable Table 1 Best-fit Model Results for the Fits to the High-resolution ALMA Data (top) and the Best-fit “Verify” Model Fit to the Lower-resolution Data (Bottom) Parameter Units General Proper 1 Proper 2 Verify Mdust 10−2M⊕ 18.70 ± 0.20 20.90 ± 1.3 21.70 ± 1.3 19.90 ± 0.20 a0 au 138.79 ± 0.12 138.93 ± 0.12 138.84 ± 0.12 138.89 ± 0.10 wr au 13.51 ± 0.29 9.3 ± 0.6 12.8 ± 1.0 15.3 ± 0.3 ef,0 % 12.56 ± 0.12 12.56 ± 0.11 12.57 ± 0.12 13.46 ± 0.14 ep % [0.0] 3.89 ± 0.16 <3.6 [0.0] ωf deg 15.2 ± 0.6 15.0 ± 0.6 15.2 ± 0.6 19.1 ± 0.6 npow ⋯ −1.75 ± 0.16 [0.0] [−1.0] −1.16 ± 0.16 h % 1.57 ± 0.13 1.56 ± 0.14 1.50 ± 0.14 [1.43] i deg 66.44 ± 0.09 66.44 ± 0.09 66.43 ± 0.09 [66.6] PA deg 336.19 ± 0.06 336.22 ± 0.05 336.20 ± 0.06 [336.4] Δχ2 ⋯ 0.0 +70.0 +21.5 ⋯ Note. The Δχ2 values are all with reference to the General model. Values presented in square brackets were fixed to the presented value during fitting. For completeness we also present the parameters determined for the verification model which only used the low-resolution data of G. M. Kennedy (2020). 6 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. residuals to those seen in Figure 2 (i.e., no residual artifacts above 4σ) and further has the lowest χ2 by 350 (versus Proper 1) and 245 (versus Proper 2). 3.5. Conclusion: An Eccentric Gradient in the Fomalhaut Debris Disk We find clear evidence that there is a gradient in the eccentricity distribution of Fomalhaut’s planetesimals, and thus in its debris dust as observed with ALMA. This finding suggests that the Fomalhaut debris disk is being shaped by the EVD effect, as discussed in Section 2. Most intriguing about the inferred power-law gradient value of the General model (npow = −1.75 ± 0.16) is its near-consistency with the classical expectation for secular interactions between an internal-planetary perturber with distantly separated outer planetesimals (for which npow = −1; see e.g., C. D. Murray & S. F. Dermott 1999, i.e., the same value that we fixed in the Proper 2 model, and fitted to the archival ALMA data). The value fitted to the newest data nevertheless appears steeper than this classical expectation (being discordant at the 4.6σ level, though we note this discrepancy disappears when we fit to the single configuration, lower-resolution mosaic of M. A. MacGregor et al. 2017, which is fully consistent with a power-law eccentricity gradient of −1). We defer further discussion to the origin of this eccentricity profile to Section 5. 4. Analysis of the Planetary System Since the inference of an eccentricity gradient in Fomal- haut’s disk has important implications for the presence and orbital properties of an internal planet interacting with the disk, we next investigate planet properties plausible with this interpretation. We first consider plausible single-planet scenarios that could be responsible for carving the disk structure in Section 4, i.e., based on the modeled properties of the disk’s inner edge and inner-edge eccentricity. In Section 5, we consider if either such planetary scenario is preferred and the plausible origins of the disk’s eccentricity. Specifically, in Section 5.1 we conduct N-body modeling of planet–disk interactions consistent with the derived planet scenarios, and in Section 5.2 we analytically investigate whether disk self- gravity (in the presence of an internal planet perturber consistent with the derived planet scenarios) may be responsible for driving a steeper gradient (than npow = −1) in to the disk. -10.00.010.0 -20.0 -10.0 0.0 10.0 20.0 R el at iv e D ec l. O ff se t [a rc se c] −20 0 20 40 60 80 100 120 -10.00.010.0 Relative RA Offset [arcsec] ep, ef ∝ a0 −20 0 20 40 60 80 100 120 -10.00.010.0 In ten sity [µ J y b eam − 1] −40 −30 −20 −10 0 10 20 30 40 -10.00.010.0 -20.0 -10.0 0.0 10.0 20.0 R el at iv e D ec l. O ff se t [a rc se c] −20 0 20 40 60 80 100 120 -10.00.010.0 Relative RA Offset [arcsec] ep, ef ∝ a−1 −20 0 20 40 60 80 100 120 -10.00.010.0 In ten sity [µ J y b eam − 1] −40 −30 −20 −10 0 10 20 30 40 Figure 3. Left: ALMA data. Center: Best-fit models which include ep (on same image scale, with top for the model fixed with npow = 0, and npow = −1 for bottom). Right: Residual emission after we subtract the best-fit models from the ALMA data. Contours are shown at the ±3σ level in the residuals and at the 5σ level for the data/model. Beams are shown in the lower-left corner of each plot. In all plots, north is up and east is left. 7 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. 4.1. Planet Properties Consistent with ALMA and JWST Observations On the assumption that a single, coplanar planet is responsible for carving the observed disk structure in the Fomalhaut debris disk (i.e., its inner-edge location, eccen- tricity, and eccentricity gradient), we consider two scenarios that can constrain such a planet’s properties.13 Recent observations of Fomalhaut with the JWST NIRCam instrument place constraints on planet masses ofMpl ≲ 200M⊕ beyond 5″ (38.3 au) from the stellar center (M. Ygouf et al. 2024). Our analysis of single-planet scenarios are consistent with such a nondetection, though we find the plausible planet mass range well below this NIRCam limit. Perhaps more important for constraining planetary properties are the JWST MIRI observa- tions (A. Gáspár et al. 2023). In these, the first evidence of an “intermediate belt” is presented, which has inner and outer edges of 83 au and 104 au, respectively, with large corresp- onding eccentricities of 0.31 and 0.265. We utilize these findings here, on the assumption that a planet could not reside within the region spanned by the intermediate belt, but instead must be located either between the ALMA-observed “main belt” and the intermediate belt, or interior to the intermediate belt. We present the single-planet properties (derived in Sections 4.1.1 and 4.1.2) in Table 2. 4.1.1. A Gap-carving Planet? We first discuss the possibility that a single planet fully carved the gap up to the disk’s inner edge directly. In this scenario, the ALMA observations can place constraints on such a planet’s semimajor axis and eccentricity, using Equation (10) of A. J. Mustill & M. C. Wyatt (2012), i.e., ( ) ( )/ /µ= ×a a a e1.8inner pl pl inner 1 5. For a planet with a mass ranging from Mpl ≈ 1–16M⊕ (or μ in the range 1.7 × 10−6 to 2.5 × 10−5 given the 2M⊙ mass of Fomalhaut, the mass range for which we justify below), this expression requires planet semimajor axes in the range apl = 109–115 au and eccentricities in the range epl = 0.18–0.20, with the innermost semimajor axis apl in this range corresponding to the highest values of epl and Mpl. We calculate these values based on the inner-edge disk semimajor axis (of ainner ∼ 125 au) and a disk eccentricity at this location (einner ≈ 0.16, estimated on an extrapolation from the disk center of the fitted values of e and npow). We justify the range of planet masses with the lower-limit corresponding to the timescale on which a given planet can plausibly secularly force its eccentricity into the disk, and at the upper limit based on direct planetesimal clearing. For example, applying Equation (15) of A. J. Mustill & M. C. Wyatt (2009) with the range of derived apl and epl and known stellar mass, semimajor axes out to 150 au require at least a μ ≈ 1.7 × 10−6 body to be secularly forced within the age of the Fomalhaut system. In addition, Equation (10) of A. J. Mustill & M. C. Wyatt (2012) suggests that planets at 109 au require masses of μ ≈ 2.5 × 10−5 to force the inner- edge disk eccentricity. Utilizing S. Morrison & R. Malhotra (2015), i.e., that planets clear material within Δa = 1.2μ0.31apl of their semimajor axis, such a body would clear material inwards to 104 au, i.e., the intermediate belt outer edge, thus we set a mass upper limit based on this location. We refer to this as the “Gap” planet in Table 2. We note that if the intermediate belt is not a long-lived planetesimal belt, planet mass constraints from JWST NIRCam (i.e., Mpl ≲ 200M⊕) would allow planet semimajor axes and eccentricities to span a slightly broader range of 100–120 au and 0.18–0.23, respectively. The planet orbital parameters of this “Gap” planet all agree with those presented for Fomalhaut b (P. Kalas et al. 2008). Although more recent observations and analysis of Fomalhaut b imply this is no longer present (A. Gaspar & G. Rieke 2020; M. Ygouf et al. 2024), this analysis suggests that a planet with properties consistent with Fomalhaut b may be present in the system. 4.1.2. A Resonant-clearing Planet? The presence of Fomalhaut’s intermediate belt presents another plausible scenario in which the gap between the main belt and the intermediate belt was cleared by an internal planet’s resonant interactions. Such a planet would need to be internal to 83 au based on the e= 0.31 intermediate belt inner edge, and have a strong resonance located between the inner edge of the main belt and the outer edge of the intermediate belt. For a planet’s (strongest) 2:1 mean motion resonance to reside within 2 au from the midpoint of these belt locations, these must span semimajor axes of 70–75 au. This semimajor axis range further requires a (narrower) range of planet masses of 7 × 10−6 ≲ μ ≲ 2.5 × 10−5 for these planets to be responsible for the eccentricity at the inner edge of the intermediate belt (i.e., we utilize once more ( ) ( )/ /µ= ×a a a e1.8inner pl pl inner 1 5, with einner = 0.31). If these bodies also follow the power-law eccentricity gradient, then these would have eccentricities of 0.38–0.41. We refer to this as the “Resonant” planet in Table 2. We note that in this resonant-clearing scenario, if such a planet was responsible for driving the eccentricity distribution of the outer belt (which we measured as ef ∝ a−1.7), then it could feasibly drive larger eccentricities in the intermediate belt given these would be much closer in their semimajor axes. By extrapolating the fitted eccentricity power law to semimajor axes of 83–104 au (i.e., those of the intermediate belt), we calculate eccentricities in the range 0.22–0.32, which overlap reasonably well with the values fitted to the inner and outer edges of the JWST observations of 0.31 and 0.265, respectively (A. Gáspár et al. 2023). While not stated explicitly in the work of A. Gáspár et al. (2023), their JWST analysis corroborates the finding that Fomalhaut hosts an eccentricity distribution that falls with semimajor axis (given the decrease in eccentricity from the intermediate belt inner to outer edge). Table 2 Derived Constraints for a Single Planet Responsible Either for Directly Carving the Inner-edge Gap of the Main Belt (“Gap”) or Clearing Via a 2:1 Resonance (“Resonant”) Scenario apl epl Mpl (au) ⋯ (μ) Gap 109–120 0.20–0.23 1.5 × 10−6–2.5 × 10−5 Resonant 70–75 0.38–0.41 7.0 × 10−6–2.5 × 10−5 13 For a discussion of the parameters one may derive for planets orbiting Fomalhaut at inclinations misaligned with the disk, we refer the reader to T. D. Pearce et al. (2021). 8 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. 5. The Origin of Fomalhaut’s Eccentric Disk 5.1. N-body Modeling of the Fomalhaut Debris Disk In this section, we investigate the plausible origins of the disk eccentricity, by considering the stability of the disk over the full 440Myr age of the star (E. E. Mamajek 2012). We base this section on the outcome of our earlier modeling and the plausible (single) planet scenarios consistent with observa- tions. In all cases, we consider a scenario with a single star (Fomalhaut A) orbited by a single (gravitating) eccentric planet, of substantially lower mass. To determine the stability of the Fomalhaut debris disk, we consider the interaction of this star–planet system with a population of test particles (the effects of the disk self-gravity are explored in Section 5.2), which sample the hypothetical birth orbital element distribu- tion of the Fomalhaut main and intermediate belts. 5.1.1. REBOUND Setup To investigate the stability of orbits in the Fomalhaut system, we use REBOUND (H. Rein & S. F. Liu 2012). We adopt the WHFAST symplectic integrator (H. Rein & D. Tamayo 2015) with a fixed timestep set to 1% of the initially innermost particles orbital period,14 using democratic heliocentric coordinates (i.e., where positions are measured relative to the star, and the corresponding canonical momenta are measured relative to the barycentre; H. Rein & D. Tamayo 2019). In total, we consider six scenarios (that investigate three planet setups, and two disk setups) that are summarized in Table 3 and described below. We explore both the “Gap” and “Resonant” planet scenarios (as per Table 2) for either initially circular or initially eccentric test particle disks. For the Resonant planet ,we consider a single mass of Mpl = 2.5 × 10−5M*, a semimajor axis of apl = 70.87 au, and an eccentricity epl = 0.4. Scenarios with a “1” in their ID include this planet. For the Gap planet, we consider two masses of Mpl = 1.0 × 10−6M* and Mpl = 1.0 × 10−7M*, both with a semimajor axis of apl = 112.9 au, and eccentricity epl = 0.215. Scenarios with either a “2” or “3” in their ID include this planet. Adopting units in which GM* = 1 and unit length equal to the planets semimajor axis, each planet has an orbital period ≈2π (strictly the orbital period is slightly shorter due to the finite planet mass). We consider two initial scenarios, either with all particles initially circular, or initially eccentric. In the case of the initially circular test particles, we distribute these with e= 0. Scenarios with a “C” in their ID relate to this initial disk setup. In the case of the initially eccentric test particles, we distribute these with ( )/=e e a apl pl 1 (i.e., consistent with our model- ing). Scenarios with an “E” in their ID relate to this initial disk setup. We setup the disk with 2 × 105 test particles uniformly distributed between ( )/a70 aupl 1 and ( )/a180 aupl 1. We note that this broad radial distribution spans all radii where Fomalhaut’s main belt and intermediate belt have been observed. In all cases, we randomly distribute the particles uniformly in semimajor axis-mean anomaly space. To avoid issues with the planetary point mass singularity, we adopt REBOUND’s gravitational softening with a gravitational soft- ening parameter set to b = 10−6[apl]. This choice of softening radius is well within the planet’s Hill sphere for all scenarios considered. 5.1.2. REBOUND Results We run simulations E1, E2, E3, C2, and C3 for the full 440Myr lifetime of Fomalhaut. We run simulation C1 for just 1 Myr because by this time the disk is already disrupted (in the sense that large sections of the disk have been broken up or cleared by interaction with the planet). We show the outcomes of these six scenarios in Figure 4, first in terms of their “conditional” mass surface density distributions (i.e., the simulation outcomes) and second to aid interpretation because their surface density distributions convolved with a specified radial profile. These convolved distributions convolve the results of the REBOUND simulations, with a mass per unit semimajor axis of ( ) ( )=M a w a a w 2 exp 2 , 8a r r 0 2 2 consistent with the modeling in Section 3. Presenting the data in this manner comes with important caveats. First, this profile choice effectively masks all emissions not imaged by ALMA. Thus, while the simulations can inform us of the stability of orbits in Fomalhaut’s intermediate belt, we restrict our investigation here purely to Fomalhaut’s main belt. By extension, this further downplays the importance of the planet’s gap carving since this hides the fact that lower-mass planets are unable to carve deep gaps in the particle Table 3 All Six REBOUND Scenarios Simulated in this Work, Starting with Disk Properties (Disk Extent and Eccentricity) and Planet Properties (Planet-stellar Mass Ratio and Semimajor Axis), the Stability of the Disk, and Comments on Any Observed Features Present ID Disk Semimajor Axis Extent Disk Eccentricity Profile Planet Mass Ratio Planet Semi- major Axis Planet Eccentricity Comments (au) (au) C1 70–180 e = 0 2.5 × 10−5 70.87 0.4 Disk disrupted C2 70–180 e = 0 1.0 × 10−6 112.9 0.215 Disk disrupted C3 70–180 e = 0 1.0 × 10−7 112.9 0.215 Disk disrupted E1 70–180 e ∝ a−1 2.5 × 10−5 70.87 0.4 Gap carved, gap in main belt, spirals E2 70–180 e ∝ a−1 1.0 × 10−6 112.9 0.215 Gap carved, weak spirals E3 70–180 e ∝ a−1 1.0 × 10−7 112.9 0.215 Gap carved, but very shallow 14 Such a timestep choice is consistent with REBOUND simulations of debris disks in the literature, e.g., J. Dong et al. (2020) and T. D. Pearce et al. (2024), and consistent with the advice provided in the REBOUND API Documentation; see comment on “Timestepping” in https://rebound.readthedocs.io/en/latest/ simulationvariables/. 9 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. https://rebound.readthedocs.io/en/latest/simulationvariables/ https://rebound.readthedocs.io/en/latest/simulationvariables/ distributions (one should thus consider both plots per scenario). Second, due to our ignorance about the initial conditions of the Fomalhaut disk, the chosen mass per orbit is, by necessity, arbitrary. The mass per orbit is not affected by secular perturbations, and so, away from resonances and close encounters, the planet has no way to modify the mass on orbits. Overall, the REBOUND simulations are thus not expected to reproduce either the ALMA and JWST observa- tions, but instead inform us about the stability of the orbits within the belt and identify features (e.g., gaps or spiral density waves) that may be expected to be observed for given planetary properties. With these caveats in mind, we find (in all cases) that scenarios initialized with circular disks do not produce stable eccentric belts. In the three cases C1, C2, and C3, we find that these all end their simulations with disrupted particle distributions (inconsistent with observations). In contrast, we find that scenarios initialized with eccentric disks (E1, E2, and E3) survive the full 440Myr simulation timescale, and result in eccentric disks consistent with those observed for Fomalhaut. In addition, both the Resonant and Gap planet scenarios carve gaps in the location between the main and intermediate belts (in the lowest mass planet case, this gap is relatively shallow in comparison to the more massive planet scenarios). One intriguing feature in the E1 scenario is the formation of a spiral density wave (due to the massive Resonant planet’s interaction with the disk). While no such features have been observed in the Fomalhaut disk as yet, deeper, higher-resolution observations could be used to investigate the plausibility of this scenario and constrain the mass of a planet with these orbital parameters. Overall, the REBOUND simulations lend weight to two important conclusions. First, the simulations suggest that Fomalhaut’s disk may need to have been formed as an eccentric disk in order to survive ∼440Myr evolution (i.e., formation within the protoplanetary disk). Second, and by extension, while a single planet could be responsible for carving the inner edge of the Fomalhaut main belt (or the gap between the intermediate and main belts), this same planet may not be capable of forcing its own eccentricity into an initially circular planetesimal belt via secular planet–disk interactions to yield an eccentric disk with Fomalhaut’s properties. In the case of the Gap planet, this body can 0 50 100 150 200 R ad iu s [a u ] Disk at t = 0: e=0 Planet: apl=71au, µ=2.5×10−5, epl=0.40 C1 Disk at t = 0: e=0 Planet: apl=113au, µ=10−6, epl=0.215 C2 Disk at t = 0: e=0 Planet: apl=113au, µ=10−7, epl=0.215 C on d . su rface M D F [cts au − 2] C3 0e+00 5e+00 1e+01 2e+01 2e+01 2e+01 0 100 200 300 Angle [deg] 0 50 100 150 200 R ad iu s [a u ] Disk at t = 0: e=0 Planet: apl=71au, µ=2.5×10−5, epl=0.40 0 100 200 300 Angle [deg] Disk at t = 0: e=0 Planet: apl=113au, µ=10−6, epl=0.215 0 100 200 300 Angle [deg] Disk at t = 0: e=0 Planet: apl=113au, µ=10−7, epl=0.215 S u rface d en sity [cts au − 2] 0e+00 1e+01 2e+01 3e+01 4e+01 0 50 100 150 200 R ad iu s [a u ] Disk at t = 0: e∝a−1 Planet: apl=71au, µ=2.5×10−5, epl=0.40 E1 Disk at t = 0: e∝a−1 Planet: apl=113au, µ=10−6, epl=0.215 E2 Disk at t = 0: e∝a−1 Planet: apl=113au, µ=10−7, epl=0.215 C on d . su rface M D F [cts au − 2] E3 0e+00 5e+00 1e+01 2e+01 2e+01 2e+01 0 100 200 300 Angle [deg] 0 50 100 150 200 R ad iu s [a u ] Disk at t = 0: e∝a−1 Planet: apl=71au, µ=2.5×10−5, epl=0.40 0 100 200 300 Angle [deg] Disk at t = 0: e∝a−1 Planet: apl=113au, µ=10−6, epl=0.215 0 100 200 300 Angle [deg] Disk at t = 0: e∝a−1 Planet: apl=113au, µ=10−7, epl=0.215 S u rface d en sity [cts au − 2] 0e+00 1e+01 2e+01 3e+01 4e+01 Figure 4. Plots for the surface density maps from the REBOUND simulations in r–f space. Upper two rows: Conditional surface mass density function and absolute surface density simulations with initially circular disk conditions. Lower two rows: Simulation outcomes with initially eccentric disk conditions presented likewise. Specific planet and disk assumptions are presented on each plot, with black-dashed lines showing the planet’s orbit. Initially circular disk conditions do not yield stable disks at 440 Myr, whereas scenarios which are initially eccentric yield stable eccentric disks at 440 Myr. The same information is tabulated in Table 3. 10 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. plausibly clear material from the main belt’s inner edge (and the intermediate belt’s outer edge) if this is situated in an initially eccentric debris belt, i.e., via close planet-planetesimal encounters. In the case of the Resonant planet, this could plausibly clear material in the same location due to the overlap with its 2:1 mean motion resonance and subsequent planete- simal ejections. 5.1.3. Comparison with Eccentric Protoplanetary Disks The conclusion that eccentric disks that survive to the stellar main sequence may need to be produced during the protoplanetary disk phase is consistent with the conclusion of G. M. Kennedy (2020), who argued this from the independent standpoint of the apparent narrowness of the debris disks of Fomalhaut and HD 202628. There are indeed a (small) number of known eccentric protoplanetary disks, such as MWC 758, HD 100546, and IRS 48 (R. Dong et al. 2018; D. Fedele et al. 2021; H. Yang et al. 2023, all of which are hosted by ∼2M⊙ stars), as well as the eccentric circumbinary disk of IRAS 04158+2805 (which is hosted by a lower-mass mid-M SpT binary; S. M. Andrews et al. 2008; E. Ragusa et al. 2021). The eccentric debris disks of HD 202628 and HD 53143 are hosted by ∼Solar-mass stars (V. Faramaz et al. 2019; M. A. MacGregor et al. 2021), and those of HD 38206 and HR 4796 likewise orbit ∼2M⊙ (A-type) stars (J. Olofsson et al. 2019; M. Booth et al. 2021). Whether this bias toward earlier-type, intermediate-mass stars is physical or observa- tional should be determined in future work. A plausible origin of Fomalhaut’s eccentric disk is that the eccentricity could have resulted from the formation of a planet (with properties consistent with those in Table 2) within the protoplanetary disk, whereby initial planet–gas disk interac- tions determined the eccentricity profile of the disk, and the eccentricity of the planet. In such a scenario, an eccentricity profile which decreases with radius is expected (e.g., J. Teys- sandier & G. I. Ogilvie 2016). This finding follows from the fact that eccentricity profiles that are either flat or increasing with radius tend to differentially precess as a result of gas pressure forces. In the case of Fomalhaut, following the dissipation of its protoplanetary disk gas (and any remnant primordial dust), subsequent planet–disk interactions with the eccentric planetesimal belt plausibly sculpted its disk’s inner and outer edges, which over the course of 440Myr resulted in the belt morphology, while maintaining its initial eccentricity distribution (or one very similar). 5.2. Additional Physics? The Role of Disk Self-gravity We have thus far explored the possible origin of Fomalhaut’s disk eccentricity assuming a stationary planet perturbing a massless debris disk. However, additional physical effects, such as the influence of a massive and self- gravitating debris disk, may affect the outcomes. As mentioned in Section 2, a planet within a debris disk induces (secular) forced planetesimal eccentricities such that ef(a) ∝ a−1 (C. D. Murray & S. F. Dermott 1999). This, however, applies specifically to a massless debris disk (Md = 0), i.e., composed of test particles. By contrast, if the disk is massive, its self-gravity can steepen the forced eccentricity profile (A. A. Sefilian 2024): namely, to as much as ef(a) ∝ a−4.5, depending on ap/ain and the mass distribution within the disk (see also, A. A. Sefilian et al. 2021). This occurs when the disk-induced apsidal precession rates of the planet and/or the planetesimals—absent in massless disk models—dominate over the planet-induced precession of planetesimal orbits. For ap ∼ ain, this condition generally translates to Md/mp ≳ 1, and would require Md/mp ≲ 1 for ap ≲ ain. Interestingly, under the same conditions on Md/mp, a similar effect arises in the disk’s vertical aspect ratio h if the planet–disk system is misaligned: rather than remaining constant with radius (as expected for Md = 0), h may instead decrease with distance from the star if Md ≠ 0 (A. A. Sefilian et al. 2025). This effect, however, is not accounted for in this study. Without presuming specific planet masses, ALMA observa- tions imply a total disk mass (i.e., including the largest, colliding planetesimals) consistent with the range inferred by the collisional modeling of A. V. Krivov & M. C. Wyatt (2021), which estimates Md = 1.8–360M⊕. If true, disk gravity could potentially explain the steep ef(a) profile indicated by our modeling. While it may be tempting to test this hypothesis using the results of A. A. Sefilian (2024), their study is a proof-of-concept based on a semianalytical frame- work that accounts for the axisymmetric component of the disk potential, but ignores the nonaxisymmetric counterpart. In principle, this limitation can be addressed using existing orbit- averaged, secular codes (e.g., the linear, N-RING code of A. A. Sefilian et al. 2023 or the more accurate GAUSS code of J. R. Touma et al. 2009) and/or direct N-body codes; however, this is more suited to a separate study (Sefilian A. A. et al. 2025, in preparation). 6. Conclusions We present a new model of Fomalhaut’s eccentric debris disk by fitting parametric debris disk models to archival ALMA observations (with a synthesized physical resolution 4–6 au). The models are developed from those introduced by E. M. Lynch & J. B. Lovell (2021) and J. B. Lovell & E. M. Lynch (2023), whereby a gradient in the forced eccentricity parameter, ( )e a anpow, is introduced to fit the known disk asymmetries: a brighter apocenter versus pericen- ter, and a broader pericenter versus apocenter. We collectively term this physical effect within disks EVD. The best-fit parameter determinations from the modeling find parameters broadly consistent with those from other models, but with the inclusion of a steep, negative eccentricity gradient as per e ∝ a−1.75 ± 0.16. By analyzing the goodness of fit, we deem that such a description provides a novel interpretation of this debris disk, and is statistically preferred to models with constant free and forced eccentricity distributions through the disk. To our knowledge, this is the first reported eccentricity gradient in a debris disk. The value we derive is very close to that expected from classical expectations of (massless) disk– planet interactions ef ∝ a−1, suggesting the gradient may be forced by a planetary perturber. Based on the ALMA modeling, and published data from JWST NIRCam and MIRI, we quantify planetary orbital and mass constraints for two coplanar single-planet scenarios consistent with observations. One scenario describes a 109–115 au planet that has directly cleared material up to the inner edge of Fomalhaut’s “main belt” (as imaged by ALMA). The second posited scenario describes a 70–75 au planet that has cleared the disk’s inner edge at its 2:1 mean motion resonance, with this planet sitting interior to the JWST-imaged 11 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. “intermediate belt.” In both cases, the implied planet mass and semimajor axis ranges are below sensitivity thresholds for existing planet detection methods. We assess the stability of the Fomalhaut debris disk by conducting N-body REBOUND modeling. We set up six sets of initial conditions, based on three planets, and two initial conditions for the disk, with this either being circular or eccentric. We show that only scenarios with an initially circular disk are stable and finalize as eccentric disks after the 440Myr simulation lifetime (matching that of the age of Fomalhaut). These findings may suggest that planet–disk interactions are primarily responsible for sculpting the disk’s morphology (i.e., its inner-edges, and as-per the JWST observations, gaps in the disk), but not its eccentricity, and thus that Fomalhaut’s eccentric ring was plausibly born eccentric. We discuss caveats with these models (e.g., the exclusion of disk self-gravity) and propose further simulations to investigate the robustness of these findings for massive debris disks. Finally, we highlight that the code used to model Fomalhaut’s surface density distribution in this work (Equations (1)–(4)) is available publicly on https://github. com/astroJLovell/eccentricDiskModels. We release this to enable our work to be reproduced easily, and such that others can utilize this to model the surface density profiles of other optically thin disks. Acknowledgments We thank the anonymous referee for their comments and suggestions, which greatly improved the clarity of our investigation. J.B.L. acknowledges the Smithsonian Institute for funding via a Submillimeter Array (SMA) Fellowship, and the North American ALMA Science Center (NAASC) for funding via an ALMA Ambassadorship. J.B.L. dedicates this study to the memory of Henry Chesters, with whom he discussed this with great interest in their final conversation, and for H.C.’s lifelong encouragement of his scientific endeavors. A.A.S. is supported by the Heising-Simons Foundation through a 51 Pegasi b Postdoctoral Fellowship. We acknowledge the operational staff and scientists involved in the collection of data presented here. Facility: ALMA. Software: This research made use of the following software packages: CASA (J. P. McMullin et al. 2007) DS9 (Smithsonian Astrophysical Observatory 2000). emcee (D. Foreman-Mackey et al. 2013); RADMC-3D (C. P. Dullemond et al. 2012); REBOUND (H. Rein & S. F. Liu 2012). Software and Third Party Data Repository Citations This paper makes use of the following ALMA data: ADS/ JAO.ALMA 2013.1.00486.S, ADS/JAO.ALMA 2015.1.00966. S, and ADS/JAO.ALMA 2017.1.01043.S. ALMA is a partner- ship of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Appendix A Data Masking We present in Figure 5 the result of fitting the General model while only masking pixels within a 200 au (projected) distance from the star, i.e., without masking by the primary beam. We show the spatial mask in the left panel of this figure (in green dashed–dotted line), and also show the primary-beam response at the 0.66 level (in orange dashed–dotted lines) which defines the inner boundary of any fits that incorporate a primary-beam mask. The regions of the two ansae between the green and orange dashed–dotted lines define the mask used in the main body of this work. It is evident that along the disk’s minor axis, the primary- beam response is significantly reduced in comparison to the disk ansae, and while the minor axes appear to be slightly better fit in comparison to the General model with a primary-beam mask, the fit is significantly worse at the apocentre, justifying the choice to fit this model only at the two disk ansae. Appendix B Model Verification We verify the modeling and fitting methodology by fitting a six-parameter model to the archival ALMA Cycle 3-only data for Fomalhaut, first presented in M. A. MacGregor et al. (2017) and subsequently in G. M. Kennedy (2020). We refer to this low-resolution model as “Verify.” We find that the parameter distributions and mean values of the fitting are consistent with those previous works (i.e., M. A. MacGregor et al. 2017; G. M. Kennedy 2020), albeit with a slightly larger (mean) forced eccentricity (which is most likely due to the difference in eccentricity profile we fit here). Such consistency provides assurance that the fitting methodology we conduct here provides reasonable uncertainties on the fitted parameters. We present in Figure 6 the posterior distribution after removing burn-in chains, the best-fit values on Table 1, and the best-fit model and residuals in Figure 7. Most interesting, is that we find a negative eccentricity gradient is likewise required with this model to fit even the lower-resolution data, with npow = −1.16 ± 0.16, and that this model setup can account for the disk’s width and brightness asymmetries (presented by G. M. Kennedy 2020). One improvement is that the model we present here fully accounts for the northwest positive emission that is still present in the apocenter residuals of G. M. Kennedy (2020). In verifying this fitting methodology, we thus find that evidence of an eccentricity gradient in Fomalhaut’s disk existed prior to its observations at higher-resolution. -10.00.010.0 -20.0 -10.0 0.0 10.0 20.0 R el at iv e D ec l. O ff se t [a rc se c] −20 0 20 40 60 80 100 120 -10.00.010.0 Relative RA Offset [arcsec] ef ∝ a−1.75 −20 0 20 40 60 80 100 120 -10.00.010.0 In ten sity [µ J y b eam − 1] −40 −30 −20 −10 0 10 20 30 40 Figure 5. Left: ALMA data. Center: Best-fit ef(a) ∝ a−1.75 model (on same image scale). Right: Residual emission after best-fit model subtraction. Contours are shown at the±3σ level in the residuals and at the 5σ level for the data/model. Beams are shown in the lower-left of each plot. In the left and right plots in amber dashed–dotted lines, we overplot the primary-beam response at the 0.66 level (and in addition at the 0.4 level on the right plot). We show in the left plot in the green dashed–dotted line the 200 au spatial mask discussed in Appendix A. 12 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. https://github.com/astroJLovell/eccentricDiskModels https://github.com/astroJLovell/eccentricDiskModels Mdust = 0.0199+0.0002 0.0002 13 8.6 13 8.8 13 9.0 13 9.2 a 0 a0 = 138.8873+0.0983 0.0993 14 .4 15 .0 15 .6 16 .2 16 .8 w r wr = 15.2977+0.3198 0.3140 0.1 30 0 0.1 32 5 0.1 35 0 0.1 37 5 0.1 40 0 e f ,0 ef, 0 = 0.1346+0.0015 0.0014 18 .0 19 .5 21 .0 f f = 19.1031+0.6276 0.6589 0.0 19 2 0.0 19 5 0.0 19 8 0.0 20 1 0.0 20 4 Mdust 1.5 0 1.2 5 1.0 0 0.7 5 n p ow 13 8.6 13 8.8 13 9.0 13 9.2 a0 14 .4 15 .0 15 .6 16 .2 16 .8 wr 0.1 30 0 0.1 32 5 0.1 35 0 0.1 37 5 0.1 40 0 ef, 0 18 .0 19 .5 21 .0 f 1.5 0 1.2 5 1.0 0 0.7 5 npow npow = 1.1651+0.1613 0.1572 Figure 6. Corner plots for the emcee chains (minus 1000 steps covering burn-in) showing the six parameters fitted to the low-resolution data of M. A. MacGregor et al. (2017) and G. M. Kennedy (2020), i.e., “Verify” model. We find comparable parameter uncertainties to those presented in G. M. Kennedy (2020) and a lower value of the rms in the residual map than the “uniform” model presented by G. M. Kennedy (2020). 13 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. Appendix C Residual Emission Maps and emcee Posterior Distributions Here, we present for the six-parameter (low-resolution) “Verify” model, the corner plots in 6 and the data, model, and residual maps in Figure 7. We present the corner plots for the General, and Proper 1 and Proper 2 models in Figures 8, 9, and 10, and the data, model, and residual maps for the Proper 1 and Proper 2 models in Figure 3. In all cases, the first 1000 steps were removed from the emcee chains (far in excess of the number of burn-in steps). -10.00.010.0 -20.0 -10.0 0.0 10.0 20.0 R el at iv e D ec l. O ff se t [a rc se c] −50 0 50 100 150 200 250 300 350 400 -10.00.010.0 Relative RA Offset [arcsec] ef ∝ a−1.16 −50 0 50 100 150 200 250 300 350 400 -10.00.010.0 In ten sity [µ J y b eam − 1] −40 −30 −20 −10 0 10 20 30 40 Figure 7. Left: ALMA data as presented in M. A. MacGregor et al. (2017) and G. M. Kennedy (2020). Center: Best-fit model for this data. Right: Residual emission (data minus best-fit model). Contours are shown at the ±3σ level in the residuals, and at the 5σ level for the data/model. Beams are shown in the lower-left corner of each plot. In all plots, north is up and east is left. The strong peak at (−4″, −3″) is the source GDC, as noted in both A. Gáspár et al. (2023) and G. M. Kennedy et al. (2023). 14 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. Mdust = 0.0187+0.0002 0.0002 13 8.4 13 8.6 13 8.8 13 9.0 13 9.2 a 0 a0 = 138.7938+0.1209 0.1193 12 .5 13 .0 13 .5 14 .0 14 .5 w r wr = 13.5075+0.2840 0.2927 0.1 22 5 0.1 25 0 0.1 27 5 0.1 30 0 e f ,0 ef, 0 = 0.1256+0.0012 0.0012 13 .5 15 .0 16 .5 18 .0 f f = 15.2347+0.6406 0.6442 2.4 2.1 1.8 1.5 1.2 n p ow npow = 1.7458+0.1577 0.1567 0.0 10 0 0.0 12 5 0.0 15 0 0.0 17 5 0.0 20 0 h h = 0.0157+0.0014 0.0013 33 5.9 5 33 6.1 0 33 6.2 5 33 6.4 0 PA PA = 336.1927+0.0568 0.0577 0.0 18 0 0.0 18 4 0.0 18 8 0.0 19 2 0.0 19 6 Mdust 66 .2 66 .4 66 .6 66 .8 in c 13 8.4 13 8.6 13 8.8 13 9.0 13 9.2 a0 12 .5 13 .0 13 .5 14 .0 14 .5 wr 0.1 22 5 0.1 25 0 0.1 27 5 0.1 30 0 ef, 0 13 .5 15 .0 16 .5 18 .0 f 2.4 2.1 1.8 1.5 1.2 npow 0.0 10 0 0.0 12 5 0.0 15 0 0.0 17 5 0.0 20 0 h 33 5.9 5 33 6.1 0 33 6.2 5 33 6.4 0 PA 66 .2 66 .4 66 .6 66 .8 inc inc = 66.4358+0.0916 0.0898 Figure 8. Corner plots for the emcee chains (minus burn-in steps) showing the nine parameters fitted in this study for the General model. All show the well-behaved features of Gaussian distributions, which are either circular or with little degeneracy (e.g., between fM,dust and wr, and between e and ωf). 15 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. Mdust = 0.0209+0.0015 0.0013 13 8.5 0 13 8.7 5 13 9.0 0 13 9.2 5 r 0 r0 = 138.9286+0.1176 0.1182 8 9 10 11 w r wr = 9.3401+0.5830 0.5634 0.1 22 0.1 24 0.1 26 0.1 28 0.1 30 e e = 0.1256+0.0011 0.0011 13 .5 15 .0 16 .5 f f = 14.9846+0.6473 0.6317 0.0 28 0.0 32 0.0 36 0.0 40 0.0 44 e p ep = 0.0389+0.0015 0.0017 0.0 10 0 0.0 12 5 0.0 15 0 0.0 17 5 0.0 20 0 h h = 0.0156+0.0014 0.0014 33 6.1 33 6.2 33 6.3 33 6.4 PA PA = 336.2189+0.0541 0.0529 0.0 17 5 0.0 20 0 0.0 22 5 0.0 25 0 0.0 27 5 Mdust 66 .15 66 .30 66 .45 66 .60 66 .75 i 13 8.5 0 13 8.7 5 13 9.0 0 13 9.2 5 r0 8 9 10 11 wr 0.1 22 0.1 24 0.1 26 0.1 28 0.1 30 e 13 .5 15 .0 16 .5 f 0.0 28 0.0 32 0.0 36 0.0 40 0.0 44 ep 0.0 10 0 0.0 12 5 0.0 15 0 0.0 17 5 0.0 20 0 h 33 6.1 33 6.2 33 6.3 33 6.4 PA 66 .15 66 .30 66 .45 66 .60 66 .75 i i = 66.4434+0.0899 0.0910 Figure 9. Corner plots for the emcee chains (minus burn-in steps) showing the nine parameters fitted for the Proper 1 model. All show the well-behaved features of Gaussian distributions, which are either circular or with little degeneracy (e.g., between fM,dust and wr, and between e and ωf). 16 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. Appendix D Model Convergence We find by fitting initial rounds of emcee (for the nine free parameter General disk model) that the burn-in stage is approximately a few hundred steps and that the autocorrelation time (τ) is≈120 steps. For the General and Proper 1 models, we measure the autocorrelation function (τ) versus step length for all chains, and find this value asymptotes (with no deviation by more than 1%) to a parameter-mean τ of ≈120. In the case of the Proper 2 model, we measure a mean τ of ≈150, and for the Verify model a mean τ of ≈60. As stated in the emcee guidance (D. Foreman-Mackey et al. 2013), convergence is achieved when the number of steps exceeds >50 × τ. As such, we run the Verify, General, and Proper 1 and Proper 2 models for 3000, 6000, 6000, and 7500 steps, respectively, and verify that convergence is indeed met in all cases. In all cases, we find after 300–500 steps that the chains are burned in. To derive the best-fit posterior distributions presented in Table 1, we removed the first half of the total steps from all chains. We used 40 walkers in all emcee runs, which is far in excess of the necessary 2× the number of free parameters i.e., 12 for the Verify model and 18 for the other three models. ORCID iDs Joshua B. Lovellaa https://orcid.org/0000-0002-4248-5443 Jay Chittidiaa https://orcid.org/0000-0002-4985-028X Mdust = 0.0217+0.0015 0.0013 13 8.5 0 13 8.7 5 13 9.0 0 13 9.2 5 r 0 r0 = 138.8482+0.1205 0.1172 10 .5 12 .0 13 .5 w r wr = 12.7650+0.6997 0.9772 0.1 22 0.1 24 0.1 26 0.1 28 0.1 30 e e = 0.1257+0.0011 0.0012 13 .5 15 .0 16 .5 f f = 15.1589+0.6316 0.6213 0.0 08 0.0 16 0.0 24 0.0 32 e p ep = 0.0202+0.0079 0.0116 0.0 10 0 0.0 12 5 0.0 15 0 0.0 17 5 h h = 0.0150+0.0013 0.0014 33 6.0 33 6.1 33 6.2 33 6.3 33 6.4 PA PA = 336.2032+0.0575 0.0542 0.0 18 0.0 21 0.0 24 0.0 27 0.0 30 Mdust 66 .2 66 .4 66 .6 66 .8 i 13 8.5 0 13 8.7 5 13 9.0 0 13 9.2 5 r0 10 .5 12 .0 13 .5 wr 0.1 22 0.1 24 0.1 26 0.1 28 0.1 30 e 13 .5 15 .0 16 .5 f 0.0 08 0.0 16 0.0 24 0.0 32 ep 0.0 10 0 0.0 12 5 0.0 15 0 0.0 17 5 h 33 6.0 33 6.1 33 6.2 33 6.3 33 6.4 PA 66 .2 66 .4 66 .6 66 .8 i i = 66.4255+0.0907 0.0899 Figure 10. Corner plots for the emcee chains (minus burn-in steps) showing the nine parameters fitted for the Proper 2 model. All show the well-behaved features of Gaussian distributions, which are either circular or with little degeneracy (e.g., between fM,dust and wr, and between e and ωf). 17 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. https://orcid.org/0000-0002-4248-5443 https://orcid.org/0000-0002-4985-028X Antranik A. Sefilianaa https://orcid.org/0000-0003- 4623-1165 Sean M. Andrewsaa https://orcid.org/0000-0003-2253-2270 Grant M. Kennedyaa https://orcid.org/0000-0001-6831-7547 Meredith MacGregoraa https://orcid.org/0000-0001- 7891-8143 David J. Wilneraa https://orcid.org/0000-0003-1526-7587 Mark C. Wyattaa https://orcid.org/0000-0001-9064-5598 References Acke, B., Min, M., Dominik, C., et al. 2012, A&A, 540, A125 Allen, R. H. 1963, Star Names. Their Lore and Meaning (New York: Dover) Andrews, S. M., Liu, M. C., Williams, J. P., & Allers, K. N. 2008, ApJ, 685, 1039 Aumann, H. H. 1985, PASP, 97, 885 Boley, A. C., Payne, M. J., Corder, S., et al. 2012, ApJL, 750, L21 Booth, M., Schulz, M., Krivov, A. V., et al. 2021, MNRAS, 500, 1604 Chittidi, J., MacGregor, M. A., Lovell, J. B., et al. 2025, ApJL, 990, L40 Dawson, R. I., & Murray-Clay, R. 2012, ApJ, 750, 43 Dohnanyi, J. S. 1969, JGR, 74, 2531 Dong, J., Dawson, R. I., Shannon, A., & Morrison, S. 2020, ApJ, 889, 47 Dong, R., Liu, S.-y., Eisner, J., et al. 2018, ApJ, 860, 124 Dullemond, C. P., Juhasz, A., Pohl, A., et al., 2012 Astrophysics Source Code Library, ascl:1202.015 Faramaz, V., Krist, J., Stapelfeldt, K. R., et al. 2019, AJ, 158, 162 Fedele, D., Toci, C., Maud, L., & Lodato, G. 2021, A&A, 651, A90 Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 Gaspar, A., & Rieke, G. 2020, PNAS, 117, 9712 Gáspár, A., Wolff, S. G., Rieke, G. H., et al. 2023, NatAs, 7, 790 Goodman, J., & Weare, J. 2010, CAMCS, 5, 65 Holland, W. S., Greaves, J. S., Dent, W. R. F., et al. 2003, ApJ, 582, 1141 Hughes, A. M., Duchêne, G., & Matthews, B. C. 2018, ARA&A, 56, 541 Kalas, P., Graham, J. R., Chiang, E., et al. 2008, Sci, 322, 1345 Kalas, P., Graham, J. R., & Clampin, M. 2005, Natur, 435, 1067 Kennedy, G. M. 2020, RSOS, 7, 200063 Kennedy, G. M., Lovell, J. B., Kalas, P., & Fitzgerald, M. P. 2023, MNRAS, 524, 2698 Krivov, A. V., & Wyatt, M. C. 2021, MNRAS, 500, 718 Kurucz, R. L. 1979, ApJS, 40, 1 Lovell, J. B., & Lynch, E. M. 2023, MNRAS, 525, L36 Lovell, J. B., Marino, S., Wyatt, M. C., et al. 2021, MNRAS, 506, 1978 Lovell, J. B., Wyatt, M. C., Kalas, P., et al. 2022, MNRAS, 517, 2546 Lynch, E. M., & Lovell, J. B. 2021, MNRAS, 510, 2538 MacGregor, M. A., Matrà, L., Kalas, P., et al. 2017, ApJ, 842, 8 MacGregor, M. A., Weinberger, A. J., Loyd, R. O. P., et al. 2021, ApJL, 911, L25 Mamajek, E. E. 2012, ApJL, 754, L20 Marino, S. 2022, arXiv:2202.03053 Matrà, L., MacGregor, M. A., Kalas, P., et al. 2017, ApJ, 842, 9 Matrà, L., Panić, O., Wyatt, M. C., & Dent, W. R. F. 2015, MNRAS, 447, 3936 Matthews, B. C., Krivov, A. V., Wyatt, M. C., Bryden, G., & Eiroa, C. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 521 McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127 Morrison, S., & Malhotra, R. 2015, ApJ, 799, 41 Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge Univ. Press) Mustill, A. J., & Wyatt, M. C. 2009, MNRAS, 399, 1403 Mustill, A. J., & Wyatt, M. C. 2012, MNRAS, 419, 3074 Olofsson, J., Milli, J., Thébault, P., et al. 2019, A&A, 630, A142 Pan, M., Nesvold, E. R., & Kuchner, M. J. 2016, ApJ, 832, 81 Pearce, T. D., Beust, H., Faramaz, V., et al. 2021, MNRAS, 503, 4767 Pearce, T. D., Krivov, A. V., Sefilian, A. A., et al. 2024, MNRAS, 527, 3876 Ragusa, E., Fasano, D., Toci, C., et al. 2021, MNRAS, 507, 1157 Rein, H., & Liu, S. F. 2012, A&A, 537, A128 Rein, H., & Tamayo, D. 2015, MNRAS, 452, 376 Rein, H., & Tamayo, D. 2019, RNAAS, 3, 16 Sefilian, A. A. 2024, ApJ, 966, 140 Sefilian, A. A., Kratter, K. M., Wyatt, M. C., et al. 2025, arXiv:2505.09578 Sefilian, A. A., Rafikov, R. R., & Wyatt, M. C. 2021, ApJ, 910, 13 Sefilian, A. A., Rafikov, R. R., & Wyatt, M. C. 2023, ApJ, 954, 100 Smithsonian Astrophysical Observatory, 2000 SAOImage DS9: A Utility for Displaying Astronomical Images in the X11 Window Environment, Astrophysics Source Code Library, ascl:0003.002 Stapelfeldt, K. R., Holmes, E. K., Chen, C., et al. 2004, ApJS, 154, 458 Teyssandier, J., & Ogilvie, G. I. 2016, MNRAS, 458, 3221 Touma, J. R., Tremaine, S., & Kazandjian, M. V. 2009, MNRAS, 394, 1085 Wyatt, M. C. 2008, ARA&A, 46, 339 Wyatt, M. C., Dermott, S. F., Telesco, C. M., et al. 1999, ApJ, 527, 918 Yang, H., Fernández-López, M., Li, Z.-Y., et al. 2023, ApJL, 948, L2 Ygouf, M., Beichman, C. A., Llop-Sayson, J., et al. 2024, AJ, 167, 26 18 The Astrophysical Journal, 990:145 (18pp), 2025 September 10 Lovell et al. https://orcid.org/0000-0003-4623-1165 https://orcid.org/0000-0003-4623-1165 https://orcid.org/0000-0003-2253-2270 https://orcid.org/0000-0001-6831-7547 https://orcid.org/0000-0001-7891-8143 https://orcid.org/0000-0001-7891-8143 https://orcid.org/0000-0003-1526-7587 https://orcid.org/0000-0001-9064-5598 https://doi.org/10.1051/0004-6361/201118581 https://ui.adsabs.harvard.edu/abs/2012A&A...540A.125A/abstract https://doi.org/10.1086/591417 https://ui.adsabs.harvard.edu/abs/2008ApJ...685.1039A/abstract https://ui.adsabs.harvard.edu/abs/2008ApJ...685.1039A/abstract https://doi.org/10.1086/131620 https://ui.adsabs.harvard.edu/abs/1985PASP...97..885A/abstract https://doi.org/10.1088/2041-8205/750/1/L21 https://ui.adsabs.harvard.edu/abs/2012ApJ...750L..21B/abstract https://doi.org/10.1093/mnras/staa3362 https://ui.adsabs.harvard.edu/abs/2021MNRAS.500.1604B/abstract https://doi.org/10.3847/2041-8213/adfadb https://doi.org/10.1088/0004-637X/750/1/43 https://ui.adsabs.harvard.edu/abs/2012ApJ...750...43D/abstract https://doi.org/10.1029/JB074i010p02531 https://ui.adsabs.harvard.edu/abs/1969JGR....74.2531D/abstract https://doi.org/10.3847/1538-4357/ab64f7 https://ui.adsabs.harvard.edu/abs/2020ApJ...889...47D/abstract https://doi.org/10.3847/1538-4357/aac6cb https://ui.adsabs.harvard.edu/abs/2018ApJ...860..124D/abstract http://www.ascl.net/1202.015 https://doi.org/10.3847/1538-3881/ab3ec1 https://ui.adsabs.harvard.edu/abs/2019AJ....158..162F/abstract https://doi.org/10.1051/0004-6361/202141278 https://ui.adsabs.harvard.edu/abs/2021A&A...651A..90F/abstract https://doi.org/10.1086/670067 https://ui.adsabs.harvard.edu/abs/2013PASP..125..306F/abstract https://ui.adsabs.harvard.edu/abs/2013PASP..125..306F/abstract https://doi.org/10.1073/pnas.1912506117 https://ui.adsabs.harvard.edu/abs/2020PNAS..117.9712G/abstract https://doi.org/10.1038/s41550-023-01962-6 https://ui.adsabs.harvard.edu/abs/2023NatAs...7..790G/abstract https://doi.org/10.2140/camcos.2010.5.65 https://ui.adsabs.harvard.edu/abs/2010CAMCS...5...65G/abstract https://doi.org/10.1086/344819 https://ui.adsabs.harvard.edu/abs/2003ApJ...582.1141H/abstract https://doi.org/10.1146/annurev-astro-081817-052035 https://ui.adsabs.harvard.edu/abs/2018ARA&A..56..541H/abstract https://doi.org/10.1126/science.1166609 https://ui.adsabs.harvard.edu/abs/2008Sci...322.1345K/abstract https://doi.org/10.1038/nature03601 https://ui.adsabs.harvard.edu/abs/2005Natur.435.1067K/abstract https://doi.org/10.1098/rsos.200063 https://ui.adsabs.harvard.edu/abs/2020RSOS....700063K/abstract https://doi.org/10.1093/mnras/stad2058 https://ui.adsabs.harvard.edu/abs/2023MNRAS.524.2698K/abstract https://ui.adsabs.harvard.edu/abs/2023MNRAS.524.2698K/abstract https://doi.org/10.1093/mnras/staa2385 https://ui.adsabs.harvard.edu/abs/2021MNRAS.500..718K/abstract https://doi.org/10.1086/190589 https://ui.adsabs.harvard.edu/abs/1979ApJS...40....1K/abstract https://doi.org/10.1093/mnrasl/slad083 https://ui.adsabs.harvard.edu/abs/2023MNRAS.525L..36L/abstract https://doi.org/10.1093/mnras/stab1678 https://ui.adsabs.harvard.edu/abs/2021MNRAS.506.1978L/abstract https://doi.org/10.1093/mnras/stac2782 https://ui.adsabs.harvard.edu/abs/2022MNRAS.517.2546L/abstract https://doi.org/10.1093/mnras/stab3566 https://ui.adsabs.harvard.edu/abs/2022MNRAS.510.2538L/abstract https://doi.org/10.3847/1538-4357/aa71ae https://ui.adsabs.harvard.edu/abs/2017ApJ...842....8M/abstract https://doi.org/10.3847/2041-8213/abf14c https://ui.adsabs.harvard.edu/abs/2021ApJ...911L..25M/abstract https://ui.adsabs.harvard.edu/abs/2021ApJ...911L..25M/abstract https://doi.org/10.1088/2041-8205/754/2/L20 https://ui.adsabs.harvard.edu/abs/2012ApJ...754L..20M/abstract http://arxiv.org/abs/2202.03053 https://doi.org/10.3847/1538-4357/aa71b4 https://ui.adsabs.harvard.edu/abs/2017ApJ...842....9M/abstract https://doi.org/10.1093/mnras/stu2619 https://ui.adsabs.harvard.edu/abs/2015MNRAS.447.3936M/abstract https://ui.adsabs.harvard.edu/abs/2015MNRAS.447.3936M/abstract https://ui.adsabs.harvard.edu/abs/2014prpl.conf..521M/abstract https://ui.adsabs.harvard.edu/abs/2007ASPC..376..127M/abstract https://doi.org/10.1088/0004-637X/799/1/41 https://ui.adsabs.harvard.edu/abs/2015ApJ...799...41M/abstract https://doi.org/10.1111/j.1365-2966.2009.15360.x https://ui.adsabs.harvard.edu/abs/2009MNRAS.399.1403M/abstract https://doi.org/10.1111/j.1365-2966.2011.19948.x https://ui.adsabs.harvard.edu/abs/2012MNRAS.419.3074M/abstract https://doi.org/10.1051/0004-6361/201935998 https://ui.adsabs.harvard.edu/abs/2019A&A...630A.142O/abstract https://doi.org/10.3847/0004-637X/832/1/81 https://ui.adsabs.harvard.edu/abs/2016ApJ...832...81P/abstract https://doi.org/10.1093/mnras/stab760 https://ui.adsabs.harvard.edu/abs/2021MNRAS.503.4767P/abstract https://doi.org/10.1093/mnras/stad3462 https://ui.adsabs.harvard.edu/abs/2024MNRAS.527.3876P/abstract https://doi.org/10.1093/mnras/stab2179 https://ui.adsabs.harvard.edu/abs/2021MNRAS.507.1157R/abstract https://doi.org/10.1051/0004-6361/201118085 https://ui.adsabs.harvard.edu/abs/2012A&A...537A.128R/abstract https://doi.org/10.1093/mnras/stv1257 https://ui.adsabs.harvard.edu/abs/2015MNRAS.452..376R/abstract https://doi.org/10.3847/2515-5172/aaff63 https://ui.adsabs.harvard.edu/abs/2019RNAAS...3...16R/abstract https://doi.org/10.3847/1538-4357/ad32d1 https://ui.adsabs.harvard.edu/abs/2024ApJ...966..140S/abstract http://arXiv.org/abs/2505.09578 https://doi.org/10.3847/1538-4357/abda46 https://ui.adsabs.harvard.edu/abs/2021ApJ...910...13S/abstract https://doi.org/10.3847/1538-4357/ace68e https://ui.adsabs.harvard.edu/abs/2023ApJ...954..100S/abstract http://www.ascl.net/0003.002 https://doi.org/10.1086/423135 https://ui.adsabs.harvard.edu/abs/2004ApJS..154..458S/abstract https://doi.org/10.1093/mnras/stw521 https://ui.adsabs.harvard.edu/abs/2016MNRAS.458.3221T/abstract https://doi.org/10.1111/j.1365-2966.2009.14409.x https://ui.adsabs.harvard.edu/abs/2009MNRAS.394.1085T/abstract https://doi.org/10.1146/annurev.astro.45.051806.110525 https://ui.adsabs.harvard.edu/abs/2008ARA&A..46..339W/abstract https://doi.org/10.1086/308093 https://ui.adsabs.harvard.edu/abs/1999ApJ...527..918W/abstract https://doi.org/10.3847/2041-8213/acccf8 https://ui.adsabs.harvard.edu/abs/2023ApJ...948L...2Y/abstract https://doi.org/10.3847/1538-3881/ad08c8 https://ui.adsabs.harvard.edu/abs/2024AJ....167...26Y/abstract 1. Introduction 2. Eccentric Velocity Divergence 3. Observations and Modeling 3.1. ALMA Observations 3.2. Model Setup 3.3. Model Fitting 3.4. Modeling Results 3.4.1. Fitting an EVD Disk Model 3.4.2. Sensitivity to Free Eccentricity? 3.4.3. Which Model is Preferred? 3.5. Conclusion: An Eccentric Gradient in the Fomalhaut Debris Disk 4. Analysis of the Planetary System 4.1. Planet Properties Consistent with ALMA and JWST Observations 4.1.1. A Gap-carving Planet? 4.1.2. A Resonant-clearing Planet? 5. The Origin of Fomalhaut’s Eccentric Disk 5.1. N-body Modeling of the Fomalhaut Debris Disk 5.1.1. REBOUND Setup 5.1.2. REBOUND Results 5.1.3. Comparison with Eccentric Protoplanetary Disks 5.2. Additional Physics? The Role of Disk Self-gravity 6. Conclusions Software and Third Party Data Repository Citations Appendix A. Data Masking Appendix B. Model Verification Appendix C. Residual Emission Maps and emcee Posterior Distributions Appendix D. Model Convergence References