Editing AxisymmetricConfigurations/PoissonEq (section)

=====In 3D=====
<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; Tor &nbsp;</td>
</tr></table>
</td></tr></table>
NOTE:  Throughout this chapter subsection, text that appears in a dark green font has been taken ''verbatim'' from [http://adsabs.harvard.edu/abs/1999ApJ...527...86C H. S. Cohl &amp; J. E. Tohline (1999)].

[[#MeshChoice|As mentioned above]], from the beginning of my research activities &#8212; following the lead of Black &amp; Bodenheimer &#8212; it has seemed reasonable to me that numerical simulations of time-evolving, rotationally flattened fluid systems should be carried out on a cylindrical, rather than cartesian, coordinate mesh.  When modeling rotationally flattened configurations, a cylindrical mesh has even seemed preferable to a ''spherical'' coordinate mesh because the "top" grid boundary (horizontal green-dashed line segment in [[#MeshChoice|Figure 1]]) can straightforwardly be dropped to a <math>~z</math>-coordinate location that is smaller than the <math>~\varpi</math>-coordinate location of the "side" grid boundary  (vertical green-dashed line segment in [[#MeshChoice|Figure 1]]), thereby reducing the number of grid cells &#8212; and, correspondingly reducing the cost of modeling the less interesting, ''vacuum'' region &#8212; outside of the fluid system.  [See, however, [[#Boss_.281980.29|Boss (1980)]] for an alternate point of view.]  At the same time, however, it has not seemed reasonable to determine the values of the potential along the (cylindrical-grid) boundary by adopting a Green's function that is expressed in terms ''spherical harmonics''. 

Over a period of approximately twenty years, off and on, my research group considered <font color="darkgreen">whether it might be advantageous in our numerical simulations to cast the Green's function in a cylindrical coordinate system.  The "familiar" expression for the cylindrical Green's function expansion can be found in variety of references (see [<b>[[Appendix/References#MF53|<font color="red">MF53</font>]]</b>]; [https://archive.org/details/ClassicalElectrodynamics2nd Jackson (1975)]</font>), and for convenience is [[#Familiar_Expression_for_the_Cylindrical_Green.27s_Function_Expansion|repeated below]].  <font color="darkgreen">It is expressible in terms of an infinite sum over the azimuthal quantum number <math>~m</math> and an infinite integral over products of Bessel functions of various orders multiplied by an exponential function. </font>  Note that [http://adsabs.harvard.edu/abs/1985ApJ...290...75V J. V. Villumsen (1985, ApJ, 290, 75 - 85)] attempted to solve the potential problem in this manner; <font color="darkgreen">he presents a technique in which each infinite integral over products of Bessel functions is evaluated numerically using a Gauss-Legendre integrator &hellip; He then emphasizes the obvious problem that, because of the infinite integrals involved, a calculation of the potential via this straightforward application of the familiar cylindrical Green's function expansion is numerically much more difficult than a calculation of the potential using a ''spherical'' Green's function expansion.</font>

<font color="red"><b>Eureka!</b></font> Via his dogged efforts and an extraordinarily in-depth investigation of this problem, [[Appendix/Ramblings/CCGF#Compact_Cylindrical_Green_Function_.28CCGF.29|in 1999 Howard S. Cohl discovered]] that, in cylindrical coordinates, the relevant Green's function can be written in a much more compact and much more practical form. Specifically, he realized that,
<table border="0" align="center" cellpadding="5">
<tr>
  <td align="right">
<math>~ \frac{1}{|\vec{x} - \vec{x}^{~'}|}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~ 
\frac{1}{\pi \sqrt{\varpi \varpi^'}} \sum_{m=-\infty}^{\infty} e^{im(\phi - \phi^')}Q_{m- 1 / 2}(\chi) \, ,
</math>
  </td>
</tr>
</table>
where,<br />
<div align="center"><math>~\chi \equiv \frac{\varpi^2 + (\varpi^')^2 + (z - z^')^2}{2\varpi \varpi^'} \, ,</math><br /><br />
[http://adsabs.harvard.edu/abs/1999ApJ...527...86C H. S. Cohl &amp; J. E. Tohline (1999)], p. 88, Eqs. (15) &amp; (16)<br />
See also: [http://adsabs.harvard.edu/abs/2007AmJPh..75..724S Selvaggi, Salon &amp; Chari (2007)] &sect;II, eq. (5)<br />
and the [https://dlmf.nist.gov/14.19#ii DLMF's definition of Toroidal Functions], <math>~Q_{m - 1 / 2}^{0}</math>
</div>
and, <math>~Q_{m - 1 / 2}</math> is the zero-order, half-(odd)integer degree, Llegendre function of the second kind &#8212; also referred to as a ''toroidal'' function of zeroth order; see [[#Toroidal_Functions|additional details, below]].  Hence, anywhere along the boundary of our cylindrical-coordinate mesh, a valid expression for the gravitational potential is,

<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>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~ 
-G \int \biggl\{ 
\frac{1}{\pi \sqrt{\varpi \varpi^'}} \sum_{m=-\infty}^{\infty} e^{im(\phi - \phi^')}Q_{m- 1 / 2}(\chi)
\biggr\}~
\rho(\varpi^',\phi^',z^') \varpi^' d\varpi^' d\phi^' dz^' 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~ 
-\frac{G}{\pi \sqrt{\varpi}} \int \rho(\varpi^',\phi^',z^') \sqrt{\varpi^'} d\varpi^' d\phi^' dz^'  
\sum_{m=0}^{\infty} \epsilon_m \cos[m(\phi - \phi^')] Q_{m- 1 / 2}(\chi) \, ,
</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. 88, Eq. (18)
  </td>
</tr>
</table>
where, <math>~\epsilon_m</math> is the Neumann factor, that is, <math>~\epsilon_0=1</math> and <math>~\epsilon_m = 2</math> for <math>~m\ge 1</math>.  Following this discovery, most of my research group's 3D numerical modeling of self-gravitating fluids has been conducted using ''Toroidal functions'' instead of ''Spherical Harmonics'' to evaluate the boundary potential on our cylindrical-coordinate meshes; see, for example, [http://adsabs.harvard.edu/abs/2002ApJS..138..121M P. M. Motl, J. E. Tohline &amp; J. Frank (2002)]; [http://adsabs.harvard.edu/abs/2005ApJ...625L.119O C. D. Ott, S. Ou, J. E. Tohline &amp; A. Burrows (2005)]; [http://adsabs.harvard.edu/abs/2006ApJ...643..381D M. C. R. D'Souza, P. M. Motl, J. E. Tohline, &amp; J. Frank (2006)]; and [http://adsabs.harvard.edu/abs/2012ApJS..199...35M D. C. Marcello &amp; J. E. Tohline (2012)].