Dust Accumulation near the Magnetospheric Truncation of Protoplanetary Discs. II. The Effects of Opacity and Thermal Evolution

Dust trapping in the global pressure bump induced by magnetospheric truncation offers a promising formation mechanism for close-in super-Earths/sub-Neptunes. These planets likely form in evolved protoplanetary discs, where the gas temperature at the expanding truncation radius become amiable to refractory solids. However, dust accumulation may alter the disc opacity such that thermal evolution is inevitable. To better understand how thermodynamics affects this planet formation pathway, we conduct a suite of local dust evolution simulations in an idealized inner disc model. Our calculations take into account self-consistent opacity-dependent temperature changes as well as dust evaporation and vapour condensation. We find that disc thermal evolution regulates dust growth and evolution, discouraging any accumulation of small particles that drives the increase of opacity and temperature. Significant retention of dust mass takes place when the disc environments allow runaway growth of large solids beyond the fragmentation barrier, where small particles are then swept up and preserved. Our results further validate dust accumulation near disc truncation as a promising mechanism to form close-in planets.

Several formation channels for such super-Earths have been proposed over the years, with a broad consensus that super-Earths are most likely formed out of materials originated from the external disc (i.e., ex situ formation) (Morbidelli & Raymond 2016, see also their Fig. 5 and references therein; Izidoro et al. 2021).The remaining key question is when and where growth takes place during inward migration of solids.If the initial growth of planetary embryos is efficient, they may form in the external disc and then migrate into the tight orbits (e.g., Ida & Lin 2010; Kley & Nelson 2012).Pebble accretion may help embryos grow much further during and after migration (e.g., Lambrechts et al. 2014;Bitsch et al. 2015;Johansen & Lambrechts 2017;Lambrechts et al. 2019).
Alternatively, rapid inward drift of dust into traps induced by pressure bumps at short-period orbits may retain solid mass fast ★ Contact e-mail: rixin@berkeley.edu† 51 Pegasi b Fellow enough to form embryos and supply their further growth.The dead zone inner boundary (DZIB) may serve this purpose, but recent disc modeling suggested that the DZIB either locates far out or may not induce a pressure bump at all (e.g., Chatterjee & Tan 2014;Hu et al. 2016Hu et al. , 2018;;Jankovic et al. 2021Jankovic et al. , 2022)).
In Li et al. (2022, hereafter Paper I), we proposed and demonstrated that the global pressure maxima induced by magnetospheric truncation in evolved discs may trap dust efficiently, providing a promising formation pathway for close-in super-Earths.Specifically, two subpathways for effective dust retention operate at different parameter regimes.When the turbulent relative velocity between dust (correlated with the Shakura-Sunyaev -parameter (Shakura & Sunyaev 1973) and the gas sound speed  s ) is relatively low and the threshold velocity for fragmentation  f is high (e.g., due to the sticky surface layer on dust at a high temperature), or when large solids can be rapidly delivered from the outer disc, a fraction of particles may grow fast and surpass the fragmentation barrier.This process triggers an accelerated accumulation of progressively larger particles, eventually reaching the typical size of planetesimals (referred to as the breakthrough scenario).Conversely, in a different scenario characterized by a high turbulent relative velocity and a low threshold velocity for fragmentation, coagulation growth beyond small grains is highly suppressed.However, when the dust supply from the outer disc is substantial enough to counteract the removal through funnel flows, the overall dust mass can still gradually rise, eventually reaching a point where gravitational instability (GI) produces planetesimals (referred to as the feedback+GI scenario).
Dust retention at the truncation radius  T , which converges with the corotation radius  Co at low accretion rates (e.g., Long et al. 2005;Bouvier et al. 2007;Romanova et al. 2008;Thanathibodee et al. 2023;Zhu et al. 2023;Pittman et al. in prep.), is of great astrophysical interest by itself even if the accumulated dust does not end up in large terrestrial planets.Dust features at these location are one of the favourable explanations proposed for "dipper" stars (e.g., Bouvier et al. 1999;Cody et al. 2014;Ansdell et al. 2016;Stauffer et al. 2017;Robinson et al. 2021;Capistrant et al. 2022).Moreover, small rocky bodies formed from the accumulated dust near  Co may serve as sources of materials that produce complex periodic variables, a new class of variable star (Bouma et al. 2023).
In Paper I, we adopted an idealized model with a static gas disc for simplicity, which allowed us to focus on a preliminary exploration of dust evolution.In this work, we take into account the dynamical thermal evolution of the gas disc caused by dust accumulation and the resulting changes in dust opacity.Our goal is to understand how opacity-driven thermal evolution affects dust evolution as well as this new planet formation pathway.
The paper is organized as follows.Section 2 first describes our dust evolution model, including our new considerations on dust opacity and disc thermal evolution.Section 3 then presents our simulation results.Finally, Section 4 summarizes our key findings and discusses the important implications.

DUST EVOLUTION MODEL WITH OPACITY AND EVAPORATION
We use the implicit dust coagulation-fragmentation code Rubble (Paper I) to model the evolution of dust size distribution, dust surface density, and the thermal evolution between dust and vapour at the dust accumulation radius  accu (i.e. the truncation-induced global pressure maximum) , in a radiative disc model initialized based on Ali-Dib et al. (2020).In this work, Rubble has been upgraded to run on graphics processing units (GPUs) by utilizing PyTorch (Paszke et al. 2019), which greatly speeds up our simulations.
In this section, we briefly reiterate our base numerical scheme to model dust size evolution.Sections 2.2 and 2.3 then describe how Rubble calculates dust opacity and solid evaporation and condensation due to temperature changes.Finally, Section 2.4 details the parameters for our simulations.

Basic Scheme
Rubble solves the Smoluchowski equation where  () ≡ d/d is the vertically integrated dust surface number density in a mass bin,  (,  ′ ,  ′′ ) is the full kernel that consists of coagulation kernel and fragmentation kernel, which elucidates the likelihood of particles of masses  ′ and  ′′ colliding with each other to produce target mass of .
To relate the total dust surface density Σ d and  (), we define the vertically integrated dust surface density distribution per logarithmic bin of grain radius (), and where where d/d log  is the vertically integrated dust surface number density in a logarithmic mass interval and is the quantity that our implicit code actually evolves.

Opacity and Energy Equation
We calculate the dust opacity from a full size distribution using the prescription from Chen et al. ( 2020), modified from the singlespecies opacity formula of Ormel (2014).For each size bin in the dust distribution with characteristic size   , the metallicity is   =   d log /Σ g where Σ g is the gas surface density.The total opacity consisting of the dust and gas opacity is (in cgs units) where  g and  denote the midplane gas density and temperature,   is the dust geometric opacity, and   is the efficiency coefficient.The first term in the equation above is the gas opacity that only dominates at high temperature (≳ 2000K or when the gas is dust-free), and the second term denotes the dust opacity that become dominant below 2000K.The gas opacity which dominates at large temperature is approximated by analytical expressions from Bell & Lin (1994).It will be dominated by dust opacity below 2000K1 .Since this opacity matters in the calculation of radiative energy flux of the gas profile (optical depth  = Σ g /2 ), we need to scale the dust cross section to that of per unit gas mass.Thus, where   is the dust internal density, and   = 2  / max , where  max = 0.29 cm/ is the peak wavelength of blackbody radiation from Wien's law.
The disc thermal evolution is then driven by the following energy equation where   = 1.7×10 8 erg g −1 K −1 is the specific heat of hydrogen gas,  + is the viscous-heating source term (Pringle 1981) at the corotation radius  Co , where the the viscosity  =  s  and Ω = 2/ ★ , where  is the gas scale height and  ★ is the stellar spin period.Ohmic heating near  Co is negligible and the stellar irradiation is also weak in the inner disc region (Garaud & Lin 2007).The last two terms on the right hand side are the cooling source terms, where  −,rad = 2 SB  4 s is the radiative cooling term and  −,d2v is the dust-to-vapour evaporation cooling term.For the first cooling term,  SB is the Stefan-Boltzmann constant and  s is the effective surface temperature, which is approximately related to the midplane temperature as where  is the optical depth and the formula above takes into account both optically thin and optically thick situations.For the second cooling term,  −,d2v is proportional to the dust-to-vapour mass conversion rate per unit area Σ d2v , and  0 is the latent heat needed during dust evaporation.

Evaporation and Condensation
To model the dust evolution that includes evaporation and condensation, we assume most solids are silica (e.g., SiO 2 ) and develop a semi-implicit scheme based on the following equations (Schoonenberg & Ormel 2017) where the dust-to-vapour conversion rate is Σ d2v =  e Σ d and the vapour-to-dust conversion rate is where Σ v is the vapour surface density and the coefficients ) are the evaporation rate and condensation rate, respectively, where  SiO 2 is the mean molecular weight for SiO 2 ,  B is the Boltzmann constant, and  eq is the saturated or equilibrium pressure given by the Clausius-Clapeyron equation where   and  eq,0 are constants depending on the species.To solve the evaporation and condensation for dust with a size distribution, we re-write Equations 9 and 10 in a discrete, semi-implicit way where  and  + 1 denote two consecutive time steps and Δ is the length of a time step.We further substitute Σ d with    d log  from each mass bin and organize the equations into a linear system where  ′ c =  c    d log .Equation 18 represents a series of equations with  = 1, • • • ,   , where   is the total number of dust species in the model.The evolution of evaporation and condensation can be then solved in the matrix format Since the condensation rate Σ v2d depends on the remaining available solid materials, we impose a tiny surface density floor for dust particles smaller than 0.01 cm such that vapour can condense back to solids when the temperature drops after a full melt-down.We adopt this idealized prescription due to the non-uniform nature of nucleation in such a turbulent disc environment and its potential enhancement by non-linear effects within local eddies (Tang et al. 2022).
Generally speaking, the timescale for thermal evolution is often much shorter than the timescale needed for dust size distribution evolution.However, our calculations are only computationally feasible with time steps based on the latter.Appendix A demonstrates that our thermal module is numerically robust against these time steps.Furthermore, particles that gain or lose too much mass may leave the original size bin and shift to another one, which is not taken into account in this work.However, we carefully examine the net effects of evaporation and condensation in our simulations and find that the consequent relative mass changes in all bins are only of the order of 1% at each time step, much less than the fractional change needed to shift bins (∼ 32%) in our setup.

Numerical Setup
Table 1 summarizes the physical and numerical parameters for our dust evolution simulations, where we adopt the same disc model outlined in Paper I as the initial conditions (see Sections 2 -4 therein for detailed descriptions).Similarly, the dust evolution is then determined by the -viscosity, the fragmentation velocity threshold  f , the dust-to-gas ratio of the accreted materials from the external disc  supp , and the maximum supplied particle size  supp,max .
In this work, we depart from the conventional approach of initializing all solids with a monodispersed distribution at the micron scale, which results in artificially elevated dust opacity.Instead, the initial dust size distribution follows the MRN power-law distribution, ranging from 10 −4 to 10 −2 cm.This choice enables us to maintain self-consistency with our disc model by normalizing the overall dust surface density to be ≲ 1% of the gas surface density and achieve a desired initial opacity of  = 1.0 cm 2 /g.For the MRN distribution, the value of  is relatively insensitive to the initial disc temperature.
The vapour surface density is initialized with the equilibrium value that can be derived from 14 given a certain temperature where  rad,0 is the initial temperature at disc midplane and  s denotes gas sound speed.During the evolution, we further assume that the accretion funnels remove vapor in the same manner as gas accretion, i.e., Σ v ≡ Σ v Σ g /Σ g .Our models take into account all the physical ingredients considered in Paper I (e.g., mass transfer, feedback effects), besides the newly included dust thermal evolution.We again only evolve all the models for 10 5 yr and consider  final ≳ 1 as the criterion for significant solid accumulation.
In this work, we focus on the temperature evolution of the gas disc following Equation 7 and maintain a constant background gas density Σ g at the accumulation site for simplicity.In reality, the gas surface density profile could be altered on the viscous timescale because the inflow of materials tends to smooth out discontinuities in the disc accretion rate , which is directly proportional to Σ g .However, this alteration occurs on a timescale significantly longer than the thermal evolution and so we simplify our model by neglecting the time evolution of  and Σ g .Furthermore, the assumption of a constant  Σ g breaks down in the limit where the vapour surface density Σ v is a significant fraction of Σ g , which however is never the case in all of our results.

RESULTS
We conduct a suite of local dust evolution simulations, incorporating self-consistent opacity-driven disc thermal evolution, with parameters listed in Table 1.Section 3.1 first presents a set of representative cases  .The evolution of (from top to bottom) dust-to-gas surface density ratio, vapour-to-gas surface density ratio, local disc temperature, and total opacity for the same models shown in Fig. 1, color-coded by  and line-stylecoded by  f .and elucidate the interplay between the collisional and thermal dust evolution.Section 3.2 then surveys disc conditions that are required to achieve effective dust retention.

Fiducial Model
Figs. 1 and 2 compare the time evolution of dust size distribution, the dust/vapour-to-gas ratio, the gas disc temperature, and the total disc opacity between our models with different  and  f , where  supp,max is fixed to 1000 cm and  supp is fixed to 0.01.
We find that almost all models experience significant thermal evolution in the very first year, before reaching a steady state or when the external dust supply makes a difference.For cases with a higher  and thus a relatively lower initial disc temperature  rad,0 ,  dust fragments easily and quickly increases the total opacity  and disc temperature , unless  f is high enough.Conversely, in cases with a lower  (i.e., higher  rad,0 ), dust particles may melt down faster than fragmentation such that the disc temperature does not change substantially, particularly with a lower  f where smaller dust dominates.
Efficient breakthrough and sweep-up growth happen in models with a larger  f , leading to significant dust accumulation.In these cases, the higher fragmentation velocity threshold promotes coagulation and allows the largest supplied particles to pile up, which eventually results in runaway growth and breakthrough.Once there are enough large particles (> 10 3.5 cm), sweep-up growth via mass transfer dominates the dust evolution and transports almost all the supplied dust mass into very large particles (Windmark et al. 2012).
Although similar dust retention processes have been seen in Paper I, their dynamic impact to the disc thermal evolution is shown for the first time in this work.We find that the initial pile up of solids near  supp,max (at  ∼ 20 yr) cause a deficit of sub-mm dust and consequently a decrease in the total opacity and a sharp drop in the disc temperature.However, a considerably lower disc temperature hinders dust loss due to fragmentation (and funnel flows) and inflates the dust surface density across the entire size distribution, which then in response leads to an upturn for  and .After breakthrough happens and most of the solid mass resides in the very large particles, the changes in  and  become subtle.Furthermore, the final disc temperature always shows a much weaker dependency on  in comparison to the initial temperature, regardless of  final and accompanied by different extent of offset in the opacity from the initial value of unity.
An accelerated external disc dust supply augments the overall reservoir of solid mass, which always facilitates dust accumulation as seen in Paper I. However, when thermal evolution is taken into consideration, this increased supply also elevates the abundance of fine dust particles that raise disc opacity.As illustrated in Figure 2, models with  supp = 0.05 consistently achieve a higher opacity and thus a higher disc temperature when compared to their counterparts in Fig. 2. Therefore, although more cases can generate substantial dust accumulation, the dust-to-gas ratios eventually converge to a plateau after ∼ 10 4 yr, representing a new balance between the thermal evolution and the dust size distribution evolution.

Threshold for significant dust accumulation
Figs. 4 and 5 present our survey results on parameters listed in Table 1 and show which runs retain substantial dust mass and which do not via the final dust-to-gas ratio  final .We find that the maximum-supplied size  supp,max and the fragmentation velocity threshold  f largely determine the dust evolution, consistent with our findings in the last section.When  supp,max = 1000 cm, substantial dust accumulation happens when  f ≳ 5.62 m s −1 if  supp = 0.01, or  f ≳ 3.16 m s −1 if  supp = 0.05.When  supp,max ≲ 300, only the corner case ( = 10 −3 ,  f = 10 m s −1 ) can efficiently retain dust.
As discussed in the previous section, for  supp,max = 1000 cm, the number of cases with  final ≳ 1 increases with  supp .However, the values of  final in these runs with  supp = 0.05 are consistently lower than those in their counterparts with  supp = 0.01, due to the higher opacity/temperature caused by the richer supply of fine dust grains.Therefore, with thermal evolution taken into account, a larger dust supply may not always yield a positive impact on all aspects of the dust accumulation process.It is possible that a even higher  supp could fully suppress dust retention, which is however beyond the scope of this work given that the choice of  supp = 0.05 is already quite large.
Compared to the previous results that neglected thermal evolution in Paper I, significant accumulation of mm-cm dust through the feedback+GI scenario (i.e., most effective at the low  high  f parameter space) becomes unviable, as small dust particles contribute most to the opacity/temperature enhancement, which in turn quench dust retention.Instead, the breakthrough scenario exhibits a better chance of accumulating solids, with a strong preference for a higher  f and a moderate preference for a higher , which provides a lower initial disc temperature.

SUMMARY AND DISCUSSIONS
In this sequel to our prior work in Paper I, we delve into the intricate dynamic synergy between the local dust distribution evolution and the local disc thermal evolution near the truncated inner boundary of protoplanetary discs.In particular, we improve our model by taking into account the self-consistent intertwined changes, including the change in disc temperature due to the change in opacity, the change in opacity due to the change in dust size distribution, and the evaporation and condensation of solids at all sizes due to the evolving disc temperature.
Our main findings are as follows: • We find that the dust thermodynamics plays an important role in dust size evolution.Since the total opacity is sensitive to the amount of small (∼ mm) dust particles, the self-consistent thermal evolution fully suppresses the feedback+GI scenario that builds a mass reservoir in small dust particles until the dust-to-gas ratio exceeds order unity.
Thus, significant dust accumulation only happens in the breakthrough scenario that overcomes the fragmentation barrier and builds a mass reservoir in the runaway growth of large rocky particles.Furthermore, the final opacity is self-regulated at the order of a few in all cases, roughly consistent with the assumed unity opacity adopted in the initial radiative disc model.
• A successful breakthrough scenario favours a larger fragmentation velocity threshold and larger particles in the accreted materials from the external disc, which makes breakthrough growth easier.However, a high dust-to-gas ratio in the supplied materials may hinder dust retention as it takes a longer time to process small dust particles in the size distribution evolution, leading to a higher opacity and therefore a hostile higher disc temperature.Kunimoto & Matthews (2020) suggested that the occurrence rate of close-in Kepler planets generally increases from early-type to late-type around FGK stars.To explore how the dust evolution scenario proposed in this work varies around stars with different masses, Fig. 6 shows the range of the accretion rate  that meets the temperature requirements of our scenario, that is, the initial disc midplane temperature  rad,0 at the truncation radius  T is above 1000 K (such that disc is truncated) and  rad,0 at the dust accumulation radius  accu is below 2000 K (such that dust can survive from sublimation), as a function of stellar properties  ★ ,  ★ , and  ★ (radius, strength of dipole magnetic field, and mass; similar to Fig. 2 in Paper I).
Under the assumptions of a fixed stellar spin period ( ★ = 8 days) and that the corotation radius  Co converges with the truncation radius  T , we find that the viable accretion rate always2 falls between ∼ 10 −9 -10 −8  ⊙ yr −1 .However, many uncertainties remain.For example,  ★ not only evolves with stellar evolution but also exhibits a very broad distribution, particularly for low-mass stars (see Fig. 4 in Reiners et al. 2022).Moreover,  ★ decreases with  ★ , but with a large scatter within each stellar mass interval.These uncertainties make it hard to pin down the range of stellar age -inferred from  ★ based on the stellar evolution model in Baraffe et al. (2015) -for the proposed dust retention scenario to happen.
Taking  ★ = 4 days as an example, Fig. 7 presents the viable accretion rate range for low-mass stars that are rapid rotators.Compared to the slow rotators shown previously in Fig. 6, these rapid ones will need to evolve to an even later stage to set off dust accumulation as the viable range of  decreases to ∼ 10 −9.5 -10 −8.5  ⊙ yr −1 .This comparison suggests that our proposed scenario for dust retention and potentially subsequent planet formation likely favours slower rotators, which is somewhat more abundant at lower stellar masses (Reiners et al. 2022;Stefánsson et al. 2023), roughly consistent with the occurrence rate findings (Kunimoto & Matthews 2020).   .The accretion rate in the stellar radius vs. stellar magnetic field strength frame for different stellar mass (from left/top to right/bottom: 0.1 ⊙ , 0.3 ⊙ , 1.0 ⊙ , 1.4 ⊙ ), assuming that the magnetospheric truncation radius converges with the corotation radius,  ⋆ = 8 days,  out = 1.25, and  = 0.01.The top axis represents the stellar age corresponding to the stellar radius on the bottom axis, as determined by the stellar evolution model in Baraffe et al. (2015).The black solid curve and the blue dash dotted curve sandwich the area where the disc temperature  rad at  T is above 10 3 K (such that disc is truncated effectively) and  rad at  accu is below 2000 K (such that refractory dust may survive), allowing potential dust trapping and retention.

(Figure 1 .
Figure 1.The evolution of vertically integrated dust surface density distribution in the particle size-time plane for all of our full-ingredient dust evolution models, with different  (from top to bottom) and  f (from left to right).In these models,  supp,max = 1000 cm,  supp = 0.01.

Figure 2
Figure2.The evolution of (from top to bottom) dust-to-gas surface density ratio, vapour-to-gas surface density ratio, local disc temperature, and total opacity for the same models shown in Fig.1, color-coded by  and line-stylecoded by  f .

Figure 4 .Figure 5 .
Figure 4. Final dust-to-gas surface density ratios  final after 10 5 years of local dust evolution for our full-ingredient dust evolution models with  supp = 0.01.Models with  supp,max = 300 cm (left) do not accumulate considerable dust mass and produce almost identical results.Models with  supp,max = 1000 cm (right) may yield runaway growth and significant solid accumulation given ideal disc conditions.An interactive version of this plot is available at https://rixinli.com/RubbleSurveyResultsII.html

T
rad @R accu = 2000K T rad @R T = 1000K Figure6.The accretion rate in the stellar radius vs. stellar magnetic field strength frame for different stellar mass (from left/top to right/bottom: 0.1 ⊙ , 0.3 ⊙ , 1.0 ⊙ , 1.4 ⊙ ), assuming that the magnetospheric truncation radius converges with the corotation radius,  ⋆ = 8 days,  out = 1.25, and  = 0.01.The top axis represents the stellar age corresponding to the stellar radius on the bottom axis, as determined by the stellar evolution model inBaraffe et al. (2015).The black solid curve and the blue dash dotted curve sandwich the area where the disc temperature  rad at  T is above 10 3 K (such that disc is truncated effectively) and  rad at  accu is below 2000 K (such that refractory dust may survive), allowing potential dust trapping and retention.

Figure 7 .
Figure 7. Similar to the top row of Fig. 6 but with the assumption of  ⋆ = 4 days.

Table 1 .
Simulation Parameters The simulation parameters (except Σ v,0 ) listed in this table are the same as those listed in Table 1 in Paper I * of  f ,  supp,max , and  supp