© The Authors 2025

1 Introduction

Conventionally, transition disks (TDs) are protoplanetary disks (PPDs) that have inner dust cavities, as indicated by a deficit in their spectral energy distribution (SED, Espaillat et al. 2014). Initially, the term transition was intended to illustrate a shift from Class II to Class III, but subsequent studies have proposed that the cavity could also represent a morphological aspect of the disk (see, e.g., the discussions in Wang & Goodman 2017; Ercolano & Pascucci 2017; van der Marel 2023).

The source of these cavities remains an open question – and proposed mechanisms are likely not exclusive (van der Marel 2023). One of the primary and most frequently explored mechanisms is gap formation by an embedded companion or planet (Marsh & Mahoney 1992; Rosotti et al. 2016; Price et al. 2018) – but only recently has this indeed been observed (e.g., Keppler et al. 2018). Another mechanism that has extensively been considered is photoevaporation, in a combination with viscous evolution (Clarke et al. 2001; Alexander et al. 2014; Ercolano & Pascucci 2017; Picogna et al. 2019). Additionally, dust migration driven by radiation pressure has been proposed as a contributing factor – although no consensus has been reached (e.g., Dominik & Dullemond 2011; Owen & Kollmeier 2019; Bi & Fung 2022).

Apart from the dust cavity, there are clear indications of the presence of gas within the dust-depleted region (Najita et al. 2007; Manara et al. 2014). While the amount of gas in the cavity is not yet certain, recent studies suggest a drop of 2–4 orders of magnitude in the gas surface density profile (van der Marel et al. 2018). Despite the uncertainty in the origin and evolution of the cavities, one recently emerging characteristic is that the accretion rate in a large number of TDs appears to be comparable to that of full disks (Manara et al. 2014; Owen 2016). Various explanations have been proposed for this phenomenon. One possibility is transonic gas accretion (Wang & Goodman 2017), where gas crosses the cavity. It has also been proposed that the high accretion rates are a result of a hidden mass reservoir in the innermost region of the disc. Models that consider this are usually also incorporating a dead zone (Gárate et al. 2022); however detailed modeling suggests that a pressure bump does not always occur at the inner side of a dead zone (Delage et al. 2022).

The majority of the models mentioned above use the αprescription (Shakura & Sunyaev 1973) for a viscously evolving disc. While the question of disk evolution and whether it is governed by turbulent transport or MHD winds is still open (see Manara et al. 2023, for a recent review), it is increasingly clear that magnetic fields play an important role in the evolution of the system: on the one hand, Ohmic resistivity (OR) and ambipolar diffusion (AD) have been shown to severely suppress the magnetorotational instability (MRI) (Bai & Stone 2011; Gressel et al. 2015), while on the other hand, magnetohydrodynamic (MHD) winds provide an effective alternative for angular momentum transport and gas accretion (Bai 2017; Béthune et al. 2017; Gressel et al. 2020).

Additionally there is the added complexity of the Hall effect (HE). The HE essentially reorganizes the magnetic field topology with vastly different outcomes depending on the magnetic field polarity, ranging from suppressing turbulence or instability (Wardle & Salmeron 2012) to introducing a new one (Hall shear instability, Kunz 2008). It has also been shown that the HE can enhance the radial and toroidal magnetic fields, thereby increasing the accretion rate through the laminar Maxwell stress (Bai & Stone 2017; Sarafidou et al. 2024).

Relatively few models of transition disks with magnetic fields exist. Early on, Combet & Ferreira (2008) presented a theoretical model of a magnetic wind-driven inner disk via what they termed a jet-emitting disk (JED), demonstrating that it can indeed sustain accretion rates similar to the outer turbulent disk – a situation reminiscent of transitional disks, as they commented upon. Later on, secular evolution models of (partially) wind-driven accretion allowed Suzuki et al. (2016) to find that a surface-density-dependent torque can indeed create and sustain a central cavity. Based on static wind models in the framework of a weakly ionized plasma, Wang & Goodman (2017) took into account the magnetic diffusivities in a TD and showed that they indeed allow for the launching of a magnetic wind – and furthermore that this would lead to transonic gas accretion. Finally, Martel & Lesur (2022) performed the first self-consistent numerical model of an accreting TD with a non ideal MHD disk wind. They confirmed that such a model cannot only maintain a cavity over thousands of orbits but also predict high accretion rates, indeed comparable with the sound speed within the cavity.

In this work, we present the first radiative non ideal MHD (RMHD) model of a TD, featuring several novel aspects: we include simplified stellar irradiation, diffuse re-radiation (Gressel 2017), as well as far ultra-violet (FUV) heating via a simple thermochemistry model (Gressel et al. 2020). We moreover feature all three non ideal MHD effects, with diffusivities calculated depending on local properties of the gas and shielded exterior ionizing sources (Gressel et al. 2015, 2020). Our goal is to examine which of the previously reported results hold under this more complex approach and to identify any new dynamical behaviors that may arise.

The remainder of this paper is organized as follows: in Section 2 we present the model equations, physical quantities and the numerical setup. Then, in Section 3, we proceed and present the results of our runs. We discuss the accretion rate (and its origin), the self organization of the flow in the cavity region, as well as the issue of magnetic flux transport. In Section 4, we discuss our results and compare them with existing work. Finally, we present our conclusions in Section 5.

2 Methods

We perform radiation MHD simulations of the inner regions of PPDs, refining the non ideal MHD disk model of Gressel et al. (2020), that featured stellar irradiation and a simplified thermochemistry. As the most notable improvement over that work, we now also include the Hall effect in the induction equation, which is implemented by means of the framework presented in Krapp et al. (2018). Additionally, since we aim to study TDs, we have introduced a cavity – that is, a drop of several orders of magnitude in the gas surface density – in our disk. Note that unlike Sarafidou et al. (2024) we do not include X-ray photoionization, in order to keep continuity with Gressel et al. (2020).

Our simulations are 2D axisymmetric with a spherical-polar coordinate system (r, θ, ϕ) – radius, colatitude, and azimuth, respectively. We are primarily interested in the inner regions of the PPD where the disk is considered to be laminar under typical conditions. Our computational domain covers an area of r ∈ (0.75, 22.5) au, θ ∈ (0, π) with a standard grid resolution of N_r × N_θ = 320 × 256 cells. All simulations were performed with the single-fluid NIRVANA-III code, which is built around a standard second-order accurate finite volume scheme (Ziegler 2004, 2016). We use a modified version of the NIRVANA-III code that has been described in detail in Gressel et al. (2020). To obtain the interface states entering the magnetohydrodynamic fluxes, we use the HLLD approximate Riemann solver by Miyoshi & Kusano (2005). The Hall term is added via operator-splitting, which moreover allows sub-stepping in time when required. The update due to the Hall electromotive force is done using the Hall Diffusion Scheme (HDS) of O’Sullivan & Downes (2006, 2007).

2.1 Equations of motion

We solve the equations of motion for mass-, momentumand total energy-density in conservation form along with the induction equation for the magnetic field, that is, $$\begin{array}{l}{\partial }_t\rho +\nabla \cdot (\rho v)=0,\\ {\partial }_t(\rho v)+\nabla \cdot \left[\rho vv+PI-BB/{\mu }_0\right]=-\rho \nabla \mathrm{\Phi },\\ {\partial }_{t}e+\nabla \cdot \left[(e+P)v-(v\cdot B)B/{\mu }_0\right]=-\rho (\nabla \mathrm{\Phi })\cdot v+\nabla \cdot \mathcal{S}+Q^{+/-},\\ {\partial }_{t}B=-\nabla \times ℰ,\end{array}$$

where the total pressure, P, is defined as the sum of the gas and the magnetic pressure, B²/2µ₀· Finally, −ρ∇Φ denotes the gravity force of the central star.

2.1.1 Non ideal MHD

In the above energy equation, we have a divergence of the Poynting flux, 𝒮 ≡ ℰ× B/µ₀, with ℰ representing the electromotive force responsible for the evolution of the magnetic field. It takes into account induction, as well as the three non ideal MHD effects, that is, $$ℰ=-v\times B+{\mu }_0\left[{\eta }_{\text{O}}J+{\eta }_{\text{H}}J\times \widehat{b}-{\eta }_{\text{A}}(J\times \widehat{b})\times \widehat{b}\right],$$(1)

where η_O , η_A and η_H denote the diffusion coefficients of Ohmic resistivity, ambipolar diffusion and Hall effect, respectively, J = (1/µ₀)∇ × B is the electric current density, $\hat{\mathbf{b}}$ the unit vector of the magnetic field, and µ₀ the magnetic permeability.

For the dynamically evolving diffusion coefficients η_O , η_H , and η_A, we use a lookup table based on Ilgner & Nelson (2006) and Desch & Turner (2015); details regarding the ionization model can be found in Section 2.5 of Gressel et al. (2020). Moreover, in order to avoid severe timestep constrains we limit η_A and η_H such that Λ_Aβ_P ≥ 10⁻³ and Λ_H β_P ≥ 10⁻², respectively¹. Additionally, we enable the Hall effect only after 30 yr, so as to allow for a less problematic initial relaxation of the magnetic field in our simulation (also see Bai 2017).

2.1.2 Radiation and thermodynamics

Thermodynamic heating cooling terms are subsumed in the energy source and sink term $$Q^{+/-}\equiv -c\rho {\kappa }_{\text{P}}\left(a_{\text{R}}{T}^4-e_{\text{r}}\right)+Q_{\text{irr}}^{+}+Q_{\text{pdr}}^{+/-},$$(2)

where e_r is the radiation energy density, a_R ≡ 4 σ/c is the radiation constant, σ is the Stefan-Boltzmann constant, κ_P is the Planck mean opacity, and c is the speed of light.

The terms ${\mathcal{Q}}_{\text{irr}}^{+}$ and ${\mathcal{Q}}_{\text{pdr}}^{+/-}$ describe the heating (and cooling) via stellar irradiation and a very simplified FUV photo chemistry, respectively. For advancing the radiation energy density, e_r, in time, we use the approached detailed in Gressel (2017) as well as Section 2.3 of Gressel et al. (2020), which essentially comprises stellar irradiation and diffusive re-radiation (including a flux limiter approach combined with a two-stream approximation to address the optical-depth transition in the vertical direction). Note that, for brevity, we have omitted the radiation force in the above momentum equation, which in general is entirely negligible in our setup. This means that radiative effects are limited to the thermodynamic aspects of our problem, not the dynamic ones.

2.2 Disk model and boundary conditions

The equilibrium structure of our models is described by a power law relationship of the locally isothermal temperature T, with the cylindrical radius, R, that is, T(R) = T₀ (R/R₀)^q. A similar relation holds for the midplane density, i.e., ρ_mid(R) = ρ₀ (R/R₀)^p, where we have chosen a temperature slope of q = −0.75 (resulting in a mildly flared disk height as a function of radius), and a density slope of p = −1.5 – these are the same values as used in Gressel et al. (2020). The initial solutions for the disk equilibrium (see Nelson et al. 2013) then become $$\rho (r)={\rho }_{\text{mid\hspace{0.17em}}}(R)\exp \left[\frac{G{M}_{\star }}{c_{\text{S}}^2(R)}\left(\frac{1}{r}-\frac{1}{R}\right)\right],$$(3) $$\tilde{\mathrm{\Omega }}(r)={\mathrm{\Omega }}_{\text{K}}(R){\left[(p+q){\left(\frac{H}{R}\right)}^2+(1+q)-q\frac{R}{r}\right]}^{1/2},$$(4)

where ${\mathrm{\Omega }}_{\text{K}}(R)=\sqrt{G{M}_{\star }}{R}^{-3/2}$ is the well-known Keplerian angular velocity, and $c_{\text{s}}^2(R)=c_{\text{s}0}^2{\left(R/R_0\right)}^q$ is the sound speed squared. The overall thickness of the disk is controlled by the parameter H(R₀) = c_s0/Ω_K(R₀), and we define H(R₀) = 0.05 R₀, which fixes the overall pressure scale height of the initial gas disc². The initial density profile was scaled to Σ₀ = 340 g cm⁻² at R₀ = 1 au, and the aspect ratio is H(R) = 0.05(R/R₀)^1/4, the mass of the central star M_⋆ = 1 M_⊙, the X-ray luminosity of the star is L_X = 10₂₉ erg s⁻¹, and the FUV luminosity, L_FUV = 0.2 L_⊙, is set to be equivalent (see Eq. (6) in Gorti et al. 2009) to a mass accretion rate onto the star of 10⁻⁸ M_⊙ yr⁻¹ . We furthermore adopt a mean molecular weight of µ = 2.4 and an adiabatic index γ = 1.4, representing a diatomic mixture of molecular hydrogen and helium gas.

In order to mimic a central cavity in our disk, we multiply the profile ρ_mid with a function f (R): $$f(R)=a_{\text{tr}}+\frac{1-a_{\text{tr}}}{2}\left[1+\tanh w_{\text{tr}}\left(R-R_{\text{tr}}\right)\right],$$(5)

with a contrast parameter, a_tr = ρ_gap/ρ_disk in the midplane, Rtr is the transition radius, where the cavity ends and the disk begins, and wtr controls the steepness of the cavity edge. For all models we have set a_tr = 10⁻⁵, wtr = 2, and R_tr = 8.

We here use the same boundary conditions as in our previous work (Sarafidou et al. 2024) – for details we refer the reader to Section 3.2.1 of Gressel et al. (2020). At radial boundaries, we apply a standard outflow condition for the momentum density, and extrapolate the mass density with the power-law slope of the equilibrium model. Moreover, standard wave-killing zones are implemented at the inner and outer domain boundary. For the magnetic field components parallel to the domain boundary, we generally enforce vanishing gradients. The remaining perpendicular component is in turn computed from the div(B) = 0 constraint. The outer radial boundary is an exception to this; here we set the azimuthal components of both the field and current to zero. Finally, the axial boundary simply enforces the demanded axisymmetry of our simulations.

2.3 Integration and averages

Vertical integration is done along a cylindrical height inside the disk, $$\overline{Q}(R,t)={\displaystyle {\int }_{-{ɀ}_d}^{+{ɀ}_d}Q}(R,ɀ,t)\text{d}ɀ$$(6)

where the integration limits ±z_d are defined as the height where β_P = 0.1, a value indicating that there is significant magnetic energy to launch a wind (Pascucci et al. 2023).

All models were run for approximately 700 years, ensuring enough orbital cycles at the transition radius to analyze the cavity dynamics, while taking into account the computational cost (we address the run-time of our simulations in Appendix A).

Time averaging, $$\langle Q(R,ɀ)\rangle =\frac{1}{T}{\displaystyle {\int }_{t_0}^{t_0+T}Q}(R,ɀ,t)\text{d}t$$(7)

when mentioned, has been done from 450 to 700 years, unless stated otherwise.

2.4 Simulation setup

The only existing MHD simulations of TDs by Martel & Lesur (2022) were restricted to the case of AD and OR. Apart from our more detailed thermodynamic and diffusivity treatment, we improve over this work by additionally including the Hall effect (HE) in our simulations. We show the profiles of the non ideal effects in Appendix B.

To compare the influence of the HE, we have performed three runs with and without the HE, but with otherwise identical input parameters: besides the disk’s power laws being the same, the three runs have also been initialized with the same vertical magnetic field of β_P = 10⁴ at the disks midplane. The non ideal MHD effects were incorporated as follows: We have two runs that include all three effects (OR,AD, and HE), one with $B\cdot \widehat{z}<0$ (Hall anti-aligned, h-aa) and the other with $B\cdot \widehat{z}>0$ (Hall aligned, h-a); the third run includes only OR and AD (Hall free, h-f ) – these effects are both invariant under reflection of the field direction, such that it suffices to study one of the field configurations in that case. The mutual orientation of the magnetic field with the rotation axis becomes significant since the HE (in the presence of a vertical field) amounts to a rotation of the horizontal field components in the presence of a horizontal current. Because the sense of rotation depends on the field orientation, this can affect the overall dynamics of the flow.

Fig. 1 Time-averaged snapshots of the gas density (top row), plasma parameter β (middle row), and toroidal magnetic field (bottom row), for: the initial model (first column), as well as the three runs with Hall anti-aligned (h-aa, second column), Hall free (h-f , third column), and Hall aligned (h-a, fourth column). The streamlines in the first row indicate the poloidal velocity vectors. The blue lines in the middle row are placed where β = 0.1 indicating the vertical integration limits ±zd. Finally, the contour lines in the bottom-row panels highlight the poloidal magnetic field lines. The initial radial location of the cavity edge is marked by a blue diamond in all panels.

3 Results

3.1 General observations

In the following, we describe some broad characteristics of our simulations. Looking at the density snapshots in the top row of Figure 1, we can see that, over the duration of the simulation, the overall extent of the cavity remains roughly stationary for all three runs. While the density contrast with the outer disk is somewhat reduced (compared with the initial setup, shown in the leftmost panel), the cavity remains clearly discernible. Notably, MHD winds originate both in the gap region and in the inner part of the outer disk (at r ≃ 9 au)³. Itis important to note that, in comparison to Sarafidou et al. (2024), we do not enforce reflection symmetry with respect to the midplane. As a result, all our models develop a slight top–bottom asymmetry, especially in the case of the h-a, where stronger horizontal fields are produced.

In the middle row of Fig. 1, we plot the evolved magnetization (again in terms of β_P) for the three cases. We can see that the evolved models features distinct levels of magnetization depending on the assumed microphysics and field orientation: In the case of the h-aa, the disk remains close to the initial β_P, with β_p ≥ 10³ in the midplane. Contrary to this, in the case of the h-a, the disk has become significantly magnetized in the most part, with β_P ~ 10² (in the midplane). The central cavity follows a similar pattern, with h-a being the most magnetized, while for the h-aa we note a small demagnetized area – this issue will be discussed in more in detail in Sect. 3.3, below.

In the bottom row of Fig. 1, we plot the toroidal magnetic field, B_ϕ . Poloidal field lines are over-plotted in terms of contour lines of RA_ϕ , where A_ϕ the toroidal component of the vector potential (that we evolve alongside the magnetic field variable; see Gressel et al. 2020, for details). In all three runs, we see symmetry breaking to various degrees. By eye, the h-aa appears to retain symmetry better than the other two cases. We moreover see, at the same time, areas devoid of field lines, as well as featuring field line accumulations. Additionally we note that in the h-a case, we have the strongest Bϕ , while for the h-aa model we have the weakest. This is linked to the difference in the disk magnetization and is because of the general tendency of the Hall effect and its mutual interplay with the Keplerian shear.

In the h-a case, there will be an enhancement of the buckling of the field lines, while on the h-aa any curvature of the field lines will be severely suppressed (Wardle & Salmeron 2012). In our models, the magnetic field lines are threading the MHD wind, thus bending and introducing a radial component to the magnetic field. This then gets sheared generating a toroidal magnetic field (see also Bai & Stone 2017; Sarafidou et al. 2024).

Fig. 2 Radially averaged accretion rates through the cavity region over time. Shaded areas indicate uncertainties due to the vertical variability.

3.2 Accretion rates in the low-density region

A central question that we aim to address with our simulations is whether the central cavity can maintain the overall accretion process onto the star, or whether it acts to isolate the outer disk instead. In Fig. 2, we accordingly compare the three models in terms of their mass throughput in the central cavity region over time; the radial accretion rate is simply defined as: $${\dot{M}}_{\text{a}}=2\pi R\overline{\rho v_{\text{R}}}$$(8)

where we average over the radial extent R ∈ [1,7] au.

As is immediately noticeable, we obtain the highest accretion rate for the h-a case. The h-f and h-aa have almost the same accretion rate in the first 200 yr. After that, though, the h-aa produces a somewhat lower M_a. More specifically, the accretion rates for the h-aa, h-f and h-a are: (2.63 ± 0.31) × 10⁻⁸, (5.13 ± 0.41) × 10⁻⁸, and (7.91 ± 0.22) × 10⁻⁸ M_⊙ yr⁻¹, respectively. As has been discussed before (see Bai & Stone 2017; Sarafidou et al. 2024), the fact that the h-a has the highest overall mass accretion rate is a direct result of the aligned Hall effect: the increased radial and toroidal magnetic components lead to a higher overall Maxwell stress, and thus a higher accretion rate.

For a more detailed analysis of the accretion rate of the disk, we turn to the equation of angular momentum conservation, in its vertically integrated form: $$\begin{array}{l}{\dot{M}}_{\text{a}}=\frac{2\pi R}{\partial /{\partial }_R\left(\tilde{\mathrm{\Omega }}{R}^2\right)}(\frac{1}{R}\frac{\partial }{\partial R}{R}^2\left[\overline{\rho v_{\varphi }^{\prime }{v}_{\text{R}}}-\frac{\overline{B_{\varphi }{B}_{\text{R}}}}{{\mu }_0}\right]\\ +{R\left(\rho v_{\varphi }^{\prime }{v}_{\text{ɀ}}-\frac{B_{\varphi }{B}_{\text{ɀ}}}{{\mu }_0}\right)|}_{-{ɀ}_d}^{+{ɀ}_d}),\end{array}$$(9)

where the mass accretion rate, Ṁ_a, is denoted $-2\pi R\overline{\rho v_{\text{R}}}$. In that expression, we have $ℛ=\overline{\rho {v^{\prime }}_{\varphi }{v}_{\text{R}}}+{R\rho {v^{\prime }}_{\varphi }{v}_{\text{z}}|}_{-z_d}^{+z_d}$ as the contribution from the (horizontal and vertical) Reynolds stress, and, moreover, $ℳ=-\overline{B_{\varphi }{B}_{\text{R}}}/{\mu }_0-R{B}_{\varphi }{B}_{\text{z}}/{{\mu }_0|}_{-z_d}^{+z_d}$ from the relevant (horizontal and vertical) Maxwell stress components, respectively. We normalize these quantities by the vertically integrated pressure $\overline{p}$ noting them as ${\alpha }_{R\varphi }^{ℛ}={ℛ}_{R\varphi }/\overline{p}$ and so on. Finally, ${v^{\prime }}_{\varphi }=v_{\varphi }-\tilde{\mathrm{\Omega }}R$ defines deviations from equilibrium.

In the top panel of Fig. 3, we show the time-averaged accretion rate over radius for the three runs, calculated from Eq. (8). Apart from the “dip” near the transition radius⁴, in all cases the accretion rates of the cavity and the outer disk are comparable. Looking specifically in the cavity, In the h-f case, the accretion rate is relatively constant, which to a lesser degree is also true for the h-aa run, where we mainly see two distinct levels of accretion for the inner and outer half of the region of interest, respectively, and a shallow ramp in-between. In contrast, for the h-a case, we observe a marked dip in the accretion rate at around r ~ 2.5–3 au – and otherwise a similar gradual step at R ≃ 4 au as for the h-aa model. The sudden drop in the accretion rate provides an indication that we might observe a substantial level of local mass accumulation just inside this radius, something that we subsequently discuss in detail in Sect. 3.3. In all cases, the accretion rate stems predominantly from the MHD wind – see the middle panel of Fig. 3. We see specifically that the surface term of the Maxwell stress is dominant compared to the in-plane stress (as shown in the bottom panel of the figure).

We highlight the specific geometry of the MHD winds in Figure 4, where we plot the magnetic field components (in cylindrical coordinates) for the three runs (top to bottom panels). For all the cases, the horizontal (purple/yellow) and vertical (green) components are roughly comparable in magnitude with each other. This has previously been highlighted in the work of Wang & Goodman (2017). As they mention, if the magnetic field components are comparable with each other, then the wind inherently has a simple geometric advantage in driving accretion compared to radial (“viscous”) stress. This is because the surface stress acts on a much larger surface area (i.e., ~ πR²) compared to the radial stress (i.e., ~ 2πR H). Their analysis further suggests that in such scenarios the accretion velocity is expected to be comparable to the speed of sound. This is further supported when looking at the accretion rates of the models of Ṁ_a ~ 3-8 × 10⁻⁸ M_⊙yr⁻¹ – values that are well within the range of expected accretion rates for full disks. Hence if the rates are comparable between full disks and transition disks but the surface density drops by at least two orders of magnitude, then the accretion velocity should be significantly increased.

We distinctly showcase the characteristic feature of nearsonic accretion in Figure 5. Looking at the accretion velocity profiles, we see that there is a clear transition from subsonic flow (v_acc/c_s ≤ 10⁻²) in the outer disk to near-sonic flow (v_acc/c_s ≥ 10⁻¹ ) in the central cavity – where specifically for R < 2.5 au we observe transonic accretion for all three runs.

Apart from R < 2.5 au, where there are non-negligible inplane Maxwell (and indeed Reynolds) stresses – see the middle and bottom panels in Fig. 3 – accretion actually happens away from the disk midplane. This is a consequence of two fundamental aspects of our simulations: a) unlike in classic wind models (Blandford & Payne 1982), the β_P = 1 line is at several scale heights above the disk, and b) the bending of field lines is suppressed by OR in the body of the disk. Regarding this topological feature of the accretion flow, we highlight that in all three runs we see that a streamer originates from the inner edge of the outer disk (at R ∼ 8 au) and then traverses the low-density cavity preferring either top or bottom (θ : 0–π/2 or θ : π/2–π), respectively. At the opposite location, an outward mass flux is typically observed – although not nearly as strong as the corresponding accretion flow through the other disk half. This can be seen in the bottom row of Fig. 5 – note the (signed) logarithmic y-axis in that plot. It is worth mentioning that there seems to be no distinction in the accretion velocity profile between the two different configurations of the HE. Moreover, the position (top/bottom) of the inward streamer is not fixed in time; it varies whenever the current sheet flips its orientation due to the non-vanishing toroidal magnetic field in the disk midplane.

Fig. 3 Time-averaged accretion rates and stress terms versus radius for the three runs. In the first panel, we show the measured accretion rates, in the middle panel we show the Reynolds (dashed) vs. the Maxwell stress (solid line), normalized by the vertically integrated pressure. In the bottom panel, we decompose the normalized Maxwell stress into its radial (dashed) and surface (solid line) components.

Fig. 4 Vertical profiles of the magnetic field components along a poloidal field line (each highlighted by a thick blue line in the bottom panels of Fig. 1, above) threading the cavity.

3.3 Emergent structures within the cavity

The time-averaged density snapshots of Figure 1 gave a clear indication that – at least overall – the density contrast between the inner cavity and the outer disk dynamically persists in all three runs. This is by no means a given. For a closer look, we provide the space-time plot of the vertically integrated density in Figure 6.

The space-time plot reveals that the gas in the cavity is quite dynamic, with burst of matter appearing (especially in the h-f case – cf. the middle panel of Fig. 1). We also note that in all cases the density of the cavity increased: It was initialized with a_tr = 10⁻⁵ , but after ∼50 yr it was partly filled-up. For the h-f case, the cavity contrast, atr, remains at a value a_tr ≃ 10⁻³ – and this sits well within the range of observational estimates for the amount of gas within the gap (van der Marel et al. 2018). For the two HE cases, though, a different behavior is observed: gas accumulates near the center of the cavity, forming a ring-like structure – this notion is strengthened by Fig. 7, where we clearly see that the two HE runs display broad peaks in their radial profiles. The one formed in the h-aa case is the smallest, spanning from (2.5–5.5) au. The ring in the h-a case is slightly larger, from (2.0–5.5) au.

Since these ring-like structures only appear in the HE runs, we explore the Hall-MHD self organization mechanism that was introduced in Kunz & Lesur (2013) and further detailed in Béthune et al. (2016). This mechanism requires the system to be in the “strong” Hall regime, characterized by the dimensionless parameter ℒ = l_H/H ≳ 0.2, where l_H = η_H/v_A is the ”Hall length”, relating the Hall coefficient η_H to the Alfvén speed $v_{\text{A}}=\left|B\right|/\sqrt{{\mu }_0\rho }$ – also cf. Kunz & Lesur (2013). Moreover, the in-plane component, ℳ_Rϕ , of the Maxwell stress tensor is characteristically found anticorrelated with the vertical magnetic field. This arises because the Maxwell stress appears in the induction equation through the Hall term. According to Béthune et al. (2016), the Hall-MHD induction equation can be expressed as $$\frac{\partial B_{ɀ}}{\partial t}\simeq {(\nabla \times [v\times B])}_{ɀ}-l_{\text{H}}\frac{\partial }{\partial R^2}{B}_R{B}_{\varphi }.$$(10)

This illustrates that – in the presence of a strong HE – the curvature of the in-plane Maxwell stress directly influences the induction equation, indicating a strong link between angular momentum transport and vertical magnetic flux: when the second derivative of ℳ_Rϕ is negative, vertical magnetic flux is pushed outward. Consequently, the vertical magnetic field B_ɀ should be anticorrelated with ℳ_Rϕ , where a local maximum in ℳ_Rϕ corresponds to a local minimum in B_ɀ .

While in the cavity we are indeed in the strong Hall regime – see Figure 8 – the radial Maxwell stress is not anti-correlated with Bz . From Figure 9, it in fact looks like they are actually correlated, instead.

The only other significant difference between the Hall and Hall-free runs, which might explain the ring formation in the former but not the latter, is observed when comparing the accretion rates from the wind to those from the radial stress, as shown in Figure 10.

In regions with ring formation, only in the Hall cases do we observe Ṁ_rad > Ṁ_wind, indicating that in those areas, more mass is replenished in the gap than is being ejected by the wind. This is remarkably similar to the wind driven instability that was described in Suriano et al. (2017); Riols & Lesur (2019): In non ideal MHD, as mass is accreted, magnetic flux diffuses away from the region of mass concentration, leading to a decrease in magnetic flux where density accumulates (see the vertically integrated profiles of density and B_ɀ in Fig. 9). This results in a dynamically weakened magnetic field, that in turn increases the accretion timescale in that region. Consequently, this creates an “accretion trap” for mass that is transported from larger radii (most apparent in the radial accretion profile of the h-a in the top panel of Fig. 3). Magnetic flux is then accumulated in the under-dense regions strengthening the wind and more effectively removing angular momentum. The reason why this is only happening in the Hall runs is currently unclear.

Fig. 5 Time-averaged snapshot of the flow Mach number over 100yr (upper panels) and normalized accretion velocity profiles (lower panels) for the three runs. On the snapshots we are overlaying the radial mass flux using a gray alpha gradient colormap. The black dashed lines indicate the disk midplane as well as the disk surface.

3.4 Magnetic topology and flux evolution

There is a broad consent that the efficiency of any accretion and/or mass ejection mechanism that relies on magnetic forces has to somehow scale with the magnetic field strength. This equally holds for torques from turbulent (laminar) Maxwell stresses, as well as for magnetic tension (pressure gradient) forces in magnetocentrifugal (pressure-driven) disk winds. To obtain a robust handle on disk evolution, one hence needs to understand the topology and evolution of the entrained magnetic flux, that is threading the disk.

In this section we thus want to look at the temporal evolution of the magnetic flux. To do that we utilize the poloidal magnetic flux function $$\psi =r\sin (\theta )A_{\varphi },$$(11)

where the azimuthal component, A_ϕ, of the magnetic vector potential is evolved alongside the induction equation. During our NIRVANA-III simulation runs, the simple equation (also cf. Sect. 2.2.1 in Gressel et al. 2020) $${\partial }_t{A}_{\varphi }=-{ℰ}_{\varphi },$$(12)

is solved, with ℰ given by Eq. (1) in Section 2. Because isocontours of the flux function Ψ describe poloidal field lines, we can use it to study the redistribution of magnetic flux. Looking at the space-time plots in Fig. 11, we see that for the h-f case (middle panel) and – to a somewhat lesser degree – for the haa case (top panel), the flux appears almost stationary over the timescale of our simulation runs. In contrast, for the h-a case, we have a different picture. We observe the flux in the outer disk (i.e., for R > 8 au) being pushed outward. Moreover, on a secular timescale, a lot of fluctuations emerge near the transition radius. At later times, this transient variability appears to have settled down, though, and gives way to a quasi-stationary state of the field evolution.

To gain a better insight into the efficiency and timescale of the transport, we use the radial motion of the vertical magnetic flux in the disk midplane as a simple proxy. To this end, we define the magnetic flux velocity (see also Bai & Stone 2017) relative to the Kepler velocity, v_K, that is $$v_{\text{B}}=\frac{{ℰ}_{\varphi }}{B_{\text{Z}}{v}_{\text{K}}}.$$(13)

In Figure 12, we plot v_B for the location of the central cavity (top panel) and, for comparison, also that of the outer regular disk (bottom panel). Comparing the overall behavior of v_B in the cavity versus that in the disk, and focusing on the time interval t > 200 yr, we note that while in the latter we mainly see outward transport, in the former we observe significant fluctuations. With the exception of the h-a case, v_B is mainly outward but sometimes inward.

Between the three runs, the h-a seems to have the fastest outward magnetic flux transport – a fact that can also be deduced from the sudden flux evacuation visible in the bottom panel of Fig. 11. The average values of v_B are shown in Table 1. We see that overall v_B is higher in the disk – in the case of the h-aa and h-f the magnetic flux of the cavity is almost unchanged.

Fig. 6 Space-time plot of the surface density for the three runs (from top to bottom: h-aa, h-f , h-a).

Fig. 7 Comparison of the vertically averaged surface density (top) and plasma parameter βP (bottom panel) over radius for the three runs. The initial profile (blue) is common for all models, the other lines are time averaged.

Fig. 8 Radial profiles of the ℒ parameter for the two HE runs. The vertical dashed black and green lines indicate the density maxima of the ring in the h-a and h-aa case, respectively. The vertical gray dashed line indicates the outer rim of the cavity.

Fig. 9 Vertically integrated profiles of the density, vertical magnetic field, and radial Maxwell stress in the cavity. We note the anticorrelation between density and Bz . The black dashed line (on the left) indicates the density maxima, while the gray dashed line (on the right) marks the end of the cavity.

4 Discussion

The results presented in the previous section demonstrate that TDs with central “cavities” – that is, with a significant depression in the gas surface density – are robust and long-lived entities. Moreover, their characteristics remain broadly stable when evolving the non ideal MHD equations over hundreds of inner orbits. This is specifically true when focusing on the effects of OR and AD only (also cf. Martel & Lesur 2022). The picture is somewhat altered when including the HE – both in its aligned and anti-aligned field configurations. In these cases, we do observe a certain tendency for a radial modulation of various quantities – such as, for instance, the surface mass density, the vertical magnetic flux, as well as the local accretion rate and the wind mass-loss rate. The overall ability of the disk to accrete material efficiently through the low-surface density region is, however, not affected when including the HE. In the contrary, in the aligned configuration, h-a, the accretion rate even gets boosted by more than 50%.

Regarding the magnetization of the disk, Martel & Lesur (2022) concluded that the cavity self regulates to β_P ~ 1–10. This is in agreement with our Hall-free run. However, for the two HE cases, h-aa, and h-a, in the regions where we have the rings, β_P naturally increases, reaching β_P ≃ 100. Another difference with Martel & Lesur (2022) arises in regards to the mass accretion rate within the cavity: in their equivalent Hallfree run it is more than three times larger than ours – that is, (~1.4 × 10⁻⁷ compared to ~4 × 10⁻⁸ M_⊙ yr⁻¹). We believe this discrepancy arises from differences between our respective models. Despite efforts to maintain comparable initial parameters, notable distinctions remain. These include the implementation of the diffusivities of the non ideal effects, the inclusion of thermochemical heating, and an overall thinner disk. Nevertheless, both values fall well within the range of typical inferred accretion rates for various observed TDs (see e.g., Espaillat et al. 2014). This is also true for the accretion rates in the Hall runs. Even though the h-a (h-aa) has the highest (lowest) rate, the values range from (2−8) × 10⁻⁸ M_⊙ yr⁻¹ – this may potentially be responsible for some of the variation in the observational scatter in Ṁ_a . Additionally, in all three runs we observe transonic accretion at small radii. This is something that is expected and, in fact, has been noted in a number of previews works: there is no clear distinction in the accretion rate of a transition disk in comparison with a full disk (e.g., Fang et al. 2013; Manara et al. 2014). As Wang & Goodman (2017) point out, a steady accretion rate combined with a several magnitude drop in density should naturally lead to (trans)sonic accretion speeds.

Magnetic flux transport is a topic that has been extensively studied (e.g., Bai & Stone 2017; Gressel et al. 2020; Lesur 2021b; Martel & Lesur 2022), yet it remains puzzling. While differences in models and simulations across these studies certainly account for some of the scatter in reported values, there is still surprisingly persistent disagreement over the transport rates – as recently highlighted by Lesur (2021b). For instance, Bai & Stone (2017) found that the h-aa configuration has the highest magnetic velocity advection v_B, while in our case, it is the h-a configuration. Gressel et al. (2020) reported higher v_B values for the Hall-free run, whereas Martel & Lesur (2022) observed inward transport. Lesur (2021b) suggested that v_B is strongly dependent to the strength of the Ohmic diffusivity – and to a lesser extent to the ambipolar diffusivity. The different implementation and treatment of the diffusivities might indeed be the case for these discrepancies – or the different evolution of the system in the cases where the setup is otherwise the same (i.e., Gressel et al. 2020, who simulated non-TD PPDs).

Similarly puzzling is the emergence of annular rings in our HE simulations. Self-organization into annular rings or so-called zonal flows is something that has been reported consistently in MHD simulations with a net magnetic field over the last decade – see Riols & Lesur (2019) for a detailed discussion. However, the mechanisms behind the emergence of these structures seem to be diverse. Bai & Stone (2014) and Bai (2015) found that the structures observed in their models are a robust outcome of the magnetorotational instability (MRI) with a net vertical magnetic flux, and these structures are further amplified by ambipolar diffusion. Kunz & Lesur (2013), Béthune et al. (2016), and Krapp et al. (2018) attributed their rings to the HE and the confinement of B_ɀ from the radial Maxwell stress. Béthune et al. (2017) deducted that the cause is ambipolar diffusion, acting anti-diffusively locally, due to the opposite signs of J_ϕ and Jϕ_,⊥ ; Suriano et al. (2018, 2019) also attributed ring formation to ambipolar diffusion, though the magnetic structure that it develops in the disk.

In this work, we noticed that in the two HE runs, the accretion rate due to the radial Maxwell stress is actually larger than that from the MHD wind (i.e., Ṁ_Rϕ > Ṁ_Rɀ), in the under-dense region before the ring is formed. As has been mentioned before, this is very similar to what was reported in Suriano et al. (2017); Riols & Lesur (2019) – however we did not collect convincing evidence as to why this should only happen in the HE runs.

Fig. 10 Accretion rate in the cavity due to the wind (Ṁɀϕ) and the radial stress (ṀRϕ – see Eq. (9)) for the three runs.

Fig. 11 Space-time evolution of the poloidal magnetic flux function ψ(θ = π/2) for the three runs (h-aa, h-f, h-a).

Fig. 12 Magnetic flux velocity in the midplane over time.

5 Conclusions

We have performed 2D global axisymmetric R-MHD simulations of TDs, specifically including the HE and a thermochemistry module (Gressel et al. 2020) for a more realistic temperature treatment⁵, as well as dynamically evolving micro-physical diffusion parameters, η_O, η_H, and η_A.

Our paper highlights results from three runs with the same initial parameters – the only difference lies in the configuration of the Hall effect: we had the Hall anti-aligned effect (h-aa) where the HE acts to suppress any inward curvature of the magnetic field lines in the poloidal plane; the Hall free (only OR and AD, termed h-f ); and the Hall aligned, where the HE enhances any concave deformation of the field lines, thus favoring the amplification of laminar fields within the disk. Our main conclusions are as follows:

In all three runs we observe a sustained cavity that retains its size, and an outer disk that behaves like a laminar disk with MHD winds. Overall we observe outward magnetic flux transport, with the transport velocity within the cavity being somewhat lower than that of the outer disk. The cavity features a similar level of magnetization for all runs, with β_P ~ 1–10. As expected from previous simulations without a cavity (Sarafidou et al. 2024), an increasing trend in magnetization is observed in the runs h-aa to h-f to h-a.

This increasing order also holds for the measured accretion rates of the cavity, with Ṁ_a ~ 3 × 10⁻⁸, 5 × 10⁻⁸, 8 × 10⁻⁸ for the h-aa, h-f , h-a case, respectively. This is due to the significantly amplified horizontal magnetic field that is the result of the h-a configuration, that in turn leads to an increased in-plane Maxwell stress. While decomposing the stresses, we find that for all cases, the magnetic Maxwell stress is significantly higher than the purely hydrodynamic Reynolds stress. Specifically, it is the vertical Maxwell stress (i.e., the stress owing to the MHD wind) that contributes the most toward the mass accretion rate. This is a result that we have already seen in full disks and that comes to support the notion that TDs behave like full disks.

A direct result of the similar accretion rates in TDs and full disks is that the gas accretion velocity is expected to be close to sonic or transonic. Indeed, for all runs, the accretion velocity within the cavity is ~c_s for r ∈ (2.5, 5.5)au and even >c_s for r < 2.5 au. We also notice that the bulk of the accretion comes from an accumulated radial mass flux (not unlike what has been termed a “streamer” in the literature) that originates from the inner radius of the outer disk (at r ~ 7.5 au) and permeates the cavity above or below the midplane, from a height of ~H/R = 0.2 (cf. also Martel & Lesur 2022, who speculate about a Rayleigh-Taylor origin of the avoidance of the midplane).

At the same time, we observe some level of self organization within the flow, where we have the formation of a ring in the cavity – but only in the case of the HE runs. We believe that this is a result of the MHD wind removing angular momentum from the gap more efficiently than in regions of high surface density – even though it is currently unclear why this did not occur in a comparable manner in the Hall free case. The formation of rings in global MHD disks with magnetic winds calls for more detailed study: It has been shown (see, e.g., Riols & Lesur 2018; Krapp et al. 2018) that the radial pressure profile resulting from the sub-Keplerian rotation of the zonal rings can trap dust. As a consequence, such features could well have implications for the accumulation of small solids – with observationally detectable features, for example, in the dust continuum emission, or more generally, for boosting planet forming mechanisms.

Acknowledgements

This work used the NIRVANA-III MHD code version 3.8, developed by Udo Ziegler at the Leibniz-Institut für Astrophysik Potsdam (AIP). Computations were performed on the TAURUS node of the local AIP cluster. We acknowledge support of the DFG (German Research Foundation), Research Unit ‘Transition discs’, funding ID 325594231. This work was co-funded⁶ by the European Union (ERC-CoG, Epoch-of-Taurus, No. 101043302).

Appendix A Assessment of the cavity longevity

Our simulations may not have run for long enough to unambiguously confirm that they have fully reached a steady state – that is, with respect to the sustained coexistence of an inner cavity and an outer disk over extended timescales. This limitation stems from the intrinsic complexity of our model and the associated computational cost. However, to address this issue, we conducted a lower-resolution, Hall-free simulation (run over 12,000 years) to assess whether the cavity can be maintained over these timescales. In the following, we briefly present the main results of this fiducial study.

Fig. A.1 Space-time plot of the vertically averaged density.

Figure A.1 presents the space-time plot of the vertically averaged density, which displays an overall steady state despite recurring high-frequency variations in the density. These results are consistent with Fig. 3 of Martel & Lesur (2022), who also observe significant density variability within the cavity. Given that their models evolved stably over longer timescales, there is no indication that our models would behave differently.

Fig. A.2 Time-averaged accretion rate.

Figure A.2 shows the accretion rate as a function of radius, averaged over the last 5,000 years of that simulation. Notably, the radial profile of the accretion rate shows a dip (with negative values) just inside the original transition radius, while it features a peak outside of that radius. The corresponding Fig. 8 in Martel & Lesur (2022) shows a similar non constant profile of the accretion rate with radius, depending on the vertical extent of the applied averaging. The radial modulation in the accretion rate suggests that the inner disk edge recedes outward at a slow speed. The space-time plot further supports this notion, revealing that indeed the cavity is expanding outward ever so slightly. Initially, this trend is also seen in Martel & Lesur (2022), where nevertheless a stationary state is obtained at even longer integration times.

The primary objective of our present work is to explore cavity dynamics and investigate the transonic accretion velocities in that region. Given the added complexity from including the Hall effect, we argue that our 33 orbits of integration time at the transition radius provide a sufficiently long timescale for this analysis to be meaningful.

Appendix B Non ideal effects in the disk

Here we showcase the non ideal MHD profiles for the cavity and the outer disk. We chose the h-aa model and we plot the Ambipolar Elsasser number ${\mathrm{\Lambda }}_{\text{A}}=v_{rm\text{\hspace{0.17em}A}}^2/\left(\mathrm{\Omega }{\eta }_{\text{A}}\right)$, the magnetic Reynolds number Rm = ΩH²/η_O and Lundquist Hall number ℒ_H = v_AH/η_H to demonstrate that indeed, our disk is dominated by ambipolar diffusion (see Fig. B.1 below). At the same time, it is important to remember that the Hall effect is dispersive, not diffusive, and can be dynamically important even in areas that are AD dominated. This is shown in Fig 8, where we see that indeed in the cavity we expect the “strong Hall regime” to be present.

Fig. B.1 Time-averaged and vertically integrated profiles of the ambipolar Elsasser number, the magnetic Reynolds number, and Lundquist Hall number. We chose these specific dimensionless quantities as they are not dependent on the magnetic field strength and thus provide a better indicator of the relative importance of each of the non ideal effects – also see Lesur (2021a).

References

All Tables

All Figures

Fig. 1 Time-averaged snapshots of the gas density (top row), plasma parameter β (middle row), and toroidal magnetic field (bottom row), for: the initial model (first column), as well as the three runs with Hall anti-aligned (h-aa, second column), Hall free (h-f , third column), and Hall aligned (h-a, fourth column). The streamlines in the first row indicate the poloidal velocity vectors. The blue lines in the middle row are placed where β = 0.1 indicating the vertical integration limits ±zd. Finally, the contour lines in the bottom-row panels highlight the poloidal magnetic field lines. The initial radial location of the cavity edge is marked by a blue diamond in all panels.

In the text

	Fig. 2 Radially averaged accretion rates through the cavity region over time. Shaded areas indicate uncertainties due to the vertical variability.
In the text

Fig. 3 Time-averaged accretion rates and stress terms versus radius for the three runs. In the first panel, we show the measured accretion rates, in the middle panel we show the Reynolds (dashed) vs. the Maxwell stress (solid line), normalized by the vertically integrated pressure. In the bottom panel, we decompose the normalized Maxwell stress into its radial (dashed) and surface (solid line) components.

In the text

	Fig. 4 Vertical profiles of the magnetic field components along a poloidal field line (each highlighted by a thick blue line in the bottom panels of Fig. 1, above) threading the cavity.
In the text

	Fig. 5 Time-averaged snapshot of the flow Mach number over 100yr (upper panels) and normalized accretion velocity profiles (lower panels) for the three runs. On the snapshots we are overlaying the radial mass flux using a gray alpha gradient colormap. The black dashed lines indicate the disk midplane as well as the disk surface.
In the text

	Fig. 6 Space-time plot of the surface density for the three runs (from top to bottom: h-aa, h-f , h-a).
In the text

	Fig. 7 Comparison of the vertically averaged surface density (top) and plasma parameter βP (bottom panel) over radius for the three runs. The initial profile (blue) is common for all models, the other lines are time averaged.
In the text

	Fig. 8 Radial profiles of the ℒ parameter for the two HE runs. The vertical dashed black and green lines indicate the density maxima of the ring in the h-a and h-aa case, respectively. The vertical gray dashed line indicates the outer rim of the cavity.
In the text

	Fig. 9 Vertically integrated profiles of the density, vertical magnetic field, and radial Maxwell stress in the cavity. We note the anticorrelation between density and Bz . The black dashed line (on the left) indicates the density maxima, while the gray dashed line (on the right) marks the end of the cavity.
In the text

	Fig. 10 Accretion rate in the cavity due to the wind (Ṁɀϕ) and the radial stress (ṀRϕ – see Eq. (9)) for the three runs.
In the text

	Fig. 11 Space-time evolution of the poloidal magnetic flux function ψ(θ = π/2) for the three runs (h-aa, h-f, h-a).
In the text

	Fig. 12 Magnetic flux velocity in the midplane over time.
In the text

	Fig. A.1 Space-time plot of the vertically averaged density.
In the text

	Fig. A.2 Time-averaged accretion rate.
In the text

	Fig. B.1 Time-averaged and vertically integrated profiles of the ambipolar Elsasser number, the magnetic Reynolds number, and Lundquist Hall number. We chose these specific dimensionless quantities as they are not dependent on the magnetic field strength and thus provide a better indicator of the relative importance of each of the non ideal effects – also see Lesur (2021a).
In the text