I'm solving the following equation
$$
\Delta \phi = 0, \quad (x, y, z) \in \Omega \equiv \mathbb{R}^3/\text{torus}
$$
with the following boundary conditions:
$$
\phi (x,y,z) \to z, \quad x, y, z \gg R
$$
$$
\vec{n} \cdot \vec{\nabla} \phi (x, y,z)|_\text{torus} = 0
$$
where the torus is given by the following parametrization
$$
\begin{aligned}
x &= (c + r \cos \theta) \cos \varphi \\
y &= (c + r \cos \theta) \sin \varphi \\
z &= r \sin \theta
\end{aligned} \quad \quad
\begin{aligned}
0 &\leq r \leq R < c \\
0 &\leq \theta, \varphi \leq 2 \pi
\end{aligned}
$$
So in other words, I'm seeking the solution for some kind of a potential $\phi$ that approaches linearly increasing function $\phi = z$ far from the torus, and so that $\vec{\nabla} \phi$ is parallel to the surface of the torus. The motivation: $\vec{\nabla} \phi$ then gives a vector field that simulates potential inviscid (irrotational) flow around the torus. Alternatively, if the torus is perfectly diamagnetic, then magnetic field outside the torus is given by $\vec{\nabla} \phi$.
Before I went on with the solution, I had a suspicion, that this is not solvable analytically (in a sense that the result will be some known function like sin, cos, exp, …). However, I still think it is expressible by some kind of spherical harmonics, or legendre polynomials, hypergeometric functions…
Secondly, I knew that these coordinates are usually not the ones that are used to solve Laplace equation (the reason might be trouble while attempting to separate variables). So I went on with the majority here and used hyperbolic coordinate chart like this:
$$
\begin{aligned}
x &= \frac{a \sinh \nu}{\cosh \nu – \cos u} \cos \phi \\
y &= \frac{a \sinh \nu}{\cosh \nu – \cos u} \sin \phi \\
z &= \frac{a \sin u}{\cosh \nu – \cos u}
\end{aligned}
$$
Skip to the end if tl;dr.
There is one slight problem here: the meaning of $a$ is not quite the same as $r$ or $c$ before, the same goes for $\nu$. I knew this would be a problem eventually, when I want to force the Neumann condition, so I expressed $\nu$ and $u$ in terms of $r$, $c$, $\theta$:
$$
a^2 = c^2 – r^2
$$
$$
\nu = \log \left( \frac{c}{r} + \sqrt{\frac{c^2}{r^2} – 1} \right)
$$
$$
\cos u = \frac{r + c \cos \theta}{c + r \cos \theta}
$$
Now this assumes that $c > r$, which is not necessarily true outside of the torus (I require $R$ to be less than $c$, that is non-intersecting torus), and it cannot be true very far from the torus. I realized I have a problem imagining what particular limiting values in each coordinate means in cartesian coordinates and how should I enforce each boundary condition. Note that $a$ and $\nu$ become complex when $c < r$. Is this reasonable?
I derived Laplacian in these coordinates, here it is:
$$
\Delta f = \frac{(\cosh \nu – \cos u)^3}{a^2 \sinh \nu} \left[ \frac{\partial}{\partial \nu} \left( \frac{\sinh \nu}{\cosh \nu – \cos u} \frac{\partial f}{\partial \nu} \right) + \right.
$$
$$
\left. + \frac{\partial}{\partial u} \left( \frac{\sinh \nu}{\cosh \nu – \cos u} \frac{\partial f}{\partial u} \right) + \frac{\partial}{\partial \phi} \left( \frac{1}{\sinh \nu (\cosh \nu – \cos u)} \frac{\partial f}{\partial \phi} \right) \right]
$$
Laplace equation $\Delta f = 0$ is not entirely separable, however the following trial solution works
$$
f (\nu, u, \phi) = \sqrt{\cosh \nu – \cos u} N (\nu) U (u) \Phi (\phi)
$$
Plugging this in and doing the standard acrobatics (divide or multiply by whatever common factor will be there with each term and don't forget to divide by the whole trial solution) will yield
$$
\sinh^2 \nu \frac{N^{\prime \prime}}{N} + \sinh \nu \cosh \nu \frac{N^\prime}{N} + \sinh^2 \nu \frac{U^{\prime \prime}}{U} + \frac{\Phi^{\prime \prime}}{\Phi} + \frac{1}{4} \sinh^2 \nu = 0
$$
It is obvious that $\Phi$ and $U$ terms should be negative perfect squares (because of the periodicity in $\phi$ and $u$)
$$
\sinh^2 \nu \frac{N^{\prime \prime}}{N} + \sinh \nu \cosh \nu \frac{N^\prime}{N} – n^2 \sinh^2 \nu – m^2 + \frac{1}{4} \sinh^2 \nu = 0
$$
so that the solutions are $\exp \pm n u$ and $\exp \pm m \phi$. However, let's assume that $m = 0$, because the problem exhibit cylindrical symmetry, so the solution will as well. That leaves us with the equation for $N$
$$
N^{\prime \prime} + \coth \nu N^\prime – \left( n^2 – \frac{1}{4} \right) N = 0
$$
The solution to this are half-integer associated Legendre functions in the argument $\coth \nu$
$$
N(\nu) = \frac{A}{\sqrt{\sinh \nu}} P_{- 1/2}^n (\coth \nu) + \frac{B}{\sqrt{\sinh \nu}} Q_{- 1/2}^n (\coth \nu)
$$
So the general solution with cylindrical symmetry can be written as
$$
f(\nu, u ,\phi) = \sqrt{\frac{\cosh \nu – \cos u}{\sinh \nu}} \times \sum_{n = 0}^\infty \left[ A_n P_{- 1/2}^n (\coth \nu) + B_n Q_{- 1/2}^n (\coth \nu) \right] \left( a_n \cos n u + b_n \sin n u \right)
$$
Note that the expression before the sum is $1/(x^2 + y^2)^{1/4}$.
I played with this expression a bit, for example, I thought that the only term needed is that with $n = 1$. For $n = 1$ I figured that $P$ and $Q$ are not real functions, so I introduced some real combinations instead
$$
f_1 (\nu) = – \frac{i \pi}{2} P_{1/2} (\coth \nu), \quad f_2 (\nu) = \frac{\pi}{2} P_{1/2} (\coth \nu) – i Q_{1/2} (\coth \nu)
$$
It seems that $f_1 (\nu) / \sqrt{\sinh \nu}$ diverges logarithmically for $\nu \to 0$, but $f_2 (\nu) / \sqrt{\sinh \nu}$ remains constant. I think that the limit $\nu \to 0$ has something to do with points far from the torus, but I'm not sure.
At this point, I can't get any further. The main reason is that I'm missing interpretation between cartesian and toroidal coordinates, namely, how do I describe points that are far from the torus? If I could do that, I would enforce the solution $\Phi$ to be equal asymptotically to $z$ in that limit. Also, I can't exactly point to some value of $\nu$ and say "this corresponds to the torus surface", let alone work out the proper equation for Neumann BC.
How can I proceed further with this? Is my effort worth anything? Thank you.
P.S.: please note that symbolic expressions don't bother me, for example solution involving a number described as "root of $f_1 (\nu) – f_2 (\nu)$" is perfectly okay. You know, that separated PDE is better than nothing, quadrature is better than separated PDE and symbolic algebraic expressions are better than quadratures.
Best Answer
Okay, so I figured it out myself. I found the article about the solutions to the Laplace's equation in toroidal coordinates. Although there seems to be no closed-form solution (expressible in finite number of terms), it can be obtained numerically. The key is to rewrite the symmetric (no dependence on $\phi$) solution as follows $$ f(\nu, u) = \sqrt{\cosh \nu - \cos u} \sum_{n = 0}^\infty \left( A_n p_n (\nu) + B_n q_n (\nu) \right) \left( a_n \cos n u + b_n \sin n u \right) \quad \quad \quad \quad (1) $$ where $$p_n (\nu) = P_{n-1/2} (\cosh \nu), \quad q_n (\nu) = Q_{n-1/2} (\cosh \nu) + \frac{i \pi}{2} P_{n-1/2} (\cosh \nu) $$ are the real combinations of half-integer Legendre functions for $\nu \geq 0$.
Now the key idea is to express the function $z \propto \sin u /(\cosh \nu - \cos u)$ in terms of (1), as this is directly the far boundary condition in case $\nu \to 0^+$ (limit $\nu \to 0^+$ corresponds to both the $z$-axis (further parametrized by $u$) and the region far from the torus).
It can be shown, that (probably hard, although I'm no mathematician, so I don't know for sure) \begin{equation} z \propto \frac{\sin u}{\cosh \nu - \cos u} = \sqrt{\cosh \nu - \cos u} \sum_{n = 1}^\infty \frac{\sqrt{32}}{\pi} \, n \, q_n (\nu) \sin n u \end{equation}
Let us now assume, that the true solution (the one satisfying the Neumann BC) is of the form \begin{equation} f (\nu, u) = \frac{\sqrt{32}}{\pi} \sqrt{\cosh \nu - \cos u} \sum_{n = 1}^\infty \Big( - n \, q_n (\nu) + A_n \, p_n (\nu) \Big) \sin n u \end{equation}
Differentiating this w.r.t. $\nu$ at $\nu = \nu_0$ corresponding to the surface of the torus and requiring this derivative to be zero yields the following infinite system of equations \begin{equation} \begin{pmatrix} \alpha_1 & - \beta_2 & 0 & 0 & 0 & \cdots \\ - \beta_1 & \alpha_2 & - \beta_3 & 0 & 0 & \cdots \\ 0 & - \beta_2 & \alpha_3 & - \beta_4 & 0 & \cdots \\ \vdots & \vdots & \ddots & \vdots & \vdots & \ddots \end{pmatrix} \begin{pmatrix} A_1 \\ A_2 \\ A_3 \\ \vdots \end{pmatrix} = \begin{pmatrix} \gamma_1 - \delta_2 \\ \gamma_2 - \delta_3 - \delta_1 \\ \gamma_3 - \delta_4 - \delta_2 \\ \vdots \end{pmatrix} \end{equation} where \begin{equation} \begin{aligned} \alpha_n &= \sinh \nu_0 \; p_n (\nu_0) + 2 \cosh \nu_0 \; p_n^\prime (\nu_0) \\ \beta_n &= p_n^\prime (\nu_0) \\ \gamma_n &= n \, \sinh \nu_0 \; q_n (\nu_0) + 2 n \, \cosh \nu_0 \; q_n^\prime (\nu_0) \\ \delta_n &= n \, q_n^\prime (\nu_0) \end{aligned} \end{equation}
This seems hard, but it turns out that the coefficients $A_n$ falls to zero rather quickly. The practical solution is to introduce some cutoff $N$ so that the system of equations becomes finite \begin{equation} \begin{pmatrix} \alpha_1 & - \beta_2 & 0 & 0 & 0 & \cdots & 0 \\ - \beta_1 & \alpha_2 & - \beta_3 & 0 & 0 & \cdots & 0 \\ 0 & - \beta_2 & \alpha_3 & - \beta_4 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 & 0 & - \beta_{N-1} & \alpha_N \end{pmatrix} \begin{pmatrix} A_1 \\ A_2 \\ A_3 \\ \vdots \\ A_N \end{pmatrix} = \begin{pmatrix} \gamma_1 - \delta_2 \\ \gamma_2 - \delta_3 - \delta_1 \\ \gamma_3 - \delta_4 - \delta_2 \\ \vdots \\ \gamma_N - \delta_{N-1} \end{pmatrix} \end{equation}
and $N$ is chosen so that the magintude of $A_N p_N (\nu_0)$ is below some threshold (this is because $0 < \nu < \nu_0$ where $\nu = 0$ corresponds to $z$ axis and/or region far from the torus and $\nu_0$ corresponds to the surface of the torus, moreover, $p_n (\nu)$ are strictly increasing functions; the value of $N$ therefore serves as a global approximation for a given $\nu_0$). The following table summarizes the maximum needed $N$ for threshold equal to $10^{-6}$ for various values of $\kappa = r/c, \; \; 0 < \kappa < 1$.
\begin{array}{| c | c | c |} \hline \kappa & \nu_0 & N \\ \hline 0.01 & 5.298 & 3 \\ 0.10 & 2.993 & 5 \\ 0.30 & 1.874 & 8 \\ 0.70 & 0.896 & 17 \\ 0.90 & 0.467 & 35 \\ 0.95 & 0.323 & 52 \\ 0.99 & 0.142 & 123 \\ \hline \end{array}
It is obvious that the larger $\kappa$ is, the more (numerical) effort is needed to calculate field to a desired precision.
The relationship between torus' major radius $c$, minor (tube) radius $r$ and parameters of hyperbolic coordinates $a$ and $\nu_0$ are $$ a = c \sqrt{1-\kappa^2}, \quad \nu_0 = \log \left( \frac{1}{\kappa} + \sqrt{\frac{1}{\kappa^2} - 1} \right) $$ where $\kappa = r/c$.
Once we have the solution $f$ we can calculate $\vec{H}$ as $\vec{\nabla} f$.
To verify all of the above, I plotted the solution for $\vec{H}$ for various $\kappa$ (= 0.1, 0.5, 0.9) and $c = 1$ in Mathematica as cuts in planes $x z$ and $x y$. It shows both field streamlines (black solid lines) and field intensity (color blue = zero, white = 1 = far field, red = maximum (different for each $\kappa$)).
$\kappa = 0.1$:
$\kappa = 0.5$:
$\kappa = 0.9$:
This should be correct, as the key features are present: $\vec{H}$ is parallel to the surface of the torus (obvious from the field lines), $\vec{H} \to \vec{e}_z$ far from the torus (obvious from the color) and $\Delta f = 0$ ($f$ is constructed in a way this holds true).
It is interesting that regions where $| \vec{H} | = 0$ (blue color, also vanishing field lines on the surface of the torus) are not exactly at the top and bottom of the torus, instead, they are shifted inward (the more the bigger $\kappa$ is).
We can also look at the magnitude of the field in three interesting points: $A: (0,0,0)$, $B: (1-\kappa, 0, 0)$, $C: (1+\kappa, 0,0)$. This is plotted below:
Blue: A, black: B, orange: C
One sanity check is $\kappa \to 0$, for which $A \to 0$, because the center doesn't "feel" the torus and $B, C \to 2$ - this corresponds to an infinitely long cylinder in perpendicular field. So In the limit $\kappa \to 0$ the system behaves like cylinder, just curled into torus, geometrically speaking (this can be seen from the image above, where in the plane $x z$ the field looks like two parallel cylinders that don't "feel" each other). For higher $\kappa$, this is no longer true. It is reasonable that the field is strongest in the point B, however A and B coincide in the limit $\kappa \to 1$ (as $B \to A$ in that limit). If I'm not mistaken, the field inside the torus diverges in the limit $\kappa \to 1$. This is due to the singularity arising when $r$ approaches $c$ (we call that a 'horn torus') in the center point $(0,0,0)$, where, in theory, the magnetic field is enforced in an infinitely small region resulting in the infinite magnitude.
I hope someone will find this useful.