Comet Evolutionary Diagram (CED)

• CED is a family of diagrammatic frameworks that represent comet evolution using empirical, dynamical, and model-derived schemata.
• The empirical CED employs mass-loss age proxies and remaining returns to classify comets into groups such as young Oort-cloud and evolved, dust-mantled objects.
• The CODE dynamical CED tracks comets through five orbital stages, integrating non-gravitational forces to predict changes in orbital energy and perihelion distance.

Searching arXiv for the specified papers to ground the encyclopedia entry in the cited literature. Search query: (Królikowska et al., 2020) OR (Ferrin et al., 10 Apr 2026) OR (Parhi et al., 30 Oct 2025) OR (Filacchione et al., 2022) OR (Ferrin et al., 2013) OR (Ferrín, 2013) The Comet Evolutionary Diagram (CED) denotes a family of diagrammatic frameworks used to represent comet evolution, but it is not a single standardized construct across the literature. In one major usage, the CED is an empirical, comparative plane whose axes are a mass-loss-based age proxy and the number of remaining returns, intended to organize Oort-cloud comets, Jupiter-family comets, asteroidal-belt comets, disintegrators, and suffocated objects. In another usage, developed for the CODE catalogue of long-period comets, the CED is a five-node dynamical track linking a comet’s previous perihelion, original barycentric orbit, present osculating orbit, future barycentric orbit, and next perihelion or escape. More recent work also uses the term for process-resolved or stage-wise schemata that encode volatile drivers, reservoir depths, and seasonal or orbital evolution (Ferrín, 2013, Królikowska et al., 2020, Parhi et al., 30 Oct 2025).

1. Terminological scope and principal variants

In the cited literature, the term “Comet Evolutionary Diagram” refers to several technically distinct representations rather than to one universally fixed diagram. The main variants can be summarized as follows.

Variant	Coordinates or nodes	Primary purpose
Empirical RR–ML-Age / RR–WB-AGE diagram	Horizontal: ML-Age or WB-AGE; vertical: RR	Comparative evolutionary classification
CODE five-stage dynamical diagram	Previous → original → osculating → future → next	Orbital-energy and perihelion-distance evolution
Model-derived activity/reservoir schema	$r$, $Q_i(r)$, depth, epoch, dominant driver	Thermal-compositional interpretation
Stage-wise mission synthesis	Orbital stage A–D with materials and observables	Seasonal and morphological evolution

The empirical RR–ML-Age form is described as a two-parameter diagnostic that places objects according to how much volatile mass they lose per apparition and how many returns they have left before either sublimating away or being choked by a dust mantle. In that framework, the CED is explicitly compared to the Hertzsprung–Russell diagram in the sense that it organizes comets into physically meaningful families and regimes (Ferrin et al., 10 Apr 2026, Ferrín, 2013).

The CODE form is different in both ontology and data product. It is a compact visualization of a long-period comet’s five standardized dynamical snapshots: a previous perihelion, the inbound original barycentric orbit at 250 au, the observed osculating orbit near the current perihelion, the outbound future barycentric orbit at 250 au, and the next perihelion or escape. Its emphasis is not comparative placement on a population plane, but the evolution of orbital energy and perihelion distance under non-gravitational outgassing, planetary scattering, and Galactic and stellar perturbations (Królikowska et al., 2020).

A third usage appears in the coupled orbital–thermal model of Oort Cloud comets. That paper does not present a named CED, but it explicitly states that its outputs provide all the ingredients to construct one, organized by activity drivers versus heliocentric distance and epoch, volatile retention and depletion by depth and size, and production rates with detectability thresholds and active-area adjustments (Parhi et al., 30 Oct 2025). A related stage-wise synthesis for 67P/Churyumov–Gerasimenko organizes composition, morphology, frost, erosion, and seasonal forcing into a structured CED anchored by Rosetta results (Filacchione et al., 2022).

2. The empirical RR–age framework

The empirical CED is built on two scalar quantities: an age proxy derived from mass loss, and the number of remaining returns implied by the current erosion rate. In the mass-loss formulation, the horizontal axis is Mass-Loss Age, ML-Age, defined as

$$\mathrm{ML\text{-}Age} = \frac{k}{M_{\mathrm{loss,tot}}},$$

with $k = 3.58 \times 10^{11}\ \mathrm{kg\cdot cy}$ and $M_{\mathrm{loss,tot}}$ the integrated mass lost over the active arc. In the water-budget formulation, the horizontal axis is Water-Budget Age,

$$\mathrm{WB\text{-}AGE}\,[\mathrm{cy}] = \frac{3.58 \times 10^{11}}{\mathrm{WB}\,[\mathrm{kg}]},$$

with the normalization chosen so that comet 28P/Neujmin 1 corresponds to $100$ comet years. Larger ML-Age or WB-AGE means smaller mass loss per orbit and therefore a more evolved object in this comparative sense (Ferrin et al., 10 Apr 2026, Ferrin et al., 2013, Ferrín, 2013).

The vertical axis is Remaining Returns, RR. Two equivalent forms are used:

$\mathrm{RR} = \frac{r}{\Delta r},$

and

$\mathrm{RR} = \frac{3\,M_\mathrm{nuc}}{\Delta M}.$

Here $r$ is nucleus radius, $Q_i(r)$0 is the shell thickness removed in one active passage, $Q_i(r)$1 is nucleus mass, and $Q_i(r)$2 is the total mass loss per return. The shell-thickness relation is

$Q_i(r)$3

or, in the water-budget formulation with dust-to-gas ratio $Q_i(r)$4,

$Q_i(r)$5

RR therefore estimates how many comparable future returns remain before the nucleus either reaches the disintegration limit, $Q_i(r)$6, or evolves toward a suffocated state in which $Q_i(r)$7 and $Q_i(r)$8 (Ferrin et al., 10 Apr 2026, Ferrin et al., 2013, Ferrín, 2013).

In this framework, the CED is plotted on logarithmic axes spanning many orders of magnitude. The left side contains young Oort-cloud comets with large mass budgets and small ML-Age. The middle contains Jupiter-family comets. The upper right contains old, dust-mantled, or “graveyard” objects, including asteroidal-belt comets and active asteroids. The lower boundary at $Q_i(r)$9 marks disintegrators such as C/2002 O4 Hönig and C/2012 S1 ISON. One paper further defines a “graveyard of comets” by $$\mathrm{ML\text{-}Age} = \frac{k}{M_{\mathrm{loss,tot}}},$$0, and another places main-belt or asteroidal-belt comets in the upper-right graveyard at large age and large RR (Ferrin et al., 2013, Ferrín, 2013).

The empirical CED is closely tied to Secular Light Curve methodology. Photometry is reduced to absolute magnitude, often using

$$\mathrm{ML\text{-}Age} = \frac{k}{M_{\mathrm{loss,tot}}},$$1

with the important caveat that for optically thick comae the nucleus is not seen and $$\mathrm{ML\text{-}Age} = \frac{k}{M_{\mathrm{loss,tot}}},$$2 is effectively undefined, so an envelope method is used rather than a phase correction to the coma-dominated flux. The envelope of the data is taken as the “truth,” because practical observing limitations bias magnitudes downward. Production-rate integration over the active interval then yields the mass budget. For dust, one widely used empirical conversion is Af$$\mathrm{ML\text{-}Age} = \frac{k}{M_{\mathrm{loss,tot}}},$$3 cm corresponding to $$\mathrm{ML\text{-}Age} = \frac{k}{M_{\mathrm{loss,tot}}},$$4, while the paper on ISON writes the equivalent scaling as $$\mathrm{ML\text{-}Age} = \frac{k}{M_{\mathrm{loss,tot}}},$$5 (Ferrin et al., 10 Apr 2026, Ferrín, 2013).

Several interpretive borders are integral to this version of the CED. A heuristic separation is drawn between sublimating-away comets, which lose ice with little residue and evolve downward toward $$\mathrm{ML\text{-}Age} = \frac{k}{M_{\mathrm{loss,tot}}},$$6, and suffocating comets, where dust fallback quenches activity and drives migration up and to the right. The “suffocation–sublimation border” is estimated in one study at $$\mathrm{ML\text{-}Age} = \frac{k}{M_{\mathrm{loss,tot}}},$$7, and a lower-right “desert” is predicted because very old, nearly exhausted comets should become undetectably small (Ferrín, 2013, Ferrin et al., 10 Apr 2026).

3. The CODE dynamical CED for long-period comets

The CODE catalogue defines a different CED centered on orbital dynamics. It records each long-period comet in five stages that cover three successive perihelion passages and form the five nodes of a polyline. Snapshot 1 is the observed osculating orbit near the current perihelion and is heliocentric. Snapshot 2 is the original barycentric orbit at 250 au inbound. Snapshot 3 is the future barycentric orbit at 250 au outbound. Snapshot 4 is the previous perihelion, obtained by long-term backward integration under external perturbations. Snapshot 5 is the next perihelion or the escape state at 120,000 au (Królikowska et al., 2020).

The original and future barycentric states isolate the net effect of planetary perturbations during the observed perihelion passage. The previous and next perihelia connect that short-timescale planetary scattering to long-timescale Galactic and stellar forcing. In practice, the nominal orbit is cloned into a swarm of 5001 virtual comets using a Monte Carlo method, and each swarm is integrated to 250 au for the original and future states. For the previous and next states, the integrations include the Galactic potential and a curated list of 643 stars or stellar systems that may approach within 4 pc over approximately $$\mathrm{ML\text{-}Age} = \frac{k}{M_{\mathrm{loss,tot}}},$$8 Myr. Integrations stop at the previous or next perihelion for returning virtual comets, or at an escape limit of 120,000 au for escaping ones (Królikowska et al., 2020).

The most common vertical axis in this CED is inverse semimajor axis, $$\mathrm{ML\text{-}Age} = \frac{k}{M_{\mathrm{loss,tot}}},$$9, plotted in units of $k = 3.58 \times 10^{11}\ \mathrm{kg\cdot cy}$0, because it is a proxy for specific orbital energy:

$k = 3.58 \times 10^{11}\ \mathrm{kg\cdot cy}$1

Higher $k = 3.58 \times 10^{11}\ \mathrm{kg\cdot cy}$2 corresponds to more tightly bound orbits; for hyperbolas, $k = 3.58 \times 10^{11}\ \mathrm{kg\cdot cy}$3 and $k = 3.58 \times 10^{11}\ \mathrm{kg\cdot cy}$4. A complementary axis is perihelion distance,

$k = 3.58 \times 10^{11}\ \mathrm{kg\cdot cy}$5

used either as the horizontal axis or as an annotation at the five nodes to show how energy and perihelion evolve together (Królikowska et al., 2020).

This formulation explicitly incorporates non-gravitational acceleration. CODE uses the Marsden–Sekanina–Yeomans formalism,

$k = 3.58 \times 10^{11}\ \mathrm{kg\cdot cy}$6

with

$k = 3.58 \times 10^{11}\ \mathrm{kg\cdot cy}$7

and

$k = 3.58 \times 10^{11}\ \mathrm{kg\cdot cy}$8

normalized so that $k = 3.58 \times 10^{11}\ \mathrm{kg\cdot cy}$9. CODE lists standard parameter sets for water-ice sublimation, typically for $M_{\mathrm{loss,tot}}$0 au, and for CO-driven activity, often for $M_{\mathrm{loss,tot}}$1 au. Asymmetric variants are used when warranted by residuals, but symmetric models are preferred unless an asymmetric model significantly improves the fit (Królikowska et al., 2020).

The dynamical classifications attached to this CED also differ from those in the RR–age literature. CODE’s sample centers on Oort-spike comets defined by purely gravitational original semimajor axis $M_{\mathrm{loss,tot}}$2 au, or equivalently

$M_{\mathrm{loss,tot}}$3

However, dynamical new versus old is classified primarily by the previous perihelion distance obtained from the virtual-comet swarm: $M_{\mathrm{loss,tot}}$4 au implies dynamically new, $M_{\mathrm{loss,tot}}$5 au dynamically old, and $M_{\mathrm{loss,tot}}$6 au uncertain. For hyperbolic cases, the five-stage scheme is adapted so that the next snapshot is taken at the 120,000 au escape limit, often with $M_{\mathrm{loss,tot}}$7 (Królikowska et al., 2020).

Uncertainty treatment is intrinsic to the CODE CED. At snapshots 2 and 3, distributions at 250 au are often approximately Gaussian, so mean and standard deviation are reported. At snapshots 4 and 5, non-Gaussian outcomes are common, so CODE provides counts of returning, escaping, and hyperbolic virtual comets and summarizes distributions by deciles. On the diagram, this permits error bars at the original and future stages and percentile boxes or ranges at the previous and next stages (Królikowska et al., 2020).

4. Thermophysical and mission-based extensions

A model-driven extension of the CED concept appears in the coupled orbital–thermal evolution of Oort Cloud comets. There the orbit $M_{\mathrm{loss,tot}}$8 is fed continuously into a thermophysical solver over $M_{\mathrm{loss,tot}}$9 Gyr, with four dynamical phases: early planetary evolution near Neptune with $$\mathrm{WB\text{-}AGE}\,[\mathrm{cy}] = \frac{3.58 \times 10^{11}}{\mathrm{WB}\,[\mathrm{kg}]},$$0 au for about $$\mathrm{WB\text{-}AGE}\,[\mathrm{cy}] = \frac{3.58 \times 10^{11}}{\mathrm{WB}\,[\mathrm{kg}]},$$1 Myr, outward migration over about $$\mathrm{WB\text{-}AGE}\,[\mathrm{cy}] = \frac{3.58 \times 10^{11}}{\mathrm{WB}\,[\mathrm{kg}]},$$2 Myr, Oort Cloud residence with $$\mathrm{WB\text{-}AGE}\,[\mathrm{cy}] = \frac{3.58 \times 10^{11}}{\mathrm{WB}\,[\mathrm{kg}]},$$3 au for about $$\mathrm{WB\text{-}AGE}\,[\mathrm{cy}] = \frac{3.58 \times 10^{11}}{\mathrm{WB}\,[\mathrm{kg}]},$$4 Gyr, and a final return inward with $$\mathrm{WB\text{-}AGE}\,[\mathrm{cy}] = \frac{3.58 \times 10^{11}}{\mathrm{WB}\,[\mathrm{kg}]},$$5 au. The nuclei are porous, spherical dust–ice mixtures with radii of 2, 10, and 50 km, and the initial volatile inventory includes CO and CO$$\mathrm{WB\text{-}AGE}\,[\mathrm{cy}] = \frac{3.58 \times 10^{11}}{\mathrm{WB}\,[\mathrm{kg}]},$$6 ices, amorphous water ice with trapped CO and CO$$\mathrm{WB\text{-}AGE}\,[\mathrm{cy}] = \frac{3.58 \times 10^{11}}{\mathrm{WB}\,[\mathrm{kg}]},$$7, and dust (Parhi et al., 30 Oct 2025).

The governing physics includes porous heat conduction, sublimation and refreezing within pores, Knudsen-regime gas transport, crystallization of amorphous ice, and a test of long-lived radiogenic heating. The surface energy balance is written as

$$\mathrm{WB\text{-}AGE}\,[\mathrm{cy}] = \frac{3.58 \times 10^{11}}{\mathrm{WB}\,[\mathrm{kg}]},$$8

and the crystallization kinetics are

$$\mathrm{WB\text{-}AGE}\,[\mathrm{cy}] = \frac{3.58 \times 10^{11}}{\mathrm{WB}\,[\mathrm{kg}]},$$9

The resulting CED content is not a single point but a structured map of dominant activity drivers, volatile reservoir accessibility, production rates, and detectability thresholds (Parhi et al., 30 Oct 2025).

Several quantitative thresholds emerge from that model. CO sublimation is the dominant large-distance driver and first exceeds detectability at $100$0 au for $100$1 km and $100$2 au for $100$3 km. CO$100$4 becomes detectable around $100$5 au in the global-average implementation, and under subsolar or limited-active-area assumptions can be detectable out to about $100$6–$100$7 au. Amorphous ice begins crystallizing around $100$8 au, and the release of trapped gases becomes the main source of outgassing around $100$9–$\mathrm{RR} = \frac{r}{\Delta r},$0 au, with peak total production reaching about $\mathrm{RR} = \frac{r}{\Delta r},$1 for the 10 km model and about $\mathrm{RR} = \frac{r}{\Delta r},$2 for the 50 km model. The model also finds that CO ice is completely depleted in the 2 km nucleus but preserved from depths of about 500 m inward in the larger nuclei, while CO$\mathrm{RR} = \frac{r}{\Delta r},$3 and amorphous H$\mathrm{RR} = \frac{r}{\Delta r},$4O are entirely preserved; radiogenic heating by long-lived isotopes is negligible (Parhi et al., 30 Oct 2025).

A mission-anchored extension appears in the structured CED for 67P/Churyumov–Gerasimenko. That synthesis uses Rosetta observations to define four stages: aphelion dormancy, inbound activation, perihelion peak activity, and outbound relaxation. The nucleus is described as very dark, red-sloped, and covered by a variable-thickness dust mantle, with geometric albedo $\mathrm{RR} = \frac{r}{\Delta r},$5 at 535 nm and $\mathrm{RR} = \frac{r}{\Delta r},$6 at 550 nm, porosity around $\mathrm{RR} = \frac{r}{\Delta r},$7–$\mathrm{RR} = \frac{r}{\Delta r},$8, patchy H$\mathrm{RR} = \frac{r}{\Delta r},$9O and transient CO$\mathrm{RR} = \frac{3\,M_\mathrm{nuc}}{\Delta M}.$0 exposures, ammoniated salts, refractory organics with aliphatic and aromatic components, and hydroxylated Mg-amorphous silicates contributing to the broad 3.2 $\mathrm{RR} = \frac{3\,M_\mathrm{nuc}}{\Delta M}.$1m absorption (Filacchione et al., 2022).

The stage-wise 67P CED is explicitly tied to heliocentric forcing and energy balance. Frost lines are given as H$\mathrm{RR} = \frac{3\,M_\mathrm{nuc}}{\Delta M}.$2O at about 3 AU, CO$\mathrm{RR} = \frac{3\,M_\mathrm{nuc}}{\Delta M}.$3 at about 14 AU, CO at about 28 AU, and N$\mathrm{RR} = \frac{3\,M_\mathrm{nuc}}{\Delta M}.$4 at about 32 AU. The instantaneous energy balance is written as

$\mathrm{RR} = \frac{3\,M_\mathrm{nuc}}{\Delta M}.$5

with sublimation flux represented by the Hertz–Knudsen relation. Observationally, the structured diagram encodes dawn frost cycles in Hapi, CO$\mathrm{RR} = \frac{3\,M_\mathrm{nuc}}{\Delta M}.$6 ice in Anhur, perihelion outbursts from cliffs and pits, net erosion of about 4 m near perihelion, and post-perihelion fallback that blankets northern terrains and restores redder colors. This suggests a CED usage in which composition, morphology, seasonal illumination asymmetry, and observable spectral or photometric changes are treated as a single evolutionary system rather than as scalar summary indices (Filacchione et al., 2022).

5. Representative placements and case studies

The empirical RR–age literature uses the CED to place specific comets within a comparative evolutionary taxonomy. For exocomet 3I/ATLAS, the Secular Light Curve and integrated production rates yield

$\mathrm{RR} = \frac{3\,M_\mathrm{nuc}}{\Delta M}.$7

The total mass loss is reported as $\mathrm{RR} = \frac{3\,M_\mathrm{nuc}}{\Delta M}.$8 kg, composed of dust $\mathrm{RR} = \frac{3\,M_\mathrm{nuc}}{\Delta M}.$9, H$r$0O $r$1, CO$r$2 $r$3, and CO $r$4, making the object strongly CO$r$5-dominated. On the log–log CED it lies in the extreme lower-left domain characteristic of young, icier comets with minimal dust mantling, colocated with the Oort-cloud group. The same work reports a photometric anomaly from $r$6 to $r$7 days before perihelion, interpreted as an eclipse consistent with a possible binary nucleus, and a maximum reduced absolute magnitude of $r$8 about 15 days after perihelion (Ferrin et al., 10 Apr 2026).

The same empirical framework places C/2012 S1 ISON at the disintegration limit. Its Secular Light Curve showed a Slope Discontinuity Event at about $r$9 AU, followed by a near-standstill lasting more than 132 days. The paper reports $Q_i(r)$00 kg, dust budget $Q_i(r)$01 kg, CO budget $Q_i(r)$02 kg, and total $Q_i(r)$03 kg, implying

$Q_i(r)$04

With $Q_i(r)$05, the corresponding diameter is $Q_i(r)$06 m. The paper treated the SDE+standstill signature as a predictor of disintegration and reports that the comet indeed fragmented at $Q_i(r)$07 AU. By contrast, C/2011 L4 Panstarrs, with an SDE at $Q_i(r)$08 AU and post-SDE brightening to $Q_i(r)$09, is placed in the Oort-cloud region on the left side of the CED (Ferrín, 2013).

Asteroidal-belt or asteroidal-belt-comet placements occupy the opposite corner. For 133P/Elst–Pizarro, the paper gives

$Q_i(r)$10

with $Q_i(r)$11 for $Q_i(r)$12 and $Q_i(r)$13 for $Q_i(r)$14. For 107P/Wilson–Harrington, $Q_i(r)$15 cy and $Q_i(r)$16 to $Q_i(r)$17. For 300163 = 2006 VW139, $Q_i(r)$18. These values place the objects firmly in the upper-right graveyard, where activity is weak, mantled, and episodically rejuvenated by decreases in perihelion distance. The same paper gives dust-mantle thicknesses of $Q_i(r)$19 m for 107P, $Q_i(r)$20 m for 133P, and $Q_i(r)$21 m for a sample of nine comets, derived from post-perihelion activation lags and a thermal-wave model (Ferrin et al., 2013).

The CODE dynamical CED yields a different type of case study. For C/2002 T7 LINEAR, the preferred NG solution “d6” gives observed osculating perihelion $Q_i(r)$22 au and

$Q_i(r)$23

at 250 au inbound. Across water-ice NG variants, $Q_i(r)$24 ranges from $Q_i(r)$25 to $Q_i(r)$26; a CO-driven model yields $Q_i(r)$27 but is judged implausible for this small-$Q_i(r)$28 comet. The previous perihelion is at approximately 150 au, making the comet dynamically new by the $Q_i(r)$29 criterion, while the future orbit is hyperbolic, with the next snapshot recorded at the 120,000 au escape limit. In the CODE CED, the track therefore moves from small positive $Q_i(r)$30 at the original stage to negative $Q_i(r)$31 at the future and next stages, visually encoding ejection by strong planetary scattering (Królikowska et al., 2020).

6. Interpretation, misconceptions, and unresolved issues

A recurrent misconception is to treat the CED as a single standard diagram. The literature does not support that usage. At minimum, the term denotes a comparative RR–age plane, a five-snapshot dynamical trajectory for long-period comets, and model-derived or stage-wise schemata that encode activity drivers and material evolution. These variants are related by their concern with comet evolution, but their variables, timescales, and inferential targets differ (Królikowska et al., 2020, Ferrin et al., 10 Apr 2026, Parhi et al., 30 Oct 2025).

A second misconception is to read ML-Age or WB-AGE as chronological age. In the empirical CED, these are explicitly proxies based on inverse mass loss, normalized in comet years by a calibration constant. They express comparative evolutionary state under current or inferred activity, not time since formation. Likewise, RR is not a deterministic lifetime; it is an estimate under the assumption that the per-apparition mass loss remains roughly constant (Ferrin et al., 10 Apr 2026, Ferrin et al., 2013).

In the CODE framework, dynamical new versus old is also frequently misunderstood. CODE does not classify long-period comets primarily by a hard threshold in $Q_i(r)$32. The operative criterion is the previous perihelion distance from the virtual-comet swarm, with $Q_i(r)$33 au marking dynamically new, $Q_i(r)$34 au dynamically old, and the intermediate range treated as uncertain. The Tisserand parameter with respect to Jupiter can aid interpretation of planetary-scattering signatures during the osculating stage,

$Q_i(r)$35

but CODE does not classify long-period comets by $Q_i(r)$36, and for LPCs it is best used only as a qualitative indicator because it is not conserved across deep planetary encounters (Królikowska et al., 2020).

The limitations of the various CEDs are explicit in the source literature. In RR–age studies, photometric apertures can under-extract coma flux, water datasets may require scale corrections by factors of about 5 to 12, and the Af$Q_i(r)$37 conversion is empirical and order-of-magnitude in nature. Derived quantities depend on adopted density, albedo, size, and production-rate integration, although the authors emphasize that the log–log diagram is “forgiving” because factor-of-two uncertainties move points little (Ferrin et al., 10 Apr 2026, Ferrín, 2013).

In CODE, previous and next perihelion outcomes depend on the completeness and accuracy of the stellar perturber list, on the adopted non-gravitational model, and on whether asymmetric $Q_i(r)$38 or data-subset solutions are needed to avoid biased osculating fits. Uncertainties become non-Gaussian at long timescales, which is why deciles rather than means and standard deviations are often reported for the previous and next stages (Królikowska et al., 2020).

The thermophysical extensions add another layer of uncertainty. The coupled orbital–thermal model is one-dimensional and spherically symmetric, anchored to a single representative dynamical history; porosity, pore size, tortuosity, thermal conductivity, dust fraction, and initial volatile abundances all shift depletion depths, refreezing profiles, onset distances, and predicted production rates. The 67P synthesis, although exceptionally rich observationally, remains a stage-wise abstraction of a nucleus whose activity is strongly controlled by local shading, morphology, and seasonal asymmetry (Parhi et al., 30 Oct 2025, Filacchione et al., 2022).

Taken together, the literature presents the CED not as one diagram but as a compact language for encoding comet evolution at different levels: comparative population state, dynamical transfer through the planetary system and the Oort cloud, and thermophysical or seasonal reorganization of surface and subsurface volatiles. A plausible implication is that future convergence among these variants will require models that connect mass-loss-based age proxies, orbit-resolved energy changes, and internally resolved volatile evolution within a single formalism.