Interior potential of a toroidal shell from pole values

TL;DR

This paper derives an analytical formula for the gravitational potential inside a toroidal shell, develops series approximations, and confirms their accuracy, advancing understanding of toroidal gravitating systems and related fields.

Contribution

It provides the first exact formula for the interior potential of a toroidal shell at the pole, along with analytical derivatives and series expansions for improved modeling.

Findings

Exact potential formula involving elliptic integrals

Series expansions accurately approximate the interior potential

Potential is non-uniform inside the toroidal cavity, unlike ellipsoids.

Abstract

We have investigated the toroidal analog of ellipsoidal shells of matter, which are of great significance in Astrophysics. The exact formula for the gravitational potential $\Psi(R,Z)$ of a shell with a circular section at the pole of toroidal coordinates is first established. It depends on the mass of the shell, its main radius and axis-ratio $e$ (i.e. core-to-main radius ratio), and involves the product of the complete elliptic integrals of the first and second kinds. Next, we show that successive partial derivatives $\partial^{n +m} \Psi/\partial_{R^n} \partial_{Z^m}$ are also accessible by analytical means at that singular point, thereby enabling the expansion of the interior potential as a bivariate series. Then, we have generated approximations at orders $0$, $1$, $2$ and $3$, corresponding to increasing accuracy. Numerical experiments confirm the great reliability of the approach, in particular for small-to-moderate axis ratios ($e^2 \lesssim 0.1$ typically). In contrast with the ellipsoidal case (Newton's theorem), the potential is not uniform inside the shell cavity as a consequence of the curvature. We explain how to construct the interior potential of toroidal shells with a thick edge (i.e. tubes), and how a core stratification can be accounted for. This is a new step towards the full description of the gravitating potential and forces of tori and rings. Applications also concern electrically-charged systems, and thus go beyond the context of gravitation.

Equations

$$\Psi (R,Z)=-2G\lambda \sqrt{\frac{a}{R}}k\mathbf{K}(k),$$

$$\mathbf{K}(k)={\int }_0^{\frac{\pi }{2}}\frac{dx}{\sqrt{1-k^2{\sin }^{2}x}}$$

$$k=\frac{2\sqrt{aR}}{\sqrt{(a+R{)}^2+(Z-z{)}^2}}\in [0,1],$$

$${\left[a-R\left(\frac{1+{k^{\prime }}^2}{k^2}\right)\right]}^2+(z-Z{)}^2={\left(\frac{2R{k}^{\prime }}{k^2}\right)}^2,$$

$$\{\begin{array}{l}{R}_c=R\frac{1+{k^{\prime }}^2}{k^2},\\ Z_c=Z,\end{array}$$

$$b=R\frac{2k^{\prime }}{k^2}.$$

$$\{\begin{array}{l}R=R_{\mathrm{p}}\frac{\sinh \eta }{\cosh \eta -\cos \zeta },\\ Z=R_{\mathrm{p}}\frac{\sin \zeta }{\cosh \eta -\cos \zeta },\end{array}$$

$${\left(a-R_{\mathrm{p}}\coth \eta \right)}^2+z^2={\left(\frac{R_{\mathrm{p}}}{\sinh \eta }\right)}^2,$$

$$\{\begin{array}{l}R\equiv R_{\mathrm{p}},\\ Z\equiv 0,\\ k^2=\frac{2\sinh \eta }{\cosh \eta +\sinh \eta },\\ \mathrm{or}\coth \eta =\frac{1+{k^{\prime }}^2}{k^2}.\end{array}$$

$$\{\begin{array}{l}a=R_c+b\cos \theta ,\\ z=b\sin \theta ,\end{array}$$

$$\Psi (R,Z)=-2G{\int }_0^{2\pi }\Sigma \sqrt{\frac{a}{R}}k\mathbf{K}(k)bd\theta ,$$

$$\psi (R_{\mathrm{p}},0)$$

$$\equiv {\Psi }_p,$$

$${\int }_0^{2\pi }\sqrt{a}d\theta$$

$$=4\sqrt{R_{\mathrm{p}}}\frac{\mathbf{E}(p)}{\sqrt{p^{\prime }}},$$

$$p^2=\frac{2e}{1+e}\in [0,1],$$

$$e=\frac{b}{R_c}\in [0,1],$$

$$\mathbf{E}(k)={\int }_0^{\pi /2}\sqrt{1-k^2{\sin }^{2}x}dx,$$

$${\Psi }_p$$

$$p=\frac{2\sqrt{k_c^{\prime }}}{1+k_c^{\prime }},$$

$$p^{\prime }=\frac{1-k_c^{\prime }}{1+k_c^{\prime }},$$

$${\Psi }_p=-16G\Sigma b\mathbf{K}(p^{\prime })\mathbf{E}(p).$$

$$M=4{\pi }^2\Sigma b{R}_c,$$

$${\Psi }_p=-\frac{GM}{R_c}\frac{4}{{\pi }^2}\mathbf{K}\left(\sqrt{\frac{1-e}{1+e}}\right)\mathbf{E}\left(\sqrt{\frac{2e}{1+e}}\right).$$

$${\Psi }_p=-8G\Sigma b\mathbf{K}(k_c)\left[2\mathbf{E}(k_c^{\prime })-{k_c}^2\mathbf{K}(k_c^{\prime })\right].$$

$$\text{u}(k)=\left(\begin{array}{c}\mathbf{K}(k)\mathbf{K}(k^{\prime })\\ \mathbf{K}(k)\mathbf{E}(k^{\prime })\\ \mathbf{E}(k)\mathbf{K}(k^{\prime })\\ \mathbf{E}(k)\mathbf{E}(k^{\prime })\end{array}\right),$$

$${\Psi }_p=-G\Sigma b{\text{f}}_{00}(k_c)\cdot \text{u}(k_c),$$

$${\text{f}}_{00}(k)=8\left(\begin{array}{c}-k^2\\ 2\\ 0\\ 0\end{array}\right).$$

$$-\frac{\partial \Psi }{\partial R}=G\lambda k\sqrt{\frac{a}{R^3}}\left[\left(1+\frac{a-R}{2a}\frac{k^2}{{k^{\prime }}^2}\right)\mathbf{E}(k)-\mathbf{K}(k)\right].$$

$$-{\frac{\partial \Psi }{\partial R}|}_{R_{\mathrm{p}}}$$

$$-\frac{G\Sigma }{\sqrt{R_{\mathrm{p}}}}\frac{{k_c}^3}{2{k_c^{\prime }}^2}\mathbf{E}(k_c)b\int \frac{d\theta }{\sqrt{a}},$$

$${\int }_0^{2\pi }\frac{d\theta }{\sqrt{a}}$$

Full text

††pagerange: Interior potential of a toroidal shell from pole values–LABEL:lastpage††pubyear: ???

1 Introduction

As elementary constituents of any model of slowly rotating star, thin homoeoids, i.e. infinitely thin ellipsoidal shells, are of major mathematical importance in Astrophysics (e.g. Binney & Tremaine, 1987). They are equipotential surfaces, as revealed already by Newton’s findings (Kellogg, 1929). In addition, the net force inside any solid ellipsoid (made of self-similar homoeoids), perpendicular to its main axis, happens to be linear with the cylindrical radius $R$, exactly as the centrifugal force does. As a consequence, any homogenous ellipsoid in rigid rotation around its main axis is a natural figure of equilibrium, as formulated by Maclaurin in the $18$th century (Chandrasekhar, 1969).

In this article, we investigate the toroidal analog of thin ellispoidal homoeoids. While tori, disks and rings are commonly observed over a vast range of scales in the Universe (around planets, thin and thick disks around forming and evolved stars, in galaxies), this topic has received almost no attention yet. Besides the physics of gravitation, various domains of science are concerned by toroidal structures : electrostatics (Andrews, 2006), plasma physics (Evangelidis, 1981; Throumoulopoulos & Tasso, 2000), nuclear physics (Wong, 1972), nano-structure physics and biology (Kuyucak et al., 1998).

The determination of gravitational attraction is a complex technical task, especially for toroids (see e.g. Cohl et al., 2000; Bannikova et al., 2011; Kondratyev & Trubitsyna, 2016; Majic, 2018). Here, we focus on the potential in the cavity of a homogeneous toroidal shell with a circular main radius and a circular core section. It is shown that not only the potential but all the partial derivatives can be expressed in terms of products of elliptic integrals of the first and second kinds at the pole of toroidal coordinates. As a consequence of this exceptional property, the interior potential can be expanded as a bivariate series of $R$ and $Z$, similar to the classical Taylor expansion. We derive the leading terms and show the performance of this approach, which is particularly good for shells with a small axis ratio $e$. We find that the interior potential is roughly linear with the cylindrical radius and weakly sensitive to the altitude $Z$, in contrast with ellipsoidal homoeoids.

The article is organized as follows. In Sect. 2, we remind the potential of a massive loop and show why the pole of toroidal coordinates plays a special role for the toroidal shell with a circular core section. The exact formula for the potential at the pole is established in Sect. 3. The corresponding acceleration is derived in Sect. 4. We then generate a two-term expansion of the potential around the pole and valid in the shell cavity. This is the aim of Sect. 5. We illustrate the method for a shell having an axis ratio of $e=0.1$. In this case, the interior potential is precise up to $4$ digits. We show in Sect. 6 how to proceed to next orders. Two additional terms are explicitly derived. The performance of the expansion versus $e\in[0,1]$ is discussed in Sect. 7. A short driver program is given in Appendix. In Sect. 8, we determine the potential at the surface of the shell and consider the case of shells with a thick edge. We also propose an empirical law for the exterior potential which is reliable for $e^{2}\ll 1$. The last section is devoted to a large discussion about the implementation of the formula, and in particular the gain with respect to a pure numerical treatment. We remind the importance of analytical solutions to understand the physics of systems hosting a massive (self-gravitating) torus. We end with a few remarks and perspectives.

2 Loop potential and the pole

2.1 The massive loop

We consider a massive circular loop of matter with radius $a$ and mass $2\pi a\lambda$, as shown in Fig. 1. The gravitational potential of this system is accessible at any point r of space by integration of the Green function $1/|\textbf{{r}}-\textbf{{r'}}|$ over the polar angle $\phi$. In cylindrical coordinates $(R,Z)$, the formula is (e.g. Kellogg, 1929; Durand, 1953)

[TABLE]

where

[TABLE]

is the complete elliptic integral of the first kind,

[TABLE]

is the modulus and $z$ is the altitude of the loop. This form assumes that the loop axis is confounded with the $Z$-axis. The potential is displayed versus $R$ and $Z$ in Fig. 2. It is logarithmically singular as $R\rightarrow a$ and $Z\rightarrow z$ (i.e. $k\rightarrow 1$) since ${\mathbf{K}}(k)\rightarrow\ln\frac{4}{k^{\prime}}$ (Byrd & Friedman, 1954; Fukushima, 2016), where $k^{\prime}=\sqrt{1-k^{2}}$ is the complementary modulus.

Any axially symmetrical system can be build by assembling co-axial loops, and the total potential is obtained by summing in Eq. (1) over the given distribution. We easily understand that the presence of the special function renders the accounting process (be discrete or continuous) somewhat cumbersome. Hence, a question arises : are there conditions and configurations under which the modulus $k$ would be a constant in Eq. (1)? If we rewrite Eq. (3) in the following form,

[TABLE]

we see that points $(a,z)$ basically belong to a circle, denoted ${\cal C}$. The centre C$(R_{c},Z_{c})$ of this circle is at

[TABLE]

and its radius is

[TABLE]

It is shown in Fig. 3. We conclude that $k$ is constant along ${\cal C}$, provided $R/R_{c}$ and $Z/R_{c}$ are fixed. We will see below which point of space is concerned. By varying $k$ in the allowed range, one gets a series of circles that are not concentric and do not have the same radius. For $k=0$, both $R_{c}$ and $b$ are infinite, and the circle is tangent to the vertical axis at the origin. Conversely, for $k=1$, the circle has null radius and $(R,Z)$ coincides with the centre C.

2.2 Link with toroidal coordinates

In axial symmetry, the toroidal coordinates $(\eta,\zeta)$ are linked to the cylindrical coordinates by

[TABLE]

where ${R_{\rm p}}>0$ corresponds to the pole (or focal ring). This point, which belongs to the radial axis, is singular in the sense that $R\rightarrow{R_{\rm p}}$ when $\eta\rightarrow\infty$. By eliminating $\zeta$ in Eq. (5), we see that lines of constant $\eta$ are non-intersecting circles surrounding the pole. In terms of cylindrical coordinates $(a,z)$, the equation of a given circle is

[TABLE]

where $({R_{\rm p}}\coth\eta,0)$ are the cylindrical coordinates of its centre and ${R_{\rm p}}/\sinh\eta$ is the radius. This is summarized in Fig. 4. A quick inspection of Eq. (4) shows that Eq. (8) is nothing but circle ${\cal C}$ (C,$b$) met above provided

[TABLE]

In summary, the circle ${\cal C}$ (C,$b$) is a particular line of constant $\eta$, which value is denoted $\eta_{c}$ in the following (the region inside the circle corresponds to $\eta>\eta_{c}$). We can therefore state that for any loop with parameters $(a,z)\in{\cal C}$, the modulus $k$ is constant at the pole of toroidal coordinates, and only at that point.

3 Potential at the pole

3.1 A new formula

We now consider a continuous concatenation of co-axial massive loops, each with radius $a$ and altitude $z$, such that $(a,z)\in{\cal C}$. This forms a toroidal shell with circular section and core radius $b$ (see again Figs. 3 and 4). In cylindrical coordinates, the equation of ${\cal C}$ is basically of the form

[TABLE]

where $\theta\in[0,2\pi]$. It follows from Eq. (1) that the total potential of the shell is given by the integral

[TABLE]

where $\Sigma$ is the surface density111The integral bounds can be modified if, instead of a full circular section, we consider an arc.. If all loops have the same mass per unit length, then $\Sigma$ is a constant.

An example of direct numerical integration of Eq. (11) is given in Fig. 5 for a homogeneous toroidal shell with axis ratio $b/R_{c}=0.1$, which corresponds to $\eta_{c}\sim 2.993$. As the quadrature scheme, we use the trapezoidal rule with $N_{\theta}=1025$ equally spaced nodes in $\theta$ (this setup is not critical). The imprint of the shell surface is clearly visible. In the exterior domain (i.e. $\eta<\eta_{c}$), the potential decreases in absolute when one moves away from the shell (i.e. $\eta\rightarrow 0$). Inside the whole cavity ($\eta>\eta_{c}$), the potential gradually increases with $R$. From the above statement (see Sect. 2), we see that $k$ (and then any function of $k$) can be taken out of the integral if $(R,Z)=({R_{\rm p}},0)$. So, Eq. (11) becomes

[TABLE]

where $a$ is given by Eq. (10a) and $k$ has been set to ${k_{c}}$. We see the remaining integral is analytical for a wide variety of angle-dependent surface densities $\Sigma(\theta)$. The result is straightforward for the homogeneous shell. We actually have

[TABLE]

where

[TABLE]

is the ratio of the core diameter to the outer radius,

[TABLE]

is the axis ratio, $p^{\prime}=\sqrt{1-p^{2}}$ is the complementary modulus, and

[TABLE]

is the complete elliptic integral of the second kind. So, Eq. (12) becomes

[TABLE]

We are not aware that this expression has already been reported in the literature. Clearly, ${\mathbf{K}}$ corresponds to the main curvature (radius $R_{c})$ while ${\mathbf{E}}$ is associated with the shell circular section (radius $b$). Note that, when $p\rightarrow 0$, $e\rightarrow 0$ while ${R_{\rm p}}\rightarrow R_{c}$. In this case, the shell tends to the loop, Eq. (17) matches Eq. (1) and its mass converges to $2\pi\lambda a$. At the opposite, for $p\rightarrow 1$, the shell has infinite radius, mass and potential.

3.2 Equivalent forms

From Eqs. (5a), (6), (14) and (15), we see that $p$ is linked to ${k_{c}}$ through the following relations

[TABLE]

and

[TABLE]

meaning that Eq. (17) can take equivalent forms by considering modulus transformations (see the Appendix A). From Eq. (74), we can eliminate ${k_{c}}$, and we find

[TABLE]

Since the mass of the toroidal shell is

[TABLE]

we can easily rewrite the potential as a function of $GM/R_{c}$ (or even $GM/{R_{\rm p}}$). This is for instance

[TABLE]

If we can eliminate ${\mathbf{E}}(p)$ in Eq. (17) by using Eq. (75), we find

[TABLE]

In order to anticipate a little bit, we introduce the four-vector

[TABLE]

whose components are plotted versus $k$ in Fig. 6. We see that Eq. (23) reads

[TABLE]

where

[TABLE]

Note that we can expand the complete elliptic integrals as a function of $k$ in extreme regimes where $k\rightarrow 0$ and $k\rightarrow 1$. Then $\Psi_{p}$ is only function of ${k_{c}}$ and $\ln{k^{\prime}_{c}}$.

4 Acceleration at the pole

We can make the same analysis for the gravitational acceleration. Since the pole belongs to the $R$-axis, the vertical acceleration is zero by symmetry. The formula for the radial acceleration caused by the loop can be determined from Eq. (1). The classical form is (Durand, 1953; Huré, 2005)

[TABLE]

To get the acceleration due to the toroidal shell, we have to multiply this expression by $\Sigma bd\theta/\lambda$, and to integrate over $\theta\in[0,2\pi]$. At the pole, any function that depends only on $k$ can be carried out of the integral, which requires to set $k={k_{c}}$ (see Sect. 2). We then find, again assuming $\Sigma$ constant :

[TABLE]

The first integral has already been met; see Eq. (13). The second one can be easily calculated. We find

[TABLE]

and so Eq. (28) becomes

[TABLE]

By using Eqs. (75) and (76), we can eliminate both ${\mathbf{K}}(p)$ and ${\mathbf{E}}(p)$. The final expression then depends only on ${\mathbf{E}}$ and ${\mathbf{K}}$ evaluated at ${k_{c}}$ and ${k^{\prime}_{c}}$, and it can be put in a form similar as Eq. (25), namely

[TABLE]

where

[TABLE]

As done for the potential, we can write this term as a function of the mass of the shell by using Eq. (21).

5 Route to the interior potential

There are three striking properties of $\Psi$ inside the cavity of the toroidal shell where $\eta\geq\eta_{c}$. First, the impact of variable $Z$ is very weak (almost not visible by eyes), especially if the core radius is small. Second, $\Psi$ is an increasing function of $R$ as already mentioned. The gravitational acceleration is oriented towards the origin of coordinates for $R\geq R_{c}-b$. Third, the variation of $\Psi$ with $R$ is very close to linear. Proofs are given in Fig. 7, which shows $\Psi$ at the equatorial plane (i.e. for $Z=0$), at $R={R_{\rm p}}$ and at $R=R_{c}$ for the shell considered in Fig. 5. The contrast with the ellipsoidal shell is therefore evident: the interior potential is not a constant but has a small positive gradient with $R$ due to the curvature around the $Z$-axis (e.g. Kellogg, 1929; Binney & Tremaine, 1987). Actually, for an observer standing at the inner edge of the shell and looking towards the origin, matter is present at relatively short seperations behind (the outer edge) and in front (opposite inner and outer edges). In contrast, at the outer edge, there is no matter behind and separations are larger. The potential well is therefore deeper at the inner edge, which is rather intuitive.

It is tempting to elaborate some kind of a fit. We could produce sets of data by varying $R/R_{c}$, $Z/R_{c}$ and $e$, but it seems more powerful to consider a bivariate series, resembling the Taylor series, i.e.

[TABLE]

since the first terms have already been calculated. Any expansion is, however, necessarily limited to the cavity and cannot be valid in the whole physical space. Actually, $\Psi$ is continuous but not differenciable at the surface of the shell. Besides, the potential must satifiy the Laplace equation in the cavity, namely

[TABLE]

which is a priori not automatic with the above form.

5.1 Zero-order approximation

We start with the crudest approximation, i.e.

[TABLE]

in which case Eq. (34) is trivially fulfilled. By constant, we mean that the potential is uniform inside the cavity. Obviously, this constant value is function of the shell parameters $\Sigma$, $b$ and ${k_{c}}$ (or $e$). We then have $C_{0}=\Psi_{p}$. We have checked this approximation by comparing $\Psi_{p}$ to values (hereafter “reference values”) obtained by the direct integration of Eq. (11); see Sect. 3. The conditions are the same as for Fig. 5. Data have been produced in the cavity only, which is easily done in toroidal coordinates by varying $\eta$ in the range $[\eta_{c},\infty]$ and $\xi\in[0,2\pi]$. The logarithmic deviation is shown in Fig. 8 as a function of $R/R_{c}$ and $Z/R_{c}$. As expected, the potential given by Eq. (35) is underestimated close to the inner edge, and overestimated close to the outer edge. We see, however, that the agreement is globally correct and even better than expected. The mean relative deviation is $\sim 10^{-2.23}$, i.e. less than one percent. We see that the precision is nominal not only around the pole but also along two directions $|\theta|\sim\pi/2$. If the core radius is decreased, the acceleration also decreases, the interior potential is flatter, and the approximation becomes even better (see Sect. 7).

5.2 First-order approximation (the curvature effect included)

In the absence of any variation in $Z$, the radial contribution must vanish in Eq. (34), which means

[TABLE]

where $C_{0}$ and $C_{1}$ are to be determined. In particular, at the pole, we have

[TABLE]

Since $\Psi_{p}$ and $-\left.\partial_{R}\Psi\right|_{{R_{\rm p}}}$ are known (see Sect. 3 and 4), this set of equations is easily solved for $C_{0}$ anf $C_{1}$, making Eq. (36) fully operational. We have compared this new approximation to reference values under the same conditions as above. Figure 9 displays the decimal logarithm of the relative error. We notice that the precision is improved by more than two orders of magnitude, with a mean deviation of $\sim 10^{-4.70}$. It is nominal around the pole, as well as in four directions $|\theta|\sim\pi/4$ and $3\pi/4$. Note that the quasi-linear behavior of $\Psi$ in the cavity is especially marked for shells with a small axis ratio. This is explained by expanding $\ln R$ around ${R_{\rm p}}$, which gives

[TABLE]

6 Next orders

6.1 Second-order approximation

The variation of $\Psi$ inside the cavity is much more complex, with a sensitivity to variable $Z$ (in fact, to $Z^{2}$ due to equatorial symmetry). If we allow for quadratic contributions, the vacuum solution is of the form

[TABLE]

where $C_{0}$, $C_{1}$ and $C_{2}$ are three quantities to be determined2. At the pole, we have

[TABLE]

and we see that a third equation is needed to fix $C_{2}$. This is done by calculating the vertical gradient of the vertical acceleration, namely

[TABLE]

Alternatively, we could also consider $\partial^{2}_{R^{2}}\Psi$ (since $\Psi$ is harmonic, derivatives with respect to $Z$ and $R$ are linked together; see below). We then reiterate the procedure described in Sects. 3 and 4. The vertical acceleration due to the circular loop is given by (Durand, 1953; Huré, 2005)

[TABLE]

and so the second derivative writes

[TABLE]

From Eq. (3), we find

[TABLE]

Besides, from Eq. (77) of Appendix B, we get

[TABLE]

where we have defined

[TABLE]

for convenience. To get $\partial^{2}_{Z^{2}}\Psi$ for the toroidal shell, we just have to multiply Eq. (43) by $\Sigma bd\theta/\lambda$ and to integrate over $\theta$. In this process, $a$ and $z$ are still given by Eq. (10). The complexity of the calculus is apparent since, at the pole, $k$ is a constant (see Sect. 2). Only terms with $a$ and $z$ are therefore retained in the integral, and we set $k={k_{c}}$. We then have

[TABLE]

where the term $\sin^{2}\theta$ comes from $(z-Z)^{2}$, evaluated at the pole. The first integral has already been employed above, and the second one is analytical, namely

[TABLE]

After some algebra and tedious calculus, we find

[TABLE]

which leads to $C_{2}$ from Eq. (41), and subsequently to $C_{0}$ and $C_{1}$ from Eq. (40). As done before, we can express ${\mathbf{K}}(p)$ and ${\mathbf{E}}(p)$ as functions of ${\mathbf{K}}({k^{\prime}_{c}})$ and ${\mathbf{E}}({k^{\prime}_{c}})$, and put the result in the form of a scalar product. The final formula is

[TABLE]

where

[TABLE]

while u is given by the Eq. (24). Again, this derivative can be expressed as a function of the mass of the shell and $R_{c}$ or ${R_{\rm p}}$. From Eqs. (31), (34) and (50), we deduce

[TABLE]

where

[TABLE]

Figure 10 displays the decimal logarithm of the relative error between Eq. (39) and reference values. The numerical setup is unchanged. We see that the deviation is now reduced by almost two orders of magnitude, with a mean value of $\sim 10^{-6.47}$. The precision is nominal in the vicinity of the pole, as well as in now six directions $|\theta|\sim\pi/6,3\pi/6$ and $4\pi/6$.

6.2 Third-order approximation

We can in principle continue the process up to the desired order, say $N+1\geq 2$. The expansion is then supplemented with a new set of terms

[TABLE]

where the coefficients $\alpha_{i,N-i}$ are found from Eq. (34), and the leading quantity $C_{2N}$ is determined from a $N$th-order partial derivative of the potential evaluated at the pole, which can be put into the form

[TABLE]

where $N+1=m+n$, where $m\geq 0$, $n\geq 0$, and $\textbf{{f}}_{nm}(k)$ is a four-vector whose components are rational functions of $k$. This is the most tricky part of the method. Here are the results for $N=2$. Including quadratic terms, the vacuum potential has the form

[TABLE]

where there are now four unknowns involved. At the pole, we have

[TABLE]

which can be solved if we can add a new equation. We see that $C_{4}$ is accessible from a cross-derivative, namely

[TABLE]

which can be calculated from Eq. (43). We actually have

[TABLE]

From Eq. (3), we get

[TABLE]

and we find (see also Appendix B)

[TABLE]

with the help of Eqs. (77) and (78). At the pole, we set $R={R_{\rm p}}$, $Z=0$, and $k={k_{c}}$. We see that two integrals involving $1/\sqrt{a^{3}}$ and $z^{2}/\sqrt{a^{5}}$ appear. Still assuming $\Sigma=$const., we have

[TABLE]

and

[TABLE]

respectiveley. So, the right-hand-side of Eq. (58) is

[TABLE]

where we have applied, as above, modulus transformations, and

[TABLE]

The $3$rd-order approximation is then ready to be used since $C_{0}$, $C_{1}$, $C_{2}$ and $C_{4}$ are found from Eqs. (57) and (58). The reader can find in the Appendix C a short driver program written in Fortran 90 to compute $\Psi$ from Eq.(56). The comparison of this new approximation with direct numerical integration is shown in Fig. 11. The conditions are the same as above. The precision is increased compared to $2$nd-order, by two more digits typically. In the cavity, the mean deviation is now $\sim 10^{-8.00}$. The pattern is similar : the smallest deviations are found around the pole and along eight directions defined by $|\theta|\sim\pi/8,3\pi/8,5\pi/8$ and $7\pi/8$.

7 Effect of the shell axis-ratio

Examples given above concern a toroidal shell with $e=0.1$ (i.e. ${k_{c}}\approx 0.9987$). It is interesting to see how all these expansions behave when varying the axis ratio. We repeated the computations for the full range, i.e. $e\in\,]0,1[$. For each parameter, potential values inside the cavity have been generated by direct computation of Eq. (11), and for the four approximations, i.e. from Eq. (36) for order [math], Eq. (36) for order $1$, Eq. (39) for order $2$ and Eq. (56) for order $3$. The log. of the relative differences, i.e. $\log|\Delta\Psi/\Psi|$, have been averaged. The results are displayed in Fig. 12. Unsurprisingly, none of the approximations is really reliable for the largest values of $e$. We remind that, when $b\rightarrow R_{c}$, the inner edge of the shell is close to the origin, while its outer edge tends to infinity, and curvature effects are important at short radii. In fact, the interior potential shows a complex, saddle-type structure. In contrast, as soon as $e\lesssim 0.3$, the four approximations are very efficient. Accuracy increases as the shell section gradually decreases. The sensitivity (i.e. slope) depends on the order too. We find

[TABLE]

8 Byproducts

8.1 Potential at the surface of the shell

Potential values on the shell (i.e. along the circle ${\cal C}$) can be estimated from any approximation derived above since the potential is continuous inside the cavity. These values can be used for instance as approximate (but accuracte) Dirichlet boundary conditions for solving the Poisson equation. There are different equivalent options, depending on the variable selected to describe the section of the shell, mainly $R(Z)$, $Z(R)$, $\theta$ or $\eta$, which can even be converted into $k$ from Eq. (9c). If we use cylindrical variables, then $Z(R)^{2}=b^{2}-(R-R_{c})^{2}$ and so we get

[TABLE]

where $R\in[R_{c}-b,R_{c}+b]$. Figure 13 shows the results obtained at orders $1$ to $3$ (order 0 is too crude to be retained). We see that the linear approximation given by Eq. (38) is quite acceptable. By varying $e$, we notice that the potential becomes flatter and flatter along ${\cal C}$, and in the whole cavity as well. This is expected because, as $e\rightarrow 0$, curvature effects become less and less pronounced. The shell becomes, locally, similar to a hollow filament, for which the interior potential is a constant by virtue of the Gauss theorem.

8.2 Shells with a thick edge, i.e. tubes

If we concatenate coaxial toroidal shells having the same main radius $R_{c}$ but a different core radius $b\in[{b_{\rm in}},{b_{\rm out}}]$, then one gets a toroidal shell with thick edge, i.e. a tube. This system is depicted in Fig. 14. Each shell has its own pole. It follows that there is no pole for the tube (otherwise, this would imply that $b/R_{c}$ is constant for each constitutive shell). The total potential is obtained from Eq. (11) by integration over $b$, i.e.

[TABLE]

where $\rho$ is the local mass density. Figure 15 shows the potential (in units of $GM/R_{c}$) computed by direct integration222We have used the trapezoidal rule as the quadrature scheme both in $\theta$ as before, and in $b$ with $N_{b}=257$ equally spaced nodes. for a homogeneous tube with parameters ${b_{\rm in}}/R_{c}=0.05$ and ${b_{\rm out}}/R_{c}=0.1$. Figure 16 is for the equatorial plane only. The mass of the tube is

[TABLE]

We observe the rounded shape of the potential inside the tube. There is no more jump in the derivative of the potential when crossing the system boundaries, as expected for any three-dimensional distribution. At short radii, the potential decreases with $R$. This continues when passing through the outer shell at $R=R_{c}-{b_{\rm out}}$. The minimum (i.e. the place where the acceleration is zero) is reached inside matter, just before crossing the inner shell at $R=R_{c}-{b_{\rm in}}$. Next, the potential increases with $R$ in the cavity, and again inside matter, and this is so up to infinity. Inside the internal cavity, we notice the quasi-linear behavior of $\Psi$ with $R$, which is a consequence of what has been observed for the shell. This is basically the expression of the superposition principle. The central cavity is actually common to all shells the tube is made of. This means that all the approximations presented in Sects. 5 and 6 are valid in this common cavity, i.e. for $R\in[R_{c}-{b_{\rm in}},R_{c}+{b_{\rm in}}]$

How does it work? For instance, if we consider order [math], we have from Eq. (35)

[TABLE]

It is important to realise that $C_{0}=\Psi_{p}$ is a function of $b$, not only through the prefactor $-G\Sigma b$, but also through ${k_{c}}$. We see that $\rho$ can also depend on $b$ if the core is stratified, and there is no difficulty there. By using Eq. (20) and changing $b=eR_{c}$ for $p$ from Eq. (14), we obtain

[TABLE]

where ${p_{\rm in}}$ and ${p_{\rm out}}$ are found from Eq. (14), and $\rho=$const. is assumed. This integral, which can also be written in terms of ${k_{c}}$ or ${k^{\prime}_{c}}$, has been found nowhere in the literature (however, see Prudnikov et al., 1988; Wan, 2012). It must therefore be estimated numerically. If we consider order $1$ for approximating $\Psi$, then we have from Eq. (36)

[TABLE]

There is no special difficulty to proceeding to next orders, which has been checked. The results for the tube fully agree with what is reported in Sect. 7 for the shell. Note that if we set ${b_{\rm in}}\rightarrow 0$, then the cavity gets smaller and smaller. The domain of validity of the above formula gets shorter and shorter. For ${b_{\rm in}}=0$, the cavity even disappears and the tube becomes a solid torus. Only the potential at $R_{c}$ becomes accessible.

8.3 An empirical law for the exterior potential

The potential outside the cavity of the toroidal shell is not accessible by the formula derived above but it is clearly needed. This would ensure a full coverage of space, which is necessary for dynamical studies (Šubr et al., 2004; Bannikova et al., 2012). We show in Fig. 16 the equatorial potential in units of $GM/R_{c}$ for the four systems considered in this article, i.e. the loop, the shell, the tube and the solid torus. We used the same radius $R_{c}$ and the same core radius ${b_{\rm out}}/R_{c}=0.1$. Since the mass densities $\lambda$, $\Sigma$ and $\rho$ have been set to unity, these systems have not the same mass. We notice that, while interior potentials are very different, exterior potentials superimpose with each other remarkably. This interesting fact has been outlined in Bannikova et al. (2011). However, the graph above can be misleading. The concordance of exterior potentials is true for small axis-ratios only and fails for $e\gtrsim 0.3$ typically. So, if the condition $e^{2}\ll 1$ is satisfyied, the potential outside the cavity of the toroidal shell is close to Eq. (1) when appropriately scaled. When coupled to one of the approximations given in the previous sections, one gets a full coverage of the physical space. In summary, by using the $1$st-order approximation for the interor solution (higher-orders can be considered instead), we have

[TABLE]

where $M$ is given by Eq. (21) and $k_{0}$ comes from Eq. (3) with $a=R_{c}$ and $z=0$. A simple program is given in the Appendix D for Eq. (73).

Besides, the knowledge of the exterior potential would enable us to determine the potential of the solid torus, which is composed of the interior potential of a tube and the exterior potential of a smaller solid torus, i.e. $\int_{0}^{b_{\rm out}}\dots db^{\prime}=\int_{0}^{b}\dots db^{\prime}+\int_{b}^{b_{\rm out}}\dots db^{\prime}$. From this point of view, an interesting and probably more sraightforward option would be to use the contour integral as the starting point (Ansorg et al., 2003; Huré et al., 2014). This question is open.

9 Discussion

9.1 About the implementation

In general, an analytical approach is more powerful than a fully numerical treatment, mainly because one captures the sensitivity of all the parameters involved on the resulting quantity. One has, however, to make sure that the implementation is not too tricky, and that the advantages in terms of precision and computing time are real. This is especially true when dealing with special functions and series as it is the case here. For instance, expanded Green functions are very often used in potential theory but their efficiency is not always satisfactory (Clement, 1974; Huré et al., 2014)

We first notice that the differences between the four approximations are indeed minor. For a given value of the shell axis-ratio $e$, only $4$ evaluations of complete elliptic integrals and a few rational functions of $k^{2}$ are needed to form $\vec{u}$ and $\vec{f}_{ij}$ respectively. These quantities, and subsequently the four constants $C_{0}$, $C_{1}$, $C_{2}$ and $C_{4}$, can be computed once for all. The potential at a single point $(R,Z)$ of the cavity is then obtained by computing only a few polynomials in $R$ and $Z$ in addition. The comparison with the direct numerical integration is then straightforward. By using the trapezoidal rule with $N_{\theta}+1$ equally spaced nodes, the number of evaluations of ${\mathbf{K}}(k)$ amounts to $N_{\theta}$ (values at $\theta=0$ and $\theta=2\pi$ are identical). This quadrature scheme being second-order accurate in the grid spacing $2\pi/N_{\theta}$, the absolute error is $E\sim 2\pi^{3}/3N_{\theta}^{2}$. So, $4$ different evaluations of the ${\mathbf{K}}(k)$ means $N_{\theta}=4$. The corresponding error is therefore of the order of unity, i.e. much larger than one percent which is the typical value obtained at order [math] (see Fig. 12 for $e\approx 0.01$). In terms of equivalent error, we need $N_{\theta}\sim\sqrt{2\pi^{3}/3E}$ to reach a given error level $E$. This gives $N_{\theta}\sim 45$ for $E=0.01$ and $N_{\theta}\sim 450$ for $E=0.0001$ (for $e\approx 0.001$). Since the computing time is mainly governed by the determination of the elliptic integrals, we see that the direct numerical integration is, by orders of magnitudes, always inferior to the analytical approach, which is what is observed in practice (see Appendix C for a possible implementation).

9.2 From shells to fluid tori

The results presented in this paper belong to the domain of classical theoretical physics. The derivation of the (vacuum) gravitational potential inside the geometrically thin massive shell, as discussed in the present paper, is a first step towards a complete solution for the potential within the matter distribution of the fluid torus. This is a longstanding and recurrent challenge in Astrophysics (Dyson, 1893b, a; Vogt & Letelier, 2009; Fukushima, 2010; Bannikova et al., 2011), and in Electrostatics as well (Belevitch & Boersma, 1983; Hernandes & Assis, 2004; Scharstein & Wilson, 2005; Majic, 2018). The knowledge of the gravity field of toroidal systems is important to derive the equilibrium structure, the shape of the surface, and the their time evolution. It is also fundamental to investigate the dynamics of particles travelling around (e.g. Šubr et al., 2004; Šubr & Karas, 2005). From this point of view, Eq. (73a), although empirical, enables such a study.

Any analytical solution readily available can be useful as a test bed for numerical approaches to more astrophysically realistic systems. Actually, we can envisage two types of objects where massive toroidal configurations are relevant. Gravitating very massive tori have been considered as a transient stage during the merger process and the associated tidal disruption event of a neutron star in the close binary system. The remnants can spread and form a toroidal structure around the tidal radius, which eventually becomes partly accreted onto the central body (presumably a black hole) and partly expelled in the form of jet or a massive outflow (e.g. Abramowicz et al., 1998; Masuda et al., 1998; Lee, 2001). Besides, gamma rays and neutrinos have been proposed to originate from the short-lived dense torus (Woosley, 1993; Janiuk, 2014). At much lower scales (a factor $10^{3}$ to $10^{6}$ typically), the formation of the Moon has probably involved a transient ring-like structure of gas, dust and peebles after the early Earth has been impacted by a large body Lock & Stewart (2017). On length-scales larger by a factor $10^{6}$ to $10^{9}$, massive tori are believed to orbit supermassive black holes in active galactic nuclei (e.g. Goodman, 2003; Karas et al., 2004). Even though the total mass of the torus is thought to be less than the central mass concentrated in SMBH in most observed systems, it has been clearly demonstrated that self-gravity must plays a significant role (Collin & Zahn, 2008; Huré, 2000).

9.3 Concluding remarks and perspectives

The article present a novel contribution to the precise construction of the gravitational potential in toroidal systems. Seen from different angles, both the analytical approaches as well as the numerical computations are rather heavy and the evaluations are always cumbersome. It is thus very useful to develop alternative approaches and simplified models based either on course grids in the numerics or fitting formulas in the analytics; here we attempt to contribute to the latter. We have shown that the potential and its successive derivatives happen to be analytical at the pole of toroidal coordinates for a toroidal shell with circular radius and circular section. On these grounds, we have proposed a new kind of expansion for the potential in the cavity of a shell. We have determined the leading terms up to order $3$. The corresponding approximations are of great precision when compared to the direct numerical integration. This is for instance $8$-digit at order $3$ for a shell with an axis ratio of $e=0.1$. Another major result concerns Newton’s theorem for the ellipdoidal homoeoid which cannot be transposed to the toroidal case. The existence of the curvature (around the $z$-axis) makes the potential well deeper at the inner edge of the shell than at the outer edge. Matter is more concentrated at short radii.

This work adresses new questions and requires further developments. While the level of accuracy reached is already very high with te $3$rd-order formula, there is a priori no limit in including more terms in the expansion of the interior potential, but an increasing effort is required to derive new partial derivatives evaluated at the pole. It would be interesting to understand the spike pattern that the error maps exhibit. This is probably the emergence of modes in $\xi$ related to the solution of the Laplace equation in toroidal coordinates with variable seperation (see Majic, 2018). This point remains to be clarified. Finally, this approach can be transposed to the case of a toroidal current density relevant in astrophysical and laboratory plasmas (Slaný et al., 2013). Actually, in the case of axial symmetry, the vector potential A for the current loop is quite similar to Eq. (1) (see Jackson, 1998; Cohl & Tohline, 1999). This would give access to the poloidal component of the magnetic field.

Acknowledgments

J.M. Huré is grateful to the Astronomical Institute in Prague where this work has been initiated during winter 2016 and fundings for the visit. A. Trova acknowledges support from the Research Training Group 1620 “Models of Gravity” funded by the German Science Foundation DFG. V. Karas thanks the Czech. Science Fundation (ref. 19-011373). We thank V. Bruneau, N. Popoff at the IMB and B. Boutin-Basillais for fruitful discussions. The anonymous referee is acknowledged for the constructive feedback.

Appendix A Modulus transformations

From Gradshteyn & Ryzhik (2007), we have

[TABLE]

[TABLE]

and

[TABLE]

Appendix B Derivatives of ${\mathbf{K}}$ and ${\mathbf{E}}$

From Gradshteyn & Ryzhik (2007), we have

[TABLE]

and

[TABLE]

Appendix C F90 program for the interior potential

Program F90drivercode ! "Interior potential of a toroidal shell from pole values" ! Hure, Trova, Karas & Lesca (2019), MNRAS ! gfortran F90drivercode.f90; ./a.out Implicit None Integer,Parameter::SP=Kind(1.00E+00),DP=Kind(1.00D+00),QP=Kind(1.00Q+00) Integer,Parameter::AP=DP Real(Kind=AP),Parameter::PI=ATAN(1._AP)4 Real(KIND=AP),Parameter::EPSMACH=Epsilon(0._AP) Real(KIND=AP)::B,RC,MASS,RP ! core radius, main radius and mass of the shell, and position of the pole Real(KIND=AP)::KMODC,KMODC2,KPRIMC,KPRIMC2 ! various modulus Real(KIND=AP)::R,Z,PSI ! cylindrical coordinates and potential value where it is estimated Real(KIND=AP)::PSIP,DPSIDRP,D2PSIDZ2P,D3PSIDZ2DRP ! partial derivatives of the potential at the pole Real(KIND=AP)::C0,C1,C2,C4 ! coefficients of the expansion, Eq.(56) Real(KIND=AP),Dimension(1:4)::UVECT ! components of the u-vector, Eq.(24) Real(KIND=AP),Dimension(1:4)::FVECT00,FVECT10,FVECT02,FVECT12 Real(KIND=AP)::ELLIPTICK,ELLIPTICKPRIM,ELLIPTICE,ELLIPTICEPRIM ! ? input parameters (properties of the shell) B=0.1_AP RC=1._AP MASS=BRCPI**24 print*,"Mass of the shell",MASS ! statements KPRIMC=(Sqrt(RC+B)-Sqrt(RC-B))/(Sqrt(RC+B)+Sqrt(RC-B)) KMODC2=Sqrt((RC-B)(RC+B))/(Sqrt(RC+B)+Sqrt(RC-B))**24 KPRIMC2=1._AP-KMODC2 KMODC=Sqrt(KMODC2) RP=RCKMODC2/(2._AP-KMODC2) print,"Radius of the pole",RP ! components of the u-vector, Eq.(24) - values of K(k=KMODC), K(k’), E(k) and E(k’) to be set ! UVECT(1:2)=ELLIPTICK*(/ELLIPTICKPRIM,ELLIPTICEPRIM/) UVECT(3:4)=ELLIPTICE*(/ELLIPTICKPRIM,ELLIPTICEPRIM/) ! components of the f_ij-vectors, Eq.(26), (32), (51) and (65) FVECT00(1:4)=(/-KMODC2,2._AP,0._AP,0._AP/)8 FVECT10(1:4)=(/KMODC2KPRIMC2,-KPRIMC22,-KMODC2,1._AP+KPRIMC2/)/KPRIMC24 FVECT02(1:4)=(/-KPRIMC2,KPRIMC2,2._AP+KPRIMC2,-(1._AP+KPRIMC2)2/)(KMODC2/KPRIMC2)2 FVECT12(1:4)=(/KPRIMC2(-4._AP+KPRIMC23-KPRIMC222),& &KPRIMC2(4._AP-KPRIMC2+KPRIMC224),8._AP-KPRIMC27-KPRIMC22*2+KPRIMC234,& &(1._AP+KPRIMC2)(-8._AP+KPRIMC211-KPRIMC2**28)/)KMODC22/KPRIMC232/3 ! ? values of R and Z where the potential is requested (must be inside the cavity!) R=RP Z=0. If ((R-RC)2+Z2-B2<0._AP) Then PSIP=-DOT_PRODUCT(FVECT00,UVECT)*B DPSIDRP=-DOT_PRODUCT(FVECT10,UVECT)*B/RP D2PSIDZ2P=-DOT_PRODUCT(FVECT02,UVECT)*B/RP2 D3PSIDZ2DRP=-DOT_PRODUCT(FVECT12,UVECT)B/RP**3 ! order 3, Eq.(56); set C4=0 for order 2, and set C2=0 for order 1, and set C1=0 for order 0 C4=-D3PSIDZ2DRP/RP/32 C2=(-D2PSIDZ2P-C4RP216)/4 C1=(DPSIDRP-C2RP2-C4RP34)RP C0=PSIP-C1LOG(RP)-C2RP2-C4*RP4 PSI=C0+C1LOG(R)+C2(R2-Z22)+C4(R4-R2Z**28+Z**4*8/3) Print ,"Potential value (3rd-order)",PSI,PSI/MASSRC Endif

Appendix D F90 program for the interior and exterior (empirical) potential

Program F90drivercode2 ! "Interior potential of a toroidal shell from pole values" ! Hure, Trova, Karas, & Lesca (2019), MNRAS ! gfortran F90drivercode2.f90; ./a.out Implicit None Integer,Parameter::SP=Kind(1.00E+00),DP=Kind(1.00D+00),QP=Kind(1.00Q+00) Integer,Parameter::AP=DP Real(Kind=AP),Parameter::PI=ATAN(1._AP)4 Real(KIND=AP),Parameter::EPSMACH=Epsilon(0._AP) Real(KIND=AP)::B,RC,MASS,RP ! core radius, main radius and mass of the shell, and position of the pole Real(KIND=AP)::KMODC,KMODC2,KPRIMC,KPRIMC2,K0 ! various modulus Real(KIND=AP)::R,Z,PSI ! cylindrical coordinates and potential value where it is estimated Real(KIND=AP)::PSIP,DPSIDRP,D2PSIDZ2P,D3PSIDZ2DRP ! partial derivatives of the potential at the pole Real(KIND=AP),Dimension(1:4)::UVECT ! components of the u-vector, Eq.(24) Real(KIND=AP),Dimension(1:4)::FVECT00,FVECT10 Real(KIND=AP)::ELLIPTICK,ELLIPTICKPRIM,ELLIPTICE,ELLIPTICEPRIM INTEGER::I,J ! ? input parameters (properties of the shell) B=0.1_AP RC=1._AP MASS=BRCPI**24 print*,"Mass of the shell",MASS ! statments KPRIMC=(Sqrt(RC+B)-Sqrt(RC-B))/(Sqrt(RC+B)+Sqrt(RC-B)) KMODC2=Sqrt((RC-B)(RC+B))/(Sqrt(RC+B)+Sqrt(RC-B))**24 KPRIMC2=1._AP-KMODC2 KMODC=Sqrt(KMODC2) RP=RCKMODC2/(2._AP-KMODC2) print,"Radius of the pole",RP ! components of the u-vector, Eq.(24) - values of K(k=KMODC), K(k’), E(k) and E(k’) to be set ! UVECT(1:2)=ELLIPTICK*(/ELLIPTICKPRIM,ELLIPTICEPRIM/) UVECT(3:4)=ELLIPTICE*(/ELLIPTICKPRIM,ELLIPTICEPRIM/) ! components of the f_ij-vectors, Eq.(26), (32), (51) and (65) FVECT00(1:4)=(/-KMODC2,2._AP,0._AP,0._AP/)8 FVECT10(1:4)=(/KMODC2KPRIMC2,-KPRIMC22,-KMODC2,1._AP+KPRIMC2/)/KPRIMC24 ! 51x51 grid on [0,2]x[-1,1], for Eq.(73) Do I=0,50 R=2._API/50 Do J=0,50 Z=2._APJ/50-1._AP If ((R-RC)2+Z2-B2<0._AP) Then ! Eq.(73b) PSI=-MASS/RC*Dot_PRODUCT(FVECT00+LOG(R/RP)*FVECT10,UVECT)/PI2/4 Else ! Eq.(73a) ! value of K(k0) to be set ! K0=Sqrt(R*RC/((R+RC)2+Z2))2 PSI=-MASS/Sqrt((R+RC)2+Z2)ELLIPTICK/PI2 Endif Print,R,Z,PSI Enddo Enddo