Analytical solutions for the evolution of MHD wind-driven accretion discs

ABSTRACT

1 INTRODUCTION

Protoplanetary discs (PPDs), often regarded as nurseries where planets may form, have been at the forefront of astrophysical research in recent decades (e.g. Armitage 2011; Andrews 2020). As our observational capabilities have progressed, we have found details within these discs, enabling the development of increasingly sophisticated theoretical models that can be rigorously tested against observational evidence. These studies have greatly advanced our understanding of the physical conditions and processes within the birthplaces of planetary systems. One of the most intriguing and complex enigmas within this domain revolves around unravelling the dominant mechanisms governing accretion in PPDs. Over the years, various mechanisms have been suggested to explain accretion disc dynamics, including gravitational instability (e.g. Lodato & Rice 2004; Rafikov 2009), magnetorotational instability (MRI; Balbus & Hawley 1991), and hydrodynamic convection (e.g. Kley, Papaloizou & Lin 1993; Held & Latter 2018). Yet, pinpointing the dominant mechanism through observational evidence continues to be a formidable task (e.g. Simon et al. 2015; Rafikov 2017; Wang & Goodman 2017; Najita & Bergin 2018; Alexander et al. 2023; Trapman et al. 2023; Wu et al. 2023).

On the other hand, a substantial body of observational data concerning PPDs has emphasized the significant impact of magnetized and photoevaporative winds, ultimately resulting in the dispersal of these discs (e.g. Pascucci et al. 2023). Notably, theoretical investigations have illuminated the connection between turbulence generation, often driven by MRI, and the concurrent production of magnetized winds (Bai & Stone 2013). These winds not only induce mass-loss but also significantly contribute to angular momentum removal. For many decades, the dominant perspective emphasized the primary role of turbulence-generating mechanisms, such as MRI, in driving the accretion process within discs. However, recent observational data have shifted this narrative, indicating that magnetohydrodynamic (MHD) winds are pivotal in driving accretion within PPDs (e.g. Weder, Mordasini & Emsenhuber 2023).

The launch of MHD winds from accretion discs has long been a subject of interest in astrophysical research. It was initially explored by Blandford & Payne (1982), who investigated the conditions required for the emergence of such winds using MHD self-similar solutions. Subsequent research expanded on this model through two primary approaches: on one hand, they explored MHD winds using analytical models (e.g. Ferreira & Pelletier 1993; Ferreira 1997), and on the other, they utilized sophisticated numerical simulations (Bai & Stone 2013). These simulations, both with and without the inclusion of non-ideal factors such as Ohmic diffusion and the Hall effect, provided deeper insights into the complexities of MHD wind ejection.

The launch of MHD winds significantly impacts the structure of the underlying disc. However, exploring this interaction through global numerical simulations poses considerable challenges. Consequently, there is a compelling motivation to refine the classical accretion paradigm, as initially formulated by Shakura & Sunyaev (1973), to study the pivotal role of MHD winds in determining the structure and dynamics of accretion discs. These models typically incorporate equations describing angular momentum extraction by MHD winds and mass-loss rates, drawing inspiration from physical principles or the findings of numerical simulations. This approach not only facilitates a relatively straightforward analysis for investigating the structure and evolution of discs influenced by MHD winds but also offers results that can be compared against existing observational evidence of PPDs. Quantities such as disc size evolution, inner mass accretion rates, and total disc mass provide valuable insights that may yield criteria for distinguishing the dominant accretion mechanisms at play within PPDs.

Another significant research domain within the context of PPDs involves investigating the dynamics of dust particles and the mechanisms governing their growth. Understanding these processes is crucial for comprehending the intricate path to planet formation and the subsequent migration of planets within the disc. Although these complex issues have been subjects of study for many years, the potential impact of MHD winds has not been widely incorporated into previous research efforts. However, with recent advancements shedding light on the substantial role of MHD winds, there is an escalating interest in exploring these phenomena. Researchers are now employing both generalized α-viscosity models (Takahashi & Muto 2018; Arakawa, Matsumoto & Honda 2021; Taki et al. 2021; Zagaria et al. 2022) and direct numerical simulations (e.g. Kimmig, Dullemond & Kley 2020; Aoyama & Bai 2023) to probe into these problems and unveil their underlying mechanisms.

The ease of constructing generalized α-viscosity models to incorporate MHD winds presents an advantage for exploring disc evolution across a wide range of input parameters. These models offer valuable physical insights into the expected behaviour of PPDs influenced by MHD winds. While the study by Suzuki, Muto & Inutsuka (2010) introduced a 1D disc model with MHD winds within the framework of the standard approach, it should be noted that their implementation of wind mass-loss was based on MHD disc simulations. Building upon this work, Suzuki et al. (2016) further extended the model to incorporate both angular momentum removal and wind mass-loss rates. A noteworthy characteristic of their evolutionary models is that disc dispersal, driven by MHD winds, begins in the inner regions and progressively extends outward as time advances.

A similar approach was undertaken by Armitage, Simon & Martin (2013) to incorporate MHD winds into their model. They used relations for disc turbulence and angular momentum removal, based on plasma parameters derived from MHD simulations. This model was later extended by Shadmehri & Ghoreyshi (2019) to include wind mass-loss rates and an in-depth study of disc evolution was conducted. Various aspects, such as disc mass, net mass accretion rate at the inner disc edge, total wind mass-loss rate, and isochrone tracks (Lodato et al. 2017), were numerically calculated for different wind strengths. The isochrone tracks derived from their disc model, incorporating MHD winds, displayed reasonable fits with observations PPDs in star-forming regions like Lupus and σ-Orion.

A steady-state adaptation of this model was implemented by Khajenabi et al. (2018), leading to analytical solutions that describe the consistent structure of PPDs with MHD winds. Their solutions also revealed an inner region characterized by a flat density profile, a consequence of the presence of MHD winds, in alignment with previous studies (e.g. Suzuki et al. 2016). These solutions have been utilized to explain the observational features of HL Tau, particularly its high accretion rate and the thinness of its dust layer (Hasegawa et al. 2017).

In parallel, several researchers have extended these disc models with MHD winds to include the dynamics of dust particles. For example, Takahashi & Muto (2018) investigated the interplay between gas and dust in young PPDs, exploring how MHD wind launch processes contribute to the formation of ring-hole structures and the emergence of pressure bumps that migrate outward over time. Meanwhile, Taki et al. (2021) found a novel growth mode leading to the coalescence of dust grains into kilometre-sized bodies. This growth phenomenon, influenced by variations in the gas pressure profile induced by the wind, gives rise to ring-like configurations of planetesimal-sized entities within the inner regions of the disc, potentially shaping the subsequent planet formation process. Additionally, the outward radial drift of pebbles within PPDs has been attributed to the effects of MHD wind launch (Arakawa et al. 2021). Remarkably, in the context of circumplanetary discs, the emergence of MHD winds has also been shown to enhance the growth rate of dust particles through mutual collisions (Shibaike & Mori 2023).

The majority of theoretical models for MHD wind-driven accretion discs, based on the generalized accretion disc model, are numerical in nature, lacking a closed analytical framework to describe disc evolution in the presence of MHD winds. However, a recent study by Tabone et al. (2022, hereafter referred to as TB22) introduced an elegant set of similarity solutions for discs featuring MHD winds, characterized by a generalized form analogous to the well-known self-similar solutions for viscous discs originally proposed by Lynden-Bell & Pringle (1974). These solutions have subsequently been applied to investigate the dynamics of dust particles within the context of MHD winds (Zagaria et al. 2022) and to determine the theoretical size of MHD wind-driven accretion discs for comparison with observational data (Trapman et al. 2022) and evolution of the distribution in the disc mass and accretion rate plane of a disc population (Somigliana et al. 2023). In addition, Chambers (2019) also provided analytical solutions for modelling the evolution of discs influenced by MHD winds. While these solutions rely on simplifying approximations, they prove invaluable for conducting population synthesis studies related to planet formation (Alessi & Pudritz 2022).

In this paper, we introduce new analytical solutions for the evolution of an accretion disc where MHD wind serves as the primary mechanism driving accretion. In this extreme configuration, inspired by observations of PPDs, as we will demonstrate, our model is integrable, allowing us to explore disc evolution as an initial value problem with a prescribed initial surface density distribution. While TB22 studied a wind-dominated disc model, their analysis was rooted in self-similar solutions, which inherently do not facilitate the tracking of disc evolution from a specified surface density distribution without further confirmation through numerical investigations. In the following section, we outline the fundamental equations of our model. Section 3 provides analytical solutions for MHD wind-driven accretion discs, using arbitrary initial surface density profiles, and extensively explores a comprehensive study of disc properties. Our findings are subsequently discussed in Section 4.

2 GENERAL FORMULATION

Our model is based on the prior studies that investigated disc evolution with MHD winds, employing parametrized relationships for both angular momentum removal by MHD wind and wind mass-loss rate (Armitage et al. 2013; Suzuki et al. 2016). TB22 introduced slightly different formulations for these relations. In this work, we adopt their parametrizations to facilitate direct comparisons. This 1D model closely resembles the model for viscous evolution introduced by Shakura & Sunyaev 1973 and Lynden-Bell & Pringle 1974. The primary evolutionary equation is as follows:

In the above equation, Σ represents the disc surface density, while c_s denotes the sound speed. The Keplerian angular velocity is given by |$\Omega = \sqrt{GM/r^3}$|⁠, where M is the mass of the central star. The model incorporates three critical input parameters: α_SS, which pertains to disc turbulence as introduced by Shakura & Sunyaev (1973); α_DW, associated with angular momentum removal by the wind, as discussed in TB22; and the magnetic lever arm parameter λ, as described by Blandford & Payne (1982) and Ferreira (1997).

Our primary equation (1) involves two unknowns: surface density and sound speed. To close the system of equations, an energy equation is necessary. However, we pursue a simplified approach by adopting a power-law function of the radial distance for the disc temperature, or equivalently, sound speed. In this context, we assume that the sound speed is |$c_{\mathrm{ s}}^{2} \propto r^{-1}$|⁠.

The mass accretion rate resulting from disc turbulence is expressed as

The mass accretion rate driven by MHD winds is formulated as

The mass-loss rate due to MHD wind is parametrized as

In the absence of MHD wind, the master equation (1) simplifies to the standard viscous evolutionary model (Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974). The second term on the right-hand side accounts for angular momentum removal by the MHD wind, which drives additional mass accretion beyond that induced by disc turbulence. The third term parametrizes the wind mass-loss rate, with larger λ values indicating reduced mass removal, essential for sustaining accretion on to the central star. Theoretical studies suggest a constraint of λ > 3/2 (Blandford & Payne 1982; Ferreira 1997).

3 ANALYTICAL SOLUTIONS

As previously mentioned, our objective is to investigate the evolution of a disc when the contribution from disc turbulence to mass accretion is negligible, especially when compared to the influence of the MHD wind. Within the context of our model, this assumption translates to setting α_SS = 0. By introducing dimensionless variables, such as |$\hat{\Sigma } = \Sigma /\Sigma _\mathrm{ c}$|⁠, |$\hat{t} = t/t_\mathrm{ c}$|⁠, and |$\hat{r} = r/r_\mathrm{ c}$|⁠, we can rewrite the master equation (1) as follows:

We highlight that Σ_c and r_c represent the characteristic surface density and radial distance, respectively. For the characteristic time-scale t_c, consistent with TB22, we adopt the accretion time-scale defined as

In this expression, c_s,c and H_c denote the sound speed and disc thickness associated with the characteristic radial distance r_c, respectively. We observe that by neglecting disc turbulence, the primary equation (1) undergoes a transformation. Specifically, it changes from a second-order partial differential equation to a first-order one, thereby simplifying its analysis.

To solve equation (5), we adopt a change of variables strategy. Upon applying the subsequent change of variables,

equation (5) simplifies, facilitating a more straightforward analysis. We then obtain

This is a first-order linear partial differential equation, commonly known as the transport or advection equation in one spatial dimension. It describes the evolution of a particular quantity, denoted by |$F(\hat{r}, \hat{t})$|⁠, which moves with a constant velocity of 1/2 in the radial direction. By employing the method of characteristics, this equation can be analytically solved, yielding the general solution:

where |$F(\hat{r} + \frac{\hat{t}}{2})$| can be any arbitrary function. However, it can be determined uniquely once the initial distribution of F is specified. Thus, the disc surface density becomes

We highlight that the fully analytical solution in equation (10) is derived without imposing specific constraints on the density profile, as is typically done in similarity methods. Once the initial surface density distribution is established, we can use solution (10) to determine F and, subsequently, deduce the disc evolution.

As an example, consider employing an initial surface density given by

where n and m are free parameters. This particular choice is inspired by the self-similar viscous evolution of a disc presented by Lynden-Bell & Pringle (1974). In our notation, their solution corresponds to n = m = 1 with Σ_c and r_c as functions of time. TB22 introduced a profile analogous to that of equation (11) for disc evolution, subsequently aiming to derive the time dependencies Σ_c(t) and r_c(t). In the scenario of extreme wind dominance, their similarity solution indicated that n = (2λ − 3)/[2(λ − 1)]. It is crucial to recognize that the self-similar approach is not inherently an initial value problem and, as such, does not incorporate the initial state adequately. Consequently, the initial surface density distribution in TB22 is contingent on the wind parameter λ.

However, the potential roles of the MHD wind are examined once the initial configuration is established, independent of the wind parameters α_DW and λ. Subsequently, we adjust these parameters to probe their influence on disc evolution. It is imperative to note that the initial state of the disc should remain uninfluenced by these wind parameters, ensuring that any alterations to them leave the initial configuration unchanged. Consequently, equation (11) presents the proposed initial surface density that aligns with our theoretical expectations and is suitable for our intended purpose. Using equations (10) and (11), we obtain

This solution characterizes the evolution of a PPD in the wind-dominated paradigm corresponding to any given exponent n. Specifically, for n = (2λ − 3)/[2(λ − 1)], our solution becomes the self-similar solution deduced by TB22 [refer to equation (27) in TB22]. This distinction implies that in the TB22 model, the initial configuration depends on the magnetic lever arm λ, while in our approach, the exponent n for the initial state is a pre-determined input. In our model, the primary input parameters encompass Σ_c, r_c, λ, t_acc,0, n, and m. The characteristic radius r_c typically ranges between 20–50 astronomical units (au).

The surface density Σ_c can be deduced from the initial disc mass by utilizing the surface density profile as presented in equation (11). This relationship can be written as

where |$I(\hat{r}_{\mathrm{in}},n,m)$| is given by

and |$\hat{r}_{\mathrm{in}} = r_{\mathrm{in}}/r_\mathrm{ c}$| and |$\hat{r}_{\mathrm{out}} = r_{\mathrm{out}}/r_\mathrm{ c}$|⁠. Notably, r_in represents the inner disc radius, and it is typically much smaller than r_c. As for the outer edge of the disc, we generally assume a large value, say, r_out = 1000 au. Furthermore, the term M_D0 represents the initial disc mass. For typical PPDs, this mass usually falls within the range 10⁻³ to 10⁻² M_⊙.

Utilizing the solution presented in equation (12), we initiate a comprehensive investigation into the characteristic evolutionary dynamics embedded within our disc model. Fig. 1 offers a view of the surface density’s evolution, depicting its progression over a range of times, which are expressed in units of t_acc,0. The initial disc mass is set at 10⁻² M_⊙ and the characteristic radius r_c is chosen to be 50 au. In addition, we set the remaining input parameters to n = m = 1 and λ = 2. In Fig. 1, the evolutionary curves for t = 0, 4t_acc,0, 10t_acc,0, and 16t_acc,0 are displayed in descending order from top to bottom.

Figure 1.

Evolution of the surface density as a function of radial distance based on our analytical solution for the wind-dominated disc with input parameters r_c = 50 au, M_D0 = 0.01 M_⊙, m = n = 1, and λ = 2. The curves, from top to bottom, correspond to time instances: t = 0, 4t_acc,0, 10t_acc,0, and 16t_acc,0, respectively.

Open in new tabDownload slide

A key feature in Fig. 1 is the consistent decrease of the surface density across all radii as time progresses. Over a span of 16t_acc,0, the reduction in surface density is significant, approximating a factor of 10 000. The reduction of the disc surface density is attributed to the dual effects of wind-driven accretion on to the star and the wind-induced mass-loss. Initially, the density distribution manifests as particularly steep, especially in the disc’s outer regions. However, as time elapses, the surface density profile adopts a more moderate incline. This moderation initiates from the innermost regions and systematically extends outwards. This inside-out manner of gas removal remains largely unaltered across varying model parameters. While the early evolution of the disc is influenced by the initially chosen configuration, over time, the disc tends toward a self-similar profile in its later stages. The characteristic time-scale t_acc,0 provides a measure of the disc’s time-scale evolution. For periods longer than t_acc,0, the disc evolution exhibits self-similar behaviour.

In the context of MHD wind mass-loss, the magnetic lever arm parameter λ serves as a measure of efficiency. An increase in λ results in a reduced rate of wind mass-loss to maintain the needed mass accretion rate on to the star. Fig. 2 displays the influence of the wind parameter λ on the evolution of the surface density. The considered input parameters align with those in Fig. 1. However, in this depiction, the surface density profile is consistently captured at 10t_acc,0, spanning varying λ values – ascending from 2 at the bottom to 7 at the top. A clear pattern is evident: as λ increases, the decline in surface density becomes less prominent, consistent with our anticipations. While the outer regions of the disc largely resist alterations from this parameter, the inner domains display marked sensitivity. However, it is worth noting that as λ exceeds an approximate value of 5, its subsequent impact on the surface density becomes considerably reduced.

Figure 2.

This figure illustrates the initial and subsequent surface density profiles across radial distances, emphasizing the role of the wind parameter λ. All the input parameters remain consistent with those specified in Fig. 1. While the profile at t = 0 remains unchanged, for t = 10t_acc,0, the surface density profiles are presented for different λ values, ranging from λ = 2 at the bottom to λ = 7 at the top.

Open in new tabDownload slide

We have also examined scenarios with varying values of n and m, which correspond to different initial surface density profiles. The observed trends in these cases align closely with the results presented thus far. While we have chosen not to detail these specific cases here, the influence of the exponent n on disc quantities, including total disc mass and the mass accretion rate, will be discussed in the subsequent figures.

Given that our model incorporates MHD wind as the dominant mechanism driving both mass accretion on to the host star and mass-loss, it is expected that the total disc mass will decrease over time. To quantify this theoretical expectation, we employ the surface density profile outlined in equation (12) to compute the total disc mass. The total disc mass is derived as follows:

where |$I(\hat{r}_{\rm in}, n, m)$| is given by equation (13) and the function |$f(\hat{r}, \hat{t})$| is

Fig. 3 presents the evolution of the total disc mass normalized by its initial value over time. The ratio M_D(t)/M_D0 is shown for two distinct initial configurations, characterized by n = 1 and n = 2. For each value of n, we also explore various values of λ as shown in the figure. The disc’s inner edge is consistently set at r_in = 0.05 au.

Figure 3.

Evolution of the normalized total disc mass over time for varying values of λ = 2, 4, 12 and exponents n = 1 and 2. The total disc mass is normalized by its initial value. The input parameters, unless otherwise mentioned, remain consistent with those presented in Fig. 1. Regardless of the chosen input parameters, the total disc mass exhibits a decline over time.

Open in new tabDownload slide

A general trend evident across all examined cases is the anticipated decline in total disc mass with time. While this decline is initially moderate, an approximate distinct point can be observed beyond which the reduction becomes sharply accelerated. This suggests a two-phase evolutionary sequence for the disc: the early phase is marked by a measured depletion of mass, whereas the subsequent phase witnesses a precipitous drop within a limited time span. Notably, similar trends have been highlighted in prior studies (e.g. Suzuki et al. 2016; Shadmehri & Ghoreyshi 2019; Zagaria et al. 2022)

A decrease in λ further amplifies the reduction in disc mass. Yet, a particular point of intrigue emerges when considering n = 2, representing a considerably steeper initial configuration. In this case, the gradual phase is shortened significantly, leading to an earlier onset of rapid decline. Given that larger n values, assuming constant initial mass, suggest smaller disc sizes, it can be inferred that steeper initial surface density profiles accelerate the disc’s depletion rate. Consequently, the lifespan of the disc is substantially reduced for steeper profiles. These observed trends diverge from the findings of TB22. While they also identified an exponential decrease in total mass, their results were not influenced by the magnetic lever arm or the slope of the initial surface density.

In Fig. 4, we explore the disc’s lifetime t_life across a range of values for both λ and n. In our analysis, we define the disc’s lifetime as the moment when its mass diminishes to a pre-defined fraction of its initial mass. It is worth noting that while we adopt a threshold fraction of 10⁻³ in this study, the general trends and behaviour of the disc’s lifetime, in relation to the input parameters, remain largely consistent irrespective of the specific choice of this threshold. Fig. 4 shows that the disc lifetime increases as the magnetic lever arm increases, but beyond λ ≃ 6, this lifetime remains nearly uniform no matter what exponent n is adopted.

Figure 4.

In this figure, we present the normalized disc lifetime t_life, scaled by t_acc,0, plotted against the wind parameter λ for various values of n as indicated. The disc lifetime is defined at the point where the disc mass falls below a specified threshold, in this instance, 10⁻³ of its initial mass.

Open in new tabDownload slide

Fig. 5 offers a comparative analysis of the total disc mass for various values of λ = 4, 8, 12, 16 relative to the baseline case of λ = 2. This comparison aims to underscore the pivotal role played by the magnetic lever arm parameter, λ, in influencing the effectiveness of mass removal via wind-driven processes. Contrasting the evolution trajectories, it becomes evident that the time-dependent ratio M_D(t, λ)/M_D(t, λ = 2) for n = 1 increases with time. Furthermore, discs associated with lower λ values experience a more pronounced wind-driven mass extraction compared to their counterparts with higher λ values. To put this into perspective, by the epoch of 20t_acc,0, the disc mass corresponding to λ = 4 is approximately double that of the λ = 2 case. Simultaneously, the disc mass associated with λ = 8 is found to be three times that of the λ = 2 case. This comparison effectively illustrates the intricate interaction between λ and the dynamics of wind-driven mass depletion in PPDs.

Figure 5.

Evolution of the mass ratio M_D(t, λ)/M_D(t, λ = 2) for n = 1. As λ increases, a larger fraction of the mass remains in the disc. Over time, this difference in retained mass becomes more pronounced.

Open in new tabDownload slide

In the framework of our model, mass reduction of the disc is primarily due to the accretion on to the central star and mass ejection through MHD winds. Now, we precisely quantify these mass-loss rates. Using equations (2) and (3), we can now determine the net mass accretion rate on to the star, denoted by |$\dot{M}_{\star }$|⁠, which is expressed as

It is important to note that both terms are evaluated at the inner disc edge. In the wind-dominated regime, however, the contribution of the first term is zero. Therefore, we obtain the following expression for |$\dot{M}_{\star }$|⁠:

Using the prescribed local mass-loss rate, as given by equation (4), we can calculate the cumulative mass-loss rate across the entirety of the disc. This rate, denoted by |$\dot{M}_{\rm W}$|⁠, is formulated through the integral of the local rate over the radial extent of the disc, spanning from its inner edge, r_in, to its outer boundary, r_out. This cumulative rate is given by the relation

Thus, we obtain

A particularly suitable parameter for assessing the influence of MHD winds within the disc is the dimensionless mass ejection-to-accretion ratio, denoted as |$\dot{M}_{\rm W}/\dot{M}_{\star }$| (e.g. Ferreira 1997). One of the notable advantages of this parameter is its observability, in addition to its independence from the disc mass. TB22 derived an expression for this ratio utilizing their self-similar solution. In a case where the wind dominates, they found this ratio as (r_c/r_in)^{1/[2(λ − 1)]} − 1. But our solution implies that

The ratio is dependent upon the input parameters and, as we will demonstrate, can vary over time.

The top plot in Fig. 6 illustrates the mass ejection-to-accretion ratio |$\dot{M}_{\rm W} (t)/ \dot{M}_{\star } (t)$| as a function of time, considering various values of λ = 2, 4, 8, 12, and 16 from the top curve to the bottom curve, respectively. It is evident that for λ values greater than 4, the ratio remains below unity throughout the entire evolution of the PPD. Conversely, for λ values less than 4, the ratio exceeds unity. In the latter case, the ratio demonstrates a non-uniform behaviour over time, which contrasts the expectations based on the self-similar solution presented by TB22. The increase in the mass ejection-to-accretion ratio over time indicates that the rate of wind mass-loss becomes greater than the mass accretion rate on to the central star.

Figure 6.

Evolution of the mass ejection-to-accretion ratio with time is shown for different input parameters. The top plot displays the ratio for various magnetic lever arm parameter values, λ = 2, 4, 8, 12, 16, highlighting the influence of λ on mass depletion dynamics. The bottom plot contrasts this by focusing on a fixed λ = 2 and varying the exponent n = 0.2, 0.5, 1, 1.4, 2, emphasizing the role of the initial disc size in the ratio’s temporal behaviour. The above trends deviate from the uniform behaviour predicted by the self-similar solution presented by TB22.

Open in new tabDownload slide

The bottom plot of Fig. 6 exhibits the evolution of the mass ejection-to-accretion ratio for a fixed λ = 2, while considering different values of the exponent n = 0.2, 0.5, 1, 1.4, 2, from the top curve to the bottom curve, respectively. Remarkably, while the ratio remains relatively uniform over time for n = 0.5, for other n values, the ratio displays substantial temporal variation. This distinctive trend and evolution of the mass ejection-to-accretion ratio are not captured by the self-similar solution presented by TB22, where a uniform ratio was obtained for all input parameters. In contrast, Fig. 6 (bottom) clearly illustrates that the ratio increases with time when the exponent n exceeds 0.5. However, for n < 0.5, the trend is reversed, and the ratio decreases with time. Our model suggests that the initial size of the disc plays a pivotal role in determining this ratio’s behaviour. In smaller discs, the ratio gradually increases over time, while in larger discs, the ratio declines.

In Fig. 7, we explore the evolution of PPDs, not only focusing on changes in the total disc mass but also variations in the mass accretion rate at the inner disc edge. Similar to previous studies, to characterize this evolution, we introduce a plane with the vertical axis represented by |$\dot{M}_{\star }(t)/\dot{M}(t=0)$| and the horizontal axis by M_D(t)/M_D(t = 0). Each snapshot in time corresponds to a distinct point on a plane |$\dot{M}_{\star }{\!-\!}M_{\rm D}$|⁠, together forming a trajectory that chronicles the disc’s progression. Although our surface density profile is analytical, the disc’s trajectory in this plane is determined numerically using equations (15) and (18).

Figure 7.

Evolutionary trajectories of PPDs in the |$\dot{M}_{\star }{\!-\!}M_{\rm D}$| plane. The top plot displays trajectories for constant n = 1 with varying magnetic lever arm parameter λ, showing diverse characteristics based on λ values. In contrast, the bottom plot, set at a constant λ = 2, explores the impact of differing exponent n values on these trajectories, highlighting the influence of the initial surface density distribution slope. As evolution unfolds, trajectory slopes from both plots exhibit a tendency for convergence.

Open in new tabDownload slide

In the top plot of Fig. 7, the trajectories of disc evolution are shown for a constant n = 1 across varying magnetic lever arm parameter λ values, specifically λ = 2, 4, 8, 12, 16, sequenced from bottom to top. Notably, TB22, in their analysis, employed their self-similar solutions to deduce these trajectories, uncovering a linear relation in wind-dominated case, i.e. |$\dot{M}_{\star } \propto M_{\rm D}$|⁠, independent of the λ parameter. However, our findings in Fig. 7 underscore that trajectories can exhibit diverse characteristics depending on the value of λ. In scenarios where wind mass removal is markedly efficient (lower λ values), the trajectory’s slope, especially in the nascent stages of disc evolution, is steeper compared to discs with diminished mass removal rates (higher λ values). Yet, as the timeline progresses, a convergence in slope is evident across the explored cases. For more substantial λ values, there is a discernible downward shift in the trajectory, suggesting a reduced inner accretion rate for equivalent total disc mass.

The bottom plot of Fig. 7 examines cases with a constant λ = 2, exploring variations in the exponent n with values n = 0.4, 1, 1.4, 2. Our analysis reveals that the trajectory within the |$\dot{M}_{\star }{\!-\!}M_{\rm D}$| plane, as the disc evolves temporally, is considerably influenced by the initial slope of the surface density distribution. As noted earlier, a rise in n for a specified total mass implies a reduced disc size. This is particularly evident when size is determined by a threshold surface density metric (e.g. Toci et al. 2023) or by the radius encompassing a certain fraction of the disc’s mass (e.g. Anderson, Adams & Calvet 2013). As a result, smaller initial disc sizes lead to steeper trajectories with a downward shift. Conversely, larger initial disc sizes produce an upward shift in the trajectory, aligning more closely with a power-law relationship. Interestingly, during the later stages of evolution, the trajectory slopes converge to the same value, independent of the adopted values for the exponent n.

We have shown corresponding curves in Fig. 7 using dimensionless physical quantities, providing a convenient framework for a more straightforward interpretation of theoretical trends. However, to facilitate a meaningful comparison with observational data, it is necessary to illustrate these trends in the |$\dot{M}_{\star }{\!-\!}M_{\mathrm{ D}}$| plane using actual physical units. It is essential to acknowledge that a comprehensive understanding of the conditions required to align these curves with observations necessitates a detailed exploration through realistic disc population synthesis, a task that lies beyond the scope of the present study.

Fig. 8 is similar to Fig. 7, with the distinction that we now present the curves in terms of actual physical units. The initial disc mass is M_D0 = 0.01 M_⊙. As previously discussed in TB22, typical values for the accretion time-scale, t_acc,0, are of the order of a million years (Myr). In this context, we adopt t_acc,0 = 1 Myr.

In both Figs 7 and 8, we can easily recognize an initial phase of disc evolution followed by a profile corresponding to self-similar evolution. The duration of the early disc phase is dependent on the wind parameter λ and the initial surface density exponent n. As demonstrated in the total disc mass evolution presented in Fig. 3, we have previously established that disc evolution commences with an early phase, followed by a more extended, long-term evolutionary phase. During the early disc phase evolution, the total disc mass undergoes only a modest reduction. However, beyond a certain time, the disc mass experiences a rapid decline.

Figure 8.

This figure is identical to Fig. 7, with the distinction that the presented curves represent the evolutionary trajectories in the |$\dot{M}_{\star }{\!-\!}M_{\mathrm{ D}}$| plane using actual physical units. These curves correspond to an ensemble of discs characterized by an initial disc mass of M_D0 = 0.01 M_⊙ and an accretion time-scale of t_acc,0 = 1 Myr.

Open in new tabDownload slide

Fig. 3 illustrates that the duration of the initial phase is of the order of the initial accretion time-scale t_acc,0. In this initial disc phase evolution, the total disc mass does not decrease rapidly, as depicted in both Figs 7 and 8 by the initial segments that deviate from the self-similar profile. In the case with n = 1 and λ = 2, the duration of the initial phase is slightly longer, though it remains shorter than the initial accretion time-scale.

For a fixed λ = 2 and different values of n, the longest duration of the initial phase belongs to n = 2. Notably, during the early phase in this specific case, the total disc mass is reduced by a factor of 10². Fig. 3 indicates that it takes a few initial time-scales to achieve this level of mass reduction. Therefore, it appears that the initial evolutionary phase is significantly shorter than the entire lifetime of the disc.

4 DISCUSSION AND CONCLUSIONS

Our analytical solutions for the evolution of MHD wind-dominated discs demonstrate that the initial disc configuration plays a pivotal role in subsequent disc evolution. When disc turbulence is parametrized by a viscous term, the disc evolves towards a state independent of the initially adopted surface density profile. In this scenario, it is as if the disc forgets its initial state due to the efficient spreading caused by viscous processes. However, in the absence of disc turbulence, the nature of the master equation changes, as we have already discussed. In such a case, the disc retains its initially chosen configuration in t < t_acc,0, but subsequently the solution converges to the self-similar one. In our study, we have utilized a commonly adopted initial surface density profile, and our results differ from those presented by TB22. It is important to note that the exact evolutionary sequence of PPDs, including the onset of MHD wind launching, remains a subject of uncertainty.

During the early stages of PPD formation, these discs tend to be more massive and susceptible to gravitational instability as the dominant mechanism for angular momentum transfer (e.g. Kratter & Lodato 2016). It remains uncertain whether, as the disc depletes its mass, it undergoes a transition into an MRI-driven accretion phase, with or without the influence of MHD winds. For a specific choice of n = 1 in our model, we have assumed that the disc starts its evolution driven by viscous mechanisms, ultimately adopting a self-similar profile similar to that described by Lynden-Bell & Pringle (1974) when the wind-dominated regime commences. It is at this stage that MHD wind becomes efficient, justifying our choice of initial configuration with n = 1. It is important to recognize that this scenario regarding the adopted initial state is speculative, and further extensive studies are warranted to explore the diverse phases of PPD evolution.

Dispersal of PPDs, whether due to MHD winds or photoevaporative winds, remains a critical area of research (e.g. Kunitomo, Suzuki & Inutsuka 2020; Pascucci et al. 2023). Kunitomo et al. (2020) investigated the impact of both MHD and photoevaporative winds by incorporating mass-loss rates associated with the winds and considering angular momentum extraction by MHD winds, employing a generalized standard disc model. Their findings reveal that, during the early stages of disc evolution, MHD winds are efficient in the inner disc region. However, as the evolution progresses, photoevaporative winds become dominant, particularly in the outer disc region and during the later phases of disc evolution. In our analytical model, we explored an extreme scenario where MHD winds are the dominant factor in disc dispersion. Within this framework, disc lifetime is quantified as the period during which the disc loses a substantial portion of its initial mass. Our findings suggest that the total disc mass decreases over time, with its evolution being intricately tied to the magnetic lever arm and the gradient of the initial density profile. This contrasts with the findings of TB22, who observed an exponential decline regardless of wind parameters and initial conditions based on their similarity solutions. In our wind-dominated solutions, smaller discs exhibit a reduced lifetime. It is important to emphasize that this characteristic trend emerges in the absence of disc turbulence, which typically promotes disc spreading.

As mentioned earlier, discriminating between the dominant accretion processes in PPDs, whether driven by viscous evolution or angular momentum removal due to MHD winds, is currently a subject of intense research (e.g. Long et al. 2022; Trapman et al. 2023). Various diagnostic approaches have been proposed (e.g. Tabone et al. 2022; Alexander et al. 2023; Somigliana et al. 2023), including the examination of disc size evolution. Typically, viscous evolutionary models tend to predict disc spreading, whereas, in wind-dominated scenarios, disc sizes often remain relatively constant or may even shrink due to mass depletion. However, a central challenge lies in precisely defining disc size for meaningful comparisons with observed data (e.g. Trapman et al. 2020; Toci et al. 2021, 2023). Recent theoretical models have integrated MHD winds and determined that they can effectively explain both the dispersal and accretion characteristics of discs (Tabone et al. 2022; Trapman et al. 2022; Zagaria et al. 2022). Our analytical solution reveals that, in addition to factors like initial disc mass, accretion time-scale, and magnetic lever arm parameter, the initial disc size can profoundly influence disc evolution in an MHD wind-dominated regime. While investigating disc size evolution using our solutions to compare with existing observed PPDs, or undertaking population disc size synthesis, is beyond the scope of this study, our analytical solutions present a computational advantage for such studies.

When MHD winds are highly efficient, and the PPDs are substantial in size, the mechanisms governing planet formation within PPDs have more time to operate. Planet formation under these conditions presents an intriguing avenue for exploration, warranting the development of more sophisticated models. In addition to disc lifetime, the surface density gradually tends to adopt flatter profiles, influencing the radial drift of dust particles and even planet migration, which depends on the surface density slope. This phenomenon leads to the emergence of local gas pressure maxima due to wind launching, capable of trapping dust particles and enhancing their growth rates (Arakawa et al. 2021; Taki et al. 2021; Zagaria et al. 2022).

Our analytical solutions serve as a valuable tool for efficiently investigating the influence of wind strength and initial disc size on grain growth within the inner disc regions. Furthermore, previous research has shown that the emergence of MHD winds can inhibit Type I planet migration in the close-in region (Ogihara, Morbidelli & Guillot 2015; Ogihara et al. 2018) or even induce outward planet migration (Kimmig et al. 2020). Although our current study does not specifically address planet migration, our solutions retain relevance for conducting investigations in this domain.

In the context of MHD wind-dominated PPDs, our investigation yielded several key findings:

• The lifespan of these discs demonstrates a positive correlation with the magnetic lever arm, provided this parameter remains below approximately 6. Importantly, this correlation appears to be independent of the initial disc size.

• For PPDs with larger values of the magnetic lever arm λ, the total disc mass remains notably higher than that of discs with smaller λ values, even after a substantial duration of around 20 characteristic time-scales, equivalent to 20 Myr, given our adoption of a characteristic time-scale of about 1 Myr.

• The ratio of mass ejection to accretion tends to increase over time in scenarios characterized by highly efficient winds, with this ratio surpassing unity. Notably, as the initial surface density profile becomes steeper, this ratio experiences a significant enhancement compared to discs that initiate their evolution with a flatter density profile.

• Furthermore, when examining the plane defined by the accretion rate at the inner disc edge and the total disc mass, we observed diverse trajectories contingent upon the values of the magnetic lever arm and the initial disc size.

These findings offer valuable insights into the intricate dynamics governing MHD wind-dominated PPDs, shedding light on the role of various parameters in shaping their evolutionary trajectories.

ACKNOWLEDGEMENTS

We express our gratitude to the referee for providing valuable insights through their insightful report, which greatly contributed to the improvement of this paper. All figures were generated with the python-based package matplotlib (Hunter 2007).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

References