Unidirectional invisibility in a PT-symmetric structure designed by topology optimization

Kei Matsushima; Yuki Noguchi; Takayuki Yamada

doi:10.1364/OL.460488

In recent years, parity–time (PT) symmetry has attracted significant attention in optics and photonics. After El-Ganainy et al. [1] theoretically demonstrated that non-Hermitian Hamiltonians have a real eigenvalue spectrum in PT-symmetric optical systems, many studies on theoretical aspects of PT symmetry, especially exceptional points [2], have been conducted, as well as experimental demonstrations [3].

Unidirectional invisibility is one of the most interesting phenomena in PT-symmetric systems. Without gain or loss, transmission and reflection must be invariant under the inversion of incident directions because of reciprocity, meaning that perfect invisibility (transmittance $T=1$ and reflectance $R=0$) should be observed for both incident directions. Therefore, reciprocity or energy conservation should be broken to achieve unidirectionality. PT symmetry allows us to break energy conservation while balancing the amount of loss and gain for a specific incident wave. Unidirectional invisibility is potentially important for designing innovative optical network systems [4].

Although unidirectional invisibility has frequently been discussed in one-dimensional two-port systems [5–9], some recent works have found that an analogous phenomenon can be observed in two-dimensional systems with infinite scattering channels. For example, Zhu et al. [10] employed transformation optics to determine a PT-symmetric material distribution to achieve the unidirectional invisibility of a perfect electric conductor (PEC). Sounas et al. [11] designed PT-symmetric metasurfaces that cloak an inside object for a specific direction. However, these designs involve continuously inhomogeneous materials, which are not always suitable for fabrication. See Ref. [12] for an outlook on the challenges involved in the design of passive and active metamaterials for invisibility cloaks.

In the present study, we demonstrate that two-dimensional unidirectional invisibility can be achieved using a PT-symmetric structure with a piecewise-constant material distribution. We use a level-set-based topology optimization algorithm [13–17] to design the cloaking structure. Topology optimization enables the determination of a PT-symmetric material design that minimizes the total scattering cross section for a given incident wave. First, we describe the cloaking model and formulation. As shown in Fig. 1, we aim to achieve unidirectional invisibility of a cylindrical PEC $\Omega _\mathrm {PEC}$ of radius $R$ against a TE-polarized plane incident wave by distributing dielectric materials with relative permittivity $\varepsilon$ and permeability $\mu =1$. We formulate the scattering problem as

(1)$$\nabla^2 u + \frac{\varepsilon \omega^2}{c^2} u = 0 \quad \text{in }\mathbb{R}^2\setminus\overline{\Omega_\mathrm{PEC}},$$

(2)$$\frac{\partial u}{\partial n}:= \boldsymbol{n}\cdot \nabla u = 0 \quad \text{on }\partial\Omega_\mathrm{PEC},$$

(3)$$\text{Sommerfeld radiation condition for} \quad u-\exp({\mathrm{i} k \boldsymbol{p}^\mathrm{in}\cdot \boldsymbol{x}}),$$

where $u$ represents the transverse component of the electric field, $\boldsymbol {n}$ is the unit normal vector, and the overline denotes the closure. It is straightforward to replace the PEC with a dielectric material. The wavenumber of the incident plane wave is given by $k=\omega /c$, where $c$ is the speed of light in vacuum and $\omega$ is the angular frequency. The propagation direction of the plane wave is denoted by the unit vector $\boldsymbol {p}^\mathrm {in}\in \mathbb {R}^2$.

Fig. 1. Level set function $\phi$ represents a quarter of a PT-symmetric structure. The material region is defined as $\{ \boldsymbol {x} \mid \phi (\boldsymbol {x})<0 \}$ in the design domain.

Download Full Size | PDF

We solve the scattering problem using a mode-matching finite element method [18], implemented in the open-source software FreeFem++ [19] with the Mmg platform [20]. To this end, we rewrite the boundary value problem as the following variational formulation:

(4)$$0 = \int_{D\setminus\overline{\Omega}_\mathrm{PEC}} \nabla \tilde{u} \cdot {\nabla} u \mathrm{d}\Omega -\frac{\omega^2}{c^2} \int_{D\setminus\overline{\Omega}_\mathrm{PEC}} \varepsilon \tilde{u}u \mathrm{d}\Omega - \int_{\partial D}\frac{\partial u}{\partial n}\tilde{u} \mathrm{d}\Gamma,$$

for all test functions $\tilde {u}$, where $D$ is a disk centered at the origin. The radius of $D$ is sufficiently large so that the disk $D$ encloses all the inhomogeneity $\varepsilon \neq 1$. The variational Eq. (4) is discretized using piecewise quadratic polynomials on a triangular mesh on $D$ and the Galerkin method. In the exterior of $D$, we write a solution as the following multipole expansion with coefficients $A_m$ and $B_m$:

(5)$$\begin{aligned} u(\boldsymbol{x}) & = \sum_{m={-}\infty}^\infty A_m J_m(k|\boldsymbol{x}|)\exp(\mathrm{i}m\theta(\boldsymbol{x}))\\ & + \sum_{m={-}\infty}^\infty B_m H^{(1)}_m(k|\boldsymbol{x}|)\exp(\mathrm{i}m\theta(\boldsymbol{x})) \quad \boldsymbol{x}\in\mathbb{R}^2\setminus\overline{D}, \end{aligned}$$

where $J_m$ and $H^{(1)}_m$ represent the $n$th-order Bessel and Hankel function of the first kind, respectively, and $\theta (\boldsymbol {x})$ denotes the angle $\theta (\boldsymbol {x})=\tan ^{-1}(x_2/x_1)$. We impose the weak continuity of $u$ and its normal derivative to couple the interior polynomial approximation with exterior multipole expansion as follows:

(6)$$0 = \int_{\partial D}\left(u|_+{-} u|_- \right)\exp(-\mathrm{i}m\theta(\boldsymbol{x})) \mathrm{d}\Gamma \quad \text{for all }m,$$

(7)$$0 = \int_{\partial D}\left(\frac{\partial u}{\partial n}\bigg|_+{-} \frac{\partial u}{\partial n}\bigg|_- \right) N_i \mathrm{d}\Gamma \quad \text{for all }i,$$

where $N_i$ represent the quadratic shape functions on the triangular mesh and $u|_\pm$ denote the trace from the exterior and interior to $\partial D$, respectively. Combining Eqs. (4)–(7), the unknown coefficients $B_m$ and nodal values on $D$ for a given incident wave $A_m$ can be obtained.

It is well known that a solution to the two-dimensional scattering problem allows the following asymptotic expression:

(8)$$\begin{aligned}u(\boldsymbol{x}) = \exp({\mathrm{i}k\boldsymbol{p}^\mathrm{in}\cdot \boldsymbol{x}}) - \mathrm{i} \frac{\exp({\mathrm{i}k|\boldsymbol{x}|}-\pi/4)}{\sqrt{|\boldsymbol{x}|}} A_\infty(\boldsymbol{x}/|\boldsymbol{x}|,\boldsymbol{p}^\mathrm{in})\\ + o(|\boldsymbol{x}|^{{-}1/2}) \quad |\boldsymbol{x}|\to\infty, \end{aligned}$$

where $A_\infty$ is called a scattering amplitude. For a given incident plane wave, invisibility is achieved when the following scattering cross section vanishes:

(9)$$\sigma(\boldsymbol{p}^\mathrm{in}) = \int_{|\hat{\boldsymbol{x}}|=1} |A_\infty(\hat{\boldsymbol{x}},\boldsymbol{p}^\mathrm{in})|^2 \mathrm{d}\hat{\boldsymbol{x}} \simeq \frac{2\pi}{N} \sum_{i=1}^N |A_\infty(\boldsymbol{p}^i,\boldsymbol{p}^\mathrm{in})|^2,$$

with $\boldsymbol {p}^i=(\cos (2\pi i/N),\sin (2\pi i/N))^T$ and sufficiently large integer $N$. We say that the system exhibits unidirectional invisibility if $\sigma (\boldsymbol {p}^\mathrm {in}) = 0$ and $\sigma (-\boldsymbol {p}^\mathrm {in}) > 0$.

The most important property of $A_\infty$ is the reciprocity, i.e., $A_\infty (\boldsymbol {p}^\mathrm {sc},\boldsymbol {p}^\mathrm {in}) = A_\infty (-\boldsymbol {p}^\mathrm {in},-\boldsymbol {p}^\mathrm {sc})$ for all unit vectors $\boldsymbol {p}^\mathrm {in},\boldsymbol {p}^\mathrm {sc}\in \mathbb {R}^2$. In addition, if the medium has neither gain nor loss, we have the following optical theorem [21]:

(10)$$\sigma(\boldsymbol{p}^\mathrm{in}) = \mathrm{Im}\,\left[ -\sqrt{\frac{8\pi c}{\omega}} A_\infty(\boldsymbol{p}^\mathrm{in},\boldsymbol{p}^\mathrm{in}) \right].$$

The reciprocity and optical theorem state that unidirectional invisibility cannot be achieved without gain or loss because

(11)$$\begin{aligned}\sigma(-\boldsymbol{p}^\mathrm{in}) & = \mathrm{Im}\,\left[ -\sqrt{\frac{8\pi c}{\omega}} A_\infty(-\boldsymbol{p}^\mathrm{in},-\boldsymbol{p}^\mathrm{in}) \right]\\ & = \mathrm{Im}\,\left[ -\sqrt{\frac{8\pi c}{\omega}} A_\infty(\boldsymbol{p}^\mathrm{in},\boldsymbol{p}^\mathrm{in}) \right] = \sigma(\boldsymbol{p}^\mathrm{in}). \end{aligned}$$

To break the scattering symmetry, we introduce material loss and gain by letting $\varepsilon$ be a complex number with a nonzero imaginary part. Assuming the time-harmonic factor $\exp ({-\mathrm {i}\omega t})$ with time $t$, we express the material loss by a positive imaginary part of $\varepsilon$. Similarly, a negative imaginary part denotes material gain.

We assume that the distribution of $\varepsilon$ satisfies the PT symmetry $\varepsilon (-\boldsymbol {x})=\bar {\varepsilon }(\boldsymbol {x})$. Such two-dimensional PT-symmetric structures have been discussed in, for example, Refs. [10,11,22]. Although we do not guarantee that the PT symmetry is strictly necessary for the unidirectionality, we follow previous studies [10,11] and limit ourselves to the PT-symmetric case. The distribution of $\varepsilon$ is determined by topology optimization. As shown in Fig. 1, we design a quarter of the PT-symmetric cloaking structure. In this design domain, we introduce a scalar function $\phi$, called a level set function, and define the material region as $\{ \boldsymbol {x} \mid \phi (\boldsymbol {x})<0 \}$ (i.e., $\phi (\boldsymbol {x})>0$ represents a vacuum). We seek a distribution of $\phi$ that minimizes the total scattering cross section $\sigma$ for the right-propagating plane wave. To this end, we iteratively update the level set function $\phi = \phi _i$, where the index $i$ denotes the current optimization step, using the following formula [23]:

(12)$$\phi_i = \phi_{i-1} + \Delta_i \left( \mathcal{T}({\phi_{i-1}}) + \tau \nabla^2 \phi_i \right),$$

where $\Delta _i>0$ denotes the step size, $\tau >0$ represents a parameter for controlling geometrical complexity, and $\mathcal {T}$ is a topological derivative of $\sigma$. The topological derivative represents a sensitivity of $\sigma$ when a small inhomogeneity $\varepsilon \neq 1$ appears at the homogeneous background [24]. We used the adjoint variable method to obtain the following explicit form of $\mathcal {T}$:

(13)$$\mathcal{T} = \frac{2\pi}{N} \sum_{i=1}^N \mathrm{Re}\,\left[\overline{A_\infty(\boldsymbol{p}^i,\boldsymbol{p}^\mathrm{in})} \delta A_\infty(\boldsymbol{p}^i,\boldsymbol{p}^\mathrm{in}) \right],$$

(14)$$\delta A_\infty(\boldsymbol{p}^i,\boldsymbol{p}^\mathrm{in}) = k^2 u(\boldsymbol{x};-\boldsymbol{p}^i)u(\boldsymbol{x};\boldsymbol{p}^\mathrm{in}) \delta \varepsilon,$$

where $\delta \varepsilon$ denotes the perturbation of $\varepsilon$ at $\boldsymbol {x}$ and $u(\boldsymbol {x};\boldsymbol {p})$ denotes the solution $u$ to the scattering problem (1)–(3) with $\boldsymbol {p}^\mathrm {in}=\boldsymbol {p}$. For details, see Refs. [13,25]. We targeted the incident direction $\boldsymbol {p}^\mathrm {in}=(1,0)^T$ with an operating frequency $\omega R/(2\pi c) = 0.75$ and performed level-set-based topology optimization to minimize $\sigma (\boldsymbol {p}^\mathrm {in})$. The lossy dielectric material is characterized by $\varepsilon =2{+}0.05\mathrm {i}$. We chose the initial level set function as $\phi _0=-\mathcal {T}_0$, where $\mathcal {T}_0$ is the topological derivative in the homogeneous background $\varepsilon =1$ [13]. During the optimization process, we plotted the objective values $\sigma$, normalized by the scattering cross section $\sigma _\mathrm {bare}$ of the bare PEC, as shown in Fig. 2. The objective values successfully decreased by almost 99% compared to the bare object until convergence.

Fig. 2. Convergence history of the objective function $\sigma$.

Download Full Size | PDF

The optimized PT-symmetric structure, shown in Fig. 3, has a striped pattern with the interval being approximately half the operating wavelength $\lambda$. We illuminated the right- and left-propagating plane waves and computed their scattering, as shown in Fig. 4. The results show that the designed PT-symmetric structure significantly suppressed the scattering of the right-propagating plane wave, whereas invisibility does not hold for the left-propagating one. Unidirectionality is also confirmed by computing the scattering amplitudes $A_\infty$. Figure 5 shows the amplitudes $|A_\infty (\boldsymbol {p}^\mathrm {sc},\boldsymbol {p}^\mathrm {in})|$ for each direction $\boldsymbol {p}^\mathrm {sc}=(\cos \theta,\sin \theta )^T$ with $0\leq \theta \leq 2\pi$. The results indicate that the designed PT-symmetric structure exhibits less scattering intensity than the bare PEC for the right-propagating incident wave, whereas strong scattering is observed for the left-propagating one. Although the scattering cross section $\sigma$ depends on the incident direction, the forward scattering intensity $A_\infty (\boldsymbol {p}^\mathrm {in},\boldsymbol {p}^\mathrm {in})$ is consistent, even if the incident direction is inverted, indicating that the PT-symmetric structure still exhibits the reciprocity $A_\infty (\boldsymbol {p}^\mathrm {sc},\boldsymbol {p}^\mathrm {in}) = A_\infty (-\boldsymbol {p}^\mathrm {in},-\boldsymbol {p}^\mathrm {sc})$.

Fig. 3. Designed PT-symmetric cloaking structure.

Download Full Size | PDF

Fig. 4. Total field $u$ when (a) right-propagating and (b) left-propagating plane waves (unit amplitude) impinge onto the designed structure.

Download Full Size | PDF

Fig. 5. Scattering amplitudes $|A_\infty |$ for (a) right-propagating and (b) left-propagating plane waves with unit amplitude. The solid and dashed lines represent the amplitudes for the PT-symmetric cloaking structure and bare PEC, respectively.

Download Full Size | PDF

To investigate scattering asymmetry, we computed the spectrum of $\sigma$ for both the right- and left-propagating plane waves. The results are presented in Fig. 6. For the right-propagating plane wave, the spectrum exhibits a rapid decline at the target frequency, with a mimimum value of $\sigma /\sigma _\mathrm {bare}=9.6\times 10^{-3}$. Figure 7 shows the incident-angle dependence of $\sigma$ at the target frequency. The cross section $\sigma$ attains its minimum value at an incident angle $\theta ^\mathrm {in}=0^{\circ }$ (right-propagating case), whereas the left-propagating wave has $\sigma /\sigma _\mathrm {bare}=0.80$, which is 83 times larger than the minimum value. Since the designed structure has gain and loss, we computed the associated absorption cross section:

(15)$${ \sigma_\mathrm{loss} ={-}\frac{1}{k}\mathrm{Im}\,\left[ \int_{\partial D} \bar{u}\frac{\partial u}{\partial n}\mathrm{d}\Gamma \right] }$$

and plotted the values in Fig. 8. For the left-propagating case, the absorption cross section is negative, i.e., the overall system produces some energy around the target frequency. On the other hand, the right-propagating wave balances the amount of loss and gain at the target frequency, which is essential for invisibility.

Fig. 6. Spectrum of the scattering cross section of the designed structure for right- and left-propagating plane waves.

Download Full Size | PDF

Fig. 7. Scattering cross section $\sigma$ for various incident angles $\theta ^\mathrm {in}$.

Download Full Size | PDF

Fig. 8. Spectrum of the absorption cross section for right- and left-propagating plane waves with unit amplitude.

Download Full Size | PDF

In conclusion, we proposed a new PT-symmetric design, to the best of our knowledge, for unidirectional invisibility of a PEC in two dimensions and showed that a piecewise homogeneous dielectric structure can exhibit unidirectionality. The cloaking structure was designed via topology optimization. We conducted numerical experiments and confirmed that the designed structure suppresses the total scattering cross section by almost 99% unidirectionally without breaking the system’s reciprocity. Since the designed structure works only in a specific narrow band, an open question is whether broadband unidirectional invisibility can be realized with a piecewise homogeneous material. Other future work is the experimental realization of the designed PT-symmetric system. The proposed design methodology will bring a new insight to the understanding of non-Hermitian physics and a technology to fabricate PT-symmetric structures.

Funding

Japan Society for the Promotion of Science (JP22K14166).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

REFERENCES

1. R. El-Ganainy, K. G. Makris, D. N. Christodoulides, and Z. H. Musslimani, Opt. Lett. 32, 2632 (2007). [CrossRef]

2. M.-A. Miri and A. Alù, Science 363, eaar7709 (2019). [CrossRef]

3. C. E. Rüter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, and D. Kip, Nat. Phys. 6, 192 (2010). [CrossRef]

4. A. Regensburger, C. Bersch, M.-A. Miri, G. Onishchukov, D. N. Christodoulides, and U. Peschel, Nature 488, 167 (2012). [CrossRef]

5. Z. Lin, H. Ramezani, T. Eichelkraut, T. Kottos, H. Cao, and D. N. Christodoulides, Phys. Rev. Lett. 106, 213901 (2011). [CrossRef]

6. L. Feng, Y. L. Xu, W. S. Fegadolli, M. H. Lu, J. E. Oliveira, V. R. Almeida, Y. F. Chen, and A. Scherer, Nat. Mater. 12, 108 (2013). [CrossRef]

7. Y. Fu, Y. Xu, and H. Chen, Opt. Express 24, 1648 (2016). [CrossRef]

8. Y. Fu, X. Zhang, Y. Xu, and H. Chen, J. Appl. Phys. 121, 094503 (2017). [CrossRef]

9. L. Yuan and Y. Y. Lu, Phys. Rev. A 100, 053805 (2019). [CrossRef]

10. X. Zhu, L. Feng, P. Zhang, X. Yin, and X. Zhang, Opt. Lett. 38, 2821 (2013). [CrossRef]

11. D. L. Sounas, R. Fleury, and A. Alù, Phys. Rev. Appl. 4, 014005 (2015). [CrossRef]

12. R. Fleury, F. Monticone, and A. Alù, Phys. Rev. Appl. 4, 037001 (2015). [CrossRef]

13. K. Nakamoto, H. Isakari, T. Takahashi, and T. Matsumoto, Mech. Eng. J. 4, 16-00268 (2017). [CrossRef]

14. G. Fujii, H. Watanabe, T. Yamada, T. Ueta, and M. Mizuno, Appl. Phys. Lett. 102, 251106 (2013). [CrossRef]

15. T. Yamada, H. Watanabe, G. Fujii, and T. Matsumoto, IEEE Trans. Magn. 49, 2073 (2013). [CrossRef]

16. G. Fujii and T. Ueta, Phys. Rev. E 94, 043301 (2016). [CrossRef]

17. N. Kishimoto, K. Izui, S. Nishiwaki, and T. Yamada, Appl. Phys. Lett. 110, 201104 (2017). [CrossRef]

18. R. J. Astley, Commun. Numer. Meth. Engng. 12, 257 (1996). [CrossRef]

19. F. Hecht, Journal of Numerical Mathematics 20, 251 (2012). [CrossRef]

20. C. Dobrzynski, ’Mmg3d: User guide’, (2012).

21. H. Ammari, B. Fitzpatrick, H. Kang, M. Ruiz, S. Yu, and H. Zhang, Mathematical and Computational Methods in Photonics and Phononics (American Mathematical Society, 2018).

22. P. Ambichl, K. G. Makris, L. Ge, Y. Chong, A. D. Stone, and S. Rotter, Phys. Rev. X 3, 041030 (2013). [CrossRef]

23. T. Yamada, K. Izui, S. Nishiwaki, and A. Takezawa, Comput. Methods Appl. Mech. Eng. 199, 2876 (2010). [CrossRef]

24. J. Sokolowski and J. P. Zolésio, Introduction to shape optimization (Springer, 1992).

25. H. Ammari and H. Kang, J. Math. Anal. Appl. 296, 190 (2004). [CrossRef]

Unidirectional invisibility in a PT-symmetric structure designed by topology optimization

Abstract

Funding

Disclosures

Data availability

REFERENCES

Data availability

Cited By

Figures (8)

Equations (15)

Optics Letters