Mouhamadou Baïdy Dia^1*, Ibrahima Sagno¹, Malick Fall², Lansana Toure¹, Mouctar Ndiaye^1
1Laboratoire de Mathématiques Appliquées (LABOMA), Faculté des Sciences de la Nature, Université Julius Nyerere de Kankan, Kankan, Guinée.
²Departement de Mathematiques, Université Gamal Abdel Nasser de Conackry, Conackry, Guinea.
DOI: 10.4236/jamp.2025.137137   PDF    HTML   XML   98 Downloads   462 Views   Citations

Abstract

The aim of this research article is to apply topological optimization gradient algorithm applied to 1D and 2D photonic charged slabs. We compute the topological gradient using min max method. We use an iterative algorithm for descending the topological gradient and with steps as the radius of circular holes, implemented under FreeFEM with the finite elements’ method.

Keywords

Photonic Crystals, Band Gap, Topological Optimization, Min-Max Principle, Electromagnetic, Elastic and Acoustic Waves, Topological Derivatives, Finite Elements, FreeFEM

Share and Cite:

1. Introduction

The year 1987 marked the discovery of periodic structures capable of trapping, confining, and guiding electromagnetic waves. Professors E. Yablonovitch [1] and S. John [2], through their respective works, independently concluded that in periodically structured media, there exist forbidden and allowed bands for electromagnetic waves. This breakthrough revolutionized technology has led to major innovations in optical fibers, compact discs, and telecommunications. This has led to a growing interest in photonic, which represents a major breakthrough in emerging technologies and information science [3]-[6].

Recent years have seen significant progress in the field of photonic crystal optimization, with a wide range of methodologies developed to tackle both structural and functional objectives. Among the most prominent approaches are level-set methods [7] [8], phase-field formulations [9] [10], and density-based techniques such as SIMP [11]. These methods have proven effective in generating complex dielectric patterns that enhance photonic performance, including band gap formation and wave guiding.

Despite these advances, most existing studies primarily focus on passive structures without accounting for the influence of external charges or electric field-induced phenomena. Furthermore, the role of boundary conditions Dirichlet vs. Neumann has rarely been explored in a comparative optimization framework.

In this work, we propose a novel perspective by investigating the topological optimization of charged photonic crystals using topological derivatives within a Mumford-Shah-type variational framework. Our approach addresses both Dirichlet and Neumann boundary conditions systematically, offering new insight into their respective impacts on algorithmic convergence and structural formation. To the best of our knowledge, such a dual-boundary analysis in the context of charged photonic systems remains largely unexplored, and we believe it provides valuable guidance for the design of efficient photonic structures.

Topological optimization provides an opportunity to obtain important information on the topology of the considered domain to optimize at least one criterion. Domain optimization is used today in many industrial environments, such as Airbus for the reduction of structures, the improvement of resistance to vibrations and many other areas of physics [12]-[14]. This gives us the idea of looking at the topological derivative, but this time using the recent work of [15]-[17]. Therefore, we study the topological derivative using the min-max method. for more information on this method the reader can consult the work of [18]. For more practical cases the reader can also consult the paper by [19], where the topological derivative of a functional linked to a linear thermoplastic problem was calculated. on the other hand in the paper by [17] the author established a practical case of the topological derivative linked to Helmholtz problems. The main objective in this article is to determine the topological derivative of the functional $J\left({\Omega }_r\right)=J\left({\Omega }_r,u_r\right)$ , where the perturbed domain ${\Omega }_r$ of $\Omega$ is defined by ${\Omega }_r=T_r\left(\Omega \right)$ or ${\Omega }_r=\Omega \backslash E_r$ depending on the derivative to be calculated but also to perform numerical simulations for both edge conditions.

The work program was as follows. Section 2 describes the modelization of photonic slabs. The application of topological derivatives to photonic slabs is described in section 3. To do this, we first establish the topological derivative for a Dirichlet condition and then for a Neumann condition using the minmax method. Section 4 establishes the numerical simulations for both the above conditions and Section 5 provides a conclusion and some possible extensions.

2. Model of Photonic Slabs

The modeling of photonic crystals aims to describe the propagation of electromagnetic waves in a medium with a periodically varying dielectric constant. Maxell’s equations govern the behavior of these waves in such structured materials.

2.1. Maxwell’s Equations for an Anisotropic Medium

Maxwell’s equations were derived for anisotropic media. The equation governing the electric field in domain $\Omega$ is given by:

$\nabla \cdot \left(ϵE\right)=\rho ,$ (1)

where $E$ is the electric field, $\rho$ represents the charge density, and $ϵ$ denotes the permittivity of the medium.

From Gauss’s theorem, we have:

${\displaystyle {\oint }_{\Sigma }E\text{d}S}=\frac{1}{{ϵ}_0}{\displaystyle {\iint }_V\rho \text{d}\tau }.$ (2)

Maxwell-Thompson’s law states that the flux of the magnetic field across a closed surface ${\Sigma }_f$ is always zero.

${\displaystyle {\iint }_{{\Sigma }_f}B\text{d}S}=0.$ (3)

Using the divergence theorem, this implies:

$\nabla \cdot B=0\text{in}\Omega .$ (4)

Because the magnetic field $H$ is related to the magnetic induction field $B$ by $B=\mu H$ , Maxwell-Faraday’s equation takes the following form.

$\nabla \times E=-\frac{\partial B}{\partial t}.$ (5)

Maxwell-Ampère’s equation is given by:

$\nabla \times B=\mu \left(j+J_D\right),$ (6)

where $\mu$ is the magnetic permeability, $j$ represents the current density ($j={\sigma }_{c}E$ , where ${\sigma }_c$ is the electrical conductivity), and $J_D=ϵ\frac{\partial E}{\partial t}$ corresponds to the displacement current.

Thus, we obtain the following system of Maxwell’s equations:

$\{\begin{array}{l}\nabla \cdot \left(ϵE\right)=\rho ,\hfill \\ \nabla \cdot B=0,\hfill \\ \nabla \times E=-\mu \frac{\partial H}{\partial t},\hfill \\ \nabla \times H=\mu \left(j+ϵ\frac{\partial E}{\partial t}\right).\hfill \end{array}$ (7)

2.2. Propagation of Electromagnetic Waves in a Good Conductor

In periodic structures such as photonic crystals, the medium is typically electrically neutral ($\rho =0$ ) and has a periodic dielectric constant $ϵ={ϵ}_0{ϵ}_r$ , where ${ϵ}_0$ and ${ϵ}_r$ are the permittivity of free space and the relative permittivity, respectively. The magnetic permeability is given by $\mu ={\mu }_0{\mu }_r$ , where ${\mu }_0$ is the permeability of free space and ${\mu }_r$ is the relative permeability.

For a good conductor, Maxwell’s equations simplify to:

$\{\begin{array}{l}\nabla \cdot \left({ϵ}_0{ϵ}_{r}E\right)=0,\hfill \\ \nabla \cdot B=0,\hfill \\ \nabla \times E=-\mu \frac{\partial H}{\partial t},\hfill \\ \nabla \times H=\mu \left(j+ϵ\frac{\partial E}{\partial t}\right).\hfill \end{array}$ (8)

Using the identity:

$\nabla \times \left(\nabla \times E\right)=\nabla \left(\nabla \cdot E\right)-{\nabla }^{2}E,$ (9)

and substituting into the above equations, we derive the wave equation:

$\text{Δ}E-ϵ\mu \frac{{\partial }^{2}E}{\partial t^2}-\mu {\sigma }_c\frac{\partial E}{\partial t}=0.$ (10)

Assuming a plane wave solution of the form $E=E_0{\text{e}}^{i\left(k\cdot x-\omega t\right)}$ , where $k$ is the wave vector, $\omega$ is the angular frequency, and $E_0$ is the wave amplitude, we rewrite the equation as:

$\text{Δ}E+\frac{{\omega }^2}{c^2}{\mu }_r{ϵ}^{\prime }E=0,$ (11)

where ${ϵ}^{\prime }={ϵ}_r+i\frac{{\sigma }_c}{\omega {ϵ}_0}$ and $c=\frac{1}{\sqrt{{\mu }_0{ϵ}_0}}$ represents the speed of light.

This equation describes the propagation of electromagnetic waves in a conducting photonic crystal, accounting for the dielectric and conductive properties of the material.

2.3. Propagation of Electromagnetic Waves in a Periodic Medium with Charge and Source

In an inhomogeneous medium with a source, we have $\rho \ne 0$ , $ϵ\ne {ϵ}_0$ , and the magnetic permeability $\mu$ satisfies $\mu ={\mu }_0{\mu }_r$ .

In this case, Maxwell’s equations become:

$\{\begin{array}{l}\text{div}\left({ϵ}_{r}E\right)=\frac{\rho }{{ϵ}_0},\\ \text{div}\left(H\right)=0,\\ \nabla \times E=-\frac{\partial H}{\partial t},\\ \nabla \times H=\mu \left(j+ϵ\frac{\partial E}{\partial t}\right).\end{array}$ (12)

Using the vector calculus identity

$\nabla \times \left(\nabla \times E\right)=\nabla \left(\nabla \cdot E\right)-\text{Δ}E,$ (13)

and combining it with (12), we obtain:

$\nabla \times \left(-\frac{\partial H}{\partial t}\right)=-\frac{\partial \nabla \times H}{\partial t}=-\frac{\partial \mu \left(j+ϵ\frac{\partial E}{\partial t}\right)}{\partial t}=-\frac{\partial \left(\mu j\right)}{\partial t}-\frac{\partial \left(\mu ϵ\frac{\partial E}{\partial t}\right)}{\partial t}.$

Simplifying, this leads to:

$\text{Δ}E-\frac{1}{ϵ}{\nabla }_x\rho -\frac{\partial \left(\mu j\right)}{\partial t}-\frac{\partial \left(\mu ϵ\frac{\partial E}{\partial t}\right)}{\partial t}=0.$

Finally, we obtain:

$\text{Δ}E-\frac{1}{ϵ}{\nabla }_x\rho -\epsilon \mu \frac{{\partial }^{2}E}{\partial t^2}-\mu {\sigma }_c\frac{\partial E}{\partial t}=0.$ (14)

By considering a plane harmonic wave solution of (14) in the form $E=u\left(x\right) {\text{e}}^{i\left(kx-wt\right)}$ , we obtain:

$\text{Δ}u\left(x\right)+\frac{w^2}{c^2}{\mu }_r{ϵ}^{\prime }u\left(x\right)=\frac{1}{ϵ}{\nabla }_x\rho {\text{e}}^{-i\left(kx-wt\right)},$ (15)

where ${ϵ}^{\prime }={ϵ}_r+i\frac{{\sigma }_c}{w{ϵ}_0}$ .

Remark 1 For the case of a monochromatic wave, we have:

$\text{Δ}u\left(x\right)+\frac{w^2}{c^2}{\mu }_r{ϵ}_{r}u\left(x\right)=\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right),$ (16)

where $k^2=\frac{w^2}{c^2}{\mu }_r{ϵ}_r$ represents the dispersion relation of the crystal.

3. Application of the Topological Derivative to Charged Photonic Slabs

In this section, we focus on the study of photonic crystals using topological optimization methods. Topological optimization is a branch of shape optimization that seeks an optimal shape by changing the topology of the initial domain. This is a rapidly developing subject in several fields.

Our research follows the studies of crystals using topological optimization methods by Ngom et al. ([20]), where they minimized a least-squares type functional and compliance for photonic and phononic crystals, respectively. Thus, in our study, we considered a shape functional that combines the two functionals used in ([20]).

Let $\Omega \subset {ℝ}^2$ (or ${ℝ}^3$ ) be an open bounded domain and $u_{\Omega }$ , a function of class ${\mathcal{C}}^2$ . We define the boundary operator of domain $\Omega$ as:

$\begin{array}{ll}{B}_{\Omega }:\hfill & \Omega \to \partial \Omega ,\hfill \\ \hfill & u_{\Omega }\to B_{\Omega }{u}_{\Omega }=\{\begin{array}{l}{u}_{\Omega },\\ \text{or},\\ \frac{\partial u_{\Omega }}{\partial n},\end{array}\hfill \end{array}$ (17)

and the partial differential equation:

$\{\begin{array}{l}\text{Δ}u\left(x\right)+\frac{w^2}{c^2}{\mu }_r{ϵ}_{r}u\left(x\right)=\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\text{in}\Omega ,\\ B_{\Omega }{u}_{\Omega }=0\text{on}\partial \Omega .\end{array}$ (18)

We then consider the functional $J$ defined by:

$J\left(\Omega \right)=j\left(u_{\Omega }\right)=\alpha {\displaystyle {\int }_{\Omega }{\left|u_{\Omega }-u_0\right|}^2\text{d}x}+\beta {\displaystyle {\int }_{\Omega }{\left|\nabla u_{\Omega }\right|}^2\text{d}x},$ (19)

where $u_{\Omega }$ is the solution of (18), and $\alpha ,\beta$ are real numbers, and $u_0$ is a function in $L^2\left(\Omega \right)$ .

Consider $x_0\in \Omega$ and let $r$ be a positive real number. We define ${\Omega }_r=\Omega \backslash \overline{E_r}$ , where $E_r=\left\{x_0+rE,E\subset \Omega \right\}$ .

Let $u_r$ , be the solution to the perturbed problem.

$\{\begin{array}{l}\text{Δ}u\left(x\right)+\frac{w^2}{c^2}{\mu }_r{ϵ}_{r}u\left(x\right)=\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\text{in}{\Omega }_r,\\ B_{\Omega }{u}_{{\Omega }_r}=0\text{on}\partial \Omega ,\\ B_{E_r}{u}_{{\Omega }_{\epsilon }}=0\text{on}\partial E_r,\end{array}$ (20)

with $k^2={ϵ}_r\frac{w^2}{c^2}=n^2\frac{w^2}{c^2}$ , where $n$ is the refractive index of the periodic medium.

We define the spaces $H_r$ and ${\tilde{H}}_r$ as follows:

$H_r=\left\{u\in H^1\left({\Omega }_{\epsilon }\right),B_{\Omega }u=0\text{on}\partial \Omega \right\},$ (21)

${\tilde{H}}_r=\left\{u\in H_r,B_{E}u=0\text{on}\partial E_r\right\}.$(22)

Consequently, we express the functional defined in (19) in the perturbed domain as:

$J\left({\Omega }_r\right)=j\left(u_{{\Omega }_r}\right)=\alpha {\displaystyle {\int }_{{\Omega }_r}{\left|u_{{\Omega }_{\epsilon }}-u_0\right|}^2\text{d}x}+\beta {\displaystyle {\int }_{{\Omega }_r}{\left|\nabla u_{{\Omega }_r}\right|}^2\text{d}x},$ (23)

where $u_{{\Omega }_r}$ is the solution of (20) and $\alpha ,\beta$ are real constants.

Thus, we determine the topological derivative $g\left(x_0\right)$ by seeking the asymptotic expansion of functional (19) in the case of Dirichlet or Neumann boundary conditions on the hole border.

In each case, we determine $g\left(x_0\right)$ and $f\left(r\right)$ such that:

$j\left(u_{{\Omega }_r}\right)-j\left(u_{\Omega }\right)=f\left(r\right)g\left(x_0\right)+o\left(f\left(r\right)\right).$ (24)

Remark 2 The expression of the topological gradient does not depend on the condition imposed on the boundary of the domain $\Omega$ , but rather on the condition considered at the boundary of the hole.

In the following theorem $N$ denotes the dimension of the workspace, $d$ the dimension of a subset $E$ of ${ℝ}^N$ , $r$ is the radius of the ball, and ${\alpha }_{N-d}$ is the volume of the unit ball in ${ℝ}^{N-d}$ .

3.1. Dirichlet Condition around the Hole

Let us associate with $r$ , $0<r\le R$ , the perturbed domain ${\Omega }_r=\Omega \backslash E_r$ , where by assumption, $\partial {\Omega }_r=\partial \Omega \cup \partial E_r$ , and $\partial \Omega \cap \partial E_r=\varnothing$ and $\partial E_r\in {\mathcal{C}}^{1,1}$ .

Let $u_{{\Omega }_r}$ be the solution for the following system:

$\{\begin{array}{l}\text{Δ}{u}_{{\Omega }_r}\left(x\right)+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_r}\left(x\right)=\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\text{in}{\Omega }_r,\\ u_{{\Omega }_r}=0\text{on}\partial \Omega ,\\ u_{{\Omega }_r}=0\text{on}\partial E_r,\end{array}$ (25)

which can be extented to $\Omega$ by intoducing the solution $E_r^0\to ℝ$ of the problem

$\begin{array}{l}\text{Δ}{u}_{{\Omega }_r}\left(x\right)+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_r}\left(x\right)\\ =\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\text{in}{E}_r^0\text{and}{u}_{{\Omega }_r}^0=u_{{\Omega }_r}\text{on}\partial E_r.\end{array}$ (26)

We suppose that ${\Omega }_r$ has two components ${\Omega }_r^m$ and ${\Omega }_r^0$ . ${\Omega }_r^m$ is the component of ${\Omega }_r$ for which $\partial \Omega$ is part of its boundary. ${\Omega }_r^0$ is the blind component of ${\Omega }_r$ whose boundary has an empty intersection with $\partial \Omega$ . The function $u_{{\Omega }_r}$ is distributed between the two components ${\Omega }_r^0$ and ${\Omega }_r^m$ as

$\begin{array}{l}{u}_{{\Omega }_r}=u_{{\Omega }_0}\text{in}{\Omega }_r^0\text{and}\\ \text{Δ}{u}_{{\Omega }_r}\left(x\right)+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_r}\left(x\right)=\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\text{in}{E}_r^0\text{in}{\Omega }_r^m,\end{array}$ (27)

$\{\begin{array}{l}{u}_{{\Omega }_r}=0\text{on}\partial \Omega ,\\ u_{{\Omega }_r}=0\text{on}\partial {\Omega }_r^m\cap \partial E_r,\end{array}$

Since $\partial E_r$ is made up of two disjoint boundary $\partial {\Omega }_r^0$ and $\partial {\Omega }_r^m\cap \partial E_r$ , we can construct an extension to $\Omega$ by defining the solution $E_r^0\leftarrow ℝ$

$\begin{array}{l}\text{Δ}{u}_{{\Omega }_r}\left(x\right)+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_r}\left(x\right)=\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\text{in}{E}_r^0\text{and}\\ u_{{\Omega }_r}^0=u_{{\Omega }_r}\text{on}\partial {\Omega }_r^m\cap \partial E_r\\ u_{{\Omega }_r}^0=u_{{\Omega }_r}\text{on}\partial {\Omega }_r^0.\end{array}$ (28)

For simplicity, we assumed that ${\Omega }_r^0$ is empty.

In the following, we also consider the functional defined in ${\Omega }_r$ , by

$J\left({\Omega }_r\right)=\beta {\displaystyle {\int }_{{\Omega }_r}{\left|\nabla u_{{\Omega }_r}\right|}^2\text{d}x}+\alpha {\displaystyle {\int }_{{\Omega }_r}{\left|u_{{\Omega }_r}-u_0\right|}^2\text{d}x}$ (29)

where $u_{{\Omega }_r}$ be the solution to the following problem

$\{\begin{array}{l}\text{Δ}{u}_{{\Omega }_r}\left(x\right)+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_r}\left(x\right)=\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\text{in}{\Omega }_r,\\ u_{{\Omega }_r}=0\text{on}\partial \Omega ,\\ u_{{\Omega }_r}=0\text{on}\partial E_r,\end{array}$ (30)

Considering a shape function $J$ defined by

$J\left(\Omega \right)=\alpha {\displaystyle {\int }_{\Omega }{\left|u_{\Omega }-u_0\right|}^2\text{d}x}+\beta {\displaystyle {\int }_{\Omega }{\left|\nabla u_{\Omega }\right|}^2\text{d}x}$ (31)

where $u_{\Omega }\in$ is solution to the variational problem

$\begin{array}{l}-{\displaystyle {\int }_{\Omega }\nabla u_{\Omega }\cdot \nabla v\text{d}x}+{\displaystyle {\int }_{\Omega }\frac{\partial u_{\Omega }}{\partial n}v\text{d}x}+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{\Omega }{u}_{\Omega }v\text{d}x}\\ ={\displaystyle {\int }_{\Omega }\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)v\text{d}x}\end{array}$ (32)

The shape functional associated the perforated domain is given by

$j\left({\chi }_r\left(x_0\right)\right)=J\left({\Omega }_r\right)=\alpha {\displaystyle {\int }_{{\Omega }_r}{\left|u_{{\Omega }_r}-u_0\right|}^2\text{d}x}+\beta {\displaystyle {\int }_{{\Omega }_r}{\left|\nabla u_{{\Omega }_r}\right|}^2\text{d}x}$ (33)

where $u_{{\text{Ω}}_r}$ is solution the variational problem

$\begin{array}{l}-{\displaystyle {\int }_{{\Omega }_r}\nabla u_{{\Omega }_r}\cdot \nabla v\text{d}x}+{\displaystyle {\int }_{{\Omega }_r}\frac{\partial u_{{\Omega }_r}}{\partial n}v\text{d}x}+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{{\Omega }_r}{u}_{{\Omega }_r}v\text{d}x}\\ ={\displaystyle {\int }_{{\Omega }_r}\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)v\text{d}x}\end{array}$ (34)

We aim to compute the topological derivative of the functional $J\left({\Omega }_r\right)$

$dJ=\underset{r\to 0}{\lim }\frac{J\left({\Omega }_r\right)-J\left(\Omega \right)}{{\alpha }_{N-d}{r}^{N-d}}.$

And for this purpose, we define the following set:

$H_r=\left\{u_{\Omega }\in H^1\left({\Omega }_r\right):u_{\Omega }=0\text{on}\partial \Omega ,u_{\Omega }=0\text{on}\partial E_r\right\}$ (35)

Our variational formulation (25) consists of finding $u_{{\text{Ω}}_r}\in H_r$ such that

$\begin{array}{l}-{\displaystyle {\int }_{{\Omega }_r}\nabla u_{{\Omega }_r}\cdot \nabla v\text{d}x}+{\displaystyle {\int }_{{\Omega }_r}\frac{\partial u_{{\Omega }_r}}{\partial n}v\text{d}x}+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{{\Omega }_r}{u}_{{\Omega }_r}v\text{d}x}\\ ={\displaystyle {\int }_{{\Omega }_r}\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)v\text{d}x}\forall v\in H\end{array}$ (36)

By taking $H_r=H$ , we have

$\begin{array}{l}-{\displaystyle {\int }_{{\Omega }_r}\nabla u_{{\Omega }_r}\cdot \nabla v\text{d}x}+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{{\Omega }_r}{u}_{{\Omega }_r}v\text{d}x}\\ ={\displaystyle {\int }_{{\Omega }_r}\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)v\text{d}x}\forall v\in H_r.\end{array}$ (37)

Thus, the Lagrangian dependent on $r$ will be written in the form :

$\begin{array}{c}L\left(r,\varphi ,\Phi \right)=\beta {\displaystyle {\int }_{{\Omega }_r}{\left|\nabla \varphi \right|}^2\text{d}x}+\alpha {\displaystyle {\int }_{{\Omega }_r}{\left|\varphi -u_0\right|}^2\text{d}x}-{\displaystyle {\int }_{{\Omega }_r}\nabla \varphi \cdot \nabla v\text{d}x}+{\displaystyle {\int }_{{\Omega }_r}\frac{\partial \varphi }{\partial n}v\text{d}x}\\ +\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{{\Omega }_r}\varphi v\text{d}x}-{\displaystyle {\int }_{{\Omega }_r}\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)v\text{d}x}\end{array}$

$J\left({\Omega }_r\right)=\underset{\varphi \in H_r}{\inf }\underset{\Phi \in H_r}{\sup }L\left(r,\varphi ,\Phi \right).$

From this, we can evaluate the derivative of the Lagrangian, dependent on $r$ , with respect to $\varphi$ .

$\begin{array}{c}{d}_{\varphi }L\left(r,\varphi ,\Phi ,{\varphi }^{\prime }\right)={\displaystyle {\int }_{{\Omega }_r}2\beta \nabla \varphi \cdot \nabla {\varphi }^{\prime }\text{d}x}+{\displaystyle {\int }_{{\Omega }_r}2\alpha \left(\varphi -u_0\right){\varphi }^{\prime }\text{d}x}\\ -{\displaystyle {\int }_{{\Omega }_r}\nabla {\varphi }^{\prime }\cdot \nabla \Phi \text{d}x}+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{{\Omega }_r}{\varphi }^{\prime }\cdot \Phi \text{d}x}.\end{array}$

Subsequently, we obtained the variational formulation of the adjoint state equation given by $d_{\varphi }L\left(0,u_{{\Omega }_0},p_0,{\varphi }^{\prime }\right)=0$ , where $u_{{\Omega }_0}=u_{{\Omega }_r}$ for $r=0$ . Find $p_0\in H_0^1\left(\Omega \right)$ such that

$\begin{array}{l}{\displaystyle {\int }_{\Omega }2\beta \nabla u_{{\Omega }_0}\cdot \nabla {\varphi }^{\prime }\text{d}x}+{\displaystyle {\int }_{\Omega }2\alpha \left(u_{{\Omega }_0}-u_0\right){\varphi }^{\prime }\text{d}x}-{\displaystyle {\int }_{\Omega }\nabla {\varphi }^{\prime }\cdot \nabla p_0\text{d}x}\\ +\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{{\Omega }_{ϵ}}{\varphi }^{\prime }\cdot p_0\text{d}x}=0.\end{array}$ (38)

And we have

${\displaystyle {\int }_{\Omega }\left[2\beta \nabla u_{{\Omega }_0}\cdot \nabla {\varphi }^{\prime }+2\alpha \left(u_{{\Omega }_0}-u_0\right){\varphi }^{\prime }-\nabla {\varphi }^{\prime }\cdot \nabla p_0+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\varphi }^{\prime }\cdot p_0\right]\text{d}x}=0.$ (39)

Next, we derive the Lagrangian with respect to $\Phi$ .

$\begin{array}{c}{d}_{\Phi }L\left(r,\varphi ,\Phi ,{\Phi }^{\prime }\right)=-{\displaystyle {\int }_{{\Omega }_r}\nabla \varphi \cdot \nabla {\Phi }^{\prime }\text{d}x}-{\displaystyle {\int }_{{\Omega }_r}\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right){\Phi }^{\prime }\text{d}x}\\ +{\displaystyle {\int }_{{\Omega }_r}\frac{w^2}{c^2}{\mu }_r{ϵ}_r\varphi \cdot {\Phi }^{\prime }\text{d}x}.\end{array}$

The initial state $u_{{\Omega }_0}=u_{\Omega }$ is a solution of $d_{\Phi }L\left(0,u_{{\Omega }_0},0,{\Phi }^{\prime }\right)=0$ $\forall$ ${\Phi }^{\prime }\in H_0^1$ and in this case, we have:

$-{\displaystyle {\int }_{\Omega }\nabla u_{{\Omega }_0}\cdot \nabla {\Phi }^{\prime }\text{d}x}-{\displaystyle {\int }_{\Omega }\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right){\Phi }^{\prime }\text{d}x}+{\displaystyle {\int }_{\Omega }\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_0}\cdot {\Phi }^{\prime }\text{d}x}=0.$

Then, we have:

${\displaystyle {\int }_{\Omega }\left[-\nabla u_{{\Omega }_0}\cdot \nabla {\Phi }^{\prime }-\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right){\Phi }^{\prime }+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_0}\cdot {\Phi }^{\prime }\right]\text{d}x}=0.$

The state $u_{{\Omega }_r}$ for all $r\ge 0$ satisfies

${\displaystyle {\int }_{{\Omega }_r}\left[-\nabla u_{{\Omega }_r}\cdot \nabla {\Phi }^{\prime }-\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right){\Phi }^{\prime }+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_r}\cdot {\Phi }^{\prime }\right]\text{d}x}=0,\forall {\Phi }^{\prime }\in H.$

In the following, we aim to determine the derivative of the Lagrangian, with respect to $r$ . To achieve this, let us first compute the quotient

$\frac{L\left(r,\varphi ,\Phi \right)-L\left(0,\varphi ,\Phi \right)}{s}.$

$\begin{array}{l}L\left(r,\varphi ,\Phi \right)-L\left(0,\varphi ,\Phi \right)\\ =\beta {\displaystyle {\int }_{{\Omega }_r}{\left|\nabla \varphi \right|}^2\text{d}x}+\alpha {\displaystyle {\int }_{{\Omega }_r}{\left|\varphi -u_0\right|}^2\text{d}x}-{\displaystyle {\int }_{{\Omega }_r}\nabla \varphi \cdot \nabla \Phi \text{d}x}+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{{\Omega }_r}\varphi \Phi \text{d}x}\\ -{\displaystyle {\int }_{{\Omega }_r}\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\Phi \text{d}x}-\left[\beta {\displaystyle {\int }_{\Omega }{\left|\nabla \varphi \right|}^2\text{d}x}+\alpha {\displaystyle {\int }_{\Omega }{\left|\varphi -u_0\right|}^2\text{d}x}\right]\\ -\left[-{\displaystyle {\int }_{\Omega }\nabla \varphi \cdot \nabla \Phi \text{d}x}+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{\Omega }\varphi \Phi \text{d}x}-{\displaystyle {\int }_{\Omega }\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\Phi \text{d}x}\right]\end{array}$

$\begin{array}{l}L\left(ϵ,\varphi ,\Phi \right)-L\left(0,\varphi ,\Phi \right)\\ =-\left[{\displaystyle {\int }_{E_r}\beta {\left|\nabla \varphi \right|}^2}+\alpha {\left|\varphi -u_0\right|}^2-\nabla \varphi \cdot \nabla \Phi -\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\Phi +\frac{w^2}{c^2}{\mu }_r{ϵ}_r\varphi \Phi \right]\text{d}x.\end{array}$

For $d=0$ , $\omega =\left\{x_0\right\}$ , $E_r=\left\{x\in {ℝ}^N:\left|x-x_0\right|\le ϵ\right\}=\overline{B}\left(x_0,r\right)$ .

$\begin{array}{c}{d}_{s}L\left(0,\varphi ,\Phi \right)=\underset{s\to 0}{\lim }-\frac{1}{\left|B\left(x_0,r\right)\right|}[{\displaystyle {\int }_{B\left(x_0,r\right)}\beta {\left|\nabla \varphi \right|}^2}+\alpha {\left|\varphi -u_0\right|}^2-\nabla \varphi \cdot \nabla \Phi \begin{array}{c}\\ \end{array}\\ -\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\Phi +M\varphi \Phi ]\text{d}x\\ =-\beta {\left|\nabla \varphi \left(x_0\right)\right|}^2-\alpha {\left|\varphi \left(x_0\right)-u_0\left(x_0\right)\right|}^2+\nabla \varphi \left(x_0\right)\cdot \nabla \Phi \left(x_0\right)\\ +\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\Phi \left(x_0\right)-\frac{w^2}{c^2}{\mu }_r{ϵ}_r\varphi \left(x_0\right)\Phi \left(x_0\right).\end{array}$

By evaluating the last equation at the point $u_{{\Omega }_0},p_0$ , we obtain:

$\begin{array}{c}{d}_{s}L\left(0,u_{{\Omega }_0},p_0\right)=-\beta {\left|\nabla u_{{\Omega }_0}\left(x_0\right)\right|}^2-\alpha {\left|u_{{\Omega }_0}\left(x_0\right)-u_0\left(x_0\right)\right|}^2+\nabla u_{{\Omega }_0}\left(x_0\right)\cdot \nabla p_0\left(x_0\right)\\ +\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)p_0\left(x_0\right)-\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_0}\left(x_0\right)p_0\left(x_0\right).\end{array}$

Hence, if $0<d\le N-1$ , we have:

$\begin{array}{l}\frac{L\left(r,\varphi ,\Phi \right)-L\left(0,\varphi ,\Phi \right)}{s}\\ =-\frac{1}{\left|E_r\right|}\left[{\displaystyle {\int }_{E_r}\beta {\left|\nabla \varphi \right|}^2}+\alpha {\left|\varphi -u_0\right|}^2-\nabla \varphi \cdot \nabla \Phi -\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\Phi +\frac{w^2}{c^2}{\mu }_r{ϵ}_r\varphi \Phi \right]\text{d}x\\ =-\frac{1}{{\alpha }_{N-d}{r}^{N-d}}\left[{\displaystyle {\int }_{E_r}\beta {\left|\nabla \varphi \right|}^2}+\alpha {\left|\varphi -u_0\right|}^2-\nabla \varphi \cdot \nabla \Phi -\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\Phi \right]\text{d}x\\ -\frac{1}{{\alpha }_{N-d}{r}^{N-d}}\left[{\displaystyle {\int }_{E_r}\frac{w^2}{c^2}{\mu }_r{ϵ}_r\varphi \Phi }\right]\text{d}x\\ \to -\left[{\displaystyle {\int }_E\beta {\left|\nabla \varphi \right|}^2}+\alpha {\left|\varphi -u_0\right|}^2-\nabla \varphi \cdot \nabla \Phi -\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\Phi +\frac{w^2}{c^2}{\mu }_r{ϵ}_r\varphi \Phi \right]\text{d}{H}^d.\end{array}$

Therefore, taking the ast result at the point $u_{{\Omega }_0},p_0$ becomes:

$\begin{array}{c}{d}_{s}L\left(0,u_{{\Omega }_0},p_0\right)=-[{\displaystyle {\int }_E\beta {\left|\nabla u_{{\Omega }_0}\right|}^2}+\alpha {\left|u_{{\Omega }_0}-u_0\right|}^2-\nabla u_{{\Omega }_0}\cdot \nabla p_0\begin{array}{c}\\ \end{array}\\ -\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)p_0+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_0}{p}_0]\text{d}{H}^d.\end{array}$

We now define $R\left(r\right)$ as

$R\left(r\right)={\displaystyle {\int }_0^1{d}_{x}L}\left(r,u_{{\Omega }_0}+\Psi \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right),p_0,\left(\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{s}\right)\right)\text{d}\Psi .$

By substituting ${\varphi }^{\prime }=\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{r}$ and $\Psi =\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{2}$ into the adjoint equation for $p_0$ , we obtain:

$\begin{array}{c}R\left(r\right)={\displaystyle {\int }_{{\Omega }_r}2\beta \nabla }\left(\frac{u_{{\Omega }_r}+u_{{\Omega }_0}}{2}\right)\cdot \nabla \left(\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{s}\right)\text{d}x\\ +{\displaystyle {\int }_{{\Omega }_r}2\alpha }\left[\left(\frac{u_{{\Omega }_r}+u_{{\Omega }_0}}{2}\right)-u_0\right]\left(\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{s}\right)\text{d}x\\ -{\displaystyle {\int }_{{\Omega }_r}\nabla \left(\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{s}\right)\cdot \nabla p_0\text{d}x}+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{{\Omega }_r}\left(\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{s}\right)\cdot p_0\text{d}x}\\ =\frac{1}{s}[{\displaystyle {\int }_{{\Omega }_r}2\beta \nabla }\left(\frac{u_{{\Omega }_r}+u_{{\Omega }_0}}{2}\right)\cdot \nabla \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\text{d}x\\ +2\alpha \left[\left(\frac{u_{{\Omega }_r}+u_{{\Omega }_0}}{2}\right)-u_0\right]\left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)]\text{d}x\\ -\frac{1}{s}\left[{\displaystyle {\int }_{{\Omega }_r}\nabla }\left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot \nabla p_0-\frac{w^2}{c^2}{\mu }_r{ϵ}_r\left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot p_0\right]\text{d}x\end{array}$

$\begin{array}{c}R\left(r\right)=\frac{1}{s}\left[{\displaystyle {\int }_{{\Omega }_r}2\beta }\left(\nabla \left(\frac{u_{{\Omega }_r}+u_{{\Omega }_0}}{2}\right)-\nabla u_{{\Omega }_0}+\nabla u_{{\Omega }_0}\right)\cdot \nabla \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\right]\text{d}x\\ +\frac{2\alpha }{s}{\displaystyle {\int }_{{\Omega }_r}\left[\left(\frac{u_{{\Omega }_r}+u_{{\Omega }_0}}{2}\right)-u_{{\Omega }_0}+u_{{\Omega }_0}-u_0\right]\left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\text{d}x}\\ -\frac{1}{s}\left[{\displaystyle {\int }_{{\Omega }_r}\nabla \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot \nabla p_0}-\frac{w^2}{c^2}{\mu }_r{ϵ}_r\left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot p_0\right]\text{d}x\\ =\frac{1}{s}\left[{\displaystyle {\int }_{{\Omega }_r}2\beta }\left(\nabla \left(\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{2}\right)+\nabla u_{{\Omega }_0}\right)\cdot \nabla \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\right]\text{d}x\\ +\frac{1}{s}\left[{\displaystyle {\int }_{{\Omega }_r}2\alpha }\left[\left(\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{2}\right)+u_{{\Omega }_0}-u_0\right]\left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\right]\text{d}x\\ -\frac{1}{s}\left[{\displaystyle {\int }_{{\Omega }_r}\nabla \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot \nabla p_0}-\frac{w^2}{c^2}{\mu }_r{ϵ}_r\left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot p_0\right]\text{d}x\end{array}$

$\begin{array}{c}R\left(r\right)={\displaystyle {\int }_{{\Omega }_r}\left(\beta {\left|\nabla \left(\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{\sqrt{s}}\right)\right|}^2+\alpha {\left|\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{\sqrt{s}}\right|}^2\right)\text{d}x}\\ +\frac{1}{s}\left[{\displaystyle {\int }_{{\Omega }_r}2\beta \nabla u_{{\Omega }_0}\nabla \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)}+2\alpha \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\left(u_{{\Omega }_0}-u_0\right)\right]\text{d}x\\ -\frac{1}{s}\left[{\displaystyle {\int }_{{\Omega }_r}\nabla \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot \nabla p_0}-\frac{w^2}{c^2}{\mu }_r{ϵ}_r\left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot p_0\right]\text{d}x\end{array}$

Thus, for all $u_{{\Omega }_0}\in H_0^1\left(\Omega \right)$ , equation (38) becomes:

$\begin{array}{l}{\displaystyle {\int }_{{\Omega }_r}\left(-\nabla u_{{\Omega }_r}\cdot \nabla v-\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)v+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_r}v\right)\text{d}x}\\ =-\left[{\displaystyle {\int }_{E_r}-\nabla u_{{\Omega }_r}\cdot \nabla v}-\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)v\right]-\left[{\displaystyle {\int }_{E_r}\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_r}v}\right]\text{d}x\\ ={\displaystyle {\int }_{\partial E_r}\frac{\partial u_{{\Omega }_0}}{\partial n}v\text{d}{H}^{N-1}},\forall v\in H_0^1\left(\Omega \right).\end{array}$

Now, considering the assumption about $E$ , $\partial E_r\in {\mathcal{C}}^{1,1}$ and

$\frac{\partial u_{{\Omega }_0}}{\partial n}=\nabla u_{{\Omega }_0}\cdot n{\Omega }_r=-\nabla u_{{\Omega }_0}\cdot \nabla dE\text{on}{E}_r.$

And in this case, we have:

$\begin{array}{l}{\displaystyle {\int }_{{\Omega }_r}\left(-\nabla u_{{\Omega }_r}\cdot \nabla v-\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)v+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_r}v\right)\text{d}x}\\ =-{\displaystyle {\int }_{\partial E_r}\nabla u_{{\Omega }_0}\cdot \nabla dEv\text{d}{H}^{N-1}}.\end{array}$ (40)

Taking the difference between equation (38) and equation (40), we obtain

${\displaystyle {\int }_{{\Omega }_r}-\nabla \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\nabla v}={\displaystyle {\int }_{\partial E_r}\nabla u_{{\Omega }_0}\cdot \nabla dEv\text{d}{H}^{N-1}}.$

The adjoint equation for $r\ge 0$ yields:

$\begin{array}{l}{\displaystyle {\int }_{{\Omega }_r}2\beta \nabla u_{{\Omega }_0}\cdot \nabla {\varphi }^{\prime }\text{d}x}+{\displaystyle {\int }_{{\Omega }_r}2\alpha \left(u_{{\Omega }_0}-u_0\right){\varphi }^{\prime }\text{d}x}-{\displaystyle {\int }_{{\Omega }_r}\nabla {\varphi }^{\prime }\cdot \nabla p_r\text{d}x}\\ +\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{{\Omega }_r}{\varphi }^{\prime }\cdot p_r\text{d}x}=0.\end{array}$ (41)

By taking ${\varphi }^{\prime }=u_{{\Omega }_r}-u_{{\Omega }_0}$ in equation (41), we have:

$\begin{array}{l}{\displaystyle {\int }_{{\Omega }_r}2\beta \nabla u_{{\Omega }_0}\cdot \nabla \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\text{d}x}+{\displaystyle {\int }_{{\Omega }_r}2\alpha \left(u_{{\Omega }_0}-u_0\right)\left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\text{d}x}\\ -{\displaystyle {\int }_{{\Omega }_r}\nabla \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot \nabla p_r}+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{{\Omega }_r}\left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot p_r\text{d}x}=0.\end{array}$

$\begin{array}{l}{\displaystyle {\int }_{{\Omega }_r}2\beta \nabla u_{{\Omega }_0}\cdot \nabla \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)+2\alpha \left(u_{{\Omega }_0}-u_0\right)\left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\text{d}x}\\ ={\displaystyle {\int }_{{\Omega }_r}\nabla \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot \nabla p_r}-\frac{w^2}{c^2}{\mu }_r{ϵ}_r\left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot p_r\text{d}x\end{array}$

The final equation for $R\left(r\right)$ becomes:

$\begin{array}{c}R\left(r\right)={\displaystyle {\int }_{{\Omega }_r}\left[\beta {\left|\nabla \left(\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{\sqrt{s}}\right)\right|}^2+\alpha {\left|\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{\sqrt{s}}\right|}^2\right]\text{d}x}\\ +\frac{1}{s}\left[{\displaystyle {\int }_{\partial {\Omega }_r^m\cap \partial E_r}\nabla u_{{\Omega }_0}\cdot \nabla d\omega p_0\text{d}{H}^{N-1}}\right]\text{d}x\\ -\frac{1}{s}{\displaystyle {\int }_{{\Omega }_r}\nabla \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot \nabla p_r}-\frac{w^2}{c^2}{\mu }_r{ϵ}_r\left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot p_r\text{d}x.\end{array}$

THEOREME Let $0\le d<N$ , and $s={\alpha }_{N-d}{r}^{N-d}$ . The topological derivative exists if and only if the following limit is satisfied:

$l=\underset{r\to 0}{\text{lim}}\left(l_0\left(r\right)+l_1\left(r\right)\right),$

exists with

$l_0\left(r\right)={\displaystyle {\int }_{{\Omega }_r}\beta {\left|\nabla \left(\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{\sqrt{s}}\right)\right|}^2+\alpha {\left|\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{\sqrt{s}}\right|}^2}$

and

$\begin{array}{c}{l}_1\left(r\right)=\frac{1}{s}\left[{\displaystyle {\int }_{\partial {\Omega }_r^m\cap \partial E_r}\nabla u_{{\Omega }_0}\cdot \nabla dE{p}_0\text{d}{H}^{N-1}}\right]\text{d}x-\frac{1}{s}{\displaystyle {\int }_{{\Omega }_r}\nabla \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot \nabla p_r}\\ -\frac{w^2}{c^2}{\mu }_r{ϵ}_r\left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot p_r\text{d}x.\end{array}$

Moreover, the topological derivative of the function is given by the expression:

$\begin{array}{c}dJ=\underset{r\to 0}{\lim }\frac{J\left({\Omega }_r\right)-J\left(\Omega \right)}{{\alpha }_{N-d}{r}^{N-d}}\\ =l-[{\displaystyle {\int }_E\beta {\left|\nabla u_{{\Omega }_0}\right|}^2}+\alpha {\left|u_{{\Omega }_0}-u_0\right|}^2-\nabla u_{{\Omega }_0}\cdot \nabla p_0\begin{array}{c}\\ \end{array}\\ -\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)p_0+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_0}{p}_0]\text{d}{H}^d.\end{array}$

where $p_0,u_{{\Omega }_0}$ are solutions of systems

${\displaystyle {\int }_{\Omega }\left[2\beta \nabla u_{{\Omega }_0}\cdot \nabla {\varphi }^{\prime }+2\alpha \left(u_{{\Omega }_0}-u_0\right){\varphi }^{\prime }-\nabla {\varphi }^{\prime }\cdot \nabla p_0+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\varphi }^{\prime }\cdot p_0\right]\text{d}x}=0.$

In particular for $d=0$ ,

$\begin{array}{c}dJ=l-\beta {\left|\nabla u_{{\Omega }_0}\left(x_0\right)\right|}^2-\alpha {\left|u_{{\Omega }_0}\left(x_0\right)-u_0\left(x_0\right)\right|}^2+\nabla u_{{\Omega }_0}\left(x_0\right)\cdot \nabla p_0\left(x_0\right)\\ +\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)p_0\left(x_0\right)-\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_0}\left(x_0\right)p_0\left(x_0\right).\end{array}$

3.2. Neumann Condition

In this section, we consider the case of a Neumann condition, and the perturbed problem thus becomes:

$\{\begin{array}{l}\Delta u_{{\Omega }_r}\left(x\right)+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_r}\left(x\right)=\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\text{in}{\Omega }_r,\\ u_{{\Omega }_r}=0\text{on}\partial \Omega ,\\ \frac{\partial u_{{\Omega }_r}}{\partial n}=0\text{on}\partial E_r.\end{array}$ (42)

First of all, the calculations were the same. The changes were simply in the spaces considered. The techniques for calculating the Lagrangian derivative with respect to the variables remain the same. For this, we do not need to go into all the details, as we did in the case of the Dirichlet condition. And we can be extended to $\Omega$ by introducing the solution $E_r^0\to ℝ$ of the problem

$\begin{array}{l}\Delta u_{{\Omega }_r}\left(x\right)+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_r}\left(x\right)=\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\text{in}{E}_r^0\text{and}\\ u_{{\Omega }_r}^0=u_{{\Omega }_r}\text{on}\partial E_r.\end{array}$ (43)

We suppose that ${\Omega }_r$ has two components ${\Omega }_r^m$ and ${\Omega }_r^0$ . ${\Omega }_r^m$ is the component of ${\Omega }_r$ for which $\partial \Omega$ is part of its boundary. ${\Omega }_r^0$ is the blind component of ${\Omega }_r$ whose boundary has an empty intersection with $\partial \Omega$ . The function $u_{{\Omega }_r}$ is distributed between the two components ${\Omega }_r^0$ and ${\Omega }_r^m$ as

$\begin{array}{l}{u}_{{\Omega }_r}=u_{{\Omega }_0}\text{in}{\Omega }_r^0\text{and}\\ \Delta u_{{\Omega }_r}\left(x\right)+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_r}\left(x\right)=\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\text{in}{E}_r^0\text{in}{\Omega }_r^m,\end{array}$ (44)

$\{\begin{array}{l}{u}_{{\Omega }_r}=0\text{on}\partial \Omega ,\\ u_{{\Omega }_r}=0\text{on}\partial {\Omega }_r^m\cap \partial E_r,\end{array}$

Since $\partial E_r$ is made up of two disjoint boundary $\partial {\Omega }_r^0$ and $\partial {\Omega }_r^m\cap \partial E_r$ , we can construct an extension to $\Omega$ by defining the solution $E_r^0\to ℝ$

$\begin{array}{l}\Delta u_{{\Omega }_r}\left(x\right)+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_r}\left(x\right)=\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\text{in}{E}_r^0\text{and}\\ u_{{\Omega }_r}^0=u_{{\Omega }_r}\text{on}\partial {\Omega }_r^m\cap \partial E_r\\ u_{{\Omega }_r}^0=u_{{\Omega }_r}\text{on}\partial {\Omega }_r^0.\end{array}$ (45)

For simplicity, we assumed that ${\Omega }_r^0$ is empty. For this purpose, we define the following set:

${\mathscr{H}}_r=\left\{u_{\Omega }\in H^1\left({\Omega }_r\right):u_{\Omega }=0\text{on}\partial \Omega \right\}$ (46)

Our variational formulation (42) consists of finding $u_{{\Omega }_r}\in {\mathscr{H}}_r$ such that

$\begin{array}{l}-{\displaystyle {\int }_{{\Omega }_r}\nabla u_{{\Omega }_r}\cdot \nabla v\text{d}x}+{\displaystyle {\int }_{\partial {\Omega }_r}\frac{\partial u_{{\Omega }_r}}{\partial n}v\text{d}\sigma }+{\displaystyle {\int }_{\partial E_r}\frac{\partial u_{{\Omega }_r}}{\partial n}v\text{d}\sigma }+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{{\Omega }_r}{u}_{{\Omega }_r}v\text{d}x}\\ ={\displaystyle {\int }_{{\Omega }_r}\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)v\text{d}x}\forall v\in {\tilde{\mathscr{H}}}_r\end{array}$(47)

where ${\tilde{\mathscr{H}}}_r$ is define by

${\tilde{\mathscr{H}}}_r=\left\{v\in H^1\left({\Omega }_r\right):v=0\text{on}\partial \Omega \right\}={\mathscr{H}}_r$

And we have

$\begin{array}{l}-{\displaystyle {\int }_{{\Omega }_r}\nabla u_{{\Omega }_r}\cdot \nabla v\text{d}x}+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{{\Omega }_r}{u}_{{\Omega }_r}v\text{d}x}\\ ={\displaystyle {\int }_{{\Omega }_r}\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)v\text{d}x}\forall v\in {\tilde{\mathscr{H}}}_r.\end{array}$(48)

Thus, the Lagrangian dependent on $r$ defined from $\left[0,R\right]\times {\mathscr{H}}_r$ to values in $ℝ$ will be written in the form:

$\begin{array}{c}L\left(r,\varphi ,\Phi \right)=\beta {\displaystyle {\int }_{{\Omega }_r}{\left|\nabla \varphi \right|}^2\text{d}x}+\alpha {\displaystyle {\int }_{{\Omega }_r}{\left|\varphi -u_0\right|}^2\text{d}x}-{\displaystyle {\int }_{{\Omega }_r}\nabla \varphi \cdot \nabla v\text{d}x}+{\displaystyle {\int }_{{\Omega }_r}\frac{\partial \varphi }{\partial n}v\text{d}x}\\ +\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{{\Omega }_r}\varphi v\text{d}x}-{\displaystyle {\int }_{{\Omega }_r}\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)v\text{d}x}\end{array}$

The derivative of the Lagrangian, dependent on $r$ , with respect to $\varphi$ is given by:

$\begin{array}{c}{d}_{\varphi }L\left(r,\varphi ,\Phi ,{\varphi }^{\prime }\right)={\displaystyle {\int }_{{\Omega }_r}2\beta \nabla \varphi \cdot \nabla {\varphi }^{\prime }\text{d}x}+{\displaystyle {\int }_{{\Omega }_r}2\alpha \left(\varphi -u_0\right){\varphi }^{\prime }\text{d}x}\\ -{\displaystyle {\int }_{{\Omega }_r}\nabla {\varphi }^{\prime }\cdot \nabla \Phi \text{d}x}+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{{\Omega }_r}{\varphi }^{\prime }\cdot \Phi \text{d}x}.\end{array}$

And the variational formulation of the adjoint state equation is given by: find $p_0\in H_0^1\left(\Omega \right)$ such that

$\begin{array}{l}{\displaystyle {\int }_{\Omega }2\beta \nabla u_{{\Omega }_0}\cdot \nabla {\varphi }^{\prime }\text{d}x}+{\displaystyle {\int }_{\Omega }2\alpha \left(u_{{\Omega }_0}-u_0\right){\varphi }^{\prime }\text{d}x}-{\displaystyle {\int }_{\Omega }\nabla {\varphi }^{\prime }\cdot \nabla p_0\text{d}x}\\ +\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\displaystyle {\int }_{{\Omega }_{ϵ}}{\varphi }^{\prime }\cdot p_0\text{d}x}=0.\end{array}$ (49)

And we have

${\displaystyle {\int }_{\Omega }\left[2\beta \nabla u_{{\Omega }_0}\cdot \nabla {\varphi }^{\prime }+2\alpha \left(u_{{\Omega }_0}-u_0\right){\varphi }^{\prime }-\nabla {\varphi }^{\prime }\cdot \nabla p_0+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\varphi }^{\prime }\cdot p_0\right]\text{d}x}=0.$ (50)

The derive of the Lagrangian with respect to $\text{Φ}$ verifies:

${\displaystyle {\int }_{\Omega }\left[-\nabla u_{{\Omega }_0}\cdot \nabla {\Phi }^{\prime }-\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right){\Phi }^{\prime }+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_0}\cdot {\Phi }^{\prime }\right]\text{d}x}=0.$

The state $u_{{\text{Ω}}_r}$ for all $r\ge 0$ satisfies

${\displaystyle {\int }_{{\Omega }_r}\left[-\nabla u_{{\Omega }_r}\cdot \nabla {\Phi }^{\prime }-\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right){\Phi }^{\prime }+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_r}\cdot {\Phi }^{\prime }\right]\text{d}x}=0,\forall {\Phi }^{\prime }\in H.$

For $d=0$ , $\omega =\left\{x_0\right\}$ , $E_r=\left\{x\in {ℝ}^N:\left|x-x_0\right|\le ϵ\right\}=\overline{B}\left(x_0,r\right)$ .

$\begin{array}{c}{d}_{s}L\left(0,\varphi ,\Phi \right)=\underset{s\to 0}{\lim }-\frac{1}{\left|B\left(x_0,r\right)\right|}[{\displaystyle {\int }_{B\left(x_0,r\right)}\beta {\left|\nabla \varphi \right|}^2}+\alpha {\left|\varphi -u_0\right|}^2-\nabla \varphi \cdot \nabla \Phi \begin{array}{c}\\ \end{array}\\ -\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\Phi +M\varphi \Phi ]\text{d}x\\ =-\beta {\left|\nabla \varphi \left(x_0\right)\right|}^2-\alpha {\left|\varphi \left(x_0\right)-u_0\left(x_0\right)\right|}^2+\nabla \varphi \left(x_0\right)\cdot \nabla \Phi \left(x_0\right)\\ +\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)\Phi \left(x_0\right)-\frac{w^2}{c^2}{\mu }_r{ϵ}_r\varphi \left(x_0\right)\Phi \left(x_0\right).\end{array}$

By evaluating the last equation at the point $u_{{\text{Ω}}_0},p_0$ , we obtain:

$\begin{array}{c}{d}_{s}L\left(0,u_{{\Omega }_0},p_0\right)=-\beta {\left|\nabla u_{{\Omega }_0}\left(x_0\right)\right|}^2-\alpha {\left|u_{{\Omega }_0}\left(x_0\right)-u_0\left(x_0\right)\right|}^2+\nabla u_{{\Omega }_0}\left(x_0\right)\cdot \nabla p_0\left(x_0\right)\\ +\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)p_0\left(x_0\right)-\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_0}\left(x_0\right)p_0\left(x_0\right).\end{array}$

Hence, if $0<d\le N-1$ , we have: Therefore, taking the ast result at the point $u_{{\Omega }_0},p_0$ becomes:

$\begin{array}{c}{d}_{s}L\left(0,u_{{\Omega }_0},p_0\right)=-[{\displaystyle {\int }_E\beta {\left|\nabla u_{{\Omega }_0}\right|}^2}+\alpha {\left|u_{{\Omega }_0}-u_0\right|}^2-\nabla u_{{\Omega }_0}\cdot \nabla p_0\begin{array}{c}\\ \end{array}\\ -\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)p_0+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_0}{p}_0]\text{d}{H}^d.\end{array}$

We now define $R\left(r\right)$ as

$R\left(r\right)={\displaystyle {\int }_0^1{d}_{x}L}\left(r,u_{{\Omega }_0}+\Psi \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right),p_0,\left(\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{s}\right)\right)\text{d}\Psi .$

And the formulas of $R\left(r\right)$ becomes:

$\begin{array}{c}R\left(r\right)={\displaystyle {\int }_{{\Omega }_r}\left[\beta {\left|\nabla \left(\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{\sqrt{s}}\right)\right|}^2+\alpha {\left|\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{\sqrt{s}}\right|}^2\right]\text{d}x}\\ +\frac{1}{s}\left[{\displaystyle {\int }_{\partial {\Omega }_r^m\cap \partial E_r}\nabla u_{{\Omega }_0}\cdot \nabla d\omega p_0\text{d}{H}^{N-1}}\right]\text{d}x\\ -\frac{1}{s}{\displaystyle {\int }_{{\Omega }_r}\nabla \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot \nabla p_r}-\frac{w^2}{c^2}{\mu }_r{ϵ}_r\left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot p_r\text{d}x.\end{array}$

THEOREME Let $0\le d<N$ , and $s={\alpha }_{N-d}{r}^{N-d}$ . The topological derivative exists if and only if the following limit is satisfied:

$l=\underset{r\to 0}{\text{lim}}\left(l_0\left(r\right)+l_1\left(r\right)\right),$

exists with

$l_0\left(r\right)={\displaystyle {\int }_{{\Omega }_r}\beta {\left|\nabla \left(\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{\sqrt{s}}\right)\right|}^2+\alpha {\left|\frac{u_{{\Omega }_r}-u_{{\Omega }_0}}{\sqrt{s}}\right|}^2}$

and

$\begin{array}{c}{l}_1\left(r\right)=\frac{1}{s}\left[{\displaystyle {\int }_{\partial {\Omega }_r^m\cap \partial E_r}\nabla u_{{\Omega }_0}\cdot \nabla dE{p}_0\text{d}{H}^{N-1}}\right]\text{d}x-\frac{1}{s}{\displaystyle {\int }_{{\Omega }_r}\nabla \left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot \nabla p_r}\\ -\frac{w^2}{c^2}{\mu }_r{ϵ}_r\left(u_{{\Omega }_r}-u_{{\Omega }_0}\right)\cdot p_r\text{d}x.\end{array}$

Moreover, the topological derivative of the function is given by the expression:

$\begin{array}{c}dJ=\underset{r\to 0}{\lim }\frac{J\left({\Omega }_r\right)-J\left(\Omega \right)}{{\alpha }_{N-d}{r}^{N-d}}\\ =l-[{\displaystyle {\int }_E\beta {\left|\nabla u_{{\Omega }_0}\right|}^2}+\alpha {\left|u_{{\Omega }_0}-u_0\right|}^2-\nabla u_{{\Omega }_0}\cdot \nabla p_0\begin{array}{c}\\ \end{array}\\ -\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)p_0+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_0}{p}_0]\text{d}{H}^d.\end{array}$

where $p_0,u_{{\Omega }_0}$ are solutions of systems

${\displaystyle {\int }_{\Omega }\left[2\beta \nabla u_{{\Omega }_0}\cdot \nabla {\varphi }^{\prime }+2\alpha \left(u_{{\Omega }_0}-u_0\right){\varphi }^{\prime }-\nabla {\varphi }^{\prime }\cdot \nabla p_0+\frac{w^2}{c^2}{\mu }_r{ϵ}_r{\varphi }^{\prime }\cdot p_0\right]\text{d}x}=0.$

In particular for $d=0$ ,

$\begin{array}{c}dJ=l-\beta {\left|\nabla u_{{\Omega }_0}\left(x_0\right)\right|}^2-\alpha {\left|u_{{\Omega }_0}\left(x_0\right)-u_0\left(x_0\right)\right|}^2+\nabla u_{{\Omega }_0}\left(x_0\right)\cdot \nabla p_0\left(x_0\right)\\ +\frac{1}{ϵ}{\nabla }_x\rho \cos \left(-\left(kx-wt\right)\right)p_0\left(x_0\right)-\frac{w^2}{c^2}{\mu }_r{ϵ}_r{u}_{{\Omega }_0}\left(x_0\right)p_0\left(x_0\right).\end{array}$

4. Numerical Simulations

In this section, we present the numerical results of the application of topological optimization to one-dimensional and two-dimensional photonic crystal problems. This method seems to be suitable for such kind models [21] [22]. To achieve this, we used the topological gradient descent method to minimize the functional (19), with $\alpha =\beta =1$ .

Indeed, this is an optimization problem of functional J under the constraint that state u is the solution of the partial differential equation (18). The existence and uniqueness of this minimization problem are ensured by the ellipticity of the bilinear form associated with the problem that corresponds to the Helmholtz equation with a source term.

Thus, we used the finite element method to solve the direct and adjoint states involved in the partial differential equation.

We then developed a topological gradient descent algorithm to obtain the optimal shape using FreeFEM ([23]).

The simulations were made under these data the light celerity $c=3\times {10}^8$ ; the pulsation $w=8\times {10}^{-8}$ , the magnetic permeability is assumed to $\mu =1$ ; and the time $T=50$ . The charge distribution follows is supposed to be Gaussian and given by: $\rho \left(x,y\right)=\text{exp}\left(-100\cdot \left({\left(x-0.5\right)}^2+{\left(y-0.5\right)}^2\right)\right)$ in 2D or $\rho \left(x,y,z\right)=\text{exp}\left(-100\cdot \left({\left(x-0.5\right)}^2+{\left(y-0.5\right)}^2+{\left(z-0.5\right)}^2\right)\right)$ in 3D.

4.1. Results in 1D Dimension Photonic Crystal with Charge

We consider an initial domain ${\Omega }_0=\left]-10;10\right[\times \left]-10;10\right[$ . Photonic crystal in one direction are characterized by a medium in which the permittivity varies periodically in one direction. In this section we assume that the permittivity function is periodic in the x direction with a period of 2.

Its expression is $\epsilon \left(x,y\right)=\{\begin{array}{ll}1,\hfill & \text{si}x\mathrm{mod}2<1\hfill \\ 18,\hfill & \text{si}x\mathrm{mod}2\ge 1\hfill \end{array}$ $\forall y\in ℝ$ with $\epsilon \left(x+2,y\right)=\epsilon \left(x,y\right)$ .

For this case, we have Figure 1, which shows the permittivity distribution in both two-dimensional and three-dimensiosnal views. We observed the periodicity of permittivity along the x-direction.

Figure 1. (a) permittivity in 1 dimension view; (b) permittivity in 3 dimensions view.

Using finite elements in FreeFEM and the above data values, Figure 2 represents the direct and adjoint states and the topological derivative with Dirichlet and Neumann conditions.

Figure 2. (a) Direct state; (b) Adjoint state; (c) Dirichlet topological derivative; (d) Neumann topological gradient.

In Figure 2, we observe the distributions of the topological gradient for both Neumann and Dirichlet boundary conditions. These distributions provide insight into the location of the gradient minima by examining the values indicated on the color bar. Using the topological gradient descent algorithm, we identify the points in ${\Omega }_0$ where the gradient reaches a minimum and insert a hole at those locations in order to perturb the topology accordingly.

1) Case of Dirichlet condition around the hole.

Inserting circular holes $B\left(x^{\text{*}},{\epsilon }_0\right)$ in the initial domain where the topological derivatives admits a minimum global at $x^{\text{*}}$ and iterate, we have

Figure 3. (a) Dirichlet topological derivative; (b) After one iteration; (c) After 245 iterations; (d) After 246 iterations.

Graphs (a) and (b) in Figure 3 illustrate the initial step of inserting a hole at the point where the topological gradient reaches its most negative value. Thus, for a Dirichlet boundary condition on the hole, we observe that after 246 iterations, the distribution of the topological gradient is zero everywhere. The descent algorithm converges, and the optimal topology of the domain is achieved, as shown in graph (d).

2) Case of a Neumann condition around the hole.

For the Neumann boundary condition, graphs (a) and (b) of Figure 4 illustrate the creation of a hole at the location where the topological gradient is the most negative. This leads to a decrease in the gradient distribution, which progressively approaches zero throughout the domain $\Omega$ .

Figure 4. (a) Neumann topological derivative after one iteration; (b) After two iterations; (c) After 443 iterations; (d) After 444 iterations. After 444 iterations, the algorithm converges. The topological gradient becomes nearly zero throughout the domain, indicating that the optimal shape minimizing the objective function has been reached.

4.2. Results in 2D Dimension Photonic Crystal with Charge and Height Near to Zero

Considering a initial domain ${\Omega }_0=\left]-1;1\right[\times \left]-1;1\right[$ .

In this section we present a quasi two dimensional photonic crystal which is characterized by a periodic permittivity $\epsilon \left(x,y\right)$ in the x and y directions.

Thus the permittivity is shown by Figure 5

Figure 5. (a) Permittivity periodicity in 2 dimension; (b) Permittivity in 3 dimension view.

With this characteristic in the charged photonic medium, we have:

Figure 6. (a) Direct state; (b) Adjoint state; (c) Dirichlet topological derivative; (d) Neumann topological gradient.

For two-dimensional photonic crystals with charge sources and unit thickness, Figure 6 shows the topological gradient plots under both Dirichlet and Neumann boundary conditions, as well as the corresponding direct and adjoint states. From these results, we observe that the location of the most negative values of the topological gradient indicates where the structure should be perturbed to improve the objective functional. The direct state reflects the wave propagation in the current structure, while the adjoint state captures the sensitivity of the objective to changes in the domain. The comparison between the two boundary conditions also highlights the influence of the choice of boundary model on the resulting optimized design.

1) Case of Dirichlet condition around the hole. Inserting circular hole $B\left(x^{\text{*}},{\epsilon }_0\right)$ in the initial domain where the topological derivatives admits a minimum global at $x^{\text{*}}$ and iterating we have the Figure 7 which shows the variation of the gradient when perturbing the topology of the initial domain.

Figure 7. (a) Dirichlet topological derivative after 1 iteration; (b) After 2 iteration; (c) After 88 iterations; (d) After 89 iterations.

Figure 7 shows that insertion a hole in the minium of the Dirichlet gradient vanishes the gradient sited at neighborhood at this minimizer. Following the procedure of the algorithm, we see after 89 iterations or insertions of holes, the gradient vanishes everywhere at the domain. We have convergence and the topological optimal design in graph (d) of Figure 7.

2) Case of a Neumann condition around the hole.

For the Neumann boundary condition, after inserting one or two holes at the points where the topological derivative reaches its minimum, we observe that the gradient distribution approaches zero.

Figure 8 shows that after 1241 iterations we get the optimal design.

Figure 8. (a) Neumann topological derivative after one iteration; (b) After two iterations; (c) After 1240 iterations; (d) After 1241 iterations.

5. Conclusions

The present work focuses on the topological optimization of one- and two-dimensional charged photonic crystal problems. We first derive the topological derivative for Dirichlet boundary conditions, followed by the case of Neumann boundary conditions. Numerical simulations are then performed and interpreted for both types of boundary conditions. The application of the topological derivative to charged photonic crystals appears to be well-suited, based on the convergence results obtained for both Dirichlet and Neumann cases.

The results show that, for the same initial domain and simulation parameters, the descent-gradient topological algorithm converges more rapidly under Dirichlet conditions than under Neumann conditions.

This observation provides a useful recommendation for choosing boundary conditions when aiming to efficiently design photonic crystals and save computational time.

In future work, we aim to couple shape and topological derivatives in the context of photonic slabs. This combined algorithm will enable the search for an optimal topology without the need to add or remove material along the domain boundary.

Finally, a comparative study between our proposed algorithm and one based on the coupling of shape and topological gradients would be essential to evaluate the practical relevance of these tools in mathematical and physical modeling.