A Primordial Origin for the Gas-rich Debris Disks around Intermediate-mass Stars

While most debris disks consist of dust with little or no gas, a fraction have significant amounts of gas detected via emission lines of CO, ionized carbon, and/or atomic oxygen. Almost all such gaseous debris disks known are around A-type stars with ages up to 50 Myr. We show, using semianalytic disk evolution modeling, that this can be understood if the gaseous debris disks are remnant protoplanetary disks that have become depleted of small grains compared to the interstellar medium. Photoelectric heating by the A stars’ far-UV (FUV) radiation is then inefficient, while the stars’ extreme-UV (EUV) and X-ray emissions are weak owing to a lack of surface convective zones capable of driving magnetic activity. In this picture, it is relatively difficult for stars outside the range of spectral types from A through early B to have such long-lived gas disks. Less-massive stars have stronger magnetic activity in the chromosphere, transition region, and corona with resulting EUV and X-ray emission, while more-massive stars have photospheres hot enough to produce strong EUV radiation. In both cases, primordial disk gas is likely to photoevaporate well before 50 Myr. These results come from 0D disk evolution models where we incorporate internal accretion stresses, MHD winds, and photoevaporation by EUV and X-ray photons with luminosities that are functions of the stellar mass and age. A key issue this work leaves open is how some disks become depleted in small dust so that FUV photoevaporation slows. Candidates include the grains’ growth, settling, radial drift, radiation force, and incorporation into planetary systems.


INTRODUCTION
A fundamental goal of planet formation research is to learn the evolutionary pathways that protoplanetary disks (PPDs) follow to become planetary systems and their associated debris disks.Classically the PPD lifetime is estimated to be no longer than 10 Myr, based on near-and mid-infrared (IR) observations which trace the inner hot dust (≲ 1-10 au; Haisch et al. 2001;Meyer et al. 2007;Hernández et al. 2007;Mamajek 2009;Ribas et al. 2014Ribas et al. , 2015) ) and ultraviolet (UV) and Hα observations tracing accreting gas (e.g., Kennedy & Kenyon 2009;Fedele et al. 2010;Sicilia-Aguilar et al. 2010).Ten Myr is thus considered a typical timescale on which gasrich PPDs evolve into gas-free debris disks.
An explanation is therefore needed for the more than 20 gaseous debris disks observed around stars older than 10 Myr (e.g., Kóspál et al. 2013;Marino et al. 2016;Lieman-Sifry et al. 2016;Moór et al. 2017;Hughes et al. 2017;Higuchi ryohei.nakatani@jpl.nasa.gov et al. 2017Higuchi ryohei.nakatani@jpl.nasa.gov et al. , 2019b,a;,a;Kral et al. 2020).Two hypotheses have been proposed for the origins of this gas.In the primordialorigin scenario, the gas is a remnant of the PPD.In the secondary-origin scenario, the gas was released more recently from volatiles in planetesimals associated with the dust in the debris disks.
While the secondary-origin scenario has been extensively examined through modeling (Kral et al. 2016(Kral et al. , 2019;;Moór et al. 2019;Matrà et al. 2019;Marino et al. 2020;Cataldi et al. 2020;Marino et al. 2022), less attention has been paid to the primordial-origin scenario (Nakatani et al. 2021;Smirnov-Pinchukov et al. 2022;Iwasaki et al. 2023).A primordial origin has been considered unlikely because of the short canonical PPD gas disk lifetime.However, hydrodynamical modeling by Nakatani et al. (2021) shows that if the PPD around an intermediate-mass star (M * ≈ 2 M ⊙ ) is depleted in submicron grains, including polycyclic aromatic hydrocarbons (PAHs), then the photoelectric heating produced when FUV photons strike grains is weak enough to slow photoevaporation so the gas could survive well past 10 Myr and into the debris disk stage -without small-grain depletion, the gas lifetimes are most likely ≲ 10 Myr regardless of stellar masses (Komaki et al. 2023).Primordial origins thus have the potential to account for the occurrence of gaseous debris disks around A stars (Lieman-Sifry et al. 2016;Moór et al. 2017;Hughes et al. 2018).
The outstanding questions include what gas lifetimes to expect quantitatively in the primordial origin scenario, how these depend on the stellar mass, and whether the results are consistent with the observed occurrence rates of gaseous debris disks.We here address these issues by constructing disk evolution models using a semi-analytic, 0D approach.
This picture of time-varying disk dispersal processes is consistent with several observational facts.For instance, class 0/I jets have much larger mass-loss rates (10 −7 -10 −6 M ⊙ yr −1 ) than photoevaporation can provide, implying the ejection is driven by MHD effects in these early stages.For another example, the luminosity of a photoevaporative wind tracer, the [Ne II] 12.8 µm low-velocity component, is higher for more evolved disks, i.e., slow accretors (∼ 10 −8 M ⊙ yr −1 ) and inner-dust-depleted disks (Pascucci et al. 2020).
Jet and wind mass-loss rates in the upper end of the range imply possibly enough material crossing the lines of sight between star and disk to block much of the stellar ultraviolet (UV) and X-ray radiation.This can slow or stop the disk's photoevaporation (for a recent review, see Pascucci et al. 2022).In 3D MHD modeling by Takasao et al. (2022), the outflow from the disk around a solar-mass star has a column density high enough to attenuate the UV and X-ray photons at accretion rates ≳ 10 −8 M ⊙ yr −1 .Overall, a fairly advanced age and low rates of accretion and outflow appear necessary before photoevaporation can dominate gas disk dispersal.

GAS LIFETIME ESTIMATE
For the primordial-origin hypothesis to be a possible scenario, one of the critical necessary conditions PPDs must meet is to retain a mass reservoir beyond 10 Myr, especially at outer radii (≳ 1-10 au) where gaseous debris disks are detected.This study's primary objective is to demonstrate its feasibility, which was previously considered implausible.Our model will be tailored for this purpose, as will described below.
We focus on the disks where PAH and small grains (≲0.01 µm) are depleted compared to the interstellar medium.These are the systems that potentially meet the necessary condition; otherwise, the disks dissipating through rapid photoelectric heating most likely have short lifetimes < 10 Myr (Gorti et al. 2009(Gorti et al. , 2015;;Komaki et al. 2021Komaki et al. , 2023)); we ignore this population.Thus, we here derive the lifetime for the sub-population of PPDs where the gas is the longestlived.The estimated lifetime is distinct from the average gas lifetime, which depends on the still-unknown probability of achieving a small-grain-depleted state.
We do not specify here the processes responsible for grain depletion.Instead, we simply assume the disks reach a smallgrain-depleted state before photoevaporation begins to dominate disk dispersal, regardless of the stellar mass.Several effects can deplete the dust, especially for intermediate-mass stars: Grain growth and radial drift reduce the local dustto-gas ratio (e.g., Sellek et al. 2020).The drift is likely especially rapid around young intermediate-mass stars (Pinilla et al. 2022).The intense radiation force from intermediatemass stars can also blow out small dust (Owen & Kollmeier 2019;Nakatani et al. 2021).
The specific level of small-grain depletion needed to inhibit FUV photoevaporation remains uncertain.Nakatani et al. (2018a,b) explored FUV photoevaporation varying disk metallicity and demonstrated that strong FUV photoevaporation demands small grain abundance exceeding ∼ 0.1-1% of the ISM level.However, these investigations focus on a 0.5 M ⊙ star with extremely high FUV luminosity, lowering oxygen and carbon abundances with metallicity, which limits the results' general applicability.The critical smallgrain abundance depends on stellar mass and FUV/X-ray luminosities, yet prior studies explored a limited parameter space.Further investigations are warranted.Nevertheless, it is likely that critical small-grain abundances exist, below which FUV photoevaporation becomes negligible.We anticipate this threshold to be around ∼ 0.1-1% of the ISM level, following Nakatani et al. (2018a,b).
Observations suggest the PAH abundance in PPDs is at least ten times less than the interstellar value (e.g., Geers et al. 2007;Oliveira et al. 2010;Vicente et al. 2013).The hybrid disk HD 141569A is depleted to a level where photoelectric heating may drive only weak photoevaporation (Li & Lunine 2003;Thi et al. 2014).Non-detection of PAHs occur with a reasonable fraction in Herbig disks, where the detection rate is ∼ 70% (e.g., Acke et al. 2010), or in T Tauri disks, where it is ≲ 10%; (e.g., Geers et al. 2006).Such non-detection hints at abundances potentially over 100 times lower than interstellar levels.This significant depletion in the disk atmosphere is also indicated from mid-IR spectra (e.g, Furlan et al. 2006Furlan et al. , 2011) ) and implied from IR scattered-light imaging (e.g., Mulders et al. 2013).These findings lead us to infer that these disks might be in the small-grain-depleted state.

Semi-Analytic Model
Our model employs a 0D approach, following the disk's mass decrease from accretion and wind mass loss in the early stage and from photoevaporation in the late stage.It is based on the analytic lifetime estimate from the UV-switch model (Clarke et al. 2001, detailed in Appendix A).In their model, the major dispersal process is viscous accretion at the early stage and switches to photoevaporation beginning at the "switching time" t 0 , which defines the transition point into the late stage.Thereafter, it disperses the remaining disk mass of M ′ disk,0 .This model returns t 0 , M ′ disk,0 , and the gas disk lifetime t life,g analytically (cf.Eqs.A5-A7) once the disk mass M disk (t) (cf.Eq.( A4)) and photoevaporation rate Ṁph (t) are given as functions of time.
The original UV-switch model uses purely viscous disks, but here we consider disks evolving under angular momentum and mass extraction through a magnetically-driven wind alongside the internal angular momentum redistribution by turbulent stresses.Observational evidence that PPDs appear to evolve by driving magnetized winds is discussed by Rafikov (2017), Manara et al. (2022), andPascucci et al. (2022).We adopt the analytic disk evolution model of Chambers (2019) for M disk (t), as set out in §B.We first derive t 0 and M ′ disk,0 by finding the time when the disk mass loss Ṁdisk (t) and photoevaporation rates Ṁph (t) equal; we will define Ṁph (t) later.t life,g is then the time till the remaining mass is photoevaporated, found by numerically integrating The resulting t life,g is insensitive to the choice of the threshold mass to define "dispersal", as M disk falls rapidly after t = t 0 .The disk's initial mass is taken as proportional to the stellar mass with M disk,0 = 0.1M * (e.g., Williams & Cieza 2011;Andrews et al. 2013, however a steeper scaling ∝ M 1.3-2.0* is found by Pascucci et al. (2016); Ansdell et al. (2016Ansdell et al. ( , 2017)); see Manara et al. (2022) for a recent review).We explore t life,g 's dependence on M disk,0 and the scaling in Appendix D.
Photoevaporation is generally driven by far-ultraviolet (FUV), extreme-ultraviolet (EUV), and X-ray radiation.However, as stated above, we here consider only cases where FUV grain photoelectric heating is negligible, and the photoevaporative winds are driven solely by the EUV and X-ray heating associated with photoionization.We also assume that the disks are isolated from massive stars, so external photoevaporation is negligible.We use the EUV photoevaporation rate estimate (1) (Hollenbach et al. 1994;Clarke et al. 2001) since it agrees well with the mass-loss rates derived by radiation hydrodynamics simulations (Nakatani et al. 2021).Based on 1+1-D modeling, EUV photoevaporation rates were thought to increase about tenfold after the disk center is cleared out, letting the stellar EUV field dominate over the diffuse EUV from recombining ions.However, 2D axisymmetric radiative transfer showed that the direct stellar component dominates even in disks without central cavities (Tanaka et al. 2013).
Additionally, radiation hydrodynamics modeling found no significant difference in photoevaporation rates with and without cavities of various sizes (Owen et al. 2010;Picogna et al. 2019;Nakatani et al. 2021).Hence, we apply the Eq.( 1) photoevaporation rates uniformly throughout disk evolution.
The stellar EUV emission rate Φ EUV can be decomposed into photospheric, chromospheric (∼ 10 4 K), transition regional (∼ 10 5 K), and coronal (∼ 10 6 K) components.We shall refer to the last three jointly as magnetic-origin EUV, Φ EUV,mag .For the photospheric EUV component, we use Kunitomo et al. (2021) Table 1 results from long-term stellar evolution calculations for 0.5 M ⊙ ≤ M * ≤ 5 M ⊙ .We estimate the magnetic-origin EUV emission rate from the X-ray luminosity L X using the Φ EUV -L X relation of Shoda & Takasao (2021), Strictly, this relation applies to solar-type stars, but we apply it to all the stars considered here (0.5 M ⊙ ≤ M * ≤ 5 M ⊙ ).Since the theoretical and observational uncertainties in Φ EUV,mag are large, we believe this estimate is a reasonable first-order approximation.
We follow Kunitomo et al. (2021) for L X time evolution, excluding the L X /L * = 10 −7 floor applied to intermediatemass stars in the original paper.We expect the floor not to apply universally, given that most Herbig stars do not present X-ray emissions detected at this level (e.g., Hamaguchi et al. 2005;Stelzer et al. 2006).Instead, we employ L X = min(10 −3.13 , 5.3 × 10 −6 Ro −2.7 )L * , where Ro represents the Rossby number (Mangeney & Praderie 1984;Noyes et al. 1984) defined as Ro = P rot /τ conv with P rot and τ conv being the rotational period and the convective turnover timescale.Kunitomo et al. (2021) set P rot to 3 days and calculate τ conv from the stellar evolution model.The rotational period can range between 1-10 days in general, but the qualitative trend in the time evolution of L X is independent of P rot (see Section 5.2 of Kunitomo et al. (2021)).Thus, the variation of P rot would have a limited impact on our results.The total EUV emission rate is Φ EUV = Φ EUV,ph + Φ EUV,mag .We omit accretiongenerated EUV and X-ray here but address their influence on t life,g in Appendix C.
Stellar/interstellar Lyman-Werner (LW) radiation-driven photoevaporation might be important for small-graindepleted disks around stars with 2 M ⊙ ≲ M * ≲ 3 M ⊙ (hereafter late intermediate-mass stars).However, our radiation hydrodynamics simulations indicate that LW photoevaporation is negligible even for late intermediate-mass stars (Nakatani et al. in prep.).Thus, we do not consider LW photoevaporation in this study.The dependence of LW photoevaporation on stellar mass and luminosity is currently unknown but warrants investigation in future research.
For X-ray photoevaporation, the mass-loss rate is still under debate.Some studies obtained large mass-loss rates (Ercolano et al. 2009;Owen et al. 2010;Picogna et al. 2019), while others found weaker mass loss (Gorti et al. 2009;Wang & Goodman 2017;Nakatani et al. 2018b;Komaki et al. 2021).The different conclusions likely originate from the adopted X-ray spectra -mass-loss rates are high if the spectrum has a certain level of soft X-rays (∼ 0.1 keV, Ercolano et al. 2009;Gorti et al. 2009;Nakatani et al. 2018b;Sellek et al. 2022), from the incorporated cooling processes -the studies with large X-ray photoevaporation rates neglect important cooling processes, and from the adopted numerical methods.To cover the uncertainty of X-ray photoevaporation rates, we estimate the lifetimes with and without X-ray photoevaporation.
For X-ray photoevaporation rates, we use the formula of Owen et al. (2012), (3) (Using the formula of Picogna et al. (2019) instead results in very similar lifetimes.)Given the disagreement over X-ray photoevaporation rates discussed above, Eq.( 3) must be considered an upper limit, so the corresponding t life,g is a lower limit.The total photoevaporation rate is max ṀEUV , ṀX (w/ X-ray photoevaporation) .
We calculate the lifetimes for two cases: (i) EUV only and (ii) EUV + X-ray with Eqs.
(1) and (3).In Appendix E, we additionally discuss how the observed spread in L X (e.g., Güdel et al. 2007) can affect t life,g .Note that using the sum instead of max ṀEUV , ṀX does not make a difference, as ṀX dominates ṀEUV significantly.Note that these lifetimes apply only for small-grain-depleted disks and not to disks where FUV photoelectric heating drives photoevaporation.

Results
To provide an overview, we start with Figure 1, the estimated t life,g as a function of M * using the Chambers (2019) model for disks dominated by slow winds.Here the accretion is predominantly driven by the wind, which also extracts a significant fraction of the disk's mass (see Section 4.3 of the original paper for details).The figure shows the corresponding spectral type at each M * determined by T eff when the stars reach the zero-age-main-sequence (ZMAS) in the table of Kunitomo et al. (2021).The relation between the spectral type and T eff is taken from Pecaut & Mamajek (2013).The table of Kunitomo et al. (2021) ends before stars with M * ≤ 1 M ⊙ reach the ZAMS.Thus, for M * = 1 M ⊙ , we used the last T eff value in the table to determine the spectral type; the spectral types of M * < 1 M ⊙ stars are undetermined, so not specified in Figure 1.
In the EUV-only case (blue), the lifetime peaks at M * = 2 M ⊙ , or equivalently spectral types around A0. Lifetimes longer than 10 Myr are found only around the late intermediate-mass stars (2 M ⊙ ≲ M * ≲ 3 M ⊙ ), corresponding to spectral types between F0 and B8.The surface convective zones of these stars disappear at ∼ 1-10 Myr, and the magnetic-origin EUV emission is already weak by the time photoevaporation comes to dominate mass loss.The results indicate that primordial gas disks relatively easily survive to the ages of young debris disks around A-and late B-type stars.
Stars with M * ≥ 4 M ⊙ (hereafter early intermediate-mass stars) produce much stronger photospheric EUV emission owing to their high T eff , leading to shorter t life,g .For stars of Solar mass and below, the magnetic-origin EUV emission remains strong throughout the disk dispersal.X-ray photoevaporation (gray line in Figure 1) reduces the lifetime by 1-2 orders of magnitude for ≲ 3 M ⊙ compared to the EUV-only case.This is because L X is kept high at ≳ 10 31 erg s −1 , corresponding to ṀX ∼ 5 × 10 −8 -10 −7 M ⊙ yr −1 , until the disk completely disperses.However, as mentioned above, the adopted X-ray photoevaporation rates might be overestimated.While X-ray photoevaporation can explain observational wind diagnostics (e.g., Weber et al. 2020;Rab et al. 2022), such high X-ray photoevaporation rates are also known to be unfavorable for explaining the observed disk mass and accretion rates with viscous disk models (Sellek et al. 2020;Alexander et al. 2023) and magnetized disk models (Weder et al. 2023); smaller mass-loss rates as EUV photoevaporation alone are preferred.However, we do not rule out the effectiveness of X-rays in driving photoevaporation, as they indeed deposit energy to the gas, but with somewhat much smaller mass-loss rates than Eq.( 3).If so, the lifetimes of the EUV models set upper limits for low-mass stars.For the early intermediate-mass stars (≥ 4 M ⊙ ), t life,g matches the EUV-only case since the convective zone disappears before 1 Myr.
We then explore the dependence on the parameters of the MHD wind-driven accretion and mass loss across the fastwind and laminar-disk regimes in Chambers (2019).In the fast-wind case, the wind mass loss is negligible, and most of the accretion is driven by the wind, while in the laminar case, wind mass loss is substantial, and the wind drives almost all the accretion (Sections 4.2 and 4.4 of the original paper).The resulting t life,g -M * relation (Figure 2) is essentially the same as Figure 1, consistently peaking at 2 M ⊙ .Small differences are that in the fast-wind model, t life,g at 1.5 M ⊙ becomes comparable to 2 M ⊙ , and in the laminar model, the peak lifetime at 2 M ⊙ is reduced to ∼ 20 Myr.However, the assumption of no wind mass loss in the fast-wind model seems difficult to reconcile with the high detection rate of [O I] λ6300 low-velocity component (∼ 80%; Nisini et al. 2018), making it an unfavored PPD evolution scenario.The stellarmass-dependent trends are consistent across a wide parameter space, as shown by the population synthesis median in Figure 2 (cf.Appendix D).Overall, varying the disk and wind parameters yields a wide range of possible lifetimes, especially near 1.5 M ⊙ .We conclude that the average lifetime peaks at early A-type stars with t life,g > 10 Myr for smallgrain-depleted disks.
We emphasize that the lifetimes derived here apply only to small-grain-depleted disks.Figures 1 and 2 neglect any short-lived disks that might be dispersed by FUV photoevaporation.To quantify the general frequency of gas disks, it would be necessary to estimate PPDs' probability of becoming depleted in small grains.This is an issue for future studies.
We attribute the stellar mass dependence of t life,g mostly to the variation with M * in the EUV and X-ray luminosities and thus the photoevaporation rates.However, the variation of MHD effects with stellar mass may also play a role.The hard X-rays from low-mass stars (≲ 1 M ⊙ ) can ionize the disk atmosphere to increase the local effective α.This enhances MHD winds' effects, shortening the disk lifetime.The soft X-rays due to shocks originating from line-driven winds of early B-and O-type stars could have a similar effect and can also directly drive soft X-ray photoevaporation.Both these effects would strengthen the gas lifetime difference between the late intermediate-mass stars and the rest.
Currently, no evidence supports the presence of companions in known gas-rich debris disks, but the possibility still remains that they could harbor unseen low-mass, close-in companions.In these scenarios, the EUV and X-ray emission may be dominated by these companions, boosting photoevaporation rates.Thus, to yield long-lived gas disks, late intermediate-mass stars should not have such companions.

COMPARISON WITH DETECTION RATE OF GASEOUS DEBRIS DISKS
We now compare the modeled lifetimes from §3 with the population of debris disks detected in the CO millimeter, [C II] 157 µm, or [O I] 63 µm lines, focusing on the dependence on host star spectral type.The detections come from Supplemental Table 1 of Hughes et al. (2018).Only 3 out of 152 sources listed in the table have ages of < 10 Myr, with all being A0 stars.The age of HD 166191 (F star) is presented as 4 Myr, but we adopt an updated value of ∼10 Myr (Potravnov et al. 2018;Su et al. 2022).As they discussed, this sample combines studies with diverse sensitivity limits and selection criteria and thus is not well-suited for determining quantitative occurrence rates.However, the spectral-type dependence resembles that in a more-uniform sampling by Moór et al. (2017).
Keeping this limitation in mind, we replot the detections compiled by Hughes et al. (2018) versus stellar type in Figure 3.We show not only the detections in CO but also those in C II and O I, since the disks may retain these species mixed with hydrogen even in the absence of CO.We note that nine of twelve young CO-detected sources have CO mass Marino et al. 2020;Moór et al. 2017Moór et al. , 2019)).The three sources below this threshold are likely better explained by secondary origins, based on the short photodissociation lifetimes for primordial CO at such low gas masses.For F stars, one of the two young COdetected sources is above this mass threshold and is early F (F2/3V); the one below the threshold is late F (F5/6V).The decreasing trend of the incidence from A to FGK stars is likely real but not due to a sensitivity bias.A comparison of debris disks with similar fractional luminosities (Tables A1  and A2 of Marino et al. 2020) indicates that only 2 out of 9 sources exhibit gas (at levels ≲ 10 −4 M ⊕ ) among FGK stars, while 11 out of 17 A-type stars display gas presence (often exceeding 10 −3 M ⊕ ).
Figure 3 shows several interesting consistencies with the modeled t life,g -M * relation.The detection rate is highest for early A-type stars (more than half the stars with gas detected are early A-type, A0-A3).Almost all the gasrich debris disks with M CO > 10 −3 M ⊕ are hosted by early A-type stars.The exception is HD 121191 (3 × 10 −3 M ⊕ , Moór et al. 2017), which is type A5 but does not necessarily disagree with our model considering the range in t life,g at M * ∼ 1.5 M ⊙ .Our model naturally explains younger (< 50 Myr) gaseous debris disks' higher numbers around Atype stars than around both lower-mass and O-and early B stars.The decreasing trend toward higher-mass stars could be due to both photoevaporation and intense CO photodissociation.For A stars, the weak photodissociating fluxes would also be helping the long survival of CO.Overall, within the limitations of the sample, the consistencies with the model are remarkable.
In contrast, the gas detected in debris disks around late Fstars appears more compatible with secondary origins, given the short gas lifetime around those low-mass stars in our model.Also, a few debris disks have multiple lines of evidence supporting a secondary origin, among them the A5-A6 star β Pic (e.g., Dent et al. 2014;Greaves et al. 2016;Iwasaki et al. 2023;Cataldi et al. 2023).We thus consider primordial and secondary origins not mutually exclusive.We stress that our hypothesized disk evolution does not necessarily rule out a secondary origin for some gaseous debris disks, even among A-type stars younger than 50 Myr.Determining the gas origins in individual sources is a challenge for the future.However, the results we present here suggest the high overall frequency of gas-rich debris disks among A-type stars can come from the primordial origin.

SUMMARY AND OUTLOOK
The substantial amounts of gas found in some debris disks could either be remnants of the primordial protoplanetary disk, or originate in the secondary release of volatiles from bodies resembling comets or asteroids.While the primordialorigin scenario was classically considered unlikely, we here demonstrated primordial gas can survive to the ages of the known gaseous debris disks, using semi-analytic models to derive lifetimes for the gas in disks that have been depleted in small grains (≲ 0.01 µm) so that FUV photoelectric heating is ineffective.The resulting lifetimes are longest around stars of 2 M ⊙ ≲ M * ≲ 3 M ⊙ , whose evolution switches off their surface convection at 1-10 Myr before the disk loses the bulk of its gas.Without a surface convective zone, these stars lack surface magnetic activity and thus have low luminosities in the EUV and X-ray bands that drive the photoevaporation of the gas.
The gas lifetime consistently exceeds 10 Myr for the late intermediate-mass stars only.The shorter lifetimes of the gas disks around lower-mass stars (≲ 1 M ⊙ ) are due to strong EUV and X-rays from the chromosphere, transition region, and corona, while those around higher-mass stars (≳ 4 M ⊙ ) are due to the high photospheric temperatures and resulting EUV luminosities.
Our model predicts a higher frequency of gaseous debris disks around A-and late B-type stars compared to types both earlier and later.This is qualitatively consistent with the observed population, motivating further investigation of the primordial-origin scenario.Longer gas disk lifetimes for A-stars could be a factor in the giant planet occurrence rate peaking at 1.7-1.9M ⊙ (Reffert et al. 2015;Wolthoff et al. 2022), though more detailed modeling of the dissipating disks would be needed to assess the possibility of gas giant formation.Likewise, our results could relate to the increasing occurrence of transition disks with stellar mass (van der Marel & Mulders 2021).Further investigations are needed to evaluate whether these transition disks might serve as precursors to gaseous debris disks.
The primordial-origin picture allows the existence around late intermediate-mass stars of long-lived gas-rich disks with very weak dust emission, which we term "phantom" disks.Finding such phantom disks would indicate that the abundant CO gas in at least some young debris disks is primordial.Searches for gas-rich, dust-poor disks especially around Aand late B-stars aged 10-50 Myr could thus help determine the range of evolutionary paths followed by planet-forming  2018) Table 1.These trends are consistent with the M * -dependence of our predicted lifetimes in Figure 1.
disks.Note that the phantoms are different from the relic disks, purely viscous, dusty disks with central cavities larger than 100 au surviving > 10 Myr around Solar-type stars owing to very weak X-ray emission (Owen et al. 2011).
The fraction of PPDs achieving the small-grain-depleted state and the critical level of depletion required to inhibit FUV photoevaporation are open questions.Further hydrodynamics investigations are needed considering a range of stellar masses and FUV/X-ray luminosities.We anticipate that PPDs with PAH/small-grain abundances below ∼ 0.1-1% of the interstellar level (see §3) may have reached this state and could be potential precursors of gas-rich debris disks around intermediate-mass stars.PPDs with PAH/smallgrain abundances above the critical level would have FUVdriven H 2 winds extend to ∼ 10-100 au with temperatures exceeding a few hundred Kelvin and mass-loss rates of ∼ 10 −9 -10 −8 M ⊙ yr −1 (Nakatani et al. 2018a,b;Komaki et al. 2021).Searching for these molecular winds could also help examine the depletion level in PPDs and very young (< 10 Myr) hybrid disks (e.g., HD 141569), including hybrid disk candidates (Iglesias et al. 2023), around weak X-ray emitters.
The point of this work is that gas disks can persist beyond 10 Myr, a necessary condition for the primordial-origin scenario.A caveat is that our 0D approach cannot unambiguously predict the detailed surface density evolution and the alignment with other observational characteristics of inner disks, such as accretion and wind mass-loss rates (Fang et al. 2018(Fang et al. , 2023)).Future studies involve 1D simulations to accurately follow the radial profile evolution and make comparisons with those observables to evaluate the proposed evolutionary pathways.Nevertheless, we expect the conclusion on the long lifetimes around intermediate-mass stars will remain consistent even in more advanced models, as the key determinant of outer gas disk lifetimes, which potentially manifest as gaseous debris disks in observations, is the photoevaporation rates set by the stellar emission rates.The speed of final-stage disk dispersal would not be significantly influenced by how the inner disks clear.Thermochemical modeling would also help understand how the CO and C masses evolve in smallgrain depleted disks and examine if the final values reproduce observations (Iwasaki et al. 2023;Cataldi et al. 2023).For a broader perspective, it is worth investigating what fraction of PPDs in general can be long-lived considering environmental factors like external photoevaporation, tidal truncation, and late infall.

ACKNOWLEDGMENTS
We are grateful to the anonymous referees for their useful comments.We thank Masanobu Kunitomo and Akimasa Kataoka for their practical and insightful comments on this study.RN is supported by the Japan Society for the Promotion of Science (JSPS), Overseas Research Fellowship.SM is supported by a Royal Society University Research Fellowship (URF-R1-221669).This research was performed in part at the Jet Propulsion Laboratory, California Institute of Technology, under contract 80NM0018D0004 with the National Aeronautics and Space Administration and with the support of the NASA Exoplanets Research Program through grant 17-XRP17_2-0081.This work was strongly benefited from the Core2disk-III residential program of Institut Pascal at Université Paris-Saclay, with the support of the program "Investissements d'avenir" ANR-11-IDEX-0003-01.
Software: Numpy (Harris et al. 2020), Matplotlib (Hunter 2007), Astropy (Astropy Collaboration et al. 2013, 2018, 2022), SciPy (Virtanen et al. 2020) Photoevaporation rates vary according to stellar evolution, where the luminosities and spectra change by orders of magnitude.Also, if Ṁph gets smaller at some point by, e.g., reduced photoelectric heating due to small-grain depletion, the lifetime is extended.This is the basic idea of our model (and Nakatani et al. 2021) to explain a long-lived gas disk in the context of the primordial origin scenario.

B. AN ANALYTIC DISK EVOLUTION MODEL WITH MHD EFFECTS
The effects of magnetohydrodynamics (MHD) winds are essential in disk evolution.They can remove not only the mass but also the angular momentum of the disk.Chambers (2019) derived an approximate analytic solution for disk evolution with MHD wind effects.Eq.( A1) is modified to where v w is the inward velocity induced by the disk wind, and Σw is the surface mass-loss rate due to MHD winds.Using a reference radius R 0 , the surface mass-loss rate is characterized by a constant parameter K as where v 0 is the inward velocity at R 0 , f w is the fraction of v 0 induced by MHD winds, v esc is the tangential velocity of the winds, and Ω is the Keplerian orbital frequency.The second equation is derived from the conservation of angular momentum.Chambers (2019) modified Eq.(B8) to make it analytically solvable (Eq.( 12) of the original paper) and found analytic solutions for the modified equation (Eq.( 36)-( 39) of the original paper).We omit to write the solutions explicitly here to avoid complexity and refer the readers to the original paper.The point of this analytic solution is that the surface density and temperature evolutions are uniquely determined for a given set of model parameters v 0 , f w , K, the temperature at R 0 : T 0 , and the surface density at R 0 : Σ 0 .
For example, f w = 0 means purely viscous disk, and choosing a high f w (≤ 1) indicates a disk that accretes mainly magnetically.The parameter K can be set to zero with a nonzero f w .The underlying assumption, in this case, is negligible mass-loss rates due to MHD winds.Small K cases are termed fast wind.On the other hand, using a high value for K (≤ 1) means that a significant fraction of the accreting gas flows out of the disk before reaching the host star.This case is called a slow wind.
The analytic formula of Chambers (2019) can overall well agree with the true solutions of Eq.(B8) regardless of the approximation taken to make the equation analytically tractable.The exception is when the wind mass-loss rate is large; the analytic formula underestimates the surface density and outer spreading at the later stage.It also fails to predict the positive slope of Σ created by the wind in the inner ∼ 1 au region.Nevertheless, the errors between the analytic formula and the true solution are small compared to the large uncertainties in the long-term evolution of disk evolution with the MHD effects.We, therefore, use the formula as a valid estimate for time-dependent disk mass evolution.

C. ACCRETION-GENERATED EMISSIONS
We have omitted accretion-generated UV and X-ray emissions from our model.This extra radiation component can dominate the magnetic and photospheric components, leading to accelerated photoevaporation and reduced t life,g .We now examine if this omission affects the study's conclusions-specifically, the peak lifetime beyond 10 Myr at M * ≈ 2 M ⊙ (cf.Figures 1 and 2).Since the EUV models potentially account for the observed higher occurrence of gas-rich debris disks around A stars, our focus remains exclusively on the EUV models.
The exact portion of accretion luminosity (L acc ≈ GM * Ṁacc /R * ) contributing to EUV radiation is uncertain.For simplicity, we assume ∼ 4% is processed into EUV emission, aligning with analogous estimation for FUV (Gorti & Hollenbach 2008) and X-ray (Hartmann et al. 2016).This can be represented as: Here, the average EUV photon energy is approximated to match that of 9000 K blackbody radiation, 14.5 eV.The total EUV emission rate is now Φ EUV = Φ EUV,ph + Φ EUV,mag + Φ EUV,acc .In practice, accretion-generated EUV radiation may terminate once the disk develops an inner cavity, when photoevaporation and accretion rates balance at ∼ 1 au (the critical radius for photoevaporation).However, as the cavityopening time is not determined unambiguously within the framework of our model, we let the accretion-generated EUV emission last until the disk dispersal; this scenario sets the lower limit for t life,g .Figure 4 shows the resulting lifetimes.The accretiongenerated EUV dominates Φ EUV,mag for the weak-and fastwind disks especially at the early time (≲ 1 Myr), while laminar disks show minimal total emission rate enhancement.Nevertheless, the photoevaporation rate remains orders of magnitude lower than the accretion and/or MHD wind massloss rates at the early stages, with no substantial impact on disk evolution.Except for fast-wind disks, the accretiongenerated component becomes negligible as the accretion rate drops in 1 Myr ≲ t ≲ 10 Myr.This occurs prior to the on- Gas disk lifetime: t life, g (Myr) G6V F1V A1V B9V B8V B5V B3V Slow wind Fast wind Laminar set of the photoevaporation-dominant epoch at t = t 0 , leading to minimal shortening of t life,g for the slow-wind and laminar disks.
In contrast, fast-wind disks can maintain relatively high accretion rates (∼ 10 −10 -10 −9 M ⊙ yr −1 ) even at t = 10 Myr, with accretion-generated EUV dominating over (or comparable to) the other components throughout dispersal.This impacts t life,g at M * = 1.5-2.5 M ⊙ , as Φ EUV,acc is significantly exceeds Φ EUV,ph + Φ EUV,mag .However, the prevalence of jets and winds observed widely (e.g., Pascucci et al. 2022) suggests that the fast-wind scenario is unlikely to be common among PPDs.It is also difficult to align with the high detection rate of [O I] λ6300 low-velocity component (Nisini et al. 2018).If termination of accretion-generated emissions due to cavity opening is considered, the peak at M * ∼ 2 M ⊙ may remain even for fast-wind disks.
Overall, incorporating the accretion-generated component can moderately reduce t life,g and lessen the stellar-mass dependence for disks with no wind mass loss.Nonetheless, not all disks would evolve as fast-wind disks, supporting the highest plausibility of long-living (> 10 Myr) gas disks around late intermediate-mass stars.Further investigation is needed to accurately determine t life,g , Φ EUV,acc , and its termination time.

D. DISK PARAMETER DEPENDENCES
We explore how t life,g depends on various disk parameters, including the initial disk mass (M disk,0 ), disk mass scaling law (M disk,0 ∝ M β * ), viscous α, and initial disk cut-off radius (r exp ).These parameters are treated as Monte Carlo variables in population synthesis to span a wide parameter space.
The fiducial v 0 values are 30 cm s −1 for slow-wind and fastwind disks and 10 cm s −1 for the laminar disk.We explore 10 cm s −1 ≤ v 0 ≤ 30 cm s −1 , as values beyond this range lead to unrealistic accretion rate variations.The wind mass-loss parameter K is varied in 0 ≤ K ≤ 1.The initial cut-off radius is tested in 10 au ≤ r exp ≤ 40 au; the fiducial is r exp = 15 au.For reference, Clarke et al. (2001) and Kunitomo et al. (2020) adopt r exp = 10 au and 30 au, respectively; the time evolution of disk radius in the model by Trapman et al. (2022) aligns with the characteristic radii of observed CO disks in Lupus and Upper Sco when r exp ≈ 20 au and M disk,0 = 0.1 M ⊙ .
Monte Carlo variables assume uniform densities in linear space, except for C disk,0 and α, which follow log-uniform sampling.We perform in total 90000 runs: 5000 parameter sets times 9 different stellar masses for each of the EUV and EUV+X-ray models.Figure 5 illustrates the resulting probability densities and cumulative probabilities for t life,g , reflecting the trends seen in Figures 1 and 2.
The lifetime is mostly unaffected by β and K.For the EUV models and EUV+X-ray models with M * ≥ 3 M ⊙ , M disk,0 has limited influence.However, for EUV+X-ray models with M * ≤ 2.5 M ⊙ , smaller initial disk masses (C disk,0 ≤ 0.1) result in t life,g ≲ 5 Myr.
Increasing r exp moderately extends t life,g .In the EUV models, t life,g tends to exceed 10 Myr for the runs with M * ≤ 1.5 M ⊙ and r exp (> 25 au).Conversely, the most runs with M * = 2-3 M ⊙ and t life,g < 10 Myr have r exp ≲ 20 au.When r exp is expanded while keeping M disk,0 constant, it decreases the initial surface density in the inner part, leading to greater mass distribution in the outer disk.This leads to slower mass loss via accretion and MHD winds.
Viscous α and v 0 also have a moderate impact on the EUV models, with smaller α and higher v 0 causing shorter lifetimes due to enhanced MHD-driven mass loss and accretion.The longest-t life,g population (> 10 2 Myr) at M * = 1.5-3M ⊙ has α > 10 −3 and v 0 < 15 cm s −1 , indicating a slow decay of the accretion rate.Unrealistically long lifetimes (∼ 10 3 yr) can occur when α is increased to ∼ 10 −2 .Hence, when the disk exhibits a characteristic of a viscous disk relatively strongly, t life,g increases significantly.
In summary, higher α and larger r exp can extend t life,g for specific cases.The other parameters exert minor impacts.Regardless of the parameter dependences, overall trends remain consistent: Late intermediate-mass stars in the EUV models consistently exhibit t life,g > 10 Myr, peaking at M * = 2 M ⊙ and declining for lower and higher masses.While the fraction of long-living (t life,g > 10 Myr) disks for M * ≤ 1 M ⊙ , most of which have r exp > 25 au, is apparently higher than observations (Figure 3), the actual CO disk lifetime is shorter than t life,g .This is because complete CO photodissociation occurs before disk dissipation, when M disk reduces to a certain threshold mass.In addition, X-ray photoevaporation likely causes mass loss at somewhat much smaller rates than Eq.(3), implying that the the lifetimes of the EUV models should be interpreted as upper limits for low-mass stars.
Thermochemistry modeling is needed to interpret these results fully, along with the potentially crucial influence of atypically weak stellar photodissociating fluxes around A stars.It is important to highlight that our findings primarily demonstrate the plausibility of gas disk survival beyond > 10 Myr, a necessary condition for supporting the primordial-origin scenario.Additionally, our model qualitatively explains the observed incidence of gaseous debris disks versus stellar mass.E. X-RAY LUMINOSITY SPREAD Observed X-ray luminosities are known to have a large spread (e.g., Güdel et al. 2007).Accounting for the variability in X-ray luminosity would introduce scattering in the derived lifetimes.In the EUV+X-ray models, amplifying X-ray luminosity results in shorter lifetimes while reducing it aligns the lifetimes with those of EUV models.Within the EUV models, the spread in luminosity could influence lifetimes for M * < 1 M ⊙ , as their lifetimes are linked to Φ EUV,mag .How-ever, combining Eqs. ( 1) and ( 2), the scaling of EUV photoevaporation rate with L 0.33 X places a relatively restrained impact on lifetimes.Furthermore, the adopted X-ray luminosity encompasses representative observational values (Fig. 6 of Kunitomo et al. (2021)).It suggests that derived lifetimes remain relatively stable even when incorporating luminosity spread with a plausible distribution of L X .Consequently, the inclusion of luminosity dispersion is unlikely to significantly alter the core finding of this study.
Taking into account the spread of L X along with its time evolution would require solving stellar evolution across a wider parameter space than covered by Kunitomo et al. (2021).

Figure 1 .
Figure1.Estimated gas disk lifetime tlife,g as a function of M * .The blue line indicates tlife,g for the EUV-only case.The gray line adds X-ray photoevaporation.The spectral type corresponding to each M * is presented on the top axis.Note that these lifetimes apply only for small-grain-depleted disks and not to disks where FUV photoelectric heating drives photoevaporation.

Figure 2 .
Figure 2. Same as Figure 1.Solid, dashed, and dot-dashed lines show lifetimes of slow-wind, fast-wind, and laminar disks, respectively.The dotted lines show the median lifetimes of the population synthesis in Appendix D.

Figure 3 .
Figure 3. Detection rates of gaseous debris disks are highest among young A-type stars.The tracers are CO millimeter emission (left), C II 157 µm (middle), and O I 63 µm (right).The light-green histograms show the all-survey count and are superposed by dark-green histograms indicating the detection count, for sources younger than 50 Myr.The light-yellow and dark-yellow histograms are sources older than 50 Myr.The horizontal axis gives the stellar type.Integrated counts are shown for M, K, and G stars.All data are from Hughes et al. (2018) Table1.These trends are consistent with the M * -dependence of our predicted lifetimes in Figure1.

Figure 4 .
Figure 4. Derived tlife,g with accretion-generated EUV (blue).The lifetimes of the fiducial runs are also plotted for comparison (light blue).

Figure 5 .
Figure 5. Resulting population.The blue and gray histograms are the probability densities for the EUV and EUV+X-ray models, respectively.Corresponding normalized cumulative probabilities are shown by the solid lines with the same colors.