Abstract

A new transform pair representing solutions to the complex Helmholtz equation in a convex 2D polygon is derived using the theory of Bessel’s functions and Green’s second identity. The derivation is a direct extension of that given by Crowdy (2015, A transform method for Laplace’s equation in multiply connected circular domains. IMA J. Appl. Math., 80, 1902–1931) for ‘Fourier–Mellin transform’ pairs associated with Laplace’s equation in various domain geometries. It is shown how the new transform pair fits into the collection of ideas known as the Fokas transform where the key step in solving any given boundary value problem is the analysis of a global relation. Here we contextualize those global relations from the point of view of ‘reciprocal theorems’, which are familiar tools in the study of the effective properties of physical systems. A survey of the many uses of this new transform approach to the complex Helmholtz equation in applications is given. This includes calculation of effective impedance in electrochemical impedance spectroscopy and in other spectroscopy methods in diffusion wave field theory, application to the 3|$\omega $| method for measuring thermal conductivity and to unsteady Stokes flow. A theoretical connection between this analysis of the global relations and Lorentz reciprocity in mathematical physics is also pointed out.

1. Introduction

This is a companion paper to that of Crowdy (2015b) on a new transform method for Laplace’s equation in multiply connected circular domains, including polygons, and to the paper of Luca & Crowdy (2018) on the biharmonic equation in that same class of domains. The present paper is concerned with the complex partial differential equation (PDE)
(1.1)
governing the 2D field variable |$\phi (x,y)$| and where, importantly, the parameter |$\sigma $| is a complex-valued constant. The domains of interest here will be restricted to polygons. The notation |$\textrm{arg}[\sigma ]$| denotes the argument of a complex number |$\sigma $|⁠. When |$\sigma =0$| (1.1) reduces to Laplace’s equation. If |$\sigma $| is purely real with |$\textrm{arg}[\sigma ] = 0$|⁠, or |$\sigma \in \mathbb{R}^+$|⁠, then (1.1) is the modified Helmholtz equation; when |$\sigma $| is real with |$\textrm{arg}[\sigma ] = \pm \pi $|⁠, or |$\sigma \in \mathbb{R}^-$|⁠, (1.1) is the Helmholtz equation. The focus of this paper is on the case strictly between these two limits when |$0 < \textrm{arg}[\sigma ] < \pi $|⁠. The PDE in (1.1) will then be referred to as the complex Helmholtz equation.

The complex Helmholtz equation (1.1) arises in periodically forced systems, ranging from electrochemical impedance spectroscopy to unsteady slow viscous flows, and in a broad area that has become known as diffusion wave field theory (Mandelis, 2013). The latter is an umbrella term embracing a variety of fields including photoacoustic spectroscopy (Rosencwaig & Gersho, 1976) and photothermal spectroscopy (Almond & Patel, 1996). Stakgold (2000) discusses the complex Helmholtz equation (1.1), for different values of |$\sigma \in \mathbb{C}$|⁠, in Section 7.12 of his monograph. Other discussions of its use in the diffraction and scattering of waves can be found in the book by Babich et al. (2008).

To elucidate how (1.1) arises in diffusion wave field theory consider the 2D diffusion equation
(1.2)
for some real-valued scalar field |$\varPhi (\textbf{x}, t)$| where |$D$| is a diffusion coefficient. Diffusion is a basic physical process and the PDE (1.2) is ubiquitous in applications. Suppose that it is required to solve (1.2) in some 2D spatial domain with a time-periodic forcing on some portion of the domain boundary having a set frequency |$\omega> 0$|⁠. Since (1.2) is linear it is natural to seek a solution, which, after initial transients, oscillates with the same frequency:
(1.3)
where |$\phi (\textbf{x})$| is a real function. On substitution into (1.2) we find
(1.4)
On the other hand, we could just as well have sought a solution of the form
(1.5)
which, on substitution into (1.2), gives
(1.6)
Now introduce the complex number |$\textrm{j} = \sqrt{-1}$| and add |$j$| times (1.6) to (1.4):
(1.7)
The common factor |$e^{\textrm{j} \omega t} $| can be dropped, leaving the complex Helmholtz equation (1.1) with
(1.8)
Since |$\omega , D> 0$| then, as a complex number, |$\textrm{arg}[\sigma ]=\pi /2$|⁠. In other words, to solve a boundary value problem for the diffusion equation (1.2) driven by a time-periodic boundary stimulus, we can pose a solution
(1.9)
and solve the complex Helmholtz equation (1.1) for |$\phi (\textbf{x})$| with the complex-valued |$\sigma $| given in (1.8). The notation |$\textrm{Re}_j[w]$| in (1.9) means that if |$w$| is some complex-valued quantity |$w = u + \textrm{j}v$| then |$\textrm{Re}_j[w] = u$|⁠.

The purpose of this paper is two-fold. First, a novel theoretical approach, based on a newly formulated ‘transform pair’, to solving the complex Helmholtz equation (1.1) in a general convex polygon is given. Second, this new method is used to tackle a range of boundary value problems arising in a variety of applications where the PDE (1.1) is relevant.

1.1 A paradigm for the new method

One application to be focussed on here, and which reduces to solving (1.7), i.e.
(1.10)
is electrochemical impedance spectroscopy, or EIS for short (MacDonald, 2006). In EIS (1.10) is relevant provided the alternating current amplitude is small (Bisquert & Roberto, 1999) with |$\phi $| being a measure of the concentration of a diffusing species, |$\omega $| is the frequency of the alternating current driving the dynamics of this species and |$D$| is a diffusion constant (Barsoukov & Macdonald, 2005). The complex-valued impedance, typically denoted by |$Z$|⁠, is an important diagnostic quantity having the dimension of Ohms that is proportional to the ratio between the concentration and the concentration flux at the stimulus boundary. Mathematically, it is proportional to the ratio between the integral of the Dirichlet data |$\phi $| and integral of the Neumann data |$\partial \phi /\partial \textbf{n}$|⁠. The quantity |$\sqrt{D/\omega }$|⁠, with the dimension of length, is called the penetration depth. By measuring the impedance of an electrochemical device physical parameters such as the electrochemical diffusivity (Huang, 2018; Huang et al., 2016), capacitance, reactance and corrosion rate can be extracted (Barsoukov & Macdonald, 2005); these quantities are important, e.g. in optimizing the design of batteries.
A simplification often made in EIS—to avoid geometrical complications associated with the domain—is to take a so-called Warburg impedance to be relevant (Bard & Faulkner, 2001; Huang, 2018; Warburg, 1899). This assumes that an AC electrochemical response is dominated by pure diffusion of the active species and that the concentration penetration depth is sufficiently small compared to the typical dimensions of the domain that, as shown in Fig. 1, a semi-infinite domain can be assumed. The Warburg impedance (Bard & Faulkner, 2001) is then given by
(1.11)
A Warburg element geometry leading to the impedance given in (1.11). The concentration penetration depth $\sqrt{D/\omega }$ is taken to be small so that a semi-infinite domain can be assumed and the influence of boundaries neglected. A concentration field is prescribed at the half-plane boundary.
Fig. 1.

A Warburg element geometry leading to the impedance given in (1.11). The concentration penetration depth |$\sqrt{D/\omega }$| is taken to be small so that a semi-infinite domain can be assumed and the influence of boundaries neglected. A concentration field is prescribed at the half-plane boundary.

where |$A_W$|⁠, measured in |$\textrm{Ohms}/\textrm{sec}^{1/2}$|⁠, is the Warburg coefficient typical of a specific experiment set-up (Bard & Faulkner, 2001; Huang, 2018). An improved approximation to the impedance is available when the finite thickness |$h$| of a diffusion layer is known as shown in Fig. 2. Two Warburg element formulas used to compute impedance are
(1.12a)
 
(1.12b)
Warburg short and open elements, with boundary conditions, corresponding to $Z_{W_S}$ (1.12) and $Z_{W_O}$ (1.12b).
Fig. 2.

Warburg short and open elements, with boundary conditions, corresponding to |$Z_{W_S}$| (1.12) and |$Z_{W_O}$| (1.12b).

where |$W_S$| and |$W_O$| denote the Warburg short and Warburg open element, respectively (Bard & Faulkner, 2001). These two cases are shown in Fig. 2. Equations (1.12) and (1.12b) are obtained by, respectively, assuming a transmissive boundary (homogeneous Dirichlet) and reflective (homogeneous Neumann) boundary on the wall opposite that of the stimulus wall. If |$S_1$| denotes the stimulus wall of this system where an AC current is imposed—in Fig. 2, |$S_1$| is the entire |$x$| axis—the impedance is defined to be
(1.13)
where |$ds$| denotes an arc-length element along the stimulus wall.

The impedance is an example of an effective parameter useful in many areas of the natural sciences. Finding effective properties in the study of heterogeneous periodic composites is a huge area (Keller, 1964; Milton, 2002). In flow problems the blockage coefficient (Crowdy, 2011) is an effective parameter quantifying the ease of passage of ideal channel flow with occlusions. In porous media flows the tortuosity factor quantifies the apparent decrease in diffusive transport resulting from convolutions of the flow paths through porous media (Cooper et al., 2017, 2016). In laminar convective heat transfer the Nusselt number is the key diagnostic of most use in engineering applications (Shah & London, 1978). In microfluidics the friction reduction properties of superhydrophobic surfaces are measured using the effective slip length (Crowdy, 2017). A mathematical device often used in finding these effective parameters is the notion of reciprocity, i.e. the use of ‘reciprocal theorems’, which is the name given to integral identities that couple together two independent fields defined over the same domain. If one of these is the field of interest, and a second ‘comparison solution’ is chosen judiciously, valuable information about effective parameters can often be obtained (Crowdy, 2017; Keller, 1964; Milton, 2002). The concept of reciprocity is also familiar in electrodynamics and electrical engineering applications, such as circuit theory, where it is also known as Lorentz, or Green’s, reciprocity (Landau & Lifshitz, 1960; Panofsky & Phillips, 1962).

To illustrate what we mean by reciprocity, and the mathematical spirit of the proposed ‘transform method for polygons’ to be presented here, a novel derivation of the Warburg impedance formula (1.12) is now given. Let |$\phi $| be the concentration field satisfying (1.10). Instead of the infinite strip the domain is taken to be a rectangle of width |$L$| and height |$h$| denoted by |${P}=[-L/2,L/2]\times [0,h]$| with boundary conditions
(1.14)
The infinite strip result corresponds to the |$L \to \infty $| limit. Each side of the rectangle is denoted by |$S_n$| for |$n=1,2,3,4$|⁠, indexed anticlockwise, with |$S_1$| being the wall stimulated by an alternating current of frequency |$\omega $|⁠. The impedance is then defined by (1.13). To determine it, the first key ingredient is Green’s second identity:
(1.15)
where |$\phi $| is the solution and |$\psi $| is any suitably differentiable field defined over |${P}$|⁠, |$\partial /\partial \textbf{n}$| denotes the outward normal derivative, |$ds$| is an element of arc length and |$dA$| denotes a surface area element of |${P}$|⁠. If both |$\phi $| and |$\psi $| are taken to satisfy (1.10), the left-hand side of (1.15) vanishes and we obtain
(1.16)
On substituting the boundary conditions (1.14) into (1.16) we obtain
(1.17)
Next we make a strategic choice of the ‘comparison solution’ |$\psi $| to be
(1.18)
where |$k \ne 0 $| is any complex constant and where we have now introduced the complex variable |$ \textrm{i} = \sqrt{-1}$| associated with the spatial variables. It is readily checked that (1.18) is a solution to (1.10). On substituting (1.18) into (1.17) we obtain
(1.19)
Since (1.18) is analytic for all |$k \ne 0$|⁠, if (1.19) is evaluated at |$k=\pm \sqrt{\sigma }$| where the left hand-side vanishes, the right-hand side must also vanish at the same points. Consequently, these two special values of the complex parameter |$k$| yield two equations involving only certain integrals along |$S_1$| and |$S_3$|⁠:
(1.20a)
 
(1.20b)
These two simultaneous algebraic equations (1.20) can be solved to find
(1.21)
leading, on use of the definition (1.13), to the impedance formula (1.12). A derivation of (1.12b) follows along similar lines.

One remarkable thing should be underlined about the above analysis: the value of the impedance was found without even solving for the concentration field |$\phi $|⁠. Instead, it arose from analysis of relations emerging from the evaluation of an integral identity for specially chosen ‘comparison solutions’. This is a familiar advantage afforded by so-called reciprocal theorems in mathematical physics.

The aim of this paper is now easy to explain: to show how to generalize the above strategy when the domain |${P}$| is an arbitrary convex polygon. More specifically, a new ‘transform pair’ will be derived relevant to fields satisfying the complex Helmholtz equation (1.1) in a general convex polygonal domain, with the relevant transform pair being ‘tailor-made’ to fit the geometry of the given polygon (Crowdy, 2016a). The derivation is an elaboration of the ideas just presented: the use of an integral identity coupled with a strategic choice of a ‘comparison solution’ to derive a set of equations—the analogues of the simple system (1.17) and later to be called the ‘global relations’—that can be solved to find valuable information on the solution. Since it is required to reconstruct the whole solution, via an inverse transform, and not just one specific integral diagnostic (such as the impedance just calculated), then we need to use an infinite family of comparison solutions, parametrized here by a complex variable |$k$|⁠, which will emerge naturally from a suitably chosen integral representation of the relevant Green’s function.

1.2 Overview of the paper

The approach to deriving the relevant transform pair will be a direct extension of one given by Crowdy (2015a,b) who derived analogous transform pairs for analytic functions defined over all the classes of domain illustrated in Fig. 3: any convex polygon, any circular-arc domain and domains having boundaries that are a mixture of straight lines and circular arcs, including multiply connected cases. As pointed out in Crowdy (2015b), any non-convex polygon can be split into a union of convex polygons so those are also included, but there are different solution representations in different sub-polygons. Such transform pairs can be used to solve boundary value problems for Laplace’s equation in all those classes of domain (Crowdy, 2015a,b). Crowdy’s work was originally motivated by mathematical challenges in quantifying hydrodynamic slip on so-called superhydrophobic surfaces where circular domains are common due to the occurrence of free surface menisci satisfying a Laplace–Young force balance requiring them to have constant curvature (i.e. they must be circular arcs). These new methods have already proven to be of great utility in that application (Crowdy, 2015b; Luca et al., 2018; Schnitzer, 2016) and others (Crowdy, 2016b; Luca & Llewellyn Smith, 2018).

Example domains to which the ‘Fourier–Mellin transforms’ presented in Crowdy (2015a,b) are applicable: a convex polygon, a lens-shaped circular-arc domain and the ‘disc-in-channel’ geometry (a multiply connected circular-arc domain).
Fig. 3.

Example domains to which the ‘Fourier–Mellin transforms’ presented in Crowdy (2015a,b) are applicable: a convex polygon, a lens-shaped circular-arc domain and the ‘disc-in-channel’ geometry (a multiply connected circular-arc domain).

Fokas & Kapaev (2003) were the first to find the transform pair for a harmonic field in a convex polygon, using a Lax pair formulation and Riemann–Hilbert methods. Analogous transform pairs for the modified Helmholtz |$\textrm{arg}[\sigma ]=0$| and Helmholtz equation |$\textrm{arg}[\sigma ]=\pi $| have similarly been derived (Fokas, 2008). That original approach has since been reappraised from various other viewpoints and this circle of mathematical ideas, which extends to nonlinear integrable equations whence it originally derives, is known collectively as the Fokas method (Fokas, 2008; Fokas & Pelloni, 2014; Fokas & Spence, 2012; Spence & Fokas, 2010). Interestingly, to the best of the authors’ knowledge, no investigators have previously examined the general complex Helmholtz equation |$0 < \textrm{arg}[\sigma ]< \pi $| from this transform perspective. That is the contribution of this paper.

For Laplace’s equation, Crowdy (2015a,b) offered a novel rederivation of the transform pair of Fokas & Kapaev (2003) using elementary considerations based in geometric function theory. This rederivation is not just a mathematical nicety but turns out to be key to unveiling the form of generalized transform pairs for Laplace’s equation in the much richer class of geometries shown in Fig. 3 including multiply connected circular domains. (Indeed, how to rederive Crowdy’s new transform pairs for multiply connected circular domains using, say, the original Lax-pair/Riemann–Hilbert approach of Fokas & Kapaev (2003) remains an open question discussed again in §5). All the ideas of Crowdy (2015a,b) have already been extended to solve boundary value problems for the biharmonic equation (Luca & Crowdy, 2018) in all the same classes of polygonal and circular domains, including multiply connected cases, shown in Fig. 3. This extends earlier work on the biharmonic equation limited to convex polygons initiated by Crowdy & Fokas (2004) and later reconsidered by Fokas & Pinotsis (2006) and Dimakos & Fokas (2015). For the complex Helmholtz equation considered here, we will be restricting attention to convex polygons. How to extend the analysis of the Helmholtz equation (1.1) defined over the general multiply connected circular domains shown in Fig. 3, even in the case where |$\sigma $| is real, is another open problem, although Spence & Fokas (2010) have obtained some results for the Helmholtz and modified Helmholtz equation in simple wedge regions.

A short summary of the relevant work of Crowdy (2015a,b) is given in §2 since an understanding of it is needed to appreciate the modifications necessary to devise the new transform pair for the complex Helmholtz equation (1.1) as done in §3. With the framework set up, §4 describes how to use it to solve a variety of problems in applications. This includes problems in EIS, diffusion wave field theory, an application to the so-called 3|$\omega $|-method, which is another experimental method for measuring material parameters, and then application to a more challenging fourth-order PDE governing unsteady slow viscous Stokes flows.

This paper is the first to write down a transform pair for the complex Helmholtz equation (1.1) even though the Helmholtz (⁠|$\textrm{arg}[\sigma ] = \pi $|⁠) and modified Helmholtz (⁠|$\textrm{arg}[\sigma ] = 0$|⁠) equations have already been treated using the Fokas methodology (Fokas, 2008; Fokas & Pelloni, 2014; Fokas & Spence, 2012). The work here not only offers new transform pairs for the complex Helmholtz equations but also gives a novel derivation of the transform pairs already known for the Helmholtz and modified Helmholtz equations.

2. A review of the ‘Fourier–Mellin transforms’ of Crowdy (2015a,b)

We first review the ‘Fourier–Mellin transform pairs’ derived by Crowdy (2015a,b). These can be used to solve boundary value problems for harmonic fields in all the domains shown in Fig. 3. This geometric function theory derivation is elementary, and all the key steps will be generalized in §3 to the complex Helmholtz equation (1.1).

2.1 Summary of the Fourier–Mellin transform pairs

Consider a bounded |$N$|-sided convex polygon |$P$| with directed straight line edges |$\lbrace S_n | n=1,...,N \rbrace $| inclined at angles |$\lbrace \chi _n | n=1,...,N \rbrace $| to the positive real axis. The triangle shown in Fig. 3 is an example with |$N=3$|⁠. Using Cauchy’s theorem and a ‘spectral representation’ of the Cauchy kernel Crowdy (2015b) derived the following transform pair to represent a function |$f(z)$| analytic in |$P$|⁠:
(2.1)
where |$L$| is the ray along the positive real axis in the |$k$|-plane—see Fig. 4—and where the spectral matrix functions  |$\rho _{mn}(k)$| satisfy the global relations  
(2.2)
The basic geometrical units used to build up more complicated polygonal and circular domains—see Fig. 6—with the corresponding ‘fundamental contours’ in the $k$ plane in the corresponding spectral representation (2.8) of the Cauchy kernel.
Fig. 4.

The basic geometrical units used to build up more complicated polygonal and circular domains—see Fig. 6—with the corresponding ‘fundamental contours’ in the |$k$| plane in the corresponding spectral representation (2.8) of the Cauchy kernel.

Only the diagonal elements of the spectral matrix appear in the inverse transform formula for |$f(z)$|⁠. A simple change of variable |$\lambda = e^{-\textrm{i} \chi _j} k$| in each of the spectral functions allows us to write the same inverse transform pair as integrals over |$N$| rotations of the fundamental contour |$L$| in a complex |$\lambda $| plane. The result above then coincides precisely with the transform pair derived by Fokas & Kapaev (2003). The appearance of the exponential terms in (2.1) reminds us of Fourier transforms.

In direct analogy, Crowdy (2015a) extended this result to circular polygons such as the lens-shaped domain of Fig. 3. While we will not discuss the complex Helmholtz equation in circular domains in this paper, it is nevertheless instructive to briefly describe how the rationale for convex polygons extends naturally to circular domains. Given a bounded convex |$N$|-sided circular polygon with directed edges |$\lbrace C_n | n=1, \ldots , N \rbrace $| that are arcs of circles with centres |$\lbrace \delta _n | n=1,...,N \rbrace $| and radii |$\lbrace q_n | n=1,...,N \rbrace $| the following transform pair can be derived to represent a function |$f(z)$| analytic in it:
(2.3)
where |$L_1, L_2$| and |$L_3$| are the contours in the |$k$|-plane making up the relevant modified fundamental contour—see Fig. 4—and where the spectral matrix functions  |$\rho _{mn}(k)$| satisfy the global relations
(2.4)
Again, only the diagonal elements of spectral matrix |$\rho _{mn}(k)$| appear in the inverse transform formula for |$f(z)$|⁠. The appearance of the power law terms |$z^k$| in (2.3) reminds us of classical Mellin transforms.
Crowdy’s geometrical derivation makes it obvious how to generalize the ideas to domains that have boundaries that are a combination of straight lines and circular arcs, including multiply connected cases such as the ‘disc-in-channel’ geometry shown in Fig. 3, which is ubiquitous in applications (e.g. the study of superhydrophobic surfaces; Crowdy, 2015a). Such a domain is the intersection of two half-planes and the exterior of a disc. The relevant transform representation of a function |$f(z)$| that is analytic and single-valued in such a domain can then be shown (Crowdy, 2015a) to be
(2.5)
where the definitions of the spectral functions are given by Crowdy (2015a) and involve integrals of the solution around the various boundary components. This formula, where both Fourier-type and Mellin-type contributions appear simultaneously, underscores the appropriateness of the designation ‘Fourier–Mellin transforms’ (Crowdy, 2015a). It would be interesting to find transform pairs for Helmholtz-type PDEs (1.1), for arbitrary |$\sigma $|⁠, in this same class of (multiply connected) circular domains. This remains an open challenge.

2.2 Spectral representation of the Cauchy kernel

Suppose a point |$z^{\prime}$| lies on some finite length slit on the real axis and |$z$| is in the upper half plane as shown on the left in Fig. 5 then it is clear from the geometry that
(2.6)
Geometrical positioning of $z$ and $z^{\prime}$ for the validity of (2.8) and (2.12), shown left and right, respectively.
Fig. 5.

Geometrical positioning of |$z$| and |$z^{\prime}$| for the validity of (2.8) and (2.12), shown left and right, respectively.

It follows, on use of elementary properties of the exponential function, that
(2.7)
or,
(2.8)
where |$L$| is the straight line contour along the positive real axis as shown in Fig. 4. It is easy to check that the contribution from the upper limit of integration vanishes for the particular choices of |$z^{\prime}$| and |$z$| to which we have restricted consideration. The left hand side of (2.8) is the Cauchy kernel appearing in Cauchy’s integral formula.
On the other hand, suppose |$z^{\prime}$| lies on some other finite length slit making angle |$\chi $| with the positive real axis and suppose that |$z$| is in the slanted half plane shown shaded on the right in Fig. 5 (the half-plane ‘to the left’ of the slit as one follows its tangent with uniform inclination angle |$\chi $|⁠). Now the transformation
(2.9)
where the (unimportant) constant |$a_0$| is shown geometrically in Fig. 5, takes the slit to the real axis, and |$z$| to the upper-half plane, and
(2.10)
Hence, on use of (2.8) with the substitutions (2.9), we can write
(2.11)
or, on cancellation of |$a_0$| and rearrangement,
(2.12)
The left hand side of (2.12) is again the Cauchy kernel and (2.12) gives a ‘spectral representation’ of it when |$z$| and |$z^{\prime}$| have the geometrical arrangement shown on the right in Fig. 5.

2.3 Convex polygon as an intersection of half-planes

The next observation is geometrical, and simple. Consider a bounded convex polygon |$P$| with |$N$| directed sides |$\lbrace S_j | j=1, \dots , N \rbrace $|⁠. Figure 6 shows an example with |$N=3$|⁠. It can be viewed as the intersection of |$N$| half-plane regions akin to those shown in Fig. 5. Any convex polygon can be constructed in this way. Any non-convex polygon can be viewed as a union of convex subpolygons.

A convex polygon $P$, here a triangle, containing an interior point $z$ as an intersection of $N=3$ half-planes with $N$ angles $\lbrace \chi _j | j=1,2,3 \rbrace $. Formula (2.12) can be used in Cauchy’s integral formula with $\chi =\chi _j$ when $z^{\prime}$ is on side $S_j$ (for $j=1,2,3$).
Fig. 6.

A convex polygon |$P$|⁠, here a triangle, containing an interior point |$z$| as an intersection of |$N=3$| half-planes with |$N$| angles |$\lbrace \chi _j | j=1,2,3 \rbrace $|⁠. Formula (2.12) can be used in Cauchy’s integral formula with |$\chi =\chi _j$| when |$z^{\prime}$| is on side |$S_j$| (for |$j=1,2,3$|⁠).

2.4 Integral representation and derivation of the transform pair

To construct the transform pair we need an integral formula for an analytic function |$f(z)$| in the polygon |$P$|⁠. This is furnished by Cauchy’s integral formula: given a function |$f(z)$| analytic in |$P$|⁠, then for |$z \in P$|⁠,
(2.13)
where |$\partial P$| denotes the boundary of |$P$|⁠. The transform pair now follows easily on combining all the ingredients prepared earlier.
First, the boundary integral (2.13) can be separated into a sum over the |$N$| sides:
(2.14)
If side |$S_j$| has inclination |$\chi _j$| then (2.12) can be used, with |$\chi \mapsto \chi _j$|⁠, to re-express the Cauchy kernel, that is |$1/(z^{\prime}-z)$|⁠, uniformly for all |$z \in P$| and |$z^{\prime}$| on the respective sides:
(2.15)
Since the integrals are uniformly convergent for |$z$| inside |$P$| and |$z^{\prime}$| on its boundary, we can reverse the order of integration to find
(2.16)
where, for integers |$m,n$| between |$1$| and |$N$|⁠, we define the spectral matrix (Crowdy, 2015a,b) to be
(2.17)
and where |${L} = [0, \infty )$| is the fundamental contour (Crowdy, 2015a,b) for straight line edges shown in Fig. 4. We have then arrived at the transform pair (2.1). The spectral matrix elements have their own analytical structure. Observe that, for any |$k \in \mathbb{C}$|⁠, and for any |$m=1,...,N$|⁠,
(2.18)
where we have used Cauchy’s theorem and the fact that |${f(z^{\prime})} e^{-\textrm{i} e^{-\textrm{i} \chi _m} k z^{\prime}} $| (for |$m=1,...,N$|⁠) is analytic inside |$P$| for any |$k \in \mathbb{C}$|⁠. The result (2.18) will be referred to as global relations.

3. A new approach to the complex Helmholtz equation

With this construction in mind, the aim now is to derive a new transform method for the inhomogeneous complex Helmholtz equation given by
(3.1)
where |$\sigma \in \mathbb{C}$|⁠. When |$\arg (\sigma ) = 0,\pi $|⁠, we obtain the familiar modified Helmholtz and Helmholtz equation, respectively. |${P}$| is an arbitrary, |$N$|-sided, convex polygon, with sides |$\lbrace S_n | n=1,...,N\rbrace $|⁠. The function |$\phi $| is twice differentiable and is subject to various boundary conditions on each side of the polygon, e.g.
(3.2a)
 
(3.2b)
 
(3.2c)
where |$\textbf{n}$| denotes the unit normal vector, and |$c$| is some constant. Throughout the derivation of the method, we assume a priori that |$\lbrace g_n \rbrace $| are sufficiently smooth and defined in such a way that there exists a unique solution to the boundary value problem we are considering. The strategy is to emulate all the steps of the construction in §2.

The main challenges are as follows. What is the integral representation that will replace use of Cauchy’s integral formula (2.13) in the foregoing construction? What formula will replace the spectral representation of the Cauchy kernel (2.12) valid in a given half-plane region?

3.1 Spectral representation of the Green’s function

The Green’s function |$\mathcal{G}(\textbf{x}^{\prime}; \textbf{x})$| for the Helmholtz equation (1.1) is defined to be the solution of
(3.3)
For |$0 \leqslant \textrm{arg}[\sigma ] < \pi $| it is well-known to be given in terms of the order-zero modified Bessel function:
(3.4)
This is an axisymmetric function, depending only on |$ |\textbf{x}^{\prime}-\textbf{x}|$|⁠, that is singular as |$\textbf{x}^{\prime} \to \textbf{x}$| and decays as |$|\textbf{x}^{\prime}| \to \infty $|⁠. The case |$\textrm{arg}[\sigma ] = \pi $|⁠, corresponding to the Helmholtz equation, is a special case for which non-decaying, but also non-singular, solutions in the far-field are admissible. These are wave-type solutions and this special case that has been treated in detail elsewhere (Fokas, 2008; Fokas & Pelloni, 2014; Fokas & Spence, 2012; Spence & Fokas, 2010) so we will not revisit it here. From our analysis we will see, however, how this interesting case is approached in the limit |$\textrm{arg}[\sigma ] \to \pi $|⁠.
Watson (1995) [see pages 182–183] has shown that an integral representation of the function |$K_0$| is
(3.5)
where |$\theta $| is a geometrical angle associated with a contour |$L_W$| joining the two essential singularities of the integrand at 0 and |$\infty $| in a complex |$\alpha $| plane as shown in Fig. 7. It is easy to see how the limits on the argument of |$u$| given in (3.5) emerge. Suppose the integration contour approaches the essential singularity of the integrand at the origin at an angle |$\theta ^{(0)}$| then for the integrand to decay (and the integral to exist) we need
(3.6)
The integration contour $L_W$ of Watson (1995) in a complex $\alpha $ plane as used in (3.5).
Fig. 7.

The integration contour |$L_W$| of Watson (1995) in a complex |$\alpha $| plane as used in (3.5).

Similarly, if the integration contour approaches the essential singularity of the integrand at infinity at an angle |$\theta ^{(\infty )}$| then for the exponential in the integrand to decay we need
(3.7)
Conditions (3.6) and (3.7) produce the same  |$\pi $|-length range of possible arguments of |$u$| provided
(3.8)
for some angle |$\theta $|⁠. This is precisely the situation for Watson’s contour choice shown in Fig. 7. Any choice of |$\theta $| in the range |$-\pi < \theta < \pi $| is allowed so there are lots of different contour choices that can, on use of Cauchy’s theorem, all be deformed into each other without altering the value of the integral. The key point for present considerations is that, for a fixed choice of |$\theta $|⁠, the representation (3.5) of |$K_0(u)$| is limited to a |$\pi $|-width range on the argument of |$u$|⁠.
On combining (3.4) and (3.5) and on letting |$u = 2 \sqrt{\sigma } X$| where we have introduced the shorthand notation |$X = |\textbf{x}^{\prime}-\textbf{x}|$|⁠, it follows that
(3.9)
The contour |$\mathcal{L}_\alpha $| joining 0 to |$\infty $| has not yet been chosen but, in view of what was presented above, it must be such that the integrand decays as its two essential singularities are approached along the contour. It is convenient to introduce the change of variables:
(3.10)
where we have introduced the natural correspondences
(3.11)
Geometrically, the change of variables (3.10) corresponds to a simple rotation and a rescaling of the integration plane so the integration contour in the |$k$| plane, which we call |$\mathcal{L}^{(\sigma )}$|⁠, remains one that joins |$k=0$| to |$k =\infty $|⁠. From (3.9) we now define the Green’s function |$G(z^{\prime};z)$| by
(3.12)
where we have used the fact that |$z-z^{\prime}=X e^{\textrm{i} \phi }$| and where |$\mathcal{L}^{(\sigma )}$| is some contour joining |$k=0$| to |$k=\infty $|⁠.
Suppose we now wish to pick this contour |$\mathcal{L}^{(\sigma )}$| in order to ensure that the integral representation is valid for the particular  |$\pi $|-width range of arguments given by
(3.13)
Geometrically, we know that points |$z$| and |$z^{\prime}$| as shown to the left of Fig. 5 satisfy such a condition, a fact we will use shortly. Let this choice of contour be called |$\mathcal{L}_0^{(\sigma )}$| and suppose it approaches the origin at an angle |$\theta ^{(0)}$| and approaches infinity at an angle |$\theta ^{(\infty )}$|⁠. For the integrand to decay at these essential singularities when |$z$| and |$z^{\prime}$| satisfy (3.13) then it is easy to see that we need
(3.14)
Such a contour |$\mathcal{L}_0^{(\sigma )}$| is shown in Fig. 8. Provided it meets the essential singularities at these angles the exact shape of such a contour is inconsequential but, for reasons to be elucidated later, we choose it to be the image of the real line |$(-\infty , \infty )$| in a parametric |$a$| plane, say, to the |$k$| plane under the transformation
(3.15)
The basic geometrical unit (shown left) used to build up a more complicated convex polygonal domain with the corresponding fundamental contour $\mathcal{L}_0^{(\sigma )}$ in the $k$ plane (shown right) for a typical $\sigma \in \mathbb{C}$, for the integral representation (3.12) of the Green’s function of the complex Helmholtz equation (3.1) to be valid. This schematic generalizes Fig. 4 to the PDE (3.1).
Fig. 8.

The basic geometrical unit (shown left) used to build up a more complicated convex polygonal domain with the corresponding fundamental contour |$\mathcal{L}_0^{(\sigma )}$| in the |$k$| plane (shown right) for a typical |$\sigma \in \mathbb{C}$|⁠, for the integral representation (3.12) of the Green’s function of the complex Helmholtz equation (3.1) to be valid. This schematic generalizes Fig. 4 to the PDE (3.1).

It is easy to verify that |$k \to |k|$| so that |$\textrm{arg}[k] \to 0$| as |$a \to \infty $|⁠, and that |$k \to \sigma /|k|$| so that |$\textrm{arg}[k] \to \arg{\sigma }$| as |$a \to -\infty $|⁠. Also easily checked is that
(3.16)
so that, on this choice of contour |$\mathcal{L}_0^{(\sigma )}$|⁠,
(3.17)
Following the designation introduced by Crowdy (2015b) for the 2D Laplacian operator, and shown in Fig. 4, we refer to |$\mathcal{L}_0^{(\sigma )}$| shown in Fig. 8 as the fundamental contour for the complex Helmholtz operator.

It is instructive to examine how this fundamental contour |$\mathcal{L}_0^{(\sigma )}$| changes with |$\arg [\sigma ]$|⁠; this is shown in Fig. 9. For the modified Helmholtz equation, or |$\arg [\sigma ]=0$|⁠, the contour is a straight line along the positive real axis. For |$0 < \textrm{arg}[\sigma ] < \pi $| the shape of the contour changes as shown by the evolution shown in Fig. 9. For the Helmholtz equation, as |$\textrm{arg}[\sigma ] \to \pi $|⁠, the curve becomes a union of three segments: a line emanating from the origin to |$-|\sqrt{\sigma }|$|⁠, the ray from |$k=|\sqrt{\sigma }|$| to infinity along the real axis and the semi-circle of radius |$|\sqrt{\sigma }|$| connecting the two segments. This corresponds to the contour for outgoing waves as discussed by Spence & Fokas (2010). If we had instead studied (1.1) with the range of |$\sigma $| with |$-\pi \leqslant \textrm{arg}[\sigma ] < 0$|—this is just the complex conjugate of the case we chose to study—then the corresponding limit would be the contour for ingoing waves in Spence & Fokas (2010). Readers are referred to Spence & Fokas (2010) for more details on the Helmholtz equation, but it is interesting to see how these two contour choices emerge from this analysis of the complex Helmholtz equation.

Evolution of the fundamental contour $\mathcal{L}_0^{(\sigma )}$ defined in (3.17) as $\arg (\sigma )$ varies. The case $\arg (\sigma )=0$ gives the straight line contour for the modified Helmholtz equation, while $\arg (\sigma )=\pi $ retrieves one of the two contours for the Helmholtz equation assigned to outgoing wave solutions in Spence & Fokas (2010).
Fig. 9.

Evolution of the fundamental contour |$\mathcal{L}_0^{(\sigma )}$| defined in (3.17) as |$\arg (\sigma )$| varies. The case |$\arg (\sigma )=0$| gives the straight line contour for the modified Helmholtz equation, while |$\arg (\sigma )=\pi $| retrieves one of the two contours for the Helmholtz equation assigned to outgoing wave solutions in Spence & Fokas (2010).

Suppose, on the other hand, that we want to pick the contour in order to ensure the integral representation for the Green’s function is valid for the a different width-|$\pi $| range of arguments given by
(3.18)
where |$0 \leqslant \chi < 2\pi $| is some given angle. Let this choice of contour be called |$\mathcal{L}_\chi ^{(\sigma )}$| and suppose it approaches the origin at an angle |$\theta ^{(0)}$| and approaches infinity at an angle |$\theta ^{(\infty )}$|⁠. Then, for the integrand to decay when |$z$| and |$z^{\prime}$| satisfy (3.13) we need
(3.19)
Such a contour |$\mathcal{L}_\chi ^{(\sigma )}$| can be chosen to be a rotation through angle |$-\chi $| of the fundamental contour |$\mathcal{L}_0^{(\sigma )}$|⁠.

At this point we wish to emphasize how our approach differs from the work of Spence & Fokas (2010) who formulate the Fokas transform for the Poisson, Helmholtz and modified Helmholtz equation. First, those authors do not consider the case of complex |$\sigma $|⁠, so the transform pair to be derived next is a new result, as is its derivation. Second, while Spence & Fokas (2010) also circumvent the Lax-pair and Riemann–Hilbert methodology of Fokas & Kapaev (2003) they derive the relevant representation of the Green’s function using a spectral representations of the Dirac-delta distribution and perform a series of manipulations (they integrate out one of the Fourier variables in a 2D Fourier transform) to arrive at the desired integral representation. In contrast, here we have exploited from the outset the connection of the Green’s function for the Helmholtz operator to the theory of (integral representations of) Bessel’s functions.

3.2 Convex polygon as an intersection of half-planes

The domain of interest is a convex polygon |$P$| so the same geometrical observation made in §2.3 still holds: any |$N$|-sided convex polygon |$P$| is the intersection of |$N$| half planes as shown in Fig. 5.

Suppose a point |$z^{\prime}$| lies on a finite length slit on the real axis, and let |$z$| be a point in the upper half-plane as shown to the left in Fig. 5. Then these two points satisfy (3.13) and we know that the representation (3.12) holds for the Green’s function with the choice |$\mathcal{L}= \mathcal{L}_0^{(\sigma )}$|⁠.

If, on the other hand, |$z^{\prime}$| lies on a finite length slit on the infinite line inclined at angle |$\chi $| to the real axis and that |$z$| lies in the right half plane to the left of this infinite line as shown on the right in Fig. 5. These two points satisfy (3.18) and we know that the representation (3.12) holds for the Green’s function with the choice |$\mathcal{L}= \mathcal{L}_\chi ^{(\sigma )}$|⁠.

3.3 Integral representation and derivation of the transform pair

It is Green’s second identity that replaces the use of the Cauchy integral formula used in §2. With the identifications (3.11), and thinking of |$\phi $| and |$\psi $| now as functions of |$z^{\prime}$| and |$\overline{z^{\prime}}$|⁠, another form of the left hand side of Green’s second identity (1.15) is
where |$dA^{\prime} = dx^{\prime} dy^{\prime}$|⁠. Using the complex form of Green’s theorem (Ablowitz & Fokas, 2003), we obtain the identity
(3.20)
We now suppose that |$\phi (z^{\prime})$| is the required solution of the complex Helmholtz equation and set
(3.21)
It follows from (3.20) that
(3.22)
This is the integral formula for the solution |$\phi (z)$| of the complex Helmholtz equation that is analogous to the expression for the analytic function |$f(z)$| by the Cauchy integral formula in (2.13).
Just as we did in (2.14), we can write (3.22) as a sum of the integrals around the |$N$| sides of the polygon:
(3.23)
On side |$S_n$|⁠, oriented at angle |$\chi _n$| to the positive real axis, we know from the representation (3.12) that
(3.24)
Each of these representations can now be substituted into (3.23):
(3.25)
The next step is to use the uniform convergence of all the integrals to justify reversing the order of integration leading to
(3.26)
where
(3.27)
The global relation can be obtained by using (3.20), splitting the contour |$\partial{P}$| into a sum around the |$N$| sides |$\lbrace S_n | n=1, \ldots , N \rbrace $|⁠, letting |$\psi (z^{\prime})=\exp (-\textrm{i}kz^{\prime}+\textrm{i}\sigma \overline{z}^{\prime}/k)$| for any |$k \ne 0$|⁠, using the governing equation (3.1) and using the definition of the spectral function (3.27) to find that
(3.28)
For the analysis of the global relation it is also useful to make a different choice of particular solution |$\psi (z^{\prime})= \exp (\textrm{i}k\overline{z^{\prime}}-\textrm{i}\frac{\sigma z^{\prime}}{k})$| in Green’s identity (3.20). On the use of the governing equation (3.1) this leads to
(3.29)
We can then define a second set of spectral functions given by
(3.30)
which satisfy a second global relation expressed as
(3.31)

An alternative way of obtaining (3.31) is by taking (3.28) and letting |$k \mapsto \frac{\sigma }{k}$|⁠. Both (3.39) and (3.41) are useful in analysing the global relations for a given boundary value problem.

3.4 Unbounded polygonal domains

For unbounded polygons the global relations will only be valid in a restricted section in the complex |$k$|-plane. We will briefly show how to find these regions.

Let |$S_{m}$| and |$S_{m+1}$| denote the sides of the polygon that extend to infinity. For each side we introduce the parametrizations
(3.32)
where |$z_j$| denotes the vertices of the polygon, and |$0\leqslant s\leqslant \infty $|⁠. We assign |$z_{m+1} = \infty $|⁠, hence the jump in index. Using the parametrizations (3.32), we express the corresponding spectral functions as
(3.33)
By examining the exponent of each of the spectral functions, we see that the functions are analytic as long as |$k$| satisfies
(3.34)
It is the overlap of these two regions that is the restricted region to which the global relation is valid.

Figure 10(a) depicts the area to which the spectral function |$\rho _m$| is valid (using |$\chi _m=0$| as an example), while Fig. 10(b) shows the area to which |$\rho _{m+1}$| is valid (using |$\chi _{m+1}=2\pi /3$| as an example). The domain of validity of the global relation is the intersection of these two domains shown in Fig. 10(c).

Figure (a) and (b) depicts the area in which the spectral function $\rho _m$ and $\rho _{m+1}$ are defined, respectively, while (c) illustrates the intersection between the two, hence giving the domain of validity of the global relation.
Fig. 10.

Figure (a) and (b) depicts the area in which the spectral function |$\rho _m$| and |$\rho _{m+1}$| are defined, respectively, while (c) illustrates the intersection between the two, hence giving the domain of validity of the global relation.

The reader will recognize that the boundaries of the regions defined by (3.34) are precisely the contours |$\mathcal{L}_{\chi _m}^{(\sigma )}$|⁠, and it is for this reason we chose these contours as our integration paths earlier. When solving particular unbounded problems, we can find explicit expressions for the spectral functions along these contours, which allows us to obtain explicit integral representations of the solutions inside the domain, respectively.

In Fig. 11 we show how the regions of analyticity change as we modify |$\arg (\sigma )$|⁠, where we used the quarter plane (⁠|$\chi _m=0$|⁠, |$\chi _{m+1}=\pi /2$|⁠) as an example.

The domain of validity of the global relation for the quarter plane ($\chi _m = 0$, $\chi _{m+1}=\pi /2$) as $\textrm{arg}[\sigma ]$ varies. For $\textrm{arg}[\sigma ]=0,\pi $ the known domains for the modified Helmholtz and Helmholtz equation, respectively, are retrieved (Fokas, 2008).
Fig. 11.

The domain of validity of the global relation for the quarter plane (⁠|$\chi _m = 0$|⁠, |$\chi _{m+1}=\pi /2$|⁠) as |$\textrm{arg}[\sigma ]$| varies. For |$\textrm{arg}[\sigma ]=0,\pi $| the known domains for the modified Helmholtz and Helmholtz equation, respectively, are retrieved (Fokas, 2008).

3.5 Other transform pairs

Since the solution |$\phi (z)$| of (3.1) is single-valued in the polygon |$P$|⁠, then there are alternative forms of the integral identity (3.20) obtainable by using the relation
(3.35)
which will hold for any single-valued function |$\psi $| in |$P$|⁠. Those modified integral representations give rise to modified transform pairs. For any given boundary value problem, it is usual to find that certain transform pairs are more convenient than others. Using the differentials
(3.36)
and (3.35), we can rewrite the right-hand side of (3.20) to obtain the following identity:
(3.37)
Based on this form of Green’s identity another choice of transform pair, derivable in exactly the same way as (3.26) and (3.27), is
(3.38a)
 
(3.38b)
with corresponding global relation, for a bounded polygon,
(3.39)
The second set of spectral functions given by
(3.40)
satisfy a second global relation, for a bounded polygon, expressed as
(3.41)
Note that the transform pair (3.38a) and (3.38b) reduces to the pair (2.1) on taking |$\phi = f(z)$| and letting |$\sigma -> 0$| and performing a sequence of integral manipulations.

4. Applications

Armed with these novel transform pairs for the complex Helmholtz equation—either (3.26)–(3.27) or (3.38a)–(3.38b)—and the relevant global relations satisfied by the spectral functions, how to use this machinery to solve a given boundary value problem is, by now, quite well established in the fast-growing body of literature on the Fokas method (Fokas, 2008; Fokas & Pelloni, 2014; Fokas & Spence, 2012). This section, together with the supplementary materials, illustrates this analysis of the global relations; its main goal, however, is to exemplify how this new transform method adds to the small supply of available analytical techniques to solve the complex Helmholtz equation (1.1) in a variety of applications.

4.1 Electrochemical impedance spectroscopy: right-isosceles triangular pores

In §1 the calculation of Warburg impedance in EIS using reciprocal theorem ideas, or integral identities, was discussed. The geometries there were simple: a semi-infinite half-plane, or finite-width infinite strips. Cooper et al. (2017) examine how the pore geometry of the electrodes can affect the impedance measurements in EIS. They develop methods to replicate early results by Keiser et al. (1976) and obtain the impedance for various idealized pore geometries, including a rectangle and triangle geometry. Those authors have built a flexible numerical software package called TauFactor (Cooper et al., 2017, 2016) that can compute impedances for quite general geometries.

The present work complements these alternative approaches and adds some new theoretical ideas to the study of this physically important class of problems. As a detailed example, we now show how to use the new transform pair to solve the boundary value problem for EIS in the triangular geometry shown in Fig. 12. The steps in the following analysis of the global relation are not new: this same geometry has been considered previously for Laplace’s equation or |$\sigma =0$| using the Fokas method by Fokas & Kapaev (2003) and for the modified Helmholtz equation |$\sigma \in \mathbb{R}^+$| (Ben-Avraham & Fokas, 2001). The following is the first analysis of a boundary value problem for the complex Helmholtz equation |$\sigma \in \mathbb{C}$| in this geometry. Many of the steps in the analysis of the global relations are similar to that in previous work for other PDEs. In §4.6 we will compare impedance formulas obtained using our new transform approach to the results given by TauFactor.

The right isosceles triangular pore geometry for EIS. The hypotenuse is the stimulus wall with a Dirichlet-type boundary condition and homogeneous Neumann conditions are imposed on the two legs.
Fig. 12.

The right isosceles triangular pore geometry for EIS. The hypotenuse is the stimulus wall with a Dirichlet-type boundary condition and homogeneous Neumann conditions are imposed on the two legs.

For this particular choice boundary conditions, relevant to EIS and shown in Fig. 12, it is possible to exploit symmetry to reduce problem to one in a rectangle where more traditional methods such as separation of variables can be used. This is because of the homogeneous Neumann-type boundary conditions on the two non-stimulated side-walls. However, all the same steps below can be generalized to the case of general inhomogeneous Neumann data on the side-walls where it is not, in general, possible to use separation of variables to solve such a problem. Our approach therefore offers genuine advantage in facilitating the solution of problems not currently solvable using other known methods. The details of this generalized analysis can be found in the supplementary materials. Incidentally, we believe our analysis there also generalizes that already presented by Fokas & Kapaev (2003) and Ben-Avraham & Fokas (2001) for other PDEs.

Let the legs of the right isosceles triangle have length |$L$| with each side labelled as in Fig. 12. The governing equation is (1.10), and the boundary conditions for the closed pore are given by
(4.1)
The goal is to compute the dimensionless impedance, which can be expressed as (Cooper et al., 2017)
(4.2)
where |$ds$| denotes the differential arc-length element along |$S_2$|⁠. Since |$\phi $| is already known along |$S_2$|⁠, we need to use the new transform method to obtain the unknown Neumann data |$\partial \phi /\partial \textbf{n}_2$| there.

4.2 Defining the spectral functions

The first step is to split the spectral functions into known and unknown parts. For this problem it turns out to be more convenient to use the transform pair (3.38a) and (3.38b) since the boundary data are either of Dirichlet or Neumann type. Hence we split the spectral functions into parts dependent on the Dirichlet and the Neumann data, respectively. On |$S_1$|⁠, e.g. we can write
(4.3)
On use of the identity
(4.4)
(4.3) becomes
(4.5)
We can perform a similar split for the spectral functions on side |$S_1$| and |$S_2$|⁠:
(4.6a)
 
(4.6b)
where we have introduced the notation
(4.7)
On |$S_2$|  |$\phi $| is given, and |$\partial \phi /\partial \textbf{n}_2$| is unknown. The function |$\rho _2(k)$| can therefore be expressed as
(4.8)
where
(4.9)
On |$S_1$| and |$S_3$| we have the opposite situation: the Neumann data are given and the Dirichlet data are unknown. Consequently the spectral functions are given by
(4.10)
where
(4.11)
In order to perform the analysis we need to use the second set of spectral functions given by (3.40). Performing a similar procedure as outlined above, the second set of spectral functions for the closed pore is given by
(4.12a)
 
(4.12b)
 
(4.12c)
where
(4.13)

4.3 Analysis of global relation

It is expedient to introduce the functions
(4.14)
On substitution of (4.8) and (4.10) into (3.39), and (4.12) into (3.41), and use of (4.14), it is found that
(4.15a)
 
(4.15b)
For determining the impedance it is the unknown boundary data on the hypotenuse |$S_2$| that is of interest. The first step is to eliminate |$\tilde{\rho }_3$| from equations (4.15):
(4.16)
Setting |$k\rightarrow -k$| in (4.16) yields
(4.17)
The function |$\tilde{\rho }_1$| can now be eliminated between (4.16) and (4.17):
(4.18)
where
(4.19)
Using (4.14) we can rewrite (4.18) as
(4.20)
A second equation is obtained by setting |$k\rightarrow \textrm{i}k$| in (4.20):
(4.21)
Finally, (4.21) can be both added to and subtracted from (4.20) to obtain the following equations:
(4.22a)
 
(4.22b)
Equations (4.22) can be solved for |$\tilde{\rho }_2(\overline{\alpha }k)\pm \tilde{\rho }_2(-\overline{\alpha }k)$|⁠:
(4.23a)
 
(4.23b)
Since |$\tilde{\rho }_2(k)$| is analytic for all |$k \ne 0$| the numerators in (4.23a) and (4.23b) must vanish when the denominators vanishes. This produces values of |$\tilde{\rho }(\alpha k)\pm \tilde{\rho }(-\alpha k)$| at special points in the |$k$|-plane:
(4.24a)
 
(4.24b)
where |$k_{1,m}\in \varSigma _1, k_{2,m}\in \varSigma _2$|⁠, and the sets |$\varSigma _j$| are defined by
(4.25a)
 
(4.25b)

4.4 Introducing a series expansion

The final step is to use (4.24) to recover the unknown boundary data |$\partial \phi /\partial \textbf{n}_2(z_2(s))$|⁠. One way to do this is to pose a series expansion of the form
(4.26)
where |$\{P_n(s)\}$| is some class of polynomials and |$\lbrace c_n \rbrace $| is a set of complex-valued coefficients. An over-determined linear system for these coefficients can then be found using (4.24), which can be solved numerically to find the unknown boundary data |$\partial \phi /\partial \textbf{n}_2(z_2(s))$|⁠. This approach is now a standard procedure in the analysis of the global relations in the Fokas method (Fokas, 2008; Fokas & Pelloni, 2014; Fokas & Spence, 2012). The supplementary materials contain details of this approach for a more general boundary value problem in the same triangular domain.
It turns out that, on making an appropriate choice of the series expansion (4.26), we can find explicit expressions for the series expansion coefficients. The relevant series expansion, known as a modified Fourier series (Iserles & Nørsett, 2008; Olver, 2009), can be expressed as
(4.27a)
 
(4.27b)
 
(4.27c)
Using (4.27), the left-hand side of (4.24) is
(4.28)
By substituting for |$G(k)$| and |$\hat{G}(k)$| in (4.19) using (4.9) and (4.13) it is found that
(4.29)
where
(4.30)
The Neumann data, or normal concentration flux, on the hypotenuse can therefore be expressed as
(4.31)
To find the impedance, we substitute for |$\phi $| and |$\partial \phi /\partial \textbf{n}_2$| in (4.2) using (4.1) and (4.31), respectively, producing the final explicit result
(4.32)
A few remarks are in order. While the modified Fourier series of the form (4.27) produces explicit expressions for the coefficients, greater accuracy is obtained by using the numerical approach based on (4.26) mentioned earlier, e.g. using Legendre or Chebyshev polynomials as the expansion functions. The supplementary materials contain a detailed comparison. Additionally, the modified Fourier series of the form (4.27) is only valid for functions that vanish at the boundary. In our case, our unknown flux satisfies this condition due to the nature of the boundary conditions. For more general problems we can split the function into two parts:
(4.33)
where |$ \check{\phi }(\pm{L}/{\sqrt{2}}) = 0$| and introduce a modified Fourier series expansion for |$\check{\phi }(s)$| instead.

4.5 Streamlines

To visualize the concentration distribution inside the domain the inverse transform (3.38a) can be used once the spectral functions have been obtained from the analysis of the global relations in the manner indicated above. Once |$\phi $| is computed in |${P}$|⁠, the concentration field |$C(x,y,t) $| at time |$t$| is
(4.34)
Figure 13 shows the concentration distribution |$C(x,y,0)$| for various values of |$\sigma $|⁠. As the magnitude of |$\sigma $| is increased, the concentration penetration depth decreases.
Contours of the concentration distribution $C(x,y,0)$ for various parameters $\sigma $. As $|\sigma |$ increases, the penetration depth gets smaller.
Fig. 13.

Contours of the concentration distribution |$C(x,y,0)$| for various parameters |$\sigma $|⁠. As |$|\sigma |$| increases, the penetration depth gets smaller.

Figure 14 and Figure 15 show contours of |$C(x,y,t)$| over one period of oscillation for |$\sigma = 2\textrm{j}$| and |$\sigma = 20\textrm{j}$|⁠. These figures show the influence of the penetration depth: in Fig. 14 where the penetration depth is larger the concentration distribution oscillates across the entire domain, while in Fig. 15 the concentration away from the stimulus wall is largely stationary.

Snapshots showing contours of the concentration distribution $C(x,y,t)$ over a period for $\sigma = 2\textrm{j}$.
Fig. 14.

Snapshots showing contours of the concentration distribution |$C(x,y,t)$| over a period for |$\sigma = 2\textrm{j}$|⁠.

Snapshots showing contours of the concentration distribution $C(x,y,t)$ over a period for $\sigma = 20\textrm{j}$.
Fig. 15.

Snapshots showing contours of the concentration distribution |$C(x,y,t)$| over a period for |$\sigma = 20\textrm{j}$|⁠.

4.6 Comparison with TauFactor

Cooper et al. (2017, 2016) have developed software called TauFactor on the MATLAB platform to compute an effective parameter called the tortuosity factor that quantifies the apparent decrease in diffusive transport resulting from convolutions of the flow paths through porous media. It can be used to compute the impedance for an arbitrary pore geometry, indeed, the software is designed to tackle materials with highly complex microstructure. Using a two colour map, where one colour represents a solid wall or obstacle, while the other represents pore space, a user can input a map into the software and it will compute the associated diffusion spectra, or impedance diagram. An analytical result such as (4.32) can serve as a numerical benchmark solution for software such as TauFactor, and vice versa. In order to verify (4.32) against TauFactor, we set |$D=1$| and |$L=1/\sqrt{2}$| giving
(4.35)
TauFactor uses a normalization of the impedance given by |$\tilde{Z}(\omega ) = ({\hat{A}D}/{\hat{L}})Z(\omega )$| where |$\hat{A}$| is the algebraic mean area accessible to diffusion, and |$\hat{L}$| is the maximum penetration distance from the stimulated surface to the longest pore path. For the right isosceles triangular path with unit length hypotenuse this normalization leads to |$\hat{A}D/\hat{L}=1$|⁠, therefore, the normalized impedance for comparison with TauFactor is simply given by (4.35). The results of this comparison are shown in Fig. 16 and show excellent agreement. Detailed information on the vertical asymptote shown in this figure can be found using a small-|$\omega $| perturbation analysis (Hauge, 2019), which reduces the problem to a set of coupled boundary value problems for Laplace’s equation where other transform methods in polygons are applicable (Crowdy, 2015a; Fokas & Kapaev, 2003).
Comparison of diffusion spectra (impedance) for a right isosceles triangular pore. Results of the TauFactor software compared with transform method result (4.35). The jagged nature of the TauFactor curve as $\omega \to 0$ curve is due to convergence issues in this limit [Cooper, private communication].
Fig. 16.

Comparison of diffusion spectra (impedance) for a right isosceles triangular pore. Results of the TauFactor software compared with transform method result (4.35). The jagged nature of the TauFactor curve as |$\omega \to 0$| curve is due to convergence issues in this limit [Cooper, private communication].

4.7 Application to the 3|$\omega $| method

Cahill & Pohl (1987) developed a method for measuring the thermal properties of substrates while minimizing the effects of black-body radiation; it has become known as the 3|$\omega $| method. The idea is to use an alternating current in a ‘heater’—used both as a heater to supply heat and as a thermocouple to measure temperature—to supply an oscillatory heat flux to a substrate, measuring the averaged temperature of the heater at various frequencies, and then using the measurements, along with a suitable solution to the governing PDE—i.e. the complex Helmholtz equation—to extract the thermal properties of the material through regression analysis. The name 3|$\omega $| arises because the method rests on measuring the amplitude of the third time harmonic induced owing to an assumed linear dependence of the heater resistance on temperature. Accurate thermal property measurement is crucial for a wide range of applications. In medicine, knowing the thermal properties of skin illuminates a patient’s hydration levels and overall dermatological health (Tian et al., 2017). The complex Helmholtz equation is at the heart of these methods: the original analysis of Cahill & Pohl (1987) uses a relatively simple solution in a half-plane, but the new techniques here provide a route to solving the same problem in much more geometrically complex geometries.

4.8 Application to diffusion wave field theory

Diffusion wave field theory covers many application areas where the diffusion equation forced by a time-periodic stimulus governs the system (Mandelis, 2000). A treatise (Mandelis, 2013) on diffusion-wave field problems focuses on this class of problems using Green’s functions as the principal tool, but the focus there is on 2D and 3D problems in infinite and semi-infinite domains. The transform method here is limited to 2D domains but applies to a very general class of polygonal geometries, and therefore adds a new tool to the armoury of available techniques in diffusion wave field theory. The new method is also based on Green’s functions, but the new ingredient is the use of a particular spectral representation, in terms of the complex variable |$k$|⁠, of the Green’s function that is relevant to convex polygons and the subsequent analysis of the global relations in this spectral |$k$| plane to obtain information on the solution.

It should be mentioned that the complex Helmholtz equation can also be used to study attenuated wave problems. There are important differences between diffusion waves and attenuated waves associated with the nature of the energy propagation mechanism (Scale & Snieder, 1999). On a related note, another possible use of our transform pair is investigation of the |$\textrm{arg}[\sigma ] \to \pi $| limit, especially in the case of unbounded domains, where the imaginary component of |$\sigma $| tends to zero and the complex Helmholtz equation reduces, in the limit, to the Helmholtz equation with its wave-like solutions.

4.9 Application to unsteady Stokes flow

An unsteady Stokes flow is one where the nonlinear inertial term in the fluid acceleration is ignored, but a time-derivative term retained. In such a case, one ends up studying the ‘linearized Navier–Stokes equations’. It is now shown how the new transform method for the complex Helmholtz operator opens up new possibilities for the analysis of such flows in polygonal domains brought about by a time-periodic stimulus. The so-called second Stokes problem is an example: it involves studying the motion of fluid in a half-plane due to an oscillatory semi-infinite wall. The fluid near the wall forms an oscillatory boundary layer, with thickness of the order of the penetration depth, called a Stokes layer (Batchelor, 1967). Unsteady Stokes flow problems in polygonal domains forced by time-periodic stimuli, or by fundamental singularities, have been treated previously using numerical methods (Chu & Kim, 2002; Clarke et al., 2005; Tuck, 1969), Fourier transforms (Chu & Kim, 2001; Tuck, 1970) and asymptotic analysis (Branicki & Moffatt, 2006; Clarke et al., 2005; Tuck, 1970).

With |$\hat{\textbf{x}}=(\hat{x},\hat{y})$|⁠, and |$\hat{\nabla }=(\partial _{\hat{x}},\partial _{\hat{y}})$|⁠, the Navier–Stokes equations for incompressible 2D Newtonian fluid flows are
(4.36)
where |$\hat{\textbf{u}}=(\hat{u},\hat{v})$| is the fluid velocity vector, |$\hat{p}$| is the fluid pressure and |$\rho $| and |$\mu $| denotes the fluid density and viscosity, respectively. On non-dimensionalization with
(4.37)
then (4.36) becomes
(4.38a)
where the dimensionless parameters are
(4.39)
where |$\nu =\mu /\rho $| is the kinematic viscosity. |$Re$| denotes the Reynolds number. If the time scale |$T$| is set by the periodic forcing of the system the parameter |$\beta $|⁠, often called a Strouhal number, controls the magnitude of inertial acceleration forces due to this time-dependent forcing relative to viscous forces.
Our interest here is in the unsteady Stokes flow equations, or the linearized Navier–Stokes equations, where |$Re\ll 1$| but where |$\beta $| is an order one quantity:
(4.40)
The curl of (4.40) leads to the vorticity equation
(4.41)
Due to the incompressibility condition (4.40), a stream function |$\psi ^{\prime}$| can be introduced where
(4.42)
If we now assume purely oscillatory solutions of the form
(4.43)
we end up with the system
(4.44)
In (4.44) both the 2D Laplacian operator and the complex Helmholtz operator appear. From the development in this paper ‘transform pairs’ for both operators in convex polygons are available opening up the possibility of a transform method to solving the unsteady Stokes equations in such domains. For unsteady Stokes flow, a layer of complication is added by the two operators being coupled: from (4.44) we obtain a fourth-order equation for |$\psi $| given by
(4.45)
This complication can be handled, however, just as transform methods for the fourth-order biharmonic equation can be formulated and solved (Crowdy & Fokas, 2004; Dimakos & Fokas, 2015; Luca & Crowdy, 2018) by viewing it as coupled problems for two harmonic fields.
We can write a solution to (4.45) as a superposition of solutions for each operator, namely,
(4.46)
where |$\phi $| satisfies Laplace’s equation, |$\omega $| is the vorticity satisfying (4.44) and |$b_j$| are arbitrary complex coefficients. In his analysis of the unsteady Stokes flow generated by an oscillating cylinder of arbitrary cross-section Tuck (1969) used a similar definition for |$\psi $| given, in dimensionless form, as
(4.47)
A benefit of this definition is that |$\phi $| becomes the harmonic conjugate to the pressure field |$p$|⁠, hence the pressure field follows from |$\phi $| using the Cauchy–Riemann equations
(4.48)

4.10 Transform method for unsteady Stokes flow

Let the polygonal domain be |${P}$| with sides |$\lbrace S_n | n=1,...,N \rbrace $| as before. From §2 we have the following integral representations for the derivatives of the harmonic field |$\phi $|⁠:
(4.49a)
 
(4.49b)
with
(4.50)
As for the vorticity, the following transform pair is found to be the most convenient:
(4.51a)
where |$\mathcal{L}_{\chi _n}$| is the contour defined earlier. Consequently, by taking (4.47) and substituting for |$\phi $| and |$\omega $| using (4.50) and (4.51a), we obtain an explicit integral representation for the stream function
(4.52)
The global relations are
(4.53a)
 
(4.53b)
It is also useful to define
(4.54)
satisfying
(4.55)
The approach to solving unsteady Stokes flow problem using the transform method derived herein is similar to the approach for the complex Helmholtz operator: the global relations (4.53) and the boundary conditions are analysed to compute the unknown boundary data. Once all the boundary data are known so are the spectral functions |$\eta _n(\alpha ), \hat{\eta }_n(\alpha )$| and |$\rho _n(k)$| and the integral representation (4.52) provides an expression for the stream function |$\psi $| inside the domain.

In the supplementary materials we use the new transform method just described to rederive a solution found by Tuck (1969) for source-driven motion into an upper half plane and then show how to generalize it to find a new solution involving a different boundary condition. This new method has also been used to retrieve the solutions of Chu & Kim (2001, 2002); those details can be found in Hauge (2019).

5. Discussion and future directions

This paper has presented a new transform pair to solve the complex Helmholtz equation in convex polygons. It can facilitate finding solutions to a wide class of boundary value problems in diffusion-wave field theory, electrochemical impedance spectroscopy and unsteady slow viscous flows in such domains. These are just a few of its potential applications.

The derivation of the new transform pair is a direct extension of one given by Crowdy (2015a,b) for analogous transform pairs for analytic functions in all domain classes illustrated in Fig. 3: any convex polygon, an arbitrary circular-arc domain, including multiply connected cases such as the important ‘disc-in-channel’ geometry. Such transform pairs can be used to solve boundary value problems for Laplace’s equation in all those classes of domain (Crowdy, 2015a,b). That work was extended to boundary value problems for biharmonic fields in Luca & Crowdy (2018). Crowdy’s extensions to circular domains were originally motivated by mathematical challenges in quantifying hydrodynamic slip on so-called superhydrophobic surfaces where circular domains are common; physically, this is due to the occurrence of free surface menisci satisfying a Laplace–Young force balance at small capillary numbers forcing them to have constant curvature. His methods have already proven to be of great utility in that application (Crowdy, 2015b; Luca et al., 2018; Schnitzer, 2016) and others (Crowdy, 2016b; Luca & Llewellyn Smith, 2018).

The most challenging part of the approach is the analysis of the global relations. This paper has given a flavour of what is involved, with more examples provided in the supplementary matter. There is, by now, a large and growing literature on analysing such global relations using a mixture of analytical and numerical techniques (Colbrook, 2020; Fokas, 2008; Fornberg & Flyer, 2011; Smitheman et al., 2010) so readers will find many other examples there. Readers are also referred to Hauge (2019) where a variety of other problems in EIS, diffusion wave field theory and the 3|$\omega $| method are treated, including how to add an inhomogeneous forcing to the right hand side of (1.1), which is quite a routine modification of what has been presented here.

As explained in the Introduction, while the term ‘global relations’ might be unfamiliar to readers not au fait with the Fokas method, these relations are closely connected to the established use of reciprocity theorems, or Lorentz reciprocity, in determining the effective properties of heterogeneous materials (Crowdy, 2017; Milton, 2002) and also useful in electrodynamics (Landau & Lifshitz, 1960; Panofsky & Phillips, 1962). We illustrated this connection, which does not appear to have been discussed previously in the literature, by carrying out a new derivation of simple Warburg impedance formulas using Green’s identity and special choices of a ‘comparison solution’ and, in particular, without ever solving for the concentration field whose effective impedance is being calculated. A well-known concept in electrical engineering is the impedance matrix, or Z-matrix, of a given system; these describe the behaviour of any linear electrical network viewed as a black box with |$N$| ports. In his derivation of ‘Fourier–Mellin’ transform pairs for |$N$|-sided polygons, or circular-arc polygons, Crowdy (2015a,b) introduces the notion of a spectral matrix with elements |$\lbrace \rho _{mn}(k) | m,n = 1, \ldots , N \rbrace $|⁠, which arises naturally in his approach; these were mentioned earlier in §2. In their studies of Laplace’s equation in a polygon Fokas & Kapaev (2003) did not introduce such a matrix, merely |$N$| spectral functions that, as shown later by Crowdy (2015a), turn out to be just the diagonal elements of the spectral matrix. In view of the connections made here between reciprocity theorems and analysis of the global relations there appears to be a tantalizing, but as yet unclear, connection between Crowdy’s spectral matrices and the concept of an impedance matrix. This matter is currently under investigation.

The complex Helmholtz equation was shown, in the Introduction, to arise naturally when studying time-periodic solutions of the 2D diffusion equation (1.2). In this case, initial transients are assumed to have decayed away. There has been much recent work (Miller & Smith, 2018; Pelloni & Smith, 2018; Sheils & Smith, 2015) on devising transform methods, again based on the Fokas method, for the 1D diffusion equation in various geometries in the unsteady case, i.e. for the solution of initial value problems. It would be interesting, to solve initial value problem for the applications discussed in this paper, to extend such approaches to the full 2D diffusion equation (1.2).

We end by reiterating two other open theoretical questions already posed earlier in this paper.

First: is there a transform pair representation for solutions of the complex Helmholtz equation—including the Helmholtz and modified Helmholtz equations—in general multiply connected circular domains? Spence & Fokas (2010) have contributed some results for the Helmholtz and modified Helmholtz equations in simple wedges, but is a more general geometric formulation available akin to the ‘Fourier–Mellin’ transform pairs in general circular domains known for the Laplace and biharmonic equation (Crowdy, 2015a,b; Luca & Crowdy, 2018)? The word ‘pair’ is emphasized here since, for theoretical interest, we also ask for an inverse transform representation, even though it is often not needed to formulate global relations that can be used to, say, solve for the Dirichlet-to-Neumann map (Colbrook, 2020; Luca et al., 2018).

Second: of theoretical interest is to ask how might the known ‘Fourier–Mellin’ transform pairs in (multiply connected) circular domains for harmonic and biharmonic fields (Crowdy, 2015a,b; Luca & Crowdy, 2018) be retrieved using a Lax-pair formulation akin to that used by Fokas & Kapaev (2003) in their original derivation of the transform pair for Laplace’s equation in a polygon? While the transform pair emerges naturally from Crowdy’s geometrical approach emulated here, much might be learned by formulating a more analytical approach via a Lax pair.

Acknowledgements

This paper is based on the Ph.D. thesis work of J.C.H. under the supervision of D.G.C. funded by the Engineering and Physical Sciences Research Council in the UK for which the financial support of a studentship is gratefully acknowledged.

References

Ablowitz
,
M. J.
and
Fokas
,
A. S.
(
2003
).
Complex Variables: Introduction and Applications
. Cambridge: Cambridge University Press.

Almond
,
D. P.
&
Patel
,
P.
(
1996
)
Photothermal Science and Techniques
.
Springer Science & Business Media
.

Babich
,
V.
,
Lyalinov
,
M.
&
Grikurov
,
V.
(
2008
)
Diffraction Theory. The Sommerfeld–Malyuzhinets Technique
, vol.
I
.
Alpha Science International
.

Bard
,
A. J.
&
Faulkner
,
L. R.
(
2001
)
Electrochemical Methods: Fundamentals and Applications
, 2nd edn.
Wiley
.

Barsoukov
,
E.
&
Macdonald
,
J. R.
(
2005
)
Impedance Spectroscopy: Theory, Experiment, and Applications
.
John Wiley & Sons, Inc
.

Batchelor
,
G. K.
(
1967
)
An Introduction to Fluid Dynamics
.
Cambridge University Press
.

Ben-Avraham
,
D.
&
Fokas
,
A. S.
(
2001
)
Solution of the modified Helmholtz equation in a triangular domain and an application to diffusion-limited coalescence
.
Phys. Rev. E
,
64
.

Bisquert
,
J.
&
Roberto
,
P.
(
1999
)
Theoretical models for ac impedance of finite diffusion layers exhibiting low frequency dispersion
.
J. Electroanal. Chem.
,
475
,
152
163
.

Branicki
,
M.
&
Moffatt
,
H. K.
(
2006
)
Evolving eddy structures in oscillatory Stokes flows in domains with sharp corners
.
J. Fluid Mech.
,
551
,
63
92
.

Cahill
,
D. G.
&
Pohl
,
R. O.
(
1987
)
Thermal conductivity of amorphous solids above the plateau
.
Phys. Rev. B
,
35
,
4067
4073
.

Chu
,
J. H.
&
Kim
,
M.
(
2001
)
Two-dimensional oscillatory Stokes flows between two parallel planes
.
Fluid Dyn. Res.
,
29
,
7
24
.

Chu
,
J. H.
&
Kim
,
M.
(
2002
)
Two-dimensional oscillatory Stokes flow in the region with a semi-infinite plate parallel to an infinite plane wall
.
Fluid Dyn. Res.
,
31
,
229
251
.

Clarke
,
R. J.
,
Cox
,
S. M.
,
Williams
,
P. M.
&
Jensen
,
O. E.
(
2005
)
The drag on a microcantilever oscillating near a wall
.
J. Fluid Mech.
,
545
,
397
426
.

Colbrook
,
M.
(
2020
)
Extending the unified transform: curvilinear polygons and variable coefficient PDEs
.
IMA J. Num. Anal.
,
40
,
976
1004
.

Cooper
,
S. J.
,
Bertei
,
A.
,
Finegan
,
D. P.
&
Brandon
,
N. P.
(
2017
)
Simulated impedance of diffusion in porous media
.
Electrochim. Acta
,
251
,
681
689
.

Cooper
,
S. J.
,
Bertei
,
A.
,
Shearing
,
P. R.
,
Kilner
,
J. A.
&
Brandon
,
N. P.
(
2016
)
TauFactor: an open-source application for calculating tortuosity factors from tomographic data
.
SoftwareX
,
5
,
203
210
.

Crowdy
,
D. G.
(
2011
)
Frictional slip lengths and blockage coefficients
.
Phys. Fluids
,
23
,
091703
.

Crowdy
,
D. G.
(
2015a
)
A transform method for Laplace’s equation in multiply connected circular domains
.
IMA J. Appl. Math.
,
80
,
1902
1931
.

Crowdy
,
D. G.
(
2015b
)
Fourier–Mellin transforms for circular domains
.
Comput. Methods Funct. Theory
,
15
,
655
687
.

Crowdy
,
D. G.
(
2016a
)
Geometry-fitted Fourier–Mellin transform pairs
.
Mathematical Analysis, Probability and Applications—Plenary Lectures
.
Springer
.

Crowdy
,
D. G.
(
2016b
)
Uniform flow past a periodic array of cylinders
.
Eur. J. Mech. B/Fluids
,
56
,
120
129
.

Crowdy
,
D. G.
(
2017
)
Perturbation analysis of subphase gas and meniscus curvature effects for longitudinal flows over superhydrophobic surfaces
.
J. Fluid Mech.
,
822
,
307
328
.

Crowdy
,
D. G.
&
Fokas
,
A. S.
(
2004
)
Explicit integral solutions for the plane elastostatic semi-strip
.
Proc. R. Soc. A
,
1285
1309
.

Dimakos
,
M.
&
Fokas
,
A. S.
(
2015
)
The Poisson and the biharmonic equations in the interior of a convex polygon
.
Stud. Appl. Math.
,
134
,
456
498
.

Fokas
,
A. S.
(
2008
)
A Unified Approach to Boundary Value Problems
, vol.
78
.
SIAM, CBMS-NSF Regional Conference Series
.

Fokas
,
A. S.
&
Kapaev
,
A. A.
(
2003
)
On a transform method for the Laplace equation in a polygon
.
IMA J. Appl. Math.
,
68
,
355
408
.

Fokas
,
A. S.
&
Pelloni
,
B.
(
2014
)
Unified Transform for Boundary Value Problems: Applications and Advances
.
SIAM
.

Fokas
,
A. S.
&
Pinotsis
,
D. A.
(
2006
)
The Dbar formalism for certain linear non-homogeneous elliptic PDEs in two dimensions
.
Eur. J. Appl. Math.
,
17
,
323
346
.

Fokas
,
A. S.
&
Spence
,
E. A.
(
2012
)
Synthesis, as opposed to separation, of variables
.
SIAM Rev.
,
54
,
291
324
.

Fornberg
,
B.
&
Flyer
,
N.
(
2011
)
A numerical implementation of Fokas boundary integral approach: Laplace’s equation on a polygonal domain
.
Proc. R. Soc. A
,
467
,
2983
3003
.

Hauge
,
J. C
. (
2019
).
The complex Helmholtz operator: a new transform approach
.
Ph.D. Thesis
.

Huang
,
J.
(
2018
)
Diffusion impedance of electroactive materials, electrolytic solutions and porous electrodes: Warburg impedance and beyond
.
Electrochim. Acta
,
281
,
170
188
.

Huang
,
J.
,
Li
,
Z.
,
Liaw
,
B. Y.
&
Zhang
,
J.
(
2016
)
Graphical analysis of electrochemical impedance spectroscopy data in Bode and Nyquist representations
.
J. Power Sources
,
309
,
82
98
.

Iserles
,
A.
&
Nørsett
,
S. P.
(
2008
)
From high oscillation to rapid approximation I: modified Fourier expansions
.
IMA J. Numer. Anal.
,
28
,
862
887
.

Keiser
,
H.
,
Beccu
,
K. D.
&
Gutjahr
,
M. A.
(
1976
)
Abschatzung der porenstruktur poroser elektroden aus impedanzmessungen
.
Electrochim. Acta
,
21
,
539
543
.

Keller
,
J.
(
1964
)
A theorem on the conductivity of a composite medium
.
J. Math. Phys.
,
5
,
548
549
.

Landau
,
L. D.
&
Lifshitz
,
E. M.
(
1960
)
Electrodynamics of Continuous Media
.
Addison-Wesley
.

Luca
,
E.
&
Crowdy
,
D. G.
(
2018
)
A transform method for the biharmonic equation in multiply connected circular domains
.
IMA J. Appl. Math.
,
83
,
942
976
.

Luca
,
E.
&
Llewellyn Smith
,
S.
(
2018
)
Stokes flow through a two-dimensional channel with a linear expansion
.
Q. J. Mech. Appl. Math
,
71
,
441
462
.

Luca
,
E.
,
Marshall
,
J. S.
&
Karamanis
,
G.
(
2018
)
Longitudinal shear flow over a bubble mattress with curved menisci: arbitrary protrusion angle and solid fraction
.
IMA J. Appl. Math.
,
83
,
917
941
.

MacDonald
,
D. D.
(
2006
)
Reflections on the history of electrochemical impedance spectroscopy
.
Electrochim. Acta
,
51
,
1376
1388
.

Mandelis
,
A.
(
2000
)
Diffusion waves and their uses
.
Phys. Today
,
53
,
29
.

Mandelis
,
A.
(
2013
)
Diffusion-Waves Fields: Mathematical Methods and Green Functions
.
Springer Science & Business Media
.

Miller
,
P. D.
&
Smith
,
D. A.
(
2018
)
The diffusion equation with nonlocal data
.
J. Math. Anal. Appl.
,
466
,
1119
1143
.

Milton
,
G.
(
2002
)
The Theory of Composites
.
Cambridge Univ. Press
.

Olver
,
S.
(
2009
)
On the convergence rate of a modified Fourier series
.
Math. Comp.
,
78
,
1629
1645
.

Panofsky
,
W. K. H.
&
Phillips
,
M.
(
1962
)
Classical Electricity and Magnetism
.
Addison-Wesley
.

Pelloni
,
B.
&
Smith
,
D. A.
(
2018
)
Nonlocal and multipoint boundary value problems for linear evolution equations
.
Stud. Appl. Math.
,
141
,
46
88
.

Rosencwaig
,
A.
&
Gersho
,
A.
(
1976
)
Theory of photoacoustic effect with solids
.
J. Appl. Phys.
,
47
.

Scale
,
J. A.
&
Snieder
,
R.
(
1999
)
What is a wave
?
Nature
,
401
,
739
.

Schnitzer
,
O.
(
2016
)
Singular effective slip length for longitudinal flow over a dense bubble mattress
.
Q. J. Mech. Appl. Math
,
1
,
052101(R)
.

Shah
,
R. K.
&
London
,
A. L.
(
1978
)
Laminar Flow Forced Convection in Ducts
.
Academic Press
.

Sheils
,
N. E.
&
Smith
,
D. A.
(
2015
)
Heat equation on a network using the Fokas method
.
J. Phys. A: Math. Theor.
,
48
,
335001
.

Smitheman
,
S. A.
,
Spence
,
E. A.
&
Fokas
,
A. S.
(
2010
)
A spectral collocation method for the Laplace and modified Helmholtz equations in a convex polygon
.
IMA J. Numer. Anal.
,
30
,
1184
1205
.

Spence
,
E. A.
&
Fokas
,
A. S.
(
2010
)
A new transform method I: domain-dependent fundamental solutions and integral representations
.
Proc. R. Soc. A
,
466
,
2259
2281
.

Stakgold
,
I.
(
2000
)
Boundary Value Problems in Mathematical Physics
, vol.
I
.
Society for Industrial and Applied Mathematics
.

Tian
,
L.
,
Li
,
Y.
,
Webb
,
R. C.
,
Krishnan
,
S.
,
Bian
,
Z.
,
Song
,
J.
,
Ning
,
X.
,
Crawford
,
K.
,
Kurniawan
,
J.
,
Bonifas
,
A.
,
Ma
,
J.
,
Liu
,
Y.
,
Xie
,
X.
,
Chen
,
J.
,
Liu
,
Y.
,
Shi
,
Z.
,
Wu
,
T.
,
Ning
,
R.
,
Li
,
D.
,
Sinha
,
S.
,
Cahill
,
D. G.
,
Huang
,
Y.
&
Rogers
,
J. A.
(
2017
)
Flexible and stretchable 3|$\omega $| sensors for thermal characterization of human skin
.
Adv. Funct. Mater.
,
1701282
,
1
9
.

Tuck
,
E. O.
(
1969
)
Calculation of unsteady flows due to small motions of cylinder in a viscous fluid
.
J. Engrg. Math.
,
3
.

Tuck
,
E. O.
(
1970
)
Unsteady flow of a viscous fluid from a source in a wall
.
J. Fluid Mech.
,
41
,
641
652
.

Warburg
,
E.
(
1899
)
Ueber das Verhalten sogenannter unpolarisirbarer Elektroden gegen Wechselstrom
.
Ann. Phys.
,
303
,
493
499
.

Watson
,
G.
(
1995
)
A Treatise on the Theory of Bessel Functions
.
Cambridge Univ. Press
.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Supplementary data