Broadband single-phase hyperbolic elastic metamaterials for super-resolution imaging

Hyperbolic metamaterials, the highly anisotropic subwavelength media, immensely widen the engineering feasibilities for wave manipulation. However, limited by the empirical structural topologies, the reported hyperbolic elastic metamaterials (HEMMs) suffer from the limitations of the relatively narrow frequency width, inflexible adjustable operating subwavelength scale and difficulty to further improve the imaging resolution. Here, we show an inverse-design strategy for HEMMs by topology optimization. We design broadband single-phase HEMMs supporting multipolar resonances at different prescribed deep-subwavelength scales, and demonstrate the super-resolution imaging for longitudinal waves. Benefiting from the extreme enhancement of the evanescent waves, an optimized HEMM at an ultra-low frequency can yield an imaging resolution of ~λ/64, representing the record in the field of elastic metamaterials. The present research provides a novel and general design methodology for exploring the HEMMs with unrevealed mechanisms and guides the ultrasonography and general biomedical applications.

where 0 x F and 0 y F represent the total induced forces on the boundaries along the x-and y-directions, respectively; ω is the angular frequency; V denotes the volume of the effective medium; ρ xx and ρ yy are the effective mass densities in the x-and y-directions, respectively; and the applied displacement field is represented by 0 x U and 0 y U . The off-diagonal elements of the mass density tensor are zero because of the orthogonal symmetry of the unit-cell's microstructure. By applying the corresponding unit eigenstates, the effective stiffness E xx can be determined from the energy equivalence between the unit-cell and the effective medium by * * ( ) is noted here that we focus on the non-propagating waves along one principal direction to obtain the hyperbolic dispersion. Therefore, the effective elastic modulus E xx or E yy should be completely generated by the induced normal motion along the x-or y-direction. In this case, the applied eigenstates are different from the work by Liu et al. [1]. However, we can also utilize the longitudinal wave modulus P=K+μ (where K is the effective bulk modulus and μ is the effective shear modulus) to characterize the whole effective behaviors concerning the longitudinal wave motion. The detailed process for the numerical determination of K and μ can be found in the work by Liu et al. [1]. Figure S1 shows the wave transmission computation of an incident elastic wave in an optimized metamaterial.

Wave transmission along the two principal directions
Let us consider a plane time-harmonic longitudinal wave along the y-direction from the left into the elastic metamaterial, where u y is the amplitude of the incident wave field, k 0 represents the wave vector in the background material (stainless steel), and  is the operating angular frequency [11].
For the sake of simplicity, the time-harmonic factor where R and T are the reflection and transmission coefficients, respectively [11]. The similar computational procedure can be used for the case of an incident wave along the x-direction.

Transmission of the propagating and evanescent elastic waves
All our optimized metamaterials show an outstanding deep-subwavelength imaging benefited from the transmission enhancement of the evanescent waves, which can convert the subwavelength information into the propagating components. In order to characterize the wave transmission [3,4], we perform the numerical analysis of the propagating and evanescent waves transmitting through a layer of the anisotropic effective medium between two semi-infinite homogeneous stainless spaces (free-space). In this paper, we only focus our attention on the hyperbolic dispersion of the longitudinal waves by neglecting the transverse waves. Due to the nearly constant value, the longitudinal wave modulus P can be utilized to describe the whole longitudinal wave motion. In this case, the anisotropic mass density dominates the formation of the hyperbolic dispersion of the elastic waves. Referring to the study reported by Zhu et al. [5], the dispersion relation for the longitudinal wave propagating in the EMMs with an anisotropic mass density can be expressed by where k x and k y are the components of the wave vector in the anisotropic effective elastic medium along the xand y-directions, respectively. The transmission coefficient T of the propagating and evanescent waves is determined as [3,4] where Z y =ρ yy /k y and Z 0y =ρ 0 /k 0y are the wave impedances; ρ 0 is the mass density of the background medium; L is the thickness of the layer in the optimized metamaterial; is the component of the wave vector in the free-space; and k 0 =2/ denotes the propagation constant of the fundamental waveguide mode [6]. If k x ≤ k 0 , the transmission coefficient T characterizes the transmission property of the propagating waves，while for k x >k 0 the corresponding waves represent the evanescent waves. For a microstructure generated from the GA, substituting the effective wave impedances into Eq. (SQ4) yields the transmission coefficient T at a fixed frequency.

Description of the objective function
Inspired by the mechanism that the dipolar resonances [1,2,5,7] can produce a negative mass density, we should search for the topology of the microstructure to induce a resonance at a certain frequency. The effective mass density ρ xx increases to the positive infinity when the operating frequency below the resonant frequency rises. However, a negative ρ xx occurs when the operating frequency above the resonant frequency increases. As a result, a negative ρ xx in a certain frequency range is obtained by the topology optimization scheme [8,9,10]. To achieve a wider negative range in the lower-frequency region, it is necessary to push down the resonant frequency and make the overall positive ρ xx smaller [7]. If we compute the effective parameters at the sampling frequencies which are equably distributed in the target operating frequency range (f min , f max ), then the most important driving force to generate a negative ρ xx in a wide frequency range is to enlarge the gap between the maximal and the minimal positive values at these sampling frequencies. Therefore, min( will be maximized to achieve this purpose. Once the resonances take places within the target frequency range, the values of ρ xx at several sampling frequencies become negative. In this case, N is introduced in Eq. (1a) to ensure the individual with a negative ρ xx to be completely superior to others, thus guiding the evolution to find out more individuals with larger N. In this way, topology optimization can explore metamaterials with a broadband negative mass density.

Description of the constraints
In our optimization formulation for the negative ρ xx , we introduce seven constraints (1b)-(1h). The corresponding physical reasons for these constraints are given as follows.

 Constraints (1b)-(1d): Warranty for the longitudinal wave motion in the y-direction
To realize the single-negativity for the hyperbolic dispersion, we have to guarantee the eigenstate motions for the positive ρ yy , E xx and E yy simultaneously when hunting for the negative ρ xx . In addition, E yy should be larger than E xx to generate the longitudinal band along the y-direction over a large frequency range in which a bandgap in the x-direction is opened. These factors will be taken into account by the constraints (1b)-(1d), which can ensure the emergence of the longitudinal wave motion along the y-direction.

 Constraints (1e): Warranty for the exclusive longitudinal wave motion in the y-direction
When ρ xx turns negative from positive, both E 11 and E 22 usually decrease, owing to the increase of the void fraction in the metamaterial. However, the multipolar resonances revealed in this paper can effectively induce not only a negative ρ xx but also a negative E xx (although very small). If E xx is positive and close to zero, however, it means that the vibration in the x-direction can be easily coupled with the excited longitudinal wave motion along the y-direction. This results in the complex hybrid vibration modes containing the longitudinal and transverse wave motions. That is, the longitudinal band along the y-direction may show the quasi-longitudinal wave motion. Therefore, we introduce the constraint (1e) to control the coupling stiffness E 12 in order to obtain an exclusive longitudinal wave motion in the y-direction.
 Constraint (1f): Warranty for a purely translational motion in the y-direction Normally, we can find the strong "local rotation" which shows the off-diagonal feature of the effective mass density and stiffness tensors beyond the classical linear elasticity theory. In particular, the topology optimization usually involves many extremely complex microstructures which can bring more serious "local rotation" problems. Accordingly, we have to introduce the constraint (1f) to make sure that the induced behavior is completely the translation motion in the y-direction when retrieving ρ yy . Otherwise, the performance of ρ yy cannot be coincident with the dispersion relation.
 Constraint (1g): Warranty for a strong anisotropy A strong anisotropy is the most important feature of the hyperbolic metamaterials. To obtain the strongly anisotropic mass densities along the two principal directions, with ρ yy as a nearly constant value when ρ xx varies from the positive infinity to the negative infinity. Because a nearly constant ρ yy indicates a nearly non-dispersive performance along the y-direction, the extremely different features in the two directions naturally contribute to the hyperbolic dispersion with a maximal anisotropy. So, the constraint (1g) is introduced here to achieve this goal. The numerical tests show that   =1.37 can effectively balance the requirements of the strong anisotropy and the large feasible solution space.

 Constraint (1h): Warranty for a sufficient stiffness
The constraint (1h) is introduced to design only such a metamaterial with a sufficient stiffness for practical manufacturing. Naturally, this constraint can also suppress the mesh-dependency problem [8,10] encountered in topology optimization.

Evolutionary history and convergence of the topology optimization strategy
To check the convergence of the present topology optimization strategy, Figure S2 presents the evolutionary history of the maximal fitness as a function of the generation number for the optimized HEMM H3 in Fig. 1(d).
The metallic structure without any void is taken as the initial "seed" structure, see snapshot of the generation 0 (G=0). From the generation G=0 to G=61, the snapshots show that the GA can quickly capture the beneficial geometry at the early evolution, i.e., thin left and right boundaries, to increase the gap between the minimal and the maximal positive ρ xx . The change from the generation G=61 to G=145 implies that the solid blocks are useful to further enlarge the gap. A significant improvement of SN takes place from the generation G=145 (SN=−0.7732) to G=209 (SN=0.7663). This implies that the geometry of the narrow boundaries with four big solid blocks and two small solid lumps in the center is effective to generate a negative ρ xx at ultra-low frequencies. When the local connections between the blocks and the boundaries get narrower from the generation G=209 (SN=0.7663) to G=229 (SN=2.6583), the frequency range with a negative ρ xx becomes wider. Then, this increasing tendency can be enhanced by removing more solid parts of the connections from the generation G=229 (SN=2.6583) to G=305 (SN=3.631). The negative ranges will increase from the generation G=305 (SN=3.631) to G=553 (SN=4.636) because of the growth of the two central lumps. At this stage, the evolution in the coarse grid begins to converge. Mapping the topology from a coarse grid into a finer one will lower SN at some extent. In the fine grid, the GA can commendably make the structural edges clearer and more smooth. From the generation G=1004 (SN=4.612) to G=2000 (SN=4.7623), the structure turns to possess the extreme narrow connections and sufficiently large solid blocks. Overall, the evolutionary history in Fig. S2 evidently demonstrates the significance of the typical structural features on the generation of a broadband negative ρ xx . Fig. 1(d)

Hyperbolic properties of the optimized HEMM H2
For completeness to demonstrate the hyperbolic dispersion for all optimized HEMMs in Fig. 1(d), we illustrate in Fig. S3 the corresponding dispersion relations, transmission properties of the longitudinal waves and the effective material parameters for H2. It is clearly seen that the different physical quantities match each other very well. The negative ρ xx with the positive E xx can accurately capture the occurrence of the bandgap along the ΓX-direction. The positive ρ yy with the positive E xx also predicts the existence of the longitudinal wave mode in the ΓY-direction. The transmission properties along the two principal directions also show that the longitudinal waves cannot propagate within the metamaterial in the ΓX-direction but can propagate along the ΓY-direction. Unlike the results in Fig. 2, the band structure in Fig. S3(a) shows that only the longitudinal wave mode exists along the ΓY-direction in the ΓX-directional bandgap range. Fig. 1(d). In addition, we illustrate in Fig. S4(a) the EFCs of the third band for H2. We can clearly observe the broadband hyperbolic dispersions. Furthermore, the EFCs become very flat at low frequencies as well. Since we only focus on the longitudinal wave propagation, the optimized metamaterials can act as the HEMMs as long as the EFCs for the longitudinal waves have the hyperbolic shapes. So the evident hyperbolic dispersions shown in Figs. 3(a), 3(b) and S3(a) validate that our proposed optimization formulation described by Eqs. (1a)-(1h) is robust for the longitudinal waves, no matter whether the transverse wave propagation exists or not. To further confirm the hyperbolic properties of H2, Figs. S4(b) and S4(c) show its subwavelength imaging at 5 kHz and 5.5 kHz, respectively. In these two cases, the imaging resolutions are as high as 0.1λ and 0.098λ, respectively. Compared with the imaging resolutions of H1 and H3, we can reasonably conclude that the frequency dominants the ability of the imaging resolutions. Furthermore, the optimized multipolar resonant microstructures in Fig. 1(d) can maintain the excellent hyperlensing effect over the sufficiently wide wavelength range. Figure S4(d) presents the longitudinal wave propagation in H2 at 5 kHz to identify its anisotropic property. It is clear that the longitudinal wave indeed propagates only along the y-direction, while no energy can come out from the metamaterials in the x-direction. In fact, the hyperbolic curves over large wave vector ranges imply that the metamaterial H2 can realize the negative refraction within the wide angle ranges. Then almost all waves can be highly focused on the top and bottom boundaries. Therefore, like H1 and H3 in Fig. 1(d), the optimized HEMM H2 also shows the broadband hyperbolic dispersion, highly anisotropy and prominent super-resolution imaging. All these benefits show the effectiveness and universality of the topology optimization scheme presented in this paper.

Optimization with impedance match constraint
A high wave transmission is very important for the wave imaging. Generally, in terms of the microstructural design, a better impedance match can improve the imaging transmission regardless of the superlens thickness. When designing the HEMMs, one common approach is to introduce a very simple relative impedance  Figure S5(a) shows the optimized HEMM which has a negative ρ xx and a negative E xx within the frequency ranges of (10.859 kHz, 19.516 kHz) and (19.604 kHz, 28.654 kHz), respectively. Both ρ yy and E yy are always positive in the same frequency ranges. This proves the occurrence of the hyperbolic dispersion in the frequency range of (10.859 kHz, 19.516 kHz) in Fig. S5(b). The optimized HEMM in Fig. S5(a) is composed of six blocks and four thick connections, which has a larger mass density and a stronger stiffness than that of H1 in Fig. 1(d). The increased stiffness gives rise to a better impedance match. We plot in Fig. S5(c) the corresponding imaging simulation at 14 kHz. Here, an improved larger imaging transmission and a noticeable higher imaging resolution (FWHM=0.138λ) are obtained, which demonstrates the positive effect of the proposed impedance match constraint. However, a better impedance match often requires the simultaneously large mass density and stiffness, leading to an incapability within the ultra-low frequency region (>=90a).