View source for Computation of Infinite Depth Green's Function

The normal incidence Green's function for infinite depth has already been published [[Squire and Dixon 2000]], [[Dixon and Squire 2001]], the solution being expressed in terms of the auxilary sine and cosine integrals, or equivalently in terms of the exponential integral (see Appendix~36). Similarly, the finite depth Green's function was derived by \can{EP03crack} and was also presented in Section~\ref{sec-fdGF} of this thesis for both normal and oblique incidence. However, at present, the required derivatives <math>g^{(n)}(x)</math> of the infinite depth Green's function for oblique incidence may only be calculated exactly when <math>x=0^\pm</math>. Formulae for these quantities were published by \can{WS02} and are also derived below in Appendix~27. Although the aforementioned derivatives may be calculated when <math>x\neq0</math> by Fast Fourier Transform, as shown in Section~\ref{sec-GFinf}, this is often not the most efficient approach, especially when the Green's function is only required to be evaluated at a few points. The following two sections outline two approaches that are more suitable for such a scenario.

The first uses a partial fractions decomposition to reduce them (the required derivatives of <math>g</math>, or more specifically of <math>g_1</math>) into the sum of five simpler transforms which each satisfy linear inhomogeneous second order ordinary differential equations (ODEs). The inhomogeneous term is the zeroth order modified Bessel function of the second kind, <math>K_0</math>, which has a standard series representation \cite[p375]{AS65}. Using a method that is similar to both the Method of Undetermined Coefficients \cite[p164ff.]{BP97} and the Method of Frobenius \cite[p262ff.]{BP97}, we search for solutions that have similar expansions to that of <math>K_0</math>, and obtain simple recurrence relations for the series coefficients. The series obtained are convergent for all values of <math>x</math>, with less terms needed in the expansion when <math>x</math> is smaller.

A second approach, which uses a contour integral to transform the inversion integral into one that may be calculated more easily by numerical quadrature, is more efficient for larger values of <math>x</math>. It is described below in Section~19, but first we will introduce the series expansion method.

Note that in both methods it is assumed in their derivations that the poles <math>\alpha_n</math> have imaginary part greater than zero. However, it will be seen that the final expressions for <math>g</math> retain the same form in the limits as any poles that need to be are allowed to return to the real line. Consequently, all results hold for any value of <math>\varpi</math>.

==Series Expansion==
As stated in the introduction to this appendix, the series representation of <math>g_1(x)</math> is obtained by breaking it up into the sum of five simpler transforms which each satisfy a different ODE. Section~3 presents this decomposition and derives the ODEs that they satisfy. Section~10 solves these ODEs by proposing series solutions of a particular type, which are similar to the series expansion of the inhomogeneous term (a modified Bessel function of the second kind).

===Ordinary Differential Equations===
From (\ref{g1hat-inf}), the Fourier transform of the <math>2n</math>-th derivative of <math>g_1</math> may be expanded by partial fractions as

<center><math>\begin{matrix}
(-ik)^{(2n)}\hat g_1(k)=-\frac1{\pi \mbox{i}}\sum^2_{m=-2}A_m\gamma_m(\mbox{i}\alpha_m)^{2n+1}\hat g_{1m}(k),
\end{matrix}</math></center>

where the <math>A_n</math> were defined in Section~\ref{sec-GFinf} as <math>A_n=-\gamma_n^3/\alpha_n(4\gamma_n^5+1)</math>, and

<center><math>\begin{matrix}
\hat g_{1n}(k)=\frac{\pi}{\kappa(\alpha_n^2-k^2)}.
\end{matrix}</math></center>

Multiplying this by <math>\alpha_n^2-k^2</math> gives the Fourier transform of the ODE

<center><math>\begin{matrix}
\big(\partial_x^2+\alpha_n^2\big)g_{1n}(x)=K_0(l|x|),
\end{matrix}</math></center>

where <math>K_0</math> is a zeroth order modified Bessel function of the second type [[AS65]]. For <math>|=Arg= (w)|<\pi</math>, it has the expansion

<center><math>\begin{matrix}
K_0(w)&=-\big(\log(w/2)+\gamma_e\big)I_0(w)+\sum_{m=1}^\infty\chi_m\frac{(w/2)^{2m}}{(m!)^2}\\
&=\sum_{m=0}^\infty a_mw^{2m}\big(d_m-\log w\big),(7)
\end{matrix}</math></center>

where <math>\gamma_e</math> is the Euler constant, <math>\chi_m=\sum_{r=1}^m1/r</math>, and <math>I_0</math> is the modified Bessel function of the first kind \citeaffixed{AS65}{also given by}:

<center><math>\begin{matrix}
I_0(w)=\sum_{m=0}^\infty\frac{(w/2)^{2m}}{(m!)^2}=\sum_{m=0}^\infty a_mw^{2m}.
\end{matrix}</math></center>

The coefficients <math>a_m</math> and <math>d_m</math> are given by <math>a_m=1/2^{2m}(m!)^2</math> and if <math>\chi_0=0</math>, <math>d_m=\log2-\gamma_e+\chi_m</math>, so the above series converges for all <math>w</math>. Thus we anticipate that series representations for the <math>g_{1n}</math> that satisfy (6) will also converge for all <math>x</math>. The <math>a_n</math> and <math>d_n</math> themselves may be calculated most efficiently by the recurrence relations <math>a_n=a_{m-1}/4m^2</math> and <math>d_m=d_{m-1}+1/m</math>, with <math>a_0=1</math> and <math>d_0=\log2-\gamma_e.</math>\\

The <math>g_{1n}</math> should also satisfy the initial conditions

<center><math>\begin{matrix}
g_{1n}(0)-q_n=g'_{1n}(0)=0,
\end{matrix}</math></center>

where the <math>q_n</math> can be evaluated analytically as shown in Section~27 below.

In addition, we know that even derivatives of <math>g_1</math> are even functions, and that odd derivatives are odd functions (since their Fourier transforms are). Consequently values of <math>g_1^{(n)}</math> for negative <math>x</math> can be generated in an obvious fashion from their values at <math>-x</math>, and so we will assume <math>x>0</math> from here on.

===Trial Series Solution===
We will solve equations (6) by using a variation of the Method of Frobenius to obtain individual series expansions for each of the <math>g_{1n}</math>. Letting <math>w=lx</math>, <math>\alpha=\alpha_n/l</math>, and <math>J(w)=l^2\big(g_{1n}(x)-q_n\cos\alpha_nx\big)</math>, <math>J</math> should now satisfy

<center><math>\begin{matrix}
J''(w)+\alpha^2J(w)=K_0(w)\quad=for <math>w>0</math>= ,(12)\\
J(0^+)=J'(0^+)=0.(13)
\end{matrix}</math></center>

We can see from their transforms that each of the <math>g_{1n}</math> are continuous, even, and have continuous first derivatives, but that they have logarithmic singularities in their second derivatives, and so we propose a trial solution for (11) of the form

<center><math>
J(w)=\sum_{m=0}^\infty\big(b_m+c_m\log w\big)w^{2m+2},
</math></center>

which clearly satisfies (13). This has first derivative

<center><math>\begin{matrix}
J'(w)=\sum_{m=0}^\infty\big((2m+2)b_m+c_m+(2m+2)c_m\log w\big)w^{2m+1},
\end{matrix}</math></center>

enabling odd derivatives of <math>g_1</math> to also be produced once the previous even one has been calculated. Proceeding to find <math>J</math>, however, we now substitute our trial solution into the left hand side of (12), giving

<center><math>\begin{matrix}
J''(w)+\alpha^2J&=2b_0+3c_0+2c_0\log w\\
&\quad+\sum_{m=1}^\infty\big( (2m+2)(2m+1)b_{m}+\alpha^2b_{m-1}+(4m+3)c_m\big)w^{2m}\\
&\quad+\sum_{m=1}^\infty\big((2m+2)(2m+1)c_m+\alpha^2c_{m-1}\big)w^{2m}\log w.
\end{matrix}</math></center>

Consequently, <math>J</math> will satisfy (11) if our <math>b_m</math> and <math>c_m</math> satisfy the recurrence relations

<center><math>\begin{matrix}
(2m+2)(2m+1)b_{m}+\alpha^2b_{m-1}+(4m+3)c_{m}=a_{m}d_{m}
\quad&=for <math>m\geq1</math>= ,(17)\\
(2m+2)(2m+1)c_m+\alpha^2c_{m-1}=-a_m\quad&=for <math>m\geq1</math>= ,(18)
\end{matrix}</math></center>

and if <math>2c_0=-a_0</math> and <math>2b_0=a_0d_0-3c_0</math> (i.e., if <math>b_0=(\log2+3/2-\gamma_e)/2</math> and <math>c_0=-1/2</math>). The latter relation (18) implies that <math>|c_m/c_{m-1}|\to0</math> as <math>m\to\infty</math>, and so the series involving those coefficients in (14) will converge for all <math>w</math>. This also implies that the <math>c_m</math> themselves tend towards zero as <math>m\to\infty</math>, and so we can also deduce from (17) that the ratio <math>|b_m/b_{m-1}|</math> will do the same for large <math>m</math>. Consequently, (14) will converge in entirety for all <math>w</math>, although it will have a branch cut along the negative real axis due to the logarithmic factor.

==Contour Integration==
An alternative means of calculating <math>g</math> is by using contour integration to transform the inverse Fourier transform into an integral which is more easily evaluated numerically. This technique is more efficient for larger values of <math>x</math>, but it turns out that the resulting integral may actually be evaluated analytically when <math>x=0</math>. Section~20 explains the contour used and outlines the numerical method of approximation of the resulting integral for general <math>x</math> values, while Section~27 derives the formula for that integral when <math>x=0</math>.

Some additional simplifications can be obtained for normal incidence. We show that the integral for general <math>x</math> reduces to an expression involving the well-known and easily evaluated exponential integral function. This also acts as a further check on the oblique incidence results because the normal incidence results can be compared to those obtained by \can{DS01}. This is done in Section~36.

===Derivation and Numerical Approximation===
For larger values of <math>x</math>, when the series expansions for the <math>g_{1n}</math> take longer to converge, we would like to transform the integral

<center><math>
g_{1n}(x)=\frac12\int_\minf^\infty\frac{\mbox{e}^{\mbox{i} kx}\mbox{d} k}{\kappa(\alpha_n^2-k^2)}
</math></center>

into one that is easier to evaluate numerically. To do this we close the contour in the upper half plane, allowing for the pole at <math>k=\alpha_n</math> and the branch cut along the line <math>k=\mbox{i}\times[l,\infty)</math>. Thus we can write

<center><math>
g_{1n}(x)=g_{3n}(x)-\frac{\epsilon_n\pi \mbox{i}}{2\alpha_n\gamma_n}\mbox{e}^{\mbox{i}\alpha_nx},
</math></center>

where <math>\epsilon_n=\mbox{sgn}\big(\mathfrak{Re}[\gamma_n]\big)</math>, and

<center><math>
g_{3n}(x)=\int_{\mbox{i} l}^{\mbox{i}\infty}\frac{\mbox{e}^{\mbox{i} kx}\mbox{d} k}{\kappa(\alpha_n^2-k^2)}
=\int_{l}^{\infty}\frac{\mbox{e}^{-\kappa x}\mbox{d}\kappa}{k(\kappa^2+\alpha_n^2)}.
</math></center>

The <math>\epsilon_n</math> results from defining the square root in <math>\kappa</math> as having a branch cut along the negative real axis. Consequently, <math>\sqrt{a_n^2+l^2}=\epsilon_n\gamma_n</math> in the residue of <math>g_{1n}</math> at <math>\alpha_n</math>. In addition, <math>\kappa</math> in the left hand integral in (22) is taken from the positive imaginary axis.

Now let us define the analogous sum to (4):

<center><math>\begin{matrix}
g_3(x)=-\frac1\pi\sum_{n=-2}^2A_n\alpha_n\gamma_ng_{3n}(x)
=-\frac1\pi\int_0^\infty\frac{k\mbox{e}^{-\kappa x}\mbox{d}\kappa}{\Lambda^2(k)k^2+1}.
\end{matrix}</math></center>

From (21), we now can write the relationship between <math>g_1</math> and <math>g_3</math>:

<center><math>\begin{matrix}
g_1(x)=g_3(x)+\frac \mbox{i}2\sum_{n=-2}^2A_n\epsilon_n\mbox{e}^{\mbox{i}\alpha_nx}.
\end{matrix}</math></center>

Adding <math>g_0</math> to <math>g_1</math> then enables us to write <math>g(x)=g_2(x)+g_3(x)</math>, where

<center><math>\begin{matrix}
g_2(x)=\mbox{i}\sum_{n=-2}^0A_n\mbox{e}^{\mbox{i}\alpha_nx}.
\end{matrix}</math></center>

This confirms the results (\ref{g2a}) and (\ref{g3a}), which were obtained by taking the limit in (\ref{g}) as <math>H\to\infty</math> directly.

<math>g_3(x)</math> as given by (23) is an integral that is more amenable to calculation by numerical quadrature than the inverse Fourier transform <math>g_1</math>. The integrand will never have any singularities on the real line and, because it is also <math>O(\mbox{e}^{-kx})</math> as <math>k\to\infty</math>, <math>g_3</math> will be infinitely differentiable for <math>x>0</math>. However, since the non-exponential factor is only <math>O(k^{-9})</math> as <math>k\to\infty</math>, its eighth derivative will be singular at the origin, as was mentioned in Section~\ref{sec-gsing}.

When <math>\theta=0</math>, <math>g_3</math> may be calculated analytically in terms of the exponential integral <math>E_1</math> as described below in Section~36. First, however, we will show how any derivative of <math>g_3</math> may be efficiently calculated for oblique incidence.

Differentiating (23) <math>n</math> times gives

<center><math>\begin{matrix}
g_3^{(n)}(x)&=-\frac{(-1)^{n}}\pi\int_l^\infty\frac{k\kappa^ne^{-\kappa x}\mbox{d}\kappa}{\Lambda^2(k)k^2+1}(25)\\
&=-\frac{(-1)^{n}}\pi\times x^{8-n}\mbox{e}^{-lx}\int_0^\infty
\frac{t^{1/2}(t+lx)^n(t+2lx)^{1/2}\mbox{e}^{-t}\mbox{d} t}{t(t+2lx)\big(t^2(t+2lx)^2+\varpi x^4\big)^2+x^{10}}(26),
\end{matrix}</math></center>

where <math>t=(\kappa-l)x</math>. This integral may now be found using generalised Gauss-Laguerre integration, using the weight function <math>t^{1/2}\mbox{e}^{-t}</math>. The remainder of the integrand will now be relatively smooth, with the exception of the <math>\sqrt{t+2lx}</math> term which has a branch cut on the line <math>(-\infty,-2lx]</math>. However, we don't need to integrate over that singularity and the larger <math>x</math> gets the further it moves away from the region of integration. Hence, the integral will converge faster for larger values of <math>x</math>. Also note that (26) shows that <math>g_3^{(n)}</math> decays exponentially as <math>x\to\infty</math>.

===Exact Calculation of Derivatives at the Origin===
Equation (22) when <math>x=0</math> is
<center><math>\begin{matrix}(28)
q_n'=q_n+\frac{\mbox{i}\epsilon_n}{2\alpha_n\gamma_n}=g_{3n}(0^+)=I(\alpha_n)+I(-\alpha_n),
\end{matrix}</math></center>
where
<center><math>\begin{matrix}(29)
I(\alpha)&=\frac1{2\mbox{i}\alpha}\int_l^\infty\frac{\mbox{d}\kappa}{k(\kappa-\mbox{i}\alpha)}
=-\frac1{\alpha}\int_0^1\frac{\mbox{d} w}{(\alpha-\mbox{i} l)w^2-(\alpha+\mbox{i} l)},
\end{matrix}</math></center>
and we have made the transformation <math>\kappa=l(1+w^2)/(1-w^2)</math>. If we also define
<math>\b</math> and <math>\gamma</math> by the relations <math>\b^2=(\alpha+\mbox{i} l)/(\alpha-\mbox{i} l)</math> and <math>\b=\gamma/(\alpha-\mbox{i} l)=(\alpha+\mbox{i} l)/\gamma</math> (so that <math>\gamma^2=\alpha^2+l^2</math>), (29) becomes
<center><math>\begin{matrix}(30)
I(\alpha)&=-\frac1{\alpha(\alpha-\mbox{i} l)}\int_0^1\frac{\mbox{d} w}{w^2-\b^2}
=-\frac1{2\alpha\gamma}\left(\int_0^1\frac{\mbox{d} w}{w-\b}
-\int_0^1\frac{\mbox{d} w}{w+\b}\mbox{i}ght)\\
&=\frac1{2\alpha\gamma}\big(\log(1+\b)-\log(1-\b)-\log\b+\log(-\b)\big).
\end{matrix}</math></center>
To simplify (30), we choose <math>\b</math> to be the square root from the upper half plane (which also implies that <math>1+\b</math> is in the upper half plane, and that <math>\log(-\b)=\log\b-\mbox{i}\pi</math>). <math>1-\b</math> will then be in the lower half-plane so the quotient <math>(1+\b)/(1-\b)</math> corresponds to an anti-clockwise rotation of <math>1+\b</math>. Potentially this could rotate <math>1+\b</math> across the branch cut in the logarithm on the negative real axis, but it can be shown explicitly that with our choice of <math>\b</math> the aforementioned quotient remains in the upper half plane. Thus, the normal logarithm rules for quotients hold, and so (30) may be written
<center><math>\begin{matrix}(31)
I(\alpha)&=\frac1{2\alpha\gamma}\big(\log\frac{1+\b}{1-\b}-\pi \mbox{i}\big)
=\frac1{2\alpha\gamma}\big(\log\frac{\gamma+\alpha+\mbox{i} l}{\gamma-\alpha-\mbox{i} l}-\pi \mbox{i}\big).
\end{matrix}</math></center>
Note that the above expression is even in <math>\gamma</math>, so it doesn't matter after all which root we choose for it. Thus,
<center><math>\begin{matrix}(32)
q_n'=\frac1{2\alpha_n\gamma_n}\big(\log\frac{\gamma_n+\alpha_n+\mbox{i} l}{\gamma_n-\alpha_n-\mbox{i} l}
-\log\frac{\gamma_n-\alpha_n+\mbox{i} l}{\gamma_n+\alpha_n-\mbox{i} l}\big)
=\frac1{2\alpha_n\gamma_n}\log\Theta_n,
\end{matrix}</math></center>
where <math>\Theta_n=(\gamma_n+\alpha_n)/(\gamma_n-\alpha_n)</math>. The latter step is justifed since if <math>\im[\alpha_n]>0</math> the argument of the first logarithm will be in the first quadrant, while the argument of the second logarithm will be in the second quadrant. This means that we can now write
<center><math>(33)
g_3^{(2n)}(0^+)=-\frac1\pi\sum_{n=-2}^2A_n\alpha_n\gamma_nq'_n
=-\frac1{2\pi}\sum_{m=-2}^2A_m(\mbox{i}\alpha_m)^{2n}\log{\Theta_m}.
</math></center>
Now, for small <math>l</math>, <math>\alpha_n\sim\epsilon_n\gamma_n\big(1-l^2/2\gamma_n+O(l^4)\big)</math>, where <math>\epsilon_n=\mbox{sgn}(\im[\gamma_n])</math>, so at first glance it seems like the above expression will become singular as <math>l\to0</math>. However, this is not the case. If we define <math>A'_n=A_n\alpha_n/\gamma_n=-\gamma_n^2/(4\gamma_n^5+1)</math>, <math>A_n\sim\epsilon_nA_n'+O(l^2)</math> and it is easily shown that for <math>n=0,1,2,3</math> they satisfy the rules
<center><math>(34)
\sum_{m=-2}^{2}A_m'\gamma_m^{2n}=0\quad=for <math>n=0,\ldots,3</math>,= 
</math></center>
which are similar to those presented by \can{SD00}. Thus as <math>l\to0</math>,
<center><math>\begin{matrix}(35)
g_3^{(2n)}(0^+)&\sim-\frac1{2\pi}\sum_{m=-2}^2A_m'\big(1+O(l^2)\big)\times
\left(\log\big( \gamma_n^2\big(1+O(l^2)\big) \big)-\log\frac{l^2}4\mbox{i}ght)\\
&\to-\frac1{2\pi}\sum_{m=-2}^2A_m'(\mbox{i}\gamma_m)^{2n}\log{\gamma^2_m}.
\end{matrix}</math></center>

===Normal Incidence Simplification===
An alternative expansion to (28) for the <math>g_{3n}</math> is

<center><math>\begin{matrix}
2\gamma_n^2g_{3n}(x)&=\int_l^\infty\left(\frac2k-\frac1{k+\mbox{i}\gamma_n}-\frac1{k-\mbox{i}\gamma_n}\mbox{i}ght)\mbox{e}^{-\kappa x}\mbox{d}\kappa\\
&=2K_0(lx)-\int_l^\infty\left(\frac1{k+\mbox{i}\gamma_n}+\frac1{k-\mbox{i}\gamma_n}\mbox{i}ght)\mbox{e}^{-\kappa x}\mbox{d}\kappa.
\end{matrix}</math></center>

Using the rules (34) to eliminate the Bessel function <math>K_0</math>, we can thus write the <math>2n</math>-th derivative of <math>g_3</math> as

<center><math>\begin{matrix}
g_3^{(2n)}(x)&=-\frac1\pi\sum_{m=-2}^2A_m'(\mbox{i}\alpha_m)^{2n}g_{3m}(x)\\
&=\frac1{2\pi}\sum_{m=-2}^2A_m'(\mbox{i}\alpha_m)^{2n}\int_l^\infty\left(\frac1{k+\mbox{i}\gamma_n}+\frac1{k-\mbox{i}\gamma_n}\mbox{i}ght)\mbox{e}^{-\kappa x}\mbox{d}\kappa.
\end{matrix}</math></center>

Now, if we let <math>\theta\to0</math>, <math>l\to0</math>, <math>\alpha_n^2\to\gamma_n^2</math>, and <math>\kappa\to k</math> on the real line, (38) becomes

<center><math>\begin{matrix}
g_3^{(2n)}(x)&=\frac1{2\pi}\sum_{m=-2}^2A_m'(\mbox{i}\gamma_m)^{2n}\left(
\mbox{e}^{\mbox{i}\gamma_mx}E_1(\mbox{i}\gamma_mx)+\mbox{e}^{-\mbox{i}\gamma_mx}E_1(-\mbox{i}\gamma_mx)\mbox{i}ght),
\end{matrix}</math></center>

where <math>E_1</math> is the well-known exponential integral function \cite[p228ff.]{AS65} defined as follows

<center><math>
E_1(w)=\int_w^\infty \frac{\mbox{e}^{-t}}t \mbox{d} t
=-\gamma-\log w-\sum_{n=1}^\infty\frac{(-w)^n}{n\,n!}
\quad=for <math>|\text{Arg= (w)|<\pi</math>}.
</math></center>

Note that, since no imaginary solutions exist for the infinite depth dispersion relation, <math>\big|=Arg= [\pm \mbox{i}\gamma_nx]\big|<\pi</math>, so the above formulae are valid. Also note that the <math>g_3^{2n}</math> will be everywhere bounded, despite the logarithmic singularities in <math>E_1</math>. This also follows from the rules (34), allowing us to deduce that taking the limit in (39) as <math>x\to0^+</math> gives

<center><math>\begin{matrix}
g_3^{(2n)}(0^+)&=-\frac1{2\pi}\sum_{m=-2}^2A_m'(\mbox{i}\gamma_m)^{2n}
\big(\log(\mbox{i}\gamma_m)+\log(-\mbox{i}\gamma_m)\big)\\
&=-\frac1\pi\sum_{m=-2}^2A_m'(\mbox{i}\gamma_m)^{2n}\log\gamma_m
-2i\big(A'_{-2}(\mbox{i}\gamma_1)^{2n}+A'_2(\mbox{i}\gamma_2)^{2n}\big)(41)\\
&=-\frac1{2\pi}\sum_{m=-2}^2A_m'(\mbox{i}\gamma_m)^{2n}\log\gamma_m^2,(42)
\end{matrix}</math></center>

which agrees with (35). Recalling that <math>A_n\to\mbox{sgn}(\mathfrak{Im}[\gamma_n])A_n'</math> as <math>l\to0</math>, adding <math>g_2^{(2n)}(0^+)</math> to (41) and using (34) to simplify the product, then gives the following result for <math>g^{(2n)}(0^+)</math> itself:

<center><math>\begin{matrix}
g^{(2n)}(0^+)&=-\frac1\pi\sum_{m=-2}^2A_m'(\mbox{i}\gamma_m)^{2n}\log\gamma_m
-2i\big(A'_{-2}(\mbox{i}\gamma_1)^{2n}+A'_2(\mbox{i}\gamma_2)^{2n}\big).
\end{matrix}</math></center>

The above result is unstable as <math>\gamma_2</math> becomes negative and so is only reliable for <math>\varpi>\varpi_\infty</math>. For a result that holds for all <math>\varpi</math>, use (42) and the ordinary definition of <math>g_2</math>.

Odd derivatives of <math>g_3</math> may also be obtained by differentiating (39) directly, again making use of (40) and also (34) to eliminate any <math>1/x</math> type singularities produced. Thus,

<center><math>\begin{matrix}
g_3^{(2n+1)}(x)&=\frac1{2\pi}\sum_{m=-2}^2A_m'(\mbox{i}\gamma_m)^{2n+1}\left(
\mbox{e}^{\mbox{i}\gamma_mx}E_1(\mbox{i}\gamma_mx)-\mbox{e}^{-\mbox{i}\gamma_mx}E_1(-\mbox{i}\gamma_mx)\mbox{i}ght),
\end{matrix}</math></center>

and so we can write the general <math>n^\mbox{th}</math> (<math>0\leq n\leq7</math>) derivative of <math>g_3</math> as

<center><math>\begin{matrix}
g_3^{(n)}(x)&=\frac1{2\pi}\sum_{m=-2}^2A_m'(\mbox{i}\gamma_m)^{n}\left(
\mbox{e}^{\mbox{i}\gamma_mx}E_1(\mbox{i}\gamma_mx)+\mbox{e}^{\mbox{i}(n\pi-\gamma_mx)}E_1(-\mbox{i}\gamma_mx)\mbox{i}ght).
\end{matrix}</math></center>

Recall that from (\ref{god}) that we already know what values <math>g^{(2n+1)}</math> will take at the origin so there is no need to take the limit in (44) as <math>x\to0^+</math>. In addition, we can also check the above results numerically by letting <math>l\to0</math> and <math>\kappa\to k</math> in (25) and using generalised Gauss-Laguerre integration (with weight function <math>k^{n+1}\mbox{e}^{-kx}</math>).

A further confirmation of (45) can be obtained by checking the <math>n=0</math> result against the equivalent result of \can{DS01}. Let us denote their nondimensional length variables as <math>x',z',\xi'</math> and <math>\zeta'</math>, which equate to <math>x/a,z/a,\xi/a</math> and <math>\zeta/a</math> respectively. In their expression (6) for the Green's function in that scheme, which we shall denote <math>G_\mbox{SD}</math>, the integral when <math>x'>\xi'</math> and <math>z',\zeta'\approx0</math> is obtained by doubling the real part of equation (14b). Doing so, and then differentiating with respect to <math>\zeta'</math> followed by <math>z'</math> gives the result

<center><math>\begin{matrix}
g(x-\xi)&=\frac1{a^2}\partial^2_{z'\zeta'}G_\mbox{SD}(x'-\xi',0,0)+\mbox{i} A_0\cos\gamma_0x\\
&=\mbox{i}\sum_{m=-2}^0A'_m\epsilon_me^{\mbox{i}\epsilon_m\gamma_mx}
+\frac1\pi\sum_{m=-2}^2A_m'g_\mbox{cos}(\epsilon_m'\gamma_mx),
\end{matrix}</math></center>

where <math>\epsilon_n'=\sgn\big(\mbox{e}[\gamma_n]\big)</math>. The cosine term needs to be added to allow for the difference between taking the Cauchy Principal value of the integral (\ref{G1}) and letting the poles <math>\pm\alpha_0</math> approach the real line from the upper or lower half plane as we did, while <math>g_\mbox{cos}</math> is the auxiliary cosine function [[AS65]] which may be written as

<center><math>
g_\mbox{cos}(w)=\int_0^\infty\frac{\cos t}{w+t}\mbox{d} t
=\frac12\big(\mbox{e}^{\mbox{i} w}E_1(\mbox{i} w)+\mbox{e}^{-\mbox{i} w}E_1(-\mbox{i} w)\big),\quad=<math>\mbox{e}[w]>0</math>= .
</math></center>

Since the first summation in (46) is obviously <math>g_2(x)</math>, the second sum must be <math>g_3(x)</math>. In light of the definition of <math>g_\mbox{cos}</math>, the result presented by \can{DS01} clearly agrees with (45).


[[Reflections on Ice]]