# Propagation of Waves from Finite Sources Arranged in Line Segments within an
Infinite Triangular Lattice

David Kapanadze
A. Razmadze Mathematical Institute, TSU,
Merab Aleksidze II Lane 2, Tbilisi 0193, Georgia
Free University of Tbilisi,
240 David Aghmashenebeli Alley, Tbilisi 0159, Georgia
<EMAIL_ADDRESS>
& Zurab Vashakidze
The University of Georgia (UG),
77a M. Kostava St., Tbilisi 0171, Georgia
Ilia Vekua Institute of Applied Mathematics (VIAM), TSU,
2 University St., Tbilisi 0186, Georgia
<EMAIL_ADDRESS><EMAIL_ADDRESS>0000-0002-5006-8615
0000-0001-8736-6213

###### Abstract

This paper examines the propagation of time-harmonic waves in a two-
dimensional triangular lattice with a lattice constant $a=1$. The sources are
positioned along line segments within the lattice. Specifically, we
investigate the discrete Helmholtz equation with a wavenumber
$k\in\left(0,2\sqrt{2}\right)$, where input data is prescribed on finite rows
or columns of lattice sites. We focus on two main questions: the efficacy of
the numerical methods employed in evaluating the Green’s function, and the
necessity of the cone condition. Consistent with a continuum theory, we employ
the notion of radiating solution and establish a unique solvability result and
Green’s representation formula using difference potentials. Finally, we
propose a numerical computation method and demonstrate its efficiency through
examples related to the propagation problems in the left-handed two-
dimensional inductor-capacitor metamaterial.

_Keywords:_ discrete Helmholtz equation $\cdot$ exterior Dirichlet problem
$\cdot$ metamaterials $\cdot$ triangular lattice model.

_MSC 2010:_ 78A45 $\cdot$ 35J05 $\cdot$ 39A14 $\cdot$ 39A60 $\cdot$ 74S20.

##  1 Introduction

We are interested in methods for solving the Dirichlet-type boundary value
problem for the discrete Helmholtz equation with uniform fixed mesh spacing.
This equation commonly arises in lattice dynamics, where the difference
equations can model wave propagation. Similar static issues can occur in
various scenarios. For instance, a study of lattice defects can be found in
Ref. [30]. Lattice dynamics traditionally finds applications in the field of
solid-state physics [22]. One can study the propagation of waves through a
periodic arrangement of interacting cells, where each cell has the same
arrangement of interacting atoms. The propagation of waves through such a
perfect lattice is a well-known topic [4]. If there is a defect in the
lattice, it will cause the waves to scatter. Wave propagation through discrete
structures remains an active area of research today.

The triangular lattice is one of the five two-dimensional Bravais lattice
types and appears naturally in applications, especially in crystals and
materials with hexagonal symmetry [4, 6]. For example, graphene, as observed
in nature, is a two-dimensional material comprising a monolayer of carbon
atoms arranged in a honeycomb structure, which can be described as a hexagonal
lattice with a two-atom basis (see [8, 25]). According to the study presented
in [18], for certain configurations, problems defined on the hexagonal lattice
can be equivalently reduced to analogous problems formulated on the triangular
lattice. Additionally, one can consider two-dimensional passive propagation
media, such as a host microstrip line network periodically loaded with series
capacitors and shunt inductors for signal processing and filtering (as
depicted in Figure 1). This type of inductor-capacitor lattice is called a
negative-refractive-index transmission-line (NRI-TL) metamaterial [7] or
simply left-handed two-dimensional metamaterial. Suppose that monochromatic
inputs are applied to finite rows/columns of lattice sites. Assume that the
number of unit cells in this slab is large enough to make it prohibitively
expensive to solve numerically for the voltage/current at every cell in the
lattice until the system reaches a steady state. As a simplifying strategy, it
can be expected that the limiting case, when the lattice is effectively
infinite, is more flexible to analysis and provides a good approximation of
the steady-state output at an exterior boundary.

The discrete Helmholtz equation is a fundamental equation in physics,
engineering, and mathematics, widely employed in various scientific
disciplines to describe wave propagation in discrete media. Its applicability
extends to fields such as electromagnetics, acoustics, optics, and quantum
mechanics. Additionally, discrete Helmholtz equations are closely related to
discrete Schrödinger equations, which naturally arise in the tight-binding
model of electrons in crystals [3, 29, 11]. Similar equations also arise in
studies involving time-harmonic elastic waves in lattice models of crystals
[5, 22, 21]. Moreover, the Helmholtz equation is closely related to the
Maxwell system (for time-harmonic fields); solutions of the scalar Helmholtz
equation are used to generate solutions of the Maxwell system (Hertz
potentials), and every component of the electric and magnetic field satisfies
an equation of Helmholtz type (see Kirsch and Hettlich [20]). Continuous
models for the Helmholtz equation with smooth boundaries have been thoroughly
studied by Colton and Kress [9] and in related works. However, the formulation
of a discrete analogue of the Rayleigh-Sommerfeld scattering theory across
different lattice structures remains an active area of research. Results
pertaining to square and triangular lattices have been previously obtained in
[14, 17].

This paper extends our investigation of the discrete Helmholtz equation and
its associated exterior problems in a two-dimensional triangular lattice with
a lattice constant $a=1$. We concentrate on two central questions: the
efficacy of the numerical technique employed to compute the lattice Green’s
function and the requirement of the cone condition. To this end, the effect of
finite sources arranged on line segments within an infinite triangular lattice
has been studied. Mathematical modelling of the propagation problem under
consideration leads us to study an exterior problem for a two-dimensional
discrete Helmholtz equation with Dirichlet boundary conditions. We use the
results from [16] and carry out our investigation without passing to the
complex wavenumber. Specifically, we employ the radiation conditions and
asymptotic estimates from [15], alongside the Rellich-Vekua type theorem
stated by Isozaki and Morioka [13]. It is worth noting that these problems can
be explored through the limiting absorption principle, cf. e.g. Sharma in Ref.
[28]. This work delves into the diffraction phenomena occurring in triangular
and hexagonal lattices caused by finite cracks and rigid constraints. We
establish the existence and uniqueness of solutions for the problem under
consideration within the interval $k\in\left(0,2\sqrt{2}\right)$. Moreover, we
derive Green’s representation formula for the solution using difference
potentials. The proposed approach aims to simplify the solution of the
prescribed boundary value problem by solving a linear system of boundary
equations. Furthermore, the solution to the original problem is represented by
a linear combination of lattice Green’s functions and the solution derived
from a reduced linear system of boundary equations. The precision of our
solution is mostly linked to the accurate computation of lattice Green’s
functions. To achieve this objective, instead of using the numerical
quadrature technique for computing lattice Green’s functions (due to the rapid
oscillatory behaviour of these functions), we adopt the approach proposed by
Berciu and Cook [2]. This method enables the computation of Green’s functions
through elementary operations alone, thus avoiding the need for integration,
and yielding physically consistent solutions. For numerical illustration, two
sample problems are considered for $k=2$. The first example features four
discrete boundary points aligned along a single line, while the second
comprises ten discrete boundary points situated on parallel line segments.
Numerical results are detailed in Section 5.

##  2 Formulation of the problem

In accordance with customary mathematical notation, we denote the sets of
integers, positive integers, real numbers, and complex numbers by the symbols
$\mathbb{Z}$, $\mathbb{N}$, $\mathbb{R}$, and $\mathbb{C}$, respectively.

Let $\left\\{\mathcal{V},\mathcal{E}\right\\}$ be a periodic simple graph that
defines a two-dimensional infinite triangular lattice $\mathfrak{T}$, where

$\mathcal{V}=\left\\{T\left(x_{1},x_{2}\right)\in\mathbb{R}^{2}:\left(x_{1},x_{2}\right)\in{\mathbb{Z}}^{2}=\mathbb{Z}\times\mathbb{Z}\right\\}\,,$
(2.1)

denotes the vertex set, while the edge set is denoted by $\mathcal{E}$. The
endpoints $\left(v,w\right)\in\mathcal{V}\times\mathcal{V}$ are considered
adjacent if and only if $\left|v-w\right|=1$. Additionally, let $T$ denote a
two-dimensional coordinate transformation, which is defined as follows:

$T\left(x_{1},x_{2}\right)=\left(x_{1}+\frac{x_{2}}{2},\frac{x_{2}\sqrt{3}}{2}\right)\,.$

Figure 1: A fragment of a two-dimensional infinite triangular lattice,
demonstrating its conceptual sense.

In Figure 1, a fragment of a two-dimensional infinite triangular lattice
$\mathfrak{T}$ is illustrated. The connection between
$x=\left({x}_{1},{x}_{2}\right)\in\mathbb{Z}^{2}$ and the Euclidean
coordinates of the vertexes (black dots) is established through the mapping
$\left({x}_{1},{x}_{2}\right)\mapsto\left({x}_{1}+{x}_{2}/2,\sqrt{3}{x}_{2}/2\right)$.
For illustrative purposes, Figure 1 additionally depicts the lattice vertex
$T\left(2,2\right)$, with Cartesian coordinates $\left(3,\sqrt{3}\right)$. It
is derived from the point $\left(2,2\right)\in\mathbb{Z}^{2}$ through the
coordinate transformation $T$. The nearest neighbour interactions on the
triangular lattice are depicted with thick blue lines. Conceptually, the
triangular lattice can be understood as a two-dimensional left-handed
metamaterial consisting of inductor-capacitor elements. Specifically, it can
be likened to a periodically loaded host transmission line with series
capacitors and shunt inductors.

We define the $6$-neighbourhood ${F}_{w}^{0}$ as the set of points
$v\in\mathcal{V}$ such that the distance between $w$ and $v$ is equal to $1$,
where $w\in\mathcal{V}$. Additionally, we define the neighbourhood ${F}_{w}$
as the union of ${F}_{w}^{0}$ and the singleton set containing $w$. Further, a
region $R\subset\mathcal{V}$ can be identified as a set consisting of points
in $\mathcal{V}$, where there exist disjoint nonempty subsets $\mathring{R}$
and $\partial R$ of $R$ satisfying the following condition:

  1. (a)

$R=\mathring{R}\cup\partial R$,

  2. (b)

if $w\in\mathring{R}$ then ${F}_{w}\subset R$,

  3. (c)

if $w\in\partial R$ then there is at least one point $v\in{F}_{w}^{0}$ such
that $v\in\mathring{R}$.

Since subsets $\mathring{R}$ and $\partial R$ are not uniquely defined by $R$,
we assume that for a given region $R$ in $\mathcal{V}$, $\mathring{R}$ and
$\partial R$ are given and fixed. We define a point $w$ as an interior or
boundary point of $R$ if $w\in\mathring{R}$ or $w\in\partial R$, respectively.
Additionally, we define a region $R\subset\mathcal{V}$ as connected if for any
$v,\tilde{v}\in R$, there exists a sequence
${w}^{\left(1\right)},\dots,{w}^{\left(n\right)}\in R$ such that
${w}^{\left(1\right)}=v$, ${w}^{\left(n\right)}=\tilde{v}$, and for all $0\leq
i\leq n-1$, $\left|{w}^{\left(i\right)}-{w}^{\left(i+1\right)}\right|=1$. By
definition, a region $R$ with a single interior point $w$ is connected and
coincides with ${F}_{w}$.

Let $\Gamma_{j}$, where $j=1,2,\dots,n$ and $n\in\mathbb{N}$, denote a finite
row of lattice sites. Consider a region $\mathring{\mathcal{W}}$ defined as
the set difference
$\mathring{\mathcal{W}}=\mathcal{V}\setminus\partial\mathcal{W}$, where
$\partial\mathcal{W}$ is the union of the sets $\Gamma_{j}$, for
$j=1,2,\dots,n$. We assume that the lattice sites $\Gamma_{j}$ are arranged in
a manner such that $\mathring{\mathcal{W}}$ is a connected region and
satisfies the cone condition, as defined in reference [16]. Lastly, we
emphasize that $\mathcal{W}$ is the union of $\mathring{\mathcal{W}}$ and
$\partial\mathcal{W}$, and that $\mathcal{V}$ and $\mathcal{W}$ are equivalent
as sets.

Given the problem and assumptions stated in the introduction, we consider a
system in which each node $w\in\mathring{\mathcal{W}}$ is connected to a
common ground plane via an inductor, and each node
$w\in\mathring{\mathcal{W}}$ is connected to its six nearest neighbours via a
capacitor. We assume that all inductances are equal to a positive constant
$L$, and all capacitances are equal to a positive constant $C$. Under these
conditions, Kirchhoff’s laws of voltage and current imply a second-order
differential equation for the voltage $U\left(w,t\right)$ across the inductor
at node $w$, as follows:

${L}{C}\frac{{\mathrm{d}}^{2}}{\mathrm{d}{t}^{2}}\left(\Delta_{d}U\left(w,t\right)\right)=U\left(w,t\right)\,.$
(2.2)

Here, the discrete Laplacian operator $\Delta_{d}$ is represented by a
$7$-point stencil and can be expressed as:

$\Delta_{d}U\left(w,t\right)=\sum_{v\in{F}_{w}^{0}}U\left(v,t\right)-6U\left(w,t\right)\,.$

It should be noted that equation (2.2) is valid for all
$w\in\mathring{\mathcal{W}}$, and the boundary of the region $\mathcal{W}$,
denoted by $\partial\mathcal{W}$, is subject to a time-dependent boundary
condition specified by a function $f\left(v\right)$, where $v$ denotes a node
on the boundary. The boundary condition is given by:

$U\left(v,t\right)=f\left(v\right){e}^{-\imath\omega t}\,.$ (2.3)

It is important to note that $\imath$ is used to denote the imaginary unit,
and $f:\partial\mathcal{W}\to\mathbb{C}$ is a given function. Under the
assumption that at time $t=0$, $U\left(w,t\right)$ and all its derivatives are
zero for all $w\in\mathring{\mathcal{W}}$, a wave propagates into the lattice
as $t$ increases due to the time-dependent boundary condition, eventually
leading the system to a steady state. In this steady state, the solution takes
the form $U\left(w,t\right)=u\left(w\right){e}^{-\imath\omega t}$, where
$u\left(w\right)$ is a complex-valued function that is constant in time. By
substituting this expression into equations (2.2) and (2.3), which correspond
to the discrete Helmholtz equation in $\mathcal{W}$, we arrive at the
following problem:

$\displaystyle\left(\Delta_{d}+k^{2}\right)u\left(w\right)$
$\displaystyle=0\,,$ in $\displaystyle\leavevmode\nobreak\
\mathring{\mathcal{W}}\,,$ (2.4a) $\displaystyle u\left(v\right)$
$\displaystyle=f\left(v\right)\,,$ on $\displaystyle\leavevmode\nobreak\
\partial\mathcal{W}\,.$ (2.4b)

The problem (2.4a)-(2.4b) involves the Laplacian operator $\Delta_{d}$ and a
parameter $k$ that is related to the frequency $\omega$ of the wave through
the equation ${k}^{2}={\left({\omega}^{2}{L}{C}\right)}^{-1}$, where $L$ and
$C$ are the inductance and capacitance of the lattice, respectively.

To enhance efficiency, we have opted to work within a simplified coordinate
system that is denoted by ${\mathbb{Z}}^{2}$, as employed in (2.1).
Specifically, we represent each node $w$ within the set $\mathcal{W}$ by a
pair of integers, $\left(x_{1},x_{2}\right)$. Within this system, every node
is associated with six neighbours: $\left(x_{1}-1,x_{2}\right)$,
$\left(x_{1}+1,x_{2}\right)$, $\left(x_{1},x_{2}-1\right)$,
$\left(x_{1},x_{2}+1\right)$, $\left(x_{1}+1,x_{2}-1\right)$, and
$\left(x_{1}-1,x_{2}+1\right)$. Consequently, we assume all definitions that
have been previously introduced with respect to content using integer
coordinates, without explicit mention.

Let us reformulate equation (2.2) within the context of the selected
coordinate system. To be precise, we set
$U\left(x,t\right)=U\left(w,t\right)$, where
$x=\left(x_{1},x_{2}\right)\in{\mathbb{Z}}^{2}$ and
$w=T\left(x_{1},x_{2}\right)\in\mathcal{W}$. With these definitions, we
obtain:

${L}{C}\frac{{\mathrm{d}}^{2}}{\mathrm{d}{t}^{2}}\left(\Delta_{d}U\left(x,t\right)\right)=U\left(x,t\right)\,.$
(2.5)

If we reformulate the boundary condition (2.3) in the coordinate system
${\mathbb{Z}}^{2}$ and follow the same reasoning as above, we obtain the
following boundary value problem of Dirichlet-type for the time-harmonic
discrete waves in
$\Omega=\mathring{\Omega}\cup\partial\Omega\subset{\mathbb{Z}}^{2}$

$\displaystyle\left(\Delta_{d}+k^{2}\right)u\left(x\right)$
$\displaystyle=0\,,$ in $\displaystyle\leavevmode\nobreak\
\mathring{\Omega}\,,$ (2.6a) $\displaystyle u\left(y\right)$
$\displaystyle=f\left(y\right)\,,$ on $\displaystyle\leavevmode\nobreak\
\partial\Omega\,,$ (2.6b)

where the discrete Laplacian operator is expressed as (its action is
represented by a stencil as illustrated in Figure 2)

$\displaystyle{\Delta}_{d}u\left(x\right)=$ $\displaystyle
u\left(x+e_{1}\right)+u\left(x-e_{1}\right)+u\left(x+e_{2}\right)+u\left(x-e_{2}\right)$
$\displaystyle+u\left(x+e_{1}-e_{2}\right)+u\left(x-e_{1}+e_{2}\right)-6u\left(x\right)\,,$

Figure 2: $7$-point stencil of the discrete Laplacian operator ${\Delta}_{d}$
in the coordinate system ${\mathbb{Z}}^{2}$.

$e_{1}=\left(1,0\right)$ and $e_{2}=\left(0,1\right)$ denote the standard
basis vectors of ${\mathbb{Z}}^{2}$, representing the directions of the
positive $x_{1}$ and $x_{2}$ axes, respectively. $\mathring{\Omega}$ denotes
the interior of the domain, and $\partial\Omega$ denotes the boundary of the
domain. The function $f:\partial\Omega\rightarrow{\mathbb{C}}$ is a given
boundary condition. Note that the coordinate system ${\mathbb{Z}}^{2}$ is used
to represent the nodes of the discrete system, and all previously introduced
definitions apply in this context.

It is established that equation (2.2) permits plane wave solutions
$U\left(x\right)={A}{e}^{-\imath\left({\xi}_{1}{x}_{1}+{\xi}_{2}{x}_{2}+{\omega}{t}\right)}$,
where $A\in\mathbb{C}$ denotes a constant, provided that the ensuing
dispersion relation holds:

${\omega}^{2}=\dfrac{1}{{4}{L}{C}\left({\sin}^{2}\frac{\xi_{1}}{2}+{\sin}^{2}\frac{\xi_{2}}{2}+{\sin}^{2}\frac{\xi_{1}-\xi_{2}}{2}\right)}\,.$

Assuming the fulfilment of the dispersion relation for
$\left(\xi_{1},\xi_{2}\right)\in{\left[-\pi,\pi\right]}^{2}$, identified as
the first Brillouin zone (see [5]). Consequently, the following holds:

${k}^{2}={4}\left({\sin}^{2}\frac{\xi_{1}}{2}+{\sin}^{2}\frac{\xi_{2}}{2}+{\sin}^{2}\frac{\xi_{1}-\xi_{2}}{2}\right)\in\left[0,9\right]\,.$

Undoubtedly, other values of $k$ are also subject to investigation, but these
cases are relatively straightforward and are not addressed within the scope of
this discussion.

As a reminder, the spectrum of the negative discrete Laplacian is known to be
absolutely continuous within the interval $\left[0,9\right]$. However, there
exists a set of exceptional points $\left\\{0,8,9\right\\}$ within this
interval where the limiting absorption principle fails, as detailed in [27].
Therefore, it is assumed that
$k\in\left(0,3\right)\setminus\left\\{2\sqrt{2}\right\\}$.

The objective of this study is to investigate the existence and uniqueness of
the function $u:\Omega\to\mathbb{C}$ that satisfies the discrete Helmholtz
equation (2.6a) with $k\in\left(0,2\sqrt{2}\right)$ and the boundary condition
(2.6b). This problem shall be denoted as _Problem_
$\mathcal{P}_{\mathrm{ext}}$ in the following discussion.

##  3 Green’s representation formula

This section primarily recollects the results from [16, 17, 15]. Let $R$ be a
region in the two-dimensional integer lattice ${\mathbb{Z}}^{2}$. For
$j=1,\dots,6$, let ${\left(\partial R\right)}_{j}$ denote the set of all
boundary points $y\in\partial R$ such that $y-{e}_{j}\in\mathring{R}$, where
$\mathring{R}$ denotes the interior of $R$. We refer to ${\left(\partial
R\right)}_{j}$ as the $j$-th side of the boundary of $R$. Note that $\partial
R$ is the union of its six sides, i.e., $\partial
R={\cup}_{j=1}^{6}{\left(\partial R\right)}_{j}$. Although a boundary point
$y$ can belong to all six sides of $R$, it will always be evident from our
subsequent arguments which side needs to be considered. Subject to this
clarification, we define the discrete derivative in the outward normal
direction ${e}_{j}$, $j=1,\dots,6$ as follows:

$\mathcal{T}u\left(y\right)=u\left(y\right)-u\left(y-{e}_{j}\right)\,,\quad
y\in(\partial R)_{j}\,,$ (3.1)

where ${e}_{3}={e}_{1}-{e}_{2}$, ${e}_{4}=-{e}_{1}$, ${e}_{5}=-{e}_{2}$ and
${e}_{6}=-{e}_{3}$.

Let ${H}_{0}=\left\\{\left(0,0\right)\right\\}$, and for each
$N\in\mathbb{N}$, define ${H}_{N}$ by the recurrence formula

${H}_{N}\coloneqq\bigcup_{x\in{H}_{N-1}}{F}_{x}\,,$

here $\mathring{H}_{N}\coloneqq{H}_{N-1}$ and ${\left(\partial
H\right)}_{N}\coloneqq{H}_{N}\setminus\mathring{H}_{N}$.

Let $R$ be a finite region, and let $R=\bigcup_{x\in\mathring{R}}{F}_{x}$ be a
representation of $R$. Then we have discrete analogues of both Green’s first
and second identities, given by:

$\sum_{x\in\mathring{R}}\left({\nabla}_{d}^{+}u\left(x\right)\cdot{\nabla}_{d}^{+}v\left(x\right)+{\nabla}_{d}^{-}u\left(x\right)\cdot{\nabla}_{d}^{-}v\left(x\right)+u\left(x\right){\Delta}_{d}v\left(x\right)\right)=\sum_{y\in\partial
R}u\left(y\right)\mathcal{T}v\left(y\right)\,,$ (3.2)

and

$\sum_{x\in\mathring{R}}\left(u\left(x\right){\Delta}_{d}v\left(x\right)-v\left(x\right){\Delta}_{d}u\left(x\right)\right)=\sum_{y\in\partial
R}\left(u\left(y\right)\mathcal{T}v\left(y\right)-v\left(y\right)\mathcal{T}u\left(y\right)\right)\,,$
(3.3)

respectively. Here, ${\nabla}_{d}^{+}$ and ${\nabla}_{d}^{-}$ are defined as
follows:

${\nabla}_{d}^{+}\coloneqq\begin{pmatrix}u\left(x+{e}_{1}\right)-u\left(x\right)\\\
u\left(x+{e}_{2}\right)-u\left(x\right)\\\
u\left(x+{e}_{3}\right)-u\left(x\right)\end{pmatrix}\,,$

and

${\nabla}_{d}^{-}\coloneqq\begin{pmatrix}u\left(x-{e}_{1}\right)-u\left(x\right)\\\
u\left(x-{e}_{2}\right)-u\left(x\right)\\\
u\left(x-{e}_{3}\right)-u\left(x\right)\end{pmatrix}\,.$

Let $\mathcal{G}\left(x;y\right)$ denote Green’s function for the discrete
Helmholtz equation (2.6a) centred at $y$ and evaluated at $x$. Then, the
function $\mathcal{G}\left(x;y\right)$ satisfies the following equation:

$\left(\Delta_{d}+k^{2}\right)\mathcal{G}\left(x;y\right)=\delta_{x,y}\,,$
(3.4)

where $\delta_{x,y}$ represents the Kronecker delta. To simplify the notation,
we denote $\mathcal{G}\left(x;0\right)$ as $\mathcal{G}\left(x\right)$. It is
important to note that $\mathcal{G}\left(x;y\right)$ can be represented as
$\mathcal{G}\left(x-y\right)$.

The lattice Green’s function $\mathcal{G}$ is a well-known concept and has
been extensively studied in the literature, see: [10, 19, 23, 31]. By
employing the discrete Fourier and the inverse Fourier transforms, we can
obtain the following result:

$\mathcal{G}\left(x\right)=\frac{1}{4\pi^{2}}\int\limits_{-\pi}^{\pi}\int\limits_{-\pi}^{\pi}\frac{e^{\imath\left(x\cdot\xi\right)}}{\sigma\left(\xi;k^{2}\right)}\,\mathrm{d}\xi\,,\quad\xi=\left(\xi_{1},\xi_{2}\right)\,,$
(3.5)

where

$\displaystyle\sigma\left(\xi;k^{2}\right)$
$\displaystyle=e^{\imath\xi_{1}}+e^{-\imath\xi_{1}}+e^{\imath\xi_{2}}+e^{-\imath\xi_{2}}+e^{\imath\xi_{1}}e^{-\imath\xi_{2}}+e^{-\imath\xi_{1}}e^{\imath\xi_{2}}-6+k^{2}$
$\displaystyle=k^{2}-6+2\cos\xi_{1}+2\cos\xi_{2}+2\cos\left(\xi_{1}-\xi_{2}\right)$
$\displaystyle=k^{2}-4\left({\sin}^{2}\frac{\xi_{1}}{2}+{\sin}^{2}\frac{\xi_{2}}{2}+{\sin}^{2}\frac{\xi_{1}-\xi_{2}}{2}\right)\,.$

The lattice Green’s function $\mathcal{G}$ is a well-established concept when
$k^{2}$ belongs to the set $\mathbb{C}\setminus\left[0,9\right]$ (cf., e.g.,
[12]). Notably, if $k^{2}\in\mathbb{C}\setminus\left[0,9\right]$, then
$\sigma\neq 0$ and $\mathcal{G}$ in (3.5) possesses a well-defined value. In
such cases, the function $\mathcal{G}\left(x\right)$ exhibits an exponential
decay as $\left|x\right|\to\infty$.

When $k\in\left(0,2\sqrt{2}\right)$, we define the lattice Green’s function as
the point-wise limit of the following expression:

$\left(R_{k^{2}+\imath\varepsilon}\delta_{x,0}\right)\left(x\right)\coloneqq\frac{1}{4\pi^{2}}\int\limits_{-\pi}^{\pi}\int\limits_{-\pi}^{\pi}\frac{e^{\imath
x\cdot\xi}\,\mathrm{d}\xi_{1}\mathrm{d}\xi_{2}}{\sigma\left(\xi;k^{2}+\imath\varepsilon\right)}\,,$
(3.6)

as $k^{2}+\imath\varepsilon$ approaches $k^{2}+\imath 0$, we denote the
resulting limit by $\mathcal{G}\left(x\right)$, i.e.,
$\mathcal{G}\left(x\right)=\left(R_{k^{2}+\imath
0}\delta_{x,0}\right)\left(x\right)$, cf., [15]. It is noteworthy that
$\mathcal{G}\left(x\right)$ satisfies the equation (3.4) and the following
equalities for all $x=\left(x_{1},x_{2}\right)\in\mathbb{Z}^{2}$ (see [17]):

$\mathcal{G}\left(x_{1},x_{2}\right)=\mathcal{G}\left(x_{2},x_{1}\right)=\mathcal{G}\left(-x_{1},-x_{2}\right)=\mathcal{G}\left(x_{1}+x_{2},-x_{2}\right)\,.$
(3.7)

Let us reconsider the expression for the lattice Green’s function (3.5) with a
different formulation. By employing the sum-to-product identities and double-
angle formula for cosine within $\sigma\left(\xi;k^{2}\right)$, we derive:

$\sigma\left(\xi;k^{2}\right)={k}^{2}-8+4\cos\frac{\xi_{1}-\xi_{2}}{2}\left(\cos\frac{\xi_{1}+\xi_{2}}{2}+\cos\frac{\xi_{1}-\xi_{2}}{2}\right)\,.$

If we introduce a change of variables as follows:

$\eta_{1}=\frac{\xi_{1}+\xi_{2}}{2}\in\left[-\pi,\pi\right]\,,\quad\eta_{2}=\frac{\xi_{1}-\xi_{2}}{2}\in\left[-\pi,\pi\right]\,,$

or equivalently,

$\xi_{1}=\eta_{1}+\eta_{2}\,,\quad\xi_{2}=\eta_{1}-\eta_{2}\,,$

we can define the function

$\widetilde{\sigma}\left(\eta;k^{2}\right)={k}^{2}-8+4\cos\eta_{2}\left(\cos\eta_{1}+\cos\eta_{2}\right)\,,\quad\eta=\left({\eta}_{1},{\eta}_{2}\right)\,.$

In this change of variables, the integration domain
$\left[-\pi,\pi\right]\times\left[-\pi,\pi\right]$ is rotated by $\pi/2$
around the origin and compressed by a factor of $1/\sqrt{2}$. Consequently,
the lattice Green’s function, denoted by (3.5), can be expressed in terms of
new variables ${\eta}_{1}$ and ${\eta}_{2}$ as follows

$\mathcal{G}\left({x}_{1},{x}_{2}\right)=\frac{1}{4\pi^{2}}\int\limits_{-\pi}^{\pi}\int\limits_{-\pi}^{\pi}\frac{e^{\imath\left(x_{1}\left(\eta_{1}+\eta_{2}\right)+x_{2}\left(\eta_{1}-\eta_{2}\right)\right)}}{\widetilde{\sigma}\left({\eta}_{1},{\eta}_{2};k^{2}\right)}\,\mathrm{d}\eta_{2}\mathrm{d}\eta_{1}\,.$

######  Theorem 3.1 (see [15]).

Consider a finite region $R$ and a function $u:R\to\mathbb{C}$. For any point,
$x\in\mathring{R}$, a discrete Green’s representation formula can be
established

$u\left(x\right)=\sum_{y\in\partial
R}\bigl{(}u\left(y\right)\mathcal{T}\mathcal{G}\left(x-y\bigr{)}-\mathcal{G}\left(x-y\right)\mathcal{T}u\left(y\right)\right)+\sum_{y\in\mathring{R}}\mathcal{G}\left(x-y\right)\left(\Delta_{d}+k^{2}\right)u\left(y\right)\,.$

Furthermore, if $u$ satisfies the discrete Helmholtz equation

$\left(\Delta_{d}+k^{2}\right)u\left(x\right)=0\,,\quad\text{in}\leavevmode\nobreak\
\mathring{R},$

then the following representation formula holds:

$u\left(x\right)=\sum_{y\in\partial
R}\bigl{(}u\left(y\right)\mathcal{T}\mathcal{G}\left(x-y\right)-\mathcal{G}\left(x-y\right)\mathcal{T}u\left(y\right)\bigr{)}\,.$
(3.8)

Recall that, here $\partial R$ refers to the boundary of a finite region $R$.

Moreover, we must apply the concept of a radiation condition to the discrete
Helmholtz operators. It should be emphasized that in the case where
$k^{2}\in\left(0,8\right)$, an extra condition at infinity is necessary (cf.,
[27]). Specifically, a function $u:\Omega\to\mathbb{C}$ satisfies the
radiation condition at infinity if (see [15])

$\left\\{\begin{aligned}
u\left(x\right)&=\mathcal{O}\left({\left|x\right|}^{-\frac{1}{2}}\right)\,,\\\
u\left(x+{e}_{j}\right)&=e^{\imath{\xi}_{j}^{*}\left(\alpha,k\right)}u\left(x\right)+\mathcal{O}\left({\left|x\right|}^{-\frac{3}{2}}\right)\,,\quad
j=1,2\,.\end{aligned}\right.$ (3.9)

Here, the term remaining after the decay is uniform in all directions
$x\,/\left|x\right|$, where $x$ is characterized by
${x}_{1}=\left|x\right|\cos\alpha$ and ${x}_{2}=\left|x\right|\sin\alpha$,
$0\leq\alpha<2\pi$. Here, ${\xi}_{j}^{*}\left(\alpha,k\right)$ is the $j$-th
coordinate of the point ${\xi}^{*}\left(\alpha,k\right)$, where
${\xi}^{*}\left(\alpha,k\right)={\xi}^{*}=\left({\xi}_{1}^{*},{\xi}_{2}^{*}\right)$
is the unique solution to the following system of equations:

$\displaystyle
2\zeta\left(\sin\xi_{1}+\sin\left(\xi_{1}-\xi_{2}\right)\right)$
$\displaystyle=\cos\alpha\,,$ $\displaystyle
2\zeta\left(\sin\xi_{2}-\sin\left(\xi_{1}-\xi_{2}\right)\right)$
$\displaystyle=\sin\alpha\,,$ $\displaystyle
k^{2}-6+2\cos\xi_{1}+2\cos\xi_{2}+2\cos\left(\xi_{1}-\xi_{2}\right)$
$\displaystyle=0\,,$

where $\zeta$ is a positive constant (cf., [15]).

######  Definition 3.2.

Let $k$ be an element of the open interval $\left(0,2\sqrt{2}\right)$.
Consider a solution $u$ to the discrete Helmholtz equation (2.6a), which is a
second-order difference equation used to model wave phenomena in a triangular
lattice. We say that $u$ is a radiating solution if it satisfies the radiation
condition (3.9), which ensures that the solution decays at infinity and
represents a physically realistic wave field.

Observe that the second condition presented in (3.9) can be formulated in a
subsequent manner:

$u\left(x\right)-u\left(x+e_{j}\right)=\left(1-e^{\imath{\xi}_{j}^{*}\left(\alpha,k\right)}\right)u\left(x\right)+\mathcal{O}\left({\left|x\right|}^{-\frac{3}{2}}\right)\,,\quad\left|x\right|\to\infty\,,$

and

$u\left(x+e_{j}\right)-u\left(x\right)=\left(1-e^{-\imath{\xi}_{j}^{*}\left(\alpha,k\right)}\right)u\left(x+e_{j}\right)+\mathcal{O}\left({\left|x\right|}^{-\frac{3}{2}}\right)\,,\quad\left|x\right|\to\infty\,.$

######  Lemma 3.3 (For the detailed proof see [16]).

Suppose that $k^{2}$ belongs to the interval $\left(0,8\right)$, and consider
a function $u$ that satisfies the radiation condition at infinity, as
expressed by equation (3.9). For any point $y$ on the boundary of $H_{N}$, we
have the following asymptotic expansion as $N$ tends to infinity:

$\mathcal{T}u\left(y\right)=\zeta\left(y,k\right)u\left(y\right)+\mathcal{O}\left({\left|y\right|}^{-\frac{3}{2}}\right)\,,\quad
N\to\infty\,,$

where $\zeta\left(y,k\right)$ is a complex-valued function satisfying
$\mathrm{\Im m}\,\zeta\left(y,k\right)>0$.

Consider a fixed point $x\in\mathbb{Z}^{2}$, and let $y$ be any point on the
boundary of $H_{N}$. The radiation conditions expressed by equation (3.9)
imply that

$\sum_{y\in\partial
H_{N}}\left(u\left(y\right)\mathcal{T}\mathcal{G}\left(x-y\right)-\mathcal{G}\left(x-y\right)\mathcal{T}u\left(y\right)\right)\to
0\,,\quad N\to\infty\,.$

Undoubtedly, as an example, it can be observed that as
$\left|y\right|\to\infty$, $\alpha\left(y-x\right)$ converges to
$\alpha=\alpha\left(y\right)$ and $\zeta\left(y-x,k\right)$ converges to
$\zeta\left(y,k\right)$, implying that for $N$ sufficiently large, the ensuing
result can be obtained (see [16]):

$\displaystyle u\left(y\right)\mathcal{T}\mathcal{G}\left(x-y\right)$
$\displaystyle-\mathcal{G}\left(x-y\right)\mathcal{T}u\left(y\right)$
$\displaystyle=u\left(y\right)\cdot\mathcal{O}\left({N}^{-\frac{3}{2}}\right)+\mathcal{G}\left(x-y\right)\cdot\mathcal{O}\left({N}^{-\frac{3}{2}}\right)=\mathcal{O}\left({N}^{-2}\right)\,.$

Moreover, by applying Theorem 3.1 to $\Omega\cap{H}_{N}$, where
$N\in\mathbb{N}$ is sufficiently large, and subsequently taking the limit as
$N\to\infty$, the ensuing Green’s formula for a radiating solution $u$ of the
discrete Helmholtz equation (2.6a) can be derived:

$u\left(x\right)=\sum_{y\in\partial\Omega}\left(u\left(y\right)\mathcal{T}\mathcal{G}\left(x-y\right)-\mathcal{G}\left(x-y\right)\mathcal{T}u\left(y\right)\right)\,.$
(3.10)

Based on the representation formula (3.10) and using the results established
in [15], it can be deduced that any radiating solution $u$ of the discrete
Helmholtz equation (2.6a) possesses the ensuing asymptotic expansion:

$u\left(x\right)=-\frac{{e}^{\imath\mu\left(\alpha,k\right)\left|x\right|}}{{\left|x\right|}^{\frac{1}{2}}}\left\\{{u}_{\infty}\left(\hat{x}\right)+\mathcal{O}\left(\frac{1}{\left|x\right|}\right)\right\\}\,,\quad\left|x\right|\to\infty\,.$
(3.11)

Furthermore, let us define
$\mu\left(\alpha,k\right)\coloneqq{\xi}^{*}\left(\alpha,k\right)\cdot\hat{x}$,
where $\hat{x}\coloneqq x\,/\left|x\right|$ and
${\xi}^{*}\left(\alpha,k\right)=\left({\xi}_{1}^{*}\left(\alpha,k\right),{\xi}_{2}^{*}\left(\alpha,k\right)\right)$.
The far-field pattern of $u$ denoted as ${u}_{\infty}\left(\hat{x}\right)$,
can be expressed through the employing of the formula $\left(12\right)$
introduced in [15].

Now, based on the previous results, we are ready to state the following
proposition:

######  Theorem 3.4.

The _Problem_ $\mathcal{P}_{\mathrm{ext}}$ possesses, at most, a single
solution that exhibits radiation properties.

  * Proof.

Demonstrating that the corresponding homogeneous problem has only the trivial
solution is a sufficient condition to establish the result.

Let $\Omega\cap{H}_{N}$ denote the intersection of the domain $\Omega$ with
the finite rectangular grid ${H}_{N}$. In references [15, 16], the discrete
version of Green’s first identity, applied in $\Omega\cap{H}_{N}$, is
expressed as follows:

$\sum_{x\in\Omega\cap{H}_{N}}\left({\nabla}_{d}^{+}u\left(x\right)\cdot{\nabla}_{d}^{+}v\left(x\right)+{\nabla}_{d}^{-}u\left(x\right)\cdot{\nabla}_{d}^{+}v\left(x\right)+u\left(x\right){\Delta}_{d}v\left(x\right)\right)=\sum_{y\in\partial{H}_{N}}u\left(y\right)\mathcal{T}v\left(y\right)\,.$
(3.12)

Setting $v\coloneqq\overline{u}$ in (3.12), we obtain the following result:

$\sum_{x\in\Omega\cap{H}_{N}}\left({\left|{\nabla}_{d}^{+}u\left(x\right)\right|}^{2}+{\left|{\nabla}_{d}^{-}u\left(x\right)\right|}^{2}-{k}^{2}{\left|u(x)\right|}^{2}\right)=\sum_{y\in\partial{H}_{N}}u\left(y\right)\mathcal{T}\overline{u}\left(y\right)\,.$
(3.13)

By applying Lemma 3.3, we can reformulate equality (3.13) in the following
manner

$\sum_{x\in\Omega\cap{H}_{N}}\left({\left|{\nabla}_{d}^{+}u\left(x\right)\right|}^{2}+{\left|{\nabla}_{d}^{-}u\left(x\right)\right|}^{2}-{k}^{2}{\left|u(x)\right|}^{2}\right)=\sum_{y\in\partial{H}_{N}}\overline{\zeta}\left(y,k\right){\left|u\left(y\right)\right|}^{2}+\mathcal{O}\left({N}^{-1}\right)\,.$

By taking the imaginary part of the final identity and then letting $N$ tend
towards infinity, we obtain the limit

$\sum_{y\in\partial{H}_{N}}{\left|u\left(y\right)\right|}^{2}\to 0\,,\quad
N\to\infty\,.$

Furthermore, employing a Rellich-type theorem [13, 14, 1], we can assert that
for sufficiently large $N$, the function $u$ vanishes outside of ${H}_{N}$.
Since $\Omega$ satisfies the cone condition, we can apply the unique
continuation property [1, Theorem 5.7], which implies that $u$ is identically
zero throughout $\Omega$. ∎

##  4 Difference potentials and the existence of a solution

For any given function $\varphi:\partial R\to\mathbb{C}$, the difference
single-layer and double-layer potentials are defined as:

$V\varphi\left(x\right)=\sum_{y\in\partial
R}\mathcal{G}\left(x-y\right)\varphi\left(y\right)\,,\quad\text{for all}\
x\in\mathbb{Z}^{2}\,,$ (4.1)

and

$W\varphi\left(x\right)=\sum_{y\in\partial
R}\bigr{(}\mathcal{T}\mathcal{G}\left(x-y\right)+\delta_{x,y}\bigl{)}\varphi\left(y\right)\,,\quad\text{for
all}\ x\in\mathbb{Z}^{2}\,,$ (4.2)

respectively. Since $\delta_{x,y}=0$ for every $x\in\mathring{R}$ and
$y\in\partial R$, expression (3.10) can be written as:

$u\left(x\right)=Wu\left(x\right)-V\left(\mathcal{T}u\right)\left(x\right)\,,\quad
x\in\mathring{\Omega}\,.$

The role of the summand $\delta_{x,y}$ is clarified by the following result.

######  Lemma 4.1.

For every $x\in\mathring{R}$, we have:

$\left({\Delta}_{d}+k^{2}\right)V\varphi\left(x\right)=0\,,\quad\text{and}\quad\left({\Delta}_{d}+k^{2}\right)W\varphi\left(x\right)=0\,.$

This indicates that both the difference single-layer potential
$V\varphi\left(x\right)$ and the double-layer potential
$W\varphi\left(x\right)$ satisfy the discrete Helmholtz equation with
wavenumber $k^{2}$.

  * Proof.

For the difference single-layer potential (4.1), we can obtain the following
result:

$\left({\Delta}_{d}+k^{2}\right)V\varphi\left(x\right)=\sum_{y\in\partial
R}\left[\left({\Delta}_{d}+k^{2}\right)\mathcal{G}\left(x-y\right)\right]\varphi\left(y\right)=\sum_{y\in\partial
R}\delta_{x,y}\varphi\left(y\right)=0\,,\quad\forall x\notin\partial R\,.$

Likewise, if $x\in\mathring{R}$ and ${F}_{x}\cap\partial R=\varnothing$, we
can establish the outcome for the difference double-layer potential (4.2).
Specifically, we have
$\left({\Delta}_{d}+k^{2}\right)\mathcal{T}\mathcal{G}\left(x-y\right)=0$ and
$\left({\Delta}_{d}+k^{2}\right)\delta_{x,y}=0$. Thus, it suffices to focus on
the scenario where $x\in\mathring{R}$ and ${F}_{x}\cap\partial
R\neq\varnothing$.

Let $y\in{F}_{x}\cap{\left(\partial R\right)}_{j}$ for $j=1,\dots,6$. This
implies that $y-{e}_{j}=x$. Applying the properties ${e}_{4}=-{e}_{1}$,
${e}_{5}=-{e}_{2}$ and ${e}_{6}=-{e}_{3}$, it follows from the definitions of
the discrete Laplace operator and Kronecker delta that:

$\displaystyle\left({\Delta}_{d}+k^{2}\right)\delta_{x,y}$
$\displaystyle=\sum_{i=1}^{3}{\left(\delta_{x+{e}_{i},y}+\delta_{x-{e}_{i},y}\right)}-\left(6-{k}^{2}\right)\delta_{x,y}$
$\displaystyle=\sum_{i=1}^{3}{\left(\delta_{x+{e}_{i},x+{e}_{j}}+\delta_{x-{e}_{i},x+{e}_{j}}\right)}-\left(6-{k}^{2}\right)\delta_{x,x+{e}_{j}}=\delta_{x+{e}_{j},x+{e}_{j}}=1\,.$

By considering the definitions of the difference double-layer potential $W$
and the operator $\mathcal{T}$ given by (4.2) and (3.1), respectively, and
employing the previous equality, we can establish the following expression for
the second equality in Lemma 4.1

$\displaystyle\left({\Delta}_{d}+k^{2}\right)\left(\mathcal{T}\mathcal{G}\left(x-y\right)+\delta_{x,y}\right)$
$\displaystyle=\left({\Delta}_{d}+k^{2}\right)\left(\mathcal{G}\left(x;y\right)-\mathcal{G}\left(x;y-{e}_{j}\right)\right)+\left({\Delta}_{d}+k^{2}\right)\delta_{x,y}$
$\displaystyle=\delta_{x,y}-\delta_{x,y-{e}_{j}}+\delta_{y+{e}_{j},y}$
$\displaystyle=0-1+1=0\,.$

Therefore, we can conclude that
$\left({\Delta}_{d}+k^{2}\right)W\varphi\left(x\right)=0$, $x\in\mathring{R}$.
∎

As a consequence of Lemma 4.1, it can be inferred that

$V\varphi\left(x\right)=\sum_{y\in\partial\Omega}\mathcal{G}\left(x-y\right)\varphi\left(y\right)\,,\quad\text{for}\
x\in\mathring{\Omega}\,,$

and

$W\varphi\left(x\right)=\sum_{y\in\partial\Omega}\left(\mathcal{T}\mathcal{G}\left(x-y\right)+\delta_{x,y}\right)\varphi\left(y\right)\,,\quad\text{for}\
x\in\mathring{\Omega}\,,$

are radiating solutions to equation (2.6a) for any function
$\varphi:\partial\Omega\to\mathbb{C}$.

Notice that, by our assumption on $\partial\Omega$ we have
$\Omega=\mathbb{Z}^{2}$. Consider a point $y_{i}$ that intersects several
sides of $\partial\Omega$. To reduce the number of numerical computations, we
shall choose and fix only one side of the boundary. Let $m$ denote the number
of points on $\partial\Omega$. We can represent $\partial\Omega$ as a sequence
of points $y_{1},y_{2},\ldots,y_{m}$ such that $y_{i}=y_{j}$ if and only if
$i=j$ for all $1\leq i,j\leq m$. This allows us to express the difference
potential $\widetilde{V}$ in (4.1) as a single summand associated with the
boundary point $y_{i}$.

Furthermore, given a function $f$ defined on $\partial\Omega$, we can form a
vector $F={\left(f_{1},\dots,f_{m}\right)}^{\top}$, where $f_{i}\coloneqq
f\left(y_{i}\right)$ for all $1\leq i\leq m$. Similarly, for an unknown
function $\varphi$ defined on $\partial\Omega$, we can write
$\Phi={\left(\varphi_{1},\dots,\varphi_{m}\right)}^{\top}$, where
$\varphi_{i}\coloneqq\varphi\left(y_{i}\right)$ for all $1\leq i\leq m$.

We seek a solution to the _Problem_ $\mathcal{P}_{\mathrm{ext}}$ in the
following form:

$u\left(x\right)=\widetilde{V}\varphi\left(x\right)=\sum_{i=1}^{m}\mathcal{G}\left(x-{y}_{i}\right){\varphi}_{i}\,,\quad
x\in\mathring{\Omega}\,.$ (4.3)

As demonstrated in the proof of Lemma 4.1, it is possible to establish that
$u$ given by (4.3) is a radiating solution to equation (2.6a). The solution
must also satisfy the boundary conditions specified in (2.6b). Consequently,
(2.6b) gives rise to the following linear system of boundary equations:

$\mathcal{H}\Phi=F\,,$ (4.4)

here

$\mathcal{H}=\begin{pmatrix}\mathcal{G}\left({y}_{1}-{y}_{1}\right)&\mathcal{G}\left({y}_{1}-{y}_{2}\right)&\mathcal{G}\left({y}_{1}-{y}_{3}\right)&\cdots&\mathcal{G}\left({y}_{1}-{y}_{m}\right)\\\
\mathcal{G}\left({y}_{2}-{y}_{1}\right)&\mathcal{G}\left({y}_{2}-{y}_{2}\right)&\mathcal{G}\left({y}_{2}-{y}_{3}\right)&\cdots&\mathcal{G}\left({y}_{2}-{y}_{m}\right)\\\
\mathcal{G}\left({y}_{3}-{y}_{1}\right)&\mathcal{G}\left({y}_{3}-{y}_{2}\right)&\mathcal{G}\left({y}_{3}-{y}_{3}\right)&\cdots&\mathcal{G}\left({y}_{3}-{y}_{m}\right)\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\mathcal{G}\left({y}_{m}-{y}_{1}\right)&\mathcal{G}\left({y}_{m}-{y}_{2}\right)&\mathcal{G}\left({y}_{m}-{y}_{3}\right)&\cdots&\mathcal{G}\left({y}_{m}-{y}_{m}\right)\end{pmatrix}\,.$

######  Lemma 4.2.

The linear system of boundary equations (4.4) admits a unique solution.

  * Proof.

By virtue of the Rouché-Capelli theorem, the unique solvability of the linear
system of boundary equations (4.4) is equivalent to the non-existence of non-
trivial solutions to the associated homogeneous system, given by

$\mathcal{H}\Phi=0\,.$ (4.5)

Suppose
$\Phi^{*}={\left({\varphi}_{1}^{*},\dots,{\varphi}_{m}^{*}\right)}^{\top}$ is
a solution to (4.5). Substituting the solution ${\varphi}^{*}$ into (4.3), we
obtain

$u\left(x\right)=\widetilde{V}{\varphi}^{*}\left(x\right)=\sum_{i=1}^{m}\mathcal{G}\left(x-{y}_{i}\right){\varphi}_{i}^{*}\,,\quad
x\in\mathring{\Omega}\,.$

Since ${\varphi}^{*}$ is a solution to the homogeneous problem (4.5) and
satisfies the radiation condition (3.9), it follows from Theorem 3.4 that
$u\equiv 0$ in $\Omega$. Moreover, at each boundary point
${y}_{i}\in\partial\Omega$, we have

$0={\left({\Delta}_{d}+{k}^{2}\right)u(x)\Big{|}}_{x={y}_{i}}=\left({\Delta}_{d}+{k}^{2}\right)\widetilde{V}{\varphi}^{*}(x)\Big{|}_{x={y}_{i}}=\sum_{i=1}^{m}\delta_{x,y_{i}}{\varphi}^{*}_{i}\Big{|}_{x={y}_{i}}={\varphi}_{i}^{*}\,,\quad\text{for
all}\ 1\leq i\leq m\,.$

Therefore, all solutions to the homogeneous system (4.5) are trivial, and the
linear system of boundary equations (4.4) has a unique solution. ∎

######  Theorem 4.3.

There exists a unique radiating solution to _Problem_
$\mathcal{P}_{\mathrm{ext}}$, which can be expressed as (4.3), where $\Phi$ is
the unique solution to the system of linear equations (4.4).

##  5 Results of numerical computations

###  5.1 Description of the procedure for computing the lattice Green’s
function

The main challenge in numerically evaluating solutions to (4.4) lies in
computing the lattice Green’s function (3.5). To address this, we employ the
method devised by Berciu and Cook [2], although various approaches have been
developed (e.g., refer to [12, 24]). Applying 8-fold symmetry, it is
sufficient to calculate the lattice Green’s function
$\mathcal{G}\left(i,j\right)$ for $i\geq j\geq 0$. Following the approach in
[2], we introduce the vectors
$\mathcal{V}_{2p}={\left(\mathcal{G}\left(2p,0\right),\mathcal{G}\left(2p-1,1\right),\ldots,\mathcal{G}\left(p,p\right)\right)}^{\top}$
and
$\mathcal{V}_{2p+1}={\left(\mathcal{G}\left(2p+1,0\right),\mathcal{G}\left(2p,1\right),\ldots,\mathcal{G}\left(p+1,p\right)\right)}^{\top}$.
These vectors collect all distinct Green’s functions
$\mathcal{G}\left(i,j\right)$ with “Manhattan distances”
$\left|i\right|+\left|j\right|$ of $2p$ and $2p+1$, respectively. For any
Manhattan distance greater than $1$, the equation

$\left(\Delta_{d}+k^{2}\right)\mathcal{G}\left(x\right)=\delta_{x,0}\,,$ (5.1)

can be expressed in matrix form as

${\gamma}_{n}\left(k\right){\mathcal{V}}_{n}={\alpha}_{n}\left(k\right){\mathcal{V}}_{n-1}+{\beta}_{n}\left(k\right){\mathcal{V}}_{n+1}\,,$
(5.2)

where the matrices ${\alpha}_{n}\left(k\right)$, ${\beta}_{n}\left(k\right)$,
and ${\gamma}_{n}\left(k\right)$ feature the property of sparsity. Notably,
only the dimensions of these matrices depend on $n$ (For details on the
dimensions of these matrices, refer to Table 2, and for the assembly of these
sparse matrices, see Appendix A). As demonstrated in [2], for any $n\geq 1$,
we obtain

${\mathcal{V}}_{n}={A}_{n}\left(k\right){\mathcal{V}}_{n-1}\,,$ (5.3)

where the matrices ${A}_{n}\left(k\right)$ are defined by inserting (5.3) into
the equation (5.2)

${A}_{n}\left(k\right)={\left[{\gamma}_{n}\left(k\right)-{\beta}_{n}\left(k\right){A}_{n+1}\right]}^{-1}{\alpha}_{n}\left(k\right)\,.$

The matrices can be computed starting from a sufficiently large $N$ with
${A}_{N+1}\left(k\right)=0$. It is important to note that an exceptional case
arises when $k=2$, leading to a failure of this procedure. In such instances,
a more refined “initial guess” is required than ${A}_{N+1}\left(k\right)=0$,
as $\mathrm{det}\left({\gamma}_{n}\left(k\right)\right)=0$, and the matrix
${\gamma}_{n}\left(k\right)-{\beta}_{n}\left(k\right){A}_{n+1}$ is non-
invertible. To address this scenario, we set
$\tilde{k}^{2}={k}^{2}+\imath\varepsilon$, where $\varepsilon$ is sufficiently
small. Assuming, for significantly large Manhattan distances
$N=\left|i\right|+\left|j\right|$ where $N+1=2p$, $p\in\mathbb{Z}$, there
exists an asymptotic dependence:
$\mathcal{G}\left(i+1,j\right)\approx\mathcal{G}\left(i,j+1\right)=\lambda\left(\tilde{k}\right)\mathcal{G}\left(i,j\right)$,
with $\left|\lambda\left(\tilde{k}\right)\right|<1$. This implies that all
propagators exhibit a uniform decrease as $N$ grows. Upon considering the
asymptotic dependence in equation (5.1), the ensuing quadratic equation is
derived: $2{\lambda}^{2}+\left(\tilde{k}^{2}-4\right)\lambda+2=0$, yielding
the immediate conclusion of the existence of a physical solution
$\left|\lambda\left(\tilde{k}\right)\right|<1$ (cf. [2]). The resulting matrix
${A}_{N+1}\left(k\right)$ has a dimension $\left(p+1\right)\times p$, and its
entries are defined as follows:
${A}_{N+1}\left(\tilde{k}\right)\mid_{i,i}={A}_{N+1}\left(\tilde{k}\right)\mid_{i,i-1}=\lambda/2$,
$i=\overline{2,p}$. Additionally,
${A}_{N+1}\left(\tilde{k}\right)\mid_{1,1}={A}_{N+1}\left(\tilde{k}\right)\mid_{p+1,p}={\lambda}$,
and all other matrix elements remain zero.

Once ${A}_{n}\left(k\right)$ are known, we can express
${\mathcal{V}}_{n}={A}_{n}\left(k\right)\cdots{A}_{1}\left(k\right){\mathcal{V}}_{0}$,
where ${\mathcal{V}}_{0}=\mathcal{G}\left(0,0\right)$. Specifically,
${\mathcal{V}}_{1}=\mathcal{G}\left(1,0\right)={A}_{1}\left(k\right)\mathcal{G}\left(0,0\right)$.
By inserting this into equation (5.1) and considering property (3.7), we
obtain the expression
$6\mathcal{G}\left(1,0\right)-\left(6-{k}^{2}\right)\mathcal{G}\left(0,0\right)=1$.
From this, we determine that
$\mathcal{G}\left(0,0\right)=1/\left[6{A}_{1}\left(k\right)-6+{k}^{2}\right]$.
This concludes the computation employing only elementary operations and devoid
of integrals. Additionally, it is important to note another meaningful
advantage of this method. The matrices ${A}_{n}\left(k\right)$ are computed by
descending from asymptotically large Manhattan distances. As they propagate
towards diminished Manhattan distances, it provides us with the physical
solution.

Let us proceed with developing another approach to determine the “initial
guess” for computing the lattice Green’s function using the described
procedure, particularly when $k$ pertains to the set of real numbers,
specifically $k\in\left(0,2\sqrt{2}\right)$. Assuming asymptotically large
Manhattan distances $\left|i\right|+\left|j\right|\gg 1$, the following
heuristic relation holds:

$\mathcal{G}\left(r,s\right)={\lambda}^{\left(r+s\right)-\left(i+j\right)}q\left({r+s},{i+j}\right){h}\mathcal{G}\left(i,j\right)\,,\quad
q\left({r+s},{i+j}\right)=\frac{i+j}{r+s}\,.$

The role of the parameters $\lambda$ and $h\neq 0$ are defined in the
subsequent discussion. According to (5.3), we express
${\mathcal{V}}_{n+1}={A}_{n+1}\left(k\right){\mathcal{V}}_{n}$, where
$n+1=2p$. Consequently, we obtain the following representations for the
vectors ${\mathcal{V}}_{n+1}$ and ${\mathcal{V}}_{n}$, viz.,
${\mathcal{V}}_{n+1}={\left(\mathcal{G}\left(2p,0\right),\mathcal{G}\left(2p-1,1\right),\ldots,\mathcal{G}\left(p,p\right)\right)}^{\top}$
and
$\mathcal{V}_{n}={\left(\mathcal{G}\left(2p-1,0\right),\mathcal{G}\left(2p-2,1\right),\ldots,\mathcal{G}\left(p,p-1\right)\right)}^{\top}$.
Let us now express the relation between ${\mathcal{V}}_{n+1}$ and
${\mathcal{V}}_{n}$ in matrix-vector form

${\begin{pmatrix}\mathcal{G}\left(2p,0\right)\\\
\mathcal{G}\left(2p-1,1\right)\\\ \mathcal{G}\left(2p-2,2\right)\\\ \vdots\\\
\mathcal{G}\left(p,p\right)\end{pmatrix}}_{\left(p+1\right)\times
1}={\left[{A}_{n+1}\left(k\right)\right]}_{\left(p+1\right)\times
p}{\begin{pmatrix}\mathcal{G}\left(2p-1,0\right)\\\
\mathcal{G}\left(2p-2,1\right)\\\ \vdots\\\
\mathcal{G}\left(p,p-1\right)\end{pmatrix}}_{p\times 1}\,.$

In our cases, we have $i=2p-\ell$ and $j=\ell$, where $\ell$ ranges from $0$
to $p$. For simplicity, we reintroduce the notation
$q\left({2p+r+s},{2p}\right)$ as $q\left(p,{r+s}\right)$, meaning
$q\left({2p+r+s},{2p}\right)=q\left(p,{r+s}\right)$. As a result, we can
easily derive

$\mathcal{G}\left(2p-\ell+r,\ell+s\right)={\lambda}^{r+s}q\left(p,{r+s}\right){h}\mathcal{G}\left(2p-\ell,\ell\right)\,,\quad
q\left(p,{r+s}\right)=\frac{2p}{2p+r+s}\,.$ (5.4)

Due to the equation (5.1), we obtain

$\displaystyle\left(6-{k}^{2}\right)\mathcal{G}\left(2p-\ell,\ell\right)$
$\displaystyle=\mathcal{G}\left(2p-\ell+1,\ell\right)+\mathcal{G}\left(2p-\ell-1,\ell\right)$
$\displaystyle+\mathcal{G}\left(2p-\ell,\ell+1\right)+\mathcal{G}\left(2p-\ell,\ell-1\right)$
$\displaystyle+\mathcal{G}\left(2p-\ell+1,\ell-1\right)+\mathcal{G}\left(2p-\ell-1,\ell+1\right)\,,$

and, considering the assumption for
$\mathcal{G}\left(2p-\ell+r,\ell+s\right)$, we arrive at the following
equation:

$a\left(p,h\right){\lambda}^{2}+b\left(k,h\right)\lambda+c\left(p,h\right)=0\,,$
(5.5)

where

$\displaystyle a\left(p,h\right)$
$\displaystyle=2q\left(p,{1}\right){h}=\frac{4p}{2p+1}{h}\,,$ $\displaystyle
b\left(k,h\right)$
$\displaystyle={k}^{2}-6+2q\left(p,{0}\right){h}={k}^{2}-6+2{h}\,,$
$\displaystyle c\left(p,h\right)$
$\displaystyle=2q\left(p,{-1}\right){h}=\frac{4p}{2p-1}{h}\,.$

The discriminant of the quadratic polynomial (5.5) is given by

$\mathcal{D}_{1}={\left({k}^{2}-6+2{h}\right)}^{2}-\frac{64{p}^{2}}{4{p}^{2}-1}{h}^{2}=-4\frac{12{p}^{2}+1}{4{p}^{2}-1}{h}^{2}+4\left({k}^{2}-6\right){h}+{\left({k}^{2}-6\right)}^{2}\,.$
(5.6)

Let us now calculate the discriminant of the quadratic polynomial (5.6), that
is

$\frac{\mathcal{D}_{2}}{4}=\frac{64{p}^{2}}{4{p}^{2}-1}{\left({k}^{2}-6\right)}^{2}\geq
0\,,\quad p\gg 1\,.$ (5.7)

The solutions to the quadratic polynomial (5.6) are expressed in the following
manner:

$\displaystyle
h=\left\\{\begin{array}[]{ll}\left(6-{k}^{2}\right){\left(\pm\frac{8p}{\sqrt{4{p}^{2}-1}}+2\right)}^{-1},&\hbox{\text{for}
$0<{k}^{2}<6$,}\\\
\left({k}^{2}-6\right){\left(\pm\frac{8p}{\sqrt{4{p}^{2}-1}}-2\right)}^{-1},&\hbox{\text{for}
$6<{k}^{2}<8$.}\end{array}\right.$ (5.8)

To establish a relationship between $\lambda$ and $h$, it is essential to
align with the behaviour of the solution to the given problem. Specifically,
wave propagation diminishes over sufficiently large Manhattan distances. In
simpler terms, all propagators decrease at a certain rate as $n$ increases. We
express this fact in the following way, observing that due to (5.4), we have
the asymptotic dependence

$\mathcal{G}\left(2p-\ell+1,\ell\right)=\mathcal{G}\left(2p-\ell,\ell+1\right)={\lambda}\frac{2p}{2p+1}{h}\mathcal{G}\left(2p-\ell,\ell\right)\,,\quad
h\neq 0\,.$

Following this point, the subsequent condition must be fulfilled

$\left|\lambda\right|\frac{2p}{2p+1}\left|{h}\right|<1\,,$ (5.9)

this implies that

$\left|h\right|<\frac{1}{\left|{\lambda}\right|}\frac{2p+1}{2p}=\frac{1}{\left|{\lambda}\right|}\left(1+\frac{1}{2p}\right)\,.$
(5.10)

In the specific case where ${k}^{2}=6$, as assumed, and with the additional
condition that $h\neq 0$, it follows from (5.6) that $\mathcal{D}_{1}<0$.
Consequently, the quadratic equation (5.5) possesses two distinct (non-real)
complex roots, which are complex conjugates of each other. In this case, we
derive

$\left|\lambda\right|=\sqrt{\frac{2p+1}{2p-1}}\,.$ (5.11)

Considering the upper bound (5.10) for $\left|\lambda\right|$, along with the
condition $h\neq 0$ and the quality (5.11), we easily obtain the result

$0<\left|{h}\right|<\sqrt{1-\frac{1}{4{p}^{2}}}\,,\quad p\gg 1\,.$

Let us now move forward to estimate $h$ in the case where
${k}^{2}\in\left(0,8\right)\setminus\left\\{6\right\\}$. In this context,
consider the scenario where the discriminant (5.6) of the quadratic polynomial
(5.5) is negative, meaning that $h$ lies outside the roots of the polynomial
(5.6). Consequently, for the polynomial (5.5), we have two distinct complex
roots ${\lambda}_{1}$ and ${\lambda}_{2}$. It is easy to show, for (5.8), that
the following estimate is valid

$\left|h\right|\geq\frac{\left|{k}^{2}-6\right|}{6}\,.$ (5.12)

The same equality is derived for $\left|\lambda\right|$ as (5.11). Considering
this result together with (5.10) and (5.12), we obtain

$\frac{\left|{k}^{2}-6\right|}{6}\leq\left|h\right|<\sqrt{1-\frac{1}{4{p}^{2}}}\,,$

providing that

$p>\frac{3}{k\sqrt{12-{k}^{2}}}\,.$

Before presenting the entries of the matrix ${A}_{n+1}\left(k\right)$, we
first provide a remark regarding $\lambda$, which is the root of the quadratic
polynomial (5.5) and depends on the values of $k$, $p$, and $h$. The entries
of the matrix ${A}_{n+1}\left(k\right)$ are denoted using the following
notation: ${a}_{i,j}^{n+1}\left(k\right)$, where $i=1,2,\ldots,p+1$ and
$j=1,2,\ldots,p$. From the matrix-vector form, we easily identify that non-
zero elements of the matrix ${A}_{n+1}\left(k\right)$ are defined as follows:

${a}_{1,1}^{n+1}\left(k\right)={a}_{p+1,p}^{n+1}\left(k\right)=\frac{2p-1}{2p}{\rho}\,,\quad{a}_{i,i-1}^{n+1}\left(k\right)={a}_{i,i}^{n+1}\left(k\right)=\frac{2p-1}{4p}{\rho}\,,\quad{\rho}=\frac{\lambda}{h}\,,\quad
i=2,3,\ldots,p\,.$

###  5.2 Numerical experiments on computing lattice Green’s functions

Prior to raising the numerical results of the given problem (2.6a)-(2.6b), it
is advisable to provide remarks regarding the numerical computations of the
lattice Green’s function (3.5). The accuracy of the calculation of the exact
solution representation (4.3) mainly depends on the numerical approximation of
the lattice Green’s function. To demonstrate the experimental convergence of
the procedure outlined in Subsec. 5.1 for the numerical computation of the
lattice Green’s function, consider ${p}_{m}={2}^{m}{p}_{0}$ with ${p}_{0}=71$,
and ${N}_{m}={2}^{m+1}{p}_{0}-1$, where $m=1,2,3,4$. For each step $m$, we
define the following vectors:
${\mathcal{V}}_{2i}^{\left(m\right)}={\left({\mathcal{G}}^{\left(m\right)}\left(2i,0\right),{\mathcal{G}}^{\left(m\right)}\left(2i-1,1\right),\ldots,{\mathcal{G}}^{\left(m\right)}\left(i,i\right)\right)}^{\top}$
and
${\mathcal{V}}_{2i+1}^{\left(m\right)}={\left({\mathcal{G}}^{\left(m\right)}\left(2i+1,0\right),{\mathcal{G}}^{\left(m\right)}\left(2i,1\right),\ldots,{\mathcal{G}}^{\left(m\right)}\left(i+1,i\right)\right)}^{\top}$,
where $i=0,1,\ldots,{p}_{m}-1$. It is important to note that, for every step
$m$ starting from $0$, there are a total of ${N}_{m}+1$ such vectors, each
with the dimension of $\left(i+1\right)\times 1$. As $m$ increases, the first
${N}_{0}+1=2{p}_{0}$ vectors are common. To collect the elements of these
common vectors, it is useful to form a matrix array by completing them.
Consider the matrix $\boldsymbol{\mathcal{V}}_{{N}_{0}}^{\left(m\right)}$ of
dimension $\left({N}_{0}+1\right)\times{p}_{0}$, with its entries denoted by
${v}_{i,j}^{\left(m\right)}$. The non-zero elements of this matrix are
precisely defined as
${v}_{i,j}^{\left(m\right)}={\mathcal{G}}^{\left(m\right)}\left(i-1,j-1\right)$,
subject to the constraints $i-j\geq 0$ and $i+j\leq{N}_{0}+2$, which
correspond to ${p}_{0}\left({p}_{0}+1\right)=5112$ entries. All other entries
of $\boldsymbol{\mathcal{V}}_{{N}_{0}}^{\left(m\right)}$ are zero. Let us
consider the element-wise norm of the matrix difference, restricted by $j\leq
i\leq-j+{N}_{0}+2$ and $j=1,2,\ldots,{p}_{0}$, defined as:

${\left\|\boldsymbol{\mathcal{V}}_{{N}_{0}}^{\left(m\right)}-\boldsymbol{\mathcal{V}}_{{N}_{0}}^{\left(m-1\right)}\right\|}_{\max}=\max\limits_{i,j}\left|{v}_{i,j}^{\left(m\right)}-{v}_{i,j}^{\left(m-1\right)}\right|\,,\quad
m=1,2,3,4\,.$

By employing analogous reasoning, we introduce the matrix
$\boldsymbol{\mathcal{V}}_{{N}_{0},{h}}^{\left(m\right)}$, wherein the non-
zero entries correspond to the lattice Green’s function computed using another
choice of the “initial guess” that depends on the parameter $h$. In our
computations, we explored various values of $\left|h\right|$, resulting in
modifications to the coefficients associated with the convergence order.
Nevertheless, the convergence order itself remains consistent with that
observed in all preceding cases. In the considered scenario, we select
$\left|h\right|=1/\left|\lambda\right|=\sqrt{(2{p}_{m}-1)/(2{p}_{m}+1)}$. The
results derived from the described computations for the wavenumber $k=2$ are
summarized in Table 1. Additionally, Figure 3 provides graphical
representations of the real part of the lattice Green’s functions for the
triangular lattice in both transformed and original coordinates, corresponding
to ${N}_{4}=2271$. For illustrative purposes, the transformed coordinates of
the triangular lattice are plotted within the range
$\left[-40,40\right]\times\left[-40,40\right]\subset{\mathbb{Z}}^{2}$.
Meanwhile, the original coordinates are confined to a subset of
${\mathbb{R}}^{2}$ through a suitable inverse transformation from
${\mathbb{Z}}^{2}$ coordinates.

It is important to highlight that all calculations were performed using the
scientific programming language GNU Octave $8.4.0$ on a standard laptop
computer equipped with an AMD Ryzen $5$ $5600$H processor with Radeon Graphics
($12$ CPUs), running at approximately $3.3$GHz, and having $8192$MB RAM. We
encountered limitations preventing further exploration (i.e., for $m\geq 5$).
Specifically, for $m=5$, the computation involves evaluating
${p}_{5}\left({p}_{5}+1\right)=5.164256e{+}6$ lattice Green’s function values,
exceeding the available memory capacity of this machine.

$m$ | ${p}_{m}$ | ${N}_{m}$ | ${\left\|\boldsymbol{\mathcal{V}}_{{N}_{0}}^{\left(m\right)}-\boldsymbol{\mathcal{V}}_{{N}_{0}}^{\left(m-1\right)}\right\|}_{\max}$ | ${\left\|\boldsymbol{\mathcal{V}}_{{N}_{0},{h}}^{\left(m\right)}-\boldsymbol{\mathcal{V}}_{{N}_{0},{h}}^{\left(m-1\right)}\right\|}_{\max}$
---|---|---|---|---
$0$ | $71$ | $141$ | $-$ | $-$
$1$ | $142$ | $283$ | $9.1546e{-}2$ | $9.0821e{-}2$
$2$ | $284$ | $567$ | $8.7608e{-}4$ | $8.2637e{-}4$
$3$ | $568$ | $1135$ | $4.1193e{-}4$ | $4.2530e{-}4$
$4$ | $1136$ | $2271$ | $1.8079e{-}4$ | $1.7728e{-}4$
Table 1: The absolute error associated with the difference between the values
of every two successive lattice Green’s functions for the case of $k=2$.

We conducted numerical experiments, presented in Table 1, to investigate the
behaviour of computed lattice Green’s functions when initiating computations
with various fixed Manhattan distances ${N}_{m}$. Our analysis focused on
observing the maximum absolute differences between lattice Green’s functions
that are common in all considered instances. The numerical results indicate a
stability property when computations start from a considerable distance. As
referred in [2], initiating computations for asymptotically large Manhattan
distances with a refined “initial guess” yields an immediate physical
solution. However, advancing further in the context of the subsequent increase
in Manhattan distance established impossible for us, given the limitations of
the machine resources.

In addition to the prescribed approach for computing lattice Green’s
functions, we have explored various numerical integration techniques,
including the composite trapezoidal rule, composite Simpson’s $1/3$ rule, and
Gaussian quadrature. These methods are employed for the numerical
approximation of the lattice Green’s functions (3.5) with
${k}^{2}+\imath\varepsilon$, which involve double integrals. It is observed
that the integrand in (3.5) exhibits fast oscillations as the magnitudes of
${x}_{1}$ and ${x}_{2}$ increase from the origin. To address this behaviour,
it becomes necessary to reduce the mesh length for both integration variables
significantly. However, this adjustment alone is insufficient to achieve the
desired precision, as the fast oscillations function is approximated.
Additionally, it becomes computationally expensive.

(a) The density plot of $\mathfrak{Re}\left(\mathcal{G}\right)$;
$\left({x}_{1},{x}_{2}\right)\in{\mathbb{Z}}^{2}$.

(b) The density plot of $\mathfrak{Re}\left(\mathcal{G}\right)$;
$\left({x}_{1},{x}_{2}\right)\in{\mathbb{R}}^{2}$.

Figure 3: The plots of the real part of the lattice Green’s functions in sub-
figures (a) and (b) are presented in both the transformed and original
coordinates of the triangular lattice, respectively, wherein $k=2$.

###  5.3 Numerical experiments for various instances

Consider _Problem_ $\mathcal{P}_{\mathrm{ext}}$ with
$\partial\Omega={\Gamma}_{1}\cup{\Gamma}_{2}$. Specifically, we have
${\Gamma}_{1}=\left\\{\left(-5,0\right),\left(-4,0\right)\right\\}$ and
${\Gamma}_{2}=\left\\{\left(4,0\right),\left(5,0\right)\right\\}$. In this
problem, two distinct cases are considered:

  * •

In the first case, referred to as the symmetric mode, we assume that the
function $f\left(y\right)\equiv 1$ is constant on both line segments. This
assumption is made for simplicity.

  * •

In the second case, referred to as the skew-symmetric mode, we assume that
$f\left(y\right)\equiv-1$ on ${\Gamma}_{1}$, while $f\left(y\right)\equiv 1$
on ${\Gamma}_{2}$.

The boundary $\partial\Omega$ consists of four points:
${y}_{1}=\left(-5,0\right)$, ${y}_{2}=\left(-4,0\right)$,
${y}_{3}=-{y}_{2}=\left(4,0\right)$ and ${y}_{4}=-{y}_{1}=\left(5,0\right)$.

The vector $\Phi={\left(\varphi_{1},\dots,\varphi_{4}\right)}^{\top}$ is the
unique solution to equation (4.4), where the coefficient matrix $\mathcal{H}$
of the system (4.4) is symmetric. Specifically, in this context, $\mathcal{H}$
is given by:

$\mathcal{H}=\begin{pmatrix}\mathcal{G}\left(0,0\right)&\mathcal{G}\left(1,0\right)&\mathcal{G}\left(9,0\right)&\mathcal{G}\left(10,0\right)\\\
\mathcal{G}\left(1,0\right)&\mathcal{G}\left(0,0\right)&\mathcal{G}\left(8,0\right)&\mathcal{G}\left(9,0\right)\\\
\mathcal{G}\left(9,0\right)&\mathcal{G}\left(8,0\right)&\mathcal{G}\left(0,0\right)&\mathcal{G}\left(1,0\right)\\\
\mathcal{G}\left(10,0\right)&\mathcal{G}\left(9,0\right)&\mathcal{G}\left(1,0\right)&\mathcal{G}\left(0,0\right)\end{pmatrix}\,.$

In order to solve the obtained system of linear equations (4.4) and find the
solution $u$, we have developed GNU Octave code implementing an efficient
method outlined previously for computing lattice Green’s functions. Notably,
these computations were completed within several minutes on a standard
personal computer. For our analyses, we fix the wavenumber at $k=2$, and
truncate the computation of lattice Green’s functions for $N=2271$.

The results of numerical evaluations are presented in Figure 4 and Figure 5.
In Figure 4, each sub-figure depicts a density plot of
$\mathfrak{Re}\left(u\right)$. Sub-figures (a) and (b) correspond to the
symmetric mode, while sub-figures (c) and (d) represent the skew-symmetric
mode in the coordinates of ${\mathbb{Z}}^{2}$ and ${\mathbb{R}}^{2}$,
respectively. Figure 5 displays the density plot of $\left|u\right|$, with
sub-figures (a) and (b) illustrating the symmetric mode, and sub-figures (c)
and (d) depicting the skew-symmetric mode in the respective coordinate
systems.

Several key features of the numerical solutions are apparent. As expected,
both $\mathfrak{Re}\left(u\right)$ and $\left|u\right|$ exhibit symmetry, as
evidenced by these figures. Moreover, we observe the interference phenomenon
of waves.

(a) Symmetric mode: the density plot of $\mathfrak{Re}\left(u\right)$.

(b) Symmetric mode: the density plot of $\mathfrak{Re}\left(u\right)$.

(c) Skew-symmetric mode: the density plot of $\mathfrak{Re}\left(u\right)$.

(d) Skew-symmetric mode: the density plot of $\mathfrak{Re}\left(u\right)$.

Figure 4: The plots shown in (a) and (c) are represented in
${\mathbb{Z}}^{2}$, whereas plots (b) and (d) are graphed on the original
coordinates of the triangular lattice, where $k=2$.

(a) Symmetric mode: the density plot of $\left|u\right|$.

(b) Symmetric mode: the density plot of $\left|u\right|$.

(c) Skew-symmetric mode: the density plot of $\left|u\right|$.

(d) Skew-symmetric mode: the density plot of $\left|u\right|$.

Figure 5: The plots displayed in (a) and (c) are presented in the coordinate
system of ${\mathbb{Z}}^{2}$, whereas plots (b) and (d) are plotted on the
initial coordinates of the triangular lattice, where the value of $k$ is equal
to $2$.

Finally, we would like to emphasize that we rely on the cone condition
requirement to prove the uniqueness result when applying the unique
continuation property. However, this requirement is not always necessary, and
we can effectively use the unique continuation property for certain
configurations. An example of such a configuration is presented in the
following problem:

Consider the problem denoted _Problem_ $\mathcal{P}_{\mathrm{ext}}$, where the
boundary $\partial\Omega$ is defined as the union of two disjoint sets:
${\Gamma}_{1}$ and ${\Gamma}_{2}$. Specifically, we fix
${\Gamma}_{1}=\left\\{\left(-3,1\right),\left(-2,1\right),\left(-1,1\right),\left(0,1\right),\left(1,1\right)\right\\}$
and
${\Gamma}_{2}=\left\\{\left(-2,-1\right),\left(-1,-1\right),\left(0,-1\right),\left(1,-1\right),\left(2,-1\right)\right\\}$.
The boundary $\partial\Omega$ comprises ten discrete points, herein
represented by ${y}_{i}$ ($i=1,2,\ldots,10$), depicted as red rhombi in Figure
6. It is assumed that the function $f\left({y}_{i}\right)\equiv 1$ remains
constant along both parallel line segments. In this case, the coefficient
matrix $\mathcal{H}$ of a linear system of boundary equations (4.4) has
dimensions $10\times 10$. Following the approach used in the previous example,
we solve this system using a standard GNU Octave routine. We set the
wavenumber to $k=2$, and truncate the computation of lattice Green’s functions
for $N=2271$.

Figure 6: Arrangement of ten discrete boundary points on parallel line
segments within a triangular lattice.

In this case, zeros (cf., end of the proof of Theorem 3.4) propagate freely
between parallel boundaries, yet the uniqueness result persists. The results
of numerical evaluations are presented in Figure 7.

(a) The density plot of $\mathfrak{Re}\left(\mathcal{G}\right)$;
$\left({x}_{1},{x}_{2}\right)\in{\mathbb{Z}}^{2}$.

(b) The density plot of $\mathfrak{Re}\left(\mathcal{G}\right)$;
$\left({x}_{1},{x}_{2}\right)\in{\mathbb{R}}^{2}$.

Figure 7: The plots of the real part of the lattice Green’s functions in sub-
figures (a) and (b) are presented in both the transformed and original
coordinates of the triangular lattice, respectively, wherein $k=2$.

##  6 Discussion

In this paper, we extend our investigation of the discrete Helmholtz equation
and its associated exterior problems within a two-dimensional triangular
lattice, focusing on two main questions: the effectiveness of the numerical
methods used to evaluate the Green’s function, including the case $k=2$, and
the necessity of the cone condition. It is important to note that the
computational accuracy of the solution representation formula (4.3) for the
discrete Helmholtz equation principally relies on the computation of lattice
Green’s functions. To illustrate these points numerically, we investigate two
sample problems. The first sample problem involves four discrete boundary
points arranged along a single line, while the second sample problem consists
of ten discrete boundary points arranged on parallel line segments.

Various numerical quadrature techniques were attempted to handle the double
integrals (3.5) for ${k}^{2}+\imath\varepsilon$. In the sense of the existence
of these integrals, instead of (3.5) the following integrals (3.6) are
considered (for details, see Section 3). However, this approach turned out to
be disadvantageous due to the rapid oscillations of the integrands as the
magnitudes of ${x}_{1}$ and ${x}_{2}$ increase significantly from the origin.
For the numerical computation of the lattice Green’s functions, we employed a
technique originated in the article [2]. Our numerical experiments demonstrate
that this approach effectively approximates the lattice Green’s functions,
showing close agreement with physical solutions. Furthermore, this approach
indicates stability properties when computations are initiated from a large
Manhattan distance (refer to Table 1).

In the second sample problem, despite not satisfying the cone condition, the
uniqueness result is still attained. Numerical experiments revealed that as
more boundary points are aligned on parallel lines (typically, the boundary
points follow the configuration illustrated in Figure 6), the determinant of
the coefficient matrix $\mathcal{H}$ of the linear system of boundary
equations (4.4) approaches zero. In the given scenario, it is observed that
$\left|\det\left(\mathcal{H}\right)\right|\approx 5.2309e{-}6$. Concurrently,
the condition number of the matrix $\mathcal{H}$, with respect to the matrix
norm induced by the (vector) Euclidean norm, is given by
${\kappa}_{2}\left(\mathcal{H}\right)={\left\|{\mathcal{H}}^{-1}\right\|}_{2}{\left\|\mathcal{H}\right\|}_{2}\approx
1.5331e{+}1$. The fact that the determinant of the coefficient matrix
$\mathcal{H}$ is of the order of ${10}^{-6}$ suggests that the matrix is
nearly singular. However, with a condition number less than $16$, it indicates
that despite the near singularity of the matrix $\mathcal{H}$, the system
(4.4) remains relatively stable and well-behaved.

Another interesting question pertains to the structure of the space of
radiating solutions. We focus our attention on the space
${\ell}_{\mathrm{R}}^{\infty}\left(\Omega\right)$, a Banach space comprising
all bounded sequences on $\Omega\subset{\mathbb{Z}}^{2}$ that satisfy the
prescribed radiation condition (3.9). Depending on the specific objectives of
the investigation, various spaces on lattices can be considered, as detailed
in [1, 26], which delves into the spectral properties of discrete Schrödinger
operators within lattice frameworks.

## Acknowledgement

The authors wish to express their gratitude to Prof. Dr. Jemal Rogava for his
helpful remarks concerning the numerical computation section of this article.

This work was supported by the Shota Rustaveli National Science Foundation of
Georgia (SRNSFG) [grant number: FR-21-301, project title: “Metamaterials with
Cracks and Wave Diffraction Problems”].

##  Appendix A Describing of assembling of sparse matrices

We introduce a set of functions for constructing sparse matrices
${\alpha}_{n}\left(k\right)$, ${\beta}_{n}\left(k\right)$, and
${\gamma}_{n}\left(k\right)$, employing a language-agnostic approach. This
procedure is delineated into two distinct algorithms: Algorithm 1 outlines the
procedure for generating these sparse matrices when $n$ is even, while
Algorithm 2 delineates the process for odd values of $n$. Specifically, the
functions Alpha2$\left(p,k\right)$, Beta2$\left(p,k\right)$, and
Gamma2$\left(p,k\right)$ in Algorithm 1 correspond to the computation of
sparse matrices ${\alpha}_{2p}\left(k\right)$, ${\beta}_{2p}\left(k\right)$,
and ${\gamma}_{2p}\left(k\right)$, respectively. Furthermore, the functions
Alpha1$\left(p,k\right)$, Beta1$\left(p,k\right)$, and
Gamma1$\left(p,k\right)$ in Algorithm 2 clarify the computations of matrices
${\alpha}_{2p+1}\left(k\right)$, ${\beta}_{2p+1}\left(k\right)$, and
${\gamma}_{2p+1}\left(k\right)$, respectively.

In Table 2, we present the size and the number of nonzero-valued elements of
these sparse matrices for each case.

Matrix | Size | Number of nonzero elements
---|---|---
${\alpha}_{2p}\left(k\right)$ | $\left(p+1\right)\times p$ | $2p$
${\beta}_{2p}\left(k\right)$ | $\left(p+1\right)\times\left(p+1\right)$ | $2p+1$
${\gamma}_{2p}\left(k\right)$ | $\left(p+1\right)\times\left(p+1\right)$ | $3p+1$
${\alpha}_{2p+1}\left(k\right)$ | $\left(p+1\right)\times\left(p+1\right)$ | $2p+1$
${\beta}_{2p+1}\left(k\right)$ | $\left(p+1\right)\times\left(p+2\right)$ | $2p+2$
${\gamma}_{2p+1}\left(k\right)$ | $\left(p+1\right)\times\left(p+1\right)$ | $3p+1$
Table 2: The dimensions of sparse matrices and their corresponding nonzero
elements in computing lattice Green’s functions. Algorithm 1 Assembly of
sparse matrices ${\alpha}_{2p}\left(k\right)$, ${\beta}_{2p}\left(k\right)$,
and ${\gamma}_{2p}\left(k\right)$.

1:$p$ (positive integer), $k$ $\triangleright$ Input parameters: $p\geq 1$ and
$k\in\left(0,2\sqrt{2}\right)$.

2:

3:function Alpha2($p,k$)

4: Initialize ${\alpha}_{2}$ as a $\left(p+1\right)\times p$ matrix with all
entries set to zero.

5: for $i\mathrel{\coloneqq}1$ to $p$ do

6: ${\alpha}_{2}\left[i,i\right]\mathrel{\coloneqq}1$ $\triangleright$ Set
diagonal entries to $1$.

7: if $i\geq 2$ then

8: ${\alpha}_{2}\left[i,i-1\right]\mathrel{\coloneqq}1$ $\triangleright$ Set
sub-diagonal entries to $1$.

9: end if

10: end for

11: ${\alpha}_{2}\left[p+1,p\right]\mathrel{\coloneqq}2$ $\triangleright$ Set
the $\left(p+1,p\right)$ entry to $2$.

12: return ${\alpha}_{2}$ $\triangleright$ Return the computed matrix
${\alpha}_{2}$.

13:end function

14:

15:function Beta2($p,k$)

16: Initialize ${\beta}_{2}$ as a $\left(p+1\right)\times\left(p+1\right)$
matrix with all entries set to zero.

17: ${\beta}_{2}\left[p+1,p+1\right]\mathrel{\coloneqq}2$ $\triangleright$ Set
the value of the bottom-right corner entry to $2$.

18:
${\beta}_{2}\left[1,2\right]\mathrel{\coloneqq}{\beta}_{2}\left[p+1,p+1\right]$
$\triangleright$ Copy the value from the bottom-right corner entry to entry
$\left(1,2\right)$.

19: for $i\mathrel{\coloneqq}1$ to $p$ do

20: ${\beta}_{2}\left[i,i\right]\mathrel{\coloneqq}1$ $\triangleright$ Set
diagonal entries to $1$.

21: if $i\geq 2$ then $\triangleright$ Check if $i$ is greater than or equal
to $2$.

22: ${\beta}_{2}\left[i,i+1\right]\mathrel{\coloneqq}1$ $\triangleright$ Set
super-diagonal entries to $1$.

23: end if

24: end for

25: return ${\beta}_{2}$ $\triangleright$ Return the computed matrix
${\beta}_{2}$.

26:end function

27:

28:function Gamma2($p,k$)

29: Initialize ${\gamma}_{2}$ as a $\left(p+1\right)\times\left(p+1\right)$
matrix with all entries set to zero.

30: ${\gamma}_{2}\left[1,2\right]\mathrel{\coloneqq}-2$ $\triangleright$ Set
the $\left(1,2\right)$ entry to $-2$.

31:
${\gamma}_{2}\left[p+1,p\right]\mathrel{\coloneqq}{\gamma}_{2}\left[1,2\right]$
$\triangleright$ Copy the value from entry $\left(1,2\right)$ to entry
$\left(p+1,p\right)$.

32: for $i\mathrel{\coloneqq}1$ to $p+1$ do

33: ${\gamma}_{2}\left[i,i\right]\mathrel{\coloneqq}6-k^{2}$ $\triangleright$
Set diagonal entries to $6-k^{2}$.

34: if $i\geq 2$ and $i\leq p$ then $\triangleright$ Set the off-diagonal
elements if applicable.

35: ${\gamma}_{2}\left[i,i+1\right]\mathrel{\coloneqq}-1$ $\triangleright$ Set
super-diagonal entries to $-1$.

36: ${\gamma}_{2}\left[i,i-1\right]\mathrel{\coloneqq}-1$ $\triangleright$ Set
sub-diagonal entries to $-1$.

37: end if

38: end for

39: return ${\gamma}_{2}$ $\triangleright$ Return the computed matrix
${\gamma}_{2}$.

40:end function

Algorithm 2 Assembly of sparse matrices ${\alpha}_{2p+1}\left(k\right)$,
${\beta}_{2p+1}\left(k\right)$, and ${\gamma}_{2p+1}\left(k\right)$.

1:$p$ (non-negative integer), $k$ $\triangleright$ Input parameters: $p\geq 0$
and $k\in\left(0,2\sqrt{2}\right)$.

2:

3:function Alpha1($p,k$)

4: Initialize ${\alpha}_{1}$ as a $\left(p+1\right)\times\left(p+1\right)$
matrix with all entries set to zero.

5: for $i\mathrel{\coloneqq}1$ to $p+1$ do

6: ${\alpha}_{1}\left[i,i\right]\mathrel{\coloneqq}1$ $\triangleright$ Set
diagonal entries to $1$.

7: if $i\geq 2$ then $\triangleright$ Set the sub-diagonal elements if
applicable.

8: ${\alpha}_{1}\left[i,i-1\right]\mathrel{\coloneqq}1$ $\triangleright$ Set
sub-diagonal entries to $1$.

9: end if

10: end for

11: return ${\alpha}_{1}$ $\triangleright$ Return the computed matrix
${\alpha}_{1}$.

12:end function

13:

14:function Beta1($p,k$)

15: Initialize ${\beta}_{1}$ as a $\left(p+1\right)\times\left(p+2\right)$
matrix with all entries set to zero.

16: ${\beta}_{1}\left[1,2\right]\mathrel{\coloneqq}2$ $\triangleright$ Set the
$\left(1,2\right)$ entry to $2$.

17: for $i\mathrel{\coloneqq}1$ to $p+1$ do

18: ${\beta}_{1}\left[i,i\right]\mathrel{\coloneqq}1$ $\triangleright$ Set
diagonal entries to $1$.

19: if $i\geq 2$ then $\triangleright$ Set the super-diagonal elements if
applicable.

20: ${\beta}_{1}\left[i,i+1\right]\mathrel{\coloneqq}1$ $\triangleright$ Set
super-diagonal entries to $1$.

21: end if

22: end for

23: return ${\beta}_{1}$ $\triangleright$ Return the computed matrix
${\beta}_{1}$.

24:end function

25:

26:function Gamma1($p,k$)

27: Initialize ${\gamma}_{1}$ as a $\left(p+1\right)\times\left(p+1\right)$
matrix with all entries set to zero.

28: if $p=0$ then $\triangleright$ Consider the case when ${\gamma}_{1}$ is a
$1\times 1$ matrix.

29: ${\gamma}_{1}\left[p+1,p+1\right]\mathrel{\coloneqq}4-k^{2}$
$\triangleright$ Set the only entry at position $\left(1,1\right)$ to
$4-k^{2}$.

30: else$\triangleright$ For non-zero values of $p$.

31: ${\gamma}_{1}\left[p+1,p+1\right]\mathrel{\coloneqq}5-k^{2}$
$\triangleright$ Set the $\left(p+1,p+1\right)$ entry to $5-k^{2}$.

32: ${\gamma}_{1}\left[1,2\right]\mathrel{\coloneqq}-2$ $\triangleright$ Set
the $\left(1,2\right)$ entry to $-2$.

33: ${\gamma}_{1}\left[p+1,p\right]\mathrel{\coloneqq}-1$ $\triangleright$ Set
the $\left(p+1,p\right)$ entry to $-1$.

34: for $i\mathrel{\coloneqq}1$ to $p$ do

35: ${\gamma}_{1}\left[i,i\right]\mathrel{\coloneqq}6-k^{2}$ $\triangleright$
Set diagonal entries to $6-k^{2}$.

36: if $i\geq 2$ then $\triangleright$ Set the off-diagonal elements if
applicable.

37: ${\gamma}_{1}\left[i,i+1\right]\mathrel{\coloneqq}-1$ $\triangleright$ Set
super-diagonal entries to $-1$.

38: ${\gamma}_{1}\left[i,i-1\right]\mathrel{\coloneqq}-1$ $\triangleright$ Set
sub-diagonal entries to $-1$.

39: end if

40: end for

41: end if

42: return ${\gamma}_{1}$ $\triangleright$ Return the computed matrix
${\gamma}_{1}$.

43:end function

## References

  * Ando et al. [2016] K. Ando, H. Isozaki, and H. Morioka.  Spectral properties of Schrödinger operators on perturbed lattices.  _Ann. Henri Poincaré_ , 17(8):2103–2171, 2016\.  ISSN 1424-0637.  doi:10.1007/s00023-015-0430-0.  URL https://doi.org/10.1007/s00023-015-0430-0.
  * Berciu and Cook [2010] M. Berciu and A. M. Cook.  Efficient computation of lattice green’s functions for models with nearest-neighbour hopping.  _EPL (Europhysics Letters)_ , 92(4):40003, Dec 2010.  ISSN 1286-4854.  doi:10.1209/0295-5075/92/40003.  URL http://dx.doi.org/10.1209/0295-5075/92/40003.
  * Bloch [1929] F. Bloch.  Über die Quantenmechanik der Elektronen in Kristallgittern.  _Zeitschrift für Physik_ , 52(7–8):555–600, July 1929.  ISSN 1434-601X.  doi:10.1007/bf01339455.  URL http://dx.doi.org/10.1007/BF01339455.
  * Born and Huang [1998] M. Born and K. Huang.  _Dynamical theory of crystal lattices_.  Oxford Classic Texts in the Physical Sciences. The Clarendon Press, Oxford University Press, New York, 1998.  ISBN 0-19-850369-5.  Reprint of the 1954 original.
  * Brillouin [1953] L. Brillouin.  _Wave propagation in periodic structures. Electric filters and crystal lattices_.  Dover Publications, Inc., New York, N.Y., 2nd edition, 1953.
  * Burke [1966] J.G. Burke.  _Origins of the Science of Crystals_.  University of California Press, 1966.  ISBN 9780520001985.  URL https://books.google.ge/books?id=qvxPbZtJu8QC.
  * Caloz and Itoh [2005] C. Caloz and T. Itoh.  _Electromagnetic Metamaterials: Transmission Line Theory and Microwave Applications: The Engineering Approach_.  Wiley, November 2005.  ISBN 9780471754329.  doi:10.1002/0471754323.  URL http://dx.doi.org/10.1002/0471754323.
  * Castro Neto et al. [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim.  The electronic properties of graphene.  _Rev. Mod. Phys._ , 81:109–162, Jan 2009.  doi:10.1103/RevModPhys.81.109.  URL https://link.aps.org/doi/10.1103/RevModPhys.81.109.
  * Colton and Kress [2019] D. Colton and R. Kress.  _Inverse acoustic and electromagnetic scattering theory_ , volume 93 of _Applied Mathematical Sciences_.  Springer, Cham, 2019.  ISBN 978-3-030-30350-1; 978-3-030-30351-8.  doi:10.1007/978-3-030-30351-8.  URL https://doi.org/10.1007/978-3-030-30351-8.  Fourth edition.
  * Economou [2006] E. N. Economou.  _Green’s functions in quantum physics_ , volume 7 of _Springer Series in Solid-State Sciences_.  Springer-Verlag, Berlin, third edition, 2006.  ISBN 978-3-540-28838-1; 3-540-28838-4.  doi:10.1007/3-540-28841-4.  URL https://doi.org/10.1007/3-540-28841-4.
  * Harrison [1989] W.A. Harrison.  _Electronic Structure and the Properties of Solids: The Physics of the Chemical Bond_.  Dover Books on Physics. Dover Publications, 1989.  ISBN 9780486660219.  URL https://books.google.ge/books?id=UPBIAwAAQBAJ.
  * Horiguchi [1972] T. Horiguchi.  Lattice Green’s functions for the triangular and honeycomb lattices.  _J. Mathematical Phys._ , 13:1411–1419, 1972.  ISSN 0022-2488.  doi:10.1063/1.1666155.  URL https://doi.org/10.1063/1.1666155.
  * Isozaki and Morioka [2014] H. Isozaki and H. Morioka.  A Rellich type theorem for discrete Schrödinger operators.  _Inverse Probl. Imaging_ , 8(2):475–489, 2014\.  ISSN 1930-8337.  doi:10.3934/ipi.2014.8.475.  URL https://doi.org/10.3934/ipi.2014.8.475.
  * Kapanadze [2018] D. Kapanadze.  Exterior diffraction problems for two-dimensional square lattice.  _Z. Angew. Math. Phys._ , 69(5):Paper No. 123, 17, 2018.  ISSN 0044-2275.  doi:10.1007/s00033-018-1019-5.  URL https://doi.org/10.1007/s00033-018-1019-5.
  * Kapanadze [2021] D. Kapanadze.  The far-field behaviour of Green’s function for a triangular lattice and radiation conditions.  _Math. Methods Appl. Sci._ , 44(17):12746–12759, 2021.  ISSN 0170-4214.  doi:10.1002/mma.7575.  URL https://doi.org/10.1002/mma.7575.
  * Kapanadze and Pesetskaya [2023a] D. Kapanadze and E. Pesetskaya.  Exterior diffraction problems for a triangular lattice.  _Math. Mech. Solids_ , 28(12):2596–2609, 2023a.  ISSN 1081-2865.  doi:10.1177/10812865231171111.  URL https://doi.org/10.1177/10812865231171111.
  * Kapanadze and Pesetskaya [2023b] D. Kapanadze and E. Pesetskaya.  Half-plane diffraction problems on a triangular lattice.  _J. Engrg. Math._ , 138:Paper No. 5, 2023b.  ISSN 0022-0833.  doi:10.1007/s10665-022-10252-5.  URL https://doi.org/10.1007/s10665-022-10252-5.
  * Kapanadze and Pesetskaya [2024] D. Kapanadze and E. Pesetskaya.  Wave propagation on hexagonal lattices.  _Georgian Math. J._ , 2024.  doi:10.1515/gmj-2024-2035.  URL https://doi.org/10.1515/gmj-2024-2035.  Ahead of Print/Just Accepted.
  * Katsura and Inawashiro [1971] S. Katsura and S. Inawashiro.  Lattice Green’s functions for the rectangular and the square lattices at arbitrary points.  _J. Mathematical Phys._ , 12:1622–1630, 1971.  ISSN 0022-2488.  doi:10.1063/1.1665785.  URL https://doi.org/10.1063/1.1665785.
  * Kirsch and Hettlich [2015] A. Kirsch and F. Hettlich.  _The mathematical theory of time-harmonic Maxwell’s equations_ , volume 190 of _Applied Mathematical Sciences_.  Springer, Cham, 2015.  ISBN 978-3-319-11085-1; 978-3-319-11086-8.  doi:10.1007/978-3-319-11086-8.  URL https://doi.org/10.1007/978-3-319-11086-8.  Expansion-, integral-, and variational methods.
  * Lifshitz and Kosevich [1966] I. M. Lifshitz and A. M. Kosevich.  The dynamics of a crystal lattice with defects.  _Reports on Progress in Physics_ , 29(1):217–254, January 1966.  ISSN 0034-4885.  doi:10.1088/0034-4885/29/1/305.  URL http://dx.doi.org/10.1088/0034-4885/29/1/305.
  * Maradudin et al. [1963] A. A. Maradudin, E. W. Montroll, and G. H. Weiss.  _Theory of lattice dynamics in the harmonic approximation_.  Solid State Physics, Supplement 3. Academic Press, New York-London, 1963\.
  * Martin [2006] P. A. Martin.  Discrete scattering theory: Green’s function for a square lattice.  _Wave Motion_ , 43(7):619–629, 2006.  ISSN 0165-2125.  doi:10.1016/j.wavemoti.2006.05.006.  URL https://doi.org/10.1016/j.wavemoti.2006.05.006.
  * Morita [1971] T. Morita.  Useful procedure for computing the lattice Green’s function—square, tetragonal, and bcc lattices.  _J. Mathematical Phys._ , 12:1744–1747, 1971.  ISSN 0022-2488.  doi:10.1063/1.1665800.  URL https://doi.org/10.1063/1.1665800.
  * Novoselov [2011] K. S. Novoselov.  Nobel lecture: Graphene: Materials in the flatland.  _Rev. Mod. Phys._ , 83:837–849, Aug 2011.  doi:10.1103/RevModPhys.83.837.  URL https://link.aps.org/doi/10.1103/RevModPhys.83.837.
  * Parra and Richard [2018] D. Parra and S. Richard.  Spectral and scattering theory for Schrödinger operators on perturbed topological crystals.  _Rev. Math. Phys._ , 30(4):1850009, 39, 2018.  ISSN 0129-055X.  doi:10.1142/S0129055X18500095.  URL https://doi.org/10.1142/S0129055X18500095.
  * Shaban and Vainberg [2001] W. Shaban and B. Vainberg.  Radiation conditions for the difference Schrödinger operators.  _Appl. Anal._ , 80(3-4):525–556, 2001.  ISSN 0003-6811.  doi:10.1080/00036810108841007.  URL https://doi.org/10.1080/00036810108841007.
  * Sharma [2016] B. L. Sharma.  Edge diffraction on triangular and hexagonal lattices: existence, uniqueness, and finite section.  _Wave Motion_ , 65:55–78, 2016.  ISSN 0165-2125.  doi:10.1016/j.wavemoti.2016.04.005.  URL https://doi.org/10.1016/j.wavemoti.2016.04.005.
  * Slater and Koster [1954] J. C. Slater and G. F. Koster.  Simplified LCAO method for the periodic potential problem.  _Physical Review_ , 94(6):1498–1524, June 1954\.  ISSN 0031-899X.  doi:10.1103/physrev.94.1498.  URL http://dx.doi.org/10.1103/PhysRev.94.1498.
  * Thomson et al. [1992] R. Thomson, S. J. Zhou, A. E. Carlsson, and V. K. Tewary.  Lattice imperfections studied by use of lattice green’s functions.  _Phys. Rev. B_ , 46:10613–10622, Nov 1992.  doi:10.1103/PhysRevB.46.10613.  URL https://link.aps.org/doi/10.1103/PhysRevB.46.10613.
  * Zemła [1995] A. Zemła.  On the fundamental solutions for the difference Helmholtz operator.  _SIAM J. Numer. Anal._ , 32(2):560–570, 1995\.  ISSN 0036-1429.  doi:10.1137/0732024.  URL https://doi.org/10.1137/0732024.