## Abstract

Acoustic metamaterials have been proposed for numerous applications including subwavelength imaging, impedance matching, and lensing. Yet, their application in compressive sensing and imaging has not been fully investigated. When metamaterials are used as resonators at certain frequencies, they can generate random radiation patterns in the transmitted waves from the transducers and received waves from a target. Compressive sensing favors such randomness inasmuch as it can increase incoherence by decreasing the amount of mutual information between any two different measurements. This study aims at assessing whether the use of resonating metamaterial unit cells in a single-layered non-optimized array between a number of ultrasound transceivers and targets can improve the sensing capacity, point-spread function of the sensing array (their beam focusing ability), and imaging performance in point-like target detection. The theoretical results are promising and can open the way for more efficient metamaterial designs with the aim of enhancing ultrasound imaging with lower number of transceivers compared to the regular systems.

## 1 Introduction

Inspired by the theoretical formulation and experimental realization of left-handed electromagnetic metamaterials [1–3], the existence of double-negative acoustic metamaterials were shown in Ref. [4]. The principles of acoustic metamaterials were founded on similar concepts discussed in electromagnetics by Veselago [5]. The introduction of acoustic metamaterials has paved the way for diverse applications such as cloaking [6–8], focusing and lensing [9–14], elastic wave manipulation [15,16], subwavelength imaging [17,18], sound absorption [19–21], canceling out aberrating layers [22], and impedance matching [23,24]. More practically, real-life applications associated with these new opportunities vary from anti-earthquake infrastructures for buildings [25,26] to penetrating bones via ultrasound waves [22].

Particularly, to enhance imaging capability, Zhu et al. proposed a holey-structured metamaterial to exceed the limit that diffraction sets on the image resolution [18]. The structure included 1600 holes made into a rigid brass block and it caused Fabry–Pérot resonances to be formed and successfully transmitted to the other side. Through simulation and experimental results, it was shown that details down to 2% of the exciting wavelength could be restored in the image from an object. Previously, a metamaterial hyperlens [13] and a single-negative metamaterial [17] were utilized to achieve resolutions of 14.7% and 13% of the used wavelength, respectively.

In this work, aiming at improving the sensing capacity, the idea of acoustic metamaterials has been combined with the concept of compressive sensing by a randomization process that enables collecting more information from the imaging domain (ID). Compressive sensing theory states that using nonadaptive linear projections, one can reconstruct a signal by samples taken at a much lower rate than what Shannon–Nyquist theorem requires [27]. One way to take advantage of compressive sensing is to create random measurements, that is, to randomize the outgoing signal from the transmitter and the incoming signals to the receivers [28]. This randomization can effectively enhance sensing mechanisms [29], which may seem unexpected at first glance. When looked into the matter more closely, however, the counter-intuitiveness of the role of randomizing the measurements in enhancing the sensing capacity can be resolved: random signals minimize the amount of mutual information between two successive measurements and thus they maximize the information that can be transferred from the imaging domain to the receivers. Utilizing metamaterials that exhibit resonant behavior at certain frequencies is one way to randomize the measurements and has been adopted in this work.

To the knowledge of the authors, the use of ultrasound metamaterials has not yet been reported for megahertz frequency imaging. This could be due to the challenges in fabricating such fine structures, since obtaining resonant frequencies in the megahertz domain requires features at micrometer scale. However, the use of metamaterials to create an ultrasound matching layer [30] and to manipulate acoustic waves [31] has been experimentally and theoretically studied, respectively, at megahertz frequencies. In the first reference, features of sizes down to 122 *μ*m were fabricated, and in the second one, employing micro gas bubbles to alter ultrasound waves were introduced. Here, the use of metamaterial resonators in increasing the sensing capacity, enhancing beam-forming ability or point-spread function (PSF) of the sensing array, and improving the imaging performance is investigated, in a theoretical, proof-of-concept fashion.

## 2 Metamaterial Design and Analysis

This section first describes the design of a resonant metamaterial unit cell to merely have at least a resonance in the working frequency band of the sensing system, without any optimization or deterministic approach. Next, the design of a randomly distributed line-array of two different sizes of the designed unit cell will be explained, which will be used in an imaging application in the upcoming section. The exact size of the unit cell, scaling factor, number of the unit cells in the array, and the distribution method are arbitrary to make the process as random as possible.

### 2.1 Metamaterial Unit Cell.

*f*

_{new}and

*f*

_{old}denote the new and old resonant frequency, respectively. As shown in Fig. 1, the square unit cell is made up of a steel central part coated by layers of rubber and aluminum, in order, in a water background. The radii of the three layers are

*r*

_{al}= 0.25

*a*,

*r*

_{rub}= 0.2

*a*, and

*r*

_{st}= 0.15

*a*, from outside to the center, with

*a*being the square unit cell side. The numerical simulations were carried out in comsol multiphysics, and the reflection (

*R*) and transmission (

*T*) coefficients were calculated on the left and right boundaries of the simulation domain according to the following formulae, respectively:

*P*

_{t}and

*P*

_{i}denote the total and incidence pressure fields, respectively;

*h*is the height of the simulation domain;

*k*

_{b}is the background wavenumber;

*f*is the guided frequency; and

*d*

_{1},

*d*

_{2}, and

*L*are geometric dimensions and are defined in Fig. 1. In the configuration shown in the figure,

*h*=

*a*= 20

*μ*m,

*d*

_{1}=

*d*

_{2}= 60

*μ*m, and

*L*= 70

*μ*m. The properties of each material, as used in the simulations, are given in Table 1. It is clear from the transmission and reflection curves that the metamaterial unit cell is working in a band-pass mode, reflecting most of the incoming waves and only passing the wave through almost completely at the resonant frequency, that is 4.5 MHz.

*homogeneity condition*[35]. Therefore, one can set the conditions on the size of a unit cell and its contents in this study as follows: $a<\lambda min/4(=53\mu m)$, 0 <

*r*

_{al}<

*a*/2, 0 <

*r*

_{rub}<

*r*

_{al}, 0 <

*r*

_{st}<

*r*

_{rub}, with $\lambda min$ being the minimum wavelength used in imaging. Now, the approach that was introduced by Fokin et al. for locally resonant acoustic metamaterials [36] is used to obtain the effective parameters. When

*R*(

*f*) and

*T*(

*f*) are available, the effective refractive index

*n*and impedance $\xi $ in the case of normal incidence can be found as [36]

*q*= 0, ±1, ±2, ±3, … is the branch number of the inverse cosine function and

*d*is the unit cell thickness in the direction of the waves. To remove the ambiguities that arise in the signs and the branch number in the above formulations, Fokin et al. set two conditions: (i) physical realization of the solution, i.e., $Re(\xi )>0$, with Re(·) showing the real part operator on a complex number, and (ii) the use of minimum metamaterial thickness for which

*q*= 0. Their approach is summarized schematically in Fig. 2. Another method to resolve the mentioned ambiguities is introduced by Szabó et al. [37], using Kramers–Kronig relationship between the real and imaginary parts of an analytical complex function such as the refractive index.

The retrieval results are shown in Fig. 3 for the designed unit cell, using both Eqs. (3) and (4) and the algorithm outlined in Fig. 2. In adopting the mentioned equations, the positive signs have been used. As can be seen, the algorithm has modified the signs of the parameters at certain frequency intervals. Furthermore, the resonant behavior observed in transmission and reflection curves also appears in the effective parameters plots, exhibiting an extremum at 4.5 MHz. This property will be employed to create the desired randomness in the compressive sensing approach, as described in the next section. It should be noted that inasmuch as the density and bulk modulus take negative values, the unit cell is of double-negative type.

### 2.2 Metamaterial Line Array.

Scaling the designed metamaterial will create unit cells of different resonant frequencies. If the original metamaterial unit has a dimension of *r*_{0} and a resonant frequency of *f*_{0}, scaling its dimensions by a factor $\sigma $ creates a metamaterial unit of dimension $\sigma r0$ and resonant frequency $f0/\sigma $. As an example, upscaling the multi-layered circle inside the designed cell by a factor of 1.2 yields a unit cell of resonant frequency 3.75 MHz. A combination of such unit cells in an array can create multiple resonant frequencies. Placing the resultant array in the way of the waves coming from the acoustic sources will randomize the wave fields that reach a target by altering their patterns at different frequencies.

In this study, two types of unit cells, one having a resonant frequency of 4.5 MHz and the other having a resonance at 3.75 MHz, are arranged randomly in a line array and placed between the sources and the targets. As mentioned, the size of the unit cells and their contents are selected so that they give one or more resonances in the operating frequency range. Moreover, the unit cells are placed side by side, without any space between them, so that they cover the entire width of the domain. It is important to note that in this work, to take advantage of the favored randomness in compressive sensing, no optimization will be done in the design of the individual unit cells or the array. Rather, unit cells with arbitrary resonances in the frequency band and their random distribution will be relied on to make the sensing matrix less coherent. Regarding the behavior of the array, the following remarks are noteworthy:

Using these cells together in a two-element array will not necessarily create isolated resonant frequencies precisely at the predicted values as shown in Fig. 4.

When the waves are not perpendicular to the metamaterial unit cell—as in the case of waves coming off a point source and reaching at the boundary of unit cells that are not directly situated below the source—the behavior of the metamaterial might change significantly. When the incident fields are not normal to the unit cell boundary, another set of formulations are required to restore the effective properties of the metamaterial cell [38]. Although these effects cannot be easily controlled, in this case, they are advantageous in increasing the randomness of the measurements.

As shown in Fig. 4, the upscaling of the base unit cell has created a transmissive response at 6.05 MHz. This can be viewed as if the *T*–*R* plot has been contracted along the frequency axis and has brought in a new resonant frequency that was outside the frequency band for the base metamaterial. The interaction between the unit cells has also generated other extermums on the plot at 2.755 MHz and 4.35 MHz that are not of reflective or transmissive nature.

## 3 Ultrasound Imaging Using Metamaterials

To assess the performance of a sensing system in the presence of resonant metamaterials as measurement randomizers, an imaging simulation is carried out with an array of ultrasound transceivers and a number of targets. In this section, the simulation domain and its parameters are introduced. Furthermore, the imaging algorithm that is utilized to retrieve the image of the targets is established.

### 3.1 Forward Modeling and Simulation Setup.

*f*denoting the frequency, $\rho (r)$ is the density of the medium, $k(r,\omega )=\omega /c$ is the wave-number with

*c*being the speed of sound, $P(r,\omega )$ is the pressure field, and $F(r,\omega )$ is the acoustic source term. The source term can be defined in various ways. If it is a monopole, and the root-mean-square (RMS) of power per unit length, that is

*P*

_{rms}, is used as excitation, then [40]

**r**

_{s}are the phase and the location of the source, in order, and $\delta (r\u2212rs)$ is the shifted Dirac delta function. The dependency of $F(r,\omega )$ on frequency is due to the two-dimensional (2D) modeling; however, in 3D modeling, this dependency does not exist [40]. Thus, in the 2D case, for the exciting source term to have a magnitude of one and no phase terms, the parameters must be chosen as $Prms(\omega )=1/(8\rho (r)\omega )$ and $\varphi s=0$.

The definition of the geometric parameters of the domain is given in Table 2, alongside with the numerical values used in this case study. Since this study is closely related with the idea of multimodal early detection of breast cancer using a mechatronic system [41,42], the domain was set up accordingly, including an acrylic interface medium between the coupling liquid and the imaging domain.

Parameter | Definition | Numerical value (mm) |
---|---|---|

w_{D} | Domain width | 6 |

h_{D} | Domain height | 6 |

t_{c} | Coupling medium height | 2 |

t_{i} | Intermediate medium height | 2 |

t_{b} | Background medium height | 2 |

d_{i} | Circular target diameter | 0.2 |

(x_{s,1}, y_{s,1}) | Source 1 location | (−0.25, 2.9) |

(x_{s,2}, y_{s,2}) | Source 2 location | (0.25, 2.9) |

(x_{t,1}, y_{t,1}) | Target 1 location | (−0.6, 2) |

(x_{t,2}, y_{t,2}) | Target 2 location | (0.8, 2) |

Parameter | Definition | Numerical value (mm) |
---|---|---|

w_{D} | Domain width | 6 |

h_{D} | Domain height | 6 |

t_{c} | Coupling medium height | 2 |

t_{i} | Intermediate medium height | 2 |

t_{b} | Background medium height | 2 |

d_{i} | Circular target diameter | 0.2 |

(x_{s,1}, y_{s,1}) | Source 1 location | (−0.25, 2.9) |

(x_{s,2}, y_{s,2}) | Source 2 location | (0.25, 2.9) |

(x_{t,1}, y_{t,1}) | Target 1 location | (−0.6, 2) |

(x_{t,2}, y_{t,2}) | Target 2 location | (0.8, 2) |

From the top of the figure to the bottom, respectively, three layers of different materials are used: (i) water in which the acoustic transceivers are placed, (ii) acrylic sheet which is assumed to be the bottom of the water container or the breast compression paddle, (iii) background medium (fibroglandular breast tissue) that encompasses a number of targets (invasive ductile carcinoma (IDC) tumors). The properties—density and speed of sound—of each of these materials, except water whose properties were given in Table 1, are listed in Table 3. The longitudinal speed of sound is given for an IDC mass in Ref. [43], but to compute the tumor’s density, the data given in Ref. [44] on the mean values of elasticity (*E*_{idc}), sound shear velocity (*c*_{sw,idc}), and the relationship between the two can be used. Having *E*_{idc} = 140.7 kPa and *c*_{sw,idc} = 6.7 m/s, the density can be calculated as $\rho idc=Eidc/3csw,idc2=1077kg/m3$ [44]. The geometry as well as the simulation were set up using the COMSOL Multiphysics with matlab module. The simulation was set up to scan the medium by a total number of 61 frequencies in the band 1–7 MHz.

### 3.2 Inverse Modeling and Imaging Algorithm.

**r**′. If the scattered field is defined as $Ps(r,\omega )=Pt(r,\omega )\u2212Pb(r,\omega )$, then Eq. (5) can be rewritten as

*s*that stands for the

*scattered field*. Thus, the solution to (12) can be found by (9), after making the relevant changes in the subscript and the source term as follows:

*X*(

**r**) takes a relatively small value of 7.47% inside the tumor and it is equal to zero everywhere else in the imaging domain. To determine

*X*(

**r**) from (13), the value of the total pressure $Pt(r,\omega )$ and scattered pressure field $Ps(r,\omega )$ are required at each point of the imaging domain. This is not feasible in most cases as the number of transceivers that can be placed in the space is limited. Yet, since the contrast value is small in the imaging domain, the total field can be estimated with the background field using the first-order Born approximation [52,53], which yields the following equation:

*n*in the imaging domain and a particular number of measurements

*m*, assuming noiseless simulations:

**is the sensing matrix, $x^$ is the unknown vector—the discretized and vectorized version of**

*A**X*(

**r**)—and

**b**is the measurement vector. In other words, $x^$ is the column-wise stacking of the values of the contrast parameter, in 2D space, at each pixel in the imaging domain. As mentioned, due to the limitation on the number of measurements, the number of unknowns is much larger than that of the known variables and thereby Eq. (15) is underdetermined, having infinitely many solutions. Consequently, a regularization method is needed to find an optimized solution for the linear system.

*i*th element of $x^$; and $\delta $ is an upper limit for the residual error $\Vert Ax^\u2212b\Vert \u21132$, with ℓ

_{2}denoting the norm-2. NESTA—short for Nesterov’s method—can be used to solve the above minimization problem. The code that implements this algorithm is available online

^{2}and it requires the input of the norm-1 regularization parameter

*μ*to optimize the solution.

### 3.3 Sensing Capacity and Point-Spread Function.

*n*

_{min}= min (

*n*,

*m*) orthogonal parallel channels as [28,55]

*P*

_{i}/

*n*

_{0}is the signal-to-noise ratio (SNR) of the

*i*th channel, with

*P*

_{i}and

*n*

_{0}being the power of the signal and the noise, respectively, and $\lambda i$ is the

*i*th nonzero singular value of the sensing matrix that can be obtained from its singular value decomposition [28]

**= (**

*U**u*

_{1}, …,

*u*

_{m}),

**= (**

*V**v*

_{1}, …,

*v*

_{n}), and $\Xi =diag$$(\lambda 1,\u2026,\lambda nmin)$. The other criterion is the beam focusing ability of the system or its PSF that can be established using the phase compensation focusing method as follows [56,57]:

*N*

_{t}is the number of transmitters,

*N*

_{f}is the number of frequencies,

*P*

_{b,ik}(

**r**,

*f*

_{k}) is the background fields due to transmitter

*i*at frequency

*f*

_{k}at location

**r**, and $\varphi p,ik$ is the phase of the background field at the location of the focus point

**r**

_{p}= (

*x*

_{p},

*y*

_{p}), that is $\u2222Pb,ik(rp,fk)$. Using this method, the acoustic beams can be focused at

**r**

_{p}, giving the highest intensity at the location of the point and fading intensity as one moves from the point. The quality of focusing is contingent upon various factors such as the number of frequencies, number of transmitters, and, as will be seen, the use of metamaterials.

## 4 Results

The results of the simulations are presented in this section and the impact of adding metamaterials between the transceivers and the targets is assessed. The sensing matrix was computed from the background pressure for each pixel of the imaging domain, at each frequency and for each transceiver. Having the sensing matrix and background pressure, sensing capacity *C* and beam focusing *BF*(**r**|**r**_{p}) are obtained immediately using Eqs. (17) and (19), for both cases where the metamaterials were present or absent. The beams were focused at the location of the tumors, as given in Fig. 5. The results are demonstrated in Figs. 6 and 7. It is clear that the addition of metamaterials has enhanced the phase-compensation focusing by narrowing down the region of the highest intensity around the focus points. Furthermore, it is evident that the metamaterials have increased the sensing capacity of the system, particularly at SNRs more than 20 dB.

The imaging results, obtained by NESTA, are shown in Fig. 8. These images are optimized by varying the NESTA algorithm parameters aimed at fewer artifacts and better localization. Also, they have been both saturated to only include the solution that is within the range 0 to −7 dB, for better demonstration of the contrast between the tumor and the background. It is clear that in the absence of the metamaterial line, the targets are not detected; rather, only their whereabouts in the *y*-direction has been revealed. On the other hand, the addition of the metamaterial line has enabled the system to detect the location of the targets—(−0.6, −2) mm and (0.8, −2) mm—both in *x*- and *y*-direction, with acceptable accuracy.

## 5 Conclusions

In this study, the use of metamaterials for increasing the sensing capacity, focusing ability, and imaging performance of ultrasound waves was discussed. Since compressive sensing technique in imaging, which enables signal retrieval at much lower rates than the Nyquist frequency, favors randomness, a metamaterial unit cell was designed at two different sizes to create random patterns in the waves at different frequencies. It was shown by three criteria—sensing capacity, beam focusing, and imaging using NESTA—that the addition of metamaterials between the transceivers and the targets can meaningfully enhance the performance of the sensing system. The design of the mostly reflective metamaterial in this work was rather preliminary and theoretical, to only exhibit transmissive resonant behavior at certain frequencies. What is more, the unit cell and array design were intentionally non-optimized to show the power of random processes in compressive sensing, which could lead to desirable results without any optimization. In future studies, it should be seen whether different metamaterial designs with different configurations will be able to improve the imaging even more (for a study of different metamaterial types and the impact of various random distributions of the unit cells, refer to [58]). Moreover, the fabrication aspects of such small, multi-layer metamaterials should be addressed.

## Footnote

## Acknowledgment

This work has been partially funded by the Department of Energy (Award DE-SC0017614) and the NSF CAREER Program (Award No. 1653671).

## Nomenclature

**x**=discretized contrast variable or the unknown vector

*C*=sensing capacity

=*A*sensing matrix

*k*(**r**, ω) =wavenumber as a function of location vector

**r**and angular frequency $\omega $- $n(f)$ =
refractive index as a function of frequency

*BF*=beam focusing

- $F(r,\omega )$ =
ultrasound source term as a function of location vector

**r**and angular frequency $\omega $ - $Gb(r,r\u2032,\omega )$ =
Green’s function as a function of transmitter location vector

**r**′, location vector at each point in the domain**r**, and angular frequency $\omega $ - $Pb(r,\omega )$ =
background pressure field as a function of location vector

**r**and angular frequency $\omega $ - $Ps(r,\omega )$ =
scattered pressure field as a function of location vector

**r**and angular frequency $\omega $ - $Pt(r,\omega )$ =
total pressure field as a function of location vector

**r**and angular frequency $\omega $ *R*(*f*) =reflection coefficient as a function of frequency

*T*(*f*) =transmission coefficient as a function of frequency

- $X(r,\omega )$ =
contrast variable a function of location vector

**r**and angular frequency $\omega $ - $\xi (f)$ =
impedance as a function of frequency

## References

*ε*and

*μ*