2DStructure/ToroidalCoordinates
Using Toroidal Coordinates to Determine the Gravitational Potential
Preface
It appears as though we have found an answer to this question posed by Trova, et al. (2012). The detailed derivations and associated scratchwork that support the summary discussion of this chapter can be found under the Appendix/Ramblings category of this H_Book and in our accompanying discussion of DysonWong toroids.
As has been known for approximately a century, at any meridionalplane coordinate location, , the gravitational potential due to an axisymmetric, infinitesimally thin ring (TR) of radius, , and total mass, , is given by the "key expression,"
and is the complete elliptic integral of the first kind. Making the coordinatename substitutions, , to match much of this chapter's variable notation, we have alternatively,
A number of research groups (e.g., Cohl & Tohline (1999), Bannikova et al. (2011), Trova, et al. (2012), Fukushima (2016)) have recognized that the gravitational potential due to any axisymmetric configuration of finite extent can be obtained by integrating over the meridionalplane surface area that is occupied by the configuration, while weighting the volume density, , according to the prescription,
where,
We, like Trova, et al. (2012), have wondered whether there is a possibility of converting this double integral into a single (line) integral. This is a particularly challenging task when this expression for the gravitational potential is couched in terms of cylindrical coordinates because the modulus of the elliptic integral is explicitly a function of both and . We have recently realized that if a switch is made from cylindrical coordinates to a toroidal coordinate system, , that is defined such that,
then the expression for the modulus of the elliptic integral becomes,
which is a function of only one coordinate — the "radial" coordinate, — and the integral expression for the gravitational potential becomes,
If the configuration's density is constant, then, as is shown below, the integral over the angular coordinate variable, , can be completed analytically. Hence, the task of evaluating the gravitational potential (both inside and outside) of a uniformdensity, axisymmetric configuration having any surface shape has been reduced a problem of carrying out a single, line integration. We specifically illustrate how this approach can be used to evaluate the gravitational potential of a uniformdensity torus with a circular crosssection and aspect ratio, . Note: As we have reviewed separately, it appears as though over forty years ago C.Y. Wong (1973) successfully employed toroidal coordinates to derive a closedform expression for the potential (inside as well as outside) of a uniformdensity torus with a circular crosssection; see, for example, his Figure 7. That is to say, he was even able to carry out the final line integral in closed form. We have not (yet) been able to demonstrate that our expression for the potential of a circularcrosssection torus is identical to his. 
Exploring the Use of Toroidal Coordinates to Determine the Gravitational Potential 

As I have studied the structure and analyzed the stability of (both selfgravitating and nonselfgravitating) toroidal configurations over the years, I have often wondered whether it might be useful to examine such systems mathematically using a toroidal — or at least a toroidallike — coordinate system. Is it possible, for example, to build an equilibrium torus for which the density distribution is onedimensional as viewed from a wellchosen toroidallike system of coordinates?
I should begin by clarifying my terminology. In volume II (p. 666) of their treatise on Methods of Theoretical Physics, Morse & Feshbach (1953; hereafter MF53) define an orthogonal toroidal coordinate system in which the Laplacian is separable.^{1} (See details, below.) It is only this system that I will refer to as the toroidal coordinate system; all other functions that trace out toroidal surfaces but that don't conform precisely to Morse & Feshbach's coordinate system will be referred to as toroidallike.
I became particularly interested in this idea while working with Howard Cohl (when he was an LSU graduate student). Howie's dissertation research uncovered a Compact Cylindrical Greens Function technique for evaluating Newtonian potentials of rotationally flattened (especially axisymmetric) configurations.^{2,3} The technique involves a multipole expansion in terms of halfintegerdegree Legendre functions of the kind — see NIST digital library discussion — where, we have discovered, the argument of this special function is a function only of the radial coordinate of Morse & Feshbach's orthogonal toroidal coordinate system — see more on this, below.
Statement of the Problem
Expression for the Axisymmetric Potential
Cohl & Tohline (1999; hereafter CT99) derive an expression for the Newtonian gravitational potential in terms of a Compact Cylindrical Green's Function expansion. They show — see, for example, their equation (31) — that when expressed in terms of cylindrical coordinates, the potential at any meridional location, and , due to an axisymmetric mass distribution, , is
where,
is a differential area element in the meridional plane, and the dimensionless argument (the modulus) of the special function, , is,
Next, following the lead of CT99, we note that according to the Abramowitz & Stegun (1965),
where, the function is the complete elliptic integral of the first kind and, for our particular problem,









Hence, we can write,
As has been explained in an accompanying set of notes, this is precisely the same expression for the gravitational potential that A. Trova, J.M. Huré and F. Hersant (2012; MNRAS, 424, 2635) used in their study of the potential of selfgravitating, axisymmetric discs.
Our objective, here, is to examine whether or not it might be advantageous to transform this expression to one in which the double integral is performed on a toroidal, rather than a cylindrical, coordinate system.
Chosen Test Mass Distribution
For purposes of illustration, we will follow the lead of Trova, Huré & Hersant (2012) — see the lefthand panel of the following figure ensemble — and seek to determine the gravitational potential, both inside and outside, of a uniformdensity, equatorialplane torus whose (pink) meridional crosssection is exactly a circle. More specifically, as illustrated in our Figure 1 — see the righthand panel of the following figure ensemble — at all azimuthal angles, a crosssection through the (pink) torus is prescribed by the familiar algebraic expression for an offcenter circle, namely,



Everywhere inside this toroidal surface we set , that is, the density is uniform with the value, .
Figure 4 extracted without modification from p. 2640 of Trova, Huré & Hersant (2012)

Our Figure 1 

Notice that another offcenter circle — this one with a purple perimeter and otherwise white, rather than pink — appears in our Figure 1 diagram. In the discussion that follows, it will be used to represent the meridionalplane crosssection of one axisymmetric surface in an MF53 toroidalcoordinate system. Here we simply point out that this "surface" is also prescribed by an algebraic expression for an offcenter circle, namely,



Toroidal Coordinates
Click here to find a supporting, detailed illustration of Toroidal Coordinates and relevant integration limits.
Properties
Here we highlight certain properties and features of the MF53 toroidal coordinate system; more details can be found in a related set of our online notes. Most importantly in the context of our discussion, if (at all azimuthal angles) the origin of the toroidal coordinate system is placed at the cylindricalcoordinate location, the pair of orthogonal coordinates, , is related to the cylindrical coordinate pair, , via the expressions,
[MF53]  Wikipedia  











An offcenter circle — such as the white circle with purple perimeter depicted in our Figure 1 diagram — is generated if a value of the "radial" coordinate, , is chosen from within the range,
and held fixed while the "angular" coordinate, , is varied over the range,
Hereafter, we will refer to this = constant circle as a "circle." A circle of radius zero and, hence, the origin of the toroidal coordinate system is associated with the upper limiting value of the radial coordinate, namely, ; as the value of is decreased monotonically, the radius of the circle (for example, the circle of radius, , in our Figure 1) steadily grows; and the radius of this circle becomes infinite at the radial coordinate's other limiting value, .
In the plane, the location of the inner and outer edges of the toroidalcoordinate surface are determined by setting (inner) and (outer). Hence,






Hence, also, the (cylindrical) radial location of the "center" of each toroidalcoordinate surface — labeled in our Figure 1 — is given by the expression,
and the surface's crosssectional radius — labeled in our Figure 1 — is given by the expression,
This last expression quantifies, and its simplicity reinforces, our earlier statement; that is, as the value of is decreased monotonically, the radius of the circle, , steadily grows. The nexttolast expression makes it clear, as well, that grows larger and, therefore, the location of the center of a circle shifts farther away from the symmetry axis as the value of is decreased. Notice that, for any offcenter circle, the ratio of these to lengths gives the value of the toroidalcoordinate system's dimensionless "radial" coordinate, that is,



Notice, furthermore, that there is a particular combination of these two lengths that is independent of , namely,



This is a manner in which one can determine the radial position, , of the origin of the toroidal coordinate system that could legitimately be associated with any particular offcenter circle, such as the white circle with a purple perimeter drawn in our Figure 1.
Connection With the Physical Problem
Earlier, we stated that the offcenter circle with purple perimeter displayed in Figure 1 is prescribed by the algebraic expression,



Let's plug in the "toroidalcoordinate" expressions for each parameter that appears on the lefthand side of this relation and see whether, after simplification, it reduces to the righthand side.
LHS 




















This, indeed, equals the righthand side of the relation, which is, . It all nicely checks out!
Next, taking a hint from the EUREKA! moment recorded in our accompanying notes, let's rewrite the function in terms of toroidal rather than cylindrical coordinates, where is the argument of the special function, , that appears in the above definition of . More specifically, let's assume that the coordinate location at which the gravitational potential is to be evaluated, , is taken to be the cylindricalcoordinate location of the origin of the toroidal coordinate system, . Given this association, we can write,



























Hence, when the function, , is rewritten in terms of the elliptic integral of the first kind, , the modulus of can be written as,









This is the key result motivating the use of a toroidal coordinate system to evaluate the gravitational potential: When expressed in an appropriately defined toroidal coordinate system, the modulus of the special function is a function of one, rather than two, spatial coordinates. This gives some hope that the integral over the second (angular) coordinate, , can be completed analytically, giving rise to an expression for the gravitational potential whose evaluation only requires numerical integration over a single (radial) coordinate, .
Finally, drawing on discussion in our accompanying set of notes, we recognize that, expressed in terms of toroidal coordinates, the differential area element in the meridional plane is,



Putting everything together, then, the (indefinite) integral expression for , expressed in terms of toroidalcoordinates, is,









Making the substitution,
,
gives,



where we have now explicitly introduced four parameters to set definite limits on the nested pair of integrations. Making the additional substitution,
,
gives,



where, written in terms of ,



Perform Angular Integration Over Test Mass Distribution
While, in principle, the pair of nested integrals must be carried out over all of space using the limits,
and
in practice, the limits should be set to reflect the volume boundaries of the mass distribution that is responsible for generating the gravitational potential. Here we present analytic expressions for these limits in the case of the (pink) toroidal test mass distribution specified above; details regarding the derivation of these expressions can be found in our accompanying set of notes.
Identifying Limits of Integration
Figure 2 

The animation shown here in Figure 2 builds upon the configuration displayed in our Figure 1, above. It shows a meridional crosssection through the selected (pink) uniformdensity, toroidal mass distribution, whose geometric properties are fully determined by specifying values for and . (For the example illustrated in Figure 2, we have specified the same values used by Trova, Huré & Hersant (2012) to construct their Figure 4, as reprinted above; namely, and .) Throughout the Figure 2 animation sequence, these two parameters have been fixed — thereby fixing the properties of the (pink) torus — and, in addition, we have fixed the location of the origin of a toroidal coordinate system — identified by the redfilled circular dot. We explicitly associate this coordinatesystem origin with the (cylindrical) coordinates of the point in space at which we choose to evaluate the gravitational potential, namely, . [For illustration purposes, in Figure 2 we have set .]
While the values of the four primary model parameters are held fixed, Figure 2 depicts in a quantitatively precise manner how the size of a = constant toroidal surface (the offcenter circle traced by the sequence of black dots) varies as the value of the radial coordinate, , is varied. In each frame of the animation sequence, the value of that was used to define the (black) circle is printed in the lowerright corner of the image; additional quantitative details associated with each animation frame can be obtained from the table titled, "Example 2" in our accompanying notes.
As the value of is varied from large values (small black circles) to smaller values (larger black circles), there is a maximum value, , at which the circle first makes contact with the (pink) equatorialplane torus, and there is a minimum value, , at which it makes its final contact. These are the limiting values of the toroidal radial coordinate to be used in the integration that produces . At all values within the parameter range,
the circle intersects the surface of the torus in two locations, defined by two different values of the associated angular coordinate, . In each frame of the animation, points of intersection are marked with small yellow diamonds; the coordinates of these points of intersection are listed in the table associated with example 2, in our accompanying notes. For each relevant value of , these are the limiting values of the toroidal angular coordinate to be used in the integration that produces . It should be realized that, at the first and final points of contact, the two values of are degenerate.
Notice in the animation that, while the origin of the selected toroidal coordinate system (the filled red dot) remains fixed, the center of the circle does not remain fixed. In order to highlight this behavior, the location of the center of the circle has been marked by a filled, lightblue square and, in keeping with the earlier Figure 1 diagram, a vertical, lightblue line connects this center to the equatorial plane of the cylindrical coordinate system.
Associated Analytic Expressions
We define the following terms that are functions only of the four principal model parameters, , and therefore can be treated as constants while carrying out the pair of nested integrals that determine :

Then,
and
which implies,
and
Also let,
in which case, and,
and define,

Figure 3 

Then, for each value of the radial coordinate, , within the range,
the limiting values of the angular coordinate are,



This provides the desired analytic expression for the two limits of integration over the angular coordinate that must be carried out in order to evaluate the key function, , in our determination of the gravitational potential at the cylindricalcoordinate location, .
We should note that at both limits — that is, at the two values of for which the functions and have been evaluated, immediately above — the ratio under the squareroot,
Hence, at both radialcoordinate limits, the two angularcoordinate limits, , are degenerate. Furthermore, we have discovered that the value of this degeneracy angle is exactly the same at both radialcoordinate limits; specifically, we find that, for our chosen testmass distribution,
This means that exactly the same constant, toroidalcoordinate "line" intersects the surface of the (pink) testmass torus at both the point of first contact and the point of final contact . In order to illustrate this, this particular constant "line" has been traced by a sequence of red dots in Figure 3. As just described, it passes smoothly through the points on the surface of the (pink) torus where the circles make first (small blackdotted circle) and final (large blackdotted circle) contact.
Integrate
Because the density is uniform throughout our test model torus, it can be pulled out of both integrals. In combination with the limits of integration just derived, we can therefore write,



Now, using WolframAlpha's online integrator, we find …
Hence, the inner integral over our toroidal system's angular coordinate gives,















This last expression exactly matches the expression found in the integral tables published by Gradshteyn & Ryzhik (1965) — specifically, relation 2.575(5) — after it is realized that the second argument of both elliptic integral functions, as published in GR65, is the squareroot of the parameter, , provided by WolframAlpha.
We now have to deal with the realization that, for our particular problem, , which means that the second argument, , of both elliptical integral functions is negative (if not zero). Following the discussion provided at mymathlib.com, we define,
and
in which case,
and we can write,






Hence,















We have therefore succeeded in deriving an expression for the gravitational potential of a uniformdensity torus that requires numerical integration over only one dimension, that is, the "radial" dimension, . The final, combined expression is,

Code Layout
Here are some suggestions related to the development of a computer program that can perform the onedimensional integral over .
STEP 1: Specify numerical values for torus parameters … varpi_t & r_t STEP 2: Associate the parameters "a" and "Z0" with the cylindrical coordinate location (R*,Z*) at which gravitational potential is to be evaluated. STEP 3: Given the numerical values for these four parameters, calculate "kappa", "C", "beta_plus", and "beta_minus" as prescribed by the following algebraic relations …

STEP 4: The (numerical values of the) limits of integration are set by the following two expressions …
and
STEP 4 (cont.): Establish a 1D grid with, say, 100 zones, assigning values of xi_1 that are equally spaced between these two limiting values; let "delta" be the grid spacing. (Perhaps use logarithmic spacing.) The midpoint value of the coordinate location of the accompanying 99 grid cells should also be determined. STEP 5: For each of the 99 gridcell coordinate values, "xi1", calculate the following parameter, or aggregateterm, values:







Special Case
In the context of the (pink torus) testmass distribution being discussed here, a review of the properties of a toroidal coordinate system highlights one particularly interesting (cylindrical) coordinate location at which the gravitational potential should be evaluated:
By placing the origin of the toroidal coordinate system at this and only this location — that is, by setting,



— the limits of integration can be specified very simply: There is only one circle that intersects the surface of the torus and it does so by, not simply intersecting, but by perfectly aligning with the surface. This perfect overlap with the surface of the torus happens for the coordinate circle associated with,
All circles larger than this one — that is, all circles corresponding to values of
— lie entirely outside of the (pink) torus and therefore need not be included in the integration to determine the gravitational potential; while all circles smaller than this one — that is, all circles corresponding to values of
— lie entirely inside of the (pink) torus. Hence, in this special case the radialcoordinate integration limits are,
and
Furthermore, because the corresponding circles fall entirely inside (or on the surface of) the torus for all values of in this latter range, the integral over the "angular" toroidal coordinate, , will have the simple limits,
and
(More practically, we will integrate from zero to , and double the result.) Hence, the inner integral over our toroidal system's angular coordinate gives,












So, we have,









where,



TO BE DONE:
This expression should be evaluated at , in which case it appears as though it should match our "special case" result for all values of the polar angle, . 
Let's see if the expression simplifies somewhat by adopting a variable,
Making this substitution, and defining,
we have,






where,



Total Mass
How do I know whether or not I have made mistakes in these derivations? Aside from the published work of Trova, Huré & Hersant (2012), how do I know what the correct answer for the gravitational potential of a uniformdensity torus is?
A degree of assurance can be drawn by carrying out a similar pair of nested integrations to determine the total mass (or volume) of the torus that is defined by our chosen testmass distribution. As can be found in many mathematics handbooks, or online, the volume of our defined torus should be,
Using Cylindrical Coordinates
Total Volume
Let's start by integrating the threedimensional volume element, , over the azimuthal angle to obtain an expression for the twodimensional differential volume element that is written in terms of the meridionalplane differential area, , as used in our definition of , above. Using the notation of MF53 (but employing the opposite sign convention from them), in cylindrical coordinates we have,












where, as before in cylindrical coordinates, . From this we obtain,






Now in order to finish the volume integration, we need the limits of integration in the meridional plane. These can be obtained from the above algebraic description of the (pink) testmass torus as an offcenter circle. Specifically, we have,


















where we have carried out the second integration using the WolframAlpha online integral calculator:
This is the answer for the volume of a torus that we expected. Good!
Volume with Cropped Top
During my development of a computer program to integrate over the volume of a circular torus while using toroidal coordinates, I decided that it would be useful to determine the volume of a torus that has a "cropped top" as illustrated in the accompanying diagram, shown here on the right. Let's evaluate this "croppedtop" volume using cylindrical coordinates so that we will know what the correct answer is when developing an integration scheme using toroidal coordinates. Specifically, let's determine the volume of the "green" portion for a specified value of .
The radial integration will be evaluated between the limits: (lower) and (upper) , where,
And the limits on the vertical integration will be: (lower) and, as before when integrating over the entire volume, (upper) With these limits, the volume integration gives,









where,












Carrying out this pair of integrations gives,



























Hence, the fractional volume is,
Analytic Expression for Green Volume  


REALITY CHECK: This should give a zero (green) volume if ; and the fractional volume should be onehalf if . In the former case, our expression gives,



as expected. And in the latter case we have,



which means that is indeed half of the total torus volume, as derived earlier.
Using Toroidal Coordinates with Special Alignment
If, instead, we use a toroidal coordinate system, we have (see p. 666 of MF53),












where, as employed above, a differential area element in the meridional plane of a toroidalcoordinate system is,



As in the case of the cylindrical coordinate system, because the torus is axisymmetric, integration over the azimuthal angular coordinate in a toroidal coordinate system gives,






Crosscheck: Wikipedia's Differential Volume Element  

where,
and, as above,
Hence, in agreement with the expression derived using the notation of MF53,
