ALMA Reveals an Eccentricity Gradient in the Fomalhaut Debris Disk

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, ef anpow, 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 440 Myr 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.

1. ALMARevealsanEccentricityGradientintheFomalhautDebrisDisk
JoshuaB.Lovell
1
aa ,ElliotM.Lynch
2
,JayChittidi
3
aa ,AntranikA.Sefilian
4
aa ,SeanM.Andrews
1
aa ,GrantM.Kennedy
5
aa ,
MeredithMacGregor
3
aa ,DavidJ.Wilner
1
aa ,andMarkC.Wyatt
6
aa
1
CenterforAstrophysics,Harvard&Smithsonian,60GardenStreet,Cambridge,MA02138-1516,USA;[email protected]
2
UnivLyon,EnsdeLyon,CNRS,CentredeRechercheAstrophysiquedeLyonUMR5574,F-69230Saint-Genis-Laval,France
3
DepartmentofPhysicsandAstronomy,JohnsHopkinsUniversity,Baltimore,MD21218,USA
4
DepartmentofAstronomyandStewardObservatory,UniversityofArizona,Tucson,AZ85721,USA
5
DepartmentofPhysics,UniversityofWarwick,CoventryCV47AL,UK
6
InstituteofAstronomy,UniversityofCambridge,MadingleyRoad,Cambridge,CB30HA,UK
Received2024September30;revised2025June10;accepted2025July8;published2025September4
Abstract
WepresentevidenceofanegativeeccentricitygradientinthedebrisdiskofthenearbyA-typemain-sequence
star,Fomalhaut.Fittingtothehigh-resolution,archivalALMA1.32mmcontinuumdataforFomalhaut(witha
synthesizedangularresolutionof0.
76×0. 55;4–6au),wepresentamodelthatdescribesthebulkpropertiesof
thedisk(semimajoraxis,width,andgeometry) anditsasymmetricmorphology.Thebest-fitmodelincorporatesa
forcedeccentricitygradientthatvarieswithsemimajoraxis,ea
f
n
pow
,ageneralizedformoftheparametric
modelsofE.M.Lynch&J.B.Lovell,withn
pow
=−1.75±0.16.Weshowthatthismodelisstatistically
preferredtomodelswithconstantforcedandfreeeccentricities.Incomparisontodiskmodelswithconstant
forcedeccentricities,negativeeccentricitygradientmodelsbroadendiskwidthsatpericenterversusapocenter,
andincreasedisksurfacedensitiesatapocenterversuspericenter,bothofwhichareseenintheFomalhautdisk,
andwhichwecollectivelyterm“EccentricVelocityDivergence.”Weproposesingle-planetarchitectures
consistentwiththemodelandinvestigatethestabilityofthediskover440Myrtoplanet–diskinteractionsvia
N-bodymodeling.WefindthatFomalhaut’sringeccentricityplausiblyformedduringtheprotoplanetarydisk
stage,withsubsequentplanet–diskinteractionsresponsibleforcarvingthediskmorphology.
UnifiedAstronomyThesaurusconcepts:Debrisdisks(363); Circumstellardisks(235); Eccentricity(441);
Planetary-diskinteractions(2204); Planetarysystemevolution(2292)
1.Introduction
Debrisdisksarecold(tensofK)dustbeltsformedvia
mutualplanetesimalcollisions(M.C.Wyatt2008; B.C.Matt-
hewsetal.2014; A.M.Hughesetal.2018; S.Marino2022).
Whilesuchcollisionsgrinddowndisksovertime,debrisdisks
cansurvivewellinto(andbeyond) thestellarmainsequence,
producingobservabledustsignaturesover10s–1000sofMyr.
Planet–diskinteractionssignificantlyalterthestructuresof
debrisdisks,e.g.,viadepleting,overpopulating,andwarping
planetesimaldistributions.Onseculartimescales,eccentric
planetsforcetheirowneccentricitiesintotheorbitsof
planetesimalsviathree-bodyinteractions(i.e.,inastar–
planet–planetesimalthree-bodysystemC.D.Murray&
S.F.Dermott1999; A.J.Mustill&M.C.Wyatt2012). In
resolvedobservationsofdisks,diagnosticssuchasthestellar
offsetfromthediskcenter,relativebrightnessasymmetries,
andrelativewidthasymmetrieshaveallbeeninterpretedas
metricstoderivediskeccentricitiesandeccentricitydistribu-
tions(M.C.Wyattetal.1999; M.Panetal.2016; E.M.Lynch
&J.B.Lovell2021; J.B.Lovell&E.M.Lynch2023).
Themostfamous,andbest-studied,eccentricdebrisdisk,
isthataroundthe7.66pc-distant,440Myrold,A-type
main-sequencestar,Fomalhaut(derivedfromtheArabic
name“famal-hūt (al-janūbī)” andtranslatestothe“mouth
ofthe(southern) whale”;R.H.Allen1963). Fomalhaut
(HD216956) isthebrightestmemberofsouthernPiscis
Austrinusconstellation,thoughitsdebrishasonlybeenknown
ofsince1985(viaanalysisofitsInfraredAstronomical
Satellite(IRAS) observations;H.H.Aumann1985). Since
then,ithasbeenresolvedathigher-resolutionwiththeHubble
SpaceTelescope(HST;inopticalscatteredlight;P.Kalas
etal.2005, 2008; A.Gaspar&G.Rieke2020), theSpitzer
SpaceTelescopeandJamesWebbSpaceTelescope(JWST;in
near-andmid-infraredthermalemission;K.R.Stapelfeldt
etal.2004; A.Gáspáretal.2023), theHerschelSpace
Telescope(infar-infraredthermalemission;B.Ackeetal.
2012), andtheJamesClerkMaxwellTelescopeandAtacama
LargeMillimeter/submillimeter Array(ALMA;in(sub)
millimeterthermalemission;W.S.Hollandetal.2003;
A.C.Boleyetal.2012;L.Matràetal.2015;M.A.MacGregor
etal.2017; L.Matràetal.2017). Mostrecently,both
G.M.Kennedy(2020) andJ.Chittidietal.(2025) presented
evidenceofawidthasymmetryinthedisk,withabroader
pericenterwidthversusapocenter.Theseobservationsdefini-
tivelyprovedtheeccentricnatureofthedebrisdisk:thestaris
offsetfromthegeometricdiskcenter(inbothscatteredlight
andthermalemission) andthediskpresentswavelength-
dependentdiskansaebrightnessasymmetries,specifically
exhibiting“pericenterglow”ininfraredemission(seee.g.,
M.C.Wyattetal.1999) and“apocenterglow”inmillimeter
emission(seee.g.,M.Panetal.2016). Bothoftheselatter
featuresareexpectedbasedonthestructureandthermal
conditionsinFomalhaut’sdiskattheobservedresolution
scales(E.M.Lynch&J.B.Lovell2021).
Here,wepresentamillimetermodelofFomalhaut’smain
(outer) debrisdiskbeltthatincorporatesaforcedeccentricity
gradientparameter,whichwefittothehigh-resolutionALMA
TheAstrophysicalJournal,990:145(18pp), 2025September10 https://doi.org/10.3847/1538-4357/adfadc
©2025.TheAuthor(s). PublishedbytheAmericanAstronomicalSociety.
aaaaaaa
Originalcontentfromthisworkmaybeusedundertheterms
oftheCreativeCommonsAttribution4.0licence. Anyfurther
distributionofthisworkmustmaintainattributiontotheauthor(s) andthetitle
ofthework,journalcitationandDOI.
1

2. imagesofJ.Chittidietal.(2025). Thismodelprovidesanovel
interpretationofthedisk’smorphology;incorporatingasteep,
negativeeccentricitygradientwithsemimajoraxis.In
Section2,weoutlineaphysicaldescriptionofthe“eccentric
velocitydivergence”(EVD) effectthatdescribesthistypeof
diskstructure.InSection3,wedescribeourmodelsetup,
fittingmethodology,andresults.InSection4,wederiveplanet
propertiesconsistentwithobservations.InSection5,we
considerthestabilityoftheseplanetaryscenarioswithN-body
modeling.WesummarizeourconclusionsinSection6.
2.EccentricVelocityDivergence
Centraltotheinvestigationthatweoutlineinthisworkis
thedependencybetweendisksurfacedensitydistributionsand
theirunderlyingeccentricitydistributions.Thishasbeen
demonstratedrecentlyintheliteraturebyE.M.Lynch&
J.B.Lovell(2021) wheresimplepower-lawforcedeccen-
tricityprofiles()eaa
f
n
pow
withfixedn
pow
valuesof0,−1,or
+1weremodeled(wherediskeccentricitygradientswere
eitherflat(∂e
f
/∂a=0),positive(∂e
f
/∂a>0),ornegative
(∂e
f
/∂a<0),respectively). Thesemodelswereshowntoalter
theresultingdisksurfacedensitiesandemissionmorphologies.
Inthiswork,weadoptthegeneralizedform()eaa
f
n
pow
(allowingfornonintegervalues) toexaminepower-law
eccentricityprofilesinthecontextoftheFomalhautdebris
disk.Weillustratetheinfluenceofthiseccentricitydistribution
inFigure1showingnormalizeddustsurfacedensitiesfora
subsetofpower-lawprofiles(0.5,0,−0.5,−1,and−2,
definedby() ()/=eae
aa
f f
n
,0 0
pow
)adoptingparameterscon-
sistentwiththeFomalhautdisk(i.e.,asemimajoraxisof
a
0
=139au,aGaussiandiskwidthofσ
r
=15.6au,anda
forcedeccentricitye
f,0
=0.127ata
0
)andanarbitrarilychosen
ω
f
=0.0.Theseplotsshowthesamebehaviorasdescribedby
E.M.Lynch&J.B.Lovell(2021), i.e.,thatsurfacedensities
indiskswithpositiveeccentricitygradientsareenhancedat
pericenter,whereasdiskswithnegativeeccentricitygradients
havesurfacedensitiesraisedatapocenter.Inaddition,these
modelsshowthatdiskwidthsdecreasewherethesearemost
dense,andbroadenwheretheseareleastdense.Forexample,
negativeeccentricitygradientprofilespreferentiallybroaden
diskwidthsatpericenter,andpreferentiallynarrowdiskwidths
atapocenter.Tosimplifydiscussionofthecombinedimpact
ondisksurfacedensitiesandwidths,wetermthiseffect
EccentricVelocityDivergence(EVD), whichsimplyarises
duetomasscontinuityineccentricdisks.
Onecanconsidertheeffectonthediskwidthanalyticallyby
calculatingtheradialdistancebetweentwoconcentricorbitsat
pericenter(Δr
peri
)andapocenter(Δr
apo
).Inthecaseoffixed
positiveeccentricities,thedefinitionsofthepericenterand
apocenterradiiresultinwidthsΔr
peri
=Δa(1 −e)and
Δr
apo
=Δa(1 +e),i.e.,Δr
apo
>Δr
peri
alwaysfore>0.
However,forthegeneralcaseofe=e(a)(retainingonlylinear-
ordertermsintheexpansionofe=e(a)),thesamederivation
findsΔr
apo
±Δr
peri
≈Δa(1±(e+a∂e/∂a)). Consequently
onecanshowthatapocenterandpericenterwidthsare
dependentoneccentricitygradients,withΔr
apo
>Δr
peri
requiring∂e/∂a ≳−e/a.Forthecaseof()eaa
f
n
pow ,one
findsthatΔr
apo
>Δr
peri
isonlytrueifn
pow
≳−1,whichcanbe
seenbycomparingthee
f
(a)∝a
−1
model(columnfourof
Figure1)withotherpresentedmodels.Overall,EVDcaninduce
similarfeaturestothoseobservedintheFomalhautdebrisdisk,
soweconsiderthepossibilityofitsoperationinthisdisk.
Power-laweccentricityprofilesarephysicallymotivated.
Forexample,intheidealizedsituationofaneccentricthree-
bodysystemwitheitheraninner-planetaryperturber(i.e.,a
star-eccentricplanet-planetesimalsystem), oranexternal-
planetaryperturber(i.e.,astar–planetesimal–eccentricplanet
system) thedisturbingfunctionpredictsanalyticalrelation-
shipsbetweenthesemimajoraxesandeccentricitiesofthe
perturbingplanetandtheplanetesimal.Tofirst-order(i.e.,
undertheconditionsofloweccentricityandwell-separated
star–planet–(massless) planetesimals), theserelationshipscan
bederivedaseithere
f
(a)∝a(forexternal-planetary
perturbations) ore
f
(a)∝a
−1
(forinternal-planetaryperturba-
tions). Thesegradientslopescansteepen,forexample,whenea
0.5
r
ea
0
r
ea
0.5
r
ea
1
r
ea
2
r
Figure1.Normalizedsurfacedensitymaps(face-onprojections,in(top)r–fspace,from0to2πand0to180au,and(bottom) x–yspace) fordifferentpower-law
eccentricityprofiles(withtheire=e
f
(a)functionalformshowninthelowerrightofeachpanel). Themodelsallhaveapse-alignedargumentsofpericenter,aswell
asradii,widths,andeccentricitiesconsistentwiththoseofFomalhaut,andanargumentofpericenterdirectedsouthortoward0or2πinr–fspace.Contoursin
increasingbrightnessscalerepresentregionsofhighersurfacedensity,setatthelevels10%,30%,50%,70%,and90%(10%levelsareshownwithwhite-dotted
contourlines).
2
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.

3. thesemimajoraxisofaperturbingplanetapproachesthoseof
itsperturbersorifplanetesimalsareinsteadmodeledwith
mass(seee.g.,A.A.Sefilianetal.2021; A.A.Sefilian2024).
Indeed,sucheccentricitydistributionsarepresentinthe
KuiperBelt,wherebythedistributionoftrans-Neptunian
Objectsinthecold(classical) belthoststeepe
f
(a)profiles
(e.g.,atthelevelofafewpercentperau;R.I.Dawson&
R.Murray-Clay2012). Moreover,theoreticalworkon
eccentricprotoplanetary(gas)diskshasshownthateccentricity
gradientsarenecessarytoavoidviolentdifferentialprecessionin()=eacons
t.
f
disks(J.Teyssandier&G.I.Ogilvie2016). If
debrisdisksformfromprimordialmaterialwithsucheccentricity
distributions,orareshapedbyeccentricplanet–diskinteractions,
itisplausiblethattheirstructuresmaybewelldescribedbythe
EVDeffect.
3.ObservationsandModeling
3.1.ALMAObservations
Forthisstudy,wefittotheimagesofFomalhautfrom
J.Chittidietal.(2025). Thatworkpresentsthecalibrationand
dataalignmentprocedureforthethreeALMAband6datasets
included(P.I.:AaronBoleyADS/JAO.ALMA#2013.1.00486.
S,P.I.:PaulKalasADS/JAO.ALMA#2015.1.00966.S, andP.
I.:MeredithMacGregorADS/JAO.ALMA#2017.1.01043.S)
andutilizesstandardCASAALMAreductionpipelines.Weuse
thesametcleanparametersforimagingasJ.Chittidietal.
(2025) withaBriggs-weightingschemeandrobustparameterof
0.5,whichproducesamosaicimagecenteredonthe2015
coordinatesofFomalhaut,withapixelsizeof0.
075anda
synthesizedcleanbeamwith0.
76×0. 55,andbeamposition
angle(PA)−87°.
4.Atthe7.66pcdistancetoFomalhaut,this
imageprovidesaphysicalresolutionof6×4au.Wepresent
theimageinFigure2(left). Thethreeobservationsthatwere
combinedtoproducethisimagewereallconductedwith
differentmosaics,sensitivities,andALMAconfigurations(and
thuseffectiveangularresolutions). Asaresult,thefinalimage
hasnon-Gaussiannoisepropertiesandprimary-beam
attenuationthatisdominatedbythetwo-pointingmosaicof
ADS/JAO.ALMA#2017.1.01043.S, whichhadphasecenters
directedtowardthediskansae.AscanbeseeninFigure2,the
dataisthereforemuchnoisierneartotheminoraxisofthedisk
incomparisontothetwoansae.
3.2.ModelSetup
Weproduceparametricmodelsofopticallythindisks
consistentwiththosepresentedinJ.B.Lovelletal.(2021),
E.M.Lynch&J.B.Lovell(2021), J.B.Lovelletal.(2022), and
J.B.Lovell&E.M.Lynch(2023), i.e.,byimagingthree-
dimensionaldustdistributionswithintheRADMC-3Dframework
(C.P.Dullemondetal.2012). Wefixthemodelwavelengthto
thatoftheband6observations(λ=1.32mm).Wefixthestellar
parameterstovaluesconsistentwithE.E.Mamajek(2012), i.e.,
theeffectivetemperature(T
�
=8500K),stellarradius
(R
�
=1.8R
⊙
),andstellarmass(M
�
=1.9M
⊙
),whichtogether
fullydefinetheKuruczstellartemplatespectra(see
R.L.Kurucz1979) fromwhichthetemperatureconditions
throughthediskarederivedbyRADMC-3D. Thefluxdensityat
thelocationofthestar(asadoptedlater) isafixedvalue
independentofthesestellarparametersandthedisktemperature
conditions.Themodelsassumethedusttohaveafixedsize
distributionbetween0.9μm–1.0cmwitha−3.5power-law
index(asperJ.S.Dohnanyi1969), andwithgraindensitiesof
ρ=3.5gcm
−3
.WeadoptasingleGaussiansemimajoraxis
distribution,basedonthe(mean) disksemimajoraxis(a
0
)and
Gaussian-width(w
r
,whichisrelatedtothediskfullwidthat
half-maximum(FWHM), by=FWHM22ln 2
w
r
).Thesur-
facedensityofthediskismodeledwitha(mean) forced
eccentricity(e
f,0
,measuredwithrespecttoa
0
),anargumentof
forcedeccentricity(ω
f
,measuredanticlockwisefromsouth), and
apower-lawindexofn
pow
(forapower-laweccentricity
distributionoftheform() ()/=eae
aa
f f
n
,0 0
pow
),asper( )
()
=
wj
aa
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
Relative Decl. Oset [arcsec]
20
0
20
40
60
80
100
120
-10.00.010.0
Relative RA Oset [arcsec]
ef/a
1:75
20
0
20
40
60
80
100
120
-10.00.010.0
Intensity [

Jy beam

1
]
40
30
20
10
0
10
20
30
40
Figure2.Left:ALMAdataaspresentedinJ.Chittidietal.(2025) whichwefittointhiswork.Center:Best-fite
f
(a)∝a
−1.75
model(onsameimagescale). Right:
Residualemissionafterwesubtractourbest-fitmodelfromtheALMAdata(left). Contoursareshownatthe±3σlevelintheresidualsandatthe5σlevelforthe
data/model. Beamsareshowninthelower-leftofeachplot.Emissionremainspresentintheminoraxisregionoftheresiduals;however,thisisdominatedby
imagingnoise(enhancedduetotheprimary-beamresponse) andbackgroundsources(suchasthe“GreatDustCloud,”seeA.Gáspáretal.2023andG.M.Kennedy
etal.2023). Northisup,eastisleft.
3
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.

4. withjthecoordinatesystemJacobiandeterminant(asalso
presentedasEquation(8)inE.M.Lynch&J.B.Lovell2021)( )
( ) ()
/
=
+
j
eeaea e
qE
1 1
1 cos 2
2
fortheeccentricanomaly(E)( )
( )
()
=
+
+
E
e
e
cos
cos
1 cos
3
f
f
andtheorbitalintersectionparameter(q,inthecaseofan
untwisteddisk,i.e.,∂ω
f
/∂a=0)( )
()
/
/
=
+
q
aea
eeaea1
. 4
Toavoidorbitalintersections,werequire|q|<1.Equation(1)
fullydescribestheplotspresentedinFigure1(griddedonr–f
andx–yspace).
7
WeadoptasingleGaussian-distributedverticaldensity
profile,withaverticalaspectratio(h=H/r) forHthe
classicalscale-heightforagivendiskradiusr,suchthatthe
verticallineardensityis()=l
H
z
H
1
2
exp
2
. 5
2
2
Thisparameterizationdefinesthetotaldensityρ=M
dust
Σl
ρ
,
whereweintroduceadustmassscalingparameter(M
dust
)
whichuniformlyscalesmodelpixelintensitiestomassunits.
Bydefault,ourmodelsareface-onandorientednorth–south
(withthepericenterdirectionsouthwards,apocenterdirection
northwards,asperFigure1).WeuseRADMC-3Dtoproject
andimagemodelsonthesky-plane.Weprovideparameters
forthePA(whichrotatesourmodelanticlockwisereferenced
tonorth), inclination(i,whichinclinesourmodelwith
referencetoanglesbetweeni=0°forface-onprojections,
andi=90°foredge-onprojections), andoffsetsintheright
ascension(ΔR.A.) anddeclination(Δdecl.), respectively
(applieduniformlytoboththestellarlocationanddiskfocus).
WepresentinFigure2(middle) onesuchimagedmodelwith
Fomalhaut’sdiskandstellarparameters(convolvedwiththe
samecleanbeamparametersastheimagedALMAdata,using
thebest-fitparametersinferredlater).
3.3.ModelFitting
Ourfittingapproachisconsistentwiththeproceduresof
J.B.Lovelletal.(2021, 2022), exceptinthisworkwefitinthe
imagedomaininsteadofthevisibilities.
8
Wefitourmodelto
theALMAdatawithemcee, an“AffineInvariantMarkov
ChainMonteCarloEnsembleSampler”
9
(J.Goodman&
J.Weare2010; D.Foreman-Mackeyetal.2013). Samplesof
theposteriordistribution(generatedbyemcee) aredeter-
minedfromtheassumedpriordistributionsandaGaussian
likelihoodfunction(( )/exp 2
2
),with()=
=
M
I I
I
l
, 6
ij
2
,1
N
data,i,j m
odel,i,j
2
data
2
pix
2
beam
i,j
pix
foranimagewithsizeN
pix
×N
pix
,intensityineachpixel(I,
foreitherthemodelordata),theimagenoise(δI
data
,seefurther
below), thebeamsize(Ω
beam
),imagepixelsize(l
pix
,wherethe
latterterm/l
pix
2
beam
reweightsχ
2
toaccountforpixel–pixel
correlation), andthefittingmask(M
i,j
).Wefittotheprimary-
beamcorrectedimage(I
data,i,j
)andmeasuretheχ
2
inthe
residualmapcorrectedbytheprimary-beamresponse,
assumingtheCASAformofALMA’sprimarybeam(as
describedine.g.,ALMAKnowledgebase,Mar2025
10
).The
ALMAimagehasanon-Gaussiannoisedistributionmost
prominentneartothediskminoraxisgiventheimagederives
fromdifferentmosaicpatternsfromdifferentALMAconfig-
urations(seeFigure2).Bycomparingfitstodatawithdifferent
maskingconditions,weproceededtofitwithamaskdefined
byacut-offintheALMAimagedprimarybeam(onlyfitting
emissionwheretheprimary-beamresponseexceeds0.66), and
aprojectedradiuscut-off(onlyfittingtoemissionwithina
projected200audistancefromthestellarcenter). Wepresent
inAppendixA(Figure5)thelocationofthefittingmaskand
theprimary-beamresponseforthedata,aswellasthe
discussionofandresultsoffittingtheGeneralmodeltothe
datawithoutincludingtheprimary-beamcut-offinthemap.
Themasksimplymultipliespixelvalueswitheither=M 1
i,j
or=M 0
i,j andisappliedtoboththedataandmodel
identically.Themaskchoicefocusesthefittingonhighsignal-
to-noiseratio(SNR) diskemission,avoidingimagingartifacts
andbackgroundsourcesnearthedisksemiminoraxes(e.g.,
suchasthe“GreatDustCloud”(GDC) asdiscussedby
G.M.Kennedyetal.2023). Wenotethatbymaskinginthe
fittingstageratherthanduringmodelgeneration,ourmodels
stillcompriseazimuthallycompletedisks.
WeadoptanuncertaintyofδI
data
=8.0μJybeam
−1
based
onthermsofpixelintensitiesindisk-freeregionsofthe
ALMAimage,whichisscaledbytheprimary-beamresponse
duringfitting,down-weightingnoisierregionsofthedata.We
findfromtestiterationsoffittingmodelstothedatawith
emceethattheresidualmapshavestandarddeviation
≈8.0μJybeam
−1
(i.e.,asexpectedgiventheimagedata)
andmean<1.0μJybeam
−1
(i.e.,verycloseto0,asexpected
forwell-fittedmodels). Toverifythefittingmethodology,we
testedthisprocedureontheCycle3-onlyALMAdataof
FomalhautasdescribedinAppendixB.
3.4.ModelingResults
3.4.1.FittinganEVDDiskModel
Wefitthestellarfluxandlocationseparatelytothediskto
reducethenumberoffreeparametersinourfitting(sincewe
simplyfitthestellarflux,weleavefixedthephysicalstellar
inputsthatderivethedisktemperatureconditions,e.g.,R
�
,M
�
,
andT
�
).Weadoptasmallaperturemaskoverthefullimageto
fitthestellarfluxandlocation,maskingallbuta10×10pixel
7
Thecodeforgeneratingeccentricringmodelsandreproducingtheseplots
isavailableviaDOI:10.5281/zenodo.16758671 andathttps://github.com/
astroJLovell/eccentricDiskModels.
8
Itisinprinciplepossibletofitsuchmodelsinthevisibilitydomain;
however,thisdemandsamuchgreatercomputationaltimeinvestment.We
notethatmodelconvergencecanbeachievedwiththisimage-basedsetupina
modest1–2days,whereasinthevisibilitydomainamuchmoreintensive1–2
weekspermodelisrequired(seeJ.Chittidietal.2025), andthustheapproach
hereenablesustoinvestigateabroadarrayofmodelsetups.
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
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.

5. box.Wefindthatthestarisbest-fittedwithoffsetsconsistent
with0.0inbothaxes(ΔR.A.=0.
01±0. 04and
Δdecl.=0.
01±0. 04,inaccordancewiththerealignment
andimagingprocedureofJ.Chittidietal.(2025), andafluxof
0.737mJy(i.e.,F
�
=0.74±0.08mJyaccountingfora10%
ALMAfluxcalibrationuncertaintyinquadraturewithour
fittingerror).
Werunaninitialfittingroutinewithemcee(with48
walkers,3000steps) toestimatetheconvergencetimescaleand
tofurtherverifythefittingresultsforthehigher-resolution
data.Thisinitialfittingroutineshowedthatseveralthousand
morestepswererequiredtoachieveconvergence,whichwe
discussindetailinAppendixD.Nevertheless,thephysical
modelparameterswemeasuredwereallconsistentwith
previousworks.Ofnote,wefittedPA=336°.
19±0°. 10,a
(Gaussian) verticalaspectratioofh=0.0157±0.0013and
ω
f
=15°. 6±0°. 6.TheGaussianverticalaspectratioisslightly
lowerthanthe0.02assumedinM.A.MacGregoretal.(2017),
andconsistentwiththatfittedbyG.M.Kennedy(2020).
11
GiventheresolutionoftheALMAdata,thelowvaluethatwe
determineforh,andthelowhconstrainedbypreviouswork,
wenotetheverticalstructureisplausiblystillunresolved.
Furthermore,inG.M.Kennedy(2020), thediskrotationwas
modeledwiththelongitudeofascendingnodeΩ=156°.
4
(measuredanticlockwisefromnorth) whichisdirectedtoward
thesouthwestansa,whereaswemodelPA=336°.
19antic-
lockwisefromnorthtoalignwiththenortheastansa(i.e.,these
aresimply180°offset,andthusfullyconsistent). Inaddition,
thedefinitionoff
ofG.M.Kennedy(2020) ismeasuredwith
respecttothelongitudeofascendingnode,whereasthe
definitionweadoptiswithrespecttosouth,i.e.,independent
ofthePA(withbothagainmeasuredanticlockwise). Assuch,
thevalueof= ±411
f
reportedinG.M.Kennedy(2020)
wouldhavebeenmeasuredas= +f
f
,i.e.,
17°.
1±1°. 0viathemodelsetupwedescribehere.Inthis
work,weconsistentlyfitvalueswithinafewdegreesofthis
valuewithsimilaruncertainties,andthuswedeemtheseto
likewiseagree.
Toreducethenumberofparameterswefitto,weproceeded
byfixingΔRAandΔDec,andfittedonlyforthediskmass
(M
disk
),semimajoraxis(a
0
),diskwidth(w
r
),(mean) forced
eccentricity(e
f,0
),argumentofforcedeccentricity(ω
f
),the
power-lawindexfortheforcedeccentricitygradient(n
pow
),the
inclination(i),thePA,andthe(Gaussian) verticalaspectratio
(h).Wetermthisthe“General”model.Wepresentin
AppendixCtheposteriordistributionoftheconvergedemcee
chains,anddiscussitsconvergencecriteriainAppendixD
afterrerunningthefittingroutinefor6000steps.Wefinda
significantdeterminationoftheforcedeccentricitygradient
power-lawindex(n
pow
=−1.75±0.16)andparametersfor
thediskmass,semimajoraxis,width,eccentricity,argumentof
pericenter,inclination,PA,andverticalaspectratiothatareall
consistentwithM.A.MacGregoretal.(2017) and(where
presented) G.M.Kennedy(2020).
Weassessthepropertiesofourresidualemissionmapto
investigateifourmodelisagoodfit.Forthecompletemasked
region,wemeasureapixelvarianceequaltothesquareofour
measuredσ
image
,apixelmeanequalto0.1×σ
image
,anda
distributionthatfollowsaGaussian-distribution.Thesevalues
implythatonlyminimalexcessemissionremainsinthefinal
residualimage(whichmaybefullyexplainedbythepresence
ofbrightpointsources,whichmaybebackgroundsources)
andthatitvariesaswouldbeexpectedifitispopulatedonly
bynoise.TheresidualmapasshowninFigure2appearsto
demonstratetheGeneralmodelisaverygoodfittothedata.
Forexample,intheansaethebroaderandfainterpericenterare
fittedwithinthenoise,andintheapocenterthepeakis
reproducedwithinthenoise.Intheresidualmapsatthetwo
ansae,weseeonlysmall(sub-beamsized) residualsthat
exceed3σ(noneofwhich≳4σ). Alongtheminoraxes,the
modelappearstobeaworsefit,giventhemodelisfainterthan
thedatainthoselocations.Someoftheselargerpeaksare
eitherbackgroundsources(i.e.,GDC,asnotedbyG.M.Ken-
nedyetal.2023) ornoisierregionsofthedata.Theresidual
emissioncoincidentwiththeapocenterouteredgeappearsto
beslightlyoversubtractedbyourmodel,whichsuggeststhat
ourwidthmaynotbequitenarrowenoughinthislocation.
12
Theseresidualfeaturesmaythereforesuggestthatasteeper
eccentricitygradientmodelmaybeneeded,buttheexisting
datadoesnotprefersuchamodel.
Similarnegativeapocenterresidualsarepresentinthe
residualmapsofG.M.Kennedy(2020) andalsothefitwe
performedtothelower-resolutiondata(discussedin
AppendixB).However,wenotethattheresidualsof
G.M.Kennedy(2020) hostinadditionastrongpositive
residual(locatedin-betweenthenegativeresidualregions),
i.e.,thosemodelscannotaccountfortheapocenterbrightness
enhancementinthedisk.Incontrast,wefindnopositive
residualsinthediskansaeinourfitstoboththehighandlow-
resolutiondata.Ourfindingsthuspresentevidencethatthe
Fomalhautdebrisdiskcanbeunderstoodintheparadigmof
EVD,withthishostingthefirstevidenceofaneccentricity
gradientinitsdebrisbelt.Wecollateourbest-fitparameter
resultsinTable1forallourfittingsetups.Weinvestigatein
Section3.4.2whetherthedataprefersamodelwithfree
eccentricityandeitheraconstantornegativeforcedeccen-
tricityprofile.
3.4.2.SensitivitytoFreeEccentricity?
Totestthepossibilitythatthepower-laweccentricity
distributionmodelfitisbiasedbynotaccountingfor
the“free”(or“proper”) eccentricityparameter(e
p
),we
produceadiskmodelwhichincludes,e
p
,e
f
,andapower-
lawdistributionintheforcedeccentricity.Wesolve
analyticallytheorbitequationforthecaseofadiskwith
bothfreeandforcedeccentricitiesfollowingea
n
pow
with
fixedvaluesofn
pow
=0,−1(andleavetofuturework
11
AlthoughwenotethattheupperlimitsoninclinationinG.M.Kennedy
(2020) werenotasstringentastheyshouldhavebeengiventhedata,wehave
tracedthistoacodingerrorthatresultedintheinclinationparametersonly
affectingradialstructure(i.e.,thediskmodelswereflat).Withthevertical
structureparameterizationcorrectlyincluded,wereranthosemodelsand
insteadyieldanupperlimitontheGaussianverticalaspectratioofh<0.03,
withonlyverysmalldifferencesforotherparameters(e.g.,e
f
=0.126±0.001
comparedto0.125±0.001previously).
12
Toinvestigatethis,weattemptedamorecomplexradialprofilefittothe
data,allowingdistinctw
r
valuesinternal/external tothesemimajoraxisin
casetheseresidualsresultedfromfittingaGaussiansurfacedensityprofileto
datathatmaynotbeGaussianinemission.Whilethefityieldedimproved
residuals(withinner-edgeandouter-edgewidthsofw
r,in
=11.2±0.7auand
w
r,in
=15.7±0.7au,respectively,asmallermeansemimajoraxisof
a=137.6±0.4au,andotherwisestatisticallyconsistentbest-fitmodel
parameters), thismodelneitherfullyremovedthesenegativeresidualfeatures,
norwastheimprovementinχ
2
statisticallysignificantversusthesimpler
Gaussianmodel.
5
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.

6. thegeneralcaseforanyvalueofn
pow
,i.e.,for=e () [( )] [( ]
)/
+eaa i e iexp
exp
f
n
f p p,0 0
pow
,forasys-
temwith(mean) forcedandfreeeccentricities(e
f,0
ande
p
,
respectively) andargumentsofforcedandfreeeccentricities
(ω
f
andω
p
,respectively). Tomodelthefreeeccentricity,we
adoptthesameapproachasJ.B.Lovell&E.M.Lynch
(2023), i.e.,welinearlysumindependent“families”oforbits
(eachcomprisingasingleargumentoffreeeccentricityω
p
)
uniformlyovertherangeof0�ω
p
<2π,whereeachfamily
comprisesadiskwithM
disk
,a
0
,w
r
,e
f,0
,ω
f
,e
p
,and(fixed) n
pow
.
J.B.Lovell&E.M.Lynch(2023) demonstratedthatthis
parametricmodelingapproachforthefreeeccentricity
distributionyieldsconsistentresultstomethodsthatinstead
modelindividualparticles(e.g.,M.A.MacGregoretal.2017
andG.M.Kennedy2020, wheren
pow
=0wasimplicitly
adopted). Giventhenegativebest-fitn
pow
thatweearlier
inferred(andthatpreviousmodelingeffortsimplicitlyhad
n
pow
=0),weconsiderhereonlythetwocasesofn
pow
=0and
n
pow
=−1.
First,byfittingasix-parametermodel(i.e.,withn
pow
fixed
to0)ande
p
allowedtofreelyvary,werunthesamemodel
fittingproceduredescribedearlierwithemcee(for3000
steps) andfindthatthebest-fitsolutionresultsin
e
p
=3.85%±0.17%,i.e.,thisrequiresasignificantfree
eccentricity.Wetermthismodel“Proper1.”Second,
rerunningthisprocessinsteadwithn
pow
fixedto−1(for
5000stepstoensureconvergence), wefindthatthebest-fit
solutionresultsine
p
=1.9%±0.9%(forwhichweplacean
upperlimitone
p
<3.59%basedonthe99.7th-percentileofall
chainsafterburn-in). Wetermthismodel“Proper2.”We
presentthemodelandresidualsofthesefittingprocessesin
Figure3,allbest-fitvaluesinTable1,andthecornerplotsfor
theposteriorsinFigures9and10.Onecanseethatwhilethe
n
pow
=−1providesareasonablefit(albeitpoorerthan
themodelwithasteepereccentricitygradientandnofree
eccentricity), themodelwithn
pow
=0andasignificantfree
eccentricityispooreratinterpretingthebroaderwidthat
pericenterandthebrightnessenhancementatapocenter
consistentwiththefindingsinJ.Chittidietal.(2025),
andalsothoseofM.A.MacGregoretal.(2017) and
G.M.Kennedy(2020). Wenotethatbydefinitionthereisa
degeneracybetweene
p
andw
r
,whichresultsinthewiderange
ofw
r
posteriordistributions.
3.4.3.WhichModelisPreferred?
InTable1,wedefineΔχ
2
astherelativedifferenceinthe
best-fitvalueχ
2
ofthe“General”modelversusthemodels
withfreeeccentricities,wheretheGeneralmodelreturnsa
lowerχ
2
thantheProper1and2modelsby70and22,
respectively.Wequantifythesedifferencesstatisticallyby
estimatingtheassociatedprobabilitythat(giveneachofthe
modelshasninefreeparameters) therelativedifferenceinχ
2
issignificantlyimprovedbycalculatingtheIncomplete
GammafunctionforeachvalueofΔχ
2
,i.e.,( )
()
/
=P
N
t tdt
1
2
exp . 7
N
par
0
2
1
2
par
Intermsofconfidenceintervals,theserespectiveΔχ
2
values
thereforeimplythattheGeneralmodelisfavoredoverthe
Proper1and2modelsby6.7σand2.6σ,respectively,andthat
theProper2modelisfavoredoverProper1by5.2σ.Wenote
thatinathirdversionofthePropermodelwhichinsteadhada
fixedeccentricitypower-lawgradientof−2,wefoundthe
Generalmodelstillhadalowerχ
2
by10.6,andthusperformed
statisticallyaswell,withapreferenceofonly1.0σ.Overall,
thesedifferencesimplythatwhilethediskmayhostsome
smallamountoffreeeccentricity,theeccentricdiskis
dominatedbythepresenceofaforcedeccentricitydistribution
whichfallswithsemimajoraxis,giventhestrongpreferenceof
theGeneralandProper2modelsovertheProper1model,and
furtherthepreferenceoftheGeneralmodelovertheProper2
model.
Althoughfittingmodelstothedatainthevisibilitydomain
requiresanorderofmagnitudelongertoconverge,wenote
thattheoverallresultwepresentisrobusttothemethodology
applied,andavisibility-domainanalysiswouldconclude
similarly.Weverifiedthisbyproducingfromeachofourthree
best-fitGeneral,Proper1andProper2models,visibilitytables
(griddedtothesameUV-spatialfrequenciesasthedata)from
theCycle5,Apocentre-onlypointingwhichcontainsthevery
highestSNRdata(seeJ.Chittidietal.2025). Bysubtracting
thesevisibilitiesfromthedatainthevisibilitydomain(after
applyingtheCASAstatwtfunctiontothedata), we
measuredtheirrelativeχ
2
values(directlyfromthevisibilities)
andimagedtheresiduals.Viathismethod,theresiduals
stronglyfavortheGeneralmodel,whichleavescomparable
Table1
Best-fitModelResultsfortheFitstotheHigh-resolutionALMAData(top)andtheBest-fit“Verify”ModelFittotheLower-resolutionData(Bottom)
Parameter Units General Proper1 Proper2 Verify
M
dust
10
−2
M
⊕
18.70±0.20 20.90±1.3 21.70±1.3 19.90±0.20
a
0
au 138.79±0.12 138.93±0.12 138.84±0.12 138.89±0.10
w
r
au 13.51±0.29 9.3±0.6 12.8±1.0 15.3±0.3
e
f,0
% 12.56±0.12 12.56±0.11 12.57±0.12 13.46±0.14
e
p
% [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
n
pow
⋯ −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
valuesareallwithreferencetotheGeneralmodel.Valuespresentedinsquarebracketswerefixedtothepresentedvalueduringfitting.For
completenesswealsopresenttheparametersdeterminedfortheverificationmodelwhichonlyusedthelow-resolutiondataofG.M.Kennedy(2020).
6
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.

7. residualstothoseseeninFigure2(i.e.,noresidualartifacts
above4σ)andfurtherhasthelowestχ
2
by350(versusProper
1)and245(versusProper2).
3.5.Conclusion:AnEccentricGradientintheFomalhaut
DebrisDisk
Wefindclearevidencethatthereisagradientinthe
eccentricitydistributionofFomalhaut’splanetesimals,and
thusinitsdebrisdustasobservedwithALMA.Thisfinding
suggeststhattheFomalhautdebrisdiskisbeingshapedbythe
EVDeffect,asdiscussedinSection2.Mostintriguingabout
theinferredpower-lawgradientvalueoftheGeneralmodel
(n
pow
=−1.75±0.16) isitsnear-consistencywiththe
classicalexpectationforsecularinteractionsbetweenan
internal-planetaryperturberwithdistantlyseparatedouter
planetesimals(forwhichn
pow
=−1;seee.g.,C.D.Murray
&S.F.Dermott1999, i.e.,thesamevaluethatwefixedinthe
Proper2model,andfittedtothearchivalALMAdata). The
valuefittedtothenewestdataneverthelessappearssteeper
thanthisclassicalexpectation(beingdiscordantatthe4.6σ
level,thoughwenotethisdiscrepancydisappearswhenwefit
tothesingleconfiguration,lower-resolutionmosaicof
M.A.MacGregoretal.2017, whichisfullyconsistentwith
apower-laweccentricitygradientof−1).Wedeferfurther
discussiontotheoriginofthiseccentricityprofiletoSection5.
4.AnalysisofthePlanetarySystem
SincetheinferenceofaneccentricitygradientinFomal-
haut’sdiskhasimportantimplicationsforthepresenceand
orbitalpropertiesofaninternalplanetinteractingwiththedisk,
wenextinvestigateplanetpropertiesplausiblewiththis
interpretation.Wefirstconsiderplausiblesingle-planet
scenariosthatcouldberesponsibleforcarvingthedisk
structureinSection4,i.e.,basedonthemodeledpropertiesof
thedisk’sinneredgeandinner-edgeeccentricity.InSection5,
weconsiderifeithersuchplanetaryscenarioispreferredand
theplausibleoriginsofthedisk’seccentricity.Specifically,in
Section5.1weconductN-bodymodelingofplanet–disk
interactionsconsistentwiththederivedplanetscenarios,andin
Section5.2weanalyticallyinvestigatewhetherdiskself-
gravity(inthepresenceofaninternalplanetperturber
consistentwiththederivedplanetscenarios) maybe
responsiblefordrivingasteepergradient(thann
pow
=−1)
intothedisk.-10.00.010.0
-20.0
-10.0
0.0
10.0
20.0
Relative Decl. Oset [arcsec]
20
0
20
40
60
80
100
120
-10.00.010.0
Relative RA Oset [arcsec]
ep,ef/a
0
20
0
20
40
60
80
100
120
-10.00.010.0
Intensity [

Jy beam

1
]
40
30
20
10
0
10
20
30
40 -10.00.010.0
-20.0
-10.0
0.0
10.0
20.0
Relative Decl. Oset [arcsec]
20
0
20
40
60
80
100
120
-10.00.010.0
Relative RA Oset [arcsec]
ep,ef/a
1
20
0
20
40
60
80
100
120
-10.00.010.0
Intensity [

Jy beam

1
]
40
30
20
10
0
10
20
30
40
Figure3.Left:ALMAdata.Center:Best-fitmodelswhichincludee
p
(onsameimagescale,withtopforthemodelfixedwithn
pow
=0,andn
pow
=−1forbottom).
Right:Residualemissionafterwesubtractthebest-fitmodelsfromtheALMAdata.Contoursareshownatthe±3σlevelintheresidualsandatthe5σlevelforthe
data/model. Beamsareshowninthelower-leftcornerofeachplot.Inallplots,northisupandeastisleft.
7
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.

8. 4.1.PlanetPropertiesConsistentwithALMAandJWST
Observations
Ontheassumptionthatasingle,coplanarplanetis
responsibleforcarvingtheobserveddiskstructureinthe
Fomalhautdebrisdisk(i.e.,itsinner-edgelocation,eccen-
tricity,andeccentricitygradient), weconsidertwoscenarios
thatcanconstrainsuchaplanet’sproperties.
13
Recent
observationsofFomalhautwiththeJWSTNIRCaminstrument
placeconstraintsonplanetmassesofM
pl
≲200M
⊕
beyond5″
(38.3au)fromthestellarcenter(M.Ygoufetal.2024). Our
analysisofsingle-planetscenariosareconsistentwithsucha
nondetection,thoughwefindtheplausibleplanetmassrange
wellbelowthisNIRCamlimit.Perhapsmoreimportantfor
constrainingplanetarypropertiesaretheJWSTMIRIobserva-
tions(A.Gáspáretal.2023). Inthese,thefirstevidenceofan
“intermediatebelt”ispresented,whichhasinnerandouter
edgesof83auand104au,respectively,withlargecorresp-
ondingeccentricitiesof0.31and0.265.Weutilizethese
findingshere,ontheassumptionthataplanetcouldnotreside
withintheregionspannedbytheintermediatebelt,butinstead
mustbelocatedeitherbetweentheALMA-observed“main
belt”andtheintermediatebelt,orinteriortotheintermediate
belt.Wepresentthesingle-planetproperties(derivedin
Sections4.1.1and4.1.2) inTable2.
4.1.1.AGap-carvingPlanet?
Wefirstdiscussthepossibilitythatasingleplanetfully
carvedthegapuptothedisk’sinneredgedirectly.Inthis
scenario,theALMAobservationscanplaceconstraintson
suchaplanet’ssemimajoraxisandeccentricity,using
Equation(10)ofA.J.Mustill&M.C.Wyatt(2012), i.e.,( ) ( )/
/
μ
= ×a aa e1.8
inner plpl
inner
15
.Foraplanetwitha
massrangingfromM
pl
≈1–16M
⊕
(orμintherange
1.7×10
−6
to2.5×10
−5
giventhe2M
⊙
massofFomalhaut,
themassrangeforwhichwejustifybelow), thisexpression
requiresplanetsemimajoraxesintherangea
pl
=109–115au
andeccentricitiesintherangee
pl
=0.18–0.20,withthe
innermostsemimajoraxisa
pl
inthisrangecorrespondingtothe
highestvaluesofe
pl
andM
pl
.Wecalculatethesevaluesbased
ontheinner-edgedisksemimajoraxis(ofa
inner
∼125au)and
adiskeccentricityatthislocation(e
inner
≈0.16,estimatedon
anextrapolationfromthediskcenterofthefittedvaluesofe
andn
pow
).
Wejustifytherangeofplanetmasseswiththelower-limit
correspondingtothetimescaleonwhichagivenplanetcan
plausiblysecularlyforceitseccentricityintothedisk,andat
theupperlimitbasedondirectplanetesimalclearing.For
example,applyingEquation(15)ofA.J.Mustill&
M.C.Wyatt(2009) withtherangeofderiveda
pl
ande
pl
andknownstellarmass,semimajoraxesoutto150aurequire
atleastaμ≈1.7×10
−6
bodytobesecularlyforcedwithin
theageoftheFomalhautsystem.Inaddition,Equation(10)of
A.J.Mustill&M.C.Wyatt(2012) suggeststhatplanetsat
109aurequiremassesofμ≈2.5×10
−5
toforcetheinner-
edgediskeccentricity.UtilizingS.Morrison&R.Malhotra
(2015), i.e.,thatplanetsclearmaterialwithinΔa=1.2μ
0.31
a
pl
oftheirsemimajoraxis,suchabodywouldclearmaterial
inwardsto104au,i.e.,theintermediatebeltouteredge,thus
wesetamassupperlimitbasedonthislocation.Wereferto
thisasthe“Gap”planetinTable2.Wenotethatifthe
intermediatebeltisnotalong-livedplanetesimalbelt,planet
massconstraintsfromJWSTNIRCam(i.e.,M
pl
≲200M
⊕
)
wouldallowplanetsemimajoraxesandeccentricitiestospana
slightlybroaderrangeof100–120auand0.18–0.23,
respectively.
Theplanetorbitalparametersofthis“Gap”planetallagree
withthosepresentedforFomalhautb(P.Kalasetal.2008).
AlthoughmorerecentobservationsandanalysisofFomalhaut
bimplythisisnolongerpresent(A.Gaspar&G.Rieke2020;
M.Ygoufetal.2024), thisanalysissuggeststhataplanetwith
propertiesconsistentwithFomalhautbmaybepresentinthe
system.
4.1.2.AResonant-clearingPlanet?
ThepresenceofFomalhaut’sintermediatebeltpresents
anotherplausiblescenarioinwhichthegapbetweenthemain
beltandtheintermediatebeltwasclearedbyaninternal
planet’sresonantinteractions.Suchaplanetwouldneedto
beinternalto83aubasedonthee=0.31intermediatebelt
inneredge,andhaveastrongresonancelocatedbetweenthe
inneredgeofthemainbeltandtheouteredgeofthe
intermediatebelt.Foraplanet’s(strongest) 2:1meanmotion
resonancetoresidewithin2aufromthemidpointofthese
beltlocations,thesemustspansemimajoraxesof70–75au.
Thissemimajoraxisrangefurtherrequiresa(narrower) range
ofplanetmassesof7×10
−6
≲μ≲2.5×10
−5
forthese
planetstoberesponsiblefortheeccentricityattheinneredge
oftheintermediatebelt(i.e.,weutilizeoncemore( ) ( )/
/
μ
= ×a aa e1.8
inner plpl
inner
15
,withe
inner
=0.31). If
thesebodiesalsofollowthepower-laweccentricitygradient,
thenthesewouldhaveeccentricitiesof0.38–0.41.Wereferto
thisasthe“Resonant”planetinTable2.
Wenotethatinthisresonant-clearingscenario,ifsucha
planetwasresponsiblefordrivingtheeccentricitydistribution
oftheouterbelt(whichwemeasuredase
f
∝a
−1.7
),thenit
couldfeasiblydrivelargereccentricitiesintheintermediate
beltgiventhesewouldbemuchcloserintheirsemimajoraxes.
Byextrapolatingthefittedeccentricitypowerlawtosemimajor
axesof83–104au(i.e.,thoseoftheintermediatebelt), we
calculateeccentricitiesintherange0.22–0.32,whichoverlap
reasonablywellwiththevaluesfittedtotheinnerandouter
edgesoftheJWSTobservationsof0.31and0.265,
respectively(A.Gáspáretal.2023). Whilenotstatedexplicitly
intheworkofA.Gáspáretal.(2023), theirJWSTanalysis
corroboratesthefindingthatFomalhauthostsaneccentricity
distributionthatfallswithsemimajoraxis(giventhedecrease
ineccentricityfromtheintermediatebeltinnertoouteredge).
Table2
DerivedConstraintsforaSinglePlanetResponsibleEitherforDirectly
CarvingtheInner-edgeGapoftheMainBelt(“Gap”) orClearingViaa2:1
Resonance(“Resonant”)
Scenario a
pl
e
pl
M
pl
(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
Foradiscussionoftheparametersonemayderiveforplanetsorbiting
Fomalhautatinclinationsmisalignedwiththedisk,wereferthereaderto
T.D.Pearceetal.(2021).
8
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.

9. 5.TheOriginofFomalhaut’sEccentricDisk
5.1.N-bodyModelingoftheFomalhautDebrisDisk
Inthissection,weinvestigatetheplausibleoriginsofthe
diskeccentricity,byconsideringthestabilityofthediskover
thefull440Myrageofthestar(E.E.Mamajek2012). We
basethissectionontheoutcomeofourearliermodelingand
theplausible(single) planetscenariosconsistentwithobserva-
tions.Inallcases,weconsiderascenariowithasinglestar
(FomalhautA)orbitedbyasingle(gravitating) eccentric
planet,ofsubstantiallylowermass.Todeterminethestability
oftheFomalhautdebrisdisk,weconsidertheinteractionof
thisstar–planetsystemwithapopulationoftestparticles(the
effectsofthediskself-gravityareexploredinSection5.2),
whichsamplethehypotheticalbirthorbitalelementdistribu-
tionoftheFomalhautmainandintermediatebelts.
5.1.1.REBOUNDSetup
ToinvestigatethestabilityoforbitsintheFomalhaut
system,weuseREBOUND(H.Rein&S.F.Liu2012). We
adopttheWHFASTsymplecticintegrator(H.Rein&
D.Tamayo2015) withafixedtimestepsetto1%ofthe
initiallyinnermostparticlesorbitalperiod,
14
usingdemocratic
heliocentriccoordinates(i.e.,wherepositionsaremeasured
relativetothestar,andthecorrespondingcanonicalmomenta
aremeasuredrelativetothebarycentre;H.Rein&
D.Tamayo2019). Intotal,weconsidersixscenarios(that
investigatethreeplanetsetups,andtwodisksetups) thatare
summarizedinTable3anddescribedbelow.
Weexploreboththe“Gap”and“Resonant”planetscenarios
(asperTable2)foreitherinitiallycircularorinitiallyeccentric
testparticledisks.FortheResonantplanet,weconsidera
singlemassofM
pl
=2.5×10
−5
M
*
,asemimajoraxisof
a
pl
=70.87au,andaneccentricitye
pl
=0.4.Scenarioswitha
“1”intheirIDincludethisplanet.FortheGapplanet,we
considertwomassesofM
pl
=1.0×10
−6
M
*
and
M
pl
=1.0×10
−7
M
*
,bothwithasemimajoraxisof
a
pl
=112.9au,andeccentricitye
pl
=0.215.Scenarioswith
eithera“2”or“3”intheirIDincludethisplanet.Adopting
unitsinwhichGM
*
=1andunitlengthequaltotheplanets
semimajoraxis,eachplanethasanorbitalperiod≈2π(strictly
theorbitalperiodisslightlyshorterduetothefinite
planetmass).
Weconsidertwoinitialscenarios,eitherwithallparticles
initiallycircular,orinitiallyeccentric.Inthecaseofthe
initiallycirculartestparticles,wedistributethesewithe=0.
Scenarioswitha“C”intheirIDrelatetothisinitialdisksetup.
Inthecaseoftheinitiallyeccentrictestparticles,wedistribute
thesewith( )/=eeaa
pl pl
1
(i.e.,consistentwithourmodel-
ing).Scenarioswithan“E”intheirIDrelatetothisinitialdisk
setup.Wesetupthediskwith2×10
5
testparticlesuniformly
distributedbetween( )/a70 au
pl
1
and( )/a180 au
pl
1 .Wenote
thatthisbroadradialdistributionspansallradiiwhere
Fomalhaut’smainbeltandintermediatebelthavebeen
observed.Inallcases,werandomlydistributetheparticles
uniformlyinsemimajoraxis-meananomalyspace.Toavoid
issueswiththeplanetarypointmasssingularity,weadopt
REBOUND’s gravitationalsofteningwithagravitationalsoft-
eningparametersettob=10
−6
[a
pl
].Thischoiceofsoftening
radiusiswellwithintheplanet’sHillsphereforallscenarios
considered.
5.1.2.REBOUNDResults
WerunsimulationsE1,E2,E3,C2,andC3forthefull
440MyrlifetimeofFomalhaut.WerunsimulationC1forjust
1Myrbecausebythistimethediskisalreadydisrupted(inthe
sensethatlargesectionsofthediskhavebeenbrokenupor
clearedbyinteractionwiththeplanet). Weshowtheoutcomes
ofthesesixscenariosinFigure4,firstintermsoftheir
“conditional”masssurfacedensitydistributions(i.e.,the
simulationoutcomes) andsecondtoaidinterpretationbecause
theirsurfacedensitydistributionsconvolvedwithaspecified
radialprofile.Theseconvolveddistributionsconvolvethe
resultsoftheREBOUNDsimulations,withamassperunit
semimajoraxisof( )
()
=M
a
w
aa
w
2
exp
2
, 8
a
r
r
0
2
2
consistentwiththemodelinginSection3.Presentingthedata
inthismannercomeswithimportantcaveats.First,thisprofile
choiceeffectivelymasksallemissionsnotimagedbyALMA.
Thus,whilethesimulationscaninformusofthestabilityof
orbitsinFomalhaut’sintermediatebelt,werestrictour
investigationherepurelytoFomalhaut’smainbelt.By
extension,thisfurtherdownplaystheimportanceofthe
planet’sgapcarvingsincethishidesthefactthatlower-mass
planetsareunabletocarvedeepgapsintheparticle
Table3
AllSixREBOUNDScenariosSimulatedinthisWork,StartingwithDiskProperties(DiskExtentandEccentricity) andPlanetProperties(Planet-stellarMassRatio
andSemimajorAxis), theStabilityoftheDisk,andCommentsonAnyObservedFeaturesPresent
ID
DiskSemimajorAxis
Extent
DiskEccentricity
Profile PlanetMassRatio
PlanetSemi-
majorAxis PlanetEccentricity Comments
(au) (au)
C1 70–180 e=0 2.5×10
−5
70.87 0.4 Diskdisrupted
C2 70–180 e=0 1.0×10
−6
112.9 0.215 Diskdisrupted
C3 70–180 e=0 1.0×10
−7
112.9 0.215 Diskdisrupted
E1 70–180 e∝a
−1
2.5×10
−5
70.87 0.4 Gapcarved,gapinmainbelt,
spirals
E2 70–180 e∝a
−1
1.0×10
−6
112.9 0.215 Gapcarved,weakspirals
E3 70–180 e∝a
−1
1.0×10
−7
112.9 0.215 Gapcarved,butveryshallow
14
SuchatimestepchoiceisconsistentwithREBOUNDsimulationsofdebris
disksintheliterature,e.g.,J.Dongetal.(2020) andT.D.Pearceetal.(2024),
andconsistentwiththeadviceprovidedintheREBOUNDAPIDocumentation;
seecommenton“Timestepping”inhttps://rebound.readthedocs.io/en/latest/
simulationvariables/.
9
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.

10. distributions(oneshouldthusconsiderbothplotsper
scenario). Second,duetoourignoranceabouttheinitial
conditionsoftheFomalhautdisk,thechosenmassperorbitis,
bynecessity,arbitrary.Themassperorbitisnotaffectedby
secularperturbations,andso,awayfromresonancesandclose
encounters,theplanethasnowaytomodifythemasson
orbits.Overall,theREBOUNDsimulationsarethusnot
expectedtoreproduceeithertheALMAandJWSTobserva-
tions,butinsteadinformusaboutthestabilityoftheorbits
withinthebeltandidentifyfeatures(e.g.,gapsorspiraldensity
waves) thatmaybeexpectedtobeobservedforgiven
planetaryproperties.
Withthesecaveatsinmind,wefind(inallcases) that
scenariosinitializedwithcirculardisksdonotproducestable
eccentricbelts.InthethreecasesC1,C2,andC3,wefindthat
theseallendtheirsimulationswithdisruptedparticle
distributions(inconsistentwithobservations). Incontrast,we
findthatscenariosinitializedwitheccentricdisks(E1,E2,and
E3)survivethefull440Myrsimulationtimescale,andresult
ineccentricdisksconsistentwiththoseobservedfor
Fomalhaut.Inaddition,boththeResonantandGapplanet
scenarioscarvegapsinthelocationbetweenthemainand
intermediatebelts(inthelowestmassplanetcase,thisgapis
relativelyshallowincomparisontothemoremassiveplanet
scenarios). OneintriguingfeatureintheE1scenarioisthe
formationofaspiraldensitywave(duetothemassive
Resonantplanet’sinteractionwiththedisk). Whilenosuch
featureshavebeenobservedintheFomalhautdiskasyet,
deeper,higher-resolutionobservationscouldbeusedto
investigatetheplausibilityofthisscenarioandconstrainthe
massofaplanetwiththeseorbitalparameters.
Overall,theREBOUNDsimulationslendweighttotwo
importantconclusions.First,thesimulationssuggestthat
Fomalhaut’sdiskmayneedtohavebeenformedasan
eccentricdiskinordertosurvive∼440Myrevolution(i.e.,
formationwithintheprotoplanetarydisk). Second,andby
extension,whileasingleplanetcouldberesponsiblefor
carvingtheinneredgeoftheFomalhautmainbelt(orthegap
betweentheintermediateandmainbelts), thissameplanet
maynotbecapableofforcingitsowneccentricityintoan
initiallycircularplanetesimalbeltviasecularplanet–disk
interactionstoyieldaneccentricdiskwithFomalhaut’s
properties.InthecaseoftheGapplanet,thisbodycan0
50
100
150
200
Radius
[
au
]
Disk att= 0:e=0
Planet:apl=71au,=2:510
5
,epl=0:40
C1 Disk att= 0:e=0
Planet:apl=113au,=10
6
,epl=0:215
C2 Disk att= 0:e=0
Planet:apl=113au,=10
7
,epl=0:215
Cond
:
surface MDF
[
ctsau

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
Radius
[
au
]
Disk att= 0:e=0
Planet:apl=71au,=2:510
5
,epl=0:40
0 100 200 300
Angle[deg]
Disk att= 0:e=0
Planet:apl=113au,=10
6
,epl=0:215
0 100 200 300
Angle[deg]
Disk att= 0:e=0
Planet:apl=113au,=10
7
,epl=0:215
Surface density
[
ctsau

2
]
0e+00
1e+01
2e+01
3e+01
4e+01 0
50
100
150
200
Radius
[
au
]
Disk att= 0:e/a
1
Planet:apl=71au,=2:510
5
,epl=0:40
E1 Disk att= 0:e/a
1
Planet:apl=113au,=10
6
,epl=0:215
E2 Disk att= 0:e/a
1
Planet:apl=113au,=10
7
,epl=0:215
Cond
:
surface MDF
[
ctsau

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
Radius
[
au
]
Disk att= 0:e/a
1
Planet:apl=71au,=2:510
5
,epl=0:40
0 100 200 300
Angle[deg]
Disk att= 0:e/a
1
Planet:apl=113au,=10
6
,epl=0:215
0 100 200 300
Angle[deg]
Disk att= 0:e/a
1
Planet:apl=113au,=10
7
,epl=0:215
Surface density
[
ctsau

2
]
0e+00
1e+01
2e+01
3e+01
4e+01
Figure4.PlotsforthesurfacedensitymapsfromtheREBOUNDsimulationsinr–fspace.Uppertworows:Conditionalsurfacemassdensityfunctionandabsolute
surfacedensitysimulationswithinitiallycirculardiskconditions.Lowertworows:Simulationoutcomeswithinitiallyeccentricdiskconditionspresentedlikewise.
Specificplanetanddiskassumptionsarepresentedoneachplot,withblack-dashedlinesshowingtheplanet’sorbit.Initiallycirculardiskconditionsdonotyield
stabledisksat440Myr,whereasscenarioswhichareinitiallyeccentricyieldstableeccentricdisksat440Myr.ThesameinformationistabulatedinTable3.
10
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.

11. plausiblyclearmaterialfromthemainbelt’sinneredge(and
theintermediatebelt’souteredge) ifthisissituatedinan
initiallyeccentricdebrisbelt,i.e.,viacloseplanet-planetesimal
encounters.InthecaseoftheResonantplanet,thiscould
plausiblyclearmaterialinthesamelocationduetotheoverlap
withits2:1meanmotionresonanceandsubsequentplanete-
simalejections.
5.1.3.ComparisonwithEccentricProtoplanetaryDisks
Theconclusionthateccentricdisksthatsurvivetothestellar
mainsequencemayneedtobeproducedduringthe
protoplanetarydiskphaseisconsistentwiththeconclusion
ofG.M.Kennedy(2020), whoarguedthisfromthe
independentstandpointoftheapparentnarrownessofthe
debrisdisksofFomalhautandHD202628.Thereareindeeda
(small) numberofknowneccentricprotoplanetarydisks,such
asMWC758,HD100546,andIRS48(R.Dongetal.2018;
D.Fedeleetal.2021; H.Yangetal.2023, allofwhichare
hostedby∼2M
⊙
stars), aswellastheeccentriccircumbinary
diskofIRAS04158+2805 (whichishostedbyalower-mass
mid-MSpTbinary;S.M.Andrewsetal.2008; E.Ragusa
etal.2021). TheeccentricdebrisdisksofHD202628and
HD53143arehostedby∼Solar-massstars(V.Faramazetal.
2019; M.A.MacGregoretal.2021), andthoseofHD38206
andHR4796likewiseorbit∼2M
⊙
(A-type) stars(J.Olofsson
etal.2019; M.Boothetal.2021). Whetherthisbiastoward
earlier-type,intermediate-massstarsisphysicalorobserva-
tionalshouldbedeterminedinfuturework.
AplausibleoriginofFomalhaut’seccentricdiskisthatthe
eccentricitycouldhaveresultedfromtheformationofaplanet
(withpropertiesconsistentwiththoseinTable2)withinthe
protoplanetarydisk,wherebyinitialplanet–gasdiskinterac-
tionsdeterminedtheeccentricityprofileofthedisk,andthe
eccentricityoftheplanet.Insuchascenario,aneccentricity
profilewhichdecreaseswithradiusisexpected(e.g.,J.Teys-
sandier&G.I.Ogilvie2016). Thisfindingfollowsfromthe
factthateccentricityprofilesthatareeitherflatorincreasing
withradiustendtodifferentiallyprecessasaresultofgas
pressureforces.InthecaseofFomalhaut,followingthe
dissipationofitsprotoplanetarydiskgas(andanyremnant
primordialdust), subsequentplanet–diskinteractionswiththe
eccentricplanetesimalbeltplausiblysculpteditsdisk’sinner
andouteredges,whichoverthecourseof440Myrresultedin
thebeltmorphology,whilemaintainingitsinitialeccentricity
distribution(oroneverysimilar).
5.2.AdditionalPhysics?TheRoleofDiskSelf-gravity
Wehavethusfarexploredthepossibleoriginof
Fomalhaut’sdiskeccentricityassumingastationaryplanet
perturbingamasslessdebrisdisk.However,additional
physicaleffects,suchastheinfluenceofamassiveandself-
gravitatingdebrisdisk,mayaffecttheoutcomes.
AsmentionedinSection2,aplanetwithinadebrisdisk
induces(secular) forcedplanetesimaleccentricitiessuchthat
e
f
(a)∝a
−1
(C.D.Murray&S.F.Dermott1999). This,
however,appliesspecificallytoamasslessdebrisdisk
(M
d
=0),i.e.,composedoftestparticles.Bycontrast,ifthe
diskismassive,itsself-gravitycansteepentheforced
eccentricityprofile(A.A.Sefilian2024): namely,toasmuch
ase
f
(a)∝a
−4.5
,dependingona
p
/a
in
andthemassdistribution
withinthedisk(seealso,A.A.Sefilianetal.2021). This
occurswhenthedisk-inducedapsidalprecessionratesofthe
planetand/or theplanetesimals—absentinmasslessdisk
models—dominateovertheplanet-inducedprecessionof
planetesimalorbits.Fora
p
∼a
in
,thisconditiongenerally
translatestoM
d
/m
p
≳1,andwouldrequireM
d
/m
p
≲1for
a
p
≲a
in
.Interestingly,underthesameconditionsonM
d
/m
p
,a
similareffectarisesinthedisk’sverticalaspectratiohifthe
planet–disksystemismisaligned:ratherthanremaining
constantwithradius(asexpectedforM
d
=0),hmayinstead
decreasewithdistancefromthestarifM
d
≠0(A.A.Sefilian
etal.2025). Thiseffect,however,isnotaccountedforinthis
study.
Withoutpresumingspecificplanetmasses,ALMAobserva-
tionsimplyatotaldiskmass(i.e.,includingthelargest,
collidingplanetesimals) consistentwiththerangeinferredby
thecollisionalmodelingofA.V.Krivov&M.C.Wyatt
(2021), whichestimatesM
d
=1.8–360M
⊕
.Iftrue,disk
gravitycouldpotentiallyexplainthesteepe
f
(a)profile
indicatedbyourmodeling.Whileitmaybetemptingtotest
thishypothesisusingtheresultsofA.A.Sefilian(2024), their
studyisaproof-of-conceptbasedonasemianalyticalframe-
workthataccountsfortheaxisymmetriccomponentofthedisk
potential,butignoresthenonaxisymmetriccounterpart.In
principle,thislimitationcanbeaddressedusingexistingorbit-
averaged,secularcodes(e.g.,thelinear,N-RINGcodeof
A.A.Sefilianetal.2023orthemoreaccurateGAUSScodeof
J.R.Toumaetal.2009) and/or directN-bodycodes;however,
thisismoresuitedtoaseparatestudy(SefilianA.A.etal.
2025,inpreparation).
6.Conclusions
WepresentanewmodelofFomalhaut’seccentricdebris
diskbyfittingparametricdebrisdiskmodelstoarchival
ALMAobservations(withasynthesizedphysicalresolution
4–6au).Themodelsaredevelopedfromthoseintroducedby
E.M.Lynch&J.B.Lovell(2021) andJ.B.Lovell&
E.M.Lynch(2023), wherebyagradientintheforced
eccentricityparameter,()eaa
n
pow
,isintroducedtofitthe
knowndiskasymmetries:abrighterapocenterversuspericen-
ter,andabroaderpericenterversusapocenter.Wecollectively
termthisphysicaleffectwithindisksEVD.Thebest-fit
parameterdeterminationsfromthemodelingfindparameters
broadlyconsistentwiththosefromothermodels,butwiththe
inclusionofasteep,negativeeccentricitygradientasper
e∝a
−1.75±0.16
.Byanalyzingthegoodnessoffit,wedeem
thatsuchadescriptionprovidesanovelinterpretationofthis
debrisdisk,andisstatisticallypreferredtomodelswith
constantfreeandforcedeccentricitydistributionsthroughthe
disk.Toourknowledge,thisisthefirstreportedeccentricity
gradientinadebrisdisk.Thevaluewederiveisverycloseto
thatexpectedfromclassicalexpectationsof(massless) disk–
planetinteractionse
f
∝a
−1
,suggestingthegradientmaybe
forcedbyaplanetaryperturber.
BasedontheALMAmodeling,andpublisheddatafrom
JWSTNIRCamandMIRI,wequantifyplanetaryorbitaland
massconstraintsfortwocoplanarsingle-planetscenarios
consistentwithobservations.Onescenariodescribesa
109–115auplanetthathasdirectlyclearedmaterialuptothe
inneredgeofFomalhaut’s“mainbelt”(asimagedbyALMA).
Thesecondpositedscenariodescribesa70–75auplanetthat
hasclearedthedisk’sinneredgeatits2:1meanmotion
resonance,withthisplanetsittinginteriortotheJWST-imaged
11
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.

12. “intermediatebelt.”Inbothcases,theimpliedplanetmassand
semimajoraxisrangesarebelowsensitivitythresholdsfor
existingplanetdetectionmethods.
WeassessthestabilityoftheFomalhautdebrisdiskby
conductingN-bodyREBOUNDmodeling.Wesetupsixsetsof
initialconditions,basedonthreeplanets,andtwoinitial
conditionsforthedisk,withthiseitherbeingcircularor
eccentric.Weshowthatonlyscenarioswithaninitially
circulardiskarestableandfinalizeaseccentricdisksafterthe
440Myrsimulationlifetime(matchingthatoftheageof
Fomalhaut). Thesefindingsmaysuggestthatplanet–disk
interactionsareprimarilyresponsibleforsculptingthedisk’s
morphology(i.e.,itsinner-edges,andas-pertheJWST
observations,gapsinthedisk), butnotitseccentricity,and
thusthatFomalhaut’seccentricringwasplausiblyborn
eccentric.Wediscusscaveatswiththesemodels(e.g.,the
exclusionofdiskself-gravity) andproposefurthersimulations
toinvestigatetherobustnessofthesefindingsformassive
debrisdisks.
Finally,wehighlightthatthecodeusedtomodel
Fomalhaut’ssurfacedensitydistributioninthiswork
(Equations(1)–(4)) isavailablepubliclyonhttps://github.
com/astroJLovell/eccentricDiskModels. Wereleasethisto
enableourworktobereproducedeasily,andsuchthatothers
canutilizethistomodelthesurfacedensityprofilesofother
opticallythindisks.
Acknowledgments
Wethanktheanonymousrefereefortheircommentsand
suggestions,whichgreatlyimprovedtheclarityofour
investigation.J.B.L.acknowledgestheSmithsonianInstitute
forfundingviaaSubmillimeterArray(SMA) Fellowship,and
theNorthAmericanALMAScienceCenter(NAASC) for
fundingviaanALMAAmbassadorship.J.B.L.dedicatesthis
studytothememoryofHenryChesters,withwhomhe
discussedthiswithgreatinterestintheirfinalconversation,
andforH.C.’slifelongencouragementofhisscientific
endeavors.A.A.S.issupportedbytheHeising-Simons
Foundationthrougha51PegasibPostdoctoralFellowship.
Weacknowledgetheoperationalstaffandscientistsinvolved
inthecollectionofdatapresentedhere.
Facility:ALMA.
Software:Thisresearchmadeuseofthefollowingsoftware
packages:CASA(J.P.McMullinetal.2007) DS9(Smithsonian
AstrophysicalObservatory2000). emcee(D.Foreman-Mackey
etal.2013); RADMC-3D(C.P.Dullemondetal.2012);
REBOUND(H.Rein&S.F.Liu2012).
SoftwareandThirdPartyDataRepositoryCitations
ThispapermakesuseofthefollowingALMAdata:ADS/
JAO.ALMA2013.1.00486.S,ADS/JAO.ALMA 2015.1.00966.
S,andADS/JAO.ALMA 2017.1.01043.S.ALMAisapartner-
shipofESO(representingitsmemberstates), NSF(USA),
andNINS(Japan), togetherwithNRC(Canada), MOSTand
ASIAA(Taiwan), andKASI(RepublicofKorea), incooperation
withtheRepublicofChile.TheJointALMAObservatoryis
operatedbyESO,AUI/NRAO, andNAOJ.TheNationalRadio
AstronomyObservatoryisafacilityoftheNationalScience
FoundationoperatedundercooperativeagreementbyAssociated
Universities,Inc.
AppendixA
DataMasking
WepresentinFigure5theresultoffittingtheGeneralmodel
whileonlymaskingpixelswithina200au(projected) distance
fromthestar,i.e.,withoutmaskingbytheprimarybeam.We
showthespatialmaskintheleftpanelofthisfigure(ingreen
dashed–dottedline),andalsoshowtheprimary-beamresponseat
the0.66level(inorangedashed–dottedlines)whichdefinesthe
innerboundaryofanyfitsthatincorporateaprimary-beammask.
Theregionsofthetwoansaebetweenthegreenandorange
dashed–dottedlinesdefinethemaskusedinthemainbodyofthis
work.Itisevidentthatalongthedisk’sminoraxis,theprimary-
beamresponseissignificantlyreducedincomparisontothedisk
ansae,andwhiletheminoraxesappeartobeslightlybetterfitin
comparisontotheGeneralmodelwithaprimary-beammask,the
fitissignificantlyworseattheapocentre,justifyingthechoiceto
fitthismodelonlyatthetwodiskansae.
AppendixB
ModelVerification
Weverifythemodelingandfittingmethodologybyfittinga
six-parametermodeltothearchivalALMACycle3-onlydata
forFomalhaut,firstpresentedinM.A.MacGregoretal.(2017)
andsubsequentlyinG.M.Kennedy(2020). Werefertothis
low-resolutionmodelas“Verify.”Wefindthattheparameter
distributionsandmeanvaluesofthefittingareconsistentwith
thosepreviousworks(i.e.,M.A.MacGregoretal.2017;
G.M.Kennedy2020), albeitwithaslightlylarger(mean)
forcedeccentricity(whichismostlikelyduetothedifferencein
eccentricityprofilewefithere). Suchconsistencyprovides
assurancethatthefittingmethodologyweconducthereprovides
reasonableuncertaintiesonthefittedparameters.Wepresentin
Figure6theposteriordistributionafterremovingburn-in
chains,thebest-fitvaluesonTable1,andthebest-fitmodeland
residualsinFigure7.Mostinteresting,isthatwefindanegative
eccentricitygradientislikewiserequiredwiththismodeltofit
eventhelower-resolutiondata,withn
pow
=−1.16±0.16,and
thatthismodelsetupcanaccountforthedisk’swidthand
brightnessasymmetries(presentedbyG.M.Kennedy2020).
Oneimprovementisthatthemodelwepresentherefully
accountsforthenorthwestpositiveemissionthatisstillpresent
intheapocenterresidualsofG.M.Kennedy(2020). In
verifyingthisfittingmethodology,wethusfindthatevidenceof
aneccentricitygradientinFomalhaut’sdiskexistedpriortoits
observationsathigher-resolution.-10.00.010.0
-20.0
-10.0
0.0
10.0
20.0
Relative Decl. Oset [arcsec]
20
0
20
40
60
80
100
120
-10.00.010.0
Relative RA Oset [arcsec]
ef/a
1:75
20
0
20
40
60
80
100
120
-10.00.010.0
Intensity [

Jy beam

1
]
40
30
20
10
0
10
20
30
40
Figure5.Left:ALMAdata.Center:Best-fite
f
(a)∝a
−1.75
model(onsame
imagescale). Right:Residualemissionafterbest-fitmodelsubtraction.Contours
areshownatthe±3σlevelintheresidualsandatthe5σlevelforthedata/model.
Beamsareshowninthelower-leftofeachplot.Intheleftandrightplotsinamber
dashed–dottedlines,weoverplottheprimary-beamresponseatthe0.66level
(andinadditionatthe0.4levelontherightplot).Weshowintheleftplotinthe
greendashed–dottedlinethe200auspatialmaskdiscussedinAppendixA.
12
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.

13. Mdust = 0.0199
+0.0002
0.0002
138.6
138.8
139.0
139.2
a0
a0 = 138.8873
+0.0983
0.0993
14.4
15.0
15.6
16.2
16.8
wr
wr = 15.2977
+0.3198
0.3140
0.1300
0.1325
0.1350
0.1375
0.1400
ef,0
ef,0 = 0.1346
+0.0015
0.0014
18.0
19.5
21.0
f
f = 19.1031
+0.6276
0.6589
0.01920.01950.01980.02010.0204
Mdust
1.50
1.25
1.00
0.75
npow
138.6138.8139.0139.2
a0
14.415.015.616.216.8
wr
0.13000.13250.13500.13750.1400
ef,0
18.0 19.5 21.0
f
1.50 1.25 1.00 0.75
npow
npow = 1.1651
+0.1613
0.1572 Figure6.Cornerplotsfortheemceechains(minus1000stepscoveringburn-in) showingthesixparametersfittedtothelow-resolutiondataofM.A.MacGregor
etal.(2017) andG.M.Kennedy(2020), i.e.,“Verify”model.WefindcomparableparameteruncertaintiestothosepresentedinG.M.Kennedy(2020) andalower
valueofthermsintheresidualmapthanthe“uniform”modelpresentedbyG.M.Kennedy(2020).
13
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.

14. AppendixC
ResidualEmissionMapsandemceePosterior
Distributions
Here,wepresentforthesix-parameter(low-resolution) “Verify”
model,thecornerplotsin6andthedata,model,andresidualmaps
inFigure7.WepresentthecornerplotsfortheGeneral,and
Proper1andProper2modelsinFigures8,9,and10,andthedata,
model,andresidualmapsfortheProper1andProper2modelsin
Figure3.Inallcases,thefirst1000stepswereremovedfromthe
emceechains(farinexcessofthenumberofburn-insteps).-10.00.010.0
-20.0
-10.0
0.0
10.0
20.0
Relative Decl. Oset [arcsec]
50
0
50
100
150
200
250
300
350
400
-10.00.010.0
Relative RA Oset [arcsec]
ef/a
1:16
50
0
50
100
150
200
250
300
350
400
-10.00.010.0
Intensity [

Jy beam

1
]
40
30
20
10
0
10
20
30
40
Figure7.Left:ALMAdataaspresentedinM.A.MacGregoretal.(2017) andG.M.Kennedy(2020). Center:Best-fitmodelforthisdata.Right:Residualemission
(dataminusbest-fitmodel). Contoursareshownatthe±3σlevelintheresiduals,andatthe5σlevelforthedata/model. Beamsareshowninthelower-leftcornerof
eachplot.Inallplots,northisupandeastisleft.Thestrongpeakat(−4″, −3″)isthesourceGDC,asnotedinbothA.Gáspáretal.(2023) andG.M.Kennedy
etal.(2023).
14
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.

15. Mdust = 0.0187
+0.0002
0.0002
138.4
138.6
138.8
139.0
139.2
a0
a0 = 138.7938
+0.1209
0.1193
12.5
13.0
13.5
14.0
14.5
wr
wr = 13.5075
+0.2840
0.2927
0.1225
0.1250
0.1275
0.1300
ef,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
npow
npow = 1.7458
+0.1577
0.1567
0.0100
0.0125
0.0150
0.0175
0.0200
h
h = 0.0157
+0.0014
0.0013
335.95
336.10
336.25
336.40
PA
PA = 336.1927
+0.0568
0.0577
0.01800.01840.01880.01920.0196
Mdust
66.2
66.4
66.6
66.8
inc
138.4138.6138.8139.0139.2
a0
12.513.013.514.014.5
wr
0.12250.12500.12750.1300
ef,0
13.515.016.518.0
f
2.4 2.1 1.8 1.5 1.2
npow
0.01000.01250.01500.01750.0200
h
335.95336.10336.25336.40
PA
66.266.466.666.8
inc
inc = 66.4358
+0.0916
0.0898 Figure8.Cornerplotsfortheemceechains(minusburn-insteps) showingthenineparametersfittedinthisstudyfortheGeneralmodel.Allshowthewell-behaved
featuresofGaussiandistributions,whichareeithercircularorwithlittledegeneracy(e.g.,betweenf
M,dust
andw
r
,andbetweeneandω
f
).
15
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.

16. Mdust = 0.0209
+0.0015
0.0013
138.50
138.75
139.00
139.25
r0
r0 = 138.9286
+0.1176
0.1182
8
9
10
11
wr
wr = 9.3401
+0.5830
0.5634
0.122
0.124
0.126
0.128
0.130
e
e = 0.1256
+0.0011
0.0011
13.5
15.0
16.5
f
f = 14.9846
+0.6473
0.6317
0.028
0.032
0.036
0.040
0.044
ep
ep = 0.0389
+0.0015
0.0017
0.0100
0.0125
0.0150
0.0175
0.0200
h
h = 0.0156
+0.0014
0.0014
336.1
336.2
336.3
336.4
PA
PA = 336.2189
+0.0541
0.0529
0.01750.02000.02250.02500.0275
Mdust
66.15
66.30
66.45
66.60
66.75
i
138.50138.75139.00139.25
r0
8 9
10 11
wr
0.1220.1240.1260.1280.130
e
13.5 15.0 16.5
f
0.0280.0320.0360.0400.044
ep
0.01000.01250.01500.01750.0200
h
336.1336.2336.3336.4
PA
66.1566.3066.4566.6066.75
i
i = 66.4434
+0.0899
0.0910 Figure9.Cornerplotsfortheemceechains(minusburn-insteps) showingthenineparametersfittedfortheProper1model.Allshowthewell-behavedfeaturesof
Gaussiandistributions,whichareeithercircularorwithlittledegeneracy(e.g.,betweenf
M,dust
andw
r
,andbetweeneandω
f
).
16
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.

17. AppendixD
ModelConvergence
Wefindbyfittinginitialroundsofemcee(fortheninefree
parameterGeneraldiskmodel) thattheburn-instageis
approximatelyafewhundredstepsandthattheautocorrelation
time(τ)is≈120steps.FortheGeneralandProper1models,we
measuretheautocorrelationfunction(τ)versussteplengthfor
allchains,andfindthisvalueasymptotes(withnodeviationby
morethan1%)toaparameter-meanτof≈120.Inthecaseof
theProper2model,wemeasureameanτof≈150,andforthe
Verifymodelameanτof≈60.Asstatedintheemcee
guidance(D.Foreman-Mackeyetal.2013), convergenceis
achievedwhenthenumberofstepsexceeds>50×τ.Assuch,
weruntheVerify,General,andProper1andProper2models
for3000,6000,6000,and7500steps,respectively,andverify
thatconvergenceisindeedmetinallcases.Inallcases,wefind
after300–500stepsthatthechainsareburnedin.Toderivethe
best-fitposteriordistributionspresentedinTable1,weremoved
thefirsthalfofthetotalstepsfromallchains.Weused40
walkersinallemceeruns,whichisfarinexcessofthe
necessary2×thenumberoffreeparametersi.e.,12forthe
Verifymodeland18fortheotherthreemodels.
ORCIDiDs
JoshuaB.Lovell
aahttps://orcid.org/0000-0002-4248-5443
JayChittidi
aahttps://orcid.org/0000-0002-4985-028XMdust = 0.0217
+0.0015
0.0013
138.50
138.75
139.00
139.25
r0
r0 = 138.8482
+0.1205
0.1172
10.5
12.0
13.5
wr
wr = 12.7650
+0.6997
0.9772
0.122
0.124
0.126
0.128
0.130
e
e = 0.1257
+0.0011
0.0012
13.5
15.0
16.5
f
f = 15.1589
+0.6316
0.6213
0.008
0.016
0.024
0.032
ep
ep = 0.0202
+0.0079
0.0116
0.0100
0.0125
0.0150
0.0175
h
h = 0.0150
+0.0013
0.0014
336.0
336.1
336.2
336.3
336.4
PA
PA = 336.2032
+0.0575
0.0542
0.0180.0210.0240.0270.030
Mdust
66.2
66.4
66.6
66.8
i
138.50138.75139.00139.25
r0
10.512.013.5
wr
0.1220.1240.1260.1280.130
e
13.5 15.0 16.5
f
0.0080.0160.0240.032
ep
0.01000.01250.01500.0175
h
336.0336.1336.2336.3336.4
PA
66.266.466.666.8
i
i = 66.4255
+0.0907
0.0899
Figure10.Cornerplotsfortheemceechains(minusburn-insteps) showingthenineparametersfittedfortheProper2model.Allshowthewell-behavedfeaturesof
Gaussiandistributions,whichareeithercircularorwithlittledegeneracy(e.g.,betweenf
M,dust
andw
r
,andbetweeneandω
f
).
17
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.

18. AntranikA.Sefilian aahttps://orcid.org/0000-0003-
4623-1165
SeanM.Andrews
aahttps://orcid.org/0000-0003-2253-2270
GrantM.Kennedy
aahttps://orcid.org/0000-0001-6831-7547
MeredithMacGregor
aahttps://orcid.org/0000-0001-
7891-8143
DavidJ.Wilner
aahttps://orcid.org/0000-0003-1526-7587
MarkC.Wyatt
aahttps://orcid.org/0000-0001-9064-5598
References
Acke,B.,Min,M.,Dominik,C.,etal.2012,A&A, 540,A125
Allen,R.H.1963,StarNames.TheirLoreandMeaning(NewYork: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.,etal.2012,ApJL, 750,L21
Booth,M.,Schulz,M.,Krivov,A.V.,etal.2021,MNRAS, 500,1604
Chittidi,J.,MacGregor,M.A.,Lovell,J.B.,etal.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.,etal.2018,ApJ,860,124
Dullemond,C.P.,Juhasz,A.,Pohl,A.,etal.,2012AstrophysicsSourceCode
Library,ascl:1202.015
Faramaz,V.,Krist,J.,Stapelfeldt,K.R.,etal.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.,etal.2023,NatAs, 7,790
Goodman,J.,&Weare,J.2010,CAMCS, 5,65
Holland,W.S.,Greaves,J.S.,Dent,W.R.F.,etal.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.,etal.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.,etal.2021,MNRAS, 506,1978
Lovell,J.B.,Wyatt,M.C.,Kalas,P.,etal.2022,MNRAS, 517,2546
Lynch,E.M.,&Lovell,J.B.2021,MNRAS, 510,2538
MacGregor,M.A.,Matrà,L.,Kalas,P.,etal.2017,ApJ,842,8
MacGregor,M.A.,Weinberger,A.J.,Loyd,R.O.P.,etal.2021,ApJL,
911,L25
Mamajek,E.E.2012,ApJL, 754,L20
Marino,S.2022,arXiv:2202.03053
Matrà,L.,MacGregor,M.A.,Kalas,P.,etal.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,
inProtostarsandPlanetsVI,ed.H.Beutheretal.(Tucson,AZ:Univ.
ArizonaPress), 521
McMullin,J.P.,Waters,B.,Schiebel,D.,Young,W.,&Golap,K.2007,in
ASPConf.Ser.376,AstronomicalDataAnalysisSoftwareandSystems
XVI,ed.R.A.Shaw,F.Hill,&D.J.Bell(SanFrancisco,CA:ASP), 127
Morrison,S.,&Malhotra,R.2015,ApJ,799,41
Murray,C.D.,&Dermott,S.F.1999,SolarSystemDynamics(Cambridge:
CambridgeUniv.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.,etal.2019,A&A, 630,A142
Pan,M.,Nesvold,E.R.,&Kuchner,M.J.2016,ApJ,832,81
Pearce,T.D.,Beust,H.,Faramaz,V.,etal.2021,MNRAS, 503,4767
Pearce,T.D.,Krivov,A.V.,Sefilian,A.A.,etal.2024,MNRAS, 527,3876
Ragusa,E.,Fasano,D.,Toci,C.,etal.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.,etal.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
SmithsonianAstrophysicalObservatory,2000SAOImageDS9:AUtilityfor
DisplayingAstronomicalImagesintheX11WindowEnvironment,
AstrophysicsSourceCodeLibrary,ascl:0003.002
Stapelfeldt,K.R.,Holmes,E.K.,Chen,C.,etal.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.,etal.1999,ApJ,527,918
Yang,H.,Fernández-López,M.,Li,Z.-Y.,etal.2023,ApJL, 948,L2
Ygouf,M.,Beichman,C.A.,Llop-Sayson,J.,etal.2024,AJ,167,26
18
TheAstrophysicalJournal,990:145(18pp), 2025September10 Lovelletal.