© The Authors 2025

1. Introduction

Protoplanetary disks (PPDs) do not evolve in isolation from their environment. Among other parameters, magnetic fields play a crucial role in understanding the evolution of PPDs and can be present in two regions of a disk: (i) strong stellar magnetic fields from the forming protostar, which is most likely partly or fully convective, interacting with the innermost disk, and (ii) large-scale magnetic fields threading the whole disk (e.g. Ohashi et al. 2018; Tang et al. 2023), whose origin is the subject of current research. Either the field is generated inside the disk by some dynamo effect (e.g. Reyes-Ruiz et al. 2004; Brandenburg et al. 1995) or the field is already present during the collapse of the protostellar cloud core (e.g. Dudorov & Khaibrakhmanov 2014). The stellar magnetic field as well as the large-scale disk field strongly influence the dynamics of an evolving PPD on a macroscopic scale (e.g. Ferreira et al. 2006; Bai & Stone 2013; Simon et al. 2013). For instance, angular momentum can be transferred between a protostar and a disk through torques of the stellar field acting on the inner disk parts. The details of such a coupled star-disk system are complex and involve the inclusion of many parameters, for example, the abundance and efficiency of accretion-powered stellar winds (APSW; Gehrig & Vorobyov 2023). Additionally, the position of a magnetically truncated inner disk radius (e.g. Hartmann et al. 2016; Steiner et al. 2021) or the launch of magnetocentrifugal winds (MCWs; Blandford & Payne 1982) are also dependent on the large-scale field. Apart from the macroscopic influence, the development of the magneto-rotational instability (MRI), which in some disk regions is believed to play a non-negligible role in driving accretion (Balbus & Hawley 1991), relies on the abundance of an initial large-scale seed field. The strength of the large-scale field is also believed to be important for determining the efficiency of MRI in PPDs (e.g. Bai & Stone 2011). One approach to investigate how MRI is altered by a magnetic field is to assume a maximally efficient MRI in a steady-state disk, which then determines the magnetic field strength for every radius in the disk (e.g. Mohanty et al. 2018; Delage et al. 2022).

The dipole field strengths of classical T Tauri stars (CTTS) range between 500 G and 1.5 kG (e.g. Johns-Krull et al. 1999; Johnstone et al. 2014) and are therefore strong enough to disrupt the disk at a certain radius, the so-called magnetic truncation radius, r_trunc (Hartmann et al. 2016). Although the details of magnetospheric accretion remain very complex to investigate (e.g. Romanova et al. 2002), it has been shown by Bessolaz et al. (2008) that under the assumption of a dipolar stellar magnetic field, the disk is magnetically truncated very close to the equipartition of magnetic pressure and thermal gas pressure (for detailed simulations of the inner disk movement during disk evolution, e.g. Steiner et al. 2021). Actual observations of large-scale magnetic fields threading a PPD remain challenging; however, magnetic fields have been detected on virtually every observationally accessible scale, including fields abundant in molecular clouds (e.g. Crutcher & Kemball 2019), which indicate that large-scale magnetic fields are likely to be present in a newly forming protostellar core and the resulting PPD. Recent observations of outflows and jets originating from the disk (e.g. Lee 2020) as well as polarisation mapping of the surroundings of T Tauri stars (e.g. Bertrang et al. 2017; Brauer et al. 2017) further strengthen the assumption of large-scale magnetic fields being abundant in PPDs. Additionally, the observation and high-precision measurement of the Zeeman effect in various emission lines has led to the calculation of an upper limit for some PPDs of the order of ∼10 mG for both the poloidal and toroidal fields at r ≈ 40 AU (e.g. Lankhaar & Teague 2023).

The interaction between the disk and the large-scale field is complex. Early 1+1D studies found that an initially purely vertical fossil field threading a thin disk is not dragged inwards efficiently due to outward-directed diffusion counterbalancing inward-directed magnetic advection very easily (e.g. Lubow et al. 1994). Taking the vertical velocity profile of the disk into account has modified this result such that a large-scale field can also be dragged inwards efficiently in a thin disk due to a larger radial mass transport and higher ionisation near the disk surface. The corresponding magnetic resistivity is then adequately averaged in the vertical direction to represent this enhanced magnetic flux transport Guilet & Ogilvie (2012, 2014). However, the source of resistivity is still a subject of investigation. Earlier studies have used turbulent resistivity to describe the non-ideal effects (e.g. Lubow et al. 1994; Guilet & Ogilvie 2012, 2014), but recent studies strongly indicate that PPDs are laminar almost everywhere but the innermost disk region (e.g. Dudorov & Khaibrakhmanov 2014; Mohanty et al. 2018; Delage et al. 2022). Laminar non-ideal diffusion can be caused by Ohmic diffusion (OD), ambipolar diffusion (AD), or the Hall effect (HE).

Apart from 1+1D models, global long-term evolution models have been conducted, including magnetic fields in 2+1D (e.g. Vorobyov et al. 2020). Such models have to deal with the numerical difficulties of conducting magnetohydrodynamic (MHD) simulations over an extended period of time. In particular, the incorporation of magnetic diffusion poses challenges to the numerical stability of a scheme. Furthermore, such models are limited with respect to the position of the inner disk boundary, since the Courant-Friedrichs-Levy (CFL) condition (Courant et al. 1928) severely limits the time step achievable in a disk simulation inside of ∼1 AU. The model of Vorobyov et al. (2020) incorporates large-scale magnetic fields but without considering non-ideal magnetic diffusion due to numerical challenges, which leads to unusually strong magnetic advection and field amplification in their models.

The previous models used to study the influence of a large-scale magnetic field fall short in the following cases: (i) modelling the long-term simulations of a disk and large-scale magnetic field including the inner disk; (ii) modelling the effects of a stellar dipole together with a large-scale field on the inner disk and the topology, especially in the inner disk; (iii) describing the effects of a MCW on the long-term evolution of a PPD. In this work, we aim to investigate the evolution of a large-scale field towards a stationary state for a disk truncated by a stellar dipole. It has been shown that the inner disk treatment is crucial in terms of long-term evolution (e.g. Vorobyov et al. 2019; Steiner et al. 2021; Gehrig et al. 2022). Therefore, a combined treatment of both stellar- and disk components has the potential to alter the disk properties in the inner disk, compared to the common assumption of a torque-free inner disk boundary. A self-consistent evolution of the disk and large-scale magnetic field offers the possibility to study the disk field topology and strength not through assuming maximum efficient MRI (e.g. Mohanty et al. 2018; Delage et al. 2022) but through inward-dragging and amplification of a more physically motivated background fossil field. Furthermore, our simulations are not confined to steady-state models and allowed us to follow the disk field evolution over time. In this work, we introduce our new model and focus on the validation of our method by comparison with current work (e.g. Dudorov & Khaibrakhmanov 2014; Mohanty et al. 2018; Delage et al. 2022). Therefore, we focus on steady-state models. By combining the stellar and disk field components, we can study the effects of various parameters on the large-scale disk field, such as the stellar rotation period, magnetic field strength, X-ray luminosity, and the accretion rate (Sect. 3). These results form the basis for investigating MCWs self-consistently with a hydrodynamic accretion disk in the subsequent studies of this series.

2. Model description

The construction of a long-term evolution simulation of a PPD including a stellar magnetic field and a large-scale magnetic field involves a model describing the hydrodynamic flow of the disk while also covering the temperature structure of the disk and the influence of magnetic fields (cf. Sect. 2.1).

2.1. Basic equations

We use TAPIR for our 1+1D simulations (Ragossnig et al. 2020; Steiner et al. 2021), which solves the hydrodynamic equations using a fully implicit time integration scheme. Hence, the CFL condition does not restrict the timestep¹. This allowed us to also include the very inner disk in our long-term hydrodynamic simulations, which is necessary to incorporate the torque of the stellar dipole and pressure gradients in our model (Steiner et al. 2021). The vertical disk extent is assumed to be hydrostatic, allowing for a vertical integration of the hydrodynamic equations Eqs. (1)–(3) (analogously to e.g. Vorobyov et al. 2020; Steiner et al. 2021). The effects of a magnetic field, B, threading the disk is included by adding the non-ideal induction equation Eq. (4),

$$\begin{array}{cc}& \frac{\partial \mathrm{\Sigma }}{\partial t}+\mathrm{\nabla }·(\mathrm{\Sigma }\mathit{u})=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\partial (\mathrm{\Sigma }\mathit{u})}{\partial t}+\mathrm{\nabla }·(\mathrm{\Sigma }\mathit{u}:\mathit{u})-\frac{B_{\mathrm{z}}\mathit{B}}{2\pi }\hfill \\ \hfill & +\mathrm{\nabla }{P}_{\mathrm{gas}}+\mathrm{\nabla }·Q+\mathrm{\Sigma }\mathrm{\nabla }\psi +H_{\mathrm{p}}\mathrm{\nabla }\left(\frac{B_{\mathrm{z}}^2}{4\pi }\right)=0,\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial (\mathrm{\Sigma }e)}{\partial t}+\mathrm{\nabla }·(\mathrm{\Sigma }\mathit{u}e)+P_{\mathrm{gas}}\mathrm{\nabla }·\mathit{u}\hfill \\ \hfill & +Q:\mathrm{\nabla }\mathit{u}-4\pi \mathrm{\Sigma }{\kappa }_{\mathrm{R}}(J-S)+{\dot{E}}_{\mathrm{rad}}=0,\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\partial \mathit{B}}{\partial t}-\mathrm{\nabla }\times (\mathit{u}\times \mathit{B})+\mathrm{\nabla }\times ({\eta }_{\mathrm{O}}\mathit{J})\hfill \\ \hfill & -\mathrm{\nabla }\times (\frac{{\eta }_{\mathrm{A}}}{{|\mathit{B}|}^2}\mathit{J}\times \mathit{B}\times \mathit{B})=0,\hfill \end{array}$$(4)

where Σ, e, and u ≡ (u_r, u_φ, 0)^T denote current density, specific inner energy density and the gas velocity restricted to the planar direction. We used an ideal equation of state for which the gas pressure can be written as P_gas = Σ e (γ − 1) with an adiabatic coefficient γ = 7/5 (e.g. Vorobyov et al. 2020) and which allows for the definition of the pressure scale height, H_p. The terms including the viscous pressure tensor, Q, describe the viscous torque enabling radial angular momentum transport and the viscous heating in Eq. (2) and Eq. (3), respectively. The disk is heated and cooled by stellar irradiation and radiative cooling through its surface on both sides of the disk, respectively. Additionally, vertical radiative transport is included in the diffusion limit by assuming local thermal equilibrium (LTE) in vertical direction. Heating, cooling and vertical radiative transport are described by the net heating- and cooling rate per unit surface in Eq. (3), Ė_rad. In our model radiative transport in radial direction is also included approximately in the diffusion limit, which is described by the term Σ κ_R (J−S), where J and S denote the zeroth moment of the radiation field and the source function, respectively (for a detailed explanation see Ragossnig et al. 2020).

The viscous-α was implemented analogously to Steiner et al. (2021). We adopted the layered-viscosity model of Gammie (1996), which reads as

$$\begin{array}{c}\hfill \alpha ={\alpha }_{\mathrm{base}}+{\alpha }_{\mathrm{surf}}+{\alpha }_{\mathrm{deep}},\end{array}$$(5)

where α_base, α_surf, and α_deep denote the base viscosity in the dead zone (DZ), the viscosity due to active MRI in the disk’s surface layer, and the α-viscosity in the deeper disk layers. The surface layer is defined on each side of the disk with a thickness corresponding to a column density of Σ₀ = 100 g cm² (measured from the disk surface inwards). If Σ > Σ₀, non-thermal ionisation cannot penetrate deeper into the disk, and a deep layer with a different viscosity α_deep is assumed:

$$\begin{array}{cc}\hfill {\alpha }_{\mathrm{surf}}(r)& ={\alpha }_{\mathrm{MRI},\mathrm{eff}}[s_{\mathrm{\Sigma }}\frac{{\mathrm{\Sigma }}_0}{\mathrm{\Sigma }(r)}+(1-s_{\mathrm{\Sigma }})],\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill {\alpha }_{\mathrm{deep}}(r)& ={\alpha }_{\mathrm{MRI},\mathrm{eff}}{s}_{\mathrm{\Sigma }}(1-\frac{{\mathrm{\Sigma }}_0}{\mathrm{\Sigma }(r)})s\left(T_0\right),\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill s_{\mathrm{\Sigma }}& =\{\begin{array}{cc}1\hfill & \forall \mathrm{\Sigma }\ge {\mathrm{\Sigma }}_0\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array},\hfill \end{array}$$(8)

where α_MRI, eff denotes the effective α due to MRI. For all radii r > r_DZ, out (where r_DZ, out denotes the outer DZ boundary), MRI operates with reduced efficiency due to AD (e.g. Delage et al. 2022). We modelled this effect by applying a reduced α_AD < α_MRI in the region r > r_DZ.out (the values for α_MRI, α_AD, and α_base used in our reference model are listed in Table 2):

$$\begin{array}{cc}\hfill {\alpha }_{\mathrm{MRI},\mathrm{eff}}& =\{\begin{array}{cc}{\alpha }_{\mathrm{MRI}}\hfill & \forall r<r_{\mathrm{DZ},\mathrm{out}}\hfill \\ {\alpha }_{\mathrm{AD}}\hfill & \forall r\ge r_{\mathrm{DZ},\mathrm{out}}\hfill \end{array}.\hfill \end{array}$$(9)

The factor s(T₀) determines if the deep layer becomes MRI-active (Steiner et al. 2021, e.g.), which would lead to episodic accretion events (e.g. Steiner et al. 2021; Cecil et al. 2024). In this work the development of episodic accretion is suppressed by setting s(T₀) = 0, which is due to the focus on stationary disks in this work (a stationary disk does not have any time-dependent features such as episodic accretion).

The net heating- or cooling rate, Ė_rad, is written as follows,

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{rad}}=\mathrm{\Gamma }-\mathrm{\Lambda },\end{array}$$(10)

where Γ and Λ describe the heating and cooling rate per unit surface, respectively. The heating rate is a combination of both viscous heating and stellar irradiation on both sides of the disk surface. Viscous heating of the gas happens mostly close to the midplane, where ρ is largest. The corresponding inner energy e is transported vertically up to the disk surface, yielding a temperature, T_surf, at the disk surface, if LTE and an optically thick disk is assumed (τ ≫ 1). With the additional assumption of the source function S being constant over the vertical extent of the disk, T_surf can be written as (e.g. Dong et al. 2016; Ragossnig et al. 2020)

$$\begin{array}{c}\hfill T_{\mathrm{surf}}^4=\frac{8{\tau }_{\mathrm{P}}{\sigma }_{\mathrm{SB}}}{1+2{\tau }_{\mathrm{P}}+\frac{3}{2}{\tau }_{\mathrm{R}}{\tau }_{\mathrm{P}}}(T_{\mathrm{gas},0}^4+T_{\mathrm{amb}}^4+\frac{L_{\star }}{4\pi r^2}{f}_{\mathrm{irr}}sin{\alpha }_{\mathrm{irr}}).\end{array}$$(11)

Here, ${\tau }_{\mathrm{R}}={\kappa }_{\mathrm{R}}{\rho }_0{H}_{\mathrm{p}}={\kappa }_{\mathrm{R}}\mathrm{\Sigma }/\sqrt{2}$ and ${\tau }_{\mathrm{P}}={\kappa }_{\mathrm{P}}\mathrm{\Sigma }/\sqrt{2}$ describe the Rosseland and Planck mean opacity, respectively. κ_R has two contributions: (i) a gaseous component κ_R, gas using the method of Ferguson et al. (2005) with the protostellar abundances of Asplund et al. (2021), and (ii) a dust-dominated component for low temperatures κ_R, dust (T_gas < 500 K). Then κ_R = κ_R, gas + f_dust κ_R, dust, where f_dust denotes the gas-to-dust ratio (Pollack et al. 1985). The Planck mean opacity is approximated with κ_P = 2.4 κ_R and is based on Nakamoto & Nakagawa (1994), Hueso & Guillot (2005). The Stefan-Boltzmann constant is denoted as σ_SB, and T_amb is the ambient temperature of the interstellar medium (ISM). The parameters f_irr and α_irr describe which percentage of the stellar irradiation is reflected at the surface and the impact angle of the irradiation at the surface (with respect to the midplane). We note that we did not include the heating of the disk surface due to stellar X-ray irradiation (cf., Güdel 2015, and also Sect. 4.4). Then Γ and Λ can be written as (e.g. Vorobyov et al. 2020)

$$\begin{array}{cc}\hfill \mathrm{\Gamma }& =\frac{8{\tau }_{\mathrm{P}}{\sigma }_{\mathrm{SB}}}{1+2{\tau }_{\mathrm{P}}+\frac{3}{2}{\tau }_{\mathrm{R}}{\tau }_{\mathrm{P}}}(\frac{L_{\star }}{4\pi r^2}{f}_{\mathrm{irr}}sin{\alpha }_{\mathrm{irr}}+T_{\mathrm{amb}}^4),\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill \mathrm{\Lambda }& =\frac{8{\tau }_{\mathrm{P}}{\sigma }_{\mathrm{SB}}}{1+2{\tau }_{\mathrm{P}}+\frac{3}{2}{\tau }_{\mathrm{R}}{\tau }_{\mathrm{P}}}{T}_{\mathrm{gas},0}^4.\hfill \end{array}$$(13)

The stellar magnetic field at the disk surface (which in our case was chosen to be the magnetic surface z_B, cf. Guilet & Ogilvie 2012, explained in detail in Sect. 2.2) was modeled as a dipole with the components in cylindrical coordinates B_⋆ = (B_r, ⋆, B_φ, ⋆, B_z, ⋆)^T (analogously to Steiner et al. 2021; Gehrig & Vorobyov 2023, but without ignoring the radial component of the dipole):

$$\begin{array}{cc}\hfill B_{\mathrm{z}},\star (r,z)& =B_{\star }{R}_{\star }^3\frac{r^2-2z^2}{{\left(\sqrt{r^2+z^2}\right)}^5},\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill B_{\mathrm{r}},\star (r,z)& =-B_{\star }{R}_{\star }^3\frac{3rz}{{\left(\sqrt{r^2+z^2}\right)}^5},\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill B_{\phi },\star & \simeq -{\alpha }_{\mathrm{cor}}{B}_{\mathrm{z}},\star (r,z)[1-{\left(\frac{{\mathrm{\Omega }}_{\star }}{\mathrm{\Omega }(r)}\right)}^{{\alpha }_{\mathrm{cor}}}],\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill {\alpha }_{\mathrm{cor}}& \equiv \{\begin{array}{cc}+1\hfill & \text{if}r<r_{\mathrm{cor}}\hfill \\ -1\hfill & \text{if}r>r_{\mathrm{cor}}\hfill \end{array},\hfill \end{array}$$(17)

where B_⋆ denotes the magnetic field strength at the stellar surface and R_⋆ the stellar radius. α_cor is a parameter which ensures that max(|B_φ, ⋆|) ≃ |B_z, ⋆|. This occurs very close to the star and at a certain distance outside the corotation radius, r_cor, and is physically motivated by the assumption that a stronger winding of the dipole would lead to buoyancy and reconnection of B_φ, ⋆ (e.g. Steiner et al. 2021; Gehrig & Vorobyov 2023).

2.2. Large-scale magnetic disk field

In this study we made the following assumptions about the large-scale disk magnetic field: (i) The magnetic field is formally averaged in azimuthal direction and temporally over the dynamical timescale at a certain radius, which means that the large-scale field in our model is axisymmetric. Additionally, the averaging results in only the macroscopic net field being considered and the microscopic fluctuations, which eventually lead to turbulence and to MRI, being modelled in the usual way using the isotropic α-viscosity (Shakura & Sunyaev 1973) and the turbulent resistivity, η_t (e.g. Guilet & Ogilvie 2012, 2014). (ii) We did not include a toroidal field component, B_φ, in the disk exterior (above the disk surface) in our model to avoid the complexity of starting a MHD wind. Hence, B_φ^s2 was set to zero by means of a Dirichlet boundary condition (i.e., B_φ^s = 0; for an in-depth analysis of the possible boundary conditions for B_φ^s cf. Khaibrakhmanov & Dudorov 2022). B_φ^s acts as an outer boundary at the disk surface for a toroidal field inside the disk. This means that B_φ inside the disk can be non-vanishing due to shearing of B_r caused by differential rotation, as long as it is connected to the exterior field at the disk surface (cf. Sect. 2.4 for a more detailed treatment). This causes additional diffusion and therefore contributes to the (vertically averaged) magnetic flux transport velocity, ${\overline{\overline{u}}}_{\psi }$ (e.g. Guilet & Ogilvie 2012; Leung & Ogilvie 2019, and Eq. (51) below). (iii) The magnetic field was assumed to be stationary in vertical direction, compared to its slow evolution in radial direction. This is due to the assumption of a geometrically thin disk (i.e., H_p/r ≪ 1; e.g. Guilet & Ogilvie 2012, for a discussion of the difference between viscous and vertical dynamical timescales in terms of an asymptotic expansion of the basic MHD equations, cf. also Sect. 2.7). (iv) No dynamo effect, which might develop in the disk (e.g. Reyes-Ruiz et al. 2004), was considered in this work. The poloidal, large-scale magnetic field, B = (B_r, 0, B_z)^T, can be represented by the magnetic flux ψ,

$$\begin{array}{c}\hfill \mathit{B}=\mathrm{\nabla }\psi \times {\mathit{e}}_{\phi },\end{array}$$(18)

where e_φ describes the azimuthal unit vector in cylindrical coordinates (e.g. Guilet & Ogilvie 2014). The induction equation Eq. (4) can then be written as follows (e.g. Ogilvie & Livio 2001),

$$\begin{array}{cc}\hfill \frac{\partial \psi (r,z)}{\partial t}+\mathit{u}·\mathrm{\nabla }\psi (r,z)& =\eta (r,z)r^2\mathrm{\nabla }·\left(\frac{1}{r^2}\mathrm{\nabla }\psi (r,z)\right),\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill \eta (r,z)& ={\eta }_{\mathrm{O}}(r,z)+{\eta }_{\mathrm{A}}(r,z)+{\eta }_{\mathrm{t}}(r,z),\hfill \end{array}$$(20)

where η is the combination of all resistivity mechanisms considered in this work, i.e. turbulent resistivity, OD, AD, which are denoted by η_t, η_O, and η_A, respectively. We note that η_A in Eq. (20) only denotes the ‘Ohm-like’ contribution. AD also has a ‘Hall-like’ contribution, which was neglected in this work due to a vanishing toroidal field in the disk exterior above the surface. This in turn causes the ‘Hall-like’ contribution to vanish (e.g. Leung & Ogilvie 2019). Assuming a thin disk allowed ψ(r, z) to be split into a radial contribution, ψ₀(r), and a small contribution, which describes the vertical dependence, ψ₁(r, z), and for which |ψ₁|≪|ψ₀| holds,

$$\begin{array}{cc}\hfill \psi (r,z)& ={\psi }_0(r)+{\psi }_1(r,z),\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill B_{\mathrm{z}}& =\frac{1}{r}\frac{\partial {\psi }_0}{\partial r}=B_{\mathrm{z},\mathrm{d}}+B_{\mathrm{z}},\star ,\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill B_{\mathrm{r}}& =-\frac{1}{r}\frac{\partial {\psi }_1}{\partial z}=B_{\mathrm{r},\mathrm{d}}+B_{\mathrm{r}},\star ,\hfill \end{array}$$(23)

where B_r and B_z denote the combined stellar and large-scale fields in radial and vertical direction, respectively. Our simulations were done in 1+1D; however Eq. (19) is height-dependent and can be evaluated at any z above the midplane, which makes it necessary to average Eq. (19) vertically. Analogously to Ogilvie & Livio (2001), Guilet & Ogilvie (2014), a sensible choice is a conductivity-weighted vertical average up to the magnetic surface, z_B (for the remainder of this work also referenced as disk surface), which is the height above the disk where magnetic and thermal gas pressure are equal (or equivalently, where the plasma beta β(z = z_B) = 1). In the case of a hydrostatically stratified disk z_B reads as

$$\begin{array}{cc}\hfill z_{\mathrm{B}}& =H_{\mathrm{p}}\sqrt{ln\left(\frac{2}{\pi }{\beta }_0\right)},\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill {\beta }_0& \equiv \frac{8\pi P_{\mathrm{gas}}}{B_{\mathrm{z}}^2},\hfill \end{array}$$(25)

where β₀ and z_B describe the plasma beta at the disk midplane and the magnetic surface, respectively. Averaging this way results in $\overline{\overline{\eta }}$ being dominated by vertical disk layers of high conductivity (low resistivity) and takes into account that the magnetic field is transported inwards more effectively in vertical disk layers of high conductivity. Applying the conductivity-weighted vertical averaging to η and the radial gas velocity, u_r, yields

$$\begin{array}{cc}\hfill \frac{1}{\overline{\overline{\eta }}(r)}& \equiv \frac{1}{z_{\mathrm{B}}}{\displaystyle {\int }_0^{z_{\mathrm{B}}}}\frac{1}{\eta (r,z)}\mathrm{dz}\hfill \\ \hfill & =\frac{1}{z_{\mathrm{B}}}{\displaystyle {\int }_0^{z_{\mathrm{B}}}}\frac{1}{{\eta }_{\mathrm{O}}(r,z)+{\eta }_{\mathrm{A}}(r,z)+{\eta }_{\mathrm{t}}(r,z)}\mathrm{dz},\hfill \end{array}$$(26)

$$\begin{array}{cc}\hfill {\overline{\overline{u}}}_{\mathrm{r}}(r)& \equiv \frac{\overline{\overline{\eta }}(r)}{z_{\mathrm{B}}}{\displaystyle {\int }_0^{z_{\mathrm{B}}}}\frac{u_{\mathrm{r}}(r,z)}{\eta (r,z)}\mathrm{dz},\hfill \end{array}$$(27)

where the double-bar notation denotes conductivity-weighted averaged quantities. Integrating Eq. (19) vertically between [ − z_B, z_B] and applying Eqs. (22)–(23) as well as using Eqs. (26)–(27), we can rewrite the magnetic flux transport equation Eq. (19) in a height-independent form (strictly following Guilet & Ogilvie 2014),

$$\begin{array}{cc}& \frac{\partial {\psi }_0(r)}{\partial t}+r{\overline{\overline{u}}}_{\psi }{B}_{\mathrm{z}}=0,\hfill \end{array}$$(28)

$$\begin{array}{cc}& {\overline{\overline{u}}}_{\psi }={\overline{\overline{u}}}_{\mathrm{r}}+\frac{\overline{\overline{\eta }}}{z_{\mathrm{B}}}\frac{1}{B_{\mathrm{z}}}(B_{\mathrm{r}}(r,z=z_{\mathrm{B}}^{-})-z_{\mathrm{B}}\frac{\partial B_{\mathrm{z}}}{\partial r}),\hfill \end{array}$$(29)

where z_B⁻ denotes the vertical height at the surface, but approached from inside the disk (analogously z_B⁺ denotes the magnetic surface approached from above the disk). This is important for modelling the connection at the disk surface between an external magnetic field and the large-scale field inside the disk Sect. 2.4 (e.g. Guilet & Ogilvie 2014). The vertical averaging performed in Eq. (29) has no description of the vertical dependency of B_r inside the disk, but only its value at the surface B_r(r, z = z_B⁻). The determination of this external field is not trivial and is described in Sect. 2.3. Additionally, averaging ${\overline{\overline{u}}}_{\mathrm{r}}$ and $\overline{\overline{\eta }}$ requires knowing the vertical profiles of u_r(r, z) and η(r, z). Assuming β ≫ 1 inside the disk and a constant viscous-α in vertical direction allows us to find an analytic expression for the vertical velocity profiles (e.g. Guilet & Ogilvie 2012; Leung & Ogilvie 2019), which we briefly summarise in Sect. 2.4). For the vertical profile of η we need the ionisation fraction, x_e, which can be caused by thermal ionisation or by various non-thermal ionisation sources such as X-rays and cosmic rays. We discuss all ionisation sources included in our model in Sect. 2.5.

2.3. External magnetic field

The radial magnetic field at the disk surface, B_r^s, is connected to the exterior field above and below the disk, which is a superposition of three components: (i) a fossil field, B_∞, formally created by currents at infinity, (ii) the stellar dipole field, B_⋆, and (iii) the field created by purely toroidal currents (due to our assumption of a poloidal field), J_φ, in the disk. Due to the thin-disk approximation used in our model we only need to describe the radial part of the large-scale field’s flux function, ψ₀ (cf. Eq. (28)). The three magnetic flux contributions from the fossil field, stellar dipole, and current density, J_φ, in the disk are denoted by ψ_∞, ψ_⋆, and ψ_d, respectively,

$$\begin{array}{cc}\hfill {\psi }_0& ={\psi }_{\infty }+{\psi }_{\star }+{\psi }_{\mathrm{d}}.\hfill \end{array}$$(30)

The Biot-Savart law was used to find an expression for the large-scale magnetic field and hence for the corresponding magnetic flux, ψ_d, of an axisymmetric current distribution, J_φ^s, by integrating the field caused by the current loops at every r over the full radial extent of the disk (e.g. Jackson 1999; Ogilvie 1997):

$$\begin{array}{cc}\hfill {\psi }_{\mathrm{d}}(r)& =\frac{1}{\pi }{\displaystyle {\int }_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}}\frac{rr\prime }{r+r\prime }\left[\frac{(2-k^2)K(k)-2E(k)}{k^2}\right]J_{\phi }^{\mathrm{s}}(dr\prime )dr\prime ,\hfill \end{array}$$(31)

$$\begin{array}{cc}\hfill k^2& =\frac{4rr\prime }{{(r+r\prime )}^2},\hfill \end{array}$$(32)

$$\begin{array}{cc}\hfill E(k)& ={\displaystyle {\int }_0^{\pi /2}}{(1-k^2{sin}^2(x))}^{1/2}dx,\hfill \end{array}$$(33)

$$\begin{array}{cc}\hfill K(k)& ={\displaystyle {\int }_0^{\pi /2}}{(1-k^2{sin}^2(x))}^{-1/2}dx,\hfill \end{array}$$(34)

where E(k) is the complete elliptic integral of the first kind and K(k) corresponds to the complete elliptic integral of the second kind (e.g. Guilet & Ogilvie 2014). Applying Stokes’ theorem to Eq. (18) allows the disk’s radial field at the surface, B_r, d^s, to be connected with J_φ^s at the disk surface (e.g. Ogilvie 1997; Guilet & Ogilvie 2014):

$$\begin{array}{cc}\hfill B_{\mathrm{r},\mathrm{d}}^{\mathrm{s}}(r)& =2J_{\phi }^{\mathrm{s}}(r).\hfill \end{array}$$(35)

We note that B_r, d^s describes the radial field corresponding to the currents in the disk only, without any contributions from stellar dipole or fossil field.

The magnetic fluxes ψ_∞ and ψ_⋆ are calculated as follows,

$$\begin{array}{cc}\hfill {\psi }_{\infty }& =\frac{r^2-R_{\star }^2}{2}{B}_{\infty },\hfill \end{array}$$(36)

$$\begin{array}{cc}\hfill {\psi }_{\star ,0}& =B_{\star }{R}_{\star }^3(\frac{1}{R_{\star }}-\frac{1}{r}),\hfill \end{array}$$(37)

which both can be readily calculated, since the stellar dipole field, B_⋆, and the fossil field, B_∞, are both parameters of our model. Since Eq. (28) describes the evolution of ψ₀ at a certain time, we can use Eq. (35) with Eq. (31) to find a relation between ψ_d and B_r, d^s and then reorganise Eq. (30) as follows,

$$\begin{array}{cc}\hfill {\psi }_{\mathrm{d}}& ={\psi }_0-{\psi }_{\infty }-{\psi }_{\star }.\hfill \end{array}$$(38)

Finally, for calculating B_r, d^s we interpreted the integral of Eq. (31) as a linear operator, ℒ, which can then be formally inverted to obtain

$$\begin{array}{c}\hfill B_{\mathrm{r},\mathrm{d}}^{\mathrm{s}}={\mathcal{L}}^{-1}{\psi }_{\mathrm{d}}.\end{array}$$(39)

Numerically, ℒ⁻¹ is represented as a matrix describing the relationship between all discretised values of B_r, d^s and ψ_d. Its construction and inversion is non-trivial and numerically expensive and is briefly discussed in Appendix A.

2.4. Analytical vertical velocity and magnetic field profile

In 1+1D, hydrodynamical models usually employ a density-weighted, vertically averaged radial and rotational velocity, ${\overline{u}}_{\mathrm{r}}(r)$ and ${\overline{u}}_{\phi }(r)$, respectively, and therefore the actual velocity structure is usually not resolved by assuming a vertically constant velocity u = (u_r(r),u_φ(r),0)^T. However, work done by Guilet & Ogilvie (2012), Leung & Ogilvie (2019) indicates that radial magnetic flux transport can significantly deviate from mass transport, if the vertical velocity structure is included. This is due to the two different vertical averaging mechanisms for mass and magnetic flux transport (cf. Sect. 2.2 and Eqs. (26)–(27)). Using that the vertical dynamical timescale is much shorter than the radial diffusion and viscous timescales, the disk structure in vertical direction appears stationary with respect to viscous evolution, and therefore is treated as such (e.g. Guilet & Ogilvie 2012; Dudorov & Khaibrakhmanov 2014; Leung & Ogilvie 2019).

Strictly following the work of Guilet & Ogilvie (2014), Leung & Ogilvie (2019), the disk can be split into two zones vertically: (i) the regions within the magnetic surface, z_B, where it is assumed that the magnetic field is passive, i.e. β ≫ 1, and (ii) a magnetically dominated region outside z_B, where in the absence of outflows the magnetic field can be modelled as force-free. It follows that in the case of a vertically constant α and with a passive large-scale field, the gas velocity profiles are parabolic in the vertical direction,

$$\begin{array}{cc}\hfill u_{\mathrm{r}}& =u_{\mathrm{r},0}(r)+u_{\mathrm{r},2}(r)\frac{z^2}{H_{\mathrm{p}}^2},\hfill \end{array}$$(40)

$$\begin{array}{cc}\hfill u_{\phi }& =u_{\phi },0(r)+u_{\phi },2(r)\frac{z^2}{H_{\mathrm{p}}^2},\hfill \end{array}$$(41)

where u_r, 0 and u_φ, 0 denote the radial and toroidal velocities at the disk midplane, respectively, whereas u_r, 2 and u_φ, 2 are functions which describe the vertical dependency of u_r and u_φ. After integration of the vertical, stationary disk profile (e.g. Guilet & Ogilvie 2012), u_r, 2 and u_φ, 2 take the following form,

$$\begin{array}{cc}\hfill u_{\mathrm{r},2}& =\frac{\alpha c_{\mathrm{s}}}{1+4{\alpha }^2}(\frac{3H_{\mathrm{p}}}{r}-5\frac{\partial H_{\mathrm{p}}}{\partial r}),\hfill \end{array}$$(42)

$$\begin{array}{cc}\hfill u_{\phi },2& =\frac{c_{\mathrm{s}}}{2}(\frac{\partial H_{\mathrm{p}}}{\partial r}-\frac{3H_{\mathrm{p}}}{2r})+\frac{{\alpha }^2{c}_{\mathrm{s}}}{1+4{\alpha }^2}(\frac{3H_{\mathrm{p}}}{r}-5\frac{\partial H_{\mathrm{p}}}{\partial r}).\hfill \end{array}$$(43)

We note that, strictly speaking, Eq. (40) and Eq. (41) can only be written as above for a vertically constant α. It was shown that a more realistic vertical profile of α, which is certainly expected in a real-world PPD, has only a minor influence on the large-scale magnetic field in steady-state (e.g. Guilet & Ogilvie 2014). Therefore, the use of those simplified profiles in this work is a viable approximation. However, the diffusion timescale at which the magnetic field evolves may be altered significantly (e.g. Leung & Ogilvie 2019), which makes it necessary to adopt a more realistic vertical velocity profile for future work for cases when the detailed time evolution matters (e.g. the time-dependent long-term disk evolution with an outflow).

The magnetic field profiles B_r(r, z) and B_φ(r, z) are dependent on the values of the exterior field components B_r^s and B_φ^s at the disk surface. As mentioned above, to avoid the complications of a magnetically induced outflow, no toroidal field at the disk surface was assumed, i.e. B_φ^s = 0. However, inside the disk (0 ≤ z ≤ z_B) a toroidal component can develop (e.g. Guilet & Ogilvie 2012). The two-zone model briefly described above has the disadvantage that the intermediate region β ≈ 1 between the passive region β ≫ 1 and the magnetically dominated disk corona β ≪ 1 is infinitesimally small, which leads to a jump between the values of B_r(r, z) and B_φ(r, z) just below (z = z_B⁻) and above (z = z_B⁺) the disk surface,

$$\begin{array}{cc}\hfill B_{\mathrm{r}}(z=z_{\mathrm{B}}^{-})& =B_{\mathrm{r}}(z=z_{\mathrm{B}}^{+})+\mathrm{\Delta }{B}_{\mathrm{r}},\hfill \end{array}$$(44)

$$\begin{array}{cc}\hfill B_{\phi }(z=z_{\mathrm{B}}^{-})& =B_{\phi }(z=z_{\mathrm{B}}^{+})+\mathrm{\Delta }{B}_{\phi },\hfill \end{array}$$(45)

where ΔB_r and ΔB_φ are the jump conditions for the radial and toroidal field, respectively. The boundary conditions for B_r and B_φ and jump conditions are read as follows (e.g. Guilet & Ogilvie 2012, for an in-depth derivation),

$$\begin{array}{cc}\hfill B_{\mathrm{r}}(z_{\mathrm{B}}^{+})& =B_{\mathrm{r}}^{\mathrm{s}}+H_{\mathrm{p}}\frac{\partial B_{\mathrm{z}}}{\partial r},\hfill \end{array}$$(46)

$$\begin{array}{cc}\hfill \mathrm{\Delta }{B}_{\mathrm{r}}& =-H_{\mathrm{p}}\pi {\frac{\partial B_{\phi }}{\partial z}|}_{z=z_{\mathrm{B}}^{-}},\hfill \end{array}$$(47)

$$\begin{array}{cc}\hfill B_{\phi }(z_{\mathrm{B}}^{+})& =0,\hfill \end{array}$$(48)

$$\begin{array}{cc}\hfill \mathrm{\Delta }{B}_{\phi }& =0.\hfill \end{array}$$(49)

We note that the vertical profile of B_φ(r, z) is indeed important to obtain the jump condition Eq. (47) for B_r. With Eqs. (46)–(49), the final form of the vertical profiles for B_r and B_φ reads as (Guilet & Ogilvie 2012; Leung & Ogilvie 2019, but using OD, ‘Ohm-like’ AD, ignoring HE, and including hydrodynamical terms)

$$\begin{array}{cc}\hfill B_{\mathrm{r}}(r,z<z_{\mathrm{B}})& =\frac{z}{1-\pi c_{\mathsf{s}}{H}_{\mathrm{p}}/(2{\eta }_{{\mathrm{z}}_{\mathrm{B}}})}[\frac{B_{\mathrm{r}}^{\mathrm{s}}}{z_{\mathrm{B}}}+\frac{\partial B_{\mathrm{z}}}{\partial r}\hfill \\ \hfill & +(5-\frac{3\pi c_{\mathsf{s}}{H}_{\mathrm{p}}}{2{\eta }_{{\mathrm{z}}_{\mathrm{B}}}})\frac{B_{\mathrm{z}}{u}_{\mathrm{r},2}{z}_{\mathrm{B}}^2}{H_{\mathrm{p}}^{2}15{\eta }_{{\mathrm{z}}_{\mathrm{B}}}}-\frac{4B_{\mathrm{z}}\pi u_{\phi },2}{3{\eta }_{{\mathrm{z}}_{\mathrm{B}}}}]\hfill \\ \hfill & -\frac{B_{\mathrm{z}}{u}_{\mathrm{r},2}{z}^3}{3{\eta }_{{\mathrm{z}}_{\mathrm{B}}}{H}_{\mathrm{p}}^2},\hfill \end{array}$$(50)

$$\begin{array}{cc}\hfill B_{\phi }(r,z<z_{\mathrm{B}})& =\frac{c_{\mathsf{s}}}{2H_{\mathrm{p}}{\eta }_{{\mathrm{z}}_{\mathrm{B}}}[1-\pi c_{\mathsf{s}}{H}_{\mathrm{p}}/(2{\eta }_{{\mathrm{z}}_{\mathrm{B}}})]}\hfill \\ \hfill & \times [\frac{1}{2}(\frac{B_{\mathrm{r}}^{\mathrm{s}}}{z_{\mathrm{B}}}+\frac{\partial B_{\mathrm{z}}}{\partial r})z(z^2-z_{\mathrm{B}}^2)\hfill \\ \hfill & +\frac{B_{\mathrm{z}}{u}_{\mathbf{r}},2z(z^2-z_{\mathrm{B}}^2)}{120{\eta }_{{\mathrm{z}}_{\mathrm{B}}}^2{H}_{\mathrm{p}}^2}[z_{\mathrm{B}}^2(3c_{\mathsf{s}}{H}_{\mathrm{p}}\pi -14{\eta }_{{\mathrm{z}}_{\mathrm{B}}})\hfill \\ \hfill & -(z^2(6c_{\mathsf{s}}{H}_{\mathrm{p}}\pi -20{\eta }_{{\mathrm{z}}_{\mathrm{B}}})]\hfill \\ \hfill & -\frac{2B_{\mathrm{z}}{u}_{\phi },2z(z_{\mathrm{B}}^2-z^2)}{3c_{\mathsf{s}}{H}_{\mathrm{p}}{\eta }_{{\mathrm{z}}_{\mathrm{B}}}}(2{\eta }_{{\mathrm{z}}_{\mathrm{B}}}-c_{\mathsf{s}}{H}_{\mathrm{p}})],\hfill \end{array}$$(51)

where η_zB = η(z = z_B). Taking the value of η at z = z_B for B_r and B_φ is justified both by taking the value of η at the transition zone (z ≈ z_B) and by the fact that the marginally stable mode allowing a single bending of the field line in radial direction is localised around the transition zone as well (e.g. Ogilvie & Livio 2001; Guilet & Ogilvie 2012; Leung & Ogilvie 2019, for a thorough analysis of the marginal stability theory). Finally, the magnetic transport velocity can be calculated by evaluating B_r at z = z_B⁻ and by using Eq. (29),

$$\begin{array}{cc}\hfill {\overline{\overline{u}}}_{\psi }& ={\overline{\overline{u}}}_{\mathrm{r}}+\frac{\overline{\overline{\eta }}}{B_{\mathrm{z}}}\frac{\pi c_{\mathsf{s}}{H}_{\mathrm{p}}/(2{\eta }_{{\mathrm{z}}_{\mathrm{B}}})}{1-\pi c_{\mathsf{s}}{H}_{\mathrm{p}}/(2{\eta }_{{\mathrm{z}}_{\mathrm{B}}})}\frac{\partial B_{\mathrm{z}}}{\partial r}\hfill \\ \hfill & +\frac{1}{1-\pi c_{\mathsf{s}}{H}_{\mathrm{p}}/(2{\eta }_{{\mathrm{z}}_{\mathrm{B}}})}[\frac{\overline{\overline{\eta }}}{z_{\mathrm{B}}}\frac{B_{\mathrm{r}}^{\mathrm{s}}}{B_{\mathrm{z}}}+(5-\frac{3\pi c_{\mathsf{s}}{H}_{\mathrm{p}}}{2{\eta }_{{\mathrm{z}}_{\mathrm{B}}}})\hfill \\ \hfill & \times \frac{z_{\mathrm{B}}^2}{H_{\mathrm{p}}^2}\frac{u_{\mathrm{r},2}\overline{\overline{\eta }}}{15{\eta }_{{\mathrm{z}}_{\mathrm{B}}}}-\frac{4\pi u_{\phi },2\overline{\overline{\eta }}}{3{\eta }_{{\mathrm{z}}_{\mathrm{B}}}}]-\frac{u_{\mathrm{r},2}{z}_{\mathrm{B}}^2\overline{\overline{\eta }}}{3H_{\mathrm{p}}^2{\eta }_{{\mathrm{z}}_{\mathrm{B}}}}·\hfill \end{array}$$(52)

We note that $\overline{\overline{\eta }}$ and η_zB are appearing multiple times in Eq. (52), which emphasises again the need for a detailed model of η(r, z) (see Sect. 2.5).

For solving Eqs. (1)–(3), the usual density-weighted vertical averages were used for the unknowns Σ, e, and u (e.g. Ragossnig et al. 2020). Therefore, using the vertical velocity profiles Eqs. (40)–(41), ${\overline{u}}_{\mathrm{r}(\phi )}$³ were obtained as follows,

$$\begin{array}{cc}\hfill {\overline{u}}_{\mathrm{r}}(\phi )& =\frac{1}{\mathrm{\Sigma }}{\displaystyle {\int }_{-\infty }^{\infty }}{u}_{\mathrm{r}}(\phi )(r,z)\rho (r,z)dz=u_{\mathrm{r}}(\phi ),0+u_{\mathrm{r}}(\phi ),2,\hfill \end{array}$$(53)

Since u_r(φ),2 can be calculated using Eqs. (42)–(43), the remaining unknowns with respect to u are u_r(φ),0.

2.5. Non-ideal resistivity model

Magnetic resistivity, η(r, z), occurs due to turbulence, in our case denoted by η_t, or laminar, non-ideal diffusion processes such as OD, η_O; AD, η_A; and diffusion due to HE, η_H. However, η_H was not considered in this work. η_t is connected to the turbulence causing the MRI, and therefore to α. The laminar diffusion processes are dependent on the ionisation fraction, x_e, and the vertical temperature structure, T_gas(r, z), and are discussed in this section.

2.5.1. Turbulent resistivity:

The assumption of η_t being caused by the same turbulence as MRI and therefore turbulent viscosity, ν, leads to the expectation that the strengths of η_t and ν are of the same order. This is described by the magnetic Prandtl number, 𝒫, which is defined as

$$\begin{array}{cc}\hfill \mathcal{P}& =\frac{\nu }{{\eta }_{\mathrm{t}}},\hfill \end{array}$$(54)

where 𝒫 is fixed in this work to the value 𝒫 = 1, which leads to

$$\begin{array}{cc}\hfill {\eta }_{\mathrm{t}}& =\nu .\hfill \end{array}$$(55)

2.5.2. Laminar resistivity:

Both η_O and η_A depend on x_e in the disk,

$$\begin{array}{cc}\hfill x_{\mathrm{e}}& =\frac{n_{\mathrm{e}}}{n_{\mathrm{tot}}}=\frac{n_{\mathrm{e}}}{n_{\mathrm{e}}+n_{\mathrm{n}}+n_{\mathrm{i}}}\approx \frac{n_{\mathrm{e}}}{n_{\mathrm{n}}},\hfill \end{array}$$(56)

where n_tot, n_e, n_n, and n_i denote the total number density, and the number densities of electrons, neutrals, and ions, respectively. By far the most abundant element in PPDs is molecular hydrogen. Therefore, the assumption of n_H2 ≈ n_n ≈ n_tot is valid (e.g. Mohanty et al. 2018). We note that η_O is defined as

$$\begin{array}{cc}\hfill {\eta }_{\mathrm{O}}& =\frac{c^2}{4\pi {\sigma }_{\mathrm{O}}},\hfill \end{array}$$(57)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{O}}& \approx \frac{e^2{n}_{\mathrm{e}}}{m_{\mathrm{e}}{⟨\sigma v⟩}_{\mathrm{en}}}=\frac{e^2{x}_{\mathrm{e}}{n}_{\mathrm{tot}}}{m_{\mathrm{e}}{⟨\sigma v⟩}_{\mathrm{en}}},\hfill \end{array}$$(58)

where σ_O denotes the Ohmic conductivity and is directly proportional to x_e (e.g. Bai & Goodman 2009). The rate coefficient of electron-neutral collisions and electron mass are described by ⟨σv⟩_en and m_e, respectively.

We note that η_A has an additional dependency on both the magnetic field strength and Σ (e.g. Lesur et al. 2014) and can be written as (again for the assumption of electrons and single-charged ions only)

$$\begin{array}{cc}\hfill {\eta }_{\mathrm{A}}& ={\beta }_{\mathrm{i}}{\beta }_{\mathrm{e}}{\eta }_{\mathrm{O}},\hfill \end{array}$$(59)

$$\begin{array}{cc}\hfill {\beta }_{\mathrm{i}(\mathrm{e})}& =\frac{eB}{m_{\mathrm{i}(\mathrm{e})}c}\frac{1}{{\gamma }_{\mathrm{i}(\mathrm{e})}{\rho }_{\mathrm{n}}}·\hfill \end{array}$$(60)

Here β_i(e) denote the Hall parameters (e.g. Mohanty et al. 2018) for single-charged ions (index ‘i’) and electrons (index ‘e’). The variables m_i, ρ_n, and γ_i(e) identify the mean ion mass, density of neutral gas and the drag coefficient for ions or electrons, respectively. With the following definitions of γ_i(e), ⟨σv⟩_en, and the rate coefficient of ion-neutral collisions, ⟨σv⟩_in (e.g. Wardle 2007; Mohanty et al. 2018),

$$\begin{array}{cc}\hfill {\gamma }_{\mathrm{i}(\mathrm{e})}& =\frac{{⟨\sigma v⟩}_{\mathrm{in}(\mathrm{en})}}{m_{\mathrm{i}(\mathrm{e})}+m_{\mathrm{n}}},\hfill \end{array}$$(61)

$$\begin{array}{cc}\hfill {⟨\sigma v⟩}_{\mathrm{in}}& =1.6\times {10}^{-9}{\mathrm{cm}}^3{s}^{-1},\hfill \end{array}$$(62)

$$\begin{array}{cc}\hfill {⟨\sigma v⟩}_{\mathrm{en}}& ={10}^{-15}\left(\frac{128k_{\mathrm{B}}{T}_{\mathrm{e}}}{9\pi m_{\mathrm{e}}}\right){\mathrm{cm}}^3{s}^{-1},\hfill \end{array}$$(63)

the laminar resistivities, η_O and η_A, can be written as (assuming that the electron temperature T_e = T_gas):

$$\begin{array}{cc}\hfill {\eta }_{\mathrm{O}}(r,z)& \approx 230\sqrt{T_{\mathrm{gas}}(r,z)}\frac{1}{x_{\mathrm{e}}(r,z)}{\mathrm{cm}}^2{s}^{-1},\hfill \end{array}$$(64)

$$\begin{array}{cc}\hfill {\eta }_{\mathrm{A}}(r,z)& \approx 4.95\times {10}^{-15}\frac{{|\mathit{B}|}^2}{\rho {(r,z)}^2{x}_{\mathrm{e}}(r,z)}{\mathrm{cm}}^2{s}^{-1}.\hfill \end{array}$$(65)

We point out that we have explicitly stated the radial and vertical dependency of T_gas, x_e, ρ, η_O, and η_A to emphasise that it is necessary to model these dependences in a 1+1D model in an approximate way. Gas in a PPD is either ionised thermally (depicted by x_e, th) or non-thermally (x_e, nth),

$$\begin{array}{cc}\hfill x_{\mathrm{e}}& =x_{\mathrm{e},\mathrm{th}}+x_{\mathrm{e},\mathrm{nth}}.\hfill \end{array}$$(66)

2.5.3. Thermal ionisation:

A high gas temperature, T_gas(r, z), causes ionising collisions between atoms, which are counteracted by recombination of those same atoms with electrons. This eventually leads to an equilibrium level of ionisation and is described by the Saha equation (e.g. Mohanty et al. 2018), in the following shown for potassium (subscript K)

$$\begin{array}{cc}\hfill \frac{n_{\mathrm{e}}{n}_{+,\mathrm{K}}}{n_{0,\mathrm{K}}}& =2.41\times {10}^{15}\sqrt{T_{\mathrm{gas}}^3}exp(-\frac{5.04\times {10}^4\mathrm{K}}{T_{\mathrm{gas}}}){\mathrm{cm}}^{-3}\hfill \end{array}$$(67)

$$\begin{array}{cc}& ={\mathcal{S}}_{\mathrm{K}}(T_{\mathrm{gas}}){\mathrm{cm}}^{-3},\hfill \end{array}$$(68)

where n_0, K and n_+, K are the number densities of neutral and ionised potassium, respectively. The function 𝒮_K(T_gas) denotes the right-hand side of Eq. (67) and is dependent on T_gas only. Following Mohanty et al. (2018), the thermal ionisation fraction, x_e, th, can then be written as,

$$\begin{array}{cc}\hfill x_{\mathrm{e},\mathrm{th}}& =\frac{-1+\sqrt{1+4x_{\mathrm{K}}(n_{\mathrm{tot}}/{\mathcal{S}}_{\mathrm{K}}(T_{\mathrm{gas}}))}}{2n_{\mathrm{tot}}/{\mathcal{S}}_{\mathrm{K}}(T_{\mathrm{gas}})},\hfill \end{array}$$(69)

where x_K is the ionisation fraction of potassium and is set to x_K = 1.97 × 10⁻⁷ (e.g. Mohanty et al. 2018). The formulation of x_e, th as in Eq. (69) covers an important edge case: In a stratified disk atmosphere, density significantly drops in the surface layers (n_tot → 0) which leads to a high x_e (cf. Eq. (56)). In the case of an equilibrium ionisation this eventually leads to x_e, th ≈ x_K.

For properly describing thermal ionisation, not only the radial gas temperature profile, but also an approximation in vertical direction T_gas(r, z) is needed. At the midplane ρ(r, z = 0) = ρ₀(r) is largest, which leads to effective viscous heating at z = 0. The resulting thermal radiation is transported radially and vertically through the disk according to Eq. (3). Additionally, the disk is both heated by stellar radiation on its surface (on both sides) and also cools radiatively across its surface (cf. Eq. (11)). Assuming LTE and using the Eddington approximation, the vertical temperature distribution can be approximated in the optical thick case (τ ≫ 1) as (e.g. Hubeny 1990),

$$\begin{array}{cc}\hfill T_{\mathrm{gas}}^4(r,z)& =T_{\mathrm{gas},0}^4-\frac{{\tau }_{\mathrm{R}}^{-}(r,z)}{{\tau }_{\mathrm{R}}}(T_{\mathrm{gas},0}^4-T_{\mathrm{surf}}^4),\hfill \end{array}$$(70)

where ${\tau }_{\mathrm{R}}^{-}(r,z)\equiv {\kappa }_{\mathrm{R}}\mathrm{\Sigma }\mathrm{erf}(z/(\sqrt{2}{H}_{\text{p}}))/\sqrt{2\pi }$ and T_gas, 0 are the optical depth from the midplane up to the height z and the midplane gas temperature, respectively. With Eq. (70) x_e can be calculated in vertical direction (cf. Eq. (69)).

Thermal ionisation is very likely to only be significant in the inner disk (r ≲ 0.2 AU), which has two reasons: (i) viscous heating is more effective in the inner regions with large Σ, as long as the disk is optically thick, and (ii) in regions very close to the star, the disk becomes optically thin and the irradiation from the star suffices for the gas to get thermally ionised, i.e. after heating by absorption.

2.5.4. Non-thermal ionisation:

The main sources for non-thermal ionisation are galactic cosmic rays and stellar X-ray as well as FUV radiation. Additionally, radioactive decay (RA) plays a role in some regions close to the disk midplane. The cosmic-ray ionisation rate, ξ_CR, is adapted from Lesur et al. (2014),

$$\begin{array}{cc}\hfill {\xi }_{\mathrm{CR}}& ={10}^{-16}exp\{-\frac{\mathrm{\Sigma }}{9.6\times {10}^2\mathrm{g}c{m}^{-2}}\}{\mathrm{s}}^{-1},\hfill \end{array}$$(71)

whereas the X-ray ionisation rate, ξ_XR, can be modelled according to a fit done by Bai & Goodman (2009), Delage et al. (2022), using an X-ray temperature of T_XR = 3 keV. X-ray photons can either directly ionise the gas, ξ_XR, dir, or indirectly ionise it through scattering, ξ_XR, sca,

$$\begin{array}{cc}\hfill {\xi }_{\mathrm{XR},\mathrm{dir}}& =\frac{L_{\mathrm{XR}}}{{10}^{29}\mathrm{erg}{\mathrm{s}}^{-1}}{\left(\frac{r}{1\mathrm{AU}}\right)}^{-2.2}\hfill \end{array}$$(72)

$$\begin{array}{cc}& {\xi }_{\mathrm{XR},\mathrm{dir},0}[exp\{-{\left(\frac{{\mathrm{\Sigma }}^{+}}{{\mathrm{\Sigma }}_{\mathrm{XR},\mathrm{dir}}}\right)}^{0.4}\}+exp\{-{\left(\frac{{\mathrm{\Sigma }}^{-}}{{\mathrm{\Sigma }}_{\mathrm{XR},\mathrm{dir}}}\right)}^{0.4}\}],\hfill \end{array}$$(73)

$$\begin{array}{cc}\hfill {\xi }_{\mathrm{XR},\mathrm{sca}}& =\frac{L_{\mathrm{XR}}}{{10}^{29}\mathrm{erg}{\mathrm{s}}^{-1}}{\left(\frac{r}{1\mathrm{AU}}\right)}^{-2.2}\hfill \\ \hfill & {\xi }_{\mathrm{XR},\mathrm{sca},0}[exp\{-{\left(\frac{{\mathrm{\Sigma }}^{+}}{{\mathrm{\Sigma }}_{\mathrm{XR},\mathrm{sca}}}\right)}^{0.65}\}+exp\{-{\left(\frac{{\mathrm{\Sigma }}^{-}}{{\mathrm{\Sigma }}_{\mathrm{XR},\mathrm{sca}}}\right)}^{0.65}\}],\hfill \end{array}$$(74)

where Σ_XR, dir = 3.5 × 10⁻³ g cm⁻² and Σ_XR, sca = 1.64 g cm⁻² are the penetration depths for direct and scattered X-ray ionisation (e.g. Delage et al. 2022), respectively. Σ⁺ = ∫_z^∞ρ(r, z) dz and Σ⁻ = Σ − Σ⁺ denote the column densities above and below a certain height z, respectively. Additionally, ξ_{XR, dir, 0} = 6 × 10⁻¹² s⁻¹ and ξ_{XR, sca, 0} = 10⁻¹⁵ s⁻¹ describe the undamped direct and scattered X-ray ionisation, respectively. Finally, a rather simple radioactive decay ionisation rate, ξ_RA, is adopted according to Lesur et al. (2014), ignoring radial variation and the effects of varying dust-to-gas ratio (e.g. Delage et al. 2022),

$$\begin{array}{cc}\hfill {\xi }_{\mathrm{RA}}& ={10}^{-19}{\mathrm{s}}^{-1}.\hfill \end{array}$$(75)

The effects of FUV radiation were neglected in this work, as the penetration depth of FUV is of order 10⁻² g cm⁻² and yields an ionisation front in the uppermost surface layer of a PPD. Therefore, FUV ionisation becomes significant when modelling MCWs or photoevaporative winds (e.g. Cecil et al. 2024), but is not likely to play an important role in modelling the resistivity for magnetic flux transport of a large-scale, purely poloidal field inside the disk. The total ionisation rate, ξ, was then written as

$$\begin{array}{cc}\hfill \xi & ={\xi }_{\mathrm{CR}}+{\xi }_{\mathrm{XR},\mathrm{dir}}+{\xi }_{\mathrm{XR},\mathrm{sca}}+{\xi }_{\mathrm{RA}}.\hfill \end{array}$$(76)

In an equilibrium state, which is necessary in assuming a steady-state in vertical direction, ionisation is balanced by radiative recombination, α_r, and recombination on dust particles, α_d, in the following way,

$$\begin{array}{cc}\hfill (1-x_{\mathrm{e},\mathrm{nth}})\xi & ={\alpha }_{\mathrm{r}}{x}_{\mathrm{e},\mathrm{nth}}^2{n}_{\mathrm{tot}}+{\alpha }_{\mathrm{d}}{x}_{\mathrm{e},\mathrm{nth}}{n}_{\mathrm{tot}},\hfill \end{array}$$(77)

where we assumed that the number density of neutrals is very close to the total number density, n_n ≈ n_tot, and that the temperature regimes of thermal and non-thermal ionisation are non-overlapping. This is a good approximation, because thermal ionisation becomes dominant above T_th, active ≳ 1000 K and has a very narrow transition temperature range. Below T_th, active non-thermal ionisation dominates. Hence, the non-thermal ionisation fraction, x_e, nth, can be calculated separately from x_e, th (e.g. Vorobyov et al. 2020). We note that dissociative recombination was not considered in this work due to the fact that radiative and dust recombination yield an upper and lower limit for x_e, nth, respectively (e.g. Dudorov & Khaibrakhmanov 2014), whereas dissociative recombination leads to an ionisation fraction between those limits.

The radiative recombination rate, α_r, can be approximated by (e.g. Spitzer 1978; Vorobyov et al. 2020)

$$\begin{array}{cc}\hfill {\alpha }_{\mathrm{r}}& =2.07\times {10}^{-11}{T}_{\mathrm{gas}}^{-1/2}{\mathrm{cm}}^3{s}^{-1},\hfill \end{array}$$(78)

whereas for the dust recombination rate, α_d, we adopt a parametrisation described by Dudorov & Khaibrakhmanov (2014), which is depicted in Table 1 and is valid for a grain size a_d = 0.1 μm and a dust-to-gas ratio f_dust = 0.01. Formally α_d is defined as

$$\begin{array}{cc}\hfill {\alpha }_{\mathrm{d}}& =X_{\mathrm{g}}⟨{\sigma }_{\mathrm{g}}{V}_{\mathrm{ig}}⟩,\hfill \end{array}$$(79)

where X_g, σ_g, and V_ig denote the number density fraction of dust grains, the grain cross-section, and ion-grain relative velocity (e.g. Dudorov & Khaibrakhmanov 2014). In the case of a fixed f_dust, the following relations hold,

$$\begin{array}{cc}\hfill X_{\mathrm{g}}& =\frac{n_{\mathrm{d}}}{n_{\mathrm{tot}}}\propto a_{\mathrm{d}}^{-3},\hfill \end{array}$$(80)

$$\begin{array}{cc}\hfill {\sigma }_{\mathrm{g}}& \propto a_{\mathrm{d}}^2,\hfill \end{array}$$(81)

which, according to Eq. (79), leads to the overall scaling of α_d:

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{d}}\propto a_{\mathrm{d}}^{-3}{a}_{\mathrm{d}}^2\propto a_{\mathrm{d}}^{-1}.\end{array}$$(82)

Therefore a different grain size, a_d, can be used by taking the values of Table 1 for a_d and then by scaling it according to Eq. (82) (e.g. Dudorov & Khaibrakhmanov 2014).

2.6. Solution procedure

Our 1+1D model aims to solve the hydrodynamical equations Eqs. (1)–(3) together with a poloidal large-scale field as described in Eq. (28) (cf. Sect. 2.1 and Sect. 2.2). The vertical structure of the disk is crucial in understanding the radial magnetic flux transport of the large-scale field. Hence, we have provided a description of the vertical velocity profile (analogously to Guilet & Ogilvie 2012; Leung & Ogilvie 2019) (cf. Sect. 2.4) and a model for the vertical ionisation fraction and resistivity (cf. Sect. 2.5. The large-scale field was assumed to be seeded by an external fossil field, which was taken to be force-free (in the absence of outflows) and which served as a boundary at the disk surface for the large-scale field in the disk interior (cf. Sect. 2.3). The vertical velocity and resistivity profiles were vertically averaged to fit into a 1+1D approach. The solving procedure followed the following steps:

1. Calculate the vertical velocity contributions, u_r, 2 and u_φ, 2 (see Eqs. (42)–(43)). The radial contributions u_r, 0 and u_φ, 0 were solved through the hydrodynamic simulation.

2. Determine the density-weighted vertical averages, ${\overline{u}}_{\mathrm{r}}$ and ${\overline{u}}_{\phi }$ (see Eq. (53)).

3. Obtain ψ_d by subtracting the magnetic flux contributions of stellar dipole and external field as in Eq. (38).

4. Calculate ℒ (cf. Eq. (31) and Eq. (39)) and inversion to obtain B_r^s (see Appendix A for details).

5. Determine x_e (see Eq. (66), Eq. (69) and Eq. (77)), the vertical temperature profile (see Eq. (70)) as well as η_O and η_A (see Eq. (64) and Eq. (65), respectively) for all radii at five different heights (z = 0, z = H_p, $z=H_{\mathrm{p}}+\frac{1}{3}(z_{\text{B}}-H_{\mathrm{p}})$, $z=H_{\mathrm{p}}+\frac{2}{3}(z_{\text{B}}-H_{\mathrm{p}})$, and z = z_B).

6. Calculate the conductivity-weighted vertically averaged $\overline{\overline{\eta }}(r)$ between −z_B ≤ z ≤ z_B. The necessary integration over η(r, z) = η_O(r, z)+η_A(r, z) is done by piece-wise numerical integration using the values for the heights z calculated in the preceding step (see Eq. (26)).

7. Determine the conductivity-weighted vertically averaged radial velocity ${\overline{\overline{u}}}_{\mathrm{r}}$ (see Eq. (27)) by using $\overline{\overline{\eta }}$.

8. Determine the magnetic flux transport velocity ${\overline{\overline{u}}}_{\psi }$ by using $\overline{\overline{\eta }}$, ${\overline{\overline{u}}}_{\mathrm{r}}$ and η(r, z = z_B) (see Eq. (52)).

9. Finally, solve the hydrodynamical equations Eqs. (1)–(3) and Eq. (28) with TAPIR (e.g. Ragossnig et al. 2020).

The most expensive operation (in the context of simulation time) is the inversion of the operator, ℒ. The details of this inversion and the steps taken to maintain adequate performance in solving the MHD equations are explained in Sect. A.

2.7. Validity of our approximations for modelling the vertical disk and magnetic field structure

We used a 1+1D disk model and assumed a geometrically thin disk (H_p/r ≪ 1), which usually means that the radial and toroidal velocity components, u_r(r) and u_φ(r), respectively, are taken to be constant in the vertical direction. However, describing the evolution of a large-scale magnetic field required modelling of magnetic advection and diffusion, which both are strongly dependent on the ionisation degree and velocity structure in vertical direction (e.g. Guilet & Ogilvie 2012; Leung & Ogilvie 2019; Delage et al. 2022, see also Sect. 2.2 and Sect. 2.4). Hence, we discuss the validity of the essential approximations for describing the vertical disk- and large-scale magnetic field structure. The key points are (for a detailed discussion cf. Appendix B):

• The focus on Class II PPDs justifies the use of the thin-disk approximation.

• A comparison of the viscous timescale, t_visc, and the dynamical timescale, t_dyn, justifies the use of stationary vertical velocity- and large-scale magnetic field profile (see also Sect. 3.2).

• The assumption of a dominant vertical magnetic field is valid for most of the disk, as long as the large-scale magnetic field remains weak, and the inclination angle does not exceed i > 30°. The validity is only questionable in the innermost disk region, where magnetic compression likely plays a role.

• The somewhat arbitrary choice of vertical layers for the calculation of x_e, as well as OD and AD in the vertical direction, is justified through a comparison with a profile with different choices of vertical layers. It can be concluded that our approach with 5 layers is a good compromise between numerical cost and reasonable accuracy.

• For describing the thermally ionised region in the inner disk a vertical temperature profile was used. It is compared to and in good agreement with 2D simulations of the vertical disk structure (e.g. Jankovic et al. 2021).

Not only the validity of our approximations, but also the comparison with former work is important. To this end we compared our reference model with the analytic profiles of Dudorov & Khaibrakhmanov (2014) and found very good agreement with their radial magnetic field profiles (see Fig. 3 in Sect. 3.3). We have also dedicated Sect. 4.1 for comparing our results with the work of various authors. We found that our results are consistent with the models of Mohanty et al. (2018) and Delage et al. (2022). Especially the calculation of the ionisation degree and resistivity profiles are in very good agreement with the work of those authors. The work of Jankovic et al. (2021) also agrees well with our model of the inner disk and the vertical temperature structure and thermally ionised region. The sum of all those comparisons strengthens the confidence that our approximate disk model is capable of reproducing the results from previous studies (e.g. Dudorov & Khaibrakhmanov 2014; Guilet & Ogilvie 2014) but also has the ability to provide additional scientific insights. We nevertheless discuss remaining fundamental limitations of our model methodology in Sect. 4.4. Some of these can be overcome in future work. However, we do not consider them to be of crucial importance in the framework of our present study that focuses on the investigation of the poloidal large-scale magnetic field topology threading a stationary disk.

3. Results

We split our results into three main categories: (i) the evolution of a background fossil field towards a steady state (see Sect. 3.1), (ii) the in-depth study of the solution obtained after the large-scale magnetic field in the disk has become stationary (see Sect. 3.3), and (iii) a parameter study of how the resulting stationary system consisting of a disk, star, and large-scale magnetic field changes for different parameters (see Sects. 3.4–4.3).

3.1. Evolution of the large-scale magnetic field towards a stationary state

For validation of the large-scale magnetic field structure obtained in this paper, it is useful to compare it to similar, already available work done on this topic (e.g. Guilet & Ogilvie 2014; Mohanty et al. 2018; Delage et al. 2022), which use a stationary magnetic field topology. A (magneto)hydrodynamic treatment of the large-scale field requires evolving the field and disk sufficiently long to eventually reach a stationary state. This is done with TAPIR in two stages for our reference model (see Table 2):

• (i)

  The first stage consists of evolving the disk with a prescribed stellar field and a prescribed large-scale disk field, analogously to Steiner et al. (2021). The stationary disk is obtained by self-consistently calculating the position of the inner disk boundary, r_in, taking into account the accretion rate, Ṁ(r_in), and the initial stellar magnetic field strength, B_⋆, 0. Additionally, we start with a radially constant, purely vertical, large-scale disk magnetic field, B_d, 0(r) = B_∞ (see Eq. (36) in Sect. 2.2), threading the disk at time t = 0. This field is superimposed with the stellar magnetic field to provide an initial star-disk magnetic field threading the disk, B₀ = B_⋆, 0 + B_d, 0 (see Sect. 2.2), and is kept constant during the first stage. This essentially means that only the disk ‘feels’ the impact of the magnetic field, whereas the magnetic field dragging due to radial mass transport is not considered during this first stage.

• (ii)

  The second stage then relaxes the constraints of a temporally constant large-scale field by allowing for the evolution of the magnetic field from the initial configuration B₀ towards a stationary field topology by interacting with the accretion disk (using the parameters of Table 2). In this phase, the disk and magnetic field are evolving as a coupled system. Starting with the stationary disk obtained from the first stage, the magnetic field is allowed to advect inward, caused by radial mass transport through the disk. Due to non-ideal MHD effects (in this study turbulent resistivity, OD, and AD are considered, see Sect. 2.5) the magnetic field advection is counteracted by diffusion, which at some point in time balances the inward-directed advection velocity, u_adv, which then yields a stationary magnetic field.

In Fig. 1, we present the large-scale magnetic field evolution with a focus on the timescales until it reaches its final stationary state. The relevant timescale describing the radial magnetic flux transport is the diffusion timescale, $t_{\mathrm{diff}}\equiv r^2/\overline{\overline{\eta }}$. As a reference, we adopt the value of t_diff, stat, corresponding to the stationary large-scale field, at the outer DZ boundary, r_DZ, out. In the very inner disk (r ≲ 0.1 AU), the stellar dipole field is dominating the large-scale field, which is represented by a low ${\overline{\overline{u}}}_{\psi }$ and a relatively high B in the inner disk for the whole disk evolution. Interestingly, despite a reduced radial mass transport, Ṁ, in the DZ (between the two vertical white dashed lines in Fig. 1), the large-scale field in the DZ is effectively advected inwards⁴. This is due to significant magnetic flux transport in the upper disk layers, which is taken into account in our model by the conductivity-weighted vertical averaged resistivity, $\overline{\overline{\eta }}$ (see Sect. 3.3), and results in a considerably larger ${\overline{\overline{u}}}_{\psi }$ compared to ${\overline{u}}_{\mathrm{r}}$ in the DZ. The sharp decrease of t_diff, stat in the vicinity r_DZ, out, together with the inward-directed magnetic field transport, leads to accumulation of magnetic flux, ψ, early on during magnetic field evolution (starting at t ∼ 10⁻¹ t_diff, stat(r_DZ, out)). This causes a steep radial gradient in ψ and therefore an increase in B_z, d just inside the outer DZ boundary, revealing a locally increased field strength, |B| (cf. colour-coding of Fig. 1, around r_DZ, out). In the outer disk, the large-scale field is not transported as effectively as in the inner region, which is due to AD being the dominant non-ideal MHD effect over the whole vertical extent in the outer disk. Eventually, this yields a larger $\overline{\overline{\eta }}$ in the outer disk and therefore a smaller ${\overline{\overline{u}}}_{\psi }$. For t ≳ 10 t_diff, stat(r_DZ, out) the large-scale field becomes essentially stationary everywhere in the disk. This means that vertically averaged magnetic advection and magnetic diffusion compensate each other, the resulting magnetic flux transport velocity, ${\overline{\overline{u}}}_{\psi }=0$ (see Eq. (29)), and the gray lines in Fig. 1 become vertical.

Fig. 1. Time evolution of the magnetic field strength over time (normalised in units of the diffusion timescale of the stationary model at the outer DZ boundary, tdiff, stat(rDZ, out)), where the colour-coding represents the magnetic field strength, |B|, at every point in space and time. The gray lines show magnetic flux transport over time, using the effective, vertically averaged magnetic field transport velocity, ${\overline{\overline{u}}}_{\psi }$. The dashed white lines denote the inner and outer boundary of the DZ, respectively, whereas the black full line denotes the diffusion timescale, tdiff, stat, of our stationary model. The inner DZ boundary is located very close to the inner disk rim.

3.2. Remarks on the validity of stationary models

In this work most of our conclusions are drawn from the analysis of stationary PPD models. Since disks are evolving over time, it is important to discuss the validity of using stationary models. In our 1+1D models we use a vertically averaged viscous-α (see Eq. (5), based on the layered viscosity model of Gammie 1996, see also Steiner et al. 2021), which implies that viscous accretion varies in vertical direction between different layers. A localised perturbation changing α(r, z) at some radius r and some height z would lead to an increased or decreased local density, ρ(r, z), and therefore would cause a deviation from a hydrostatic vertical equilibrium at r. Additionally, the vertically averaged α would change at r, which over time would lead to a bump or gap in Σ(r). The development of such gaps or bumps would occur on the local viscous timescale, t_visc, whereas the return to a local vertical hydrostatic equilibrium would happen on the local dynamical timescale, t_dyn (e.g. Delage et al. 2022). Therefore, an approximate stationary model can be understood as a time-averaged disk, where such local perturbations are smoothed out under the assumption that the relaxation to a hydrostatic equilibrium happens much faster than the development of a bump or gap in Σ, or t_dyn(r)≪t_visc(r).

Additionally, we have to discuss the validity of a stationary large-scale magnetic field, which can be argued similarly in comparing the timescale of global magnetic field evolution to the timescale at which such perturbations typically occur. It has been argued by Dudorov & Khaibrakhmanov (2014) that the magnetic field in vertical direction also changes on t_dyn. In contrast, the large-scale magnetic field in radial direction evolves on the local diffusion timescale, t_diff, which requires t_dyn ≪ t_diff for a valid temporal averaging of local large-scale magnetic field perturbations.

Finally, PPDs typically disperse after ∼1 − 4 Myrs (e.g. Delage et al. 2022; Ben et al. 2025), which raises the question if such disks live long enough to allow the disk and large-scale magnetic field to evolve towards a stationary state (i.e. steady-state accretion with Ṁ = const.). This means that we also have to verify that the global evolution timescales of both disk and large-scale magnetic field, t_visc and t_diff, respectively, are shorter or of the order of typical disk lifetimes (e.g. Delage et al. 2022). In Fig. 2 we can see that the conditions, t_dyn ≪ t_visc and t_dyn ≪ t_diff, are both fulfilled everywhere in the disk, which allowed us temporally average over local perturbations of the vertical disk structure and vertical large-scale magnetic field profile. Furthermore, Fig. 2 shows that t_diff and t_visc are indeed shorter or of the order of a typical PPD lifetime, which justifies the assumption that the disk is likely to establish (near) steady-state accretion and therefore stationarity over most of its radial extent before it dissolves (e.g. Delage et al. 2022, for a similar argumentation).

Fig. 2. Timescales of our reference model (cf. Table 2) relevant to the long-term evolution of a PPD and a large-scale magnetic field. The solid blue, orange, and green lines correspond to the viscous timescale, tvisc; the dynamical timescale, tdyn; and the magnetic diffusion timescale, tdiff, respectively. The vertical, dashed black lines correspond to the corotation radius, Rcor, and the outer DZ boundary, rDZ, out. Finally, the gray area marks the region between 1 Myr and 4 Myrs and depicts the expected lifetime of a PPD.

3.3. The reference model

In this section, we describe the details of our reference disk model, especially the large-scale disk magnetic field topology. The field is strongly dependent on the resistivity profile, $\overline{\overline{\eta }}(r,z)$, which in turn is dependent on the ionisation fraction, x_e(r, z). Therefore, we investigated in more detail the vertical and radial profile of x_e and η, which then can be used to construct the conductivity-weighted vertical resistivity average, $\overline{\overline{\eta }}$, counteracting advection.

3.3.1. Large-scale magnetic field profile:

We present and analyse the stationary magnetic field structure. As already discussed in Sect. 3.1, the field in the inner disk is dominated by the stellar dipole field (see Eqs. (14)–(16)). Since B_z, ⋆ ∝ r⁻³, the vertical magnetic field, B_z, in the inner disk also follows this power law (see (a) in Fig. 3). The field evolution in the outer disk is strongly influenced by AD, which according to an analytical model from (Dudorov & Khaibrakhmanov 2014), should lead to a dependency B_z ∝ r^−5/8 in the case of an irradiated disk, which we can verify with our model (see (a) in Fig. 3). Amplification of the fossil field in the outer region is roughly B_z/B_z(t = 0)≈10, whereas in the inner disk, due to the dominant stellar field component, almost no amplification is observed. However, the disk field is advected inwards and a radial field component at the disk surface, B_r^s, develops, which leads to an inclination of the large-scale field with respect to the disk normal (see (b) in Fig. 3). The value of 30° shown in (b) of Fig. 3 corresponds to the critical field inclination, for which cold MCWs can be launched (e.g. Ogilvie 1997; Ogilvie & Livio 2001, and Sect. 4.2).

Fig. 3. Poloidal stationary magnetic field profile of our reference model (see Table 2). Panel (a) shows both the vertical and radial magnetic field component at the surface, Bz and Brs, respectively. Panel (b) depicts the ratio of Brs/Bz, which is a measure for the inclination of the field lines, whereas panel (c) shows the magnetic plasma beta β0 at the disk midplane. The vertical dashed line correspond to the corotation radius, Rcor, and rDZ, out, respectively. The dotted, gray lines in panels (a) and (c) describe the corresponding values at t = 0, whereas the horizontal dotted line in panel (b) denotes a surface angle of 30° with respect to the midplane normal.

Inside the DZ close to r_DZ, out, the radial field, B_r^s, decreases significantly compared to its value outside the DZ. This can be explained by comparing the viscous timescale, t_visc, to t_diff, stat: If t_diff, stat is significantly shorter than t_visc, then the large-scale field diffuses outwards before any significant inward-directed advection takes place. As a consequence, field lines are only weakly inclined and B_r^s stays small compared to B_z (see zone around t_diff, DZ in Fig. 4). Only very close to the inner rim t_diff, stat/t_visc ≈ 1, which means that there B_r^s ≈ B_z (cf. spike in B_r^s close to the inner rim in (b) of Fig. 3 and in Fig. 4). This steep increase in t_diff, stat close to r_in is a result of three effects, which occur in a very narrow radial region close to the inner rim: (i) At the transition from optically thin to optically thick gas a dip in T_gas (for details see below paragraph about the temperature profile and (b) of Fig. 6 leads to a reduced x_e,th and consequently to a larger $\overline{\overline{\eta }}$, which is equivalent to a shorter t_diff, stat. (ii) However, just slightly further inwards, the stellar irradiation starts to dominate the disk temperature, which leads to rising T_gas (see (b) of Fig. 6) and a steep decrease of $\overline{\overline{\eta }}$ once more (or equivalently to a sharp increase of t_diff, stat, see Fig. 4) (iii) Closest to r_in, $\overline{\overline{\eta }}$ rises again (equivalently t_diff, stat becomes shorter) slightly towards r_in which is due to small Σ and relatively large |B| provided by the stellar dipole field (we recall that η_A ∝ |B|² ρ⁻², also see Eq. (59)).

Fig. 4. Comparison of the viscous timescale, tvisc, and the diffusion timescale, (tdiff, stat). The vertical dashed lines close to the inner rim and at ∼4 AU represent Rcor and rDZ, out, respectively.

The magnetic plasma beta at the midplane, β₀, describes the relative importance of the large-scale field in the disk. We start the time evolution with a fossil field, which corresponds to β₀(r = r_out) = 10² at the outer disk rim, r_out. In the AD-dominated outer disk, the model of Dudorov & Khaibrakhmanov (2014) predicts a radial dependency β₀ ∝ r^−3/2 in the case of an irradiated, optically thick disk, which our results also reflect (see (c) in Fig. 3).

3.3.2. Resistivity profile:

The ionisation fraction profile in radial and vertical direction, x_e(r, z), is necessary to calculate the corresponding OD and AD profiles, η_O and η_A, respectively. Both thermal and non-thermal ionisation are taken into account (see Sect. 2.5) and show that the innermost region (r ≤ 0.044 AU) is optically thin due to low Σ and the gas temperature, T_gas, is dominated by stellar irradiation. Hence, the very inner disk has a high x_e; however, at ∼0.044 AU, the disk becomes sufficiently dense to turn optically thick for stellar radiation and the primary heating process transitions from stellar irradiation to viscous heating. Additionally, the disk has a DZ between 0.06 AU ≲ r ≲ 5 AU, where the inner part of the DZ close to the midplane r ≲ 0.2 AU is subject to thermal ionisation due to viscous heating (cf. region of large ionisation in the inner disk region close to the midplane in (a) of Fig. 5). Above z ≈ H_p, thermal ionisation is not effective anymore, since the disk is hydrostatically stratified and ρ above H_p becomes very small. Finally, the upper disk is penetrated mainly by X-rays and cosmic rays, yielding an expected significant x_e,nth in the surface layer.

Fig. 5. For better understanding of the regimes of thermal and non-thermal ionisation, the radial and vertical profiles of (a) ionisation fraction, xe; (b) the contribution of OD, ηO; (c) the contribution of AD, ηA; and (d) the combined resistivity, η(r, z) = ηO(r, z)+ηA(r, z), are plotted for our reference model, which is explained in Sect. 3.1 (cf. Table 2). The dashed gray lines in (a)–(d) correspond to specific heights above the disk midplane (z = 0, z = Hp, $z=H_{\mathrm{p}}+\frac{1}{3}(z_{\mathrm{B}}-H_{\mathrm{p}})$, $z=H_{\mathrm{p}}+\frac{2}{3}(z_{\mathrm{B}}-H_{\mathrm{p}})$, and z = zB). In (e) the combined resistivity is plotted for those specific heights, as well as the resulting vertically averaged resistivity $\overline{\overline{\eta }}$ (thick blue line). For comparison the turbulent resistivity, ηt, corresponding to a Prandtl number of 1 is also plotted, showing that ${\eta }_{\mathrm{t}}\ll \overline{\overline{\eta }}$ almost everywhere except for the very inner region.

Ohmic diffusion, η_O, is only significant inside the DZ, which is also in agreement with predictions made by other authors (e.g. Bai & Stone 2013; Delage et al. 2022). Interestingly, the DZ also spans above (and below) the thermally ionised region in the inner disk, which means that the very weakly ionised gas is shifted to intermediate heights, where thermal ionisation cannot be maintained due to low gas density and non-thermal ionisation sources cannot penetrate deep enough to ionise the gas. The disk surface layer is sufficiently well ionised by non-thermal processes, such that η_O ∝ x_e⁻¹ is very inefficient in the disk surface (see (b) of Fig. 5).

Ambipolar diffusion, η_A, is effectively quenched in the thermally ionised disk region, which is due to both high gas density, ρ, and large x_e (we recall that η_A ∝ |B|² ρ⁻² x_e⁻¹, see Eq. (65)). Inside the DZ, η_A is also negligible compared to η_O due to its strong dependency on ρ. In the surface layer as well as in the outer disk; however, AD becomes the dominant diffusion process. The magnetic field strength, |B|, is roughly independent of height, z, and therefore the relative importance of |B| compared to ρ increases with z (see (c) of Fig. 5). The combined resistivity profile η = η_O + η_A is plotted in panel (d) of Fig. 5. The regions of low η are the thermally ionised disk region, a thin layer at intermediate height z inside the DZ and in the outer disk close to the midplane. Recalling Eq. (26), the vertical averaging process of $\overline{\overline{\eta }}$ favours layers of low η, as in those regions magnetic flux transport, ${\overline{u}}_{\psi }$, is more efficient. As explained in Sect. 2.5 it is important for our model to vertically average the resistivity profile, $\overline{\overline{\eta }}$, due to constraint of a 1+1D model. Therefore, we evaluated (η_O + η_A)(r, z) at five different heights, including at the midplane, the pressure scale height, and the magnetic surface: z = 0, z = H_p, and z_B, respectively. Additionally, for $\overline{\overline{\eta }}$ to also incorporate the effects of the DZ above the thermally ionised region and of the disk layer of comparatively low resistivity between the DZ and the surface layer, $\overline{\overline{\eta }}$ is evaluated at two additional heights, $z=H_{\mathrm{p}}+\frac{1}{3}(z_{\mathrm{B}}-H_{\mathrm{p}})$ and $z=H_{\mathrm{p}}+\frac{2}{3}(z_{\mathrm{B}}-H_{\mathrm{p}})$ to adequately describe that in different radial sections different vertical layers contribute dominantly to $\overline{\overline{\eta }}$ (see (e) of Fig. 5):

• In the thermally ionised region (r ≲ 0.3 AU), η is dominated by the vertical disk layers close to the disk midplane (the dashed orange line and dash-dotted green line are closest to the blue full line for r ≲ 0.3 AU, see panel (e) of Fig. 5).

• In the region inside the DZ, but further out than the thermally ionised region (0.3 AU ≲ r ≤ r_DZ, out), η is dominated by a region above the vertical extent of the DZ but deep enough below the disk surface for non-thermal radiation not ionising the gas. In this intermediate region, AD is the dominant non-ideal diffusion process (the purple dash-dotted line is very close to the full blue line in this region, see panel (e) of Fig. 5).

• Outside of the DZ (r ≥ r_DZ, out), the disk is dominated by AD, which is least effective close to the disk midplane and therefore most dominant for calculating the average $\overline{\overline{\eta }}$ (the dashed orange line and dash-dotted green line are very close to the blue full line, see panel (e) of Fig. 5).

For comparison, the turbulent resistivity, η_t, resulting from a magnetic Prandtl number, 𝒫 = 1, is also included in the vertical average of $\overline{\overline{\eta }}$ and is taken to be vertically constant, which means that η_t is added at every height to η(r, z). It has been pointed out by other authors (e.g. Dudorov & Khaibrakhmanov 2014) that η_t is only important in the innermost hot disk with high x_e (cf. the magenta line in panel (e) of Fig. 5).

3.3.3. Temperature profile:

The disk temperature profile, T_gas(r), is important for determining the thermal gas pressure, P_gas, as well as the pressure scale height, H_p, and the magnetic surface, z_B. All of those quantities are important to get the magnetic flux transport and the large-scale magnetic field topology, (B_r^s(r), 0,  B_z(r))^T. The optically thick regions of the disk are dominated by viscous heating, whereas T_gas in the optically thin regions is determined by stellar irradiation and the ambient temperature (for implementation details, see Ragossnig et al. 2020; Steiner et al. 2021).

We highlight one important temperature feature close to the inner edge of the disk, resulting in a distinct dip in the temperature of the disk (see (b) of Fig. 6). At the transition between the optically thick (white) and the optically thin (grey) region, we compare the individual energy contributions. The pressure work rate, P_gas∇⋅u, has a negative sign, which means that pressure aids the gas to be pushed inwards due to −ΔP_gas < 0 in this area. Additionally, the net energy advected by the radial gas flow into a disk annulus just outside R_cor is positive, because of the stellar torques acting on the inner disk. However, for r < R_cor the stellar dipole brakes the gas and gas starts to fall radially faster towards the star. This leads to Ė_adv < 0, because more material exits than enters any disk annulus at r < R_cor (see (a) of Fig. 6). The combined effect of pressure work and advection at r ≈ R_cor, P_gas ∇ · u + Ė_adv < 0, leads to a net drain of inner energy e from the disk annuli in this area and is equivalent to a lower T_gas, 0. For our reference model, this leads to a drop below T_gas, 0 < 1500 K (see (b) of Fig. 6). At this temperature, some dust can still survive (e.g. Ferguson et al. 2005).

Fig. 6. Panel (a), energy contributions to the overall inner energy budget for every disk annulus at radius r. Q : ∇u, Γ − Λ, Pgas ∇⋅u, 4 π Σ κR(J − S), and Ėadv, denote energy rates per unit volume due to viscous heating, heating or cooling over the surface, work done by gas pressure per unit time, radial radiation flux, and net inner energy advection, respectively (see Eq. (3)). Panel (b), plot of the midplane gas temperature, Tgas, 0. The vertical dashed line in all panels denotes the corotation radius, Rcor, and the gray area marks the optically thin region in the inner disk.