A framework of shape optimisation based on the isogeometric boundary element method toward designing thin-silicon photovoltaic devices

Abstract

We propose a gradient-based optimisation framework to design the shape of the interfaces in periodic layered structures, which can model photovoltaic devices, optical gratings, etc., according to a desired electromagnetic property. To this end, we first develop an isogeometric boundary element method (IGBEM) to analyse singly periodic boundary value problems for 2D Helmholtz equation in layered media. By nature, the IGBEM is not only numerically accurate but also suitable for the shape optimisation in comparison with the conventional boundary- and domain-type solvers. Next, we derive the shape derivative (or shape sensitivity) of the objective function in terms of the magnetic field strength or energy absorption rate on the basis of the adjoint variable method. The shape derivative can be computed by solving the primal and adjoint problems with the IGBEM. Subsequently, we map our shape optimisation problems to nonlinear programming problems in order to exploit a general-purpose software for the latter problems. After verifying our optimisation framework rigorously by comparing with the exact solutions, we demonstrate a shape optimisation of the silicon-metal interface in a model of thin-silicon photovoltaic devices: it is crucial to determine the interface so that the energy of the incident light can be confined in the silicon layer as much as possible when the thickness of the layer is unconventionally small; one micrometre in our model. The optimal shape achieved 8.6 times higher absorption rate than the reference (flat) shape. This simulation shows that the proposed framework is capable of making a fundamental contribution to analysing and designing photovoltaic devices as well as photonic and plasmonic devices that consist of singly periodic layers.

Appendices

Appendix A: Evaluation of the quasi-periodic Green’s function by Kummer’s transformation

In this appendix, we will denote the quasi-periodic Green’s function $G_i$ and the wavenumber $k_i$ in Eq. 8 as $G$ and k, respectively, for concise. Also, we denote $\varvec{x}$ and $\varvec{y}\in \mathbb {R}^2$ as $(x_1,x_2)^T$ and $(y_1,y_2)^T$, respectively. In addition, we consider the case that k is real, because $G$ converges rapidly if k is not real.

Applying Kummer’s transformation [88] to $G$ leads to the two series $G^{1}$ and $G^{2}$ that can converge more rapidly than $G$, i.e.

$$\begin{aligned} G^{1}(\varvec{x}-\varvec{y}):= & {} \frac{1}{2\pi }\sum _{m=-\infty }^{\infty }K_0(k\left| \varvec{x}-(\varvec{y}+ mL\varvec{e}_x)\right| )e^{\mathrm {i}\beta L m}, \end{aligned}$$

(42a)

$$\begin{aligned} G^{2}(\varvec{x}-\varvec{y}):= & {} \frac{1}{2L}\sum _{m=-\infty }^{\infty }\left( \frac{e^{-\sqrt{\xi _m^2-k^2}|x_2-y_2|}}{\sqrt{\xi _m^2-k^2}} - \frac{e^{-\sqrt{\xi _m^2+k^2}|x_2-y_2|}}{\sqrt{\xi _m^2+k^2}} \right) e^{\mathrm {i}\xi _m(x_1-y_1)}, \end{aligned}$$

(42b)

where $G\equiv G^{1}+G^{2}$ holds, $K_0$ denotes the modified Bessel function of the second kind of order zero and $\xi _m:=\beta +\frac{2\pi m}{L}$. Note that, regarding any real wavenumber k, we evaluate $\sqrt{\xi _m^2-k^2}$ as follows:

$$\begin{aligned} \sqrt{\xi _m^2-k^2}= {\left\{ \begin{array}{ll} \sqrt{\xi _m^2-k^2} &{} ~\mathrm{for}~ \xi _m^2> k^2, \\ -\mathrm {i}\sqrt{k^2-\xi _m^2} &{} ~\mathrm{for}~ \xi _m^2 < k^2. \end{array}\right. } \end{aligned}$$

(43)

Note also that the term of $m=0$ and $-1$ in Eq. 42a can be singular under the configuration of the present IGBEM, as mentioned in Subsubsect. 2.3.5.

The convergence rate of the series $G^1$ is $\mathcal {O}(m^{-\frac{1}{2}}e^{-m})$, while that of $G^2$ is $\mathcal {O}(m^{-1}e^{-m})$ if $|x_2-y_2|>0$ and $\mathcal {O}(m^{-3})$ otherwise. These rates are genuinely better than that of $G$, i.e. $\mathcal {O}(m^{-1})$.

We need to truncate the infinite series of Eq. 42a and b. To this end, we compare the value of $G^1$ or $G^2$ evaluated for $\left| m\right| \le p$ for a certain number p with the previous value for $\left| m\right| \le p-1$. If the change of the former value relative to the latter one is smaller than $10^{-10}$, we employ the former value. Otherwise, we increase the number p by one and repeat the evaluation and comparison.

Regarding the double layer kernel, i.e. $\frac{\partial G(\varvec{x}-\varvec{y})}{\partial n}$, the normal derivative of $G^1$ and $G^2$ is computed as follows:

$$\begin{aligned}&\frac{\partial G^1(\varvec{x}-\varvec{y})}{\partial n(\varvec{y})} = \frac{k}{2\pi }\sum _{m=-\infty }^{\infty }K_1(k|\varvec{x}-(\varvec{y}+ mL\varvec{e}_x)|)\frac{(\varvec{x}-(\varvec{y}+ mL\varvec{e}_x))\cdot \varvec{n}(\varvec{y})}{|\varvec{x}-(\varvec{y}+ mL\varvec{e}_x)|}e^{\mathrm {i}\beta L m},\end{aligned}$$

(44a)

$$\begin{aligned}&\frac{\partial G^2(\varvec{x}-\varvec{y})}{\partial n(\varvec{y})} = \frac{1}{2L}\sum _{m=-\infty }^{\infty }\left( \left( -\mathrm {i}\xi _m n_1(\varvec{y})+\mathrm {sgn}(x_2-y_2)\sqrt{\xi _m^2-k^2}n_2(\varvec{y})\right) \frac{e^{-\sqrt{\xi _m^2-k^2}|x_2-y_2|}}{\sqrt{\xi _m^2-k^2}}\right. \nonumber \\&\left. -\left( -\mathrm {i}\xi _m n_1(\varvec{y})+\mathrm {sgn}(x_2-y_2)\sqrt{\xi _m^2+k^2}n_2(\varvec{y})\right) \frac{e^{-\sqrt{\xi _m^2+k^2}|x_2-y_2|}}{\sqrt{\xi _m^2+k^2}} \right) e^{\mathrm {i}\xi _m(x_1-y_1)}, \end{aligned}$$

(44b)

where $K_1$ denotes the modified Bessel function of the second kind of order one and $\mathrm {sgn}(f)$ denotes the sign of f. The way to truncate the infinite series in Eq. 44 is the same as $G^1$ and $G^2$.

Appendix B: B-spline basis function of degree 2

Letting $\{t_0,\ldots ,t_{n+2}\}$ be the knots, the B-spline basis function of degree 2, denoted by $N_\nu ^2$, and its derivative are given as follows:

$$\begin{aligned} N^2_\nu (t) = {\left\{ \begin{array}{ll} 0 &{} t< t_\nu \\ \frac{(t-t_\nu )^2}{(t_{\nu +2}-t_\nu )(t_{\nu +1}-t_\nu )} &{} t_\nu \le t< t_{\nu +1} \\ \frac{(t-t_\nu )(t_{\nu +2}-t)}{(t_{\nu +2}-t_\nu )(t_{\nu +2}-t_{\nu +1})} + \frac{(t_{\nu +3} - t)(t - t_{\nu +1})}{(t_{\nu +3} - t_{\nu +1})(t_{\nu +2} - t_{\nu +1})} &{} t_{\nu +1} \le t< t_{\nu +2}\\ \frac{(t_{\nu +3} - t)^2}{(t_{\nu +3} - t_{\nu +1})(t_{\nu +3} - t_{\nu +2})} &{} t_{\nu +2} \le t < t_{\nu +3} \\ 0 &{} t \ge t_{\nu +3},\\ \end{array}\right. } \end{aligned}$$

$$\begin{aligned} \frac{\mathrm {d}N^2_\nu (t)}{\mathrm {d}t} = {\left\{ \begin{array}{ll} 0 &{} t< t_\nu \\ \frac{2(t-t_\nu )}{(t_{\nu +2}-t_\nu )(t_{\nu +1}-t_\nu )} &{} t_\nu \le t< t_{\nu +1} \\ \frac{-2t+(t_\nu +t_{\nu +2})}{(t_{\nu +2}-t_\nu )(t_{\nu +2}-t_{\nu +1})} + \frac{-2t+(t_{\nu +1}+t_{\nu +3})}{(t_{\nu +3} - t_{\nu +1})(t_{\nu +2} - t_{\nu +1})} &{} t_{\nu +1} \le t< t_{\nu +2}\\ \frac{-2(t_{\nu +3} - t)}{(t_{\nu +3} - t_{\nu +1})(t_{\nu +3} - t_{\nu +2})} &{} t_{\nu +2} \le t < t_{\nu +3} \\ 0 &{} t \ge t_{\nu +3}. \end{array}\right. } \end{aligned}$$

Note that, when the denominator of a particular term in the RHSs vanishes, the term can be considered as zero.

Appendix C: Derivation of shape derivatives in Eqs. (28) and (31)

We derive shape derivative, i.e. the sensitivity of an objective function $\mathcal {J}$ with respect to the shape of $\Gamma$, in the framework of the adjoint variable method. To this end, we consider the following Lagrangian $\mathcal {L}$:

$$\begin{aligned} \mathcal {L}(\Gamma ;u,\lambda ):=\mathcal {J}(\Gamma ;u) + \sum _{i=0}^{d-1} \mathop {\mathrm {Re}}\left[ \int _{\Omega _i} \frac{\lambda ^*}{\varepsilon _i}( {{{\Delta }}}u + k_i^2 u) \, \mathrm {d}\Omega \right] , \end{aligned}$$

(45)

where the Lagrange multiplier $\lambda$ is an arbitrary function and the operator ${}^*$ denotes the complex conjugate. As long as u satisfies Eq. 2, $\mathcal {L}$ is equivalent to $\mathcal {J}$. Hence, we may maximise or minimise $\mathcal {L}$ instead of $\mathcal {J}$. Note that dividing $\lambda$ by $\varepsilon _i({}\ne 0)$ is not essential but necessary to obtain the adjoint problem in the same fashion as the primal problem. In what follows, we denote the integral over $\Omega _i$ in the RHS of Eq. 45 by $\mathcal {G}$, i.e.

$$\begin{aligned} \mathcal {G}(\Omega _i;u,\lambda ):=\int _{\Omega _i} \frac{\lambda ^*}{\varepsilon _i}( {{{\Delta }}}u + k_i^2 u) \, \mathrm {d}\Omega . \end{aligned}$$

This can be rewritten by the Green’s first identity as follows:

$$\begin{aligned} \mathcal {G}(\Omega _i;u,\lambda ) = \int _{\Omega _i} \frac{1}{\varepsilon _i} \left( \nabla \lambda ^* \cdot \nabla u - k_i^2 \lambda ^* u \right) \mathrm {d}\Omega - \int _{\partial \Omega _i} \lambda ^* q \, \mathrm {d}\Gamma . \end{aligned}$$

(46)

Next, we slightly move (oneperturbate) every point $\varvec{x}$ in $\Omega$ by $\epsilon \varvec{v}(\varvec{x})$. Here, a real number $\epsilon$ is infinitesimal and represents the degree of the onedesign perturbation (variation) and the vector $\epsilon \varvec{v}(\varvec{x})$ represents the onedesign velocity at $\varvec{x}$. Then, the domain $\Omega _i$ after onedesign perturbation can be represented by

$$\begin{aligned} \tilde{\Omega }_i := \{\tilde{\varvec{x}} | \tilde{\varvec{x}} := \varvec{x} + \epsilon \varvec{v}(\varvec{x}), \,\varvec{x} \in \Omega _i\}. \end{aligned}$$

(47)

Here and hereafter, $\tilde{f}$ stands for the quantity f after onedesign perturbation. We now give some assumptions for $\varvec{v}$ on $\Gamma _{-1}$, $\Gamma _{d}$, $\Gamma _i^\mathrm{L}$ and $\Gamma _i^\mathrm{R}$. Regarding the infinite boundaries $\Gamma _{-1}$ and $\Gamma _{d-1}$, we may suppose that they are fixed, i.e.

$$\begin{aligned} \varvec{v}\equiv \varvec{0}\quad \text {on } \Gamma _{-1} \,\,\text {and}\,\, {\Gamma _{d-1}}. \end{aligned}$$

(48)

Meanwhile, we assume that the periodic boundaries $\Gamma _i^\mathrm{L}$ and $\Gamma _i^\mathrm{R}$ can move only in the $+y$-direction so that the period L is maintained, i.e.

$$\begin{aligned} \varvec{v}\cdot \varvec{n}=0\quad \text {on } \Gamma _i^\mathrm{L}\cup \Gamma _i^\mathrm{R}, \end{aligned}$$

(49)

where the normal vector $\varvec{n}$ is identical to $(-1,0)^T$ and $(1,0)^T$ on $\Gamma _i^\mathrm{L}$ and $\Gamma _i^\mathrm{R}$, respectively. Moreover, the periodicity of the geometry in the x-direction imposes that the onedesign velocity at $\varvec{x}\in \Gamma _i^\mathrm{L}$ and $\varvec{x}+L\varvec{e}_x\in \Gamma _i^\mathrm{R}$ must be same, i.e.

$$\begin{aligned} \varvec{v}(\varvec{x}+L\varvec{e}_x) = \varvec{v}(\varvec{x})\quad ~\mathrm{for}~ \varvec{x}\in \Gamma _i^\mathrm{L}. \end{aligned}$$

(50)

We will define the shape derivative of $\mathcal {L}$ from its variation due to the infinitesimal displacement $\epsilon \varvec{v}$. We suppose that $\tilde{u}$ is the solution of the primal problem with respect to $\tilde{\Omega }$. Then, the Lagrangian after onedesign perturbation can be obtained from Eq. 45 as

$$\begin{aligned} \mathcal {L}(\tilde{\Gamma }; \tilde{u}, \lambda ) = \mathcal {J}(\tilde{\Gamma };\tilde{u}) + \sum _{i=0}^{d-1} \mathop {\mathrm {Re}}\left[ \mathcal {G}(\tilde{\Omega }_i; \tilde{u}, \lambda ) \right] , \end{aligned}$$

(51)

where

$$\begin{aligned} \mathcal {G}(\tilde{\Omega }_i; \tilde{u}, \lambda ) = \int _{\tilde{\Omega }_i} \frac{1}{\varepsilon _i} \left( \tilde{\nabla } \lambda ^* \cdot \tilde{\nabla } \tilde{u} - k_i^2 \lambda ^* \tilde{u} \right) \mathrm {d}\tilde{\Omega } -\int _{\partial \tilde{\Omega }_i} \lambda ^* \tilde{q} \, \mathrm {d}\tilde{\Gamma } \end{aligned}$$

(52)

follows from Eq. 46. In Eq. 51, it is allowed to use the same $\lambda$ as that in Eq. 45 because the variation of $\mathcal {J}$, which we actually desire to evaluate, is irrelevant to the choice of $\lambda$ in Eq. 51 as long as $\tilde{u}$ satisfies Eq. 2. Subtracting Eq. 45 from Eq. 51, we can obtain the variation of $\mathcal {L}$, i.e.

$$\begin{aligned} \delta \mathcal {L}:= \mathcal {L}(\tilde{\Gamma }; \tilde{u}, \lambda ) - \mathcal {L}(\Gamma ; u, \lambda ) = \delta \mathcal {J}+ \sum _{i=0}^{d-1} \mathop {\mathrm {Re}}\left[ \delta \mathcal {G}_i \right] , \end{aligned}$$

(53)

where $\delta \mathcal {J}:=\mathcal {J}(\tilde{\Gamma };\tilde{u})-\mathcal {J}(\Gamma ;u)$ and $\delta \mathcal {G}_i:=\mathcal {G}(\tilde{\Omega }_i;\tilde{u},\lambda )-\mathcal {G}(\Omega _i;u,\lambda )$. Here, since $\delta \mathcal {J}\rightarrow 0$ as $\epsilon \rightarrow 0$, the variation $\delta \mathcal {J}$ can be expanded as the polynomial of $\epsilon$ without the term of $\epsilon ^0$, i.e.

$$\begin{aligned} \delta \mathcal {J}= \epsilon \mathcal {S}(\Gamma ) + \mathcal {O}(\epsilon ^2) \Longleftrightarrow \mathcal {J}(\tilde{\Gamma }) = \mathcal {J}(\Gamma ) + \epsilon \mathcal {S}(\Gamma ) + \mathcal {O}(\epsilon ^2), \end{aligned}$$

(54)

where $\mathcal {S}$ denotes the coefficient of $\epsilon$ and is defined as the (first-order) shape derivative of $\mathcal {J}$, i.e.

$$\begin{aligned} \mathcal {S}(\Gamma ):=\lim _{\epsilon \rightarrow 0}\frac{\delta \mathcal {J}}{\epsilon }=\lim _{\epsilon \rightarrow 0}\frac{\delta \mathcal {L}}{\epsilon }=\lim _{\epsilon \rightarrow 0}\frac{\delta \mathcal {J}+ \sum _{i=0}^{d-1} \mathop {\mathrm {Re}}\left[ \delta \mathcal {G}_i \right] }{\epsilon }. \end{aligned}$$

(55)

To express $\mathcal {S}$ explicitly, we need to expand $\mathcal {J}(\tilde{\Gamma };\tilde{u})$ and $\mathcal {G}(\tilde{\Omega }_i; \tilde{u}, \lambda )$ in the RHS of Eq. 51 with respect to the infinitesimal number $\epsilon$. First of all, $\tilde{u}$ and $\tilde{q}$ can be expanded in terms of $\epsilon$ as follows:

$$\begin{aligned} \tilde{u}(\tilde{\varvec{x}}) = u(\varvec{x}) + \epsilon \dot{u}(\varvec{x}) + \mathcal {O}( \epsilon ^2),\quad \tilde{q}(\tilde{\varvec{x}}) = q(\varvec{x}) + \epsilon \dot{q}(\varvec{x}) + \mathcal {O}( \epsilon ^2), \end{aligned}$$

(56)

where

$$\begin{aligned} \dot{u}(\varvec{x}):=\left. \frac{\mathrm {d}u(\tilde{\varvec{x}})}{\mathrm {d}\epsilon } \right| _{\epsilon =0},\quad \dot{q}(\varvec{x}):=\left. \frac{\mathrm {d}q(\tilde{\varvec{x}})}{\mathrm {d}\epsilon } \right| _{\epsilon =0}. \end{aligned}$$

Accordingly, applying Eq. 56 to $\tilde{u}_j(\tilde{\varvec{x}})=\tilde{u}_{j+1}(\tilde{\varvec{x}})$ and $\tilde{q}_j(\tilde{\varvec{x}})=-\tilde{q}_{j+1}(\tilde{\varvec{x}})$ obtained from the boundary conditions in Eq. 3 and then gathering the terms of $\varepsilon$, we can obtain the following boundary conditions in terms of $\dot{u}$ and $\dot{q}$:

$$\begin{aligned} \dot{u}_j = \dot{u}_{j+1}(=:\dot{u}),\quad \dot{q}_j = -\dot{q}_{j+1}(=:\dot{q})\quad \text { on }\Gamma _j. \end{aligned}$$

(57)

On the other hand, substituting Eq. 56 into Eq. 5 where $u(\varvec{x})$ is replaced with $\tilde{u}(\tilde{\varvec{x}})$ yields the following quasi-periodic condition for $\dot{u}$:

$$\begin{aligned} \partial _x^n \dot{u}(\varvec{x}+ L\varvec{e}_x) = e^{\mathrm {i}\beta L} \partial _x^n \dot{u}(\varvec{x}),\quad \partial _y^n \dot{u}(\varvec{x}+ L\varvec{e}_x) = e^{\mathrm {i}\beta L} \partial _y^n \dot{u}(\varvec{x}). \end{aligned}$$

(58)

Next, we will write down $\delta \mathcal {G}_i$, using the quantities without tilde. To this end, we utilise the following formulae [83]:

$$\begin{aligned}&\mathrm {d}\tilde{\Omega } = \det \varvec{F}\,\mathrm {d}\Omega = (1 + \epsilon \nabla \cdot \varvec{v})\,\mathrm {d}\Omega + \mathcal {O}(\epsilon ^2), \end{aligned}$$

(59a)

$$\begin{aligned}&\mathrm {d}\tilde{\Gamma } = \det \varvec{F}\left| \varvec{F}^{-T} \varvec{n}\right| \,\mathrm {d}\Gamma = \left[ 1 + \epsilon \bigl ( \nabla \cdot \varvec{v} - \varvec{n} \cdot (\nabla \varvec{v}) \varvec{n} \bigr )\right] \, \mathrm {d}\Gamma + \mathcal {O}(\epsilon ^2), \end{aligned}$$

(59b)

where $\varvec{F}\in \mathbb {R}^{2\times 2}$ represents the deformation gradient defined by $F_{ij}:=\frac{\partial \tilde{x}_i}{\partial x_j}$ and $\varvec{n}$ denotes the outward normal vector of $\partial \Omega _i$ in this context. Now, by substituting Eq. 56 and Eq. 59a into Eq. 52, we can express $\delta \mathcal {G}_i$ as follows:

$$\begin{aligned} \delta \mathcal {G}_i= & {} \epsilon \, \Bigg [ \int _{\partial \Omega _i} \left( \varphi ^*\dot{u} - \lambda ^* \dot{q}\right) \, \mathrm {d}\Gamma \nonumber \\&- \int _{\Omega _i} \frac{1}{\varepsilon _i} \left( {{{\Delta }}}\lambda ^* + k_i^2 \lambda ^* \right) \dot{u} \, \mathrm {d}\Omega \nonumber \\&+ \int _{\partial \Omega _i} \frac{1}{\varepsilon _i} \left( \nabla \lambda ^* \cdot \nabla u - k_i^2 \lambda ^* u \right) \varvec{v} \cdot \varvec{n} \, \mathrm {d}\Gamma \nonumber \\&- \int _{\partial \Omega _i} \left( q (\varvec{v} \cdot \nabla \lambda ^*) + \varphi ^* (\varvec{v} \cdot \nabla u) \right) \, \mathrm {d}\Gamma \nonumber \\&- \int _{\partial \Omega _i} \lambda ^* q \left( \nabla \cdot \varvec{v} - \varvec{n} \cdot (\nabla \varvec{v}) \varvec{n} \right) \, \mathrm {d}\Gamma \nonumber \\&+ \int _{\Omega _i} \frac{1}{\varepsilon _i} \left( {{{\Delta }}}\lambda ^* + k_i^2 \lambda ^*\right) (\varvec{v} \cdot \nabla u) \, \mathrm {d}\Omega \Bigg ] + \mathcal {O}(\epsilon ^2), \end{aligned}$$

(60)

where the Gauss’s divergence theorem was used and $\varphi ^*$ represents the flux of $\lambda ^*$, i.e. $\varphi ^*:=\frac{1}{\varepsilon _i}\frac{\partial \lambda ^*}{\partial n_i}$; note that the relative permittivity $\varepsilon _i$ is not conjugated.

Meanwhile, the expression of $\delta \mathcal {J}$ (thus, $\mathcal {S}$) depends on the choice of $\mathcal {J}$. In this paper, we consider two cases of $\mathcal {J}$ in the subsequent sections.

1.1 Case I: magnetic field strength

Let us consider the objective function $\mathcal {J}^\mathrm{I}$ in Eq. 27. From Eq. 56, we can obtain

$$\begin{aligned} \frac{\left| \tilde{u}(\tilde{\varvec{x}})\right| ^2}{2} - \frac{\left| u(\varvec{x})\right| ^2}{2}= & {} \frac{(u^*(\varvec{x}) + \epsilon \dot{u}^*(\varvec{x}))(u(\varvec{x}) + \epsilon \dot{u}(\varvec{x}))}{2} - \frac{u^*(\varvec{x}) u(\varvec{x})}{2} + \mathcal {O}(\epsilon ^2)\\= & {} \epsilon \mathop {\mathrm {Re}}[ u^*(\varvec{x}) \dot{u}(\varvec{x})] + \mathcal {O}(\epsilon ^2). \end{aligned}$$

Thus, we can compute the variation of $\mathcal {J}^\mathrm{I}$ as follows:

$$\begin{aligned} \delta \mathcal {J}^\mathrm{I}:= & {} \mathcal {J}^\mathrm{I}(\tilde{\Gamma }; \tilde{u}) - \mathcal {J}^\mathrm{I}(\Gamma ; u) \nonumber \\= & {} \epsilon \sum _{i=0}^{d-1} \sum _{m=0}^{M_i-1} \mathop {\mathrm {Re}}\left[ \int _{\Omega _i} u^*(\varvec{x}) \dot{u}(\varvec{x}) \updelta (\varvec{x} - \varvec{z}_m^{(i)}) \, \mathrm {d}\Omega \right] + \mathcal {O}(\epsilon ^2), \end{aligned}$$

(61)

where we assume that each observation point belongs to the same domain before and after deformation.

By substituting Eqs. 60 and 61 into 55, the shape derivative of $\mathcal {J}^\mathrm{I}$, denoted by $\mathcal {S}^\mathrm{I}$, can be written as follows:

$$\begin{aligned} \mathcal {S}^\mathrm{I}(\Gamma )= & {} \mathop {\mathrm {Re}}\sum _{i=0}^{d-1} \Bigg [ \int _{\partial \Omega _i} \left( \varphi ^* \dot{u} - \lambda ^* \dot{q} \right) \, \mathrm {d}\Gamma \nonumber \\&+ \int _{\partial \Omega _i} \frac{1}{\varepsilon _i}\left( \nabla \lambda ^* \cdot \nabla u - k_i^2 \lambda ^* u\right) \varvec{v} \cdot \varvec{n} \, \mathrm {d}\Gamma \nonumber \\&- \int _{\partial \Omega _i} \frac{1}{\varepsilon _i} \left( q (\varvec{v} \cdot \nabla \lambda ^*) + \varphi ^* (\varvec{v} \cdot \nabla u) \right) \, \mathrm {d}\Gamma \nonumber \\&- \int _{\partial \Omega _i} \lambda ^* q \left( \nabla \cdot \varvec{v} - \varvec{n} \cdot (\nabla \varvec{v} \varvec{n}) \right) \, \mathrm {d}\Gamma \nonumber \\&+ \int _{\Omega _i} \frac{1}{\varepsilon _i} \left( {{{\Delta }}}\lambda ^* + k_i^2 \lambda ^*\right) (\varvec{v} \cdot \nabla u) \, \mathrm {d}\Omega \nonumber \\&- \int _{\Omega _i} \frac{1}{\varepsilon _i} \left( {{{\Delta }}}\lambda ^* + k_i^2 \lambda ^* - \varepsilon _i \sum _{m=0}^{M_i-1} u^* \updelta (\varvec{x} - \varvec{z}_m^{(i)}) \right) \dot{u} \, \mathrm {d}\Omega \Bigg ]. \end{aligned}$$

(62)

Now, we will re-express the boundary integrals in terms of $\partial \Omega _i$ in Eq. 62 with those in terms of $\Gamma _j$. To this end, we will show that $\Gamma _i^\mathrm{R}$, $\Gamma _i^\mathrm{L}$, $\Gamma _{-1}$ and $\Gamma _{d-1}$ do not actually contribute to Eq. 62 by giving the quasi-periodic condition, i.e.

$$\begin{aligned} \partial _x^n \lambda ^*(\varvec{x}+ L\varvec{e}_x) = e^{-\mathrm {i}\beta L} \partial _x^n \lambda ^*(\varvec{x}),\quad \partial _y^n \lambda ^*(\varvec{x}+ L\varvec{e}_x) = e^{-\mathrm {i}\beta L} \partial _y^n \lambda ^*(\varvec{x}), \end{aligned}$$

(63)

and the radiation condition, i.e.

$$\begin{aligned} \lambda ^*(x,y)= {\left\{ \begin{array}{ll} \sum _n C_n e^{\mathrm {i}\phi _n x} e^{\mathrm {i}\psi _n (y-y_0)} &{} ~\mathrm{for}~ y>y_0,\\ \sum _n D_n e^{\mathrm {i}\phi _n x} e^{\mathrm {i}\psi _n (-y-y_0)} &{} ~\mathrm{for}~ y<-y_0, \end{array}\right. } \end{aligned}$$

(64)

where $C_n$ and $D_n$ are complex numbers, $\phi _n:=-\beta +\frac{2\pi n}{L}$ and $\psi _n:=\sqrt{k^2-\phi _n^2}$. Similarly to the radiation condition for $u^\mathrm{sca}$ in Eq. 6, we consider $n\in \mathbb {Z}$ such as $k^2\ge \phi _n^2$. Note that, unlikely to $\beta$ for u, the phase difference is $-\beta$ for $\lambda ^*$ in both Eqs. 63 and 64. On the other hand, the quasi-periodic conditions of u in Eq. 5, $\dot{u}$ in Eq. 58 and $\lambda ^*$ in Eq. 63 yield those of q, $\dot{q}$ and $\varphi ^*$, respectively, i.e.

$$\begin{aligned}&q(\varvec{x}+L\varvec{e}_x)=-e^{\mathrm {i}\beta L}q(\varvec{x})\quad ~\mathrm{for}~ \varvec{x}\in \Gamma _i^\mathrm{L},\end{aligned}$$

(65a)

$$\begin{aligned}&\dot{q}(\varvec{x}+L\varvec{e}_x)=-e^{\mathrm {i}\beta L}\dot{q}(\varvec{x})\quad ~\mathrm{for}~ \varvec{x}\in \Gamma _i^\mathrm{L},\end{aligned}$$

(65b)

$$\begin{aligned}&\varphi ^*(\varvec{x}+L\varvec{e}_x)=-e^{-\mathrm {i}\beta L}\varphi ^*(\varvec{x})\quad ~\mathrm{for}~ \varvec{x}\in \Gamma _i^\mathrm{L}, \end{aligned}$$

(65c)

where the negative sign in the RHSs contemplates the fact that the normal vector at $\varvec{x}\in \Gamma _i^\mathrm{L}$ is opposite to that at the corresponding point $\varvec{x}+L\varvec{e}_x\in \Gamma _i^\mathrm{R}$. First, from Eqs. 58, 63, 65b and 65c, the integral over $\Gamma _i^\mathrm{R}$ is cancelled with that over $\Gamma _i^\mathrm{L}$ in the first boundary term of Eq. 62. Second, the integrals over $\Gamma _i^\mathrm{L}$ and $\Gamma _i^\mathrm{R}$ in the second boundary term vanish because of Eq. 49. Third, the integrals over $\Gamma _i^\mathrm{L}$ and $\Gamma _i^\mathrm{R}$ in the third boundary term are cancelled because of the identity of $\varvec{v}$ in Eq. 50 and the quasi-periodic conditions of $\nabla u$ from Eq. 5, $\nabla \lambda ^*$ from Eq. 63, q in Eq. 65a and $\dot{q}$ in Eq. 65b. Fourth, the integrals over $\Gamma _i^\mathrm{L}$ and $\Gamma _i^\mathrm{R}$ in the fourth boundary term are cancelled.¹⁶ Fifth, because of the radiation conditions for u in Eq. 6 and $\lambda ^*$ in Eq. 64, the integrals over the infinite boundaries $\Gamma _{-1}$ and $\Gamma _{d-1}$ in the first boundary term vanish (see the proof in 8.1). Sixth and last, the integrals over $\Gamma _{-1}$ and $\Gamma _{d-1}$ in the second, third and fourth boundary terms vanish because of the assumption in Eq. 48. Therefore, we can regard $\partial \Omega _i$ in Eq. 62 as follows:

$$\begin{aligned} \partial \Omega _i={\left\{ \begin{array}{ll} \Gamma _0 &{} ~\mathrm{for}~ i=0,\\ \Gamma _{d-2} &{} ~\mathrm{for}~ i=d-1,\\ \Gamma _{i-1}\cup \Gamma _i &{} \mathrm{otherwise}. \end{array}\right. } \end{aligned}$$

Hence, by considering the both sides of each boundary, we can rewrite Eq. 62 in the following boundary-centric expression:

$$\begin{aligned}\mathcal {S}^\mathrm{I}(\Gamma )= & {} \mathop {\mathrm {Re}}\sum _{j=0}^{d-2} \Bigg [ \int _{\Gamma _j} \left\{ \left( \varphi _j^* + \varphi _{j+1}^* \right) \dot{u} - \left( \lambda _j^* - \lambda _{j+1}^* \right) \dot{q} \right\} \, \mathrm {d}\Gamma \nonumber \\&+ \int _{\Gamma _j} \left\{ \left( \frac{1}{\varepsilon _j} \nabla \lambda _j^* \cdot \nabla u_j - \frac{1}{\varepsilon _{j+1}} \nabla \lambda _{j+1}^* \cdot \nabla u_{j+1} \right) - \left( \frac{k_j^2}{\varepsilon _j} \lambda _j^* - \frac{k_{j+1}^2}{\varepsilon _{j+1}} \lambda _{j+1}^* \right) u \right\} \varvec{v} \cdot \varvec{n} \, \mathrm {d}\Gamma \nonumber \\&- \int _{\Gamma _j} \left\{ q \varvec{v} \cdot ( \nabla \lambda _j^* - \nabla \lambda _{j+1}^* ) + \left( \varphi _j^* (\varvec{v} \cdot \nabla u_j) + \varphi _{j+1}^* (\varvec{v} \cdot \nabla u_{j+1}) \right) \right\} \mathrm {d}\Gamma \nonumber \\&- \int _{\Gamma _j} (\lambda _j^* - \lambda _{j+1}^* ) q \left( \nabla \cdot \varvec{v} - \varvec{n} \cdot (\nabla \varvec{v} \varvec{n}) \right) \, \mathrm {d}\Gamma \Bigg ] \nonumber \\&+\mathop {\mathrm {Re}}\sum _{i=0}^{d-1} \Bigg [\int _{\Omega _i} \frac{1}{\varepsilon _i} ({{{\Delta }}}\lambda ^* + k_i^2 \lambda ^*) ( \varvec{v} \cdot \nabla u) \, \mathrm {d}\Omega \nonumber \\&- \int _{\Omega _i} \frac{1}{\varepsilon _i} \left( {{{\Delta }}}\lambda ^* + k_i^2 \lambda ^* - \varepsilon _i \sum _{m=0}^{M_i-1} u^* \updelta (\varvec{x} - \varvec{z}_m^{(i)}) \right) \dot{u} \, \mathrm {d}\Omega \Bigg ] \end{aligned}$$

(66)

where the boundary conditions Eqs. 3 and 57 were considered and thus the symbols u, q, $\dot{u}$ and $\dot{q}$ (without subscript) have appeared on $\Gamma _j$. Regarding the symbol $\varvec{n}$ on $\Gamma _j$, recall Fig. 3 (right). Note that we define $\nabla u_j(\varvec{x}):=\displaystyle \lim _{\varvec{x'}\in \Omega _{j}\rightarrow \varvec{x}}\nabla u(\varvec{x'})$ and $\nabla u_{j+1}(\varvec{x}):=\displaystyle \lim _{\varvec{x'}\in \Omega _{j+1}\rightarrow \varvec{x}}\nabla u(\varvec{x'})$ for $\varvec{x}\in \Gamma _j$ and these gradients are different from each other in general. Similarly, $\nabla \lambda _j$, $\nabla \lambda _{j+1}$, $\varphi _j^*$ and $\varphi _{j+1}^*$ are defined as the limits to $\Gamma _j$ from $\Omega _j$ or $\Omega _{j+1}$.

In accordance with the concept of the adjoint variable method, we give some constraints to $\lambda ^*$ so that we can avoid calculating the derivatives with respect to $\epsilon$, i.e. $\dot{u}$ and $\dot{q}$, in the first and last terms of Eq. 66. First, we request that $\lambda$ should be the solution of the following inhomogeneous Helmholtz equation:

$$\begin{aligned} {{{\Delta }}}\lambda ^*(\varvec{x}) + k_i^2 \lambda ^*(\varvec{x}) = \varepsilon _i \sum _{m=0}^{M_i-1} u^*(\varvec{x}) \updelta (\varvec{x} - \varvec{z}_m^{(i)})\quad ~\mathrm{for}~ \varvec{x}\in \Omega _i. \end{aligned}$$

(67)

Second, we request that $\lambda$ should satisfy the following boundary conditions:

$$\begin{aligned} \lambda _j^* = \lambda _{j+1}^*(=:\lambda ^*),\quad \varphi _j^*=-\varphi _{j+1}^*(=:\varphi ^*)\quad \text {on }\Gamma _j. \end{aligned}$$

(68)

Then, by substituting Eq. 67 into the two domain integrals and Eq. 68 into the first and fourth boundary integrals, we can write Eq. 66 as follows:

$$\begin{aligned} \mathcal {S}^\mathrm{I}(\Gamma )= & {} \mathop {\mathrm {Re}}\sum _{j=0}^{d-2} \Bigg [ \int _{\Gamma _j} \left( \frac{1}{\varepsilon _j} \nabla \lambda _j^* \cdot \nabla u_j - \frac{1}{\varepsilon _{j+1}} \nabla \lambda _{j+1}^* \cdot \nabla u_{j+1} \right) \varvec{v} \cdot \varvec{n} \, \mathrm {d}\Gamma \nonumber \\&- \int _{\Gamma _j} q \varvec{v} \cdot ( \nabla \lambda _j^* - \nabla \lambda _{j+1}^* ) \, \mathrm {d}\Gamma - \int _{\Gamma _j} \varphi ^* \varvec{v} \cdot (\nabla u_j - \nabla u_{j+1}) \, \mathrm {d}\Gamma \Bigg ], \end{aligned}$$

(69)

where the first domain integral in Eq. 66, which is reduced to $u^*(\varvec{z}_m^{(i)})\varvec{v}(\varvec{z}_m^{(i)})\cdot \nabla u(\varvec{z}_m^{(i)})$ due to Eq. 67 and the property of the Dirac’s delta function, vanished because of the assumption that the observation points are fixed, i.e. $\varvec{v}(\varvec{z}_m^{(i)})\equiv \varvec{0}$. In addition, $\frac{k_j^2}{\varepsilon _j}-\frac{k_{j+1}^2}{\varepsilon _{j+1}}\equiv 0$ was used.

We can show that any deformation in tangential direction does not contribute to the shape derivative. Namely, if we set $\varvec{v}=\sigma \varvec{s}$, where $\varvec{s}$ is the unit tangential vector to $\Gamma _j$ and $\sigma$ is any real number, we can derive $\mathcal {S}^\mathrm{I}\equiv 0$. As a matter of fact, since $\varvec{v}\cdot \varvec{n}=0$ holds, the first integral in Eq. 69 vanishes. In addition, the boundary conditions $u_j=u_{j+1}$ in Eq. 3 and $\lambda _j^*=\lambda _{j+1}^*$ in Eq. 68 yield the continuities of the tangential derivatives, i.e. $\varvec{s}\cdot \nabla u_j=\varvec{s}\cdot \nabla u_{j+1}$ and $\varvec{s}\cdot \nabla \lambda _j^*=\varvec{s}\cdot \nabla \lambda _{j+1}^*$ on $\Gamma _j$, respectively, from which the second and third integrals vanish. Consequently, we can write $\varvec{v}$ as $(\varvec{v}\cdot \varvec{n})\varvec{n}$ in Eq. 69.

Now, we represent $\nabla u_j$ and $\nabla \lambda _j^*$ in Eq. 69 with their normal and tangential components, i.e.

$$\begin{aligned} \nabla u_j = \frac{\partial u_j}{\partial n} \varvec{n} + \frac{\partial u_j}{\partial s} \varvec{s},\quad \nabla \lambda _j^* = \frac{\partial \lambda _j^*}{\partial n} \varvec{n} + \frac{\partial \lambda _j^*}{\partial s} \varvec{s}. \end{aligned}$$

Then, the first term in Eq. 69 can be reduced as:

$$\begin{aligned} \int _{\Gamma _j} \left( (\varepsilon _{j} - \varepsilon _{j+1}) \varphi ^* q + \left( \frac{1}{\varepsilon _{j}} - \frac{1}{\varepsilon _{j+1}}\right) \frac{\partial \lambda ^*}{\partial s} \frac{\partial u}{\partial s} \right) \varvec{v} \cdot \varvec{n} \, \mathrm {d}\Gamma , \end{aligned}$$

while the second and third integrals can be expressed as:

$$\begin{aligned}&- \int _{\Gamma _j} q (\varvec{v}\cdot \varvec{n})\varvec{n}\cdot \left( \frac{\partial \lambda _j^*}{\partial n} \varvec{n} + \frac{\partial \lambda _j^*}{\partial s} \varvec{s} - \frac{\partial \lambda _{j+1}^*}{\partial n} \varvec{n} - \frac{\partial \lambda _{j+1}^*}{\partial s} \varvec{s} \right) \, \mathrm {d}\Gamma = - \int _{\Gamma _j} \left( \varepsilon _j - \varepsilon _{j+1} \right) \varphi ^* q \varvec{v}\cdot \varvec{n}\, \mathrm {d}\Gamma ,\\&- \int _{\Gamma _j} \varphi ^* (\varvec{v}\cdot \varvec{n})\varvec{n}\cdot \left( \frac{\partial u_j}{\partial n} \varvec{n} + \frac{\partial u_j}{\partial s} \varvec{s} - \frac{\partial u_{j+1}}{\partial n} \varvec{n} - \frac{\partial u_{j+1}}{\partial s} \varvec{s} \right) \, \mathrm {d}\Gamma = - \int _{\Gamma _j} \left( \varepsilon _j - \varepsilon _{j+1} \right) \varphi ^* q \varvec{v}\cdot \varvec{n}\, \mathrm {d}\Gamma . \end{aligned}$$

In these expressions, the boundary conditions in Eqs. 3 and 68 were used. Consequently, we can obtain the shape derivative $\mathcal {S}^\mathrm{I}$ in Eq. 28.

1.2 Appendix C.2: Case II: energy absorption rate

We consider the objective function $\mathcal {J}^\mathrm{II}$ in Eq. 30. Using Eqs. 56 and 59b, we can obtain the variation of $\mathcal {J}^\mathrm{II}$ as follows:

$$\begin{aligned} \delta \mathcal {J}^\mathrm{II}:= & {} \mathcal {J}^\mathrm{II}(\tilde{\Gamma };\tilde{u},\tilde{q})-\mathcal {J}^\mathrm{II}(\Gamma ;u,q)\nonumber \\= & {} \kappa \epsilon \mathop {\mathrm {Re}}\Bigg [ \int _{\partial \Omega _\alpha } \mathrm {i}(\dot{u} q^* - u^* \dot{q}) \, \mathrm {d}\Gamma - \int _{\partial \Omega _\alpha } \mathrm {i}u^* q \big (\nabla \cdot \varvec{v} - \varvec{n} \cdot (\nabla \varvec{v} \varvec{n}) \big ) \, \mathrm {d}\Gamma \Bigg ] + \mathcal {O}(\epsilon ^2). \end{aligned}$$

(70)

Successively, we can compute the shape derivative, denoted by $\mathcal {S}^\mathrm{II}$, by substituting $\delta \mathcal {J}^\mathrm{II}$ in Eq. 70 and $\delta \mathcal {G}_i$ in Eq. 60 into the definition of the shape derivative in Eq. 55. The way of the derivation is similar to that of the previous case. In this case, it is requested that the (conjugated) Lagrange multiplier $\lambda ^*$ should be the solution of the following adjoint problem:

$$\begin{aligned}&{{{\Delta }}}\lambda ^* + k_i^2 \lambda ^* = 0 \quad \text {in} \Omega _i, \end{aligned}$$

(71a)

$$\begin{aligned}&\lambda _{\alpha -1}^* = \lambda _\alpha ^* + \mathrm {i}\kappa u_\alpha ^* \, (=: \lambda ^*), \quad \varphi _{\alpha -1}^* = - \varphi _\alpha ^* - \mathrm {i}\kappa q_\alpha ^* \, (=: \varphi ^*) \quad \text {on } \Gamma _{\alpha -1}, \end{aligned}$$

(71b)

$$\begin{aligned}&\lambda _\alpha ^* = \lambda _{\alpha +1}^* - \mathrm {i}\kappa u_\alpha ^*\, (=: \lambda ^*), \quad \varphi _\alpha ^* = - \varphi _{\alpha +1}^* - \mathrm {i}\kappa q_\alpha ^*\, (=: \varphi ^*) \quad \text {on } \Gamma _{\alpha }, \end{aligned}$$

(71c)

$$\begin{aligned}&\lambda _j^* = \lambda _{j+1}^* \, (=: \lambda ^*),\quad \varphi _j^* = - \varphi _{j+1}^* \, (=: \varphi ^*) \quad \text {on } \Gamma _{j} for j\ne \alpha -1,\alpha , \end{aligned}$$

(71d)

$$\begin{aligned}&\lambda ^*(\varvec{x} + L \varvec{e}_x) = e^{-\mathrm {i}\beta L} \lambda ^*(\varvec{x}),\quad \varphi ^*(\varvec{x} + L \varvec{e}_x) = e^{-\mathrm {i}\beta L} \varphi ^*(\varvec{x})\quad \text {on }\Gamma _j^\mathrm{L}, \end{aligned}$$

(71e)

$$\begin{aligned}&\lambda ^*(x,y)= {\left\{ \begin{array}{ll} \sum _n C_n e^{\mathrm {i}\phi _n x} e^{\mathrm {i}\psi _n (y-y_0)} &{} ~\mathrm{for}~ y>y_0,\\ \sum _n D_n e^{\mathrm {i}\phi _n x} e^{\mathrm {i}\psi _n (-y-y_0)} &{} ~\mathrm{for}~ y<-y_0, \end{array}\right. } \end{aligned}$$

(71f)

where $u_\alpha$ and $q_\alpha$ [which are defined in Eq. 3] in the boundary conditions Eqs. 71b and 71c are the solutions obtained from the primal problem. Finally, we can obtain $\mathcal {S}^\mathrm{II}$ in Eq. 31.

Appendix D: Evaluation of the integral of $\varphi ^*\dot{u}-\lambda ^*\dot{q}$ over the infinite boundaries $\Gamma _{-1}$ and $\Gamma _{d-1}$ in Eq. (62)

Let us consider the case of $y>y_0$. Similarly to Eq. 56, the scattered field $\tilde{u}^\mathrm{sca}$ for the deformed domain $\tilde{\Omega }_0$ can be expanded in terms of $\epsilon$ as follows:

$$\begin{aligned} \tilde{u}^\mathrm{sca}(\tilde{\varvec{x}}) = u^\mathrm{sca}(\varvec{x}) + \epsilon \dot{u}^\mathrm{sca}(\varvec{x}) + \mathcal {O}( \epsilon ^2). \end{aligned}$$

Here, $\dot{u}^\mathrm{sca}(\varvec{x}):=\left. \frac{\mathrm {d}u^\mathrm{sca}(\tilde{\varvec{x}})}{\mathrm {d}\epsilon }\right| _{\epsilon =0}$ can be rewritten from the radiation condition of $u^\mathrm{sca}$ in Eq. 6 as follows:

$$\begin{aligned} \dot{u}^\mathrm{sca}(x,y)=\left. \frac{\mathrm {d}}{\mathrm {d}\epsilon }\sum _m A_m e^{\mathrm {i}\xi _m (x+\epsilon v_x)} e^{\mathrm {i}\eta _m (y+\epsilon v_y-y_0)}\right| _{\epsilon =0} =\sum _m A_m \mathrm {i}(\xi _m v_x + \eta _m v_y) e^{\mathrm {i}\xi _m x}e^{\mathrm {i}\eta _m (y-y_0)}, \end{aligned}$$

where $\epsilon \varvec{v}=(\epsilon v_x,\epsilon v_y)^T$ is the infinitesimal vector defined in Eq. 47. By comparing with the radiation condition of $u^\mathrm{sca}$ in Eq. 6, this equation means that $\dot{u}^\mathrm{sca}$ also satisfies the radiation condition. Here, $\dot{u}^\mathrm{sca}$ is equivalent to $\dot{u}(\equiv \dot{u}^\mathrm{inc}+\dot{u}^\mathrm{sca})$ because the incident field is irrelevant to the deformation, i.e. $\tilde{u}^\mathrm{inc}(\tilde{\varvec{x}})\equiv u^\mathrm{inc}(\varvec{x})\Leftrightarrow \dot{u}^\mathrm{inc}\equiv 0$. Therefore, $\dot{u}$ satisfies the following radiation condition:

$$\begin{aligned} \dot{u}(x,y)=\sum _m E_m e^{\mathrm {i}\xi _m x}e^{\mathrm {i}\eta _m (y-y_0)}, \end{aligned}$$

(72)

where $E_m:=A_m \mathrm {i}(\xi _m v_x + \eta _m v_y)$.

With using the radiation condition for $\dot{u}$ in Eq. 72 as well as $\varphi ^*$ in Eq. 64, we can write down the integrand of the first boundary term in Eq. 62 regarding the upper infinite boundary $\Gamma _{-1}$, where the outward normal is $(0,1)^T$, as follows:

$$\begin{aligned} \varphi ^*\dot{u}-\lambda ^*\dot{q} =\frac{1}{\varepsilon _0}\frac{\partial \lambda ^*}{\partial y}\dot{u}-\lambda ^*\frac{1}{\varepsilon _0}\frac{\partial \dot{u}}{\partial y}=\frac{1}{\varepsilon _0}\sum _n \sum _m C_n E_m e^{\mathrm {i}\frac{2\pi (m+n)}{L}x} e^{\mathrm {i}(\psi _n + \eta _m)(y-y_0)} \mathrm {i}(\psi _n-\eta _m). \end{aligned}$$

Then, the terms such as $n+m\ne 0$ vanish when they are integrated over $\Gamma _{-1}$ (i.e. $x\in [0,L]$ for any $y>y_0$). Meanwhile, the terms such as $n+m=0$ also vanish because $n+m=0$ is followed by $\phi _n=\phi _{-m}=-\beta -\frac{2\pi m}{L}=-\xi _m$ and thus $\psi _n=\sqrt{k^2-(-\xi _m)^2}=\eta _m$. Therefore, the integral of $\varphi ^*\dot{u}-\lambda ^*\dot{q}$ over $\Gamma _{-1}$ vanishes and does not contribute to $\mathcal {S}^\mathrm{I}$ in Eq. 62.

Similarly, the integral of $\varphi ^*\dot{u}-\lambda ^*\dot{q}$ over the lower infinite boundary $\Gamma _{d-1}$ results in zero by starting with $y<-y_0$ instead of $y>y_0$.

Appendix E: Energy absorption rate $\mathcal {A}$

The time-averaged energy per unit time that absorbed in a closed domain $\Omega _\alpha$ (where $\alpha \ne 0$ and $n-1$) is defined by

$$\begin{aligned} e_\mathrm{absorp} := - \int _{\Omega _\alpha } \mathrm {div} \langle {\varvec{S}}\rangle \, \mathrm {d}\Omega = - \int _{\partial \Omega _\alpha } \langle {\varvec{S}}\rangle \cdot \varvec{n} \, \mathrm {d}\Gamma , \end{aligned}$$

(73)

where the Gauss’s divergence theorem was used, $\varvec{n}$ denotes the unit outward normal to the boundary of $\Omega _\alpha$, $\varvec{S}(\varvec{x}, t):=\varvec{E}(\varvec{x}, t)\times \varvec{H}(\varvec{x}, t)$ denotes the (real-valued) Poynting vector, and $\langle {\varvec{S}}\rangle$ denotes the time-average of $\varvec{S}$ over the period $T:=\frac{2\pi }{\omega }$, i.e. $\langle {\varvec{S}(\varvec{x})}\rangle := \frac{1}{T} \int _0^T \varvec{S}(\varvec{x},t) \mathrm {d}t$. Since we suppose the TM-mode, where the magnetic field $\varvec{H}$ is given by Eq. 1, the corresponding electric field $\varvec{E}$ in $\Omega _\alpha$ is represented as:

$$\begin{aligned} \varvec{E}(\varvec{x},t)=\left( \mathop {\mathrm {Re}}\left[ \frac{-1}{\mathrm {i}\omega \upepsilon _0\varepsilon _\alpha }\frac{\partial u(\varvec{x})}{\partial x_2}e^{-\mathrm {i}\omega t}\right] ,\mathop {\mathrm {Re}}\left[ \frac{1}{\mathrm {i}\omega \upepsilon _0\varepsilon _\alpha }\frac{\partial u(\varvec{x})}{\partial x_1}e^{-\mathrm {i}\omega t}\right] ,0\right) ^T \equiv \mathop {\mathrm {Re}}\left[ \hat{\varvec{E}}(\varvec{x})e^{-\mathrm {i}\omega t}\right] , \end{aligned}$$

where $\upepsilon _0$ denotes the permittivity in vacuum ($\approx$air)¹⁷ and $\hat{\varvec{E}}(\varvec{x}):=\left( \frac{-1}{\mathrm {i}\omega \upepsilon _0\varepsilon }\frac{\partial u(\varvec{x})}{\partial x_2},\frac{1}{\mathrm {i}\omega \upepsilon _0\epsilon }\frac{\partial u(\varvec{x})}{\partial x_1},0\right) ^T\in \mathbb {C}^{3}$. Similarly to $\hat{\varvec{E}}$, we define $\hat{\varvec{H}}(\varvec{x}):=\left( 0,0,u(\varvec{x})\right) ^T\in \mathbb {C}^{3}$ so that $\varvec{H}(\varvec{x},t)\equiv \mathop {\mathrm {Re}}[\hat{\varvec{H}}(\varvec{x})e^{-\mathrm {i}\omega t}]$ holds. Subsequently, we obtain $\langle {\varvec{S}}\rangle = \frac{1}{2}\left( \mathop {\mathrm {Re}}[\hat{\varvec{E}}]\times \mathop {\mathrm {Re}}[\hat{\varvec{H}}]+\mathfrak {I}[\hat{\varvec{E}}]\times \mathfrak {I}[\hat{\varvec{H}}]\right)$. By substituting this into Eq. 73 and doing some calculus by means of the Maxwell equations, $e_\mathrm{absorp}$ can be expressed as follows:

$$\begin{aligned} e_\mathrm{absorp} = - \frac{\mathrm {i}}{4 \omega \upepsilon _0} \int _{\partial \Omega _\alpha } (u q^* - u^* q) \mathrm {d}\Gamma , \end{aligned}$$

(74)

where $q:=\frac{1}{\varepsilon _\alpha }\frac{\partial u}{\partial n}$, which was defined in Eq. 7. Note that, since u satisfies the quasi-periodicity in Eq. 5 and, thus, $uq^*-u^*q(\equiv u\frac{\mathrm {d}u^*}{\mathrm {d}x}-u^*\frac{\mathrm {d}u}{\mathrm {d}x})$ on $\Gamma _\alpha ^\mathrm{R}$ is identical to that on $\Gamma _\alpha ^\mathrm{L}$ except the sign, the integral over $\Gamma _\alpha ^\mathrm{L}\cup \Gamma _\alpha ^\mathrm{R}$ vanishes. Therefore, the integral over $\partial \Omega _\alpha$ is equivalent to that over $\Gamma _{\alpha -1}\cup \Gamma _\alpha$ in net.

To normalise $e_\mathrm{absorp}$ in Eq. 74, we consider the time-averaged energy of the incident planewave $u^\mathrm{inc}$ (see Eq. 4) through the effective length of $L\sin \theta$ per unit time, i.e.

$$\begin{aligned} e_\mathrm{inc}:=\frac{1}{2}\sqrt{\frac{{\mu }_0}{\upepsilon _0}}\left| a^\mathrm{inc}\right| ^2 L\sin \theta , \end{aligned}$$

where ${\mu }_0$ denotes the permeability in vacuum.

Consequently, we define the (normalised) energy absorption rate $\mathcal {A}$ as:

$$\begin{aligned} \mathcal {A}:=\frac{e_\mathrm{absorp}}{e_\mathrm{inc}}=\frac{\mathrm {i}\kappa }{2} \int _{\partial \Omega _\alpha } (u q^* - u^* q)\mathrm {d}\Gamma = \kappa \mathfrak {I}\left[ \int _{\partial \Omega _\alpha } u q^* \mathrm {d}\Gamma \right] \end{aligned}$$

(75)

where $\kappa :=-\frac{1}{kL\sin \theta \left| a^\mathrm{inc}\right| ^2}$. Note that $k:=\frac{\omega }{c}$ denotes the wavenumber in vacuum, which corresponds to $k_0$, i.e. the wavenumber in $\Omega _0$ (to which the incident wave is given), in the main text. Note also that $0\le \mathcal {A}\le 1$ holds if $\mathfrak {I}\varepsilon \ge 0$ because of the conservation of the energy.

Appendix F: Lists of optimal control points

Table 1 shows the control points that can maximise the objective function $\mathcal {A}$ in the first and second tries of the demonstration described in Subsect. 4.3.

Table 1 Design variables (i.e. coordinates of the sixteen control points) with side constraints $(x_\nu ,y_\nu )^T(\equiv \varvec{p}_\nu ^{(1)})$ obtained by the first try (left) and second try (right)