© ESO, 2014

1. Introduction

Self-gravity in a disc orbiting a central star is significant if it is massive enough or/and very extended (Abramowicz et al. 1984; Mohanty et al. 2013). It can influence not only its rotation curve and trigger instabilities (including large-scale patterns and fragmentation), but also the accumulation of matter around the mid-plane. The perpendicular component of the field governs – together with the temperature – the disc thickness, internal pressure, density, etc. In a seminal paper, Paczynski (1978a) has studied the impact of vertical self-gravity on the hydrostatic equilibrium of a Keplerian disc. He assumed that the gradients of the potential are mainly vertical, making the Poisson equation fully integrable by analytical means, like for plane-parallel sheets (Goldreich & Lynden-Bell 1965). With this hypothesis, the field is linear with the surface density, which enables treating discs as concentric rings, independent of each other. Paczynski’s approximation was subsequently considered by many authors in the context of star formation and accretion in active galactic nuclei through “one zone” and bi-dimensional computations (Paczynski 1978b; Sakimoto & Coroniti 1981; Shore & White 1982; Cannizzo et al. 1982; Shlosman & Begelman 1987; Lantian & Xiaoci 1990; Cannizzo 1992; Cannizzo & Reiff 1992; Fukue & Sakamoto 1992; Huré et al. 1994, 2001; Liu et al. 1994, 1995; Mineshige & Umemura 1996, 1997; Mineshige et al. 1997; Pandey & Gupta 1998; Bardou et al. 1998; Huré 1998, 2000; Duschl et al. 2000; Vollmer & Beckert 2003; Mayer & Duschl 2005; Khajenabi & Duffy 2008). Depending on the viscosity prescription (Shakura & Sunyaev 1973; Pringle 1981), models predict the rise of the density along with the disc flattening and even pinching, which conditions are suited to the birth of global gravitational instabilities (Armitage 2010).

Paczynski’s approximation is expected to break down in certain zones where gradients of density and thickness are strong. These are typically disc edges, and possibly internal gaps and transition regions where matter experiences intrinsic changes (e.g., variation of opacity, or mean molecular weight; see Pickett et al. 1997; Huré 2000; Nelson et al. 2000). Like any open boundary, edges are important components of discs for geometrical and dynamical reasons (significant area, boundary layer, reservoir of angular momentum, escape of radiation, contraction/expansion into the ambient medium, etc.). In many models and grid-based hydrodynamical simulations, edges are often located outside the computational box, and are assumed not to be influential (Shakura & Sunyaev 1973; Collin-Souffrin & Dumont 1990; Dubrulle 1992), which requires specific boundary conditions (e.g. Paardekooper et al. 2011). Unfortunately, there is no easy way to determine with precision how the gravitational field varies at edges, except from a three-dimensional numerical approach. In this paper, we reconsider Paczynski’s formula for vertical self-gravity. Our main objective is to extend its domain of validity, in particular at the outer edge (i.e., large radii) where the field gradually falls and self-gravity can be significant. It is important to improve the understanding of those regions where most observations come from. We also caution that Paczynski’s formula is never valid in the common context of astrophysical discs, which questions models made previously (see references above). We reduce the technical difficulties by making a few main assumptions, namely: i) axial and mid-plane symmetries; ii) low aspect ratio; and iii) constant density profile in the direction perpendicular to the mid-plane. This framework is therefore well-suited to “vertically averaged, geometrically thin” discs. Thus, we derive an extended formula, discuss its limits and accuracy. In particular, by considering a few realistic density profiles, we conclude that only the local surface density is critical, not the vertical stratification. We then briefly discuss how the extended formula is expected to impact hydrostatic equilibrium. Clearly, it would be worth exploring its consequences by performing a detailed disc model, which may be the aim of a future article.

The paper is organised as follows: in Sect. 2, we recall the expression for the vertical field in an infinite plane-parallel medium. The monopole approximation is fairly irrelevant for cases of astrophysical interest. We show through a few concrete examples that Paczynski’s approximation fails i) close to disc edges; and ii) for centrally condensed discs. The estimate of the vertical field begins in Sect. 3 from the exact the integral expression valid for vertically homogeneous discs. We first simplify the axially symmetric kernel (which involves complete elliptic integrals) with a three-term fit (typical accuracy of 1%). We then use this fit to derive a formula appropriate for geometrically thin discs and valid for a wide range of radial density profiles and disc shapes. Next, we derive a criterion stating the condition of non-local contributions (i.e., the case of centrally condensed discs). We show that this result is compatible with Paczynski’s formula in the limit of infinite, plane-parallel sheets. In Sect. 4, we discuss edge effects. We consider power-law disc models, give a zero-order correction to Paczynski’s formula, check the assumption of vertical homogeneity, and stress to what extent self-gravity influences hydrostatic equilibrium and disc flaring. The last section is devoted to concluding remarks.

2. Local vs. non-local character of self-gravity

2.1. Plane-parallel sheets and Paczynski’s approximation

In cylindrical coordinates (R,θ,Z), the Poisson equation linking the gravitational potential ψ to the mass density ρ reads $\frac{\mathrm{1}}{\mathit{R}}\frac{\mathit{\partial }}{\mathit{\partial R}}\left(\mathit{R}\frac{\mathit{\partial \psi }}{\mathit{\partial R}}\right)\mathrm{+}\frac{\mathrm{1}}{{\mathit{R}}^{\mathrm{2}}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\psi }}{\mathit{\partial }{\mathit{\theta }}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\psi }}{\mathit{\partial }{\mathit{Z}}^{\mathrm{2}}}\mathrm{=}\mathrm{4}\mathit{\pi G\rho }\mathit{.}$(1)When the mass density is constant in any (R,θ)-plane, the radial and azimuthal gradients of the potential are strictly zero everywhere. The medium then consists of plane-parallel and homogeneous sheets (Goldreich & Lynden-Bell 1965). The configuration is depicted in Fig. 1a. Equation (1) can then be integrated in the vertical direction from the mid-plane Z = 0, leading to ${\frac{\mathit{\partial \psi }}{\mathit{\partial Z}}|}_{\mathit{Z}}\mathrm{-}{\frac{\mathit{\partial \psi }}{\mathit{\partial Z}}|}_{\mathrm{0}}\mathrm{=}\mathrm{4}\mathit{\pi G}\mathrm{\Sigma }\mathrm{\left(}\mathit{Z}\mathrm{\right)}\mathit{,}$(2)where $\mathrm{\Sigma }\mathrm{\left(}\mathit{Z}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{Z}}\mathit{\rho }\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{d}\mathit{z}$(3)is the cumulative surface density, i.e. the mass of gas per unit area from the midplane to Z. If the medium is symmetric with respect to the mid-plane, then ∂ψ/∂Z|₀ = 0, and the vertical component of gravity g_z = −∂_Zψ reads ${\mathit{g}}_{\mathit{Z}}\mathrm{\left(}\mathit{Z}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{4}\mathit{\pi G}\mathrm{\Sigma }\mathrm{\left(}\mathit{Z}\mathrm{\right)}\mathit{,}$(4)which result can also be established from the Gauss theorem (Armitage 2010). As long as the following conditions (hereafter “Paczynski’s approximation”; Paczynski 1978a)

Fig. 1 In the infinite, plane-parallel sheet model a), the magnitude of vertical self-gravity gZ is fully determined from the Poisson equation (or Gauss theorem) and is linear with the cumulative surface density. In real systems b), this is only correct in magnitude (Paczynski’s approximation): steep density gradients or high concentrations of mass in the inner disc produce significant deviations to this law.

Fig. 2 Vertical component gZ of self-gravity of a homogeneous disc computed from Eq. (14) compared with Paczynski’s formula versus the radius for Z = h (i.e., at the top of the disc). The disc parameters are ain/aout = 0.05, mass density ρ = 1, and semi-thickness h = 0.1aout (model 1 in Table 1). When the radial extent of the disc increases, gz tends to Paczynski’s formula, as indicated by the arrows, but edge effects are still present.

$\{\begin{array}{c}\\ \\ \\ \end{array}$(5)are satisfied, Eq.(2) roughly holds, but the field is only approximate, i.e. ${\mathit{g}}_{\mathit{Z}}\mathrm{\left(}\mathit{Z}\mathrm{\right)}\mathrm{\approx }\mathrm{-}\mathrm{4}\mathit{\pi G}\mathrm{\Sigma }\mathrm{\left(}\mathit{Z}\mathrm{\right)}\mathrm{\equiv }{\mathit{g}}_{\mathit{Z}}^{\mathrm{Pacz}\mathit{.}}\mathit{.}$(6)In fact, this expression is not very accurate in discs, which can generally not be considered as plane-parallel and homogeneous sheets. Discs are finite in size and are inhomogeneous. A simple proof is shown in Fig. 2 where we compare the “exact” vertical field¹ at the top of a thin, finite-size disc with constant density and constant semi-thickness h (disc model 1) with the value given by Eq. (6). The deviations are the strongest at the edges. Paczynski’s formula overestimates the field by a factor 2 typically close to the outer edge: there is no matter beyond the outer edge, which weakens the field. Obviously, the deviation is reduced at intermediate radii (i.e. far from the edges) when the disc axis ratio a_in/a_out decreases, but edge effects are always present. These concern a few local thicknesses in radius. Figure 3 shows the ratio $\frac{{\mathit{g}}_{\mathit{Z}}^{\mathrm{Pacz}\mathit{.}}}{{\mathit{g}}_{\mathit{Z}}}\mathrm{\equiv }\frac{\mathrm{1}}{\mathit{f}}$(7)for the three disc models listed in Table 1, including the fully homogeneous case considered before. Thus, f < 1 means that Paczynski’s formula overestimates the field. We see that Eq. (6) always overestimates the field at large radii.

2.2. Centrally condensed discs

Fig. 3 Component ratio 1/f defined by Eq. (7) for the three disc models defined in Table 1 including the case shown in Fig. 2 (the disc parameters are the same). When f < 1, Paczynski’s formula overestimates the field.

For centrally condensed discs (like model 3), a non-local effect appears. If the mass contained in the inner disc is large enough, then the outer disc feels the field of a central condensation, and so ${\mathit{g}}_{\mathit{Z}}\mathrm{\approx }\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{inner}\hspace{0.17em}\mathrm{disc}}}{{\mathit{R}}^{\mathrm{2}}}\mathrm{\times }\frac{\mathit{Z}}{\mathit{R}}\mathit{,}$(8)where M_inner disc is close to the total disc mass. Since this non-local field is not limited in magnitude, it can surpass the local contribution and Paczynski’s value as well, resulting in f > 1 in Eq. (7). To check this point, we considered model 3, but for different power-law density profiles, namely $\mathit{\rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{\propto }{\mathit{a}}^{\mathit{s}}\mathit{,}$(9)where the exponent s is a constant. Figure 4 displays 1/f versus R and for a continuum of the power-law exponents s ∈ [1, −4], which should include most cases of astrophysical interest. Here, the disc aspect ratio $\mathit{ϵ}\mathrm{=}\frac{\mathit{h}\mathrm{\left(}\mathit{a}\mathrm{\right)}}{\mathit{a}}$(10)is set to 0.1, and the axis ratio is the same as for Fig. 2. We clearly see that for s ≲ −2, the gravity of the inner disc cannot be ignored. This depends on ϵ and a_in/a_out.

Fig. 4 Component ratio 1/f defined by Eq. (7) for the same parameters as for model 3 defined in Table 1, but with different power-law exponents s in the range [−4,1] (contour lines every 0.1).

We conclude that the presence of edges makes the use of Paczynski’s approximation incorrect typically by a factor 0.5−2. Two contributions can play a role: the local disc, and the inner regions if these contain a large part of the disc mass. It would be erroneous to consider this range not wide enough to require a more detailed analysis. As quoted above, the hydrostatic equilibrium of a disc, through the equation $\begin{array}{ccc}{\mathit{g}}_{\mathit{Z}}\mathrm{=}\frac{\mathrm{1}}{\mathit{\rho }}{\mathit{\partial }}_{\mathit{Z}}\mathit{P}\mathrm{\approx }\mathrm{-}\frac{\mathrm{1}}{\mathit{h}}\mathrm{\times }{\frac{\mathit{P}}{\mathit{\rho }}|}_{\mathit{Z}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}}\mathit{,}& & \end{array}$(11)where P is the pressure, is very sensitive to gravity. We see that an error by a factor 2 in the field results in a factor $\sqrt{\mathrm{2}}$ on the gas temperature.

3. Vertical gravity of vertically homogeneous discs

The component we wish to calculate and compare with Eq. (6) is determined from the integral: ${\mathit{g}}_{\mathit{Z}}\mathrm{\left(}{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathit{G}{\$}_{\mathrm{disc}}\frac{\mathit{\rho }\mathrm{\left(}{{r}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{(}{r}\mathrm{-}{{r}}^{\mathrm{\prime }}\mathrm{)}\mathit{.}{{e}}_{\mathit{Z}}\mathit{a}\mathrm{d}\mathit{a}\mathrm{d}{\mathit{\theta }}^{\mathrm{\prime }}\mathrm{d}\mathit{z}}{\mathrm{|}{r}\mathrm{-}{{r}}^{\mathrm{\prime }}{\mathrm{|}}^{\mathrm{3}}}\mathit{,}$(12)where r′(a,θ′,z) refers to matter and e_Z is the unit vector. This integral is too complex to be estimated without additional hypothesis. A series expansion of the kernel, which leads to a component of the form $\mathrm{-}\mathit{Gm}\mathrm{\left(}\mathit{R}\mathrm{\right)}\frac{\mathit{Z}}{{\mathit{R}}^{\mathrm{3}}}\mathrm{+}\mathit{...}\mathit{,}$(13)where m(R) is the mass of gas from the inner edge of the disc to radius R, could be envisaged but it reproduces the field inaccurately, especially for an exponent power-law of s ≳ −4 (Mineshige & Umemura 1997; Faghei 2012). This is illustrated in Fig. 5, which displays the exact component g_z at the top of the disc divided by Eq. (13) at the same altitude versus the radius, for flared power-law discs with exponents ∈[−5, −1]. We see that this “monopole-type” approximation is valid only for very condensed discs with s ≲ −5. The above integral approach is therefore fully justified.

For axial and mid-plane symmetries, and for ∂_zρ = 0, the integrations over z and θ′ in Eq. (12) can be performed analytically. We have (Durand 1953; Huré 2005) ${\mathit{g}}_{\mathit{Z}}\mathrm{=}\mathrm{-}\mathrm{2}\mathit{G}{\mathrm{\int }}_{\mathit{a}}\mathit{\rho }\sqrt{\frac{\mathit{a}}{\mathit{R}}}\mathrm{\left[}{\mathit{k}}_{\mathrm{+}}\mathit{K}\mathrm{\right(}{\mathit{k}}_{\mathrm{+}}\mathrm{)}\mathrm{-}{\mathit{k}}_{\mathrm{-}}\mathit{K}\mathrm{(}{\mathit{k}}_{\mathrm{-}}\mathrm{\left)}\mathrm{\right]}\mathrm{d}\mathit{a}\mathit{,}$(14)where ${\mathit{k}}_{\mathrm{\pm }}\mathrm{=}\frac{\mathrm{2}\sqrt{\mathit{aR}}}{\sqrt{\mathrm{(}\mathit{a}\mathrm{+}\mathit{R}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{Z}\mathrm{\mp }\mathit{h}{\mathrm{)}}^{\mathrm{2}}}}\mathrm{\in }\mathrm{\left[}\mathrm{0}\mathit{,}\mathrm{1}\mathrm{\right]}$(15)is the modulus of the complete elliptic integral of the first kind K(k), and ± h is the altitude of the top/bottom surface of the disc. Both ρ and h can depend on the radius a. We note that the integrand in Eq. (14) is zero at the mid-plane (i.e., for Z = 0), regardless of a and R because k₊ = k₋. Off the mid-plane, the function kK(k) exhibits a peak at a = R; this peak rises as Z → ± h and is even infinite when Z = ± h, which corresponds to k = 1 (the elliptic integral is logarithmically singular; see the appendix). This peak lends weight to the mass density distributed around R and is therefore responsible for the local character of self-gravity. This function also has broad wings, which are important if ρ is large at radii a ≪ R or a ≫ R. This is a non-local effect.

Fig. 5 Ratio of the vertical component of the disc gravity to the monopole value versus the radius for model 3, but with different power law exponents s of the density. Same edges as for Fig. 2.

Fig. 6 Function kK(k) + lnk′ versus k and its linear fit by Eq. (17). The absolute error is also shown (dotted line).

3.1. Three-term fit of the kernel

To estimate Eq. (14) by analytical means, we must replace the special function by a series. Expanding K(k) gives at the lowest order (Gradshteyn & Ryzhik 2007) $\{\begin{array}{c}\\ \\ \\ \end{array}$(16)where ${\mathit{k}}^{\mathrm{\prime }}\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}{\mathit{k}}^{\mathrm{2}}}$ is the complementary modulus. These two asymptotic forms do not connect well at intermediate modulus, however. The five-term and nine-term fits given in Abramowitz & Stegun (1964) are not very practical and are too accurate for the present purpose. We have noticed that the function kK(k) + lnk′, when plotted versus k, is very close to linear. This is shown in Fig. 6. We therefore considered the following three-term fit: $\mathit{kK}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathit{C}}_{\mathrm{0}}\mathrm{+}{\mathit{C}}_{\mathrm{1}}\mathit{k}\mathrm{-}\ln {\mathit{k}}^{\mathrm{\prime }}\mathrm{+}\mathrm{err}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}$(17)where err(k) is the error. The result of a least-squares regression gives C₀ = 0.0246 and C₁ = 1.3371. Figure 6 shows the associated errors. The area below the curve is remarkably well conserved because we have $\begin{array}{ccc}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\mathrm{err}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{d}\mathit{k}\mathrm{=}\mathrm{1}\mathrm{-}\left({\mathit{C}}_{\mathrm{0}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{C}}_{\mathrm{1}}\mathrm{-}\ln \mathrm{2}\mathrm{+}\mathrm{1}\right)\mathrm{\approx }\mathrm{2}\mathrm{\times }{\mathrm{10}}^{-6}\mathit{.}& & \end{array}$(18)This approximation is certainly not perfect², but it minimises the subsequent analytical treatment and provides good results since | err(k)| ≲ 0.015 in the whole domain k ∈ [0,1]. Using Eq. (17), the integrand in Eq. (14) is ${\mathit{k}}_{\mathrm{+}}\mathit{K}\mathrm{\left(}{\mathit{k}}_{\mathrm{+}}\mathrm{\right)}\mathrm{-}{\mathit{k}}_{\mathrm{-}}\mathit{K}\mathrm{\left(}{\mathit{k}}_{\mathrm{-}}\mathrm{\right)}\mathrm{\approx }\mathrm{-}{\mathit{C}}_{\mathrm{1}}\mathrm{(}{\mathit{k}}_{\mathrm{+}}\mathrm{-}{\mathit{k}}_{\mathrm{-}}\mathrm{)}\mathrm{-}\ln \frac{{\mathit{k}}_{\mathrm{+}}^{\mathrm{\prime }}}{{\mathit{k}}_{\mathrm{-}}^{\mathrm{\prime }}}\mathit{,}$(19)with an accuracy of 1% or less. This is proved in Fig. 7, which gives the relative error made on the vertical component using this three-term fit instead of the exact kernel. The conditions are the same as for Fig. 4. Note that coefficient C₀ does not play a role in the following.

Fig. 7 Relative deviation between the vertical component gZ computed from the three-term fit of the function kK(k) and from the exact kernel kK(k). The conditions are the same as for Fig. 4 (contour dotted lines every 0.001).

3.2. Application to geometrically thin discs

At this level, we can estimate the vertical field for any type of vertically homogeneous medium by inserting the three-term fit in Eq. (14). We find ${\mathit{g}}_{\mathit{Z}}\mathrm{=}\mathrm{+}\mathrm{2}\mathit{G}{\mathrm{\int }}_{\mathit{a}}\mathit{\rho }\sqrt{\frac{\mathit{a}}{\mathit{R}}}\left[{\mathit{C}}_{\mathrm{1}}\mathrm{(}{\mathit{k}}_{\mathrm{+}}\mathrm{-}{\mathit{k}}_{\mathrm{-}}\mathrm{)}\mathrm{+}\ln \frac{{\mathit{k}}_{\mathrm{+}}^{\mathrm{\prime }}}{{\mathit{k}}_{\mathrm{-}}^{\mathrm{\prime }}}\right]\mathrm{d}\mathit{a}\mathit{,}$(20)which only avoids manipulating special functions (see the appendix for the numerical evaluation). We can further simplify the formula by assuming $\mathrm{(}\mathit{Z}\mathrm{\mp }\mathit{h}{\mathrm{)}}^{\mathrm{2}}\mathrm{\ll }\mathrm{(}\mathit{a}\mathrm{+}\mathit{R}{\mathrm{)}}^{\mathrm{2}}\mathit{,}$(21)which is suited to discs such that ϵ² ≪ 1 and close neighbourhood. Equation (20) can then be written in the form ${\mathit{g}}_{\mathit{Z}}\mathrm{\approx }\mathrm{-}\mathit{G}{\mathrm{\int }}_{\mathit{a}}\mathit{\rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}{\mathit{w}}_{\mathrm{1}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{d}\mathit{a}\mathrm{-}\mathit{G}{\mathrm{\int }}_{\mathit{a}}\mathit{\rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}{\mathit{w}}_{\mathrm{2}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{d}\mathit{a}\mathit{,}$(22)where $\{\begin{array}{c}\\ \\ \\ \end{array}$(23)and $\mathit{m}\mathrm{=}\frac{\mathrm{2}\sqrt{\mathit{aR}}}{\mathrm{(}\mathit{a}\mathrm{+}\mathit{R}\mathrm{)}}\mathrm{\in }\mathrm{\left[}\mathrm{0}\mathit{,}\mathrm{1}\mathrm{\right]}\mathit{.}$(24)Figure 8 displays w₁ and w₂ versus a/a_out in a typical case. We see that w₁ dominates in absolute close to the outer edge. This is true for any altitude between the mid-plane and the surface Z = ± h. For a ≪ R, however, w₁ and w₂ are similar in magnitude, meaning that w₂ can probably not be neglected for centrally condensed discs.

Fig. 8 Two weighting functions w1 and w2 plotted versus the normalised radius a/aout for model 2 (see Table 1), ϵ = 0.1, $\mathit{Z}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}\mathit{h}$ and R/aout = 0.75.

3.3. Local vs. non-local contributions: a criterion

As stated above, the local contribution of self-gravity is due to the peak of the function w₁(a) at a = R, while the possible role of the inner disc is due to the left wing of w₁. We can derive a limit between these two regimes by equating the local and non-local contributions. To avoid the possible singularity at a = R, we integrate w₁(a) over one semi-thickness leftward and rightward to R, and then divide by 2h. To compute this radial average, we assume that ρ(a) and h(a) do not change significantly in the integration range a ∈ [R − h(a),R + h(a)]. After some calculation, we find $\begin{array}{ccc}\frac{\mathrm{1}}{\mathrm{2}\mathit{h}}{\mathrm{\int }}_{\mathrm{-}\mathit{h}}^{\mathit{h}}{\mathit{w}}_{\mathrm{1}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{d}\mathit{a}\mathrm{=}\ln \mathrm{4}\mathrm{+}\mathrm{4}\hspace{0.17em}\mathrm{atan}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\equiv }\mathit{A}\mathrm{\approx }\mathrm{3.241}\mathit{,}& & \end{array}$(25)at Z = h. For a ≪ R, still assuming weak radial gradients for ρ(a) and h(a), we have $\frac{\mathrm{1}}{\mathrm{2}\mathit{h}}{\mathrm{\int }}_{\mathrm{-}\mathit{h}}^{\mathit{h}}{\mathit{w}}_{\mathrm{1}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{d}\mathit{a}\mathrm{\approx }{\mathit{w}}_{\mathrm{2}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathit{.}$(26)Provided ρ(a) decreases monotonically outward, the local contribution exceeds the non-local (i.e., inner disc) contribution by a factor α if $\mathit{A\rho }\mathrm{\left(}\mathit{R}\mathrm{\right)}\gtrsim \mathit{\alpha \rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}{\mathit{w}}_{\mathrm{1}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathit{,}$(27)which leads to the criterion $\frac{\mathit{A}}{\mathrm{4}\mathit{\alpha }}\frac{\mathit{\rho }\mathrm{\left(}\mathit{R}\mathrm{\right)}}{\mathit{\rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}}\mathrm{\ge }\sqrt{\frac{\mathit{a}}{\mathit{R}}}\frac{\mathit{h}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathit{h}\mathrm{\left(}\mathit{a}\mathrm{\right)}}{{\mathit{R}}^{\mathrm{2}}}\mathrm{·}$(28)A safe criterion is then obtained with α ~ 10 (see Sect. 4 for power-law discs).

3.4. Recovering Paczynski’s formula?

To verify that Eq. (22) is compatible with Paczynski’s formula, we must include appropriate assumptions, namely ρ(a) = cst and h(a) = cst. Next, we omit w₂, which is much lower than w₁ in magnitude. We then have ${\mathit{g}}_{\mathit{Z}}\mathrm{\approx }\mathrm{-}\mathit{G\rho }\mathrm{\left(}\mathit{R}\mathrm{\right)}{\mathrm{\int }}_{\mathit{a}}{\mathit{w}}_{\mathrm{1}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{d}\mathit{a}\mathit{.}$(29)Setting $\sqrt{\frac{\mathit{a}}{\mathit{R}}}\mathrm{\approx }\mathrm{1}\mathrm{+}\frac{\mathit{a}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{R}}{\mathrm{2}\mathit{R}}$, the order zero term is $\begin{array}{ccc}\mathrm{\int }{\mathit{w}}_{\mathrm{1}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{d}\mathit{a}& & \mathrm{\approx }\mathrm{\int }\ln \frac{\mathrm{(}\mathit{a}\mathrm{-}\mathit{R}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{Z}\mathrm{-}\mathit{h}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{(}\mathit{a}\mathrm{-}\mathit{R}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{Z}\mathrm{+}\mathit{h}{\mathrm{)}}^{\mathrm{2}}}\mathrm{d}\mathit{a}\\ & \mathrm{=}& \mathrm{(}\mathit{a}\mathrm{-}\mathit{R}\mathrm{)}\ln \frac{\mathrm{(}\mathit{a}\mathrm{-}\mathit{R}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{Z}\mathrm{-}\mathit{h}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{(}\mathit{a}\mathrm{-}\mathit{R}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{Z}\mathrm{+}\mathit{h}{\mathrm{)}}^{\mathrm{2}}}\\ & & \mathrm{+}\mathrm{2}\left[\mathrm{(}\mathit{Z}\mathrm{-}\mathit{h}\mathrm{)}\hspace{0.17em}\mathrm{atan}\frac{\mathit{a}\mathrm{-}\mathit{R}}{\mathit{Z}\mathrm{-}\mathit{h}}\mathrm{-}\mathrm{(}\mathit{Z}\mathrm{+}\mathit{h}\mathrm{)}\hspace{0.17em}\mathrm{atan}\frac{\mathit{a}\mathrm{-}\mathit{R}}{\mathit{Z}\mathrm{+}\mathit{h}}\right]\mathrm{·}\end{array}$(30)We now take R − nh and R + nh as integral bounds, where n is a positive number. The contribution of matter located leftward of R is $\begin{array}{ccc}{\mathrm{\int }}_{\mathit{a}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{R}}^{\mathit{a}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{R}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{nh}}{\mathit{w}}_{\mathrm{1}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{d}\mathit{a}\mathrm{\approx }\hspace{0.17em}& & \mathrm{4}\mathit{Z}\frac{\mathit{n}}{\mathrm{1}\mathrm{+}{\mathit{n}}^{\mathrm{2}}}\mathrm{-}\mathrm{2}\mathrm{(}\mathit{h}\mathrm{-}\mathit{Z}\mathrm{)}\hspace{0.17em}\mathrm{atan}\frac{\mathit{nh}}{\mathit{h}\mathrm{-}\mathit{Z}}\\ & & \mathrm{+}\mathrm{2}\mathrm{(}\mathit{h}\mathrm{+}\mathit{Z}\mathrm{)}\hspace{0.17em}\mathrm{atan}\frac{\mathit{nh}}{\mathit{h}\mathrm{+}\mathit{Z}}\mathrm{·}\end{array}$(31)Assuming n ≫ 1, this gives ${\mathrm{\int }}_{\mathit{a}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{R}}^{\mathit{a}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{R}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathit{nh}}{\mathit{w}}_{\mathrm{1}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{d}\mathit{a}\mathrm{\approx }\mathrm{2}\mathit{\pi Z}\mathit{.}$(32)Since the integrand in Eq. (30) is odd in a − R, the right-side contribution is the same. We finally find ${\mathit{g}}_{\mathit{Z}}\mathrm{\approx }\mathrm{-}\mathit{G\rho }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{\times }\mathrm{2}\mathit{\pi Z}\mathrm{\times }\mathrm{2}\mathit{,}$(33)which is equivalent to Eq. (6) for a vertically homogeneous slab. This demonstration clearly explains why the presence of an edge suppresses one wing in w₁(a), which diminishes the vertical field typically by a factor ~2.

3.5. Is the vertical stratification critical?

The results presented so far are only valid for vertically homogeneous systems. To determine the sensitivity to this hypothesis, which may seem too restrictive for applications, we considered two density profiles: $\mathit{\rho }\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\times }\{\begin{array}{c}\\ \\ \\ \end{array}$(34)The first one is a classical Gaussian distribution, typical of vertically isothermal discs, and the second one its quadratic version, which, in contrast, has a vanishing density at the surface. Both have the same cumulative surface density at z = h, i.e. ρ₀h, so the comparison of the vertical fields with the homogeneous case is straightforward. The results are displayed in Fig. 9 for model 1, and in Fig. 10 for model 3 (with a zoom at the outer edge). The deviations caused by a non-uniform vertical stratification are of the order of one percent or weaker. Note that in the Gaussian case, matter located at z ≥ h (~5% of the total mass per unit area) has almost no influence. We conclude that all formulae derived here assuming ∂_zρ = 0 at any radius are relatively robust.

Fig. 9 Same legend and same conditions as for Fig. 3 (disc model 1), but for the Gaussian and quadratic profiles. Deviations with respect to the vertically homogeneous case are globally of the order of one percent or weaker.

Fig. 10 Same as for Fig. 9 but for model 3.

4. Measuring edge effects

4.1. Zero-order correction to Paczynski’s formula

It is difficult to deduce an algebraic formula for g_Z at the top of the disc, even when the local contribution is dominant. However, we can derive the order zero by considering the development above. If we define the variable $\mathit{\varpi }\mathrm{=}\frac{\mathrm{|}\mathit{a}\mathrm{-}\mathit{R}\mathrm{|}}{\mathrm{2}\mathit{h}\mathrm{\left(}\mathit{R}\mathrm{\right)}}\mathit{,}$(35)which measures the distance to the actual radius R in units of local thickness 2h(R), then Eq. (30) gives $\mathrm{\int }{\mathit{w}}_{\mathrm{1}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{d}\mathit{a}\mathrm{\approx }\mathrm{4}\mathit{h}\mathrm{\left(}\mathit{R}\mathrm{\right)}\left[\mathit{\varpi }\ln \frac{\mathit{\varpi }}{\sqrt{\mathrm{1}\mathrm{+}{\mathit{\varpi }}^{\mathrm{2}}}}\mathrm{-}\hspace{0.17em}\mathrm{atan}\mathit{\varpi }\right]\mathit{,}$(36)where we set Z = h(R). The vertical field at the surface of the disc is then given for R ∈ [a_in,a_out] by the relation ${\mathit{g}}_{\mathit{Z}}\mathrm{\approx }{\mathit{g}}_{\mathit{Z}}^{\mathrm{Pacz}\mathit{.}}\mathrm{\times }{\mathit{f}}_{\mathrm{edges}}\mathit{,}$(37)where

$\begin{array}{ccc}& & {\mathit{f}}_{\mathrm{edges}}\mathrm{=}{\mathit{f}}_{\mathrm{in}}\mathrm{+}{\mathit{f}}_{\mathrm{out}}\mathit{,}\\ & & \{\begin{array}{c}\\ & \\ & \end{array}\end{array}$and ϖ_in is for a_in, and ϖ_out is for a_out. Clearly, f_edges is a function of the radius R. Figure 11 shows the accuracy of this zero-order approximation (i.e. the relative error made on f_edges), in the same conditions as Fig. 4. We see that the error does not exceed 10% ~ ϵ in the domain where the local contribution dominates (i.e. for s ≳ s₀). If the mass stored in the inner discs does not play a role, we have $\mathit{\pi }{\mathit{f}}_{\mathrm{in}}\mathrm{\approx }\frac{\mathit{\pi }}{\mathrm{2}}$, and so ${\mathit{f}}_{\mathrm{edges}}\mathrm{\approx }\frac{\mathrm{1}}{\mathrm{2}}\mathrm{+}\hspace{0.17em}\mathrm{atan}{\mathit{\varpi }}_{\mathrm{out}}\mathrm{-}{\mathit{\varpi }}_{\mathrm{out}}\ln \frac{{\mathit{\varpi }}_{\mathrm{out}}}{\sqrt{\mathrm{1}\mathrm{+}{\mathit{\varpi }}_{\mathrm{out}}^{\mathrm{2}}}}\mathrm{·}$(40)

Fig. 11 Relative error on the vertical field when computed from Eq. (37) (i.e. the zero-order extension of Paczynski’s formula) for the same parameters as for model 3 defined in Table 1, but for different power-law exponents s in the range [−4,1] (contour lines every 0.01).

4.2. Criterion for power-law discs

For power-law discs where ρ(a) is defined by Eq. (9) and which have a constant aspect ratio ϵ, Eq. (28) becomes $\mathit{s}\gtrsim \mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\mathrm{+}\frac{\ln \frac{\mathrm{4}\mathit{\alpha }{\mathit{ϵ}}^{\mathrm{2}}}{\mathit{\pi }}}{\ln \frac{\mathit{R}}{{\mathit{a}}_{\mathrm{in}}}}\mathrm{·}$(41)For the conditions of Fig. 8, the criterion is s ≳ −2.26 for α = 10, i.e., the local contribution is higher than the inner edge contribution by a factor 10. This is s ≳ −3.11 for α = 1. We see in Fig. 12, where we have plotted w₁(a) × ρ(a) versus a, that the criterion is reliable. To render the criterion robust and independent of R, we set R = a_out and α = 10. Then, the criterion reads $\begin{array}{ccc}\mathit{s}\mathrm{\ge }\mathrm{-}\mathrm{1.5}\mathrm{+}\frac{\mathrm{2}\log \mathit{ϵ}\mathrm{+}\mathrm{1.1}}{\log {\mathit{a}}_{\mathrm{out}}\mathit{/}{\mathit{a}}_{\mathrm{in}}}\mathrm{·}& & \end{array}$(42)

Fig. 12 Function ρ(a) × w1(a) plotted versus the normalised radius a/aout in the same condition as for Fig. 8, but for two different exponents: s = −2.26 (which corresponds to α = 10) and for s = −3.11 (which corresponds to α = 1).

4.3. Are self-gravitating discs still flaring?

How is the disc shape eventually modified in the self-gravitating region? It is not possible to give a firm answer to this question without considering a detailed disc model including balance and transport equations. As quoted in the introduction, models usually show that the thickness is reduced when Paczynski’s formula is applied. The fact that g_Z/Σ ∝ f_edges decreases close to the edge implies a thickening of the disc. Figure 13 compares the local flaring angle ϵ(R) in a disc close to the outer edge when using Paczynski’s formula and when using our new expression. In this simple example, the temperature ~P/ρ is kept constant and we just compare h/R resulting from Eq. (11) with two different formulae for g_Z. We also include the contribution of a central star such that the disc mass is 10% of the central mass. We recognize the disc pinching, typical of Paczysnki’s formula. This property is not found by considering the edge effect. On the contrary, the disc is expected to exhibit a certain flaring, less than in a Keplerian case, but still a flaring. Obvisouly, a detailed disc model is required to take into account any decrease of Σ and T outward.

Fig. 13 Flaring angle ϵ(R) close to the disc outer edge for various models of vertical field including a central star ten times more massive than the disc. The disc parameters other than the semi-thickness h are the same as for Fig. 2 (see text).

5. Conclusion

We have revisited Paczynski’s formula for vertical self-gravity to extend its domain of validity. The primary aim was to provide a reliable analytical tool to better investigate the structure of outer parts of discs, typically through “one-zone” models. We demonstrated that the classical approximation is erroneous typically by a factor 0.5−2 in the conditions of common use, which can have a noticeable influence on the disc equilibrium. For very centrally condensed discs, it can even gravely underestimate gravity. Our formula works for a wide family of radial profiles and shapes, and its accuracy is better than 1% in any case. Depending on the disc thickness (small or large), different formulae can be implemented in models. Interestingly enough, the details of the vertical density structure are not important, only the surface density at the top of the disc is critical. For geometrically thin discs, however, the vertical field is, at zero order at least, fully algebraic. If necessary, high orders can be considered.

There are a few interesting consequences for the decrease of vertical self-gravity close to the disc edge. As noted, the disc could not be pinched as previously concluded (see references in the introduction), but could still flare. Obviously, this claim must be adjusted with the necessary transition region from the disc to the ambient/external medium through a decrease of the surface density and temperature. Nevertheless, the edge effect may render the disc thicker than thought, and subsequently less dense. This raises two interesting questions. Can the outer disc still be illuminated by a central star if present (e.g. Dullemond & Dominik 2004)? Can the decrease of the density postpone the triggering of global instabilities and affect Toomre’s criterion? There may be a range of disc masses, high enough to be subject to vertical self-gravity, but low enough to prevent global instabilities.

Appendix: Estimation of ${{}^{\mathrm{\int }}}_{\mathit{a}}\sqrt{\frac{\mathit{a}}{\mathit{R}}}\mathit{\rho }\mathrm{\left(}\mathit{a}\mathit{)}\mathit{kK}\mathrm{\right(}\mathit{k}\mathrm{)}\mathrm{d}\mathit{a}$ through integral splitting

We can numerically estimate Eq. (14) with accuracy by using integral splitting (Pierens & Huré 2004; Huré 2005). In the present case, we extracted the log. singularity from the function K, leaving its regular part, denoted K^∗. Then, we decomposed the integrand as follows:

$\begin{array}{ccc}& & \sqrt{\frac{\mathit{a}}{\mathit{R}}}\mathit{\rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathit{kK}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\sqrt{\frac{\mathit{a}}{\mathit{R}}}\mathit{\rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathit{k}\left[{\mathit{K}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{-}\ln {\mathit{k}}^{\mathrm{\prime }}\right]\\ & & \mathrm{=}\sqrt{\frac{\mathit{a}}{\mathit{R}}}\mathit{\rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathit{k}{\mathit{K}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{-}\sqrt{\frac{\mathit{a}}{\mathit{R}}}\mathit{k\rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}\ln {\mathit{k}}^{\mathrm{\prime }}\\ & & \mathrm{=}\sqrt{\frac{\mathit{a}}{\mathit{R}}}\mathit{\rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathit{k}{\mathit{K}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{-}\left[\sqrt{\frac{\mathit{a}}{\mathit{R}}}\mathit{k\rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{-}\mathit{\rho }\mathrm{\left(}\mathit{R}\mathrm{\right)}\right]\ln {\mathit{k}}^{\mathrm{\prime }}\mathrm{-}\mathit{\rho }\mathrm{\left(}\mathit{R}\mathrm{\right)}\ln {\mathit{k}}^{\mathrm{\prime }}\\ & & \mathrm{=}\sqrt{\frac{\mathit{a}}{\mathit{R}}}\mathit{\rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathit{k}{\mathit{K}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{-}\left[\sqrt{\frac{\mathit{a}}{\mathit{R}}}\mathit{k\rho }\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{-}\mathit{\rho }\mathrm{\left(}\mathit{R}\mathrm{\right)}\right]\ln {\mathit{k}}^{\mathrm{\prime }}\\ & & \mathrm{-}\mathit{\rho }\mathrm{\left(}\mathit{R}\mathrm{\right)}\begin{array}{c}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\ln {\mathit{k}}^{\mathrm{\prime }}\mathrm{-}\ln \sqrt{\frac{\mathrm{(}\mathit{a}\mathrm{-}\mathit{R}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{[}\mathit{Z}\mathrm{-}\mathit{h}\mathrm{\left(}\mathit{R}\mathrm{\right)}{\mathrm{]}}^{\mathrm{2}}}{\mathrm{(}\mathit{a}\mathrm{+}\mathit{R}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{[}\mathit{Z}\mathrm{-}\mathit{h}\mathrm{\left(}\mathit{R}\mathrm{\right)}{\mathrm{]}}^{\mathrm{2}}}}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\\ & & \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }\mathrm{\left(}\mathit{R}\mathrm{\right)}\ln \frac{\mathrm{(}\mathit{a}\mathrm{-}\mathit{R}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{[}\mathit{Z}\mathrm{-}\mathit{h}\mathrm{\left(}\mathit{R}\mathrm{\right)}{\mathrm{]}}^{\mathrm{2}}}{\mathrm{(}\mathit{a}\mathrm{+}\mathit{R}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{[}\mathit{Z}\mathrm{-}\mathit{h}\mathrm{\left(}\mathit{R}\mathrm{\right)}{\mathrm{]}}^{\mathrm{2}}}\mathrm{·}\end{array}$

We see that the first three terms are always finite, even when k = 1. The last term is possibly singular, but it contains terms of the form cst × ln(x² + cst), which are integrable by analytical means. This method also works to estimate Eqs. (20) and (22).

Acknowledgments

It is a pleasure to thank the referee, W.J. Duschl, for valuable comments and suggestions, in particular about vertical stratification.

References

All Tables

All Figures

Fig. 1 In the infinite, plane-parallel sheet model a), the magnitude of vertical self-gravity gZ is fully determined from the Poisson equation (or Gauss theorem) and is linear with the cumulative surface density. In real systems b), this is only correct in magnitude (Paczynski’s approximation): steep density gradients or high concentrations of mass in the inner disc produce significant deviations to this law.

In the text

Fig. 2 Vertical component gZ of self-gravity of a homogeneous disc computed from Eq. (14) compared with Paczynski’s formula versus the radius for Z = h (i.e., at the top of the disc). The disc parameters are ain/aout = 0.05, mass density ρ = 1, and semi-thickness h = 0.1aout (model 1 in Table 1). When the radial extent of the disc increases, gz tends to Paczynski’s formula, as indicated by the arrows, but edge effects are still present.

In the text

	Fig. 3 Component ratio 1/f defined by Eq. (7) for the three disc models defined in Table 1 including the case shown in Fig. 2 (the disc parameters are the same). When f < 1, Paczynski’s formula overestimates the field.
In the text

	Fig. 4 Component ratio 1/f defined by Eq. (7) for the same parameters as for model 3 defined in Table 1, but with different power-law exponents s in the range [−4,1] (contour lines every 0.1).
In the text

	Fig. 5 Ratio of the vertical component of the disc gravity to the monopole value versus the radius for model 3, but with different power law exponents s of the density. Same edges as for Fig. 2.
In the text

	Fig. 6 Function kK(k) + lnk′ versus k and its linear fit by Eq. (17). The absolute error is also shown (dotted line).
In the text

	Fig. 7 Relative deviation between the vertical component gZ computed from the three-term fit of the function kK(k) and from the exact kernel kK(k). The conditions are the same as for Fig. 4 (contour dotted lines every 0.001).
In the text

	Fig. 8 Two weighting functions w1 and w2 plotted versus the normalised radius a/aout for model 2 (see Table 1), ϵ = 0.1, $\mathit{Z}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}\mathit{h}$ and R/aout = 0.75.
In the text

	Fig. 9 Same legend and same conditions as for Fig. 3 (disc model 1), but for the Gaussian and quadratic profiles. Deviations with respect to the vertically homogeneous case are globally of the order of one percent or weaker.
In the text

	Fig. 10 Same as for Fig. 9 but for model 3.
In the text

	Fig. 11 Relative error on the vertical field when computed from Eq. (37) (i.e. the zero-order extension of Paczynski’s formula) for the same parameters as for model 3 defined in Table 1, but for different power-law exponents s in the range [−4,1] (contour lines every 0.01).
In the text

	Fig. 12 Function ρ(a) × w1(a) plotted versus the normalised radius a/aout in the same condition as for Fig. 8, but for two different exponents: s = −2.26 (which corresponds to α = 10) and for s = −3.11 (which corresponds to α = 1).
In the text

	Fig. 13 Flaring angle ϵ(R) close to the disc outer edge for various models of vertical field including a central star ten times more massive than the disc. The disc parameters other than the semi-thickness h are the same as for Fig. 2 (see text).
In the text