1. Introduction
This paper began with a quest for solvable analytical models for fire whirls (Figure 1). The phenomenology is complicated, combining fluid flow and combustion; for an extensive review, see [Reference Tohidi, Gollner and Xiao18]. Mathematical models can be built, but, inevitably, they are complicated. For early models, see [Reference Byram and Martin5, Reference Morton12]. For a much more sophisticated model, including compressibility, viscosity and unsteadiness, see [Reference Forbes and Walters7].
We start from a scalar nonlinear partial differential equation (PDE) for a stream function
$\psi (r,z)$
, as derived by Battaglia et al. [Reference Battaglia, Rehm and Baum3]; here,
$r$
and
$z$
are cylindrical polar coordinates. After linearisation, we obtain a second-order elliptic PDE containing three (dimensionless) parameters,
$A$
,
$B$
and
$C$
;
$A$
is related to the circulation around the whirl,
$B$
is related to the head (or Bernoulli’s constant), and
$C$
is related to temperature changes. In the absence of thermal effects,
$C=0$
and then we recover models used for steady axisymmetric flows of an incompressible inviscid fluid [Reference Batchelor2, Section 7.5], [Reference Saffman16, Section 3.13]. In what follows, we only consider the PDE with
$C\neq 0$
; although the PDE obtained is linear and homogeneous, we do not have any non-trivial solutions – the method of separation of variables does not succeed.
Image of a typical firewhirl. The image is courtesy of U.S. Department of Fish and Wildlife Services.
In the absence of exact solutions, we construct Frobenius-type expansions, both near the ground at
$z=0$
(Section 4.1) and near the axis at
$r=0$
(Section 4.2). These may be regarded as local expansions, but they have limited use if we want to apply boundary conditions at
$z=0$
and
$z=\infty$
, for example. Therefore, we are interested in asymptotic approximations. Motivated by classic approximations for high-frequency diffraction problems, we seek approximations based on the assumption that the parameter
$A$
is large (Section 5.1). As with diffraction theory (and long-distance sound propagation through an inhomogeneous ocean [Reference Brekhovskikh and Godin4, Section 7.6]), the result is a parabolic PDE (see (5.3) below), but it is not separable; as such, it is similar to a PDE known as Popov’s equation [Reference Popov15]. We develop two approaches for solving this equation. In the first, we assume that the parameter
$B$
is also large (Section 5.2). In the second, we solve the PDE using Laplace transforms with respect to a radial coordinate (Section 5.3), inspired by Ockendon and Tew’s analysis of Popov’s equation using Fourier transforms [Reference Ockendon and Tew14]. In both approaches, we have to solve a first-order PDE using standard methods (as outlined in Section B). The large-
$B$
approximation is fairly simple (see (5.7) below), but the more general solution (which does not require
$B$
to be large) is much more complicated, requiring the inverse Laplace transform of an expression involving rational combinations of Airy functions (see (5.10) below). Evidently, further work is needed before this solution can be used and tested.
2. The basic nonlinear model
Batchelor [Reference Batchelor2, Section 7.5] and Saffman [Reference Saffman16, Section 3.13] give clear discussions of steady axisymmetric flow of an incompressible inviscid fluid. The continuity equation is satisfied using a stream function
$\psi$
, with the velocity components given by
\begin{equation} u=v_r=-\frac {1}{r}\,\frac {\partial \psi }{\partial z} \quad \mbox{and}\quad w=v_z=\frac {1}{r}\,\frac {\partial \psi }{\partial r}, \end{equation}
where
$r$
and
$z$
are cylindrical polar coordinates. Then, Batchelor [Reference Batchelor2, equation (7.5.11)] and Saffman [Reference Saffman16, equation (3.13.9)] show that
$\psi (r,z)$
satisfies the nonlinear PDE
\begin{equation} \frac {\partial ^2 \psi }{\partial r^2}-\frac {1}{r}\,\frac {\partial \psi }{\partial r}+\frac {\partial ^2 \psi }{\partial z^2}+F(\psi )=0 \end{equation}
with
\begin{equation} F(\psi )=\Gamma \frac {\textrm {d} \Gamma }{\textrm {d} \psi } - r^2\frac {\textrm {d} H}{\textrm {d} \psi }. \end{equation}
Here,
$\Gamma$
and
$H$
are given functions of
$\psi$
. In the context of magnetohydrodynamics, (2.2) is known as the Grad–Shafranov equation [Reference Landau and Lifshitz10, equation (68.20)], [Reference Atanasiu, Günter, Lackner and Miron1]. In [Reference Llewellyn Smith and Hattori11, equation (12)], the authors have a third term in (2.2b
), proportional to
$r^4$
.
To introduce buoyancy-driven convection, suppose that the
$z$
-axis points upward from the ground at
$z=0$
. Then (2.2b
) has a third term,
$zr^2\,\textrm {d}\Theta /\textrm {d}\psi$
, where
$\Theta =(T-T_0)/T_0$
,
$T$
is the temperature, and
$T_0$
is the ambient temperature [Reference Gledzer and Makarov8, equation (1.4)], [Reference Battaglia, Rehm and Baum3, equation (10)]. The resulting PDE for
$\psi$
will be of interest to us.
3. The basic linear model
To make analytical progress, we assume that
$F(\psi )$
is linear in
$\psi$
and take
\begin{equation*} F(\psi )=( A^2 -B^2r^2+Czr^2)\psi , \end{equation*}
where
$A$
,
$B$
and
$C$
are constants. Dimensionally, given a length scale
$L$
,
$A$
,
$B$
and
$C$
have dimensions
$L^{-1}$
,
$L^{-2}$
and
$L^{-5}$
, respectively. Therefore, to make (2.2a
) dimensionless (ignoring the dimensions of
$\psi$
), we send
$r/L\to r$
,
$z/L\to z$
,
$AL\to A$
,
$BL^2\to B$
and
$CL^5\to C$
, giving (see [Reference Gledzer and Makarov8, equation (1.5)])
\begin{equation} \frac {\partial ^2 \psi }{\partial z^2}+ \frac {\partial ^2 \psi }{\partial r^2}-\frac {1}{r}\,\frac {\partial \psi }{\partial r}+(A^2-B^2r^2+Czr^2)\psi =0. \end{equation}
Now,
$r$
,
$z$
,
$A$
,
$B$
and
$C$
are dimensionless.
To eliminate the term in
${\partial \psi }/{\partial r}$
, put
$\rho =r^2$
. Hence, (3.1) becomes
\begin{equation} \frac {\partial ^2 \psi }{\partial z^2}+4\rho \frac {\partial ^2 \psi }{\partial \rho ^2} +(A^2-B^2\rho +Cz\rho )\psi =0. \end{equation}
Although this is a linear homogeneous PDE, it cannot be solved by the method of separation of variables when
$C\neq 0$
. There are numerous papers on solving (3.2), when
$C=0$
, leading to solutions involving various known special functions; we do not consider these further and concentrate on the non-separable PDE.
For a natural boundary condition, we may specify the vertical velocity
$w(r,z)$
at the ground,
$z=0$
. For example, take
$w(r,0)=U\,\textrm {e}^{-\frac {1}{2}\gamma r^2}$
, where
$U$
and
$\gamma$
are constants. Then, from (2.1),
$\psi (\rho ,0)=-(U/\gamma )\,\textrm {e}^{-\frac {1}{2}\gamma \rho }$
. The choice
$\gamma =B$
, giving
\begin{equation} \psi (\rho ,0)=-(U/B)\,\textrm {e}^{-\frac {1}{2}B\rho }, \end{equation}
fits with later approximations (see (5.7)) and is consistent with a rough estimate from (3.2) when
$\rho \gg 1$
.
4. Frobenius-type expansions
We are unaware of any exact solutions of (3.2) (with
$C\neq 0$
). For this reason, we first attempt to construct solutions using Frobenius expansions. We start in Section 4.1, looking for solutions near the ground at
$z=0$
. Solutions near the axis at
$r=0$
are constructed in Section 4.2.
As (3.2) is a second-order elliptic PDE, we anticipate applying one boundary condition at
$z=0$
(such as (3.3)), but the solutions constructed below involve two arbitrary functions. This difficulty can only be resolved by adjusting the solution so as to satisfy a desired condition for large
$z$
. A similar difficulty arises with the solutions obtained in powers of
$r$
.
4.1. Solution near
$z=0$
We are interested in solving (3.2) for
$\rho \gt 0$
and
$z\gt 0$
, with an appropriate boundary condition at
$z=0$
. Near
$z=0$
, we can try writing
\begin{equation} \psi (\rho ,z)=\sum _{n=0}^\infty z^{n+\delta } f_n(\rho )\quad \mbox{with}\quad f_0(\rho )\not \equiv 0. \end{equation}
We substitute in (3.2). We have
\begin{align*} 4\rho \,\frac {\partial ^2 \psi }{\partial \rho ^2} +(A^2-B^2\rho )\psi &= \sum _{n=0}^\infty z^{n+\delta } ({\mathcal P} f_n)(\rho ), \quad C\rho z \psi =C\rho \sum _{n=1}^\infty z^{n+\delta } f_{n-1}(\rho ),\\ \frac {\partial ^2 \psi }{\partial z^2} & =\sum _{n=0}^\infty z^{n+\delta -2}(n+\delta -1)(n+\delta )f_{n}(\rho )\\ & =\delta (\delta -1)f_0z^{\delta -2}+\delta (\delta +1)f_1z^{\delta -1} +\sum _{n=0}^\infty z^{n+\delta }(n+\delta +2)(n+\delta +1)f_{n+2}, \end{align*}
where
$({\mathcal P} f)(\rho )=4\rho f''(\rho )+(A^2-B^2\rho )f(\rho )$
.
As the term in
$z^{\delta -2}$
must vanish, we obtain
$\delta =0$
or
$\delta =1$
. The choice
$\delta =0$
also removes the term in
$\delta ^{-1}$
(so that
$f_1$
is arbitrary), followed by
$2f_2+{\mathcal P} f_0=0$
and then
${\mathcal P} f_n+C\rho f_{n-1}+(n+2)(n+1)f_{n+2}=0$
for
$n=1,2,\ldots$
. The solution has the form
\begin{equation} \psi (\rho , z)=f_0(\rho )+zf_1(\rho ) +z^2f_2(\rho )+\cdots\kern-1pt , \end{equation}
in which both
$f_0$
and
$f_1$
are arbitrary; they determine
$f_n$
for
$n\geq 2$
.
Evidently, from (4.2),
$\psi (\rho ,0)=f_0(\rho )$
and
$({\partial \psi }/{\partial z})(\rho ,0)=f_1(\rho )$
, and these are related to the velocity at
$z=0$
; see (2.1). In particular, specifying
$f_0$
is equivalent to specifying the vertical velocity at
$z=0$
; see (3.3). Then,
$f_1$
must be related, somehow, to the solution for large
$z$
. (It turns out that the choice
$\delta =1$
is of no interest because it requires that
$f_1\equiv 0$
.)
4.2. Solution near the axis,
$r=0$
Let us look for a solution in the form
\begin{equation*} \psi (\rho ,z)=\sum _{n=0}^\infty \rho ^ng_n(z)\quad \mbox{with}\quad g_0(z)\not \equiv 0. \end{equation*}
We substitute in (3.2). We have
\begin{align*} \frac {\partial ^2 \psi }{\partial z^2}+A^2\psi & =\sum _{n=0}^\infty \rho ^n\big\{g_n''(z)+A^2g_n(z)\big\}, \\ 4\rho \,\frac {\partial ^2 \psi }{\partial \rho ^2} +\rho ( Cz -B^2)\psi & = \sum _{n=1}^\infty \rho ^n\big\{4n(n+1)g_{n+1}+(Cz-B^2)g_{n-1}\big\}. \end{align*}
At the lowest order (
$n=0$
), we have
$g_0''(z)+A^2g_0(z)=0$
, whence
\begin{equation} \psi (0,z)=g_0(z)=\psi _0\,\textrm {e}^{\pm \textrm {i} Az}, \end{equation}
where
$\psi _0$
is a constant. (This formula is consistent with the ansatz (5.1) below). Higher orders give
\begin{equation*} g_n''+A^2g_n+ 4n(n+1)g_{n+1}+(Cz-B^2)g_{n-1}=0,\quad n=1,2,\ldots . \end{equation*}
This can be used to construct
$g_{n+1}(z)$
from
$g_n(z)$
and
$g_{n-1}(z)$
. Thus, we choose
$g_1(z)$
and then proceed; the choice
$g_1=0$
seems to be admissible but may not be interesting.
5. Asymptotic approximations
In this section, we seek approximations to
$\psi$
, assuming that the (dimensionless) coefficient
$A$
is large (
$|A|\gg 1$
). This strategy is motivated by methods coming from high-frequency diffraction, in which a (dimensionless) wavenumber
$k$
is assumed to be large,
$k\gg 1$
. (We shall cite relevant literature below.) The result is a PDE with some resemblance to an equation known as Popov’s equation. To solve our equation, we offer two approaches. In the first, we assume further that the parameter
$B$
is also large (
$|B|\gg 1$
); see Section 5.2. The second approach (which does not assume that
$B$
is large) is more complicated; it uses Laplace transforms with respect to
$\rho$
(Section 5.3).
5.1. A Popov-type equation
Start by writing
\begin{equation} \psi (\rho , z)=\Psi (\rho , z)\,\textrm {e}^{\textrm {i} Az}. \end{equation}
We have
\begin{equation*} \frac {\partial \psi }{\partial z}=\left (\frac {\partial \Psi }{\partial z}+\textrm {i} A\Psi \right )\textrm {e}^{\textrm {i} Az},\quad \frac {\partial ^2 \psi }{\partial z^2}=\left (\frac {\partial ^2 \Psi }{\partial z^2} +2\textrm {i} A\frac {\partial \Psi }{\partial z}-A^2\Psi \right )\textrm {e}^{\textrm {i} Az}. \end{equation*}
Hence, (3.2) becomes (the terms in
$A^2$
cancel)
\begin{equation} \frac {\partial ^2 \Psi }{\partial z^2} +2\textrm {i} A\frac {\partial \Psi }{\partial z} +4\rho \frac {\partial ^2 \Psi }{\partial \rho ^2}-\rho (B^2-Cz )\Psi =0. \end{equation}
Approximate solutions of (5.2) for large
$z$
are constructed in Appendix A using a similarity variable.
If
$|A|\gg 1$
, we may justify discarding the first term in (5.2) when compared to the second term, leaving the approximate equation
\begin{equation} 2\textrm {i} A\frac {\partial \Psi }{\partial z} +4\rho \frac {\partial ^2 \Psi }{\partial \rho ^2}-\rho (B^2-Cz )\Psi =0. \end{equation}
This is a parabolic PDE, reminiscent of the approach introduced by Leontovich and Fock in the 1940s; see, for example, [Reference Kiselev and Zlobina9, equation (4)] or [Reference Ockendon and Tew13, equation (2.5)]. However, the last term
$C\rho z\Psi$
means that (5.3) is closer to a PDE introduced by Mikhail Popov in 1979 [Reference Popov15, equation (1.8)]. Since then Popov’s equation has been studied by several authors [Reference Kiselev and Zlobina9, Reference Ockendon and Tew14, Reference Smyshlyaev and Kamotski17]. In fact, Kiselev and Zlobina [Reference Kiselev and Zlobina9, p. 101] remark that the non-separability of Popov’s equation ‘makes its study an extremely difficult task’!
Actually, Popov’s equation is not exactly the same as (5.3); his PDE has the form
\begin{equation} 2\textrm {i} A\frac {\partial U}{\partial z}+\frac {\partial ^2 U}{\partial y^2}+CyzU=0. \end{equation}
We see that (5.3) contains the additional parameter
$B$
, and this can be exploited to obtain approximate solutions. Just as the ansatz (5.1) enabled us to discard the term in
${\partial ^2 \Psi }/{\partial z^2}$
from (5.2) for
$|A|\gg 1$
, a similar ansatz for
$\Psi$
will enable us to discard the term in
${\partial ^2 \Psi }/{\partial \rho ^2}$
from (5.3) for
$|B|\gg 1$
. The result of doing this will be a first-order PDE which we can solve explicitly; see Section 5.2.
In [Reference Ockendon and Tew14, Section 4.1], the authors obtain a first-order PDE from Popov’s equation by taking a Fourier transform with respect to
$y$
. In Section 5.3, we show that we can do the same for (5.3) (without assuming
$|B|\gg 1$
) using a Laplace transform with respect to
$\rho$
. The resulting PDE can be solved exactly in terms of Airy functions.
5.2. A further approximation assuming
$|B|\gg 1$
Writing
\begin{equation} \Psi (\rho , z)=\Phi (\rho , z)\,\textrm {e}^{-\frac {1}{2}B\rho }, \end{equation}
we have
\begin{equation*} \frac {\partial \Psi }{\partial \rho }=\left (\frac {\partial \Phi }{\partial \rho }-\frac {B}{2}\Phi \right )\textrm {e}^{-\frac {1}{2}B\rho }, \quad \frac {\partial ^2 \Psi }{\partial \rho ^2}=\left (\frac {\partial ^2 \Phi }{\partial \rho ^2} -B\frac {\partial \Phi }{\partial \rho }+\frac {B^2}{4}\Phi \right )\textrm {e}^{-\frac {1}{2}B\rho }. \end{equation*}
Hence, (5.3) becomes (the terms in
$B^2$
cancel)
\begin{equation*} 2\textrm {i} A\frac {\partial \Phi }{\partial z} +4\rho \left (\frac {\partial ^2 \Phi }{\partial \rho ^2} -B\frac {\partial \Phi }{\partial \rho }\right ) +C\rho z\Phi =0. \end{equation*}
Then, for
$|B|\gg 1$
, we may justify discarding the second-derivative term in comparison with the next term, leaving the approximate equation
\begin{equation} 2\textrm {i} A\frac {\partial \Phi }{\partial z} -4B \rho \frac {\partial \Phi }{\partial \rho }+C\rho z\Phi =0. \end{equation}
This is a first-order linear PDE for
$\Phi$
. It can be solved explicitly (see Appendix B.1 for details). Having done this, we obtain the approximation (assuming that
$|A|\gg 1$
and
$|B|\gg 1$
)
\begin{equation} \psi (r,z)=\textrm {e}^{\textrm {i} Az-\frac {1}{2}Br^2}\,\Phi \end{equation}
where
\begin{equation} \Phi (r,z) ={\mathcal F}(\xi )\exp {\left \{\frac {Cr^2}{4B}\!\left (z+\frac {\textrm {i} A}{2B}\right ) \! \right \}}, \end{equation}
$\xi = r^2 \,\textrm {e}^{-2\textrm {i} Bz/A}$
and
$\mathcal F$
is an arbitrary smooth function. In particular,
\begin{equation*} \psi (r,0)={\mathcal F}(r^2)\exp {\left \{ \left (-\frac {1}{2}B+ \frac {\textrm {i} AC}{8B^2}\right )r^2 \right \}} \quad \mbox{and}\quad \psi (0,z)={\mathcal F}(0)\, \textrm {e}^{\textrm {i} Az}. \end{equation*}
Matching with the boundary condition (3.3) gives
\begin{equation} {\mathcal F}(\rho )=-\frac {U}{B} \exp {\left ( - \frac {\textrm {i} AC}{8B^2}\,\rho \right )}. \end{equation}
This determines the function
$\mathcal F$
for positive real arguments.
If we accept the approximation (5.7), we can then calculate
${\partial \psi }/{\partial z}$
at
$z=0$
, and this can be used to define
$f_1(\rho )$
in the Frobenius expansion (4.2). Differentiating (5.7) and making use of (5.8), we obtain the simple result
$f_1(\rho )=\textrm {i} A f_0(\rho )$
.
5.3. Solving (5.3) using Laplace transforms
Consider (5.3), derived assuming that
$|A|\gg 1$
. In order to obtain a first-order PDE, we compute the Laplace transform of
$\Psi (\rho ,z)$
with respect to
$\rho$
, defining
\begin{equation*} \widetilde {\Psi }(s,z)=\int _0^\infty \Psi (\rho ,z)\,\textrm {e}^{-s\rho }\,\textrm {d}\rho . \end{equation*}
Now, for a given function
$f$
, its Laplace transform is
$F(s)={\mathcal L}\{f\}=\int _0^\infty f(t)\,\textrm {e}^{-st}\,\textrm {d} t$
. Then
\begin{equation*} {\mathcal L}\{f''\}=s^2F(s)-sf(0)-f'(0), {\mathcal L}\{tf(t)\}=-F'(s), {\mathcal L}\{tf''(t)\}=-s^2F'(s)-2sF(s)+f(0). \end{equation*}
Using these results, we take the Laplace transform of (5.3) giving
\begin{equation} 2\textrm {i} A\frac {\partial \widetilde {\Psi }}{\partial z} +(B^2-Cz -4s^2)\frac {\partial \widetilde {\Psi }}{\partial s}=8s\widetilde {\Psi }-4\psi _0, \end{equation}
where
$\psi _0=\Psi (0,z)$
, a constant, after combining (4.3) with (5.1). This is a first-order linear PDE for
$\widetilde {\Psi }$
. It can be solved explicitly (see Appendix B.2 for details). Having done this, we obtain the approximation (assuming that
$|A|\gg 1$
)
$\psi (r,z)=\textrm {e}^{\textrm {i} Az}\Psi (\rho , z)$
where the Laplace transform
$\widetilde {\Psi }(s,z)$
of the stream function
$\Psi (\rho ,z)$
is given by
\begin{equation} \widetilde {\Psi }(s,z)=\frac {2\textrm {i}}{A}\psi _0\, [v(z)]^{2}\int _0^z \frac {\textrm {d}\zeta }{[v(\zeta )]^2}+{\mathcal F}(\xi ) [v(z)]^{2}, \end{equation}
with
$v(z)=\mbox{Ai}(\varphi (z))+\xi (z,s)\,\mbox{Bi}(\varphi (z))$
,
\begin{equation*} \xi (z,s)=\frac {s\,\mbox{Ai}(\varphi (z))-E\,\mbox{Ai}'(\varphi (z))}{E\,\mbox{Bi}'(\varphi (z))-s\,\mbox{Bi}(\varphi (z))} \end{equation*}
$E=\frac {1}{2}\textrm {i} (AC)^{1/3}$
and
$\varphi (z) =(AC)^{-2/3}(z C -B^2)$
. As before,
$\mathcal F$
is an arbitrary smooth function.
6. Discussion and conclusions
In this paper, we have studied a very simple model of fire whirls. The model consists of a linear second-order elliptic PDE for a stream function
$\psi (r,z)$
, containing three dimensionless parameters,
$A$
,
$B$
and
$C$
, where
$r$
and
$z$
are positive cylindrical polar coordinates. When
$C=0$
(no thermal effects), the PDE can be solved by the method of separation of variables (as has been known for many years). This method fails when
$C\neq 0$
, and this is the case of interest here.
Two Frobenius-type solutions are constructed, one for small
$z$
(near the ground) and one for small
$r$
(near the axis of rotation). Approximations for large
$z$
are also obtained, using a similarity variable.
Then, motivated by strategies coming from high-frequency diffraction, we assume that the parameter
$A$
is large. Doing this leads to a PDE that has similarities with Popov’s equation but also some differences. One difference involves the parameter
$B$
: if we also assume that
$B$
is large, we obtain a first-order PDE that we solve by standard methods. When
$B$
is not assumed to be large, we obtain a more complicated but solvable first-order PDE, using a Laplace transform with respect to a radial coordinate. The solution of this PDE is found; it involves Airy functions. We have not attempted to invert
$\widetilde {\Psi }$
in order to obtain
$\Psi$
itself.
The calculations presented herein constitute a first step. Clearly, the basic model PDE is very simple and could be improved. On the other hand, encountering a non-separable PDE of Popov type presents mathematical difficulties and challenges.
Competing interests
None.
Appendix A. Solving (5.2) approximately for large
$z$
For large
$z$
, we can introduce a similarity variable
$\tau =\rho z^\beta$
[Reference Battaglia, Rehm and Baum3, Section II.D]. We will choose the parameter
$\beta$
below. Eliminate
$\rho$
in favour of
$\tau$
using
\begin{equation*} \frac {\partial \Psi }{\partial \rho }=z^\beta \frac {\partial \Psi }{\partial \tau } \quad \mbox{and}\quad \frac {\partial ^2 \Psi }{\partial \rho ^2}=z^{2\beta }\frac {\partial ^2 \Psi }{\partial \tau ^2}. \end{equation*}
Substituting in (5.2) gives
\begin{equation} \frac {\partial ^2 \Psi }{\partial z^2} +2\textrm {i} A\frac {\partial \Psi }{\partial z} +4\tau z^\beta \frac {\partial ^2 \Psi }{\partial \tau ^2}-B^2 \tau z^{-\beta } \Psi +C\tau z^{1-\beta }\Psi =0. \end{equation}
We are interested in approximations for large
$z$
. To obtain a non-trivial result as
$z\to \infty$
, we balance the third and last terms in (A.1). Thus, take
$\beta =1-\beta$
giving
$\beta =\frac {1}{2}$
and
$\tau =\rho \sqrt {z}$
. Then, dividing by
$\tau \sqrt {z}$
,
\begin{equation} 4 \frac {\partial ^2 \Psi }{\partial \tau ^2} +\frac {1}{\tau \sqrt {z}}\frac {\partial ^2 \Psi }{\partial z^2} +\frac {2\textrm {i} A}{\tau \sqrt {z}}\frac {\partial \Psi }{\partial z} +C\Psi -\frac {B^2}{z}\Psi =0. \end{equation}
To solve this equation, try
$\Psi (\tau ,z)=\Psi _0(\tau )+z^{-\gamma }\Psi _1(\tau )+\cdots$
, where
$\gamma \gt 0$
will be chosen later. Substituting in (A.2) gives, at leading order,
$4\Psi _0''+C\Psi _0=0$
, whence
\begin{equation} \Psi _0(\tau )=\textrm {e}^{\pm \mu \tau }\quad \mbox{with}\quad \mu ^2=-C/4 \quad \mathrm{and}\quad \tau =r^2\sqrt {z}. \end{equation}
Thus, we obtain the large-
$z$
approximation
$\psi (r,z)\simeq \exp {\{\textrm {i} A z \pm \mu r^2\sqrt {z}\,\}}$
. At next order, the dominant terms are
$z^{-\gamma }$
and
$z^{-1/2}$
; balancing these gives
$\gamma =\frac {1}{2}$
, and then
$\Psi _1(\tau )$
can be calculated, if desired.
Appendix B. Solving first-order linear PDEs
In Section 5.2, we obtained (5.6) whereas we obtained (5.9) in Section 5.3. These first-order PDEs can be written in standard form as
\begin{equation} a(x,y)\frac {\partial u}{\partial x}+b(x,y)\frac {\partial u}{\partial y}=q(x,y, u), \end{equation}
and they can be solved using a standard method. The first step is to find the characteristics by solving
$\textrm {d} y/\textrm {d} x=b/a$
; the solution can be written as
$\xi (x,y)=c$
, where
$c$
is an arbitrary constant. We note that
$\xi$
solves the homogeneous form of (B.1),
\begin{equation} a(x,y)\frac {\partial \xi }{\partial x}+b(x,y)\frac {\partial \xi }{\partial y}=0. \end{equation}
Next, change the independent variables in (B.1) from
$x$
and
$y$
to
$\xi$
and
$\eta$
, where
$\eta$
is unspecified at present. Application of the chain rule, together with (B.2), yields
\begin{equation} \left (a\frac {\partial \eta }{\partial x}+b\frac {\partial \eta }{\partial y}\right )\frac {\partial u}{\partial \eta }=q. \end{equation}
Finally, a judicious choice for
$\eta$
may allow us to integrate this equation with respect to
$\eta$
.
B.1. Solving (5.6)
Let us start with (5.6). We write it as (B.1), in which
$x=z$
,
$y=\rho$
,
$u=\Phi$
,
$q=u$
,
$a=-2\textrm {i} A/(Cxy)$
and
$b=4B/(Cx)$
. The characteristics are determined by solving
\begin{equation*} \frac {\textrm {d} y}{\textrm {d} x}=\frac {b}{a}=\textrm {i}\lambda y \quad \mbox{with}\quad \lambda =\frac {2B}{A}. \end{equation*}
The general solution is
$y=c\,\textrm {e}^{\textrm {i}\lambda x}$
, that is,
$\xi (x,y)=c$
with
$\xi =y\,\textrm {e}^{-\textrm {i}\lambda x}$
. So, the general solution of (B.1) with
$q\equiv 0$
is
$u=f(\xi )$
, where
$f$
is an arbitrary smooth function.
Next, choose
$\eta =x^2$
in (B.3); as
$q=u$
, we obtain
\begin{equation*} \frac {\partial u}{\partial \eta }=\Delta \,u \quad \mbox{with}\quad \Delta =\frac {1}{2xa}=\frac {\textrm {i} C y}{4A}= \frac {\textrm {i} C \xi }{4A}\,\textrm {e}^{\textrm {i} \lambda \sqrt {\eta }}. \end{equation*}
Integrating with respect to
$\eta$
,
\begin{align} u&=f(\xi )\exp {\left \{\frac {\textrm {i} C \xi }{4A}\int ^\eta \!\textrm {e}^{\textrm {i}\lambda \sqrt {t}}\,\textrm {d} t \right \}} =f(\xi )\exp {\left \{ \frac {C \xi \textrm {e}^{\textrm {i}\lambda x}}{2A\lambda }\left (x+\frac {\textrm {i}}{\lambda }\right ) \!\right \}} \nonumber \\[4pt] &=f(\xi )\exp {\left \{ \frac {C y}{4B}\!\left (x+\frac {\textrm {i}}{\lambda }\right ) \!\right \}}. \end{align}
One can check that this expression for
$u$
does solve (B.1).
B.2. Solving (5.9)
Let us write (5.9) in standard form as (B.1), in which
$x=z$
,
$y=s$
,
$u=\widetilde {\Psi }$
,
$a=2\textrm {i} A$
,
$b=B^2-Cx-4y^2$
and
$q=8yu-4\psi _0$
. The characteristics are determined by solving
\begin{equation*} \frac {\textrm {d} y}{\textrm {d} x}=\frac {b}{a}=\frac {2\textrm {i}}{A}\,y^2 +\frac {B^2-Cx}{2\textrm {i} A}. \end{equation*}
This is a Riccati equation [Reference Davis6, p. 59]: it is of the form
$y'(x)+yQ(x)+y^2R(x)=P(x)$
with
$Q=0$
,
$R(x)=R_0$
and
$P(x)=P_0+P_1x$
, where
$P_0=B^2/(2\textrm {i} A)$
,
$P_1=-C/(2\textrm {i} A)$
and
$R_0=-2\textrm {i}/A$
are constants. Making the substitution
$y(x)=v'(x)/(Rv)$
results in a linear second-order equation for
$v(x)$
,
\begin{equation*} Rv''+(QR-R')v'-PR^2v=0, \end{equation*}
which simplifies to
$v''(x)-R_0(P_0+P_1x)v(x)=0$
; this can be solved in terms of Airy functions,
$\mbox{Ai}$
and
$\mbox{Bi}$
. Thus, for example, write
$v(x)=\mbox{Ai}(\varphi (x))$
with
$\varphi (x)=q_0+q_1x$
. Hence
$v''(x)=q_1^2\mbox{Ai}''(\varphi )=q_1^2\varphi \mbox{Ai}(\varphi )$
, using
$\mbox{Ai}''(w)=w\mbox{Ai}(w)$
, and we obtain
\begin{equation*} 0=v''(x)-R_0(P_0+P_1x)v(x)=\big\{q_1^2(q_0+q_1x)-R_0(P_0+P_1x)\big\}\mbox{Ai}(\varphi ). \end{equation*}
Matching powers of
$x$
gives
$q_1^3=R_0P_1$
and
$q_0q_1^2=R_0P_0$
, whence
\begin{equation*} q_1=(R_0P_1)^{1/3}= (C/A^2)^{1/3} \quad \mbox{and}\quad q_0=R_0P_0/q_1^2=-B^2 (AC)^{-2/3}. \end{equation*}
Therefore,
$v(x)=c_1\mbox{Ai}(\varphi (x))+c_2\mbox{Bi}(\varphi (x))$
,
$v'=q_1\{c_1\mbox{Ai}'(\varphi )+c_2\mbox{Bi}'(\varphi )\}$
and
\begin{equation} y(x)=\frac {v'(x)}{R_0v(x)}=E \,\frac {\mbox{Ai}'(\varphi )+c\,\mbox{Bi}'(\varphi ) }{\mbox{Ai}(\varphi )+c\,\mbox{Bi}(\varphi ) } \end{equation}
with
$E=q_1/R_0=\frac {1}{2}\textrm {i} (AC)^{1/3}$
and
$c=c_2/c_1$
. Solving for
$c$
, we obtain
$\xi (x,y)=c$
where
\begin{equation*} \xi (x,y)=\frac {y\,\mbox{Ai}(\varphi )-E\,\mbox{Ai}'(\varphi )}{E\,\mbox{Bi}'(\varphi )-y\,\mbox{Bi}(\varphi )}. \end{equation*}
A tedious calculation verifies that
$\xi (x,y)$
satisfies (B.2), as expected.
As before, change variables to
$\xi$
and
$\eta$
, but with
$\eta =x$
, giving (B.3), which reduces to
\begin{equation} \frac {\partial u}{\partial \eta }=\frac {q}{a}=\frac {1}{2\textrm {i} A}(8yu-4\psi _0). \end{equation}
Now, from (B.5),
$y=\frac {1}{2}\textrm {i} A v'(\eta )/v(\eta )$
and
$c$
can be replaced by
$\xi$
. Hence (B.6) becomes
\begin{equation*} \frac {\partial u}{\partial \eta } -2 \frac {v'(\eta )}{v(\eta )}u=\frac {2\textrm {i}}{A}\psi _0, \end{equation*}
a first-order linear equation with integrating factor
$[v(\eta )]^{-2}$
; its use gives
\begin{equation*} \frac {\partial }{\partial \eta }\left \{ \frac {u}{ [v(\eta )]^{2}}\right \}=\frac {2\textrm {i}\,\psi _0}{A[v(\eta )]^2}. \end{equation*}
Integrating with respect to
$\eta$
,
\begin{equation} u=\frac {2\textrm {i}}{A}\psi _0\, [v(\eta )]^{2}\int _0^\eta \frac {\textrm {d}\zeta }{[v(\zeta )]^2}+f(\xi ) [v(\eta )]^{2} \end{equation}
where
$f$
is an arbitrary function.
Open access
