*Article* **The Modified Helmholtz Equation on a Regular Hexagon—The Symmetric Dirichlet Problem**

#### **Konstantinos Kalimeris 1,\* and Athanassios S. Fokas 1,2,3**


Received: 9 June 2020; Accepted: 27 June 2020; Published: 28 July 2020

**Abstract:** Using the unified transform, also known as the Fokas method, we analyse the modified Helmholtz equation in the regular hexagon with symmetric Dirichlet boundary conditions; namely, the boundary value problem where the trace of the solution is given by the same function on each side of the hexagon. We show that if this function is odd, then this problem can be solved in closed form; numerical verification is also provided.

**Keywords:** unified transform; modified Helmholtz equation; global relation

#### **1. Introduction**

We analyse the modified Helmholtz equation in a regular hexagon using the unified transform, also known as the Fokas method. This method was introduced by one of the authors [1], for analysing integrable nonlinear partial differential equations (PDEs) [2]. Later, it was realized that it also yields novel results for linear evolution PDEs [3]; results in this direction are obtained by several authors [4–10]. Furthermore, it yields new integral representations for the solution of linear elliptic PDEs in polygonal domains [11], which in the case of simple domains can be used to obtain the analytical solution of several problems which apparently cannot be solved by the standard methods [12,13]. Recently, researchers utilised the integral representations provided by the Fokas method for the local and global wellposedness analysis of Korteweg-de Vries and nonlinear Schrödinger type PDEs [14–18], as well as for studying problems from control theory [19].

The Fokas method is based on two basic ingredients:


For linear PDEs, the Fokas method involves the following:


• The explicit solution of the problem is derived by determining the contribution of the unknown boundary values to the integral representation. This can be achieved by using the global relation, as well as equations obtained from the global relation through certain invariant transformations.

The global relation has had important analytical and numerical implications: first, it has led to novel analytical formulations of a variety of important physical problems from water waves [20–26] to three-dimensional layer scattering [27]. Second, it has led to the development of new techniques for the Laplace, modified Helmholtz, Helmholtz, biharmonic equations, both analytical [28–35] and numerical [36–47].

The above analytical solutions are given in terms of infinite series; this is to be contrasted to other techniques based on the eigenvalues of the Laplace operator that yield the solution as a bi-infinite series. The eigenvalues of the Laplace operator for the Dirichlet, Neumann and Robin problems in the interior of an equilateral triangle were first obtained by Lamé in 1833 [48]; these results have also been derived using the Fokas method [49]. Completeness for the associated expansions for the Dirichlet and Neumann problems was obtained in [50–53] using group theoretic techniques. McCartin rederived these results [54,55] and studied the connection of the eigen-structure of the equilateral triangle with that of the regular hexagon [56]. The above remarks indicate that the existing literature is based on an implicit way for deriving the solution of specific BVPs of the regular hexagon in terms of bi-infinite series. This is to be contrasted with our work which presents a direct approach for deriving explicit integral representations of the solution of a special BVP on the regular hexagon; the extension of the current methodology to more general problems is under investigation.

#### *Organisation of the Paper*

In Section 2 we implement the four steps discussed above for solving the symmetric Dirichlet problem of the modified Helmholtz equation in a regular hexagon. The main achievement of this work is presented in Section 3 and concerns the fourth step: our analysis yields the solution for the case of odd symmetric Dirichlet data in the closed form (34). We study the case of even symmetric data in Section 4, where we derive the expression (37); this expression in addition to known terms also involves an unknown term. In Section 5, Figures 1 and 2 depict the numerical verification of the main result of Section 3; also, Figures 7 and 8 indicate that the unknown term in the expression (37) is exponentially small in the high frequency limit, and hence this result provides an excellent approximation for this physically significant limit.

#### **2. The Basic Elements**

The equation investigated here is the modified Helmholtz equation in the interior of the regular hexagon, *D*, namely,

$$q\_{xx} + q\_{yy} - 4\beta^2 q = 0, \ (\mathbf{x}, y) \in D,\tag{1}$$

where *q*(*x*, *y*) is a real valued function and *β* > 0.

Using complex coordinates,

$$z = \mathfrak{x} + i\mathfrak{y}\_{\prime} \qquad \mathfrak{z} = \mathfrak{x} - i\mathfrak{y}\_{\prime}$$

Equation (1) becomes

$$
\rho\_{zz} - \beta^2 q = 0.\tag{2}
$$

*2.1. The Global Relation and the Integral Representation of the Solution in the Interior of a Convex Polygon*

We first derive the global relation:

The formal adjoint also satisfies the modified Helmholtz equation

$$
\ddot{q}\_{\mathbb{Z}^2} - \beta^2 \ddot{q} = 0. \tag{3}
$$

Multiplying Equation (2) by *q*˜, Equation (3) by *q* and subtracting, we find

$$
\dagger q\_{\overline{z}\overline{z}} - q \overline{q}\_{\overline{z}\overline{z}} = 0,
\tag{4}
$$

or equivalently

$$\frac{\partial}{\partial z} \left( \vec{q} q\_{\vec{z}} - \vec{q}\_{\vec{z}} q \right) + \frac{\partial}{\partial \vec{z}} \left( q \vec{q}\_{\vec{z}} - q\_{\vec{z}} \vec{q} \right) = 0. \tag{5}$$

Using in (5) the special solution *q*˜ = *e* <sup>−</sup>*iβ*(*kz*<sup>−</sup> *<sup>z</sup>*¯ *<sup>k</sup>* ) and employing Green's theorem, we obtain

$$\int\_{\partial\Omega} W(z,\bar{z},k) = 0, \qquad k \in \mathbb{C}\_{\prime} \tag{6}$$

where *W* is defined by

$$\mathcal{W}(z,\bar{z},k) = \varepsilon^{-i\beta\left(kz-\frac{\pi}{k}\right)} \left[ \left( q\_z + ik\beta q \right) dz - \left( q\_{\bar{z}} + \frac{\beta}{i\bar{k}} q \right) d\bar{z} \right], \qquad k \in \mathbb{C}. \tag{7}$$

Suppose that Ω is the polygon defined via the points *z*1, *z*2, ... , *zn*, *zn*+<sup>1</sup> = *z*1. Then (6) gives the following global relation for the modified Helmholtz in this polygon:

$$\sum\_{j=1}^{n} \mathfrak{q}\_{j}(k) = 0, \qquad k \in \mathbb{C}, \tag{8}$$

where % *q*ˆ*j*(*k*) &*n* <sup>1</sup> are defined by

$$\hat{q}\_{j}(k) = \int\_{z\_{j}}^{z\_{j+1}} e^{-i\beta\left(kz - \frac{\alpha}{k}\right)} \left[ \left( q\_{\bar{z}} + ik\beta q \right) dz - \left( q\_{\bar{z}} + \frac{\beta}{i\bar{k}} q \right) d\bar{z} \right], \ k \in \mathbb{C}, \tag{9}$$

or equivalently (in local coordinates) by

$$\mathfrak{d}\_{\vec{j}}(k) = \int\_{z\_{\vec{j}}}^{z\_{\vec{j}+1}} e^{-i\beta(kz - \frac{\pi}{k})} \left[ iq\_{N}^{(j)}(s) + i\beta \left( \frac{1}{k} \frac{d\vec{z}}{ds} + k \frac{dz}{ds} \right) q^{(j)}(s) \right] ds, \ k \in \mathbb{C},$$

$$j = 1, \ldots, n. \tag{10}$$

In Equation (10) we have used the identity

$$q\_z dz - q\_{\overline{z}} d\overline{z} = i q\_N ds\_{\prime}$$

where *s* is the arclength on the boundary *z*(*s*) = *x*(*s*) + *iy*(*s*) of the polygon and *qN* denotes the derivative in the outward normal direction to the boundary of the polygon.

In order to derive the integral representation of the solution one has to implement the spectral analysis of the differential form

$$d\left[e^{-i\beta\left(kz-\frac{\bar{z}}{k}\right)}\mu(z,k)\right] = \mathcal{W}(z,\bar{z},k), \qquad k \in \mathbb{C}.\tag{11}$$

This procedure yields the following theorem, proven in [6]:

**Theorem 1.** *Let* Ω *be the interior of a convex closed polygon in the complex z-plane, with corners z*1, ... , *zn*, *zn*+<sup>1</sup> ≡ *z*1*. Assume that there exists a solution q*(*z*, *z*¯) *of the modified Helmholtz equation, i.e., of Equation (2) with β* > 0*, valid on* Ω*, and suppose that this solution has sufficient smoothness on the boundary of the polygon.*

*Then, q can be expressed in the form*

$$q(z, \bar{z}) = \frac{1}{4\pi i} \sum\_{j=1}^{n} \int\_{l\_j} e^{j\beta \left(kz - \frac{\pi}{k}\right)} \hat{q}\_j(k) \frac{dk}{k},\tag{12}$$

*where* {*q*ˆ*j*(*k*)}*<sup>n</sup>* <sup>1</sup> *are defined by* (10)*, and* {*lj*}*<sup>n</sup>* <sup>1</sup> *are the rays in the complex k-plane*

$$l\_j = \{ k \in \mathbb{C} : \arg k = -\arg(z\_{j+1} - z\_j) \}, \qquad j = 1, \dots, n$$

*oriented from zero to infinity.*

Observe that the solution given in (12) is given in terms of {*q*ˆ*j*}*<sup>n</sup>* <sup>1</sup> which involve integral transforms of both *q* and *qN* on the boundary, i.e., both known and unknown functions.

#### *2.2. The Dirichlet Problem on a Regular Hexagon*

Let *<sup>D</sup>* <sup>⊂</sup> <sup>C</sup> be the interior of a regular hexagon with vertices {*zj*}<sup>6</sup> 1,

$$z\_1 = \frac{l\sqrt{3}}{2} - i\frac{l}{2} = l e^{\frac{-j\pi}{6}} \qquad \text{and} \qquad z\_j = \omega^{j-1} z\_{1\prime} \tag{13}$$

where *l* is the length of the side and *ω* = *e iπ* <sup>3</sup> . The sides {(*zj*, *zj*+1)}<sup>6</sup> <sup>1</sup>, *z*<sup>7</sup> ≡ *z*<sup>1</sup> will be referred to as sides {(*j*)}<sup>6</sup> 1.

For the sides {(*j*)}<sup>6</sup> <sup>1</sup> the following parametrizations will be used:

$$z\_1(\mathbf{s}) = \frac{l\sqrt{3}}{2} + i\mathbf{s}, \quad z\_j(\mathbf{s}) = \left(\frac{l\sqrt{3}}{2} + i\mathbf{s}\right)\omega^{j-1}, \qquad \mathbf{s} \in \left[-\frac{l}{2}, \frac{l}{2}\right].$$

The general Dirichlet problem can be uniquely decomposed to 6 simpler Dirichlet problems, by employing the decomposition

$$q^{(j)}(s) = \sum\_{i=1}^{6} \omega^{(j-1)(i-1)} g\_i(s), \qquad j = 1, \dots, 6, \qquad s \in \left[ -\frac{l}{2}, \frac{l}{2} \right];$$

indeed the determinant of the matrix ! *ω*(*j*−1)(*i*−1) *i*,*j*=1,...,6" is non-zero (Its value is 216 = 63, and for the general case Det ! *ω*(*j*−1)(*i*−1) *i*,*j*=1,...,*n* " = *i* 2−*n*(*n*+1) <sup>2</sup> *n<sup>n</sup>*/2).

The existence and uniqueness of the solution of the modified Helmholtz equation shows that it is sufficient to solve each one of the above Dirichlet problems. The first of them is the symmetric Dirichlet problem, where the value *g*1(*s*) = *d*(*s*) is prescribed on each side. This symmetric problem is analysed in the next section.

#### *2.3. The Symmetric Dirichlet Problem*

The problem analysed in this subsection is the symmetric Dirichlet problem for the modified Helmholtz equation in the regular hexagon (Ω ≡ *D*). Let *d*(*s*) be a real function with sufficient smoothness and compatibility at the vertices of the hexagon, i.e., *d l* 2 = *d* − *l* 2 . We prescribe the boundary conditions

$$q^{(j)}(s) = d(s), \ s \in \left[ -\frac{l}{2}, \frac{l}{2} \right], \ j = 1, \ldots, 6.$$

The above 'symmetry' property also holds for the Neumann boundary values. This fact is the consequence of the following three observations:


$$\frac{\partial q(z)}{\partial z} dz = \frac{\partial q(\omega z)}{\partial z} dz = \frac{\partial(\omega z)}{\partial z} \frac{\partial q(\omega z)}{\partial(\omega z)} \frac{1}{\omega} d(\omega z) = \frac{\partial q(\omega z)}{\partial(\omega z)} d(\omega z).$$

• Evaluating the above differential form on each side we obtain

$$q\_z dz = \frac{1}{2} \left( \dot{q}^{(j)}(s) + i q\_N^{(j)}(s) \right) ds = \frac{1}{2} \left( d'(s) + i q\_N^{(j)}(s) \right) ds\_\prime$$

where the second equality is a direct consequence of the fact that the Dirichlet data are invariant under this rotation.

Thus,

$$q\_N^{(j)}(s) = \mu(s), \ s \in \left[ -\frac{l}{2}, \frac{l}{2} \right], \ j = 1, \ldots, 6.$$

Applying the parametrization of the regular hexagon on Equation (10) we obtain:

$$\mathfrak{q}\_1(k) = \mathfrak{q}(k), \qquad \mathfrak{q}\_j(k) = \mathfrak{q}\left(\omega^{j-1}k\right), \qquad j = 1, \ldots, 6,\tag{14}$$

with

$$
\hat{q}(k) = E(-ik)[i\mathcal{U}(k) + D(k)],\tag{15}
$$

where *E*(*k*), *D*(*k*) and *U*(*k*) are defined by

$$\begin{aligned} E(k) &= e^{\mathcal{S}(k+\frac{1}{k})\frac{l\sqrt{d}}{2}}, \\ D(k) &= \mathcal{S}\left(\frac{1}{k} - k\right) \int\_{-\frac{1}{2}}^{\frac{l}{2}} e^{\mathcal{S}(k+\frac{1}{k})s} d(s) ds, \\ \mathcal{U}(k) &= \int\_{-\frac{1}{2}}^{\frac{l}{2}} e^{\mathcal{S}(k+\frac{1}{k})s} u(s) ds, \quad &k \in \mathbb{C}. \end{aligned} \tag{16}$$

The function *D*(*k*) is known, whereas the unknown function *U*(*k*) contains the unknown Neumann boundary value *u*(*s*) = *qN*.

Using (15), the global relation (8) takes the form

$$\begin{aligned} E(-ik)\mathcal{U}(k) + E(-i\omega k)\mathcal{U}(\omega k) + E(-i\omega^2 k)\mathcal{U}(\omega^2 k) \\ + E(ik)\mathcal{U}(-k) + E(i\omega k)\mathcal{U}(-\omega k) + E(i\omega^2 k)\mathcal{U}(-\omega^2 k) = i\mathcal{G}(k), \quad k \in \mathbb{C}, \end{aligned} \tag{17}$$

where the known function *G*(*k*) is defined by

$$G(k) = \sum\_{j=1}^{6} E\left(-i\omega^{j-1}k\right)D\left(\omega^{j-1}k\right), \qquad k \in \mathbb{C}.\tag{18}$$

The integral representation (12) of the solution takes the form

$$q(z,\overline{z}) = \frac{1}{4\pi i} \sum\_{j=1}^{6} \int\_{l\_j} e^{i\beta\left(kz - \frac{z}{k}\right)} \mathbb{E}\left(-i\omega^{j-1}k\right) \left[\mathcal{D}\left(\omega^{j-1}k\right) + i\mathcal{U}\left(\omega^{j-1}k\right)\right] \frac{dk}{k},\tag{19}$$

where {*lj*}<sup>6</sup> <sup>1</sup> are the rays defined by

$$l\_j = \left\{ k \in \mathbb{C} : \arg k = \frac{11 - 2j}{6}\pi \right\}, \qquad j = 1, \ldots, 6,\tag{20}$$

oriented from zero to infinity. The principal arguments of {*l*1, *<sup>l</sup>*2, *<sup>l</sup>*3, *<sup>l</sup>*4, *<sup>l</sup>*5, *<sup>l</sup>*6} are <sup>3</sup>*<sup>π</sup>* 2 , 7*π* 6 , 5*π* 6 , *π* 2 , *π* 6 , 11*π* 6 , respectively.

Since the function *d*(*s*) can be uniquely written as a sum of an odd and an even function, we will only consider two particular cases:


The solution and the Neumann boundary values inherit the analogous properties:


#### **3. Derivation of the Solution for the Symmetric Odd Case**

In what follows we will show that the contribution of the unknown functions % *U ωj*−1*k* &<sup>6</sup> <sup>1</sup> to the solution representation (19) can be computed explicitly.

Applying the condition *U*(−*k*) = −*U*(*k*) in (17) we obtain the equation

$$
\Delta(ik)\mathcal{U}(k) + \Delta(i\omega k)\mathcal{U}(\omega k) + \Delta(i\omega^2 k)\mathcal{U}(\omega^2 k) = -i\mathcal{G}(k), \qquad k \in \mathbb{C} \tag{21}
$$

where *G*(*k*) is given in (18) and Δ(*k*) is defined by

$$
\Delta(k) = E(k) - E(-k).
$$

Solving (21) for *U*(*k*) and substituting the resulting expression in (15) we find

$$q(k) = E(-ik)D(k) + \frac{E(-ik)G(k)}{\Delta(ik)}$$

$$+ i[E(-ik)E(-i\omega k) - E(-ik)E(i\omega k)]\frac{l\varPi(\omega k)}{\Delta(ik)}$$

$$+ i[E(-ik)E(-i\omega^2 k) - E(-ik)E(i\omega^2 k)]\frac{l\varPi(\omega^2 k)}{\Delta(ik)}.\tag{22}$$

The functions *q*ˆ*j*(*k*) can be obtained from (22) by replacing *k* with *ωj*−1*k* for *j* = 1, . . . , 6.

Regarding the integral representation of the solution, we restrict our attention to the first integral of (19), namely the integral along *l*<sup>1</sup> (the negative imaginary axis).

Let

$$\mathcal{P} = e^{i\beta \left(kz - \frac{\pi}{2}\right)\_\cdot}$$

Solving (21) for *U*(*k*) and substituting the resulting expression in the first integral of (19) we find that the known part of this integral is given by the expression

$$F\_1 = \frac{1}{4\pi i} \int\_{I\_1} \mathcal{P}E(-ik) \left[ D(k) + \frac{G(k)}{\Delta(ik)} \right] \frac{dk}{k}. \tag{23}$$

The unknown part involves the functions *U*(*ωk*) and *U*(*ω*2*k*) and is given by

$$\begin{split} \mathbb{C}\_{1} &= \frac{1}{4\pi} \int\_{I\_{1}} \mathcal{P}\left[E(-ik)E(-i\omega k)\frac{\mathcal{U}(\omega k)}{\Delta(ik)} + E(-ik)E(-i\omega^{2}k)\frac{\mathcal{U}(\omega^{2}k)}{\Delta(ik)}\right] \frac{dk}{k} \\ &- \frac{1}{4\pi} \int\_{I\_{1}} \mathcal{P}\left[E(-ik)E(i\omega k)\frac{\mathcal{U}(\omega k)}{\Delta(ik)} + E(-ik)E(i\omega^{2}k)\frac{\mathcal{U}(\omega^{2}k)}{\Delta(ik)}\right] \frac{dk}{k} .\end{split}$$

In what follows we will show that the contribution of the unknown functions, namely of the sum ∑6 <sup>1</sup> *Cj*, can be computed in terms of the given boundary conditions.

The first integral in the rhs of *C*<sup>1</sup> can be deformed from *l*<sup>1</sup> to *l* <sup>1</sup>, where *l* <sup>1</sup> is a ray with <sup>7</sup>*<sup>π</sup>* <sup>6</sup> <sup>≤</sup> arg *<sup>k</sup>* <sup>≤</sup> <sup>3</sup>*<sup>π</sup>* <sup>2</sup> ; choosing *l* <sup>1</sup> ≡ *l*<sup>2</sup> we obtain

$$\mathbf{C}\_{1} = \mathbf{\hat{C}}\_{1} + \mathbf{\hat{C}}\_{1} \,\tag{24}$$

where

$$\hat{\mathcal{C}}\_{1} = \frac{1}{4\pi} \int\_{l\_{2}} \mathcal{P}\left[E(-ik)E(-i\omega k)\frac{\mathcal{U}(\omega k)}{\Delta(ik)} + E(-ik)E(-i\omega^{2}k)\frac{\mathcal{U}(\omega^{2}k)}{\Delta(ik)}\right] \frac{dk}{k}$$

and

$$\mathcal{L}\_1 = -\frac{1}{4\pi} \int\_{l\_1} \mathcal{P}\left[E(-ik)E(i\omega k)\frac{\mathcal{U}(\omega k)}{\Delta(ik)} + E(-ik)E(i\omega^2 k)\frac{\mathcal{U}(\omega^2 k)}{\Delta(ik)}\right] \frac{dk}{k}.$$

The above deformation is justified, since it can be shown that the integrand of *C*ˆ <sup>1</sup> is bounded and analytic in the region where arg *<sup>k</sup>* <sup>∈</sup> <sup>7</sup>*<sup>π</sup>* <sup>6</sup> , <sup>3</sup>*<sup>π</sup>* 2 : letting *a* = *e<sup>i</sup> <sup>π</sup>* <sup>6</sup> , we can rewrite the first term of the integrand of *C*ˆ <sup>1</sup> in the form

$$\mathcal{P}E^{-\frac{2}{\sqrt{5}}}(iak)\frac{E^{\frac{2}{\sqrt{5}}}(iak)E(-ik)E(-i\omega k)E^{\frac{1}{\sqrt{5}}}(\omega k)}{\Delta(ik)}E^{-\frac{1}{\sqrt{5}}}(\omega k)\mathcal{U}(\omega k).$$

We observe the following:


Indeed, if *<sup>z</sup>* <sup>∈</sup> *<sup>D</sup>*, then <sup>5</sup>*<sup>π</sup>* <sup>≤</sup> arg(*<sup>z</sup>* <sup>−</sup> *<sup>z</sup>*2) <sup>≤</sup> <sup>3</sup>*<sup>π</sup>* . Thus, if <sup>7</sup>*<sup>π</sup>* <sup>≤</sup> arg *<sup>k</sup>* <sup>≤</sup> <sup>3</sup>*<sup>π</sup>* , it follows that *π* ≤ arg[*k*(*z* − *z*2)] ≤ 3*π*. Hence, Re% *ik*(*z* − *z*2) & ≤ 0.

Therefore, the exponentials *eiβk*(*z*−*z*2) and *e β ik* (*z*¯−*z*¯2) are bounded.

• The function *<sup>E</sup>*<sup>−</sup> <sup>√</sup><sup>1</sup> <sup>3</sup> (*ωk*)*U*(*ωk*) is bounded and analytic for arg *<sup>k</sup>* <sup>∈</sup> [ <sup>7</sup>*<sup>π</sup>* <sup>6</sup> , <sup>13</sup>*<sup>π</sup>* <sup>6</sup> ], namely in the region where Re(*ωk*) ≥ 0.

Indeed, this expression involves the exponentials *eβωk*(*s*<sup>−</sup> *<sup>l</sup>* <sup>2</sup> ) and *e<sup>β</sup>* <sup>1</sup> *<sup>ω</sup><sup>k</sup>* (*s*<sup>−</sup> *<sup>l</sup>* 2 ) , which are bounded in this region, since *<sup>s</sup>* <sup>≤</sup> *<sup>l</sup>* 2 .

• The function

$$\frac{E^{\frac{2}{\sqrt{3}}}(iak)E(-ik)E(-i\omega k)E^{\frac{1}{\sqrt{3}}}(\omega k)}{\Delta(ik)} = \frac{E^{\frac{1}{\sqrt{3}}}(k)}{\Delta(ik)}\lambda$$

is bounded and analytic for arg *<sup>k</sup>* <sup>∈</sup> [ <sup>7</sup>*<sup>π</sup>* <sup>6</sup> , <sup>3</sup>*<sup>π</sup>* 2 ].

Indeed, since *k* is at the lower half plane, then

$$\frac{E^{\frac{1}{\sqrt{3}}}(k)}{\Delta(ik)} \sim \frac{E^{\frac{1}{\sqrt{3}}}(k)}{E(ik)} = E^{-\frac{2}{\sqrt{3}}}(\omega^2 k), \qquad k \to \infty,$$

which is bounded if Re *ω*2*k* ≥ 0.

If arg *<sup>k</sup>* <sup>∈</sup> [ <sup>7</sup>*<sup>π</sup>* <sup>6</sup> , <sup>3</sup>*<sup>π</sup>* <sup>2</sup> ], then arg *ω*2*k* <sup>∈</sup> [ <sup>11</sup>*<sup>π</sup>* <sup>6</sup> , <sup>13</sup>*<sup>π</sup>* <sup>6</sup> ], which yields Re *ω*2*k* > 0.

Similar considerations apply to the second term of the integrand of *C*ˆ 1; this term can be rewritten in the form

$$\mathcal{P}E^{-\frac{2}{\sqrt{3}}}(iak)\frac{E^{\frac{2}{\sqrt{3}}}(iak)E(-ik)E(-i\omega^2k)E^{\frac{1}{\sqrt{3}}}(\omega^2k)}{\Delta(ik)}E^{-\frac{1}{\sqrt{3}}}(\omega^2k)\mathcal{U}(\omega^2k).$$

We observe the following:


$$\frac{E^{\frac{2}{\sqrt{3}}}(iak)E(-ik)E(-i\omega^2k)E^{\frac{1}{\sqrt{3}}}(\omega^2k)}{\Delta(ik)} \sim 1, \qquad k \to \infty.$$

Thus, it is bounded and analytic for arg *<sup>k</sup>* <sup>∈</sup> [ <sup>7</sup>*<sup>π</sup>* <sup>6</sup> , <sup>3</sup>*<sup>π</sup>* 2 ].

Using the underlined symmetries, we can express the integral representation of the solution in the form

$$q = \sum\_{j=1}^{6} F\_j + \sum\_{j=1}^{6} \mathbb{C}\_j = \sum\_{j=1}^{6} F\_j + \sum\_{j=1}^{6} \left(\mathbb{C}\_j + \mathbb{C}\_j\right),\tag{25}$$

where *Fj* and *Cj* are given by applying in (23) and (24) the following rotations:

*<sup>k</sup>* <sup>→</sup> *<sup>ω</sup>j*−1*k*, *<sup>l</sup>*<sup>1</sup> <sup>→</sup> *lj*, *<sup>l</sup>*<sup>2</sup> <sup>→</sup> *lj*+1, *<sup>j</sup>* <sup>=</sup> 2, . . . , 6; *<sup>l</sup>*<sup>7</sup> :<sup>=</sup> *<sup>l</sup>*1.

We define *<sup>C</sup>*,*<sup>j</sup>* <sup>=</sup> *<sup>C</sup>*<sup>ˆ</sup> *<sup>j</sup>*−<sup>1</sup> <sup>+</sup> *<sup>C</sup>*<sup>ˇ</sup> *<sup>j</sup>*, *j* = 1, ... , 6, where we employ the notation *C*ˆ <sup>0</sup> = *C*ˆ 6. Then, we rewrite the expression in (25) in the form

$$q = \sum\_{j=1}^{6} F\_j + \sum\_{j=0}^{5} \mathcal{C}\_j + \sum\_{j=1}^{6} \mathcal{C}\_j = \sum\_{j=1}^{6} F\_j + \sum\_{j=1}^{6} \left(\mathcal{C}\_{j-1} + \mathcal{C}\_j\right) = \sum\_{j=1}^{6} F\_j + \sum\_{j=1}^{6} \mathcal{C}\_j. \tag{26}$$

Thus, it is sufficient to compute the contribution {*C*,*j*}<sup>6</sup> <sup>1</sup>. In this direction we find (via rotation) that

$$\mathcal{L}\_2 = -\frac{1}{4\pi} \int\_{l\_2} \mathcal{P}\left[E(-i\omega k)E(i\omega^2 k)\frac{\mathcal{U}(\omega^2 k)}{\Delta(i\omega k)} + E(-i\omega k)E(i\omega^3 k)\frac{\mathcal{U}(\omega^3 k)}{\Delta(i\omega k)}\right] \frac{dk}{k}.$$

Thus

$$\begin{split} \tilde{C}\_{2} &= \tilde{C}\_{1} + \tilde{C}\_{2} \\ &= \frac{1}{4\pi} \int\_{l\_{2}} \mathcal{P}\left[E(-ik)E(-i\omega k)\frac{\mathcal{U}(\omega k)}{\Delta(ik)} + E(-ik)E(-i\omega^{2}k)\frac{\mathcal{U}(\omega^{2}k)}{\Delta(ik)}\right] \frac{dk}{k} \\ &- \frac{1}{4\pi} \int\_{l\_{2}} \mathcal{P}\left[E(-i\omega k)E(i\omega^{2}k)\frac{\mathcal{U}(\omega^{2}k)}{\Delta(i\omega k)} + E(-i\omega k)E(i\omega^{3}k)\frac{\mathcal{U}(\omega^{3}k)}{\Delta(i\omega k)}\right] \frac{dk}{k} .\end{split}$$

Using that *<sup>ω</sup>*<sup>3</sup> <sup>=</sup> <sup>−</sup>1 and *<sup>U</sup>*(−*k*) = <sup>−</sup>*U*(*k*) the above expression is simplified to

$$\widetilde{\mathbf{C}}\_{2} = \frac{1}{4\pi} \int\_{l\_{2}} \mathcal{P}E(-ik)E(-i\omega k) \frac{\Delta(ik)\mathcal{U}(k) + \Delta(i\omega k)\mathcal{U}(\omega k) + \Delta(i\omega^{2}k)\mathcal{U}(\omega^{2}k)}{\Delta(ik)\Delta(i\omega k)} \frac{dk}{k}.\tag{27}$$

Employing the global relation (21) we obtain

$$\check{C}\_2 = \frac{1}{4\pi i} \int\_{l\_2} \mathcal{P}E(-ik)E(-i\omega k) \frac{G(k)}{\Delta(ik)\Delta(i\omega k)} \frac{dk}{k} . \tag{28}$$

In summary, the solution takes the form

$$q = \sum\_{j=1}^{6} F\_j + \sum\_{j=1}^{6} \widetilde{\mathcal{C}}\_{j\star} \tag{29}$$

where *Fj* is defined by

$$F\_{\bar{\jmath}} = \frac{1}{4\pi i} \int\_{l\_{\bar{\jmath}}} \mathcal{P}E\left(-i\omega^{\bar{\jmath}-1}k\right) \left[D\left(\omega^{\bar{\jmath}-1}k\right) + \frac{G\left(\omega^{\bar{\jmath}-1}k\right)}{\Delta\left(i\omega^{\bar{\jmath}-1}k\right)}\right] \frac{dk}{k} \tag{30}$$

and *<sup>C</sup>*,*<sup>j</sup>* is defined by

$$\widetilde{C}\_{j} = \frac{1}{4\pi i} \int\_{l\_{j}} \mathcal{P}E(-i\omega^{j-2}k)E(-i\omega^{j-1}k) \frac{G(\omega^{j-2}k)}{\Delta(i\omega^{j-1}k)\Delta(i\omega^{j-2}k)} \frac{dk}{k}.\tag{31}$$

Note also that the integrals of *<sup>C</sup>*,*<sup>j</sup>* can be deformed on a sector of angle <sup>2</sup>*<sup>π</sup>* <sup>3</sup> . For example, in *<sup>C</sup>*,<sup>2</sup> the ray *<sup>l</sup>*<sup>2</sup> can be deformed in a ray *l* <sup>2</sup> in the sector arg *<sup>k</sup>* <sup>∈</sup> (*π*, <sup>5</sup>*<sup>π</sup>* <sup>3</sup> ); analogous results are valid for the remaining {*C*,*j*}<sup>6</sup> 1.

Observing that *G*(*ωk*) = *G*(*k*), Equation (29) can be further simplified to

$$q = \frac{1}{4\pi i} \sum\_{j=1}^{6} \int\_{l\_j} \mathcal{P}\left[E(-i\omega^{j-1}k)D(\omega^{j-1}k) + \frac{E(-i\omega^j k)G(\omega^{j-1}k)}{\Delta(i\omega^{j-1}k)\Delta(i\omega^{j-2}k)}\right] \frac{dk}{k}.\tag{32}$$

In order to write the integral representation in a more compact form we make the change of variables *<sup>k</sup>* <sup>→</sup> *<sup>ω</sup>*1−*<sup>j</sup> <sup>k</sup>* in the integrals in *Fj* and *<sup>C</sup>*,*j*. In this procedure:


$$\text{i.3. } \quad \text{the exponent } \mathcal{P} = \epsilon^{i\beta \left( kz - \frac{z}{k} \right)} \text{ becomes } \epsilon^{i\beta \left( \omega^{1-j} kz - \frac{z}{\omega^{1-j} k} \right)} \text{.}$$

4. the remaining integrands are equal to the corresponding integrands in *<sup>F</sup>*<sup>1</sup> and *<sup>C</sup>*,1.

Thus, we obtain

$$q = \frac{1}{4\pi i} \int\_{l\_1} \mathcal{T} \left[ E(-ik)D(k) - \frac{E(-i\omega k)}{\Delta(ik)\Delta(i\omega^2 k)} G(k) \right] \frac{dk}{k}.\tag{33}$$

where

$$\mathcal{T} = \sum\_{j=1}^{6} e^{i\beta \left(\omega^{1-j}k\overline{z} - \frac{\pi}{\omega^{1-j}k}\right)}.$$

We make the change of variables *k* → −*ik* in the integrand of (33), so that the contour of integration transforms from the negative imaginary axis *l*<sup>1</sup> to the real imaginary axis, and we summarize the above result in the form of a proposition.

**Proposition 1.** *Let q satisfy the modified Helmholtz Equation* (2) *in the interior of a regular hexagon defined in* (13)*. Assume that on each side of this hexagon an odd symmetric Dirichlet boundary condition is prescribed, namely,*

$$q^{(j)}(s) = d(s), \qquad s \in \left[ -\frac{l}{2}, \frac{l}{2} \right], \ j = 1, \ldots, 6,$$

*with d*(−*s*) = <sup>−</sup>*d*(*s*) *and d* − *l* 2 = *d l* 2 = 0*.* *The solution q can be computed in closed form:*

$$q(z,\overline{z}) = \frac{1}{4\pi i} \int\_0^\infty R(k, z, \overline{z}) \left[ E(-k)D(-ik) - \frac{E(-\omega k)}{\Delta(k)\Delta(\omega^2 k)} G(-ik) \right] \frac{dk}{k},\tag{34}$$

*where R*(*k*, *z*, *z*¯), *D*(*k*), *E*(*k*), *G*(*k*), Δ(*k*) *are defined as follows:*

$$\begin{aligned} R(k, z, z) &= \sum\_{j=1}^{6} e^{\beta \left(\omega^{1-j} kz + \frac{z}{\omega^{1-j} k}\right)} \\ E(k) &= e^{\beta(k + \frac{1}{k}) \frac{l \sqrt{3}}{2}}, \qquad D(k) = \beta \left(\frac{1}{k} - k\right) \int\_{-\frac{1}{2}}^{\frac{l}{2}} e^{\beta(k + \frac{1}{k})s} d(s) ds, \\ G(k) &= \sum\_{j=1}^{6} E\left(-i\omega^{j-1} k\right) D\left(\omega^{j-1} k\right), \qquad \Delta(k) = E(k) - E(-k), \qquad k \in \mathbb{C}. \end{aligned}$$

#### **4. The Symmetric Even Case**

Applying the condition *U*(−*k*) = *U*(*k*) in (17) we obtain the following equation

$$
\Delta^{+}(ik)\mathcal{U}(k) + \Delta^{+}(i\omega k)\mathcal{U}(\omega k) + \Delta^{+}(i\omega^{2}k)\mathcal{U}(\omega^{2}k) = iG(k), \qquad k \in \mathbb{C}, \tag{35}
$$

where

$$
\Delta^+(k) = E(k) + E(-k)
$$

and *G*(*k*) is known and given in (18).

Following the same stems used in Section 3 we derive the analogue of (28), which yields the following formula for *<sup>C</sup>*,2:

$$\begin{split} \breve{C}\_{2} &= \frac{1}{4i\pi} \int\_{l\_{2}} \mathcal{P}E(-ik)E(-i\omega k) \frac{G(k)}{\Delta^{+}(ik)\Delta^{+}(i\omega k)} \frac{dk}{k} \\ &+ \frac{1}{2\pi} \int\_{l\_{2}} \mathcal{P} \frac{\mathcal{U}(\omega^{2}k)}{\Delta^{+}(ik)\Delta^{+}(i\omega k)} \frac{dk}{k} .\end{split} \tag{36}$$

where in addition to the known part which involves *G*(*k*), there now exists an unknown part which involves *U*(*ω*2*k*).

Thus, the analogue of (29) now takes the form

$$q = \sum\_{j=1}^{6} F\_j + \sum\_{j=1}^{6} A\_j + \sum\_{j=1}^{6} B\_{j\prime} \tag{37}$$

where *Fj* is known function defined by

$$F\_{\bar{\jmath}} = \frac{1}{4\pi i} \int\_{l\_{\bar{\jmath}}} \mathcal{P}E(-i\omega^{\bar{\jmath}-1}k) \left[ D(\omega^{\bar{\jmath}-1}k) + \frac{G(\omega^{\bar{\jmath}-1}k)}{\Delta^{+}(i\omega^{\bar{\jmath}-1}k)} \right] \frac{dk}{k} \,\tag{38}$$

*Aj* is also known and defined by

$$A\_{\bar{j}} = \frac{1}{4\pi i} \int\_{l\_{\bar{j}}} \mathcal{P}E(-i\omega^{j-2}k) \mathcal{E}(-i\omega^{j-1}k) \frac{G(\omega^{j-2}k)}{\Delta^{+}(i\omega^{j-1}k)\Delta^{+}(i\omega^{j-2}k)} \frac{dk}{k},\tag{39}$$

whereas *Bj* is the unknown function defined by

$$B\_{\vec{j}} = \frac{1}{2\pi} \int\_{\vec{l}\_{\vec{j}}} \mathcal{P} \frac{\mathcal{U}(\omega^{\vec{j}}k)}{\Delta^{+}(i\omega^{\vec{j}-1}k)\Delta^{+}(i\omega^{\vec{j}-2}k)} \frac{dk}{k} . \tag{40}$$

It can be shown that each of *Bj* decays exponentially fast as *β* → ∞. The rigorous proof of this statement will be presented elsewhere. In the next section, this fact will be indicated via the numerical evaluation of each of the terms appearing in Equation (37).

#### **5. Illustration of the Results**

#### *5.1. Odd Case*

Below we depict the solution obtained by (34) for various choices of the Dirichlet datum *d*(*s*) and of the parameter *β*. At all the examples we have fixed the length of the side of the hexagon *l* = 2.

For the first example we employ the Dirichlet datum *d*(*s*) = sin(*πs*) and the parameter *β* = 1; see Figure 1.

**Figure 1.** The solution *q* given by (34) for *d*(*s*) = sin(*πs*), *l* = 2 and *β* = 1.

We also depict the deviation of *d*(*s*) from the function obtained by the integral representation (34) evaluated at the side of the hexagon, namely at *x* = *<sup>l</sup>* √3 <sup>2</sup> <sup>=</sup> <sup>√</sup><sup>3</sup> and *<sup>y</sup>* <sup>=</sup> *<sup>s</sup>* <sup>∈</sup> - − *l* 2 , *l* 2 . ≡ [−1, 1]; see Figure 2.

**Figure 2.** The deviation of *q* (given by (34)) from the actual Dirichlet datum *d*(*s*) evaluated at the side of the hexagon; here we employ *d*(*s*) = sin(*πs*), *l* = 2 and *β* = 1.

For the second example we employ the Dirichlet datum *d*(*s*) = sin(*πs*) and the parameter *β* = 1/5; see Figure 3.

**Figure 3.** The solution *q* given by (34) for *d*(*s*) = sin(*πs*), *l* = 2 and *β* = 1/5.

For the third example we employ the Dirichlet datum *d*(*s*) = sin(2*πs*) and the parameter *β* = 1; see Figure 4.

**Figure 4.** The solution *q* given by (34) for *d*(*s*) = sin(2*πs*), *l* = 2 and *β* = 1.

#### *5.2. Even Case*

In this case we employ the Dirichlet datum *d*(*s*) = cos *<sup>π</sup>* 2 *s* and the parameter *β* = 1 at the known part of the rhs of the formula (37), namely the expression

$$\sum\_{j=1}^{6} F\_j + \sum\_{j=1}^{6} A\_{j\prime} \tag{41}$$

where *Fj* and *Aj* are given by (38) and (39), respectively; see Figure 5.

**Figure 5.** The known part of the solution *q* given by (41) for *d*(*s*) = cos *<sup>π</sup>* 2 *s* , *l* = 2 and *β* = 1.

We also depict the deviation of *d*(*s*) from the above expression evaluated at the side of the hexagon, namely at *<sup>x</sup>* <sup>=</sup> <sup>√</sup><sup>3</sup> and *<sup>y</sup>* <sup>=</sup> *<sup>s</sup>* <sup>∈</sup> [−1, 1]. This is equal to the contribution <sup>∑</sup><sup>6</sup> *<sup>j</sup>*=<sup>1</sup> *Bj*, with *Bj* given by (40); see Figure 6.

**Figure 6.** The deviation of the known part of the solution *q* given by (41) from the actual Dirichlet datum *d*(*s*) = cos *<sup>π</sup>* 2 *s* , *l* = 2 and *β* = 1, evaluated at the side of the hexagon.

Furthermore, we depict the latter contribution for the different values of *β* = <sup>1</sup> 4 , 1 <sup>2</sup> , 1, 2, 4, where it is clearly shown that the error decreases drastically with the increase of *β*; see Figure 7. We observe exponential decay for *z* = *zj*, *j* = 1, ... , 6: in Figure 8 we depict the deviation from the actual Dirichlet data for three points on side (1) of the hexagon, namely *α*<sup>1</sup> = √3, 0 , *α*<sup>2</sup> = √3, <sup>3</sup> 10 , *α*<sup>3</sup> = √3, <sup>9</sup> 10 , with *β* in the intervals *I*<sup>1</sup> = [1, 8], *I*<sup>2</sup> = [1, 10], *I*<sup>3</sup> = [1, 58], respectively.

**Figure 7.** The deviation of the known part of the solution *q* given by (41) from the actual Dirichlet datum *d*(*s*) = cos *<sup>π</sup>* 2 *s* and *l* = 2, evaluated at the side of the hexagon. This deviation is depicted for the different values of *β* = <sup>1</sup> 4 , 1 <sup>2</sup> , 1, 2, 4, and it decreases drastically with the increase of *β*.

**Figure 8.** The deviation of the known part of the solution *q* given by (41) from the actual Dirichlet datum *d*(*s*), evaluated at three points of side (1) of the hexagon, namely *α*<sup>1</sup> = √3, 0 in red, *α*<sup>2</sup> = √3, <sup>3</sup> 10 in blue, *α*<sup>3</sup> = √3, <sup>9</sup> 10 in black. The deviation is depicted against *β* and it displays exponential decay.

#### **6. Conclusions**

In this work we have presented the explicit solution of a particular boundary value problem for the modified Helmholtz equation in a regular hexagon: we have solved the case where the same Dirichlet datum *d*(*s*) is prescribed in all sides of the hexagon, and this function is odd. This explicit solution is described in Proposition 1. We have also obtained an approximate analytical representation for the solution for the case that *d*(*s*) is even. The exact representation is given by Equation (37), where the terms *Fj* and *Aj* are given in terms of *d*(*s*), but the terms *Bj* involve the unknown Neumann boundary value. However, these terms are exponentially small as *β* → ∞. Thus, for the case of large *β*, Equation (37) provides the solution to this problem with an exponentially small error. The above analytical results were verified numerically in Section 5. The rigorous investigation on the analytical and numerical accuracy of the latter approximate formula will be presented in future work.

It should be noted that the arbitrary Dirichlet problem can be decomposed into 6 separate and simpler Dirichlet BVPs, which are defined in Section 2.3; the first of these BVPs is the symmetric Dirichlet problem. The analysis of the remaining problems is a work in progress.

**Author Contributions:** Conceptualization, K.K. and A.S.F.; methodology, K.K. and A.S.F.; validation, K.K. and A.S.F.; formal analysis, K.K. and A.S.F.; investigation, K.K. and A.S.F.; writing–original draft preparation, K.K. and A.S.F.; writing–review and editing, K.K. and A.S.F.; visualization, K.K. and A.S.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** A.S.F. was supported by EPSRC, UK, via a senior fellowship.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
