Half-plane diffraction problems on a triangular lattice

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. For example, close-packed planes occur frequently in some kinds of crystals in the form of triangular lattice [2, 3]. Besides, one can consider two-dimensional passive propagation media, a host microstrip line network periodically loaded with series capacitors and shunt inductors for signal processing and filtering (as shown in Figure 1). This type of inductor-capacitor lattice is referred to a negative-refractive-index transmission-line (NRI-TL) metamaterial [6] or, simply, left-handed 2D metamaterial.

Motivated by applications of recent interest related to analog circuits, crystalline materials, and metamaterials, we investigate thin-slit diffraction problems for triangular lattices. The peculiar structure allows us to consider the problems on the semi-infinite triangular lattice. For this lattice, we study Dirichlet problems for the two-dimensional discrete Helmholtz equation in a half-plane (mathematically formulated in the next section). When one is interested in the analysis of regular processes in which waves corresponding to the microstructural scales can be neglected, then the continuum limit of corresponding equations can be investigated. In this case, we arrive at the famous problem in applied mathematics – the half-plane diffraction problem. It is well known that the classical continuum model of wave diffraction can be considered only as the slowly-varying approximation of a discrete or structured material. Nowadays the understanding of nanostructure and microstructure phenomena in modern materials and composites in critical conditions is of crucial importance. It turns out that a discrete structure provides an effective way to describe microstructural processes and influences within, cf., e.g., [5, 6, 13, 17]. Therefore, there is an obvious necessity to avoid continuum limits and instead directly analyze the discrete Helmholtz diffraction problems. Although the similar problem for square lattice has been studied in [1], see also [15], its extension to a triangular lattice model is not direct.

Our main interest lies in the investigation of the problems with wave number $k$ within the pass-band which is an arbitrary non-zero complex number in general. In this paper, we use the well approved approach to provide new results for a real wave number $k\in(0,3)\backslash\{2\sqrt{2}\}$ without passing to a complex wave number. In order to have the unique solvability results and an exact representation formula for the solution, we use the notion of the radiating solution [14] and the method of images to construct the Dirichlet Green’s function for the half-plane. Notice that the cases with the wave numbers within the stop-band is mathematically more simple, they do not need the radiation conditions and is not considered here. Finally, we propose a method for numerical calculation. The efficiency of our approach is demonstrated in an example related to the propagation of wave fronts in metamaterials through two small openings.

Following the customary notation in mathematics, let $\mathbb{Z}$, $\mathbb{Z}^{+}$, $\mathbb{Z}^{-}$, $\mathbb{N}$, $\mathbb{R}$, and $\mathbb{C}$ denote the set of integers, positive integers, negative integers, non-negative integers, real numbers and complex numbers, respectively. We denote by $e_{1}=(1,0)$, $e_{2}=(0,1)$ the standard base of $\mathbb{Z}^{2}$ ($=\mathbb{Z}\times\mathbb{Z}$).

Consider a two-dimensional infinite triangular lattice $\mathfrak{T}$ defined as a periodic simple graph $\{\mathcal{V},\mathcal{E}\}$, where

is a vertex set, and $\mathcal{E}$ is an edge set, whose endpoints $(a,b)\in\mathcal{V}\times\mathcal{V}$ are adjacent points, i.e., $\mid a-b\mid=1$, cf. Figure 1. The time-harmonic discrete waves in $\mathfrak{T}$ can be described by solutions of the discrete Helmholtz equation

Here, $\Delta_{d}$ denotes the discrete (a 7-point) Laplacian defined as follows

where $g$ is a function with bounded support, i.e., $g\in C_{0}(\mathbb{Z}^{2})$.

Let $\Omega=\mathring{\Omega}\cup\Gamma$ be a discrete half-plane in $\mathbb{Z}^{2}$, where $\mathring{\Omega}:=\mathbb{Z}\times\mathbb{\mathbb{Z}^{+}}$ and $\Gamma:=\partial\Omega=\{(x_{1},0):x_{1}\in\mathbb{Z}\}$. We consider the Dirichlet poblem for the discrete Helmholtz equation on the semi-infinite lattice:

Here, $f$ is a given function supported on a finite subset of $\Gamma$.

Thus, we are interested in studying the problem of the existence and uniqueness of a function $u:\Omega\to\mathbb{C}$ such that $u$ satisfies the discrete Helmoltz equation (2a) with $k\in(0,3)\backslash\{2\sqrt{2}\}$ and the boundary condition (2b). From now on we will refer to this problem as Problem $\mathcal{P}_{\mathrm{D}}$.

Denote by $\mathcal{G}(x;y)$ the Green’s function for the discrete Helmholtz equation (1) centered at $y$ and evaluated at $x$. Then the function $\mathcal{G}(x;y)$ satisfies

where $\delta_{x,y}$ is the Kronecker delta. For brevity, we use the notation $\mathcal{G}(x)$ for $\mathcal{G}(x;0)$. Notice that $\mathcal{G}(x;y)=\mathcal{G}(x-y)$.

Using the discrete Fourier transform and the inverse Fourier transform we get

where

The lattice Green’s function $\mathcal{G}$ is quite well known when $k^{2}\in\mathbb{C}\backslash[0,9]$ (cf., e.g., [16]). Notice that if $k^{2}\in\mathbb{C}\backslash[0,9]$ then $\sigma\neq 0$ and, consequently, $\mathcal{G}$ in (4) is well defined. In this case $\mathcal{G}(x)$ decays exponentially when $\mid x\mid\to\infty$. Moreover, we have

for all $x=(x_{1},x_{2})\in\mathbb{Z}^{2}$.

Let us show $\mathcal{G}(x_{1},x_{2})=\mathcal{G}(x_{1}+x_{2},-x_{2})$ while other identities are trivial. According to (4), we have

with $\eta_{1}=\xi_{1},\,\eta_{2}=\xi_{1}-\xi_{2}$. The following equality can be easily obtained by changing the variable $\eta_{2}$ to $\eta_{2}-2\pi$

The factor $e^{\iota x_{2}2\pi}$ is equal to $1$ since $x_{2}\in\mathbb{Z}$. Similar arguments give us

Taking into account the last two equalities, we get

In the case $k^{2}\in(0,9)\backslash\{8\}$ the expression (4) is understood as follows: we replace $k^{2}$ by $k^{2}+\iota\varepsilon$ with $0<\varepsilon\ll 1$ and let $\varepsilon\to 0$ at the end of the calculation, cf. [14]. Thus, we define the lattice Green’s function for $k\in(0,3)\backslash\{2\sqrt{2}\}$ as a pointwise limit of

as $k^{2}+\iota\varepsilon\to k^{2}+\iota 0$ and denote it again by $\mathcal{G}(x)$, i.e., $\mathcal{G}(x)=(R_{\lambda+\iota 0}\delta_{x,0})(x)$. Clearly, $\mathcal{G}(x)$ is a solution to equation (3) and satisfies equalities (6).

In order to simplify further arguments let us introduce following vectors:

For any point $x\in\mathbb{Z}^{2}$ we define the 6-neighbourhood $F^{0}_{x}$ as the set of points $\{x+e_{j}:j=1,...,6\}$ and the neighbourhood $F_{x}$ as $F^{0}_{x}\bigcup\{x\}$. We say that $R\subset\mathbb{Z}^{2}$ is a region if there exist disjoint nonempty subsets $\mathring{R}$ and $\partial R$ of $R$ such that

• (a)

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

• (b)

  if $x\in\mathring{R}$ then $F_{x}\subset R$,

• (c)

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

Clearly, the subsets $\mathring{R}$ and $\partial R$ are not defined uniquely by $R$, but henceforth, it will be always assumed that $\mathring{R}$ and $\partial R$ are also given and fixed for a given region $R$ in $\mathbb{Z}^{2}$. We also say that $x$ is an interior (boundary) point of $R$ if $x\in\mathring{R}$ ($x\in\partial R$).

Denote by $(\partial R)_{j}$, $j=1,...,6$, a set of all boundary points $y\in\partial R$ such that $y-e_{j}\in\mathring{R}$ and call it the sides of the boundary $R$. Note that a boundary point $x$ can simultaneously belong to all six sides of $R$. Thus, these sides may overlap each other. Clearly, $\partial R$ is the union of its six sides, $\partial R=\cup_{j=1}^{6}(\partial R)_{j}$.

As it was already mentioned above, $y\in\partial R$ may be a point of intersection of several sides of $\partial R$. However, in our arguments presented below it will be always clear which side is needed to be considered. Under this condition, we define the discrete derivative in the outward normal direction $e_{j}$, $j=1,\dots,6$,

Let us introduce the following set $H_{0}=\{(0,0)\}$ and then define $H_{N}$, $N\in\mathbb{N}$, with the help of recurrence formula

with $\mathring{H}_{N}:=H_{N-1}$ and $(\partial H)_{N}:=H_{N}\backslash\mathring{H}_{N}$.

Let $R$ be a finite region. Representing $R=\bigcup_{x\in\mathring{R}}{F_{x}}$, we can easily derive a discrete analogue of Green’s second identity

Now let us give a definition of a radiating solution on the discrete half-plane. First, we consider the case $k^{2}\in(0,8)$. We say that $u:\Omega\to\mathbb{C}$ satisfies the radiation condition at infinity if

with the remaining term decaying uniformly in all directions $x/\mid x\mid$, where $x$ is characterized as $x_{1}=\mid x\mid\cos\alpha$, $x_{2}=\mid x\mid\sin\alpha$, $0\leq\alpha\leq\pi$. Here, $\xi^{*}_{j}(\alpha,k)$ is the $j$th coordinate of the point $\xi^{*}(\alpha,k)$.

For $k^{2}\in(8,9)$ we say that $u:\Omega\to\mathbb{C}$ satisfies the radiation condition at infinity if it can be represented as a sum of functions $u_{1}$ and $u_{2}$ such that each of them satisfies

with the remaining term decaying uniformly in all directions $x/\mid x\mid$, where $x$ is characterized as $x_{1}=\mid x\mid\cos\alpha$, $x_{2}=\mid x\mid\sin\alpha$, $0\leq\alpha\leq\pi$. Here, $\xi^{*,i}_{j}(\alpha,k)$ is the $j$th coordinate of the point $\xi^{*,i}(\alpha,k)$, $i=1,2$.

It is shown in [14] that $\mathcal{G}(x)$ satisfies the radiation conditions introduced above. Now we are ready to introduce a notion of a radiating solution to the discrete Helmholtz equation (1).

We represent the Dirichlet Green’s function for the half-plane as follows:

with $x=(x_{1},x_{2})\in\Omega$ and $\hat{x}=(x_{1}+x_{2},-x_{2})$. Notice that $\mathcal{G}^{+}(x;y)$ can be represented equivalently as

Indeed, using the property (6), we have

Since $\widehat{(y_{1},0)}=(y_{1},0)$, we immediately get $\mathcal{G}^{+}(x;y)=0$ for $y\in\Gamma$. Moreover,

for all $x\in\mathring{\Omega}$ and $y\in\Omega$. From (10) we see that $\mathcal{G}^{+}(x;y)$ satisfies the radiation condition. Indeed, for any fixed $y\in\Gamma$ the angle $\alpha$ from the radiation conditions (8) and (9) is the same for $\mathcal{G}(x;y)=\mathcal{G}(x-y)$ and $\mathcal{G}(x;\hat{y})=\mathcal{G}(x-\hat{y})$ when $\mid x\mid\to\infty$.

• Proof.  Denote by $\tilde{u}$ the extension of $u$ to $\mathbb{Z}^{2}$ by zeros, i.e, $\tilde{u}(x)=u(x)$ if $x\in\Omega$ and $\tilde{u}(x)=0$ if $x\in\mathbb{Z}\times\mathbb{Z}^{-}$. Then, for any finite region $H_{N}$, $N\in\mathbb{N}$, we apply Green’s second identity (7) where we take $v(y)=\mathcal{G}^{+}(x;y)$. Here, note that $\mathcal{G}^{+}$ is a fundamental solution for $y\in\mathring{\Omega}$, and $\tilde{u}(y)=0$ for $y\in\mathbb{Z}\times\mathbb{Z}^{-}$. Thus, $u(y)(\Delta_{d}+k^{2})v(y)$ disappears in the following identity

  $$u(y)\Delta_{d}v(y)-v(y)\Delta_{d}u(y)=u(y)(\Delta_{d}+k^{2})v(y)-v(y)(\Delta_{d}+k^{2})u(y)$$

  as far as $y\in\mathbb{Z}\times\mathbb{Z}^{-}$. For $x\in\mathring{H}_{N}$, we get

  	$\displaystyle\tilde{u}(x)$	$\displaystyle=\sum_{y\in\mathring{H}_{N}}\tilde{u}(y)\delta_{x,y}$
  		$\displaystyle=\sum_{y\in\partial H_{N}}\big{(}\tilde{u}(y)T\mathcal{G}^{+}(x;y)-\mathcal{G}^{+}(x;y)T\tilde{u}(y)\big{)}+\sum_{y\in\mathring{H}_{N}}(\mathcal{G}^{+}(x;y)(\Delta_{d}+k^{2})\tilde{u}(y)).$

  Passing to the limit $N\to\infty$, we use exactly the same arguments as in the proof of Theorem 5.2 from [14] and, consequently, we get

  $$\tilde{u}(x)=\sum_{y\in\mathbb{Z}^{2}}\mathcal{G}^{+}(x;y)(\Delta_{d}+k^{2})\tilde{u}(y).$$

  For the function $u$, we obtain

  $$u(x)=\sum_{y\in\Gamma}\big{(}u(y)T\mathcal{G}^{+}(x;y)-\mathcal{G}^{+}(x;y)Tu(y)\big{)}+\sum_{y\in\mathring{\Omega}}\mathcal{G}^{+}(x;y)(\Delta_{d}+k^{2})u(y).$$

  Taking into the account that $\mathcal{G}^{+}(x;y)=0$ when $y\in\Gamma$, for a solution $u$ to the discrete Helmholtz equation we get the following quality

  $$u(x)=\sum_{y\in\Gamma}u(y)T\mathcal{G}^{+}(x;y),\quad x\in\mathring{\Omega}.$$

  $\Box$

Now we are ready to prove the unique solvability result for the discrete Helmholtz equation on the semi-infinite triangular lattice.

• Proof.  To prove the uniqueness result, it is sufficient to show that the corresponding homogeneous problem has only the trivial solution. Let $u$ be a radiating solution to the homogeneous problem $\mathcal{P}_{D}$. Then Theorem 4.2 immediately implies $u(x)=0$ for all $x\in\mathring{\Omega}$.

  To show the existence results, let us first check that $u(x)$ satisfies the boundary condition (2b). Since $\mathcal{G}^{+}(x;y)=0$ for all $x\in\Gamma$ and any $y$ we get

  $$u(x_{1},0)=\sum_{y_{1}\in\mathbb{R}}\delta_{x_{1},y_{1}}\delta_{0,0}f(y)=f(x_{1}).$$

  Now let us check that $u(x)$ satisfies the Helmholtz equation (2a). For $x_{2}>1$ all terms $(\Delta_{d}+k^{2})\delta_{x,y}$, $(\Delta_{d}+k^{2})\mathcal{G}^{+}(x;y+e_{2})=\delta_{x,y+e_{2}}$ and $(\Delta_{d}+k^{2})\mathcal{G}^{+}(x;y+e_{2}-e_{1})=\delta_{x,y+e_{2}-e_{1}}$ are equal to zero. Consequently, $(\Delta_{d}+k^{2})u(x)=0$ for points $x=(x_{1},x_{2})$ with $x_{2}>1$. It remains to consider the case $x_{2}=1$. By the direct calculation we have

  	$\displaystyle(\Delta_{d}+k^{2})\delta_{x,y}f(y)$	$\displaystyle=$	$\displaystyle f(x_{1})+f(x_{1}+1),$
  	$\displaystyle(\Delta_{d}+k^{2})\mathcal{G}^{+}(x;y+e_{2})f(y)$	$\displaystyle=$	$\displaystyle f(x_{1}),$
  	$\displaystyle(\Delta_{d}+k^{2})\mathcal{G}^{+}(x;y+e_{2}-e_{1})f(y)$	$\displaystyle=$	$\displaystyle f(x_{1}+1),$

  and as a result $(\Delta_{d}+k^{2})u(x)=0$. Thus, $u(x)$ is a solution to (2a). Since $\mathcal{G}^{+}$ satisfies the radiation condition and $f$ is supported on the finite subset of $\Gamma$ it follows that $u(x)$ is the unique radiating solution to the problem (2a), (2b). $\Box$

The main issue for numerical evaluation of the solution (12) is to compute the lattice Green’s function. For this purpose we apply the method developed in [4]. Using 8-fold symmetry, we need only to compute the lattice Green’s function $\mathcal{G}(i,j)$ with $i\geq j\geq 0$. Following to [4], let us introduce the vectors $\mathcal{V}_{2p}=(\mathcal{G}(2p,0),\mathcal{G}(2p-1,1),\dots,\mathcal{G}(p,p))^{\top}$ and $\mathcal{V}_{2p+1}=(\mathcal{G}(2p+1,0),\mathcal{G}(2p,1),\dots,\mathcal{G}(p+1,p))^{\top}$ that collect all distinct Green’s functions $\mathcal{G}(i,j)$ with “Manhattan distances” $\mid i\mid+\mid j\mid$ of $2p$ and $2p+1$, respectively. For any Manhattan distance larger than 1, equation

can be written in the matrix form $\gamma_{n}(k)\mathcal{V}_{n}=\alpha_{n}(k)\mathcal{V}_{n-1}+\beta_{n}(k)\mathcal{V}_{n+1}$ where $\alpha_{n}(k)$, $\beta_{n}(k)$ and $\gamma_{n}(k)$ are sparse matrices (cf., Appendix A). Notice that only the dimensions of these matrices depend on $n$. It is shown in [4] that, for any $n\geq 1$, we have

where the matrices $A_{n}(k)$ are defined by the following recurrence formula

They can be computed starting from a sufficiently large $N$ with $A_{N+1}(k)=0$. Here, it is worth mentioning, that for $k=2$ we need to choose a better “initial guess” than $A_{N+1}(k)=0$, since in this case $\mathrm{det}\gamma_{n}(k)=0$ and the matrix $\gamma_{n}(k)-\beta_{n}(k)A_{n+1}$ is not invertible.

Once $A_{n}(k)$ are known, we have $\mathcal{V}_{n}=A_{n}(k)\dots A_{1}(k)\mathcal{V}_{0}$, where $\mathcal{V}_{0}=\mathcal{G}(0,0)$. In particular, $\mathcal{V}_{1}=\mathcal{G}(1,0)=A_{1}(k)\mathcal{G}(0,0)$ which, together with $6\mathcal{G}(1,0)-(6-k^{2})\mathcal{G}(0,0)=1$, gives $\mathcal{G}(0,0)=1/[6A_{1}(k)-6+k^{2}]$. This completes the calculation of the Green’s function using elementary operations and no integrals. Notice also one more important advantage of this method. The $A_{n}(k)$ matrices are calculated coming down from asymptotically large Manhattan distances. As they are propagated towards smaller Manhattan distances, it definitely gives us the physical solution.

Finally, we demonstrate our theoretical and numerical approaches for the diffraction problem on triangular lattice. For this purpose we take $k=\sqrt{2}$ and consider the wave $U(x_{1},x_{2},t)=\exp(\iota x_{2}\pi/3-\iota\omega t)$ on the semi-infinite lattice $\mathbb{Z}\times\mathbb{Z}^{-}$ which encounters an obstacle at $x_{2}=0$ with two small openings formed by four nodes $\{(-11,0),(-10,0),(10,0),(11,0)\}$. It is easy to check that $\exp(\iota x_{2}\pi/3)$ satisfies the discrete Helmholtz equation. Thus, our goal is to evaluate numerically a radiating solution of Problem $\mathcal{P}_{\mathrm{D}}$ where, for the given data, we take $f(\pm 10)=f(\pm 11)=1$ and $f(m)=0$, $m\in\mathbb{Z}^{+}\backslash\{-11,-10,10,11\}$.

The results of numerical evaluations are plotted in Figure 2 in original coordinates of $\mathfrak{T}$, where (a) shows the density plots of the Green’s function $\mathrm{\Re e}\,\mathcal{G}$, (b) shows the density plots of $\cos(\iota x_{2}\pi/3)$ on the lower half-plane and $\mathrm{\Re e}\,u$ for the case $k=\sqrt{2}$ on the upper half-plane, and (c) presents the graph of $\mathrm{Re}\,u(\cdot,25\sqrt{3})$. Some key features of the numerical solution can be immediately observed, namely, as expected, the symmetry of $\mathrm{\Re e}\,u$ and contributions from the wavefront that create a variable intensity.

In this paper, we have constructed the discrete scattering theory for the two-dimensional discrete Helmholtz equation with a real wave number $k\in(0,3)\backslash\{2\sqrt{2}\}$ for the semi-infinite triangular lattices. The main objective was to prove the unique solvability result and derive an exact formula for the solution to a Dirichlet problem. For simplicity we restricted ourselves to compact boundary data. For non-compact boundary data we should require some extra conditions at infinity, cf. [15].

Similarly to the continuum theory, we used the notion of the radiating solution for the continuous Helmholtz equation

Recall that a solution $u$ to the equation (13) is called radiating if it satisfies the Sommerfeld radiation condition

for $\mid x\mid\to\infty$ uniformly in all directions $x/\mid x\mid$. There are evident similarities and differences between conditions (8), (9) and (14). For example, we see the same decay at infinity, but, in contrast to (14), the second condition of (8) or (9) is anisotropic: the factor $e^{\iota\xi^{*}_{j}(\alpha,k)}$ is not a constant. Moreover, it is well known that the radiating solutions $u$ to the continuous Helmholtz equation automatically satisfy the Sommerfeld’s finiteness condition, however this fact is not shown for the discrete case so far. Therefore, the first condition of (8) or (9) is included in our definition.

In the present paper, the problems under consideration have an infinite boundary. Within the continuum framework it is well-known that, in general, when the surface is unbounded, we cannot neglect the contribution of that surface waves at infinity. In this case, the Sommerfeld radiation condition is no longer appropriate and a proper modification is needed. To the best of the authors’ knowledge, different radiation conditions are provided only for the half-plane (and locally perturbed half-plane) problems, cf. [7, 8, 9, 10, 11, 12], and finding radiation conditions for arbitrary wedge-shaped regions remains open. In case of the square lattice we have proposed sufficient conditions for the given boundary data at infinity, which ensures to have an unique radiation solution to the corresponding problem, cf. [15]. Here, in order to avoid further technical difficulties, we restricted ourselves with compactly supported data on the boundary.

Finally, let us note that comparing the results of problems in the continuous framework and results obtained in the continuum limit of discrete problems deserves the high interest of scientists, however, it is beyond the scope of this paper.

This work was supported by Shota Rustaveli National Science Foundation of Georgia (SRNSFG) [FR-21-301]

The sparse matrices $\alpha_{n}(k)$, $\beta_{n}(k)$ and $\gamma_{n}(k)$ are defined as follows: if $n=2p$ then $\alpha_{2p}(k)$ is a $(p+1)\times p$ matrix such that $\alpha_{2p}(k)\mid_{i,i}=1$, $i=\overline{1,p}$, $\alpha_{2p}(k)\mid_{i,i-1}=1$, $i=\overline{2,p}$, while $\alpha_{2p}(k)\mid_{p+1,p}=2$, and all other matrix elements are zero. The $\beta_{2p}(k)$ is a $(p+1)\times(p+1)$ matrix such that $\beta_{2p}(k)\mid_{i,i}=1$, $i=\overline{1,p}$, $\beta_{2p}(k)\mid_{i,i+1}=1$, $i=\overline{2,p}$, while $\beta_{2p}(k)\mid_{p+1,p+1}=\beta_{2p}(k)\mid_{1,2}=2$, and all other matrix elements are zero. The $\gamma_{2p}(k)$ is a $(p+1)\times(p+1)$ matrix such that $\gamma_{2p}(k)\mid_{i,i}=6-k^{2}$, $i=\overline{1,p+1}$, $\gamma_{2p}(k)\mid_{i,i+1}=\gamma_{2p}(k)\mid_{i,i-1}=-1$, $i=\overline{2,p}$, and $\gamma_{2p}(k)\mid_{1,2}=\gamma_{2p}(k)\mid_{p+1,p}=-2$.

If $n=2p+1$ then $\alpha_{2p+1}(k)$ is a $(p+1)\times(p+1)$ matrix such that $\alpha_{2p+1}(k)\mid_{i,i}=1$, $i=\overline{1,p+1}$, $\alpha_{2p+1}(k)\mid_{i,i-1}=1$, $i=\overline{2,p+1}$, and all other matrix elements are zero. The $\beta_{2p+1}(k)$ is a $(p+1)\times(p+2)$ matrix such that $\beta_{2p+1}(k)\mid_{i,i}=1$, $i=\overline{1,p+1}$, $\beta_{2p+1}(k)\mid_{i,i+1}=1$, $i=\overline{2,p+1}$, while $\beta_{2p+1}(k)\mid_{1,2}=2$, and all other matrix elements are zero. The $\gamma_{2p}(k)$ is a $(p+1)\times(p+1)$ matrix such that $\gamma_{2p+1}(k)\mid_{i,i}=6-k^{2}$, $i=\overline{1,p}$, $\gamma_{2p+1}(k)\mid_{p+1,p+1}=5-k^{2}$, while $\gamma_{2p+1}(k)\mid_{i,i+1}=-1$, $i=\overline{2,p}$, $\gamma_{2p+1}(k)\mid_{i,i-1}=-1$, $i=\overline{2,p+1}$, and $\gamma_{2p+1}(k)\mid_{1,2}=\gamma_{2p+1}(k)\mid_{p+1,p}=-2$. Finally, $\gamma_{1}(k)$ is a $1\times 1$ matrixs with an element $4-k^{2}$.