One
To explain how allium works, let’s start (as I did), with a sphere
with constant sound speed \(c\). The general
solutions of the Helmholtz equation are the spherical Bessel
functions \(j_\ell(\omega r/c)\) and \(y_\ell(\omega
r/c)\), where \(\ell\) is the angular degree. The solutions that satisfy
the inner boundary conditions are functions of the first kind. The
outer boundary condition is only satisfied at particular values of the
frequency \(\omega\). These values are eigenvalues or
eigenfrequencies and the solutions at these frequencies are
eigenfunctions.
Because the Helmholtz equation is linear, the amplitude of the
solution is arbitrary, so I (and allium) choose a solution with
amplitude 1. i.e., \(f(r)=j_\ell(\omega r/c)\).
Two
For a given \(\ell\), consider the case of two sound speeds, i.e.
\[\begin{split}c(r) = \left\{ \begin{array}{lcr}
c_1 & \mathrm{if} & 0 < r < R_1 \\
c_2 & \mathrm{if} & R_1 < r < R_2 \\
\end{array} \right.\end{split}\]
Because the centre must still be regular, the inner solution is still
just \(f_1(r)=j_\ell(\omega r/c_1)\). Let’s call the outer solution
\[f_2(r) = \alpha_2 j_\ell(\omega r/c_2) + \beta_2 y_\ell(\omega r/c_2)\]
You can think of the inner solution similarly, with \(\alpha_1=1\)
and \(\beta_1=0\).
The solution and its first derivative must be continuous at the
interface \(R_1\), i.e.
\[\begin{split}f_1(R_1) = f_2(R_1) & \Rightarrow
j_\ell(\omega R_1/c_1) &= \alpha_2 j_\ell(\omega R_1/c_2) + \beta_2 y_\ell(\omega R_1/c_2) \\
f'_1(R_1) = f'_2(R_1) & \Rightarrow
j'_\ell(\omega R_1/c_1) &= \alpha_2 j'_\ell(\omega R_1/c_2) + \beta_2 y'_\ell(\omega R_1/c_2) \\\end{split}\]
which we can rewrite as a linear equation
\[\begin{split}\left(\begin{array}{cc}
j_\ell(\omega R_1/c_2) & y_\ell(\omega R_1/c_2) \\
j'_\ell(\omega R_1/c_2) & y'_\ell(\omega R_1/c_2)
\end{array}\right)
\left(\begin{array}{c}
\alpha_2 \\
\beta_2
\end{array}\right)=
\left(\begin{array}{c}
j_\ell(\omega R_1/c_1) \\
j'_\ell(\omega R_1/c_1)
\end{array}\right)\end{split}\]
But this only gives us coefficients for some wavefunction at a
frequency \(\omega\). To be an eigenfunction, the solution must
also satisfy the outer boundary condition. To find such a solution,
allium uses \(f(R_2)\) as a residual for which to find a root
(between user-specified bounds), which then returns an eigenfrequency.
Many
Now consider the case of an arbitrary number of layers \(n\). i.e.,
\[\begin{split}c(r) = \left\{ \begin{array}{lcr}
c_1 & \mathrm{if} & 0 < r < R_1 \\
c_2 & \mathrm{if} & R_1 < r < R_2 \\
\ldots \\
c_n & \mathrm{if} & R_{n-1} < r < R_n \\
\end{array} \right.\end{split}\]
In each extra layer, we have a solution
\[f_i(r)=\alpha_i j_\ell(\omega r/c_i) + \beta_i y_\ell(\omega r/c_i)\]
To solve for the extra coefficients, we use additional constraints
from the continuity of \(f\) and \(f'\) at each \(R_i\),
which give
\[\begin{split}\alpha_i j_\ell(\omega R_{i-1}/c_i)+\beta_i y_\ell(\omega R_{i-1}/c_i)
&=\alpha_{i-1}j_\ell(\omega R_{i-1}/c_{i-1})+\beta_{i-1}y_\ell(\omega R_{i-1}/c_{i-1}) \\
\alpha_i j'_\ell(\omega R_{i-1}/c_i)+\beta_i y'_\ell(\omega R_{i-1}/c_i)
&=\alpha_{i-1}j'_\ell(\omega R_{i-1}/c_{i-1})+\beta'_{i-1}y_\ell(\omega R_{i-1}/c_{i-1})\end{split}\]
for \(2\leq i\leq n\). This can again be turned into a linear
equation of the form \(Ax=b\), where, for brevity, we denote
\[\begin{split}\mathcal{J}_{i,j} &= j_\ell(\omega R_i/c_j) \\
\mathcal{J}'_{i,j} &= j'_\ell(\omega R_i/c_j) \\
\mathcal{Y}_{i,j} &= y_\ell(\omega R_i/c_j) \\
\mathcal{Y}'_{i,j} &= y'_\ell(\omega R_i/c_j) \\\end{split}\]
in which case the terms in the matrix equation are
\[\begin{split}x &= (\alpha_2, \beta_2, \alpha_3, \beta_3, \ldots, \alpha_n, \beta_n) \\
b &= (\mathcal{J}_{1,1}, \mathcal{J}'_{1,1}, 0, 0, \ldots, 0, 0) \\
A &= \left(\begin{array}{ccccccccccc}
\mathcal{J}_{1,2} & \mathcal{Y}_{1,2} \\
\mathcal{J}'_{1,2} & \mathcal{Y}'_{1,2} \\
-\mathcal{J}_{2,2} & -\mathcal{Y}_{2,2} & \mathcal{J}_{2,3} & \mathcal{Y}_{2,3} \\
-\mathcal{J}'_{2,2} & -\mathcal{Y}'_{2,2} & \mathcal{J}'_{2,3} & \mathcal{Y}'_{2,3} \\
& & -\mathcal{J}_{3,3} & -\mathcal{Y}_{3,3} & \mathcal{J}_{3,4} & \mathcal{Y}_{3,4} \\
& & -\mathcal{J}'_{3,3} & -\mathcal{Y}'_{3,3} & \mathcal{J}'_{3,4} & \mathcal{Y}'_{3,4} \\
& & & & & & \ddots \\
& & & & & & & -\mathcal{J}_{n-1,n-1} & -\mathcal{Y}_{n-1,n-1} & \mathcal{J}_{n-1,n} & \mathcal{Y}_{n-1,n} \\
& & & & & & & -\mathcal{J}'_{n-1,n-1} & -\mathcal{Y}'_{n-1,n-1} & \mathcal{J}'_{n-1,n} & \mathcal{Y}'_{n-1,n} \\
\end{array}\right)\end{split}\]
This is the matrix problem that Sphere.solve_coefficients()
solves.