Editing AxisymmetricConfigurations/PoissonEq (section)

===Determining Values of the Potential on the Mesh Boundary===
Let's determine the potential, <math>~\Phi_B</math>, at all points along the ''boundary'' of the cylindrical coordinate mesh by evaluating Table 1's ''integral representation'' of the Poisson equation. 

====Using a Spherical Harmonic Expansion====
=====Full Three-Dimensional Generality=====

Following the lead of [http://adsabs.harvard.edu/abs/1975ApJ...199..619B Black &amp; Bodenheimer (1975)], we will insert into this integral relation the Green's function expression for <math>~|\vec{x}^{~'}- \vec{x} |^{-1} </math> as given in terms of ''Spherical Harmonics'', <math>~Y_{\ell m}</math>, which in turn can be written in terms of ''Associated Legendre Functions.''  [[#Ylm|Table 2, below]], provides the primary details.

Written in the context of a spherical coordinate system 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>~ 
-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 
\sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{+\ell} \frac{4\pi}{2\ell+1} \biggl[  \frac{r_<^\ell}{r_>^{\ell+1}} \biggr] Y_{\ell m}^*(\theta^', \phi^') Y_{\ell m}(\theta,\phi)
~\rho(r^', \theta^', \phi^') d^3x^' 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~ 
-4\pi G  
\sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{+\ell} \frac{Y_{\ell m}(\theta,\phi)}{(2\ell+1)} \biggl[ \frac{1}{r^{\ell+1}}\int_0^r (r^')^\ell  Y_{\ell m}^*(\theta^', \phi^')
~\rho(r^', \theta^', \phi^') d^3x^' 
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
&nbsp;
  </td>
  <td align="left">
&nbsp; &nbsp; &nbsp; &nbsp; <math>~
+ r^\ell \int_r^\infty (r^')^{-(\ell+1)}  Y_{\ell m}^*(\theta^', \phi^')
~\rho(r^', \theta^', \phi^') d^3x^' \biggr] \, .
</math>
  </td>
</tr>
<tr>
  <td align="center" colspan="3">
[<b>[[Appendix/References#BT87|<font color="red">BT87</font>]]</b>], p. 66, Eq. (2-122)
  </td>
</tr>
</table>

If the distance from the origin, <math>~r</math>, of a boundary point (''i.e.'', any point lying along the dashed green lines in Figure 1) is greater than the distance from the origin, <math>~r^'</math>, of ''all'' of the (pink) mass elements, the second integral can be completely ignored.  We have, then, an expression that will henceforth be referred to as,
<!-- OLD COMMENT
In making this last step, we have moved the radial distance with its associated exponent, <math>~(r_>)^{-(\ell+1)}</math>, outside of the mass integral.  At the same time, we have left the alternate radial distance along with its associated exponent, <math>~(r_<)^\ell</math>, inside the integral and, accordingly, have labeled it with a "prime" to emphasize its association with the integral.  This has been done under the assumption that ''every'' (pink) mass element (tagged by a "primed" coordinate) lies closer to the coordinate origin than ''every'' point on the boundary (dashed green lines).  The integral must be split into two parts with the locations (inside or outside of the integral) of <math>~r_></math> and <math>~r_<</math> swapped in the second part of the integral if, in any case, the point on the boundary lies closer to the coordinate origin than ''any'' (pink) mass element(s).
-->
<div align="center" id="PotentialA">
<table border="0" align="center">
<tr>
<td align="center" colspan="3">
<font color="#770000">'''Form A of the Boundary Potential'''</font><br />
</td>
</tr>
<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{Y_{\ell m}(\theta,\phi)}{(2\ell+1)} \biggl[ \frac{1}{r^{\ell+1}}\int_0^r (r^')^\ell  Y_{\ell m}^*(\theta^', \phi^')
~\rho(r^', \theta^', \phi^') d^3x^' \biggr] \, .
</math>
  </td>
</tr>
<tr>
  <td align="center" colspan="3">
[https://archive.org/details/ClassicalElectrodynamics2nd Jackson (1975)], p. 137, Eq. (4.2)<br />
[<b>[[Appendix/References#BLRY07|<font color="red">BLRY07</font>]]</b>], p. 238, Eqs. (7.53) - (7.54)
  </td>
</tr>
</table>
</div>

Rewriting this expression for <math>~\Phi_B</math> in terms of cylindrical coordinates &#8212; which aligns with our chosen grid coordinate system &#8212; and admitting that in practice our summation over the index, <math>~\ell</math>, cannot extend to infinity, we have,

<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>~ 
-4\pi G 
\sum_{\ell=0}^{\ell_\mathrm{max}} \sum_{m=-\ell}^{+\ell} \frac{Y_{\ell m}}{(2\ell+1)} \biggl[ \varpi^2 + z^2 \biggr]^{-(\ell+1)/2} \int  Y_{\ell m}^* \biggl[ (\varpi^')^2 + (z^')^2 \biggr]^{\ell/2} 
~\rho(\varpi^', \phi^', z^') d^3x^' \, .
</math>
  </td>
</tr>
</table>

Note that, as a consequence of assuming that our configurations have equatorial-plane symmetry, the weighted integral over the mass distribution necessarily goes to zero anytime the sum of the two indexes, <math>~(\ell + m)</math>, is an odd number.  This is because, in each of these situations &#8212; again, see [[#Ylm|Table 2, below]] for details &#8212; the <math>~Y_{\ell m}</math> includes an overall factor of <math>~\cos\theta</math>, which necessarily switches signs between the two hemispheres.  After setting <math>~\ell_\mathrm{max}=4</math> &#8212; and dropping all terms in the summation for which the index sum, <math>~(\ell + m)</math>, is odd &#8212; this expression becomes precisely the relation that was used to determine the ''boundary'' values of the gravitational potential in our earliest set of simulations; see, for example, [http://adsabs.harvard.edu/abs/1978PhDT.........6T Tohline (1978)], [http://adsabs.harvard.edu/abs/1980ApJ...235..866T Tohline (1980)], and [http://adsabs.harvard.edu/abs/1980ApJ...242..209B Bodenheimer, Tohline, &amp; Black (1980)].

=====Simplification for 2D, Axisymmetric Systems=====
It is easy to show that this last expression for <math>~\Phi_B</math> &#8212; which has been used in our 3D simulations &#8212; is a ''generalization'' of the expression for <math>~\Phi_B</math> that was employed by [http://adsabs.harvard.edu/abs/1975ApJ...199..619B Black &amp; Bodenheimer (1975)] for 2D, axisymmetric simulations.  In axisymmetric systems, by definition, physical variables exhibit no variation in the azimuthal coordinate direction.  Hence, in the expression for <math>~\Phi_B</math>: 
* the azimuthal coordinate, <math>~\phi</math>, need not appear explicitly as an independent variable; 
* the index, <math>~m</math>, must be set to zero, so there is no summation over this index; and,
* every surviving spherical harmonic can be written more simply in terms of a Legendre function, namely,
<div align="center">
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>~Y_{\ell m} \rightarrow Y_{\ell 0}</math>
  </td>
  <td align="center">
<math>~=</math>
  </td>
  <td align="left">
<math>~\sqrt{\frac{(2\ell+1 )}{4\pi}} P_\ell(\chi) \, ,</math>
  </td>
  <td align="center">&nbsp; &nbsp; &nbsp; &nbsp; where,</td>
  <td align="left">
<math>~\chi \equiv \frac{z}{(\varpi^2 + z^2)^{1 / 2}}  \, .</math>
  </td>
</tr>
</table>
</div>
Note that the argument, <math>~\chi</math>, is still the spherical-coordinate expression, <math>~\cos\theta</math>, but here it has been written in terms of cylindrical coordinates. We have, therefore,

<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>~ 
- G 
\sum_{\ell=0}^{\ell_\mathrm{max}}  P_\ell(\chi) \biggl[ \varpi^2 + z^2 \biggr]^{-(\ell+1)/2} \int  P_\ell(\chi^') \biggl[ (\varpi^')^2 + (z^')^2 \biggr]^{\ell/2} 
~\rho(\varpi^', z^') d^3x^' \, ,
</math>
  </td>
</tr>
</table>
where, <math>~d^3x^' = 2\pi \varpi^' d\varpi^' dz^'</math>.  This is precisely the same as equation (5) from [http://adsabs.harvard.edu/abs/1975ApJ...199..619B Black &amp; Bodenheimer (1975)]; see also, equations (8) and (9) in [http://adsabs.harvard.edu/abs/1999ApJ...527...86C Cohl &amp; Tohline (1999)].

====Using Toroidal Functions====

=====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)].

=====For Axisymmetric Systems=====

As was done, [[#Simplification_for_2D.2C_Axisymmetric_Systems|above]], in the context of our discussion of a spherical-harmonics-based expression for the boundary potential, let's consider how this toroidal-function-based expression for the boundary potential can be simplified when examining 2D (axisymmetric) rather than fully 3D systems.  In axisymmetric systems, by definition, physical variables exhibit no variation in the azimuthal coordinate direction.  Hence, in the expression for <math>~\Phi_B</math>: 
* the azimuthal coordinate, <math>~\phi</math>, need not appear explicitly as an independent variable, so the integral over this angular coordinate immediately gives, <math>~2\pi</math>; 
* the index, <math>~m</math>, must be set to zero, so there is no summation over this index; and,
* drawing from the [[#Toroidal_Functions|additional details provided in Table 5, below]] the single surviving toroidal function is, <math>~Q_{-1 / 2}</math>, which can be rewritten in terms of the complete elliptic integral of the first kind.

As a result of these simplifications, anywhere along the boundary of our cylindrical-coordinate mesh, a valid expression for the ''axisymmetric'' gravitational potential is,

<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^'}  
Q_{- 1 / 2}(\chi) 
</math>
  </td>
</tr>
<tr>
  <td align="right">
&nbsp;
  </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>~\equiv</math>
  </td>
  <td align="left">
<math>~\biggl[ \frac{2}{1+\chi} \biggr]^{1 / 2} =
\biggl[ \frac{4\varpi \varpi^'}{(\varpi + \varpi^')^2 + (z - z^')^2} \biggr]^{1 / 2}
\, .
</math>
  </td>
</tr>
</table>
Actually, this expression for the potential is not only valid along the outer boundaries of the computational mesh, but ''anywhere'' inside or outside of the mass distribution.