Editing AxisymmetricConfigurations/PoissonEq (section)

==Other Approaches==

===Deupree (1974)===

<table border="0" cellpadding="5" width="100%" align="center"><tr><td align="left">
<table border="1" cellpadding="5" align="left"><tr>
<td align="center" bgcolor="red">&nbsp; 2D &nbsp;</td>
<td align="center" bgcolor="lightgreen">&nbsp; Sph &nbsp;</td>
<td align="center" bgcolor="yellow">&nbsp; Tor &nbsp;</td>
</tr></table>
</td></tr></table>

<table border="0" cellpadding="3" align="center" width="75%">
<tr><td align="left">
<font color="#009999">
"The fully nonlinear, nonradial, nonadiabatic calculation of stellar oscillations has not as yet been attempted by anyone, to the best of the author's knowledge.  Only Deupree (1974b, 1975) so far seems to have taken some steps in this direction.  He has carried out numerical calculations of the nonlinear axisymmetric oscillations in the adiabatic ([http://adsabs.harvard.edu/abs/1974ApJ...194..393D Deupree 1974b]) and nonadiabatic ([http://adsabs.harvard.edu/abs/1975ApJ...198..419D Deupree 1975]) approximations, with apparently encouraging results."</font>
</td></tr>
<tr><td align="right">
&#8212; Drawn from &sect;3.5 of the review article by [http://adsabs.harvard.edu/abs/1976ARA%26A..14..247C J. P. Cox (1976, ARAA, 14, 247 - 273)] 
</td></tr></table>

Apparently, [http://adsabs.harvard.edu/abs/1974ApJ...194..393D R. G. Deupree (1974, ApJ, 194, 393 - 402)] was the first astronomer to employ self-gravitating, numerical hydrodynamic techniques to model ''nonlinear, nonradial stellar pulsations.''  He chose to carry out his simulations on a spherical coordinate mesh with his earliest simulations being restricted to the examination of 2D (axisymmetric) configurations.  The outer boundary of his computational grid was identified as (see his &sect;IIb) <font color="darkgreen">a spherical surface completely exterior to the star.</font>

====His Derived Expression for the Boundary Potential====
In designing an algorithm to solve the Poisson equation, [http://adsabs.harvard.edu/abs/1974ApJ...194..393D Deupree (1974)] considered <font color="darkgreen">the major difficulty [to be] the evaluation of the potential boundary conditions.</font>  He decided to determine the boundary potential via an evaluation of the ''integral representation'' of the Poisson equation; <font color="darkgreen">once the boundary conditions [were] specified, the potential inside the star [was] evaluated by [employing] the [http://adsabs.harvard.edu/abs/1964ApJ...139..306H Henyey method]</font> to solve the ''differential representation'' of the Poisson equation. The assumption of azimuthal symmetry meant that, for example, the density was specified at various meridional-plane locations, <math>~\rho(r,\theta)</math>, with <font color="darkgreen">each zone [being] considered as a uniform density circular ring.</font>  He argued that <font color="darkgreen">the potential of [each such] ring at any point on the spherical boundary [could] easily be evaluated by</font> treating it as an infinitesimally thin ring of mass, <math>~\delta M = 2\pi a \rho ~\delta A</math> &#8212; where, <math>~\delta A</math> <font color="darkgreen">is the cross-sectional area of the [grid] zone, and <math>~a</math> is the radius of [that] ring</font> &#8212; then drawing upon the analytic analysis described by [https://archive.org/details/foundationsofpot033485mbp Kellogg (1929)] or, equivalently, by [https://www.amazon.com/Theory-Potential-W-D-Macmillan/dp/0486604861/ref=sr_1_2?s=books&ie=UTF8&qid=1503444466&sr=1-2&keywords=the+theory+of+the+potential MacMillan (1958; originally 1930)].  Then, <font color="darkgreen">the total potential at the boundary is the sum of the potentials from all of the [meridional-plane] zones.</font>

----

<table border="0" align="center" cellpadding="5" width="80%"><tr><td align="left">
Note that, over the years, a number of other research groups have adopted this same approach to evaluate the gravitational potential of axisymmetric mass distributions &#8212; that is, summing up the potential contributions due to an ensemble of "infinitesimally thin rings."  But none appear to have recognized Deupree's earlier, pioneering investigation.  See, for example, the brief reviews that we have written regarding the related investigations by:  [[#Stahler_.281983.29|Stahler (1983)]], [[Apps/DysonWongTori#Bannikova_et_al._.282011.29|Bannikova et al. (2011)]], and [[Apps/DysonWongTori#Fukushima_.282016.29|Fukushima (2016)]].
</td></tr></table>

----

More specifically, equation (16) from [http://adsabs.harvard.edu/abs/1974ApJ...194..393D Deupree (1974)] states that the differential contribution to the boundary potential due to each infinitesimally thin ring is,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\delta\Phi_B</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
- \frac{2G(2\pi a\rho ~\delta A)}{\pi p} \int_0^{\pi/2} \frac{d\psi}{[\cos^2\psi + (q/p)^2\sin^2\psi]^{1 / 2}} \, ,
</math>
  </td>
</tr>
</table>
</div>

where, <math>~q</math> is the shortest distance between the ring and the boundary point, and <math>~p</math> is the longest distance.   (Note that, in order to match our conventions, we have inserted the negative sign along with a leading factor of <math>~G</math>.)  Referencing our accompanying [[Apps/DysonWongTori#DeupreeReference|detailed analysis of the potential due to a thin ring]], and adopting the variable mappings, <math>~p \leftrightarrow \rho_1</math> and <math>~q \leftrightarrow \rho_2</math>, we see that Deupree's expression is indeed identical to the expression for the potential derived by [https://www.amazon.com/Theory-Potential-W-D-Macmillan/dp/0486604861/ref=sr_1_2?s=books&ie=UTF8&qid=1503444466&sr=1-2&keywords=the+theory+of+the+potential MacMillan (1958)].
<!--
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\Phi</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- \biggl( \frac{2GM}{\pi} \biggr) \int_0^{\pi/2} \frac{d\omega}{[\rho_1^2 \cos^2\omega + \rho_2^2\sin^2\omega]^{1 / 2}} \, .</math>
  </td>
</tr>
</table>
</div>
-->
Recognizing, [[Apps/DysonWongTori#RingPotential|as did MacMillan]], that the definite integral in this expression is related to the complete elliptic integral of the first kind, and introducing the ratio of lengths, <math>~c \equiv p/q</math>, Deupree's expression for the (differential contribution to the) potential can be rewritten as,

<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\delta\Phi_B</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~-  \frac{2G (2\pi a\rho ~\delta A) c}{\pi p} \int_0^{\pi/2} \frac{d\psi}{ \sqrt{1 - k^2 \sin^2\psi }} </math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~-  \biggl[ \frac{2G (2\pi a\rho ~\delta A) c}{\pi p} \biggr] K(k)  \, ,</math>
  </td>
</tr>
</table>
where, <math>~K(k)</math> is the [http://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html complete elliptic integral of the first kind] for the modulus, <math>~k = \sqrt{1-c^2}</math>.

====Comparison With Other Related Derivations====
Now, given that Deupree chose to construct and evolve his models using a spherical coordinate system, he would have specified the relevant lengths, <math>~p</math> and <math>~q</math>, and each ring's differential cross-section, <math>~\delta A</math>, in terms of spherical coordinates.  In an effort to more clearly illustrate the connection between Deupree's expression for the (differential contribution to the) boundary potential and the expression for the boundary potential that we have [[#For_Axisymmetric_Systems|derived above for axisymmetric systems]], we will insert expressions for these terms that apply, instead, to a cylindrical-coordinate mesh.

Following the same line of reasoning as has been presented in our [[Apps/DysonWongTori#CylindricalLocation|accompanying discussion of MacMillan's work]], if the meridional-plane locations of the infinitesimally thin ring and the desired point on the boundary are, respectively, <math>~(\varpi^',z^')</math> and <math>~(\varpi,z)</math>, we see that,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~k = \sqrt{1 - c^2}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~ \biggl[\frac{4\varpi \varpi^'}{(\varpi + \varpi^')^2 + (z - z^')^2 }\biggr]^{1 / 2} \, ,</math>
  </td>
</tr>
</table>
</div>
and,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\frac{c}{p}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\biggl[(\varpi + \varpi^')^2 + (z - z^')^2  \biggr]^{- 1 / 2} = \frac{k}{\sqrt{4\varpi \varpi^'}} \, .</math>
  </td>
</tr>
</table>
</div>
Hence &#8212; after acknowledging that, in cylindrical coordinates, the radius of each "infinitesimally thin ring" is, <math>~a = \varpi^'</math>, and the differential cross-section of each ring is, <math>~\delta A = \delta\varpi^' \delta z^'</math> &#8212; Deupree's expression for the (differential contribution to the) potential may be rewritten as,

<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\delta\Phi_B(\varpi,z)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- \frac{G (2\pi a \rho ~\delta A) }{\pi } \biggl[ \frac{k }{\sqrt{\varpi \varpi^'}} \biggr]  K(k) </math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- \frac{2G (\delta M) }{\pi } \cdot \frac{K(k) }{ \sqrt{(\varpi + a)^2 + (z - z^')^2 }}  \, ,  </math>
  </td>
</tr>
</table>
or it may be rewritten as,

<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\delta\Phi_B(\varpi,z)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- \frac{G (2\pi a \rho ~\delta A) }{\pi } \biggl[ \frac{k }{\sqrt{\varpi \varpi^'}} \biggr]  K(k) </math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- \frac{2G }{\sqrt{\varpi}} \biggl[ \delta\varpi^' \delta z^' \rho(\varpi^', z^') \sqrt{\varpi^'} k K(k) \biggr] \, .</math>
  </td>
</tr>
</table>

Notice that the first of these two rewritten expressions aligns perfectly with our "[[Appendix/EquationTemplates#Other_Equations_with_Assigned_Templates|key equation]]" that gives the gravitational potential of an axisymmetric torus in the thin ring (TR) approximation, namely,
<table border="0" align="center" cellpadding="10"><tr><td align="center">
{{ Math/EQ TRApproximation }}
</td>
<td align="center" rowspan="2">[[File:FlatColorContoursCropped.png|225px|link=Apps/DysonWongTori#ThinRingContours]]</td>
</tr></table>
(See our [[Apps/DysonWongTori#ThinRingContours|accompanying discussion]] for more information on the meridional-plane contour plot that is displayed to the right of this equation.)  Next, referring back to the expression that was [[#For_Axisymmetric_Systems|derived above for axisymmetric systems from a toroidal-function-based Green's function]], namely,

<table border="0" align="center">
<tr>
  <td align="right">
<math>~ \Phi_B(\varpi,z)\biggr|_{2D}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~ 
-\frac{2G}{\sqrt{\varpi}} \int\limits_{\varpi^'} \int\limits_{z^'} d\varpi^' dz^' \rho(\varpi^',z^') \sqrt{\varpi^'}  
\mu K(\mu) \, , 
</math>
  </td>
</tr>
<tr>
  <td align="center" colspan="3">
[http://adsabs.harvard.edu/abs/1999ApJ...527...86C H. S. Cohl &amp; J. E. Tohline (1999)], p. 89, Eqs. (31) &amp; (32)
  </td>
</tr>
</table>
where,
<table border="0" align="center">
<tr>
  <td align="right">
<math>~\mu </math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\biggl[ \frac{4\varpi \varpi^'}{(\varpi + \varpi^')^2 + (z - z^')^2} \biggr]^{1 / 2}
\, ,
</math>
  </td>
</tr>
</table>
we see that the second of these rewritten expressions for Deupree's <math>~\delta\Phi_B</math> aligns perfectly with our derived expression for the differential contribution to the potential of any axisymmetric mass distribution.  It is therefore fair to say that the expression that Deupree used to determine the gravitational potential along the boundary of his modeled configurations is derivable from a 3D Green's function that is written in terms of ''toroidal functions''.

===Cook (1977)===

<table border="0" cellpadding="5" width="100%" align="center"><tr><td align="left">
<table border="1" cellpadding="5" align="left"><tr>
<td align="center" bgcolor="red">&nbsp; 3D &nbsp;</td>
<td align="center" bgcolor="lightgreen">&nbsp; Cyl &nbsp;</td>
<td align="center" bgcolor="yellow">&nbsp; Cyl &nbsp;</td>
</tr></table>
</td></tr></table>

Key results from [https://www.osti.gov/servlets/purl/7294335 T. L. Cook's (1977) doctoral dissertation] &#8212; titled, ''Three-Dimensional Dynamics of Protostellar Evolution'' &#8212; were published as [http://adsabs.harvard.edu/abs/1978ApJ...225.1005C T. L. Cook &amp; F. H. Harlow (1978, ApJ, 225, 1005 - 1020)]. Cook's three-dimensional hydrodynamic simulations were conducted on a cylindrical-coordinate mesh.  While few details of the technique used to solve the Poisson equation are provided in the ApJ article, &sect;II.B of [https://www.osti.gov/servlets/purl/7294335 Cook's (1977) dissertation] explains that values of the gravitational potential across the ''interior'' regions of the mesh were obtained by solving the ''differential representation'' of the Poisson equation, subject to the boundary conditions at each boundary point, calculated by "performing a numerical integration over all mass points" using the ''integral representation'' of the Poisson equation, namely,
<table border="0" align="center">
<tr>
  <td align="right">
<math>~ \Phi_B(\varpi,\phi,z)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~ 
-G \int \frac{1}{|\vec{x}^{~'} - \vec{x}|} ~\rho(\vec{x}^{~'}) d^3x^' \, ,
</math>
  </td>
</tr>
</table>
with,
<table border="0" cellpadding="5" align="center">
<tr>
  <td align="right">
<math>~\frac{1}{|\vec{x}^{~'} - \vec{x}|}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\biggl[ \varpi^2 + (\varpi^')^2 - 2\varpi \varpi^' \cos|\phi - \phi^'| + (z - z^')^2 \biggr]^{- 1 / 2} \, .
</math>
  </td>
</tr>
<tr>
  <td align="center" colspan="3">
[https://www.osti.gov/servlets/purl/7294335 Cook (1977)], p. 15, Eq. (II-16)<br />
See also: [http://adsabs.harvard.edu/abs/2007AmJPh..75..724S Selvaggi, Salon &amp; Chari (2007)] &sect;II, eq. (2)
  </td>
</tr>
</table>

(Note that the leading factor of the gravitational constant, <math>~G</math>, does not appear explicitly in Cook's equation II-16, although it should.)  Presumably the coordinate locations of the boundary cells, <math>~(\varpi,\phi,z)</math>, were in no case coincident with the coordinate locations of any of the ''interior'' grid cells, <math>~(\varpi^',\phi^',z^')</math>, so there was no danger that the integration would encounter a singularity as a consequence of the distance, <math>~|\vec{x}^{~'} - \vec{x}|</math>, being zero.

===Norman and Wilson (1978)===

<table border="0" cellpadding="5" width="100%" align="center"><tr><td align="left">
<table border="1" cellpadding="5" align="left"><tr>
<td align="center" bgcolor="red">&nbsp; 3D &nbsp;</td>
<td align="center" bgcolor="lightgreen">&nbsp; Cyl &nbsp;</td>
<td align="center" bgcolor="yellow">&nbsp; Sph &nbsp;</td>
</tr></table>
</td></tr></table>

Following the discussion presented in &sect;4.1 of [https://archive.org/details/ClassicalElectrodynamics2nd Jackson (1975)], it is sometimes useful to rewrite [[#PotentialA|''Form A'' of the boundary potential]] as,
<table border="0" align="center">
<tr>
  <td align="right">
<math>~ \Phi_B(r,\theta,\phi)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~ 
-4\pi G  
\sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{+\ell} \frac{q_{\ell m}}{r^{\ell+1}} \frac{Y_{\ell m}(\theta,\phi)}{(2\ell+1)}  \, ,
</math>
  </td>
</tr>
</table>

<span id="MultipoleMoments">where we have introduced what is commonly referred to as the,
<table border="0" cellpadding="5" align="center">

<tr>
<td align="center" colspan="3">
<font color="#770000">'''Multipole Moments of the Mass Distribution'''</font>
</td>
</tr>
<tr>
  <td align="right">
<math>~q_{\ell m}</math>
  </td>
  <td align="center">
<math>~\equiv</math>
  </td>
  <td align="left">
<math>~
\int (r^')^\ell  Y_{\ell m}^*(\theta^', \phi^')
~\rho(r^', \theta^', \phi^') d^3x^' 
\, .
</math>
  </td>
</tr>
<tr>
  <td align="center" colspan="3">
[https://archive.org/details/ClassicalElectrodynamics2nd Jackson (1975)], p. 137, Eq. (4.3)
  </td>
</tr>
</table>

When written explicitly in terms of cartesian coordinates &#8212; see [[#Ylm|Table 2, below]], for each of the relevant <math>~Y_{\ell m}</math> expressions &#8212; the first few of these moments have the functional representations [[#Multipole_Moments_of_the_Mass_Distribution|derived below]].  For the cases that correspond to positive values of the index, <math>~m</math>, the set of ''multipole moment'' expressions that have been included in our [[#qlm|Table 3 summary]] exactly matches the set of expressions presented as equations (4.4), (4.5), and (4.6) in [https://archive.org/details/ClassicalElectrodynamics2nd Jackson (1975)].  (Jackson's equation 4.7 explains how to map these expressions to the cases corresponding to negative values of the index, <math>~m</math>.)  

<span id="PhiSeriesExpansion">Let's now look at various terms in the summed expression for the boundary potential</span> with each term expressing the contribution for a separate value of the index, <math>~\ell</math>.  Isolating the first three terms from all the rest, for example, we have,

<table border="0" align="center">
<tr>
  <td align="right">
<math>~ \Phi_B(r,\theta,\phi)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~ \Phi_B(r,\theta,\phi) \biggr|_{\ell=0} +~ \Phi_B(r,\theta,\phi)\biggr|_{\ell=1} + ~\Phi_B(r,\theta,\phi)\biggr|_{\ell=2} 
-~4\pi G  
\sum_{\ell=3}^{\infty} \sum_{m=-\ell}^{+\ell} \frac{q_{\ell m}}{r^{\ell+1}} \frac{Y_{\ell m}(\theta,\phi)}{(2\ell+1)}  \, .
</math>
  </td>
</tr>
</table>
We recognize, first, that as we consider boundary points that lie farther and farther away from the mass distribution, the magnitude of the first <math>(\ell = 0)</math> term drops off as <math>~r^{-1}</math> &#8212; the expected behavior of the potential outside of a point mass; the second <math>~(\ell=1)</math> term drops off as <math>~r^{-3}</math>; and the third <math>~(\ell=2)</math> term drops off as <math>~r^{-5}</math>.  Below, we have evaluated in more detail the behavior of these first three terms.  This evaluation gives us the,
<table border="0" align="center">
<tr>
<td align="center" colspan="3">
<font color="#770000">'''Boundary Potential Written in Terms of Multipole Moments '''</font>
</td>
</tr>
<tr>
  <td align="right">
<math>~ \Phi_B(r,\theta,\phi)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~ 
~-~~\frac{GM}{r}
~~- ~~\frac{GM}{r^3} \biggl[ \vec{x} \cdot \vec{x}_\mathrm{com} \biggr]
~~-~~\frac{G}{2r^5} \sum_{i=1}^3 \sum_{j=1}^3 Q_{i,j} \biggl[ x_i x_j \biggr]
~~-~~4\pi G  \sum_{\ell=3}^{\infty} \sum_{m=-\ell}^{+\ell} \frac{q_{\ell m}}{r^{\ell+1}} \frac{Y_{\ell m}(\theta,\phi)}{(2\ell+1)}  \, ,
</math>
  </td>
</tr>
</table>
where, as defined below, <math>~\vec{x}_\mathrm{com}</math> is the [[#Second_Term|center-of-mass location]], and <math>~Q_{i,j} </math> is the [[#QuadrupoleMomentTensor|traceless quadrupole moment tensor]].

If a modeled mass-density distribution, <math>~\rho( \vec{x}^{~'})</math>, has been configured such that the center-of-mass of the system coincides with the origin of the coordinate system &#8212; that is, if it has been configured such that <math>~\vec{x}_\mathrm{com} = 0</math> &#8212; then the second term in this series can be set to zero.  If, in addition, all terms having <math>~\ell \ge 3</math> are ignored because their values drop off rapidly with distance &#8212; specifically, inverse distance to the <math>~(2\ell + 1)</math> power &#8212; then a reasonably good approximation for the potential on the boundary of the modeled system is given by the expression,
<table border="0" align="center">
<tr>
  <td align="right">
<math>~ \Phi_B</math>
  </td>
  <td align="center">
<math>~\approx</math>
  </td>
  <td align="left">
<math>~ 
-\frac{GM}{r}
~-\frac{G}{2r^5} \sum_{i=1}^3 \sum_{j=1}^3 Q_{i,j} \biggl[ x_i x_j \biggr] \, .
</math>
  </td>
</tr>
</table>

This is precisely the relation that was adopted by [http://adsabs.harvard.edu/abs/1978ApJ...224..497N M. L. Norman &amp; J. R. Wilson (1978, ApJ, 224, 497 - 511)] when, in the context of star formation, they modeled the ''Fragmentation of Isothermal Rings''  &#8212; see specifically their equation (18).

===Boss (1980)===

<table border="0" cellpadding="5" width="100%" align="center"><tr><td align="left">
<table border="1" cellpadding="5" align="left"><tr>
<td align="center" bgcolor="red">&nbsp; 3D &nbsp;</td>
<td align="center" bgcolor="lightgreen">&nbsp; Sph &nbsp;</td>
<td align="center" bgcolor="yellow">&nbsp; Sph &nbsp;</td>
</tr></table>
</td></tr></table>

Building upon equations (5) and (7) from [http://adsabs.harvard.edu/abs/1980ApJ...236..619B A. P. Boss (1980, ApJ, 236, 619 - 627)] &#8212; see also Jackson's equation (3.58) &#8212; we can write,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\Phi(r,\theta,\phi)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\sum_{\ell=0}^\infty \sum_{m = -\ell}^{\ell} \Phi_{\ell m} ( r ) Y_{\ell m}(\theta,\phi) \, ,
</math>
  </td>
</tr>
</table>
</div>
and,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\rho(r,\theta,\phi)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\sum_{\ell=0}^\infty \sum_{m = -\ell}^{\ell} \rho_{\ell m} ( r ) Y_{\ell m}(\theta,\phi) \, ,
</math>
  </td>
</tr>
</table>
</div>
where the coefficients in the series expansion of <math>~\rho</math> can each be obtained from the known spatial density distribution via the integral expression,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\rho_{\ell m}(r)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\int d\Omega Y_{\ell m}^*(\theta,\phi) \rho(r,\theta,\phi) \, ,
</math>
  </td>
</tr>
</table>
</div>
and the differential solid angle is,
<div align="center">
<math>~d\Omega \equiv \sin\theta d\theta d\phi \, .</math>
</div>

===Stahler (1983)===

<table border="0" cellpadding="5" width="100%" align="center"><tr><td align="left">
<table border="1" cellpadding="5" align="left"><tr>
<td align="center" bgcolor="red">&nbsp; 2D &nbsp;</td>
<td align="center" bgcolor="yellow">&nbsp; Tor &nbsp;</td>
</tr></table>
</td></tr></table>

[http://adsabs.harvard.edu/abs/1983ApJ...268..155S Stahler (1983a)] used a [[AxisymmetricConfigurations/HSCF#Introduction|self-consistent-field]] technique to construct equilibrium sequences of rotationally flattened, isothermal gas clouds.  At each iteration step, the method that he adopted to evaluate the gravitational potential along the outer boundary of his computational mesh was essentially the same as the method used by [http://adsabs.harvard.edu/abs/1974ApJ...194..393D Deupree (1974)] &#8212; see our [[#Deupree_.281974.29|above description]] &#8212; to model the time-dependent behavior of pulsating stars.  It appears as though Stahler was unaware of Deupree's earlier development of this method to evaluate the boundary potential, as Deupree's (1974) paper is not among Stahler's list of references.  

<!--
In the chapter of this H_Book that focuses on a discussion of [[Apps/DysonWongTori|Dyson-Wong tori]], we have included the expression for the [[Apps/DysonWongTori#RingPotential|gravitational potential of a thin ring]] of mass, <math>~M</math>, that passes through the meridional plane at coordinate location, <math>~(\varpi^', z^') = (a, 0)</math>, as derived, for example, by [https://archive.org/details/foundationsofpot033485mbp O. D. Kellogg (1929)] and by [https://www.amazon.com/Theory-Potential-W-D-Macmillan/dp/0486604861/ref=sr_1_2?s=books&ie=UTF8&qid=1503444466&sr=1-2&keywords=the+theory+of+the+potential W. D. MacMillan (1958; originally, 1930)], namely,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\Phi(\varpi, z)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
- \biggl[ \frac{2GMc}{\pi\rho_1}\biggr] K(k)
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
- \biggl[ \frac{2GM}{\pi } \biggr]\frac{1}{\sqrt{(\varpi+a)^2 + z^2}} \times K(k) \, ,
</math>
  </td>
</tr>
</table>
</div>
where,

<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~k</math>
 </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\biggl[ \frac{4a\varpi }{(\varpi+a)^2 + z^2} \biggr]^{1 / 2}  \, .
</math>
  </td>
</tr>
</table>
</div>
-->
In describing [http://adsabs.harvard.edu/abs/1983ApJ...268..155S Stahler's (1983a)] method, we will first draw upon our "[[Appendix/EquationTemplates#Other_Equations_with_Assigned_Templates|key equation]]" that gives the gravitational potential of an axisymmetric torus in the thin ring (TR) approximation, namely,
<table border="0" align="center" cellpadding="10"><tr><td align="center">
{{ Math/EQ TRApproximation }}
</td>
<td align="center" rowspan="2">[[File:FlatColorContoursCropped.png|225px|link=Apps/DysonWongTori#ThinRingContours]]</td>
</tr></table>
(See our [[Apps/DysonWongTori#ThinRingContours|accompanying discussion]] for more information on the meridional-plane contour plot that is displayed to the right of this equation.)  Stahler has argued that a reasonably good approximation to the gravitational potential due to any extended axisymmetric mass distribution can be obtained by adding up the contributions due to many ''thin rings'' &#8212; <math>~\delta M(\varpi^', z^')</math>  being the appropriate differential mass contributed by each ring element &#8212; positioned at various meridional coordinate locations throughout the mass distribution.  According to his independent derivation, the differential contribution to the potential, <math>~\delta\Phi_g(\varpi, z)</math>, due to each differential mass element is (see his equation 11 and the explanatory text that follows it):
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\delta\Phi_g(\varpi,z)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
- \biggl[\frac{2G}{\pi \varpi^'}\biggr] \frac{\delta M}{[(\alpha + 1)^2 + \beta^2]^{1 / 2}}
\times K\biggl\{ \biggl[ \frac{4\alpha}{(\alpha+1)^2 + \beta^2} \biggr]^{1 / 2} \biggr\}
</math>
  </td>
</tr>
<!--
<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
- \biggl[\frac{2G}{\pi }\biggr] \frac{\delta M}{[(\varpi^' \alpha + \varpi^')^2 + (\varpi^' \beta)^2]^{1 / 2}}
\times K\biggl\{ \biggl[ \frac{4\alpha (\varpi^')^2}{(\varpi^' \alpha+\varpi^')^2 + (\varpi^' \beta)^2} \biggr]^{1 / 2} \biggr\}
</math>
  </td>
</tr>
-->
<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
- \biggl[\frac{2G}{\pi }\biggr] \frac{\delta M}{[(\varpi + \varpi^')^2 + (z^' - z)^2]^{1 / 2}}
\times K\biggl\{ \biggl[ \frac{4\varpi^' \varpi}{(\varpi +\varpi^')^2 + (z^' - z)^2} \biggr]^{1 / 2} \biggr\} \, .
</math>
  </td>
</tr>
</table>
</div>
Stahler's expression for each ''thin ring'' contribution is clearly the same as our "key equation" expression for <math>~\Phi_\mathrm{TR}</math> if the individual ring being considered cuts through the meridional plane at <math>~(\varpi^', z^') = (a, 0)</math>.

In the context of our [[#For_Axisymmetric_Systems|above discussion of a Green's function expression written in terms of toroidal functions]], we have shown that the exact integral expression for the gravitational potential due to any axisymmetric mass-density distribution, <math>~\rho(\varpi^', z^')</math>, is,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\Phi(\varpi, z)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- \frac{2G}{\varpi^{1 / 2}} \int\int (\varpi^')^{1 / 2} \mu K(\mu) \rho(\varpi^', z^') d\varpi^' dz^' \, ,</math>
  </td>
</tr>
</table>
</div>
where, 
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\mu^2</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~
\biggl[\frac{4\varpi^' \varpi}{(\varpi + \varpi^')^2 + (z^' - z)^2} \biggr] \, .
</math>
  </td>
</tr>
</table>
</div>
Recognizing that, for axisymmetric structures, the differential mass element is, <math>~dM^' = 2\pi \rho(\varpi^', z^') \varpi^' d\varpi^' dz^'</math>, this integral expression may be rewritten as,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~\Phi(\varpi, z)</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- \frac{2G}{\varpi^{1 / 2}} \int\int (\varpi^')^{1 / 2} \mu K(\mu) \biggl[ \frac{dM^'}{2\pi \varpi^'} \biggr]  </math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- \frac{G}{\pi} \int\int \biggl[ \frac{1}{\varpi^'\varpi}\biggr]^{1 / 2} \biggl[\frac{4\varpi^' \varpi}{(\varpi + \varpi^')^2 + (z^' - z)^2} \biggr]^{1 / 2} K(\mu) dM^'  </math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~- \frac{2G}{\pi} \int\int  \biggl[\frac{1}{(\varpi + \varpi^')^2 + (z^' - z)^2} \biggr]^{1 / 2} K(\mu) dM^'  \, .</math>
  </td>
</tr>
</table>
</div>
We see that our expression for the differential contribution to the potential exactly matches [http://adsabs.harvard.edu/abs/1983ApJ...268..155S Stahler's (1983a)].  It is therefore fair to say that Stahler's expression for the gravitational potential is derivable from a 3D Green's function that is written in terms of toroidal functions.  We note that, in his study of axisymmetric systems, Stahler made the decision to evaluate the gravitational potential both inside as well as outside of the mass distribution using the same integral expression.