© The Authors 2025

1 Introduction

Star formation occurs in dense cores within molecular clouds. Low-mass protostars grow in mass by accreting gas from the collapsing molecular cloud. Since the contracting cloud contains a substantial amount of angular momentum, the infalling gas forms an accretion disk around the protostar (Shu et al. 1987). If a reliable viscosity-generating mechanism exists to transport angular momentum outward in the disk, the disk’s material spirals inward, while gravitational energy is being transformed to radiation and internal energy. Mechanisms that can provide a viscosity in the disk are the gravitational instability (GI) (e.g., Lodato & Rice 2004) and the magneto-rotational instability (MRI) (Balbus & Hawley 1998).

Despite this broad picture, a detailed physical understanding of the accretion process during the proto-stellar stage, when the forming star is still embedded within and accreting from its parent dense core, remains elusive. In the simplest model of star formation, the collapse of a singular isothermal sphere, mass accretes onto the protostar at a constant rate of ~2 × 10⁻⁶ M_⊙ yr⁻¹ (Shu et al. 1987). However, as Kenyon et al. (1990) first discussed, most protostars have luminosities that are significantly lower than the accretion luminosity expected from accretion at this rate. This luminosity problem has been highlighted by Spitzer Space Telescope (Werner et al. 2004) observations of nearby star-forming regions, which reveal a significant population of protostars with luminosities lower than these theoretical predictions (e.g., Evans et al. 2009; Dunham et al. 2015).

There is growing evidence that episodic rather than continuous mass accretion is the solution to the luminosity problem, with long periods of very low accretion punctuated by short bursts of rapid accretion (resulting in luminosity outbursts) (Audard et al. 2014; Herczeg et al. 2017; Mercer & Stamatellos 2020; Fischer et al. 2023). The luminosities during the periods of low accretion predicted by simulations (Dunham & Vorobyov 2012; Stamatellos et al. 2012; Vorobyov & Basu 2015) do indeed match the observed proto-stellar luminosity distribution and can alleviate the luminosity problem.

An open question in star formation is whether the formation and evolutionary stages of substellar mass, that is, brown dwarfs (BDs) (0.013 M_⊙ ≤ M ≤ 0.08 M_⊙) and planetary-mass objects (≤0.013 M_⊙) are comparable to those of low-mass stars. In order to unravel this mystery, a combination of theoretical models and observations of young BDs in their most embedded stages of their formation (comparable to the Class 0/I stages of low-mass star formation, here called proto-BDs) is essential. Despite the fact that numerical simulation studies have provided insight into the various brown dwarf formation processes, for example, whether it is a star-like formation via gravitational core collapse (Machida et al. 2009) or alternative mechanisms of formation via disk fragmentation or as ejected embryos (Goodwin & Whitworth 2007; Stamatellos & Whitworth 2009; Bate 2012; Mercer & Stamatellos 2017; Mercer & Stamatellos 2020), the number of observed proto-BDs candidates remains small.

Observations of Class 0/I brown dwarf candidates (Palau et al. 2014; Morata et al. 2015; Riaz et al. 2015; Riaz et al. 2024) suggest that the formation of these systems can be explained by a scaled-down version of star formation. In more recent years, telescopes such as Atacama Large Millimeter/submillimeter Array (ALMA) have allowed us to spatially resolve the inner few au scales of these very low luminosity objects (VeLLOs) and constrain theoretical models (Lee et al. 2018). Observations have revealed features, such as infalling envelopes, rotationally supported pseudo-disks, and outflow-driving areas (Riaz et al. 2017, 2019; Riaz & Machida 2021), that are similar to those found in proto-stellar systems. The bolometric luminosity of these objects is generally ≤0.1 L_⊙, while mean accretion rates range from ~10⁻⁶ to ~10⁻⁸ M_⊙ yr⁻¹ (Riaz & Bally 2021).

Observations of the formation of proto-BDs are quite challenging. The physical scales of early-stage Class 0/I brown dwarfs are estimated to be at least 10 times smaller than those of low-mass protostars (Machida et al. 2009; Riaz et al. 2019). Furthermore, the main accretion phase, that is, the duration of the Class 0 and Class I phases combined, significantly differs between brown dwarfs and low-mass stars. The envelope dissipation timescale is predicted to be quite brief, ≲10⁴ yr, in the formation of brown dwarfs, according to simulations performed by Machida et al. (2009); Machida & Matsumoto (2012). In contrast, protostars with ages of less than 0.1 Myr have been shown observable signs of infall from an envelope (Takakuwa et al. 2012; Davidson et al. 2011), while Dunham et al. (2015) estimate lifetimes of 0.13–0.26 Myr for the Class 0 phase and 0.27–0.52 Myr for the Class I phase.

In most previous theoretical studies, a continuous gas accretion is assumed onto the proto-BD. However, in the inner few au, the gas can be hot enough to be thermally ionized and couple with magnetic fields activating the MRI (e.g., Zhu et al. 2009, 2010). As a result, the accretion can be episodic, switching between periods of enhanced and quiescent accretion when the MRI is active or not, respectively, just as it happens on solar-mass stars. Episodic accretion may also be due to infalling fragments that formed due to GI (e.g., Vorobyov & Basu 2005, 2010, 2015), or due to thermal instability in the inner disk region (Bell et al. 1995).

Stamatellos et al. (2011); Stamatellos et al. (2012) examine the role of episodic accretion on disk fragmentation around young solar-mass stars, regulating the number of low-mass stars and BDs forming (see also Mercer & Stamatellos 2017). Lomax et al. (2014, 2015) show that episodic accretion is necessary for reproducing the initial mass function (IMF) of stellar clusters, as continuous accretion leads to a lack of low-mass stars and BDs. MacFarlane & Stamatellos (2017) and Vorobyov et al. (2020) study the effect of episodic outbursts on the structures of pro- tostellar disks. MacFarlane et al. (2019a,b) discuss the expected observational signatures of episodic accretion/outbursts around young stars.

The role of episodic accretion during the formation of proto-BD remains mostly unexplored. Observations of multiple shocked emission knots at the extended jet (HH1165) of the very low-luminosity Class I candidate proto-BD Mayrit 1701117 (Riaz et al. 2017, 2019) indicate that episodic accretion events may also happen during brown dwarf formation.

In this paper, we use the Smoothed Particle Hydrodynamics (SPH) code SEREN (Hubber et al. 2011) to simulate the effect of episodic accretion/episodic outbursts during brown dwarf formation via the gravitational collapse of a very low-mass pre-BD core. We probe the evolution of the proto-BD until the end of the main accretion phase (Machida et al. 2009), examining the role of the timescale of outbursts and the interval between them. We discuss the observational implications of our models and compare our results with the proto-brown dwarf observations of Riaz & Machida (2021). Additionally, we compare the results of simulations of BD formation with those of low-mass star formation.

2 Computational method

We use the Smoothed Particle Hydrodynamics code SEREN (Hubber et al. 2011) to treat the gas hydro- and thermodynamics. The code uses a 2nd-order Runge-Kutta integration approach, multiple particle time-steps, and an octal tree to compute gravitational forces and find particle neighbors. In order to decrease artificial shear viscosity, the code employs timedependent artificial viscosity (Morris & Monaghan 1997) with parameters α_min = 0.1, β = 2α, and a Balsara switch (Balsara 1995).

The Stamatellos et al. (2007) technique is used to treat the chemical and radiative processes that control gas temperature (see also Forgan et al. 2009; Mercer et al. 2018). The technique accounts for radiative cooling/heating, viscous heating, background radiation field heating, and compressional heating. The equation of state takes into account the rotational and vibrational degrees of freedom of molecular hydrogen, as well as the various chemical states of hydrogen and helium (Boley et al. 2007). It is assumed that the gas is a mixture of 70% hydrogen and 30% helium. The dust and gas opacity is set to κ(ρ, T) = κ_0ρ^αT^β, where κ₀, α, β are constants (Bell et al. 1995) that depend on the physical processes causing the opacity, such as dust sublimation, ice mantle melting, molecular contributions, and H⁻ contributions.

The self-gravitating gas dynamics, energy equation, and related radiation transport are treated explicitly up to densities of ρ = 10⁻⁹ g cm⁻³ . It is assumed that a gravitationally confined condensation with ρ > 10⁻⁹ g cm⁻³ will eventually develop into a bound object (protostar or proto-BD). At this point, a sink (Bate et al. 1995), with a radius of 1 au, replaces the condensation to avoid very short time-steps. Only the sink’s gravity and luminosity affects the rest of the computational domain. Gas, that flows within the sink radius and is bound to it, is added to the sink’s mass. A forming proto-BD will provide radiative feedback to the surrounding environment, setting the pseudo-background radiation field, that is, a minimum temperature floor under which the gas cannot cool radiatively (for details see Stamatellos et al. 2007; Mercer & Stamatellos 2017) that is set to $$T_{\text{BGR}}^4(\mathit{r})={(10K)}^4+\frac{L}{16\pi {\sigma }_{SB}{\left|\mathit{r}-{\mathit{r}}_{BD}\right|}^2},$$(1)

where L and r_sink are the luminosity and the position of the proto-BD, respectively. The luminosity at the early stages that we model is due to gas accretion onto the object, and is given by $$L=f\frac{G{M}_{\text{BD}}{\dot{M}}_{BD}}{R_{\text{BD}}},$$(2)

where M_BD is the mass of the proto-BD, R_BD its radius and Ṁ_BD the accretion rate onto it. f = 0.75 is the fraction of the accretion energy that is radiated away at the photosphere of the proto-BD, rather than being expended, for instance, in driving jets and/or winds (Shang et al. 2007; Offner et al. 2010; Stamatellos et al. 2012). We assume R_BD = 3 R_⊙ which is the typical radius of a young forming object (Palla & Stahler 1993; Bate 2011). The background temperature profile is variable (see Eqs. (1) and (2)). It depends on the accretion rate onto the proto-BD and on the distance from it, T_BGR ∝ r⁻², which is the anticipated radial temperature profile on an envelope around a luminosity source.

3 Model setup

3.1 A semi-analytical model of episodic accretion

We use the model of Stamatellos et al. (2011) to include the effect of episodic accretion (resulting in episodic radiative feedback, i.e., in luminosity outbursts) on a forming proto-BD. The model uses a sub-grid model to describe the region within a sink SPH particle. Here, we briefly describe the basic aspects of the model that are needed to interpret the results of our simulations. More details of the model are given in Stamatellos et al. (2011); Stamatellos et al. (2012). We assume that the MRI develops in the inner region (<1 AU) around BDs the same way it does around solar-mass stars. It is therefore assumed that the density and temperature conditions very close to the BDs are appropriate for this to happen (Bae et al. 2014).

The model assumes that the sink consists of the proto- BD (M_BD) and an inner accretion disk (IAD; M_IAD) which is notionally the part of the disk that is not resolved by the SPH simulation. Gas that is accreted onto the sink, first is deposited on the IAD and then spirals toward the proto-BD according to the local disk viscosity. When the viscosity is low, gas piles up in the IAD until the temperature rises to the point where thermal ionization couples the matter to the magnetic field. The temperature at which this occurs is set to T_MRI ~ 1400 K. At this point, the MRI is triggered, moves angular momentum outward and the IAD gas spirals inward and onto the central proto-BD. The accretion rate onto the proto-BD is $${\dot{M}}_{\text{BD}}={\dot{M}}_{\text{BGR}}+{\dot{M}}_{\text{MRI}},$$(3)

where Ṁ_BGR is the quiescent accretion rate from the lAD to the proto-BD, whereas Ṁ_MRI is a significantly greater accretion rate that only occurs when the MRI operates.

The MRI is activated when enough mass has accumulated in the IAD, thus $$M_{\text{IAD}}>M_{\text{MRI}}~{\dot{M}}_{\text{MRI}}\text{Δ}{t}_{\text{MRI}},$$(4)

where $${\dot{M}}_{\text{MRI}}~5\times {10}^{-4}{\text{M}}_{\odot }{\text{yr}}^{-1}\left(\frac{{\alpha }_{\text{MRI}}}{0.1}\right).$$(5)

α_MRI is the effective Shakura-Sunyayev parameter (Shakura & Sunyaev 1973) for the MRI-provided viscosity, and Δt_MRI is the duration of the episodic event $$\begin{array}{c}\text{Δ}{t}_{\text{MRI}}\sim 0.25\text{\hspace{0.17em}kyr\hspace{0.17em}}\left(\frac{0.1}{{\alpha }_{\text{MRI}}}\right){\left(\frac{M_{\text{BD}}}{0.2{\text{M}}_{\odot }}\right)}^{2/3}\\ \times {\left(\frac{{\dot{M}}_{\text{IAD}}}{{10}^{-5}{\text{M}}_{\odot }y{r}^{-1}}\right)}^{1/9},\end{array}$$(6)

where Ṁ_IAD is the accretion rate onto the sink and onto the IAD.

Using Equations (4), (5) and (6), the condition for the MRI activation becomes $$M_{\text{IAD}}>0.13{\text{M}}_{\odot }{\left(\frac{M_{\text{BD}}}{0.2{\text{M}}_{\odot }}\right)}^{2/3}{\left(\frac{{\dot{M}}_{\text{IAD}}}{{10}^{-5}{\text{M}}_{\odot }{\text{yr}}^{-1}}\right)}^{1/9}.$$(7)

To match observations of FU Ori-type outbursts the accretion rate during the outburst is set to exponentially drop as $$\begin{array}{c}{\dot{M}}_{\text{MRI}}=7.9\times {10}^{-4}{\text{M}}_{\odot }y{r}^{-1}\left(\frac{{\alpha }_{\text{MRI}}}{0.1}\right)\\ \times \exp \left(-\frac{t-t_0}{\text{Δ}{t}_{\text{MRI}}}\right),t_0<t<\text{Δ}{t}_{\text{MRI}},\end{array}$$(8)

where t₀ marks the start of the event. After the IAD empties, the time taken to re-accumulate enough gas for another episodic accretion event is Δt_acc ~ M_MRI /Ṁ_IAD, thus $$\text{Δ}{t}_{acc}=13\text{kyr}{\left(\frac{M_{\text{BD}}}{0.2{\text{M}}_{\odot }}\right)}^{2/3}{\left(\frac{{\dot{M}}_{\text{IAD}}}{{10}^{-5}{\text{M}}_{\odot }{\text{yr}}^{-1}}\right)}^{-8/9}.$$(9)

There are two free parameters in the model: (i) the viscosity parameter provided by the MRI, α_MRI , and (ii) the quiescent accretion rate, Ṁ_BGR. The first, regulates the intensity and duration of the outburst; raising α_MRI results in a shorter (Eq. (6)), more intense outburst (Eq. (2)). Ṁ_BGR takes over as the dominant mode of gas accretion onto the proto-BD when the MRI is inactive (Eq. (3)). The values of α_MRI and Ṁ_BGR are uncertain.

3.2 BD formation – initial conditions

We simulate the collapse of a low-mass molecular cloud core so that the central object forming is a BD rather than a star. The initial density profile of the cloud is $$\rho (r)=\frac{{\rho }_{\text{c}}}{1+{\left(r/R_{\text{kernel }}\right)}^2},$$(10)

where ρ_c = 2 × 10⁻¹⁴ g cm⁻³ is the central density and R_kernel = 21 au is the radius of the core where the density is roughly uniform. The core extends to R_core = 500 au and its total mass is M_core = 0.087 M_⊙. The initial ratio of the rotational to gravitational energy is β_rot = U_rot/|U_grav| = 0.01. The initial temperature of the gas is T_gas = 10 K and thus the initial ratio of thermal to gravitational energy is α_thermal = U_thermal /|U_grav | = 0.41. We assume that the cloud core does not have any turbulence.

The cloud core is represented by 3 × 10⁵ SPH particles, thus each particle has mass m_particle = 3 × 10⁻⁷ M_⊙. The minimum mass that can be resolved is therefore M_min = N_neibm_particle = 1.8 × 10⁻⁵ M. This is two orders of magnitude smaller than the opacity limit for fragmentation, that is, the theoretical minimum mass for star formation, ~10⁻³ M_⊙ (e.g. Whitworth & Stamatellos 2006), thus the self-gravitating condensations are well resolved.

We perform seven simulations with the same cloud initial conditions. In six of them, we vary the parameters aMRI and Ṁ_BGR, which control the episodic accretion/episodic outburst events. In the remaining one, the episodic accretion model is inactive, hence the accretion and the associated radiative feedback are continuous. Observations and simulations of proto- planetary disks around low-mass stars suggest α_MRI ∈ [0.01, 0.4] (King et al. 2007; Isella et al. 2009; Bae et al. 2014; Kadam et al. 2020; Jankovic et al. 2021). In this work, we choose α_MRI = [0.01, 0.1, 0.3]. The accretion rate onto Class0/I proto-BD candidates (Palau et al. 2014; Riaz et al. 2015; Lee et al. 2018; Riaz & Bally 2021) are thought to be Ṁ_BGR ∈ [10⁻⁶, 10⁻⁹] M_⊙ yr⁻¹, thus we chose four values in this range. In Table 1, we summarize the initial conditions of our simulations.

4 Episodic accretion during brown dwarf formation

All simulations run for t = 100 kyr. In Fig. 1, we present the results of the simulations for different α_MRI (runs 1–3), and in Fig. 2 the results of the simulations for different Ṁ_BGR (runs 2, 4–6). In both figures, we also include the results of the simulation with continuous accretion, where the episodic accretion model is inactive (run 7). It is important to note that the proto-BD and the IAD are not defined for run 7.

The collapsing cloud core forms a sink at t ~ 1.85 kyr, hence a proto-BD with an inner accretion disk is acquired. The sink’s initial mass is M_sink,1.85kyr = 0.0026 M_⊙ At first, the sink accretes vigorously, with Ṁ_IAD ~ 7.8 × 10⁻⁵ M_⊙ yr⁻¹, but the accretion rate decreases as the available cloud mass is reduced. Following the formation of a sink, gas begins to accumulate in the IAD before accreting onto the proto-BD. When MRI is inactive, the forming proto-BD accumulates mass at a constant Ṁ_BGR rate. Once the MRI is triggered, Ṁ_BD rises up by ~4–6 orders of magnitude, and the IAD deposes intensively its mass on to the forming proto-BD (see Figs. 1 and 2). After the end of the event, the IAD fills up until the next episodic accretion event.

Accretion on to the proto-BD releases energy (Eq. (2)). When the MRI is inactive, the forming object’s luminosity remains constant and ≤0.1 L_⊙. However, during episodic accretion events, the luminosity dramatically increases and then declines exponentially until the event ends. The young proto-BD exhibits outbursts during the first few kyr of its life, but no outbursts occur after t = 20 kyr, because the majority of the cloud core mass has already been accreted onto the forming object. The exception is run 5, where the last outburst occurs at t ~ 80 kyr. It is therefore expected that outbursts from BDs happen within a relatively short timescale, making their discovery difficult.

In the continuous accretion case (run 7) the evolution of gas accretion is different. Initially, the sink accretes at similar rates as in other models, and the accretion rate gradually declines as the available cloud mass decreases. However, in runs with episodic accretion (1–6) the accretion rate onto the sink shows a significant drop after the last outburst, and during the subsequent evolution (except for run 5). Furthermore, in the continuous accretion run the object displays a higher luminosity than models 1–6, a pattern that is temporarily disrupted when MRI is active and an outburst occurs. Figure 3 illustrates two snapshots of the disks around the forming BD for run 2 and run 7, respectively. For continuous accretion, the higher luminosity results in a hotter and thus thicker disc. Nevertheless, the total energy emitted during each run is on the same order of magnitude, E_total ~ 10⁴⁵ erg. Episodic accretion distributes the total emitted energy into short bursts, allowing for long periods of low luminosities of proto-brown dwarfs.

Radiative feedback directly affects the accretion onto the sink, because it heats the gas, leading to a thicker disc. Consequently, the surface density drops (Fig. 3) and accretion is delayed or even temporary suppressed (e.g. Mercer & Stamatellos 2017). This is evident in Fig. 2 (2nd row, inset), where the accretion rate onto the IAD/sink drops during each outburst, but after the radiative feedback is terminated, the gas cools down quickly and the accretion rate increases again. For a longer episodic feedback, the accretion onto the sink decreases for a longer period of time. As a result, the duration and the time between successive outbursts affects the sink’s growth during the early stage of formation. In the next subsections, we explore the effect of the models’ parameters on the duration and the time intervals between successive outbursts, and the overall evolution of the BD.

Fig. 1 Luminosity of the central object, accretion rate onto the sink and proto-BD, and the masses of the proto-BD, IAD, and sink (i.e., proto-BD + lAD mass) (top to bottom row) for runs 1–3 (left to right column); see also Table 1. The black line shows the luminosity, the accretion rate onto the sink and its mass for run 7.

Fig. 2 Luminosity of the central object, accretion rate onto the sink and proto-BD, and the masses of the proto-BD, IAD and sink (i.e., proto-BD + lAD mass) for runs 2, 4–6 (top to bottom row); see also Table 1. The black line shows the luminosity, the accretion rate onto the sink and its mass for run 7.

4.1 The effect of increasing a_MRI

In Figure 4, we present the duration of successive outbursts for different viscosity parameter α_MRI values during the first 20 kyr of the simulation, as no episodic event occurs after ~17 kyr. The duration of the episodic events decreases, with increasing αMRI. As briefly explained in Sect. 3, when the MRI is active, the viscosity parameter α_MRI regulates how quickly matter accretes onto the proto-BD; the greater α_MRI, the more quickly angular momentum can be transferred outward through the IAD and the faster matter accretes onto the proto-BD. For α_MRI = 0.01 the duration of the shortest is ~156 yr; while, for α_MRI = 0.3, the duration of the lengthiest event is just ~20 yr. Consequently, the duration of the respective outbursts decreases as α_MRI increases (see also Eq. (6)).

In Figure 5, we plot the time intervals between successive outbursts. For a given α_MRI value, these intervals increase as the system evolves. This is because the mass of the proto-BD (M_BD) increases over time, while the accretion rate onto the sink (Ṁ_IAD) decreases (Eq. (9)).

Additionally, with increasing α_MRI , we observe two trends: (1) the total number of episodic events increases, and (2) the time intervals between successive outbursts decrease. Higher α_MRI values are associated with shorter outbursts and, therefore, reduced periods of radiative feedback. This reduction allows the accretion rate onto the IAD/sink (Ṁ_IAD) to decrease initially but recover over shorter timescales. As a result, the IAD can re-accumulate gas more rapidly, leading to subsequent bursts on shorter timescales. For instance, the time intervals between the first and the second events are ~1432 yr and ~369 yr, for α_MRI = 0.01 and α_MRI = 0.3, respectively. Therefore, shorter outbursts lead to shorter time intervals between successive outbursts, ultimately resulting in a greater number of events.

During an episodic accretion event, the accretion rate onto the proto-BD increases. Higher values of α_MRI result to a more significant increase of the accretion rate onto the proto-BD (see 2nd row in Fig. 1 and Eq. (8)). For α_MRI = 0.01 the accretion rate onto the proto-BD reaches a value of ~9 × 10⁻⁵ M_⊙ yr⁻¹, while for α_MRI = 0.3, maximizes at ~2 × 10⁻² M_⊙ yr⁻¹. Despite the fact, that α_MRI determines how quickly matter accretes onto the proto-BD, the amount of mass that is delivered onto the proto- BD is not directly dependent on α_MRI (Eq. (7)); nevertheless, a shorter event duration results to a more intense outburst. Hence, for greater, α_MRI the outburst brightness is also higher (Eq. (2)). For α_MRI = 0.01 the luminosity of the proto-BD maximizes at L_BD ~ 18 L_⊙ and lasts ~533 yr (see Fig. 1). For α_MRI = 0.3 the luminosity reaches up to L_BD ~ 462 L_⊙, but the duration of the event is < 20 yr, and therefore statistically difficult to be observed.

Fig. 3 Cross section of the density on the x − ɀ plane showing the disk (edge-on view) at t = 15 kyr, which forms around the proto-BD in run 2 (episodic accretion; top) and in run 7 (continuous accretion; bottom).

Fig. 4 Effect of increasing αMRI on the duration of episodic outbursts for runs 1–3; see also Table 1.

Fig. 5 Effect of increasing αMRI on the time intervals between successive episodic outbursts for runs 1–3; see also Table 1.

4.2 The effect of increasing Ṁ_BGR

The quiescent accretion rate, Ṁ_BGR, regulates the growth of the sink and the amount of energy released, when MRI is inactive (see Fig. 2 and Eq. (3)); higher Ṁ_BGR values, correspond to higher brightness released by the proto-BD. In Figure 2, we see that the sink accretion rates, the IAD and the proto-BD mass converge independently of the Ṁ_BGR value for runs 2, 4 and 6, except for run 5.

In Figures 6 and 7, we present the duration of successive outbursts and the time intervals between consecutive outbursts for runs 2, 4, 5 and 6, respectively. The duration of the respective outbursts varies by ≤15 yr for varied Ṁ_BGR, and so the duration of the outbursts is mainly independent of Ṁ_BGR . The time intervals between successive outbursts vary also by ≤15 yr for varied Ṁ_BGR, except for run 5. For the latter, the third outburst, takes place at ~2.8 kyr, and the subsequent evolution differs from the rest runs (see subset in Fig. 7). In runs 2, 4 and 6, four episodic accretion events occur no later than ~20 kyr, while in run 5 we encounter five events with the last one taking place at t ~ 80 kyr.

Overall, the duration of episodic outbursts is largely independent of Ṁ_BGR , while this parameter may affect the time intervals between episodic outbursts. However, this seems to be indirect and at least one order of magnitude less important than the effect of α_MRI . In the next sections, we constrain ourselves to the models for different α_MRI values, runs 1–3, and we explore the implications of episodic accretion in the formation of BDs.

Fig. 6 Effect of increasing ṀBGR on the duration of episodic outbursts for runs 2, 4–6; see also Table 1.

Fig. 7 Effect of increasing ṀBGR on the time intervals between successive episodic outbursts for runs 2, 4–6; see also Table 1.

4.3 Episodic accretion during the main accretion phase of BD formation

We investigate the implications of episodic accretion in the very early stages of BD formation. We focus on the timescale of the main accretion phase, which corresponds to the Class 0/I stages of low-mass star formation. We define the main accretion phase as the timescale on which the accretion rate onto the sink has decreased by ~3 orders of magnitude, similar to Machida et al. (2009). In Figure 8, we present the results for runs 1–3, 7.

The duration of the main accretion phase increases, with increasing α_MRI, as seen in Fig. 8. Higher α_MRI values correspond to shorter but more intense events (see also Figs. 6 and 1). Furthermore, the time intervals between successive events shorten as α_MRI increases, while the total number of events increases (see Fig. 7). As a result, we argue that, despite the larger number of high-accretion events, the presence of less significant but longer ones, combined with prolonged but less intense radiative feedback, leads to shorter timescales for the main accretion phase in BD formation.

Additionally, Fig. 8 shows that, regardless of the α_MRI values, the main accretion phase of the episodic accretion models 1–3 is shorter than in the continuous accretion model (run 7). In Figure 1, we see that the sink particles in runs 1–3 gain mass more rapidly (at the expense of the envelope) than the sink in run 7 during the early stages of formation, Δt ≤ 20 kyr. During these early episodic events, M_sink ≪ M_core, and after the radiative feedback from an event ends, the accretion rate onto the sink rises again to values similar to those in the continuous accretion model. However, following the final event, a clear divergence occurs. By this point, the majority of the available mass has been accreted onto the sink, while the remaining low-mass core is heated by radiative feedback, causing a significant drop in the accretion rate. Notably, the accretion rate in model 3 rises slowly afterwards, suggesting that there is still a non-negligible amount of gas in the core to be accreted onto the proto-BD. We conclude that episodic accretion during BD formation leads to a shorter timescale for the main accretion phase compared to continuous accretion.

In Figure 9, we present sink masses at the end of the main accretion phases as a function of the time of the last episodic event, for runs 1–3. We additionally show the sink mass at the end of the main accretion phase for run 7 and the results of Machida et al. (2009). The sink particles accumulate more mass when the final episodic event occurs later. Interestingly, this result is also consistent with runs 4–6. Finally, in the continuous accretion model (run 7) the sink accretes more mass than in the episodic accretion models (runs 1–3). In Table 2, we summarize our results. We conclude that episodic accretion accompanied by radiative outbursts results to lower overall gas accretion onto the proto-BD than continuous accretion.

From Figure 9, we see an approximately linear correlation between the sink mass at the end of the main accretion phase and the time of the last outburst. The later the final outburst occurs, the more massive the sink at the end of the main accretion phase.

We note that Machida et al. (2009) find a shorter duration for the main accretion phase (~10 kyr) but they investigate a system that is different to ours, including different physics: a rotating magnetized compact low-mass cloud (M_core = 0.22M_⊙). Therefore, such differences are to be expected; their result is used for reference only.

Fig. 8 Duration of the main accretion phase for runs 1–3, 7; see also Table 1. The black horizontal dashed line corresponds to run 7. The red horizontal dashed line corresponds to the results of Machida et al. (2009).

Fig. 9 Sink mass at the end of the main accretion phase versus time of the last outburst for runs 1–3, 7; see also Table 1. The horizontal dashed black line corresponds to run 7 (continuous accretion). The red horizontal dashed line corresponds to the results of Machida et al. (2009).

5 Comparison of outbursting proto-BDs with outbursting low-mass protostars

5.1 Low-mass star formation – initial conditions

We investigate the impact of episodic accretion during the early stages of low-mass star formation. We simulate the collapse of a more massive molecular cloud core so that the forming central object is a protostar, rather than a proto-BD. The cloud core’s density profile is given by Eq. (10), we use ρ_c = 2 × 10⁻¹⁶ g cm⁻³, R_kernel = 215 au, and R_core = 5500 au. Thus, the cloud’s core total mass is M_core = 1.0 M_⊙. The initial ratio of the rotational to gravitational energy is β_rot = U_rot/|U_grav| = 0.01. The initial temperature of the gas is T_gas = 10 K and thus the initial ratio of thermal to gravitational energy is α_thermal = U_thermal/|U_grav| = 0.39.

The cloud core is represented by 3 × 10⁵ SPH particles, with each particle having mass m_particle = 3.35 × 10⁻⁶ M_⊙. The minimum resolved mass now is M_min = N_neibm_particle = 2.012 × 10⁻⁴ M_⊙, and the self-gravitating condensations are well resolved (e.g. Whitworth & Stamatellos 2006).

We preform four simulations of the collapse of this more massive cloud. For runs 8–10, we choose Ṁ_BGR = 10⁻⁷ M_⊙ yr⁻¹ and a_MRI = 0.01, 0.1, 0.3, respectively. Run 11 corresponds to a model where accretion is continuous. In Table 3, we summarize the initial conditions for the simulations related to low-mass star formation.

5.2 Outbursting low-mass protostars

All simulations run for t = 1.6 Myr. In Figure 10, we present the results. We observe qualitatively similar behavior for protostars as seen for proto-BDs (see Sect. 4), though the timescales differ. The collapsing core forms a sink at t ~ 15.26 kyr, after which a protostar with an inner accretion disk is acquired. The sink’s initial mass is M_{sink,15.26kyr} = 0.0048 M_⊙. The young protostars undergo 7, 10, and 15 episodic accretion events, accompanied by outbursts, for a_MRI = 0.01, 0.1, 0.3, respectively. The majority of these events occur within the first 150 kyr after the sink formation, with the final two and three events happening later for a_MRI = 0.1 and 0.3, respectively.

In all simulations, the sink initially accretes vigorously with Ṁ_IAD ~ 4 × 10⁻⁵ M_⊙ yr⁻¹, with Ṁ_IAD gradually decreasing over time. Each outburst is marked by a temporary drop in the sink’s accretion rate due to a decrease in the disk surface density (see Fig. 11). During the quiescent accretion phases, the luminosity mostly remains within the range of [0.1–1] L_⊙. However, during episodic accretion events, it spikes up. Higher values of a_MRI result in more intense episodic outbursts, a trend similar to that observed in forming BD (see Fig. 1).

The continuous accretion model (run 11) differs from the episodic accretion models (runs 8–10) during low-mass star formation, in a manner similar to its behaviour during BD formation. Initially, the sink accretes at comparable rates across all simulations, but the accretion rate gradually declines as the available cloud core mass is depleted. Run 11 exhibits higher luminosity than models 8–10, except during outbursts, resulting in a hotter and thus also thicker disc, as can be seen in Fig. 11. Nevertheless, the total energy emitted during each run is approximately the same, E_total ~ 10⁴⁷ erg. Thus, we conclude that episodic accretion in low-mass star formation distributes the total emitted energy in a manner analogous to BD formation.

In Figures 12 and 13, we present the time duration and time intervals between successive outbursts during the first 1.1 Myr of the simulation, as no episodic event occurs afterwards. Higher values of a_MRI lead to shorter outbursts and shorter time intervals between successive events during low-mass star formation. For example, the shortest event lasts ~144 yr for α_MRI = 0.01, whereas the longest outburst for α_MRI = 0.3 lasts only~ 107 yr. Additionally, for α_MRI = 0.01, the time interval between the third and the fourth event is ~4217 yr. However, for α_MRI = 0.1 and α_MRI = 0.3, the time intervals between the third and the fourth events are ~1316 and ~328 yr, respectively. Thus, we observe a consistent trend for outbursts from low-mass stars and BDs, where higher α_MRI values are associated with shorter event durations and intervals between successive events.

Fig. 10 Luminosity of the central object, accretion rate onto the sink and protostar, and the masses of the protostar, IAD, and sink (i.e., protostar + IAD mass) (top to bottom row) for runs 8–10 (left to right column); see also Table 3. The black line shows the luminosity, the accretion rate onto the sink and its mass for run 11.

Fig. 11 Cross section of the density on the x − ɀ plane showing the disk (edge-on view) at t = 15 kyr, which forms around the protostar in run 9 (episodic accretion; top) and in run 11 (continuous accretion; bottom).

5.3 Episodic accretion during the main accretion phase of low-mass star formation

In Figure 14, we illustrate the duration of the main accretion phase as a function of α_MRI for the forming protostars. To ensure consistency in our comparison between the formation of proto- BDs and low-mass protostars, we apply the same criteria to define the main accretion phase as used in the proto-BD simulations. For α_MRI = 0.1, the main accretion phase is the longest, while for α_MRI = 0.01, the shortest. Additionally, across all values of α_MRI, the episodic accretion models (runs 8–10) exhibit a longer main accretion phase than the continuous accretion model (run 11), which is the opposite of what seen in the BD formation simulations (see Fig. 8).

In Figure 15, we show the sink mass at the end of the main accretion phases in relation to the time of the last episodic event, for runs 8–10. Notably, sink particles in the episodic accretion models accumulate more mass when the final episodic event occurs earlier. Moreover, in the continuous accretion simulation (run 11) the sink accretes more mass overall compared to the sink in the episodic accretion simulations (runs 8–10). In Table 4 we summarize our results. We conclude that episodic accretion, accompanied by radiative outbursts, leads to lower gas accretion than continuous accretion, similar to what is observed in BD formation.

A notable difference between proto-BDs and low-mass stars can be seen by comparing the sink accretion rates in the episodic and continuous cases (black vs orange lines) in Figs. 1 and 10 (middle row). In contrast to BD formation, for low-mass star formation, following the last outburst, the accretion rate onto the sink swiftly converges into the one of the continuous accretion case. We see that the effect of episodic outbursts is more prominent in BD-mass cloud cores than in more massive cloud cores, presumably due to the lower gravitational potential energy of the BD-mass cores.

Fig. 12 Effect of increasing αMRI on the duration of episodic outbursts for runs 8–10; see also Table 3.

Fig. 13 Effect of increasing αMRI on the time intervals between successive episodic outbursts for runs 8–10; see also Table 3.

Fig. 14 Duration of the main accretion phase for runs 8–11; see also Table 3. The black horizontal dashed line corresponds to run 11 (continuous accretion).

Fig. 15 Sink mass at the end of the main accretion phase versus time of the last outburst for runs 8–11; see also Table 3. The black horizontal dashed line corresponds to run 11 (continuous accretion).

6 Discussion

The MRI viscosity α_MRI is a free parameter in our models, but its value is uncertain. Therefore, observations of outburst- ing proto-BDs and protostars are needed in order to constrain it. In Figures 1 and 10, the models for α_MRI = 0.3 result to maximum outburst luminosities reaching ~10³ L_⊙ and ~10⁴ L_⊙, respectively. Vorobyov et al. (2021) presented the main characteristics of observed FUors and FUor-like objects (see Vorobyov et al. 2021, Table A.1.); since there are no observed outbursts with ≥10³ L_⊙, values of α_MRI ≥ 0.3 are probably unlikely. Values of α_MRI ~ 0.1 result predominantly to luminosities ~10¹ L_⊙, ~10² L_⊙ and outbursts duration on the order of ~10¹ yr, ~10² yr for proto-BDs and protostars, respectively, providing a better agreement with the observations. We note though that if outbursts with luminosities ≥ 10³ L_⊙ do occur, then their duration is < 20 yr (see Figs. 4 and 12), and therefore statistically difficult to be observed.

Riaz & Bally (2021) presented a near-infrared study of accretion and outflow activity in six proto-BDs. The bolometric luminosity, of their targets is in the range ~0.03–0.09 L_⊙ and was used to derive mean accretion rates, which vary from ~10⁻⁶ to ~10⁻⁸ M_⊙ yr⁻¹. In the brown dwarf formation simulations, we notice that the last episodic event separates the evolution into high-accretion and low-accretion phases, notably for models 1 and 2, see Fig. 1. This pattern is not apparent in low-mass star formation simulations, as the initial cloud core is more massive. During the high accretion phase, the accretion rate onto the proto-BD, Ṁ_IAD, ranges from ∼10⁻⁵ to ∼10⁻⁶ M⊙ yr⁻¹, whereas during low-accretion from ∼10⁻⁷ to ∼10⁻⁹ M⊙ yr⁻¹. The high- accretion period is short (≤10–17 kyr), making it statistically more difficult to observe. In contrast, it is plausible that we will frequently observe low-accretion rates, since low-accretion periods last longer. Machida et al. (2009) reach to a qualitatively similar conclusion.

6.1 ISO-OPH 200: a possible outbursting proto-BD

A detailed physical and chemical modeling of the ALMA observations for ISO-Oph 200 was conducted by Riaz & Machida (2021). Their observational and model-derived signals are compatible with the predictions of the core collapse model for BD formation (Machida et al. 2009) and they classify this object as a Class 0 proto-BD.

They identify an infalling cloud core that merges onto a circumstellar pseudo-disk and an inner Keplerian pseudo-disk, and they deduced two different outflows from the object; a wideangled low-velocity CO molecular outflow (∼1000 au), and a high-velocity (unresolved) [Fe II] atomic jet, with the former preceding the latter. The spatial scale of the CO outflow is comparable with the size of the cloud/pseudo-disk regions, resulting in the most extended CO molecular outflow identified to date.

Riaz & Machida (2021) estimate a kinematic age for the system, namely, the time since the system began to collapse to form the brown dwarf, of around 6 ± 1 kyr, and a dynamical age for the CO outflow, that is, the time since the outflow was launched by the proto-BD, at roughly 0.6 ± 0.1 kyr.

The presence of two different flows may indicate variable outflow activity. Additionally, the CO outflow’s very young dynamical age, in comparison to the proto-BD system’s kinematic age, implies a recent accretion event resulting in an outburst. This is corroborated by the ALMA CO(2–1) spectrum, which shows a sudden rise in the flux by a factor of ∼20 over just two velocity channels, and the observed spectral shape of a P Cygni profile which is typically seen during accretion bursts (Riaz & Machida 2021).

Our simulations support this argument. More specifically, the kinematic age of the system is consistent with an episodic event occurring in the first few thousand years from the formation of the proto-BD. Furthermore, according to the simulations presented here, the dynamical age of the outflow (∼0.6 kyr) is comparable to what is expected for an outburst duration for αMRI ≲ 0.01. From Fig. 4 we see that the duration of the outbursts for α_MRI ∼ 0.01 is on the order of a hundred years, while smaller α_MRI will result in even longer outbursts. On the other hand, higher α_MRI values lead to outbursts lasting only a few decades, thus they are not consistent with the dynamical age of the observed outflow.

7 Conclusions

We capitalized the semi-analytic models of Zhu et al. (2009, 2010), which have been incorporated into hydrodynamic simulations of star formation by Stamatellos et al. (2011) to study episodic outbursts from proto-BDs. The model has two free parameters: the MRI viscosity parameter α_MRI , which regulates how quickly gas flows onto the proto-BD while the MRI is active, and M_BGR, which reflects the proto-BD’s quiescent accretion rate when the MRI is dormant. These parameters are weakly restricted and may differ among proto-BDs.

We performed seven hydrodynamical simulations, all starting with an initial cloud mass of M_cloud = 0.09 M_⊙, leading to the formation of objects within the brown dwarf mass range. These simulations included one continuous accretion model and six episodic accretion models, where we varied the free parameters, α_MRI and M_BGR . In addition, we ran four additional simulations using a more massive cloud, M_cloud = 1.0 M_⊙, to investigate how episodic accretion affects brown dwarf formation and to examine the key differences and similarities compared to low-mass star formation. Our main conclusions are given below:

1. The formation timescales for proto-BDs are at least one order of magnitude shorter than for protostars;

2. During continuous and episodic accretion, the total amount of emitted energy is similar. However, episodic accretion redistributes the total emitted energy, allowing for long periods of low accretion luminosities and short luminosity outbursts. This behaviour is the same for forming proto-BDs and low-mass protostars;

3. Almost all BD outbursts occur within 20 kyr from the BD formation, whereas the majority of outbursts from low-mass stars occur within 150 kyr. Therefore, BD outbursts are more rare and hence more difficult to observe than outbursts from low-mass stars;

4. In forming brown dwarfs, episodic accretion results to shorter main accretion phase than continuous accretion. However, the opposite is true in forming low-mass stars. Additionally, higher α_MRI leads to a longer main accretion phase for brown dwarfs, but no such trend appears for low-mass stars;

5. Episodic accretion accelerates mass accretion onto to the forming object during the very early stages of its formation. However, both for brown dwarfs and low-mass stars, we find that episodic accretion results to less massive objects than continuous accretion;

6. We find a correlation between the object’s mass at the end of the main accretion phase and the timing of the last episodic accretion event. For brown dwarfs, the later this final event occurs, the more massive the forming object is. However, for low-mass stars, we find the opposite trend: the later the final event, the less massive the forming object is;

7. The effect of episodic outbursts is more prominent in BD- mass cloud cores than in more massive cloud cores. For brown dwarf forming cloud cores, the last outburst essentially splits the evolution into a short high-accretion and a much longer low-accretion phase.

Acknowledgements

We would like to thank the anonymous referee for their constructive review that has helped to improve the paper. BR acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) – Projekt number RI-2919/2-3. DS acknowledges support from STFC grant ST/Y002741/1.

References

All Tables

All Figures

	Fig. 1 Luminosity of the central object, accretion rate onto the sink and proto-BD, and the masses of the proto-BD, IAD, and sink (i.e., proto-BD + lAD mass) (top to bottom row) for runs 1–3 (left to right column); see also Table 1. The black line shows the luminosity, the accretion rate onto the sink and its mass for run 7.
In the text

	Fig. 2 Luminosity of the central object, accretion rate onto the sink and proto-BD, and the masses of the proto-BD, IAD and sink (i.e., proto-BD + lAD mass) for runs 2, 4–6 (top to bottom row); see also Table 1. The black line shows the luminosity, the accretion rate onto the sink and its mass for run 7.
In the text

	Fig. 3 Cross section of the density on the x − ɀ plane showing the disk (edge-on view) at t = 15 kyr, which forms around the proto-BD in run 2 (episodic accretion; top) and in run 7 (continuous accretion; bottom).
In the text

	Fig. 4 Effect of increasing αMRI on the duration of episodic outbursts for runs 1–3; see also Table 1.
In the text

	Fig. 5 Effect of increasing αMRI on the time intervals between successive episodic outbursts for runs 1–3; see also Table 1.
In the text

	Fig. 6 Effect of increasing ṀBGR on the duration of episodic outbursts for runs 2, 4–6; see also Table 1.
In the text

	Fig. 7 Effect of increasing ṀBGR on the time intervals between successive episodic outbursts for runs 2, 4–6; see also Table 1.
In the text

	Fig. 8 Duration of the main accretion phase for runs 1–3, 7; see also Table 1. The black horizontal dashed line corresponds to run 7. The red horizontal dashed line corresponds to the results of Machida et al. (2009).
In the text

	Fig. 9 Sink mass at the end of the main accretion phase versus time of the last outburst for runs 1–3, 7; see also Table 1. The horizontal dashed black line corresponds to run 7 (continuous accretion). The red horizontal dashed line corresponds to the results of Machida et al. (2009).
In the text

	Fig. 10 Luminosity of the central object, accretion rate onto the sink and protostar, and the masses of the protostar, IAD, and sink (i.e., protostar + IAD mass) (top to bottom row) for runs 8–10 (left to right column); see also Table 3. The black line shows the luminosity, the accretion rate onto the sink and its mass for run 11.
In the text

	Fig. 11 Cross section of the density on the x − ɀ plane showing the disk (edge-on view) at t = 15 kyr, which forms around the protostar in run 9 (episodic accretion; top) and in run 11 (continuous accretion; bottom).
In the text

	Fig. 12 Effect of increasing αMRI on the duration of episodic outbursts for runs 8–10; see also Table 3.
In the text

	Fig. 13 Effect of increasing αMRI on the time intervals between successive episodic outbursts for runs 8–10; see also Table 3.
In the text

	Fig. 14 Duration of the main accretion phase for runs 8–11; see also Table 3. The black horizontal dashed line corresponds to run 11 (continuous accretion).
In the text

	Fig. 15 Sink mass at the end of the main accretion phase versus time of the last outburst for runs 8–11; see also Table 3. The black horizontal dashed line corresponds to run 11 (continuous accretion).
In the text