Synchronous oscillations for a coupled cell-bulk ODE–PDE model with localized cells on ${\mathbb {R}}^2$

Abstract

In many micro- and macro-scale systems, collective dynamics occurs from the coupling of small spatially segregated, but dynamically active, units through a bulk diffusion field. This bulk diffusion field, which is both produced and sensed by the active units, can trigger and then synchronize oscillatory dynamics associated with the individual units. In this context, we analyze diffusion-induced synchrony for a class of cell-bulk ODE–PDE system in ${\mathbb {R}}^2$ that has two spatially segregated dynamically active circular cells of small radius. By using strong localized perturbation theory in the limit of small cell radius, we calculate the steady-state solution and formulate the linear stability problem. For Sel’kov intracellular reaction kinetics, we analyze how the effect of bulk diffusion can trigger, via a Hopf bifurcation, either in-phase or anti-phase intracellular oscillations that would otherwise not occur for cells that are uncoupled from the bulk medium. Phase diagrams in parameter space where these oscillations occur are presented, and the theoretical results from the linear stability theory are validated from full numerical simulations of the ODE–PDE system. In addition, the two-cell case is extended to study the onset of synchronous oscillatory instabilities associated with an infinite hexagonal arrangement of small identical cells in ${\mathbb {R}}^2$ with Sel’kov intracellular kinetics. This analysis for the hexagonal cell pattern relies on determining a new, computationally efficient, explicit formula for the regular part of a certain periodic reduced-wave Green’s function.

Appendices

Appendices

The quasi-periodic reduced-wave Green’s function

We develop a rapidly converging expansion for the Bloch Green’s function $G_\mathrm{{b}}(\mathbf{x})$ and its regular part $R_\mathrm{{b}}({\pmb k};z)$ as defined by (4.6). To do so, we extend the analysis in §5.1 of [20] for the case $z=1$, as was motivated by the methodology in [19], to allow for complex-valued z with $\text{ Re }(z)>0$. Then, by setting $z=1+\tau \lambda$ we obtain $R_\mathrm{{b}}({\pmb k};1+\tau \lambda )$, as needed in (4.5). Moreover, we identify that $R_p\equiv R_\mathrm{{b}}(\pmb {0};1)$ in (4.3) and (4.4).

We begin by representing the solution to (4.6) as the sum of free-space Green’s functions

$$\begin{aligned} G_\mathrm{{b}}(\mathbf{x}) = \sum _{\pmb l\in \Lambda } G_{\mathrm{free}}(\mathbf{x}+\pmb l)\, \mathrm{{e}}^{\mathrm{{i}}\pmb k\cdot \pmb l}, \end{aligned}$$

(A.1)

which ensures that the quasi-periodicity condition in (4.6) is satisfied. Then, to analyze (A.1) we use the Poisson summation formula which converts a sum over the hexagonal lattice $\Lambda$ to a sum over the reciprocal lattice $\Lambda ^{\star }$ of (4.2). In the notation of [19], we have (see Proposition 2.1 of [19])

$$\begin{aligned} \sum _{\pmb l\in \Lambda } f(\mathbf{x}+ \pmb l)\,\mathrm{{e}}^{\mathrm{{i}}\pmb k\cdot \pmb l} = \frac{1}{|\Omega |} \sum _{\pmb d\in \Lambda ^{\star }} {\hat{f}}(2\pi \pmb d - \pmb k) \,\mathrm{{e}}^{\mathrm{{i}}\mathbf{x}\cdot (2\pi \pmb d-\pmb k)}, \quad \mathbf{x}, \pmb k\in {\mathbb {R}}^2. \end{aligned}$$

(A.2)

Here $|\Omega |$ is the area of the primitive cell of the lattice, while ${\hat{f}}$ is the Fourier transform of f, defined in ${\mathbb {R}}^2$ by

$$\begin{aligned} {\hat{f}}(\pmb p) \equiv \int _{{\mathbb {R}}^2} f(\mathbf{x})\,\mathrm{{e}}^{-\mathrm{{i}}\mathbf{x}\cdot \pmb p}\, \mathrm{d}\mathbf{x}, \qquad f(\mathbf{x}) = \frac{1}{4\pi ^2}\int _{{\mathbb {R}}^2} {\hat{f}}(\pmb p) \,\mathrm{{e}}^{\mathrm{{i}}\pmb p\cdot \mathbf{x}}\, \mathrm{d}\pmb p. \end{aligned}$$

(A.3)

Upon applying (A.2) to (A.1) we obtain that the sum over the reciprocal lattice consists of free-space Green’s functions in the Fourier domain. By taking the Fourier transform of the PDE $\Delta G_{\mathrm{free}} - z D^{-1} G_{\mathrm{free}}=-\delta (\mathbf{x})$ for the free-space Green’s function, we obtain that ${\hat{G}}_{\mathrm{free}} (\pmb p) = {\hat{G}}_{\mathrm{free}} (|\pmb p|)$, where

$$\begin{aligned} {\hat{G}}_{\mathrm{free}}(\rho ) = \frac{1}{\rho ^2+\frac{z}{D}} \qquad \text{ with } \quad \rho \equiv |\pmb p|. \end{aligned}$$

(A.4)

Then, from (A.2) and (A.1), and by using $|\Omega |=1$, we obtain that

$$\begin{aligned} G_\mathrm{{b}}(\mathbf{x}) = \frac{1}{|\Omega |}\sum _{\pmb d\in \Lambda ^{\star }} {\hat{G}}_{\mathrm{free}}(2\pi \pmb d - \pmb k)\,\mathrm{{e}}^{\mathrm{{i}}\mathbf{x}\cdot (2\pi \pmb d-\pmb k)} = \sum _{\pmb d\in \Lambda ^{\star }} \frac{\mathrm{{e}}^{\mathrm{{i}}\mathbf{x}\cdot (2\pi \pmb d-\pmb k)} }{|2\pi \pmb d - \pmb k|^2+\frac{z}{D}}. \end{aligned}$$

(A.5)

In order to obtain a rapidly converging infinite series representation for $G_\mathrm{{b}}(\mathbf{x})$, we introduce the decomposition

$$\begin{aligned} {\hat{G}}_{\mathrm{free}}(2\pi \pmb d-\pmb k) = \alpha (2\pi \pmb d-\pmb k,\eta )\, {\hat{G}}_{\mathrm{free}}(2\pi \pmb d-\pmb k) + \left( 1-\alpha (2\pi \pmb d-\pmb k,\eta )\right) \, {\hat{G}}_{\mathrm{free}}(2\pi \pmb d-\pmb k), \end{aligned}$$

(A.6)

where the function $\alpha (2\pi \pmb d-\pmb k,\eta )$ is defined by

$$\begin{aligned} \alpha (2\pi \pmb d-\pmb k,\eta ) = \exp \!\left( -\frac{|2\pi \pmb d-\pmb k|^2+\frac{z}{D}}{\eta ^2}\right) . \end{aligned}$$

(A.7)

Here $\eta >0$ is a real-valued cut-off parameter, which is specified below. Since $\text{ Re }(z)>0$, we readily calculate that

$$\begin{aligned} \lim _{\eta \rightarrow 0}\alpha (2\pi \pmb d-\pmb k,\eta ) = 0; \qquad \lim _{\eta \rightarrow \infty }\alpha (2\pi \pmb d-\pmb k,\eta ) = 1. \end{aligned}$$

With this choice for $\alpha$, the sum over $\pmb d\in \Lambda ^{\star }$ in (A.5) of the first set of terms in (A.6), as labeled by

$$\begin{aligned} G_{\mathrm{fourier}}(\mathbf{x}) \equiv \sum _{\pmb d\in \Lambda ^{\star }} \exp \!\left( -\frac{|2\pi \pmb d-\pmb k|^2+ \frac{z}{D}}{\eta ^2}\right) \, \frac{\,\mathrm{{e}}^{\mathrm{{i}}\mathbf{x}\cdot (2\pi \pmb d-\pmb k)} }{|2\pi \pmb d - \pmb k|^2+\frac{z}{D}}, \end{aligned}$$

(A.8)

converges absolutely when $\text{ Re }(z)>0$.

Next, we calculate the lattice sum in (A.5) over the second set of terms in (A.6) by using the inverse transform (A.3) after first writing $(1-\alpha )\,{\hat{G}}_{\mathrm{free}}$ as an integral. To do so, we define $\rho \equiv |2\pi \pmb d-\pmb k|$, so that from (A.7) and (A.4)

$$\begin{aligned} \left( 1- \alpha (2\pi \pmb d-\pmb k,\eta )\right) \, {\hat{G}}_{\mathrm{free}}(2\pi \pmb d-\pmb k)= & {} \frac{1}{\rho ^2 + \frac{z}{D}} \left( 1-\exp \!\left( -\frac{\rho ^2+\frac{z}{D}}{\eta ^2}\right) \right) \nonumber \\= & {} {\hat{\chi }}(\rho ) \equiv \int _{\log \eta }^\infty \!\! 2\,\mathrm{{e}}^{-\left( \rho ^2+\frac{z}{D}\right) \,\mathrm{{e}}^{-2s}-2s}\, \mathrm{{d}}s. \end{aligned}$$

(A.9)

In recognizing the middle term in (A.9) as the integral ${\hat{\chi }}(\rho )$ we used the easily verified definite integral

$$\begin{aligned} 2\!\!\int _{\log \eta }^{\infty }\!\! \mathrm{{e}}^{-\left( \rho ^2 + \frac{z}{D}\right) \mathrm{{e}}^{-2s}-2s}\,\mathrm{{d}}s =\frac{\mathrm{{e}}^{-\left( \rho ^2 + \frac{z}{D}\right) \mathrm{{e}}^{-2s}}}{\rho ^2+ \frac{z}{D} } \Big \vert _{s=\log \eta }^{s=\infty } = \frac{1}{\rho ^2 + \frac{z}{D}} \left( 1-\exp \!\left( -\frac{\rho ^2+\frac{z}{D}}{\eta ^2}\right) \right) . \end{aligned}$$

(A.10)

Next, we take the inverse Fourier transform of (A.9). To do so, we use two key facts. Firstly, the inverse Fourier transform of a radially symmetric function ${\hat{f}}(\rho )$ is the inverse Hankel transform of order zero (cf. [35]), so that $f(r)=(2\pi )^{-1}\int _0^\infty {\hat{f}}(\rho )\,J_0(\rho r)\,\rho \, \mathrm{{d}}\rho$. Secondly, from [35], we recall the well-known inverse Hankel transform

$$\begin{aligned} \int _0^\infty \mathrm{{e}}^{-\rho ^2\,\mathrm{{e}}^{-2s}} \, \rho \, J_0(\rho r)\, \mathrm{{d}}\rho = \frac{1}{2}\,\mathrm{{e}}^{2s-{r^2\,\mathrm{{e}}^{2s}/4}}. \end{aligned}$$

In this way, we calculate using the definition of ${\hat{\chi }}(\rho )$ in (A.9) that

$$\begin{aligned} \chi (r)\equiv \frac{1}{2\pi }\int _0^\infty {\hat{\chi }}(\rho )\,J_0(\rho r)\,\rho \, \mathrm{{d}}\rho&=\frac{1}{\pi } \!\!\int ^{\infty }_{\log \eta }\!\! \mathrm{{e}}^{-2s - z \mathrm{{e}}^{-2s}/D} \left( \int _0^\infty \mathrm{{e}}^{-\rho ^2 \mathrm{{e}}^{-2s}}\, \rho \, J_0(\rho r) \, \mathrm{{d}}\rho \,\right) \,\mathrm{{d}}s\\&=\frac{1}{2\pi } \!\!\int ^{\infty }_{\log \eta }\!\! \mathrm{{e}}^{-2s - z \mathrm{{e}}^{-2s}/D} \,\mathrm{{e}}^{2s- \frac{r^2}{4}\,\mathrm{{e}}^{2s}} \,\mathrm{{d}}s = \frac{1}{2\pi } \!\!\int ^{\infty }_{\log \eta }\!\! \,\mathrm{{e}}^{-\left( \frac{r^2}{4}\,\mathrm{{e}}^{2s} + \frac{z}{D}\,\mathrm{{e}}^{-2s}\right) }\,\mathrm{{d}}s. \end{aligned}$$

In the notation of [19], we then define $F_{\mathrm{sing}}(|\mathbf{x}|)$ as

$$\begin{aligned} F_{\mathrm{sing}}(|\mathbf{x}|)\equiv \frac{1}{2\pi } \!\!\int ^{\infty }_{\log \eta }\!\! \,\mathrm{{e}}^{-\left( \frac{|\mathbf{x}|^2}{4}\,\mathrm{{e}}^{2s} + \frac{z}{D}\,\mathrm{{e}}^{-2s}\right) }\,\mathrm{{d}}s. \end{aligned}$$

(A.11)

Therefore, by using the Poisson summation formula (A.2) to calculate lattice sum in (A.5) over the second set of terms in (A.6), we obtain

$$\begin{aligned} G_{\mathrm{spatial}}(\mathbf{x}) \equiv F_{\mathrm{sing}}(|\mathbf{x}|) + \sum _{\genfrac{}{}{0.0pt}{}{\pmb l\in \Lambda }{\pmb l\ne \mathbf{0}}}\mathrm{{e}}^{\mathrm{{i}}\pmb k\cdot \pmb l}\, F_{\mathrm{sing}}(|\mathbf{x}+ \pmb l|). \end{aligned}$$

(A.12)

In summary, the Bloch Green’s function for the reduced-wave operator in the spatial domain, satisfying (4.6), is $G_\mathrm{{b}}(\mathbf{x})=G_{\mathrm{fourier}}(\mathbf{x}) + G_{\mathrm{spatial}}(\mathbf{x})$, representing the sum of (A.8) and (A.12). In this way, we have

$$\begin{aligned} G_\mathrm{{b}}(\mathbf{x}) = \sum _{\pmb d\in \Lambda ^{\star }} \frac{ \mathrm{{e}}^{-\frac{|2\pi \pmb d-\pmb k|^2+\frac{z}{D}}{\eta ^2}}\, \,\mathrm{{e}}^{\mathrm{{i}}\mathbf{x}\cdot (2\pi \pmb d-\pmb k)} }{|2\pi \pmb d - \pmb k|^2+ \frac{z}{D}} + F_{\mathrm{sing}}(|\mathbf{x}|) + \sum _{\genfrac{}{}{0.0pt}{}{\pmb l\in \Lambda }{\pmb l\ne \mathbf{0}}} \mathrm{{e}}^{\mathrm{{i}}\pmb k\cdot \pmb l}\, F_{\mathrm{sing}}(|\mathbf{x}+ \pmb l|), \end{aligned}$$

(A.13)

where $F_{\mathrm{sing}}(|\mathbf{x}|)$ is defined in (A.11). In terms of the Ewald cut-off parameter $\eta$, we observe from (A.8) that $G_{\mathrm{fourier}}(\mathbf{x}) \rightarrow 0$ as $\eta \rightarrow 0$, while from (A.12) and (A.11) we get that $G_{\mathrm{spatial}}(\mathbf{x}) \rightarrow 0$ as $\eta \rightarrow \infty$.

The next step in the analysis is to isolate the logarithmic singularity in (A.13) as $\mathbf{x}\rightarrow \mathbf{0}$, so as to identify the regular part $R_\mathrm{{b}}({\pmb k};z)$ defined by

$$\begin{aligned} R_\mathrm{{b}}({\pmb k};z) \equiv \lim _{\mathbf{x}\rightarrow \mathbf{0}} \left( G_\mathrm{{b}}(\mathbf{x}) + \frac{1}{2\pi } \log |\mathbf{x}| \right) . \end{aligned}$$

(A.14)

From (A.8) we observe that $G_{\mathrm{fourier}}(\mathbf{0})$ is finite for all $|2\pi \pmb d -\pmb k|$ and $0<\eta <\infty$ when $\text{ Re }(z)>0$. Likewise, for the last set of terms in (A.13), we can take the limit $\mathbf{x}\rightarrow \mathbf{0}$ to conclude that $\sum _{\genfrac{}{}{0.0pt}{}{\pmb l\in \Lambda }{\pmb l\ne \mathbf{0}}}\mathrm{{e}}^{\mathrm{{i}}\pmb k\cdot \pmb l}\, F_{\mathrm{sing}}(|\pmb l|)$ is also finite. As a result, the logarithmic singularity in $G_\mathrm{{b}}$ as $\mathbf{x}\rightarrow 0$ must arise from the middle term, $F_{\mathrm{sing}}(r)$, in (A.13), where $r\equiv |\mathbf{x}|$.

To analyze $F_{\mathrm{sing}}(r)$ when $\text{ Re }(z)>0$, we write $z=|z|\mathrm{{e}}^{\mathrm{{i}}\theta }$, where $\theta \equiv \text{ Arg }z$ satisfies $|\theta |<{\pi /2}$. In the integrand of (A.11) we introduce the new variable t by

$$\begin{aligned} s=-\frac{1}{4}\log \left( \frac{r^2 D}{4|z|}\right) + \frac{t}{2} \end{aligned}$$

(A.15)

so that (A.11) becomes

$$\begin{aligned} F_{\mathrm{sing}}(r) = \frac{1}{4\pi } \int _{\beta }^{\infty } \exp \left( - r \sqrt{\frac{z}{D}} \cosh \left( t-{\mathrm{{i}}\theta /2}\right) \right) \, \mathrm{{d}}t, \qquad \text{ where } \qquad \beta \equiv \log \left( \frac{\eta ^2 r \sqrt{D}}{2\sqrt{|z|}}\right) . \end{aligned}$$

(A.16)

Here $\theta \equiv \text{ Arg }\,z$ and we have specified the principal branch for $\sqrt{z}$. We then split the integration range in (A.16) as

$$\begin{aligned} F_{\mathrm{sing}}(r) = \frac{1}{4\pi } \int _{-\infty }^{\infty } \exp \left( - r \sqrt{\frac{z}{D}} \cosh \left( t-{\mathrm{{i}}\theta /2}\right) \right) \, \mathrm{{d}}t - \frac{1}{4\pi } \int _{-\infty }^{\beta } \exp \left( - r \sqrt{\frac{z}{D}} \cosh \left( t-{\mathrm{{i}}\theta /2}\right) \right) \, \mathrm{{d}}t. \end{aligned}$$

(A.17)

To calculate the first term in (A.17), we use a contour integration over the box-shaped contour $-b\le \text{ Re }(t)\le b$ and $0\le \text{ Im }(t)<{\pi /2}$ in the complex t-plane. Since there are no residues within the contour, and we have exponential decay on the sides of the box as $b\rightarrow \infty$, we can replace $t-{\mathrm{{i}}\theta /2}$ by t in (A.17) and then use symmetry to obtain

$$\begin{aligned} \frac{1}{4\pi } \int _{-\infty }^{\infty } \exp \left( - r \sqrt{\frac{z}{D}} \cosh \left( t-{\mathrm{{i}}\theta /2}\right) \right) \, \mathrm{{d}}t= & {} \frac{1}{2\pi } \int _{0}^{\infty } \exp \left( - r \sqrt{\frac{z}{D}} \cosh \left( t\right) \right) \, \mathrm{{d}}t\nonumber \\= & {} \frac{1}{2\pi } \int _{1}^{\infty } \frac{\mathrm{{e}}^{-r \xi \sqrt{{z/D}}}}{\sqrt{\xi ^2-1}} \, \mathrm{{d}}\xi , \end{aligned}$$

(A.18)

where the last equality in (A.18) follows from the substitution $\xi =\cosh (t)$. The last integral in (A.18) is identified by using the well-known integral representation $K_{0}(\omega ) = \int _{1}^{\infty } {\mathrm{{e}}^{-\omega \xi }/\sqrt{\xi ^2-1}} \, \mathrm{{d}}\xi$ for the modified Bessel function of the second kind of order zero, which is valid for $|\text{ Arg }\,\omega |<{\pi /2}$. In this way, we obtain from (A.18) and (A.17) that

$$\begin{aligned} F_{\mathrm{sing}}(r) = \frac{1}{2\pi } K_0\left( r \sqrt{ \frac{z}{D}} \right) -\frac{J(r;z)}{4\pi }, \quad \text{ where } \quad J(r;z) \equiv \int _{-\infty }^{\beta } \exp \left( - r \sqrt{\frac{z}{D}} \cosh \left( t-{\mathrm{{i}}\theta /2}\right) \right) \, \mathrm{{d}}t. \end{aligned}$$

(A.19)

Next, we provide a more tractable representation for J(r; z). We introduce a new variable $\xi$ defined by

$$\begin{aligned} \xi = \mathrm{{e}}^{-(t-\beta )} = \frac{\eta ^2 r \sqrt{D}}{2\sqrt{|z|}} \mathrm{{e}}^{-t} \end{aligned}$$

(A.20)

so that $\xi =1$ when $t=\beta$ and $\xi \rightarrow +\infty$ as $t\rightarrow -\infty$. In terms of $\xi$, the exponential in the integral J is

$$\begin{aligned} -r \sqrt{\frac{z}{D}} \cosh \left( t-{\mathrm{{i}}\theta /2}\right) = - \frac{z \xi }{\eta ^2 D} \left( 1 + \frac{\eta ^4 r^2 D}{4 z\xi ^2} \right) . \end{aligned}$$

(A.21)

By using (A.21) and $\mathrm{{d}}t=-{\mathrm{{d}}\xi /\xi }$ within the integral J in (A.19), we obtain that (A.19) transforms to

$$\begin{aligned} F_{\mathrm{sing}}(r) = \frac{1}{2\pi } K_0\left( r \sqrt{ \frac{z}{D}} \right) -\frac{J(r;z)}{4\pi }, \qquad \text{ where } \qquad J(r;z) \equiv \int _{1}^{\infty } \frac{\mathrm{{e}}^{-z\xi /(\eta ^2 D)}}{\xi } \mathrm{{e}}^{-r^2\eta ^2/(4\xi )} \, \mathrm{{d}}\xi . \end{aligned}$$

(A.22)

Finally, we use (A.22) to extract the local behavior of $F_{\mathrm{sing}}(r)$ as $r\rightarrow 0$. As $r\rightarrow 0$, we calculate for J(r; z) that

$$\begin{aligned} J = E_1\left( \frac{z}{\eta ^2 D} \right) + {{\mathcal {O}}}(r^2) \quad \quad \text{ when } \quad \text{ Re }(z)>0, \end{aligned}$$

(A.23)

where $E_1(w)\equiv \int _{1}^{\infty } \xi ^{-1} \mathrm{{e}}^{-w\xi }\, \mathrm{{d}}\xi$, for $\text{ Re }(w)>0$, is the well-known exponential integral. In addition, in (A.22) we use $K_0(w) = -\log {w}+\log {2}-\gamma _e+o(1)$ as $|w|\rightarrow 0$, where $\gamma _e$ is Euler’s constant. In this way, (A.22) and (A.23) yield

$$\begin{aligned} F_{\mathrm{sing}}(r) = -\frac{1}{2\pi } \log {r} + \frac{1}{2\pi } \left[ \log (2\sqrt{D}) - \gamma _e - \log (\sqrt{z})\right] - \frac{1}{4\pi } E_1\left( \frac{z}{\eta ^2 D} \right) + o(1) \qquad \text{ as } \quad r\rightarrow 0. \end{aligned}$$

(A.24)

Finally, by substituting (A.13) in (A.14), and using (A.24), we obtain the result for $R_\mathrm{{b}}({\pmb k};z)$ given in (4.7).

Formulation on the Wigner–Seitz cell

In this appendix, we provide a key result for $R_\mathrm{{b}}({\pmb k};z)$. However, before doing so, we first must obtain a more refined description of the fundamental Wigner–Seitz (FWS) cell, as was discussed in §2.2 of [31]. For a general lattice, there are eight nearest neighbor lattice points to $\mathbf{x}=0$ given by the set

$$\begin{aligned} P\equiv \lbrace { \, m{\pmb l}_1+ n{\pmb l}_2 \, \vert \,\, m\in \lbrace {0,1,-1\rbrace }, \,\, n\in \lbrace {0,1,-1\rbrace },\,\, (m,n)\ne 0 \rbrace }. \end{aligned}$$

(B.1)

For each (vector) point ${\pmb P}_i\in P$, for $i=1,\ldots ,8$, the Bragg line $L_i$ is defined as the line that crosses the point ${{\pmb P}_i/2}$ orthogonally to ${\pmb P}_i$. The unit outer normal to $L_i$ is labeled by ${\pmb \eta }_i\equiv {{\pmb P}_i/|{\pmb P}_i|}$. The convex hull generated by these Bragg lines is the FWS cell $\Omega$. Specifically, for the hexagonal lattice (4.1) its boundary $\partial \Omega$ is the union of exactly six Bragg lines, and the centers of the Bragg lines generating $\partial \Omega$ are re-indexed as ${{\pmb P}_i/2}$ for $i=1,\ldots ,6$. The boundary $\partial \Omega$ of $\Omega$ is the union of the re-indexed Bragg lines $L_i$ for $i=1,\ldots ,6$, and is parameterized segment-wise as

$$\begin{aligned} \partial \Omega = \left\{ \mathbf{x}\in \bigcup _i \left\{ \frac{{\pmb P}_i}{2} + t {\pmb \eta }_{i}^{\perp }\right\} \,\, \Big \vert \,\,\, -t_i\le t\le t_i, \,\,\, i=1,\ldots ,6\right\} , \end{aligned}$$

(B.2)

where $2t_i$ is the length of $L_i$. Here ${\pmb \eta }_{i}^{\perp }$ is the direction perpendicular to ${\pmb P}_i$, and is, therefore, tangent to $L_i$. From this construction, Bragg lines on $\partial \Omega$ must come in pairs. In particular, if ${\pmb P}$ is a neighbor of $0$ and the Bragg line crossing ${{\pmb P}/2}$ lies on $\partial \Omega$, it follows by symmetry that the Bragg line crossing ${-{\pmb P}/2}$ must also lie on $\partial \Omega$.

Next, we reformulate the PDE (4.6) for the reduced-wave Bloch Green’s function $G_\mathrm{{b}}$ on ${\mathbb {R}}^2$ to an equivalent PDE on the FWS $\Omega$. This is done by imposing a boundary operator $\mathcal{{P}}_k$ on $\partial \Omega$ that incorporates the quasi-periodic condition in (4.6). This equivalent PDE is

$$\begin{aligned} \Delta G_\mathrm{{b}} -\frac{z}{D} G_\mathrm{{b}} = -\delta (\mathbf{x}), \quad \mathbf{x}\in \Omega ; \qquad G_\mathrm{{b}} \in \mathcal{{P}}_k, \quad \mathbf{x}\in \partial \Omega ; \qquad R_\mathrm{{b}}({\pmb k};z) \equiv \lim _{\mathbf{x}\rightarrow \mathbf{0}} \left( G_\mathrm{{b}}(\mathbf{x}) + \frac{1}{2\pi } \log |\mathbf{x}| \right) , \end{aligned}$$

(B.3)

where ${{\pmb k}/(2\pi )}\in \Omega _\mathrm{{b}}$. In (B.3), the boundary operator is defined by

$$\begin{aligned} \mathcal{{P}}_{k} \equiv \left\{ u \, \Big \vert \, \left( \begin{array}{c} u(\mathbf{x}_{i1}) \\ \partial _n u(\mathbf{x}_{i1}) \end{array} \right) = \mathrm{{e}}^{-\mathrm{{i}}{\pmb k}\cdot {{\pmb l}_i}} \left( \begin{array}{c} u(\mathbf{x}_{i2}) \\ \partial _n u(\mathbf{x}_{i2}) \end{array} \right) , \,\, \forall \, \mathbf{x}_{i1}\in L_i, \,\, \forall \,\mathbf{x}_{i2}\in L_{-i}, \,\, {{\pmb l}}_i \in \Lambda , \,\, i=1,\ldots ,3 \right\} . \end{aligned}$$

(B.4)

Here $L_{i}$ and $L_{-i}$ denote two parallel Bragg lines on opposite sides of $\partial \Omega$ for $i=1,\ldots ,3$, while $\mathbf{x}_{i1}\in L_i$ and $\mathbf{x}_{i2}\in L_{-i}$ are any two opposing points on these Bragg lines. In this notation, periodic boundary conditions refer to $G_\mathrm{{b}}\in \mathcal{{P}}_{0}$.

With this reformulation of (4.6) to the PDE (B.3) on the FWS cell $\Omega$, with unit area $|\Omega |=1$, we now derive a result is needed in Sect. 4.2. Firstly, we obtain for the periodic problem, with $G_\mathrm{{b}}\in \mathcal{{P}}_{0}$ on $\partial \Omega$, that $G_\mathrm{{b}}={{\mathcal {O}}}(D)$ when $D\gg 1$. By expanding in powers of D, we readily derive for ${\pmb k}={\pmb 0}$ that

$$\begin{aligned} G_\mathrm{{b}} = \frac{D}{1+\tau \lambda } + G_{p0} + {{\mathcal {O}}}(D^{-1}), \end{aligned}$$

(B.5)

where $G_{p0}(\mathbf{x})$ is the periodic source-neutral Green’s function satisfying

$$\begin{aligned} \Delta G_{p0} = \frac{1}{|\Omega |} - \delta (\mathbf{x}), \quad \mathbf{x}\in \Omega ; \qquad G_\mathrm{{B}0} \in \mathcal{{P}}_0, \quad \mathbf{x}\in \partial \Omega ; \qquad R_{p0} \equiv \lim _{\mathbf{x}\rightarrow \mathbf{0}} \left( G_{p0}(\mathbf{x}) + \frac{1}{2\pi } \log |\mathbf{x}| \right) , \end{aligned}$$

(B.6)

and normalized by $\int _{\Omega } G_{p0} \, \mathrm{{d}}\mathbf{x}= 0$. From Theorem 1 of [36], $R_{p0}\approx -0.21027$ is given explicitly by

$$\begin{aligned} R_{p0}=-\frac{1}{2\pi }\ln (2\pi ) -\frac{1}{2\pi } \ln \Big \vert \sqrt{ \sin \left( \frac{\pi }{3}\right) } \, \mathrm{{e}}^{\pi \mathrm{{i}} \xi /6} \prod _{n=1}^{\infty } \left( 1- \mathrm{{e}}^{2\pi \mathrm{{i}} n\xi } \right) ^2\Big \vert , \qquad \text{ where } \qquad \xi \equiv \mathrm{{e}}^{\mathrm{{i}}\pi /3}. \end{aligned}$$

(B.7)

By taking the regular part of (B.5) we obtain the two-term result for $R_\mathrm{{b}}({\pmb 0};1+\tau \lambda )$ in (4.12), which is valid for $D\gg 1$.