Magnetic field transport in geometrically thick discs: multi-dimensional effects on the field strength and inclination angle

We theoretically investigate the magnetic flux transport in geometrically thick accretion discs which may form around black holes. We utilize a two-dimensional (2D) kinematic mean-field model for poloidal field transport which is governed by both inward advection and outward diffusion of the field. Assuming a steady state, we analytically show that the multi-dimensional effects prevent the field accumulation toward the centre and reduce the field inclination angle. We also numerically investigate the radial profile of the field strength and the inclination angle for two geometrically thick discs for which (quasi-)analytic solutions exist: radiatively inefficient accretion flows (RIAFs) and super-Eddington accretion flows. We develop a 2D kinematic mean-field code and perform simulations of flux transport to study the multi-dimensional effects. The numerical simulations are consistent with our analytical prediction. We also discuss a condition for the external field strength that RIAF can be a magnetically arrested disc. This study could be important for understanding the origin of a large-scale magnetic field that drives jets and disc winds around black holes.


INTRODUCTION
To understand the origin of magnetically driven outflows, it is crucial to elucidate the distribution of large-scale poloidal magnetic fields, including their strength and inclination (e.g., Contopoulos & Lovelace 1994;Kudoh & Shibata 1997;Fukumura et al. 2010;Jacquemin-Ide et al. 2019;Dihingia et al. 2021).This is because the large-scale magnetic fields neither dissipate through local magnetic diffusion nor vanish when accreted on to black holes.The disc acquires magnetic fields from the accretion flow from outside the disc during its growth (e.g., Lovelace 1976;Hawley et al. 2015;Takasao et al. 2022), and the imported magnetic fields distribute within the disc according to the disc's advection and effective magnetic diffusion (e.g., Lubow et al. 1994;Beckwith et al. 2009).In addition to the magnetic flux transport, disc dynamo could also be important for the generation of the disc poloidal field (e.g., Stepanovs et al. 2014).However, the effectiveness of dynamos seems to depend on the initial magnetic field strength and disc thickness and the importance remains elusive (e.g., Hogg & Reynolds 2018;Liska et al. 2020).Recent AGN observations suggest that Radio Loud AGNs with jets have stronger magnetic fields compared to Radio Quiet AGNs without jets (Lopez-Rodriguez et al. 2023) at a scale of ≳ 5 pc, emphasizing the importance of magnetic flux transport.
The amount of magnetic flux brought into the central black hole is crucial for the driving of relativistic jets (e.g., Blandford & Znajek 1977;Tchekhovskoy et al. 2011).It is believed that when a rotating black hole accumulates a significant amount of poloidal magnetic fields, it forms a magnetically arrested disc (MAD) that can drive powerful jets (e.g., Narayan et al. 2003;Tchekhovskoy et al. 2011; ★ E-mail: ryoya@astro-osaka.jpMcKinney et al. 2012).Supermassive black holes like Sgr A* exhibit recurrent flares (GRAVITY Collaboration et al. 2018), and recent simulation studies suggest that flares can be driven by magnetic reconnection in the inner disc if it is in the MAD state (Dexter et al. 2020;Porth et al. 2021;Ripperda et al. 2022) (a similar flare mechanism is also discussed in the context of accretion of protostars (Takasao et al. 2019)).Recent observations of X-ray binary (XRB) systems suggest the existence of MAD (Zdziarski et al. 2022), and future X-ray spectroscopic observations will find evidence of MAD (e.g., Inoue 2023).Understanding the conditions for MAD manifestation is crucial for revealing the origin of jets.
However, the mechanism determining the disc magnetic field distribution remains poorly understood, making it challenging to discuss the magnetization of black holes.The transport of global poloidal magnetic fields is one of the challenging problems in astrophysics, primarily due to the need to simultaneously track a wide range of scales.The typical size of the X-ray binary disc is considered to be ∼ 10 5 − 10 6  g ( g =  / 2 is the gravitational radius) (e.g., Alfonso-Garzón et al. 2018;Hynes et al. 2019).Ressler et al. (2020a,b) performed three-dimensional magnetohydrodynamic (3D MHD) simulations covering up to seven orders of magnitude in scale using a zoom-in technique.Such an approach is powerful when the backreaction from smaller scales to larger scales is negligible.As the smaller scales can affect the larger scales during the system's evolution, there is still a need to solve a wide range of spacetime simultaneously.
3D MHD simulations that cover a wide range of scales are numerically challenging.For this reason, a one-dimensional kinematic mean-field model (1D model), which integrates the disc structure in the vertical direction, has been widely used (Lubow et al. 1994;Okuzumi et al. 2014;Guilet & Ogilvie 2014).Such a model is use-ful for modelling geometrically thin discs such as a standard disc (Shakura & Sunyaev 1973).Lubow et al. (1994) demonstrated for geometrically thin discs that the following dimensionless parameter  is key for determining the efficiency of inward magnetic flux transport and the inclination of the magnetic field: where  m represents the effective magnetic Prandtl number due to turbulence, and ℎ is the aspect ratio of the disc.
It is believed that jet-related discs are thick discs in many cases (e.g., Cao 2011;Begelman & Armitage 2014;Dhang et al. 2023).Although 1D models are a powerful tool for geometrically thin discs, a simple application to geometrically thick discs could be problematic because the model ignores a detailed vertical structure.There have been theoretical studies examining the multi-dimensional effects of discs.Guilet & Ogilvie (2012) studied flux transport by considering the vertical structure of the discs, but their models adopt the thin disc approximation (see also Lovelace et al. 2009).In addition, the radial field distribution (or the radial gradient of the field strength), which is intrinsically a result of magnetic flux transport, is treated as an input parameter.Li & Cao (2021) investigated flux transport in vertically structured discs using their 2D axisymmetric model.However, the impact of the geometrical thickness on the disc field structure remains unclear.
In this study, we analytically and numerically investigate the distribution of poloidal magnetic fields achieved in geometrically thick discs around black holes.We examine the magnetic field distribution in two types of geometrically thick discs: radiatively inefficient accretion flows (RIAF; e.g., Narayan & Yi 1994;Yuan & Narayan 2014) and slim discs (e.g., Abramowicz et al. 1988).The structure of this paper is as follows.Section 2 provides an overview of past 1D models regarding magnetic flux transport.We consider the vertical structure of the disc's internal magnetic field, which was partly ignored in the 1D approximation.We analytically derive the parameter  eff applicable to thick discs.We then introduce 1D and 2D kinematic mean field models used for the numerical validation of  eff .We also present the disc models used for validation.Section 3 presents detailed results of the magnetic field distribution, focusing on the steady state.In Section 4, our results are compared with previous studies with a particular focus on the vertical magnetic field distribution.We also discuss the onset condition of MAD in a RIAF.Section 5 provides a summary of our results.

Overview of kinematic 1D mean-field model
Here we introduce a 1D model formulated by Lubow et al. (1994).The evolution of the axisymmetric poloidal magnetic field in the disc (− <  < , where  ≈  S /Ω,  S is sound velocity, Ω is rotational angular velocity) is calculated in cylindrical coordinates (, , ).The aspect ratio of the disc is denoted as The disc is assumed to be geometrically thin (ℎ ≪ 1), and the magnetic field within the disc is approximated as follows: where  ,surf =   (,  = ) and  ,mid =   (,  = 0).The induction equation governs the time evolution of the poloidal magnetic field.In an axisymmetric configuration, it is helpful to utilize a flux function (, ) defined by B =   e  +   e  = ∇ × (e  /).The following relations hold: The induction equation for (, ) is given by equation ( 10) of Lubow et al. (1994): where   is the radial advection velocity and  is the magnetic diffusion coefficient, which depend on a disc model (see Section 2.2).The azimuthal current density   is given by The 1D model averages equation ( 7) in the vertical direction with weighting by 1/ over the range of − <  <  (e.g., Ogilvie & Livio 2001;Okuzumi et al. 2014), which yields where 1 Here,  0 = (,  = 0) is the flux function on the equatorial plane of the disc, and  is assumed to be approximately constant in the vertical direction.  is the azimuthal component of the surface current defined as We determine   from Biot-Savart's law using  0 (see e.g., Lubow et al. 1994).Using equation ( 8), we also get This is the equation to calculate  ,surf in the 1D models.
The second term in the right-hand side of the equation ( 13) includes multi-dimensional effects as it depends on the vertical structure and the disc thickness (e.g., Okuzumi et al. 2014).Using equation (4), we find where We will examine the validity of equation ( 14) for geometrically thick discs by performing two-dimensional numerical calculations (see Section 3).When | ,surf / ,mid | ≫ /, we can further simplify the equation ( 13) as follows: (e.g., Lubow et al. 1994;Guilet & Ogilvie 2014).Under the assumption, the radial component of the magnetic field on the disc surface can be expressed as We note that  ,mid is unaffected by the choice of the calculation methods of  ,surf as long as the accretion velocity and the resistivity do not depend on  ,surf (see the induction equation 9).  is calculated exactly from Biot-Savart's law.
In the steady state, the following equation holds true based on equations (9): Substituting equation ( 18) into equation ( 19) yields the following equation: is a dimensionless parameter defined as Equation ( 20) indicates that  −1 is a measure of the inclination of the magnetic field at the disc surface.The situation of  < 1 corresponds to the case where advection is dominant.In such a case, the magnetic flux is efficiently transported to the centre, and the inclination of the magnetic field at the disc surface is also large.In contrast,  > 1 corresponds to the case where diffusion is dominant.The magnetic flux transport toward the centre is inefficient, and the magnetic field is more vertical at the disc surface.Thus,  serves as an indicator of the efficiency of magnetic flux transport in the 1D model (see Section 4 of Lubow et al. 1994).

Evaluation of the multi-dimensional effect
The 1D model ignores a detailed vertical magnetic field structure inside the disc by assuming that the disc is thin.The assumption breaks down for geometrically thick discs.As a result, the vertical distribution of the magnetic field can deviate from a linear distribution, as represented by equations (3) and (4).Here, we analytically evaluate the effect of the thickness of the disc on the vertical distribution of the magnetic field inside the disc.We expand the magnetic field in the disc as follows: where Substituting equations ( 22) and ( 23) into the induction equation for a steady state, we obtain the following: After some calculations, we get the following equation: where When the vertical dependence on  and   is sufficiently small such that  = O ( 0 ) and   = O ( 0 ), equation ( 26) requires that the following two relations must hold true regardless of :

𝐵
(2) Equation ( 31) is an extension of equation ( 20).By comparing the two equations, we can rewrite equation (31) as where eff is similar to  but includes multi-dimensional effects due to the thickness of the disc and magnetic field gradient.One can easily find that  eff ≈  in the limit of a thin disc (ℎ ≪ 1).When the magnetic field strength decreases monotonically outward from the centre,    0 < 0, and thus  eff > .The difference between them becomes more significant as ℎ increases.For a thick disc (ℎ ∼ 1),  eff > , which indicates that the multi-dimensional effects effectively increase magnetic diffusivity.Jet-related discs are commonly thick discs, but the multi-dimensional effects suppress the magnetic flux accumulation.Considering these results, careful investigations of flux transport are necessary for thick discs.
The assumption of the linear distribution (equations 3 and 4) will break down when multi-dimensional effects are important.To show this, we evaluate the magnitude of higher-order components of the magnetic field.By substituting equations ( 22) and ( 23) into the divergence-free condition ∇ • B = 0, we obtain where Equation ( 35) requires that where we have used equation (31) to eliminate (1) and defined Equation ( 37) is a direct comparison between the 0th and 2nd order components of   .If (2) / (0)  ≈ 1, the linear distribution approximation given by equation ( 4) is invalid.Since (2) / (0)  ∝ ℎ 2 , higher-order components should be more significant for thicker discs.

Geometrically thick disc models
To highlight the impact of disc thickness on magnetic flux transport, we perform numerical calculations for geometrically thick discs.We compare the results of the 1D and axisymmetric 2D models to examine the limitation of the assumptions in the 1D model.
Our disc models are the following two: the self-similar solution of the radiatively inefficient accretion flow (RIAF; Narayan & Yi 1994) and the quasi-analytical solution of the super-Eddington accretion flow (SEAF, slim disc and standard disc; Watarai 2006).They are representative models of geometrically thick discs around black holes.Both are advection-dominated accretion flows (ADAF), but RIAF and SEAF correspond to optically thin and thick discs, respectively.The following scaling laws hold for both models: where  is the -parameter (Shakura & Sunyaev 1973).Σ and v are the surface density and density-weighted average advection velocity defined as respectively. is the advection parameter, defined as the ratio of the advective cooling rate per unit area to the viscous heating rate per unit area, i.e., We have  = 1 in cases of no radiative cooling.The RIAF solution is a self-similar solution obtained when  is a constant value throughout the disc.The aspect ratio ℎ is an increasing function of  .Assuming isothermal in the vertical direction, the density is written as where  = √︁ Π/Σ/Ω.This density structure is used only for visualization.A different semi-analytical solution has been proposed recently (Xu 2023).The application of such other disc models will be our future work.
Watarai ( 2006) constructed the SEAF solution by assuming that the disc radiation is blackbody radiation and the radiation pressure dominates over the gas pressure ( rad ≫  gas ).The radius dependence of  is expressed as where  is a dimensionless constant of order unity (see equation 24 of Watarai 2006)1 .Also,  ≡ (/  )( /  Edd ) −1 ( S = 2 / 2 and  Edd =  Edd / 2 are the Schwarzschild radius and the Eddington accretion rate, respectively). takes a value in the range 0 ≤  ≤ 1. SEAF has the characteristic radius that determines the slim region.Photons travel diffusively in the disc because of their high gas density.When the accretion time becomes shorter than the diffusion time of photons, the photons are trapped inside the disc.The photon trapping radius can be expressed as  trap ≈ /(/)  ( is the optical depth).Within the photon trapping radius, the cooling effect of radiation becomes inefficient, causing the disc to expand and form a slim disc (e.g., Abramowicz et al. 1988;Watarai 2006).In terms of the accretion rate, the photon trapping radius can be written as For  ≲  trap ,  → 1 (  ∝  0 ) and the advection cooling is dominant.Such a geometrically and optically thick disc is commonly referred to as a "slim disc".For  ≳  trap ,  ∝  −2 and the radiative cooling dominates the advection cooling.Therefore, the outer disc can be described as a standard disc.The quasi-analytical solution by Watarai ( 2006) has also been reproduced in recent 2D radiation hydrodynamic simulations (Kitaki et al. 2021;Yoshioka et al. 2022).
The density structure of SEAF is introduced only for better visualization of the results.The temperature of SEAF has a non-isothermal distribution given by  (, ) =  mid ()(1 −  2 / 2 ), and it is assumed that the polytropic relation  ∝  Γ holds.Therefore, the density distribution is given by where  = 1/(Γ − 1) (polytropic index) and  = (2 + 3) √︁ Π/Σ/Ω.SEAF has a photon trapping radius (approximately 85 S ) and a radius where the "cold" (gas-pressure-dominated) standard disc appears (approximately 500 S ), resulting in a significantly varying aspect ratio in the radial direction.
We describe the fiducial sets of parameters in this study.Both models assume a black hole mass of 10 ⊙ ,  = 0.01, and a specific heat ratio of Γ = 4/3.Regarding , we refer to the results of past 3D MHD simulations on MRI turbulent discs, which show  = O (0.01) − O (0.1) (e.g., Hawley et al. 2013;Suzuki & Inutsuka 2014;Takasao et al. 2018).We adopt the mass accretion rates of 10 −3  Edd and 100  Edd for RIAF and SEAF, respectively.The disc structure of RIAF remains unchanged regardless of the mass accretion rate, while the mass accretion rate affects the photon trapping radius in SEAF.For the RIAF model, we adopt a fiducial value of  = 1.The aspect ratios ℎ for these parameter values are ℎ ≈ 0.53 and ≈ 3 for the RIAF and slim disc in SEAF, respectively.We also have  ≈ 1.77 for the SEAF solution (equation 46).
We characterise the density-weighted magnetic diffusivity using the magnetic Prandtl number where the density-weighted viscosity coefficient ν is given by Here, Ω =   / and Π ≡ ∫ ∞ −∞  represent the disc rotation angular velocity and surface pressure, respectively, both of which are as we have already used  (equation 21), we instead use  to denote the constant to avoid confusion.
provided by the disc model.As the relation   ≈ ν/ holds in a viscous disc, we obtain equation (1).

Numerical methods
To directly investigate the multi-dimensional effects discussed in Section 2.1.2,we perform magnetic flux transport calculations using a two-dimensional axisymmetric model in spherical coordinates (, , ) (hereinafter referred to as the 2D model), and compare the results with those obtained using a 1D model.We describe the computational methods for the 1D and 2D models below.

1D model
We solve equation ( 9) and use the equation ( 18) to obtain  ,surf .
The surface current density   in the diffusion term (the second term of equation 9) is calculated in the same way as in Lubow et al. (1994) (see Section 3 of Lubow et al. 1994), but we set  = 10 −3 in our calculations, whereas Lubow et al. (1994) used their equation ( 27).We calculate the advection term (the first term of equation 9) using second-order upwind differencing with the Monotonized Central (MC) interpolation and the MUSCL method (van Leer 1997) (first-order upwind differencing was used in Lubow et al. 1994).We use second-order strong stability preserving Runge-Kutta methods (SSPRK, Gottlieb et al. 2009) for time integration.The time step is determined by the CFL condition: where Δ  =  +1/2 −  −1/2 (  denotes the position of the -th cell center), and the CFL number is  CFL = 0.4.The logarithmically spaced cells are used, and the cell number   is set to match the spatial resolution of the corresponding 2D models.
The discs are initially threaded by a uniform magnetic field, (,  = 0) = (1/2) 2  0 ( 0 is the uniform imposed field strength).The field evolves in response to the advection and effective magnetic diffusion. *  and  * are given according to the disc models and are assumed to be constants with time.We ignore the back reaction of the Lorentz force on the disc flows.We also assume that   and  are uniform in the vertical direction of the disc and set  *  = v and  * = η =  m ν for simplicity.The velocity generally varies with height in response to the coupling between the magnetic field and the plasma, as suggested by 3D MHD simulations (e.g., Bai & Stone 2013;Lesur et al. 2013;Suzuki & Inutsuka 2014;Takasao et al. 2018).However, there are no established analytic methods that selfconsistently describe both the radial and vertical velocity structures (for the analytic study of the local disc, see, e.g., Guilet & Ogilvie 2012).For this reason, we adopt the vertically uniform velocity in this study.
As the inner boundary condition, we impose that   = const.This condition approximates the boundary condition of the 2D models.We note that a simple outflow boundary ( = const or   = 0) behaves as a sink of the field (see also Appendix A).
For comparisons with 2D models, we define the radius of an effective disc outer edge inside the numerical domain,  trunc .We set  *  =  * = 0 for  >  trunc .

2D model
The time evolution equation for the poloidal magnetic field is solved using the flux function  in the 2D spherical coordinate system (, ).
The relationship between  and the components of the magnetic field is given as follows: Here,  satisfies B = ∇ × (e  / sin ).Therefore, the induction equation is given by The first and second terms of equation ( 56) are the advection terms, and the third and fourth terms are the diffusion terms.Note that our code can handle the time evolution of the magnetic field, unlike the 2D model of Li & Cao (2021).
The velocity field (  ,   ) is given by a disc model within the disc and is set to 0 outside the disc.Assuming that the disc surface is located at a height of scale height  from the equatorial plane, we set where (, ) = ( cos ,  sin ) and Δ = 0.1.We then set   and  as The magnetic diffusion coefficient  is assumed to be constant in the vertical direction, including above the disc (i.e., (, ) = η() =  m ν()), and its value at each radius  is determined by the disc model.Although this study focuses only on steady-state solutions, the steady-state is independent of the method used to determine the magnetic diffusion coefficient outside the disc.To prevent numerical instability near the poles, we set (, ) = 0 for the physical domain mesh that borders the  direction boundaries.We also set   and  to zero for  >  trunc , respectively, to prevent magnetic field inflow/outflow from the external boundary for  <  trunc .Setting the velocity field to zero and the magnetic diffusion coefficient to a finite value above the disc ensures that the magnetic field distribution above the disc in the steady state is a potential field (i.e.,   = 0).This is consistent with the requirements of a 1D model.The magnetic field distribution within and above the disc in the steady state does not depend on the magnitude of the magnetic diffusion coefficient above the disc.
We solve the advection terms with the second-order upwind scheme using the MC limiter and the diffusion terms with secondorder central differencing, achieving second-order spatial accuracy.We also integrate the time evolution with second-order SSPRK same as the 1D model.The time step is determined by the CFL condition ( CFL = 0.4), given by where the subscripts  and  denote the labels of the radial and latitudinal locations of the cell centre.The domain is descretized such that Δ  =  +1/2 −  −1/2 and Δ  =  +1/2 −   −1/2 .We set a uniform vertical magnetic field for the initial condition, similar to the 1D model (i.e.(, ,  = 0) = (1/2) 0 ( sin ) 2 ).
We set the boundary conditions for  at the inner and outer boundaries in the  direction such that    = const (see Appendix A for details).Additionally, we set reflective boundary conditions for  in the  direction.
We adopt a logarithmic mesh structure in the  direction with the same number of meshes as the  direction of the 1D model.For the  direction, we use a compressed non-uniform mesh so that the grid spacing at the equator is 2.5 times smaller than that at the poles.Table 1 summarises the information about the domain size and resolution.As shown in Appendix B, our numerical resolution is sufficient to obtain the converged results.Our code for the 2D models utilizes the framework of the publicly available MHD code, Athena++ (Stone et al. 2020).

Overview
We first review the results for  m = 1.The dependence on  m is investigated in Section 3.2.2.

Overview of steady-state solutions
In this study, we focus only on the steady state.For a better comparison between 1D and 2D models, the results of the 2D model are displayed in cylindrical coordinates.The solid lines in Fig. 1 represent the poloidal magnetic field structures of the RIAF and SEAF obtained from the 2D model.The disc density distribution is also indicated by colour for better visualization of the disc structure.The disc aspect ratios ℎ are approximately 0.5 for RIAF and 3 for the slim disc region in SEAF.The upper panel represents the global magnetic field distribution inside  <  trunc , while the lower panel represents the magnetic field distribution focused on the inner region (approximately the slim region in SEAF).Furthermore, in this section, we assume a fiducial magnetic Prandtl number of  m = 1.
The timescale until the magnetic flux transport reaches a steady state can be estimated by the relaxation timescale defined as  relax () ≡ / * (e.g., Lubow et al. 1994;Lovelace et al. 2009).The relaxation timescale at  =  trunc is approximately  relax ≈ 1.5 × 10 5  g / for RIAF and  relax ≈ 3.5 × 10 8  g / for SEAF.The state shown in we find that the distribution of  ,mid is more diffusive in the 2D model.
The 1D model exhibits a steeper  ,mid distribution near the inner boundary than the 2D model.This steep structure is seen in previous studies (e.g., Lubow et al. 1994;Okuzumi et al. 2014;Guilet & Ogilvie 2014).The 1D model does not allow the magnetic field to go inside the inner boundary, which results in the accumulation of the field there.On the other hand, in the 2D model, the magnetic field can spread around the polar regions (Fig. 1), which leads to a smoother distribution.Fig. 3 illustrates the differences in the magnetic field distribution between the 1D and 2D models near the inner boundary.Fig. 2 indicates that the 1D model can significantly overestimate the magnetic field strength near the central object.This result is important when discussing the injection of a magnetic field into the central object.

Importance of multi-dimensional effects
We analyse the 2D results to study the multi-dimensional effects.Fig. 4 show the radial distributions of  and  eff and compare them Table 1. min and  max are the minimum and maximum spherical radii of the calculation domain of 2D models, respectively.  is the number of cells in the  direction of 2D models, while   is that in the  direction of 1D models.We set   =   .  is the number of cells in the  direction of 2D models.(/Δ  ) min and (/Δ ) max are the minimum and maximum resolutions within the disc, respectively.The aspect ratios of cells, Δ/Δ, at the pole and equatorial plane are also listed. are obtained by fitting the magnetic field distribution with a 15th-order function in terms of  (= /).In the case of RIAF, we find   /  | = ≈  −1 eff except near the inner boundary (see black dotted and red lines).In the case of SEAF, we find a greater mismatch between   /  | = and  −1 eff .As we will see later, the mismatch is caused by higher-order components in the magnetic fields.Nevertheless,  −1 eff provides a much better approximation for   /  | = (red line) than the classical estimate  −1 (blue line) over a wide range of radius.We note that  −1 significantly overestimates the inclination in both disc models, which emphasizes the impact of the disc thickness.
The result of  eff >  indicates that the effective magnetic diffusivity is enhanced by multi-dimensional effects.The multidimensional effects are prominent only when  m ≲ 1.We will see this in Section 3.2.2.
The field inclination is a key factor for driving the magnetocentrifugal wind (Blandford & Payne 1982, hereafter BP wind).If the field inclination is larger than 30 • , the disc may drive the BP wind.The 1D model predicts that the BP wind can appear in both RIAF and SEAF.However, in the 2D model, the condition for the BP wind is not satisfied for  m = 1.
In Section 2.1.2,we show that the multi-dimensional effects of the disc thickness lead to nonlinear field profiles in the disc vertical direction. in equations ( 22) and ( 23), respectively).In the case of RIAF, the 1D approximation almost holds around the equatorial plane, but   deviates from the dashed line near the disc surfaces.The breakdown of the 1D approximation is more pronounced in SEAF.The higher-order components of the magnetic fields are significant even around the equatorial plane.SEAF shows larger deviations than RIAF because it has a larger aspect ratio ℎ (RIAF has ℎ ≈ 0.5, while the slim region of SEAF has ℎ ≈ 3).
The nonlinearity is examined by looking at  and the prediction of equation ( 37).One can see that equation ( 37) holds well at all radii except in the very vicinity of the inner boundary.This result demonstrates the validity of the discussion in Section 2.1.2,particularly the prediction of equation ( 37) that (2) As mentioned in Section 2.1.1,the 1D model can partially take into account the multi-dimensional effects using equation ( 14).We examine the applicability of this 1D approach through the comparison with the 2D model.Fig. 9 compares the different predictions of the surface field,  ,surf .The blue dashed lines denote the results of the 1D model using the relation  ,surf =   /2 (equation 18).The blue dotted lines indicate the results of the 1D model using the relation  ,surf =   /2 + ℎ   0  ,mid (equation 14), which takes into account the field gradient.The red solid lines show the results of the 2D models.In the case of RIAF, the 1D prediction of equation ( 14) agrees well with the 2D result except near the boundaries.In the case of SEAF, the 1D prediction matches the 2D result only outside the slim region.We note that  ,surf based on equation ( 14) is negative in the inner region.Fig. 8 compares the poloidal field structures of SEAF obtained by the 1D and 2D models.It is shown that the field near the centre is inclined inward, which is unrealistic.In Appendix C, we demonstrate for the SEAF model that the unphysical structure emerges because the 1D model cannot handle the multidimensional  effects correctly even if the radial gradient of  ,mid is considered in the calculation of  ,surf .
Fig. 8 shows that in the slim region (i.e.inside the photon trapping radius), the field structure of the 1D model significantly deviates from that of the 2D model, although they match in the standard disk region (outside  trap ).This is because higher-order components in the magnetic fields are critically important to determine  ,surf in the slim region.For example, the fitting to our 2D result shows 72, 4.30, −7.93, 8.45 for  = 3, 5, 7, 9, respectively, at  = 20  .
The following analysis demonstrates an example of the interplay among high-order terms.Using equations ( 31), (32), and (37), we can express (3)  as follows:  At  = 20  (inside the slim region), ℎ ≈ 2.5,   ≈ 1, and  −1 * ≈ 1.5. −1 * +  0 ≈ 0.3.In addition, the gradient of the highorder term    2 is approximately −1.4 and plays a role in increasing the absolute magnitude of   | ≈ 1.7, which shows the significance of the third-order component and explains the mismatch of the field structures of the two models (Fig. 8).In the standard disc region, however, the high-order term is unimportant mainly owing to its thickness and weak radial dependence of This explains the good agreement between the field structures of the two models in that region.

Dependence on aspect ratio
As we are interested in the flux transport in thick discs, we study the dependence of the solution on the disc aspect ratio ℎ, with a particular focus on RIAF.In the RIAF solutions, when we  , which denotes the relative importance of the second-order term.The orange solid lines show the numerical results, and the blue dashed lines denote the analytic prediction (equation 37).The results are for  m = 1.set  = 1, 0.01, 0.001, the disc aspect ratios are approximately ℎ ≈ 0.5, 0.1, 0.03, respectively.We use different mesh numbers for discs with different thicknesses to resolve the disc structure.Namely, (  ,   ) = (128, 80), (360, 200), (800, 460) for the cases  = 1, 0.01, 0.001, respectively.Again, we align the number of mesh points in the  direction of the 1D model with the number of mesh points in the  direction of the 2D model.Fig. 10 presents the results with different ℎ.The upper panel displays  ,mid .For the thickest model (ℎ ≈ 0.5), the results differ significantly between the two models.However, the 2D model shows results close to the 1D model for ℎ ≲ 0.1 except near the inner boundary, which suggests that the assumption of the linear field distribution is valid for such thin discs.The result that a thinner disc has a weaker magnetic field is also expected because a thinner disc shows a larger  for a given  m .
The lower panel of Fig. 10 shows the inclination of the field at the disc surface.Again, the discrepancy between the 1D and 2D models is significant for the thick disc but is insignificant for thinner discs (ℎ ≲ 0.1).Therefore, the higher-order components of the magnetic field are unimportant for thin discs, and  −1 eff ≈  −1 is a good measure of the field inclination (equation 20).

Dependence on effective magnetic Prandtl number
3D MHD shearing box simulations suggest that the magnetic Prandtl number  m is of order unity (Guan & Gammie 2009;Lesur & Longaretti 2009;Fromang & Stone 2009;Käpylä et al. 2020).In geometrically thick discs, the magnetic Prandtl number and  are related as  ∼ (ℎ m ) −1 ∼  −1 m (see equation 1), which implies that the efficiency of magnetic flux transport strongly depends on the local quantity  m .Considering this, we investigate the  m dependence of the solutions by performing a parameter search.
Figs. 11 and 12 compare the results with different  m .Fig. 11 shows the difference in the field structure between the models with  m = 2 (left) and 1 (right).The figure shows that the discs with larger  m accumulate magnetic fields more efficiently.It is found that even a twofold increase in the magnetic Prandtl number leads to significant changes in the magnetic field.Fig. 12 quantitatively compares the distributions of  ,mid for different  m .When comparing the cases of  m = 0.3 and  m = 3.0, there is an order of magnitude difference in the magnetic field strength near the inner boundary for RIAF and about four orders of magnitude difference for SEAF.Therefore, the distribution of the disc field strongly depends on the local properties of effective viscosity and magnetic diffusion arising from turbulence and other factors.
Since RIAF is a self-similar disc model,  ,mid follows a powerlaw dependence on .We investigate the dependence of the powerlaw index on  m .Fig. 13 illustrates the relationship between  m and the power-law index, ln ,mid /ln.From the figure, we observe that the behaviour of the power-law index is similar between the 1D and 2D models.Furthermore, the data points for the 1D model match the analytical solution derived by Okuzumi et al. (2014) (dashed line).The 2D model exhibits a smaller power-law index than the 1D model due to the multi-dimensional effect.When  m = 1, the power-law indexes in the 2D model and the analytical solution by Okuzumi et al. (2014) are    0 ≈ 0.55 and    0 ≈ 0.69, respectively.Therefore, they agree with approximately 80% accuracy.For  m ≲ 1, the powerlaw index is proportional to  m , while in the advection-dominated regime of  m ≫ 1, the power-law index approaches −2.In the limit of advection dominance, the magnetic field distribution asymptotically follows   ∝   −2 , as analytically derived by Okuzumi et al. (2014).Since our RIAF model possesses  = const., the asymptotic behaviour is consistent with the analytical estimation.By fitting the data in Fig. 13, we derive an approximate expression for the power-law index as a function of  m : Figure 10.Dependence of the field structure on the disc aspect ratio ℎ for the RIAF case.The coloured solid lines show the results for 2D models, and the dashed lines display the 1D results.The grey line in the upper panel denotes the initial profile, while the grey line in the lower panel denotes the field inclination angle corresponding to 30 • .We set  m = 1.The inclination for the 1D results is calculated using equation ( 18).tanh fixed to −2 as  m → ∞ (solid lines in Fig. 13): The maximum error of this fitting function is approximately 7% for the 2D model and 13% for the 1D model.By using the results for the 2D model (equation 64) and equation (34), we can obtain the relationship between  eff and  m as follows: The higher-order components of the magnetic field also depend on  m .Using equations ( 37) and ( 64), we have / (0) ∝  2 m and strongly depends on  m .However, since the magnitude of   64) and ( 65).The dashed black line represents the analytical solution for the power-law index in the 1D model (Okuzumi et al. 2014).10 6 , respectively.For comparison, giving the typical parameters of RIAF in this study (ℎ = 0.5,  = 0.01), the inclination of the magnetic field   /  | = for  mid = 10 3 and 10 6 becomes approximately 16 m and 47 m , respectively.However, our results indicate that   /  | = ≈ 0.36 m when  m ≲ 1 (equation 66).Therefore, equation ( 68) predicts larger values by approximately an order of magnitude than ours.
We consider that the differences between our study and Guilet & Ogilvie (2012) originate from the following two points.The first point is that they approximate the disc thickness as negligible.However, as we have shown in Section 3.1.2,the multi-dimensional effects due to the disc thickness enhance the effective magnetic diffusivity.The second point is the behaviour of the accretion flow.In their model, the vertical range of the accretion flow increases as the magnetic field weakens.Therefore, a weaker magnetic field results in an effectively thicker disc (ℎ is larger, and  is smaller).On the other hand, our model assumes that the accretion flow exists only within the vertical range of || < .Future detailed investigations using 3D MHD simulations of accretion flows will enable us to model a realistic velocity structure.Lovelace et al. (2009) argued that  m ≳ 2.7 is required for the disc to blow the BP wind.Although their constraint is similar to ours ( m ≳ 2 from Fig. 14), their argument is based on the thin disc approximation.In fact, the field inclination (their equation 13) is essentially the same as equation ( 20) which ignores the disc thickness effect.On the other hand, they solved not only the induction equation but also the momentum equation to take into account the outflow motion.Including the back-reaction of the outflow is our future work.
We have ignored the vertical structure of the accretion speed in this study just for simplicity.However, near the disc surfaces, efficient angular momentum loss due to the magnetic field can lead to a very fast accretion, as seen in MHD simulations (e.g., Matsumoto et al. 1996;Beckwith et al. 2009;Zhu & Stone 2018;Takasao et al. 2018Takasao et al. , 2019;;Mishra et al. 2020;Jacquemin-Ide et al. 2021).Such a coronal accretion will enhance the efficiency of the field transport toward the We briefly compare our 2D model with the model of Li & Cao (2021).Their model is based on steady-state equations, but our model can solve the time evolution of the system.Therefore, our model enables us to study how the disc magnetic field builds up to drive BP winds, for example.Our 2D model can accept more complicated distributions of the physical quantities and can also handle coronal accretion.To demonstrate this, we perform the flux transport calculation of a model of Li & Cao (2021) using our 2D code.Fig. 15 shows the result for the model with a scale height of  = 0.05 and a height of the corona surface at  ℎ = 0.2598 (see the left figure of Fig. 7 in their paper).In this model, the accretion speed drastically changes with height, and the advection speed in the corona (the region between the dotted and dashed lines in Fig. 15) reaches a maximum of approximately 35 times the equatorial plane velocity.With the coronal accretion, magnetic flux is more strongly dragged toward the centre, as expected (compare Figs. 1 and 15).The magnetic field distribution we obtained is generally consistent with their results.The grey-shaded areas correspond to the parameter spaces that satisfy all three conditions.Fig. 16 illustrates the parameter space of  ext that can achieve MAD in XRBs as a function of  m .The grey dashed line denotes the condition for  m given by equation ( 80).Namely,  m ≳ 2.8 for the given parameter set.The grey region indicates the range of  ext that can satisfy the possible three MAD conditions ( ext ≥  ext,1 ,  ext ≥  ext,2 , and equation 80).The top panel of Fig. 16 shows the result for a small MAD radius ( MAD = 2 g ).This case is a limiting case of MAD in the sense that the MAD size is just above the BH size.We find that the constraint from the plasma  is slightly stronger than but is similar to the constraint from the MAD parameter ( ext,1 ∼  ext,2 ).Therefore,  BH ∼ 40 would be a good criterion for MAD formation.The bottom panel of Fig. 16 displays the case with  MAD = 10 g .In this case,  ext,1 <  ext,2 in the allowed range of  m .Therefore, when the MAD size is much larger than the BH size, the constraint from the plasma  seems to be more important than that from the MAD parameter.
We note that the present calculation does not consider coronal accretion.If coronal accretion exists, the effective Prandtl number will increase.Therefore, we expect that the disc can reach the MAD state with a weaker external magnetic field.

SUMMARY
To understand the mechanism determining the magnetic field distribution in geometrically thick accretion discs around black holes, we analytically and numerically investigated magnetic flux transport in RIAF and SEAF.For the numerical study, we developed a time-dependent magnetic flux transport model in the axisymmetric spherical coordinate system (referred to as the 2D model in the text).Below, we summarize our findings.
It has been revealed that the magnetic field behaves diffusively due to the multi-dimensional effects arising from the thickness of the disc.The multi-dimensional effects of the thickness can be evaluated by the dimensionless parameter  eff (Sections 2.1.2and 3.1.2).In the conventional 1D models, these multi-dimensional effects cannot be accounted for, which may lead to an overestimation of the magnetic field strength and inclination (Section 3.1.2).However, when the disc is sufficiently thin (ℎ ≲ 0.1), the multi-dimensional effects are small, and the conventional 1D models are adequate for modelling (Section 3.2.1).
We demonstrated that the magnetic field distribution in a thick disc strongly depends on the magnetic Prandtl number  m (Section 3.2.2).This result points out the importance of local physical quantities in determining the global distribution of the magnetic field.
To investigate the driving condition for BP winds, we studied the dependence of the inclination of the magnetic field on the disc aspect ratio (ℎ) and magnetic Prandtl number ( m ) (Section 3.2.1 and 3.2.2).As thicker discs can drag the poloidal field more efficiently toward the centre, they form a more inclined field.In addition, discs with a higher  m can also do so.However, our 2D model revealed that the classical 1D model tends to overestimate the field inclination because of the lack of multi-dimensional effects.Therefore, in reality, the radial extent of the outflow driving regions could be narrower than predicted by the 1D model.
We investigated the condition for a RIAF to be MAD (Section 4.2).We derived the lower limit for the poloidal field strength at the radius of the outer edge of RIAF for a given  m .The diagram would be helpful for interpreting the 3D GRMHD simulations and observations.

Figure 1 .
Figure 1.The poloidal field structures in steady states (solid lines).The left and right panels show the results of RIAF and SEAF, respectively.The bottom panels are zoom-in images of the top panels.In the case of SEAF, the zoomed-in region approximately corresponds to the slim disc region.The figure shows the density distribution by colour to indicate the disc structure.The pressure scale height is shown as dashed-dotted lines.The results are for  m = 1.

Figure 2 .
Figure 2. Comparison between the 1D and 2D models.The red solid lines denote the 2D results, while the blue dashed lines indicate the 1D results.The top and bottom panels show the radial distribution of  ,mid and  ,surf , respectively.The horizontal grey lines in the top panels show the initial profile of  ,mid .
Fig.5compares the field inclination in the 1D and 2D models.As expected from Fig.4, the inclination is smaller in the 2D models.

Figure 3 .
Figure 3. Schematic pictures of the magnetic field structures just around the inner boundary in the 2D (top) and 1D (bottom) models.Note that in the 1D model, the poloidal field accumulates more significantly at the inner edge of the disc and has a larger inclination than in the 2D model.
Fig. 6 examines the distributions of   and   at  = 20 S in the 2D model.For SEAF, the profile in the slim region is shown.The solid lines represent the numerical results, and the dashed lines represent the lowest-order terms obtained by fitting the magnetic field distribution with a 15th-order function in terms of  (= /) (corresponding to

Figure 4 .
Figure 4. Comparison of the field inclination with the parameters  −1 and  −1 eff .The dashed lines denote the field inclination defined by the linear components of the field,  (1)  / (0)  .The dotted lines indicate the field inclination   /  | = .The red lines denote  −1 eff , while the blue lines show  −1 .The results are for 2D models with  m = 1.

Figure 5 .
Figure 5.Comparison of   /  | = between the 1D and 2D models.The red solid lines denote the results for the 2D models, while the blue dashed lines show the results for the 1D models.The red lines in this figure correspond to the dotted lines in Fig. 4. The horizontal grey lines denote the inclination angle of 30 • , above which the BP wind can blow.The results are for  m = 1.

Figure 6 .
Figure 6.The vertical profiles of   and   at  = 20  in the 2D models.The orange and blue solid lines denote   and   , respectively.The dashed lines indicate the linear terms of the field.The vertical black lines denote the height of the scale height,  =  ( = 1).For SEAF, the radius of  = 20  is inside the slim region.

Figure 8 .
Figure 8.The poloidal field structures of SEAF obtained by the 1D (blue lines) and 2D (red lines) models.The dashed line shows the photon trapping radius ( trap ≈ 85 S ).

Figure 9 .
Figure 9.Comparison of the 2D models to the 1D models with different approximation levels.The red lines denote the field profiles obtained by the 2D models.The blue lines are for the 1D models.The dotted blue lines show the results based on equation (14), while the dashed blue lines display the results based on equation (18).The latter ignores the effect of the radial gradient of the vertical field.The results are for  m = 1.
mid ln ≈ −1.93 tanh(0.279m ) for 2D model (62) 1.85 tanh(0.390m ) for 1D model (63) The slight deviation of the asymptotic value for  m → ∞ from −2 arises from the limited range of data used in the fitting, which is restricted to  m ≲ 20.The maximum errors of this fitting function are approximately 5 and 7% for the 2D and 1D models, respectively.Additionally, considering the analytical solution by Okuzumi et al. (2014), we performed fitting using a function with the coefficient of

Figure 13 .
Figure13.The relationship between  m and the power law index of  ,mid (−ln ,mid /ln) in RIAF (  = 1).The red circles and blue triangles denote the results of the 2D and 1D models, respectively.The red and blue solid lines indicate the fitting functions for the data of the 2D and 1D models, respectively.The fitting functions are shown in equations (64) and (65).The dashed black line represents the analytical solution for the power-law index in the 1D model(Okuzumi et al. 2014).

Figure 14 .
Figure 14.The  m dependence of the field inclination   /  | = .The grey lines denote the inclination angle of 30 • .

Figure 16 .
Figure16.The range of  ext that can achieve MAD in XRBs as a function of  m .The parameters are shown in the panels.The top panel shows the case in which the MAD radius  MAD is just 2 g , while the bottom panel corresponds to the case with a larger MAD radius,  MAD = 10 g .The red and blue lines denote the lower limits given by the constraints of the MAD parameter  BH and the plasma , respectively (see equations 76 and 78).The dashed vertical lines indicate the lower limit for  m (see the description about equation 80).The grey-shaded areas correspond to the parameter spaces that satisfy all three conditions.