Introduction

Topological insulators (TIs) have been extensively studied in condensed matter physics1,2,3 but are also moving forward at fast pace in classical systems such as photonics4,5,6,7, acoustics8,9,10,11,12,13,14, and mechanics15,16,17,18,19. A key feature of the topological insulators is the existence of robust edge states that are immune to backscattering from disorder or interface variations. Conventionally, the edge states were predicted by the bulk-edge correspondence principle, but recently, the concept of TIs has been generalized to higher-order topological insulators (HOTIs), at which the conventional bulk-boundary correspondence does not apply20,21,22,23,24. Beyond theoretical advances, classical lattices have shown to be a fruitful base since numerous experimental verifications have seen the light of day25 including photonic26,27,28 and phononic crystal systems29,30,31,32,33,34,35,36.

Parity-time \({\mathcal{P}}{\mathcal{T}}\) symmetry on the other hand, describes the invariance of a non-Hermitian system that can have real eigenvalue spectra despite its complex entities. Those systems are a special kind of physical configurations invariant upon the combined parity \({\mathcal{P}}\) and time reversal \({\mathcal{T}}\) operations, which have been realized in artificial structures hosting balanced gain and loss constituents. Systems with \({\mathcal{P}}{\mathcal{T}}\)-symmetry have been extensively studied in optics and acoustics enabling unusual wave propagation characteristics such as unidirectional invisibility and focusing, as well as complex topological valley-Hall effects37,38,39,40,41,42.

In this article, we use a multiple scattering theory (MST) to calculate the complex properties of \({\mathcal{P}}{\mathcal{T}}\)-symmetric second-order topological sonic crystals. Generally speaking, non-Hermitian first-order TIs also have a bulk-boundary correspondence that however has to be modified significantly43,44. Most recently, non-Hermitian HOTIs, as a whole, have been reported owning an even more complex bulk-boundary correspondence45,46,47.

This paper is a follow up to our short letter48, in which all results have been numerically obtained by a commercial finite element method (FEM) solver. Beyond the excellent agreement between the two techniques permitting us to study topological corner, edge and bulk states, as well as their spatial pressure profiles, we deliberately employ the MST to elaborate the robustness of non-Hermitian corner states in the presence of interstitial cylinder-defects.

Results

Multiple scattering theory

We consider a 2D arrangement of N-cylinders located at Rα, where α = 1, 2...N. These cylinders, made of fluid A with density ρa and speed of sound ca, are embedded in a fluid background with density ρb and speed of sound cb. If an external field Pin, with angular frequency ω, impinges the cluster of cylinders, the total scattered field, at a given position r in polar coordinates (r, θ), will be given by:

$$P^{{\mathrm{sc}}}(r,\theta ) = \mathop {\sum}\limits_{\alpha = 1}^N {\mathop {\sum}\limits_q {(A_{{\alpha}} )_q} } H_q(kr_{{\alpha}} )e^{{{i}}q\theta _{{\alpha}} },$$
(1)

where Hq is the qth order Hankel function of the first kind, (rα, θα) are the polar coordinates of the αth cylinder in the reference frame rα = r − Rα, k = ω/cb is the wavenumber and (Aα)s are the coefficients of the scattered field to be determined. The total incident field at the αth cylinder can be expressed either as a linear combination of Bessel functions of the qth order (in a more general approach),

$$P_{\mathrm{in}}({{r}}_\alpha ,{{\theta}} _{{\alpha}} ) = \mathop {\sum}\limits_{{q}} {(B_\alpha )_{{q}}} J_{{q}}(kr_\alpha )e^{{{i}}q\theta _{\mathrm{\alpha}} },$$
(2)

or as a sum of the external field Pin and the field scattered by all cylinders except the αth one. If we rewrite the latter on a α basis by means of the Graft’s addition theorem49, we obtain the following expression:

$$\begin{array}{l}P_\alpha ^{{\mathrm{in}}}(r_\alpha ,\theta _\alpha ) \hskip -0pt = \mathop {\sum}\limits_q {(A_\alpha ^0)_q} J_q(kr_\alpha )e^{iq\theta _\alpha } \\ + \mathop {\sum}\limits_s {\mathop {\sum}\limits_q {\mathop {\sum}\limits_{\beta \ne \alpha } {(A_\beta )_s} } } H_{q - s}(kR_{\alpha \beta })e^{i(s - q)\theta _{\alpha \beta }}J_q(kr_\alpha )e^{iq\theta _\alpha }.\end{array}$$
(3)

Reformulating the equations into matrix form, coefficients (Bα)q and (Aα)q of the incident and scattered field, respectively, are related by means of the T-matrix \((A_\alpha )_q = \mathop {\sum}\nolimits_s {(T_\alpha )_{qs}} (B_\alpha )_s\). Thus, once we determine the value of the T-matrix, as well as setting the amplitude \((A_\alpha ^0)_q\) of the incident field, we are able to compute the scattered amplitude (Aα)q, by solving the system of linear equations, resulting from substituting Eq. (2) into Eq. (3) and introducing the value of the T-matrix therein. The structure that we are going to study is made of fluid cylinders of circular cross section of radius Ra, whose T-matrix reads:

$$T_q = - \frac{{\rho _{q}J_{q}^{\prime} (kR_a) - J_q(kR_a)}}{{\rho _qH_q^\prime (kR_a) - H_q(kR_a)}} \\ \rho _q = \frac{{\rho _ac_aJ_q(k_aR_a)}}{{\rho _bc_bJ_q^\prime (k_aR_a)}}.$$
(4)

The fluid cylinders are irradiated by an external point source \(P_{{\mathrm{in}}}(r,\theta ) = \mathop {\sum}\nolimits_s {H_s} (kr_s)e^{is\theta _s}\) located at Rs = r − rs. Rewriting this expression on an α basis allows one to determine the amplitude of the incident field as follows system, \((A_q^0) = H_{s - q}(kR_s^\prime )e^{i(s - q)\theta _s^\prime }\).

Corner states

We begin with a 2D regular Hermitian square sonic crystal with lattice constant a, whose unit cell is composed of four isotropic fluid cylinders embedded in air with an adjacent center-to-center distance D/a = 0.5. When adjusting the ratio from D/a < 0.5 (shrunk) to D/a > 0.5 (expanded), a topological phase transition characterized by the vectorial Zak phase P emerge. This phase is defined as an integration of the non-Abelian Berry connection Amn over the first Brillouin Zone (BZ):

$${\boldsymbol{P}} = \mathop {\iint}\nolimits_{1^{{\mathrm{st}}}{\mathrm{BZ}}} {{\mathrm{d}}k_1{\mathrm{d}}k_2{\mathrm{Tr}}({{\bf{A}}}),} \quad {{\bf{A}}}_{{\bf{mn}}} = i\langle u_m|\partial _{\mathbf{k}}|u_n\rangle$$
(5)

where |um〉 is the Bloch function for the mth band and m, n run over occupied bands. Due to the C4v symmetry, the previous expression can be rewritten as50,51

$$P_i = \pi \left( {\mathop {\sum}\limits_n {q_i^n} \,{\mathrm{mod}}2} \right), \quad ( - 1)^{q_i^n} = \frac{{\eta _n(X_i)}}{{\eta _n({\it{\Gamma }})}}$$
(6)

where i represents the directions x or y, Xi denotes the X(Y) point of the first BZ and ηn stands for the parity of the nth band at the high symmetry point. In a physical context, the Zak phase is related to 2D wave polarization (σi with i = x, y) through σi = Pi/2π. When the origin is chosen to be the center of the metamolecule the Zak phase has only two possible values, 0 or π, with corresponding polarization quantized as 0 and 1/252,53. In the shrunk scenario, when D/a < 0.5 the 2D Zak phase equals (0, 0) with σ = (0, 0), indicating a trivial state. In contrary, when the metamolecule is expanded beyond D/a > 0.5 the Zak phase equals (π, π) with \(\sigma = (\frac{1}{2},\frac{1}{2})\), indicating a non-trivial state. The coexistence of the non-trivial dipole moments in both directions generate a corner confining state according to the edge-corner correspondence22,48,54. As a result, a corner state is sustained at the 90-degree corners between the trivial and non-trivial regions having non-zero Zak phases in both directions along the interfaces.

Resting on the aforementioned fundamental theory, we construct a finite concentric square crystal (CSC) shown in Fig. 1a, which is based on an inner sonic crystal of size 8a × 8a constituting the topological non-trivial region with D/a = 0.714 that is enclosed by a topological trivial region of width 5a and D/a = 0.280. As simulations show, indeed we are able to select other ranges of geometrical parameters (see Supplementary Note 1). The lattice constant of the both sonic crystals are taken to be a = 1 cm and the radii of the cylinders are r = 0.12a. The mass density and the bulk modulus of the background correspond to air, ρ0 = 1.21 kg/m3 and K0 = 1.4 ⋅ 105 Pa, whereas the mass density and the bulk modulus of the cylinders are ρ = ρ0 and K = K0/9, respectively. Our aim is to employ the MST to obtain the topological corner states of a two-dimensional second-order topological insulator. In doing this, we first compare the eigenfrequencies calculated with FEMs with the pressure spectrum peaks calculated with a MST. Later one, the MST is also employed to compute pressure field distribution of the corner states. Within the spectral region of the bulk band gap of the periodic crystal48, by means of FEM simulations we compute the eigenfrequencies at their corresponding state number of the Hermitian CSC as represented by dots in Fig. 1b. Specifically, we characterize these dots as bulk (black), edge (green) and corner states (red). When conducting the numerical experiment based on the MST we place an acoustic point source at a relevant position of the CSC to computationally be able to measure spectrally the response of the specific confined states. In other words, as corner states confine in the corner, we place the source and probe its excitation near one of the four corners. Similarly, as schematically represented by dots (sources) and stars (measurement point) in Fig. 1a, we excite and detect the topological edge and bulk states near their origins of the CSC. The spectral pressure responses obtained by the MST and the topological eingestates from FEM simulations shown in Fig. 1b display a remarkably good agreement. The consistency between the spectral locations obtained by the two distinct methods permits us to faithfully conduct a robust study solely employed by semi-analytical means. Before that, however, we focus on the spatial characteristics of the corner states shown in Fig. 1c. The MST pressure maps visualize equivalent intensity confinement in all of the corners associated to the four possible corner states computed in the spectrum, Fig. 1b. In addition, our predictions display some leakage at the corners enabling sound to penetrate at certain distance into the trivial region. By adjusting its width, the topological confinement can be controlled. However, increasing it would also increase the number of cylinders included in the structure. With a width of 5a the total amount of cylinders forming the CSC is 1296, already setting high demands in terms of computation time.

Fig. 1
figure 1

Hermitian concentric square crystal. a Schematic representation of the crystal made of 8a × 8a metamolecule arrays defining the non-trivial region with D/a = 0.714 enclosed by a trivial region with D/a = 0.280 and thickness 5a. Circles represent the position of the source exciting the bulk (black), the edge (green) and the corner states (red), which correspondingly are indicated by stars. b Comparison between the real eigenfrequencies obtained by finite element method simulations (dots) with the normalized pressure of the excited states from the multiple scattering theory computation. c Pressure maps of all corner states. Inserts illustrate the excitation point

Full size image

We now study the evolution of the corner states while a total amount of nine interstitial fluid cylindrical defects are gradually introduced. As can be seen in the inset of Fig. 2a, the cylindrical defects are symmetrically located near the upper interface between the non-trivial and the trivial region, and their radii r = 0.15a are slightly larger than the radii of the cylinders forming the CSC, although we chose the same acoustic parameters. As the intensity is the same of all the corner states, a complete study is done by only analyzing the evolution of one of them. Thus, the source is placed in the center of the CSC (red dot in Fig. 2a) and the detection point (red star) is chosen to be located at the upper-right corner between the trivial and the non trivial region. Corner states are exited in narrow frequency window centered at 12.50 kHz as displayed in Fig. 2b. The results show a slight pressure increase when adding defect up to ndefects = 3. Beyond this, a smooth frequency shift together with a slight pressure reduction is predicted when ndefects is further increased. For some selected number of interstitial defects as labeled and highlighted by white dots in Fig. 2c, we plot the spatial pressure field maps at their corresponding spectral peaks, clearly displaying the resilience of the corner states in response to a point-like excitation.

Fig. 2
figure 2

Study of the topological protection of corner states against perturbations. a Schematic representation of the locations of the interstitial defects embedded into the non trivial part of the Hermitian concentric square crystal. b Spectral dependence of the pressure at the upper-right corner when the source is placed in the crystal center with a steadily increased number of defects. c Mapping the pressure distribution of the corner states in the presence of 0, 3, and 7 defects. In all simulations the source has been placed in the crystal center

Full size image

We now extend the study of HOTIs to non-Hermitian CSCs by introducing acoustic gain and loss into the metamolcules by employing active acoustic metamaterials. Thus, the complex mass densities of the cylinders are defined as ρgain = (1 + iβ)ρ0 for the gain components and ρloss = (1 − iβ)ρ0 for the lossy ones, while the bulk modulus remain the same as in the Hermitian case, K = K0/9. In order to have an equal amount of loss and gain as measured by the non-Hermiticity parameter β, we construct two types of metamolecules, each having two gain and two loss cylinders. The first non-Hermitian arrangement of the CSC introduces symmetrical gain and loss units diagonally distributed as illustrated in Fig. 3a. Notice that gain (loss) cylinders are situated at the diagonal (off diagonal) sites within the trivial region of the CSC and vice versa for the non-trivial part. Figure 3b shows the eigenfrequencies of the CSC containing diagonal non-hermiticity with Γ = 0.03. Similar to the Hermitian case, the spectral locations of the edge, bulk and corner states obtained by means of the aforementioned numerical experiments agree very well with the eigenfrequencies calculated by FEM. The positions of the sources and measurement points considered for the MST predictions are the same as those considered in the Hermitian CSC as illustrated in Fig. 3a. By comparing Figs. 1b and 3b, we can conclude that the real eigenfrequencies are almost identical in both cases. However, the pressure distribution of the corner states shown in Fig. 3c differs from the one in Fig. 1c, since the pressure is not the same in all of the corners. The states at the off-diagonal corners confine more energy (amplifying) than those at the diagonal sites (attenuating). In particular, the pressure of the attenuating corner states is 60% in comparison to the one confined via the amplifying corner state as Fig. 3d shows. Concerning the study of the robustness in this non-Hermitian configuration, we conduct the same defect study as the one described for the Hermitian case. But now, due to the fact that the pressure of the corner states, as discussed earlier, is asymmetric, the defect analysis has to be undertaken for both the amplifying and the attenuating corner states. Hence, the topological robustness study can be done by focusing only on the upper two corners as depicted in the insets of Fig. 4a, b. Again, the analysis comprises in steady increase of the number of interstitial defect cylinders ndefects. Both corner states display a remarkable durability against defects, though unlike the Hermitian case, the present non-Hermitian corner states display slightly differing shifts and pressure drops when ndefects increases as seen in the two spectra in Fig. 4. The pressure field maps in Fig. 4c also confirms the pertinacity against defects where opposite corners appear with slightly different pressure values.

Fig. 3
figure 3

Concentric square crystal (CSC) with diagonal non-Hermiticity. a Magnified view of the CSC displaying the arrangement of the trivial and nontrivial regions. b Comparison between the real eigenfrequencies obtained by finite element method simulations (dots) with the normalized pressure of the excited states from the multiple scattering theory computation. c Pressure maps of all corner states. Inserts illustrate the excitation point. d Pressure difference between the amplifying and attenuating upper corner states

Full size image
Fig. 4
figure 4

Pressure spectra with growing number of interstitial defects computed at (a) the upper-right and (b) the upper-left corner of the concentric square crystal with diagonal non-Hermiticity. c Mapping the pressure distribution of the corner states in the presence of 0, 3 and 7 defects

Full size image

The second non-Hermitian arrangement of the CSC is composed of parallel symmetrically gain and loss units as the magnification of Fig. 5a shows. In this scenario we obtain a mirror-like symmetry between gain and loss components that are embedded into the trivial and the non-trivial regions. The pressure spectra shown in Fig. 5b again display a perfect agreement between MST and FEM. Unlike the previously discussed corner states sustained in the diagonal non-hermiticity CSC, whose attenuating (amplifying) states confine at the diagonal (off-diagonal) corners, the parallel CSC on the contrary, binds the amplifying modes at the two right corners whereas the attenuating ones at the two left corners, as Fig. 5c shows. Furthermore, the pressure difference between the attenuating and the amplifying corner states as predicted in Fig. 5d appears to be remarkably pronounced.

Fig. 5
figure 5

Concentric square crystal (CSC) with parallel non-Hermiticity. a Magnified view of the CSC displaying the arrangement of the trivial and nontrivial regions. b Comparison between the real eigenfrequencies obtained by finite element method simulations (dots) with the normalized pressure of the excited states from the multiple scattering theory computation. c Pressure maps of all corner states. Inserts illustrate the excitation point. d Pressure difference between the amplifying and attenuating upper corner states

Full size image

From this follows that we should expect a highly selective response upon introducing defects when comparing the corner states situated at the left and right.

Indeed, the pressure spectra at the upper-right corner in Fig. 6a display a remarkable peak when compared to the opposing corner state computed in Fig. 6b. Curiously, the left corner state displays a pressure enhancement with growing interstitial defect cylinders up to ndefects = 3 above which, we witness the aforementioned shift and pressure drop, however, as seen in Fig. 6c all corner states appear clearly excited in real space of the CSC.

Fig. 6
figure 6

Pressure spectra with growing number of interstitial defects computed at (a) the upper-left and (b) the upper-right corner of the concentric square crystal with parallel non-Hermiticity. c Mapping the pressure distribution of the corner states in the presence of 0, 3 and 7 defects

Full size image

Discussion

In conclusion, we have conducted numerical studies based on a multiple scattering theory of non-Hermitian sonic second-order topological insulators. The predicted spectral locations of the corner states in response to an acoustic point source agree very well with eigenfrequencies computations by a FEM. Our simulations further confirm the topological robustness of these corner states in the presence of interstitial cylindrical crystal defects. Interestingly, depending on the geometrical arrangement of the non-Hermitian components we are able to excite corner states of variable acoustic energy. Hence, our powerful semi-numerical predictions can reliably pave the way for the design of new acoustic devices where topological protection in combination with lossy and amplifying ingredients are in great demand.

Methods

Multiple scattering theory

The semi-analytical method used to compute the pressure distribution and the resonance spectra was developed by means of the multiple scattering technique, whose formal definition has been derived in the previous result section.

Finite element method

Numerically computed eigenmodes of the SOTI were implemented using COMSOL Multiphysics, a finite-element analysis and solver software. The simulations were performed in the pressure acoustic module. The largest mesh element size was less than a tenth of the shortest incident wavelength. Plane wave radiation conditions were imposed on the exterior of the air domain to eliminate interference from the reflected waves.