The Modiﬁed Helmholtz Equation on a Regular Hexagon—The Symmetric Dirichlet Problem

: Using the uniﬁed transform, also known as the Fokas method, we analyse the modiﬁed 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 veriﬁcation is also provided.


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][5][6][7][8][9][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][15][16][17][18], as well as for studying problems from control theory [19].
The Fokas method is based on two basic ingredients: (1) a global relation, which is an algebraic equation that involves certain transforms of all (known and unknown) boundary values. (2) an integral representation of the solution, which involves transforms of all boundary values.
For linear PDEs, the Fokas method involves the following: • Given a PDE, define its formal adjoint and construct a one parameter family of solutions of this equation.

•
By employing the given PDE and its adjoint, obtain a one parameter family of equations in conservation form. This family, together with Green's theorem, yield the global relation.

•
The above family also gives rise to a certain closed differential form. The spectral analysis of this form gives rise to a scalar Riemann-Hilbert problem, which consequently yields an integral representation of the solution. This representation involves integral transforms of all the boundary values, and since some of them are not prescribed as boundary conditions, this form of solution is not yet effective.

•
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 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][51][52][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.

The Basic Elements
The equation investigated here is the modified Helmholtz equation in the interior of the regular hexagon, D, namely, where q(x, y) is a real valued function and β > 0. Using complex coordinates,

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 equatioñ Multiplying Equation (2) byq, Equation (3) by q and subtracting, we find Using in (5) the special solutionq = e −iβ(kz−z k ) and employing Green's theorem, we obtain where W is defined by Suppose that Ω is the polygon defined via the points z 1 , z 2 , . . . , z n , z n+1 = z 1 . Then (6) gives the following global relation for the modified Helmholtz in this polygon: where q j (k) n 1 are defined bŷ or equivalently (in local coordinates) bŷ In Equation (10) we have used the identity where s is the arclength on the boundary z(s) = x(s) + iy(s) of the polygon and q N 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 This procedure yields the following theorem, proven in [6]: Let Ω be the interior of a convex closed polygon in the complex z-plane, with corners z 1 , . . . , z n , z n+1 ≡ 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 where {q j (k)} n 1 are defined by (10), and {l j } n 1 are the rays in the complex k-plane oriented from zero to infinity.
Observe that the solution given in (12) is given in terms of {q j } n 1 which involve integral transforms of both q and q N on the boundary, i.e., both known and unknown functions.

The Dirichlet Problem on a Regular Hexagon
Let D ⊂ C be the interior of a regular hexagon with vertices {z j } 6 1 , where l is the length of the side and ω = e iπ 3 . The sides {(z j , z j+1 )} 6 1 , z 7 ≡ z 1 will be referred to as sides {(j)} 6 1 . For the sides {(j)} 6 1 the following parametrizations will be used: The general Dirichlet problem can be uniquely decomposed to 6 simpler Dirichlet problems, by employing the decomposition indeed the determinant of the matrix ω (j−1)(i−1) i,j=1,...,6 is non-zero (Its value is 216 = 6 3 , and for the general case Det 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.

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 The above 'symmetry' property also holds for the Neumann boundary values. This fact is the consequence of the following three observations: • The modified Helmholtz operator namely under rotation of 2π/3. Since the Dirichlet data are invariant under this rotation, then the (unique) solution q(z,z) of the Helmholtz equation is also invariant under this rotation.

•
If q is invariant under this transformation, then the differential form q z dz is also invariant under the transformation z → ωz: • Evaluating the above differential form on each side we obtain where the second equality is a direct consequence of the fact that the Dirichlet data are invariant under this rotation. Thus, Applying the parametrization of the regular hexagon on Equation (10) we obtain: where E(k), D(k) and U(k) are defined by The function D(k) is known, whereas the unknown function U(k) contains the unknown Neumann boundary value u(s) = q N .
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:

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 6 1 to the solution representation (19) can be computed explicitly.
Applying the condition U(−k) = −U(k) in (17) we obtain the equation where G(k) is given in (18) and ∆(k) is defined by Solving (21) for U(k) and substituting the resulting expression in (15) we find The functionsq 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 1 (the negative imaginary axis).
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 The unknown part involves the functions U(ωk) and U(ω 2 k) and is given by In what follows we will show that the contribution of the unknown functions, namely of the sum ∑ 6 1 C j , can be computed in terms of the given boundary conditions. The first integral in the rhs of C 1 can be deformed from l 1 to l 1 , where l 1 is a ray with 7π 6 ≤ arg k ≤ 3π 2 ; choosing l 1 ≡ l 2 we obtain The above deformation is justified, since it can be shown that the integrand ofĈ 1 is bounded and analytic in the region where arg k ∈ 7π 6 , 3π 2 : letting a = e i π 6 , we can rewrite the first term of the integrand ofĈ 1 in the form We observe the following: • The zeros of ∆(ik) occur when ik + 1 ik ∈ e −i π 2 R, thus k ∈ R.

•
The function E is bounded and analytic for arg k ∈ [ 7π 6 , 3π 2 ]. Indeed, since k is at the lower half plane, then which is bounded if Re ω 2 k ≥ 0.
Similar considerations apply to the second term of the integrand ofĈ 1 ; this term can be rewritten in the form We observe the following: is bounded and analytic for arg k ∈ [ 5π 6 , 11π 6 ], namely in the region where Re(ω 2 k) ≥ 0.

•
In the lower half plane Thus, it is bounded and analytic for arg k ∈ [ 7π 6 , 3π 2 ]. Using the underlined symmetries, we can express the integral representation of the solution in the form where F j and C j are given by applying in (23) and (24) the following rotations: k → ω j−1 k, l 1 → l j , l 2 → l j+1 , j = 2, . . . , 6; l 7 := l 1 .
We define C j =Ĉ j−1 +Č j , j = 1, . . . , 6, where we employ the notationĈ 0 =Ĉ 6 . Then, we rewrite the expression in (25) in the form Thus, it is sufficient to compute the contribution { C j } 6 1 . In this direction we find (via rotation) thať Thus Using that ω 3 = −1 and U(−k) = −U(k) the above expression is simplified to Employing the global relation (21) we obtain In summary, the solution takes the form where F j is defined by and C j is defined by Note also that the integrals of C j can be deformed on a sector of angle 2π 3 . For example, in C 2 the ray l 2 can be deformed in a ray l 2 in the sector arg k ∈ (π, 5π 3 ); analogous results are valid for the remaining { C j } 6 1 .
Observing that G(ωk) = G(k), Equation (29) can be further simplified to In order to write the integral representation in a more compact form we make the change of variables k → ω 1−j k in the integrals in F j and C j . In this procedure: 1.
the fraction dk k remains invariant; 2.
the rays l j become l 1 ; 3. the exponent P = e iβ(kz−z k ) becomes e iβ ω 1−j kz−z the remaining integrands are equal to the corresponding integrands in F 1 and C 1 .
Thus, we obtain where 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 1 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, with d(−s) = −d(s) and d − l The solution q can be computed in closed form: where R(k, z,z), D(k), E(k), G(k), ∆(k) are defined as follows:

The Symmetric Even Case
Applying the condition U(−k) = U(k) in (17) we obtain the following equation where 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 C 2 : 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 where F j is known function defined by A j is also known and defined by whereas B j is the unknown function defined by It can be shown that each of B j 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).

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.  For the second example we employ the Dirichlet datum d(s) = sin(πs) and the parameter β = 1/5; see Figure 3. For the third example we employ the Dirichlet datum d(s) = sin(2πs) and the parameter β = 1; see Figure 4.

Even Case
In this case we employ the Dirichlet datum d(s) = cos π 2 s and the parameter β = 1 at the known part of the rhs of the formula (37), namely the expression where F j and A j are given by (38) and (39), respectively; see Figure 5. We also depict the deviation of d(s) from the above expression evaluated at the side of the hexagon, namely at x = √ 3 and y = s ∈ [−1, 1]. This is equal to the contribution ∑ 6 j=1 B j , with B j given by (40); see Figure 6. Furthermore, we depict the latter contribution for the different values of β = 1 4 , 1 2 , 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 = z j , 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 α 1 = √ 3, 0 , α 2 = √ 3, 3 10 , α 3 = √ 3, 9 10 , with β in the intervals I 1 = [1,8], I 2 = [1, 10], I 3 = [1, 58], respectively.  in blue, α 3 = √ 3, 9 10 in black. The deviation is depicted against β and it displays exponential decay.

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 F j and A j are given in terms of d(s), but the terms B j 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.