Over the past few years, acoustic gradient index metasurfaces (GIMs) have been actively studied for the numerous wave control capabilities that they facilitate. Previous research, however, has primarily focused on GIMs that operate in the audible frequency range, due to the difficulties in fabricating such intricate structures at the millimeter and submillimeter scales, for ultrasonic applications. In this work, we design, fabricate, and experimentally demonstrate the working of a hybrid resonant acoustic gradient index metasurface for airborne ultrasound at 40 kHz. The fabrication of such a GIM is made possible by projection microstereolithography, an emerging additive manufacturing technique. Numerical simulations were conducted to verify the metasurface design, and experiments were performed to corroborate these simulations. The stronger dissipation associated with airborne ultrasound is highlighted in this paper. The experimental demonstration of such a metasurface for airborne ultrasound could further its prospects as a candidate for miniaturized acoustic devices.

The emergence of metasurfaces1–11 has bolstered the interest of various engineering communities in compact wave-based devices. Several studies over the past decade have thus proposed subwavelength structures that facilitate passive wavefront manipulation.12–15 In acoustics, this not only has unfolded the realization of unconventional designs16–18 but also has revisited classical structures which could now be redesigned to provide unusual functionalities.4,14,15 In this context, gradient index metasurfaces (GIMs) have been ardently examined to propose several features such as focusing,17,19 bending,14 retro-reflection,20,21 and holographic rendering.22 Additionally, some recent works have examined and embraced the inherent dissipation23–25 in these GIMs, as not just a loss but an avenue to more unconventional applications such as asymmetric transmission7,26 and wide-angle absorption.27 A hybrid design consisting of shunted Helmholtz resonators (HRs) and a straight channel is a frequent candidate for the building blocks of these structures, due to their relative simplicity and high sound transmission. The fabrication of such structures has been enabled by conventional 3D printing techniques such as fused deposition modeling16 or stereo-lithography.28,29 The precision of these manufacturing methods, however, has rendered their implementation in the ultrasound range rather scarce.6,30 Metasurface-based designs operating at ultrasound frequencies would be advantageous for a variety of applications such as levitation,6,31 sensing, and imaging,32 owing to the better precision brought about by the smaller wavelengths.

On the other hand, ongoing developments in the field of mechanical metamaterials are accompanied by the rapid advances of manufacturing approaches that help put forward artificial materials with exceptional mechanical properties. Such materials possess complex three dimensional micro- and nanoarchitectures and hence require sophisticated manufacturing capabilities. Large area projection microstereolithography (PμSL),33,34 in this regard, is an emerging additive manufacturing technique that is capable of fabricating samples with high structural complexity and feature sizes ranging from a few micrometers to tens of centimeters—a characteristic that could benefit the fabrication of acoustic metasurfaces/metamaterials at ultrasound frequencies. In this work, we design, fabricate, and experimentally demonstrate the performance of a hybrid resonant acoustic gradient index metasurface (HRAGIM) that operates at 40 kHz. The fabrication of the designed prototype is made possible by the aforementioned additive manufacturing technique. The role of thermoviscous dissipation in its elementary units is discussed, and full wave simulations that incorporate these losses are put forward. The results from these simulations are then validated via experimental measurements obtained from a two-dimensional experimental platform. The challenges at ultrasound frequencies such as stronger dissipation are highlighted.

The HRAGIM has elementary units which consist of a series of four subwavelength Helmholtz resonators (HRs) and an open channel. While the shunted resonators provide the reactance to shift the phase of the incident wave, the open channels serve as subwavelength slits that enhance the rate of transmission due to Fabry–Pérot resonance. By varying the parameter h1, shown in Fig. 1(a), a full range of phase shifts from 0 to 2π can be obtained. Numerical simulations were first performed using the pressure acoustics module on COMSOL Multiphysics 5.3a, a commercial finite element package. The unit cell has a fixed height of h = 0.1λ and a width of w = 0.4λ, where λ is the wavelength of the incident wave. This was done in order to obtain values of h1, whose phase shifts (ϕ) are equally spaced. As will be discussed later, the material used for fabrication was a high stiffness UV curable acrylate polymer, and it is thus safe to assume that the walls of the resonators are acoustically hard since their impedance is much greater than that of air. The phase gradient, ξ, is then engineered,35,36 by arranging the appropriate unit cells in the required pattern where

ξ=dϕsdx=(sinθtsinθi)k0nG.
(1)

Here, θt and θi are the angles of refraction and incidence, respectively, k0 is the wavenumber at the operating frequency, n is the order of diffraction that is associated with the grating, and G is the reciprocal lattice vector. In this case, ξ = 2π/γ = 97.2 (2π rad m−1), where the array period, γ, of the metasurface is 10.28 mm. The critical angle of incidence for this metasurface can hence be calculated to be θc=sin1(1ξ/ko)=9.6°. As explained in previous works,12,13n vanishes for angles of incidence smaller than the critical angle (θi < θc). In this letter, we consider only the case of normal incidence, and hence, the second term in Eq. (1) can be neglected (n is the grating induced diffraction order and must not be confused with p, the overall diffraction that results from the interplay of the phase gradient and the period grating7). It should be noted that such a GIM for airborne ultrasound can also be formed utilizing other types of unit cells.16–18 In this paper, however, we employ the above-mentioned design due to the versatility of the associated hybrid resonances and the firm theoretical framework that previous research in such structures has offered. Here, it is worthwhile mentioning that the metasurface in this paper is an adaptation from the works of Li et al.,7,14,15 where a HRAGIM was designed to operate at 3.43 kHz. The solid black lines and the dotted red lines in Figs. 1(b) and 1(c) are from pressure acoustics simulations for unit cell designs that operate at 40 kHz and 3.43 kHz, respectively. This illustrates that the dimensions of the structure that was previously designed to operate at a lower frequency can be scaled down to work for airborne ultrasound. However, it is pivotal to examine the prospective challenges that such a metasurface could face, while operating at higher frequencies.

FIG. 1.

(a) The elementary unit of the HRAGIM, where the incident sound propagates in the positive x direction. To operate at 40 kHz, the dimensions of the unit are w = 3.43 mm, h = 0.86 mm, w2 = 0.75 mm, and h2 = 0.086 mm. The values of h1 are chosen to be 0.13 mm, 0.15 mm, 0.19 mm, 0.23 mm, 0.32 mm, and 0.48 mm to achieve a discretized phase distribution. (b) and (c) show the numerical prediction of the transmission coefficient and phase shift, respectively, through the unit cell as a function of h1/h. The solid black lines and dotted red lines indicate the results at 40 kHz and 3.4 kHz, respectively, for the case without dissipation. Similarly, the black and red dashed lines are the results shown when thermoviscous dissipation is included. (d) shows the diffractive efficiencies for the metasurfaces (with loss) operating at 3.4 kHz and 40 kHz, in comparison to that of the lossless GIM.

FIG. 1.

(a) The elementary unit of the HRAGIM, where the incident sound propagates in the positive x direction. To operate at 40 kHz, the dimensions of the unit are w = 3.43 mm, h = 0.86 mm, w2 = 0.75 mm, and h2 = 0.086 mm. The values of h1 are chosen to be 0.13 mm, 0.15 mm, 0.19 mm, 0.23 mm, 0.32 mm, and 0.48 mm to achieve a discretized phase distribution. (b) and (c) show the numerical prediction of the transmission coefficient and phase shift, respectively, through the unit cell as a function of h1/h. The solid black lines and dotted red lines indicate the results at 40 kHz and 3.4 kHz, respectively, for the case without dissipation. Similarly, the black and red dashed lines are the results shown when thermoviscous dissipation is included. (d) shows the diffractive efficiencies for the metasurfaces (with loss) operating at 3.4 kHz and 40 kHz, in comparison to that of the lossless GIM.

Close modal

Foremost, as discussed in prior literature,24,25 the energy dissipation in such structures is due to thermoviscosity. This is characterized by Δ=δ/h1, where h1 is the slit's transverse dimension [as seen in Fig. 1(a)] and δ=2ν/ω is the viscous boundary layer thickness. Here, ν = 1.45 × 10−5 m2/s is the kinematic viscosity of air, ω is the angular frequency, and Δω. Hence, the largest value of Δ at 3.43 kHz is only 3.67%, while it is 12.53% for the same value of h1/h (scaled) at 40 kHz. This was further examined by performing numerical simulations that take thermoviscosity into account. No-slip and isothermal conditions were imposed on the solid boundaries. The dashed black and red lines in Fig. 1(b) compare the results from these calculations at 40 kHz with those carried out on the hybrid resonant structure designed to operate at 3.43 kHz. It can be seen that with the incorporation of themoviscous effects, the decrease in the transmission coefficient is as expected, much larger in the case of the 40 kHz design than that of 3.43 kHz, as seen in Fig. 1(b). The low rate of transmission for smaller values of h1/h is also as anticipated. In addition to the boundary layer effects, the attenuation of sound, αs, in free space due to the thermal and viscous losses in air is directly proportional to the square of the frequency, f2. In the audible range, this attenuation is rather negligible and can be ignored. At higher frequencies, however, it becomes more important to take this attenuation into account. In our case, at 40 kHz, αs ≈ 0.2 Np/m (Ref. 37), and its inclusion in our simulations (both the unit cell and the upcoming full wave) was still found to have a negligible effect. This is because of the much larger loss from the boundary layer effects. It is hence important to note that although HRAGIMs can be scaled as a function of frequency, the effects of losses (Δ and αs) are more profound at higher frequencies and are not readily scalable. Furthermore, Fig. 1(c) shows that thermoviscous dissipation has a relatively weak influence on the phase shift at both 3.43 kHz and 40 kHz. This is true as dissipation dominates the resistance part of impedance more than the reactance.24 This was further verified by means of full wave thermoviscous acoustics simulations (shown in the supplementary material), through which the diffractive efficiency of the HRAGIM was calculated and compared with the lossless and 3.43 kHz cases. This is shown in Fig. 1(d), where it can be seen that p = +1 is the dominant diffractive mode in all cases. However, with the inclusion of loss, the p = 0 mode is enhanced, and the amplitude of the +1 mode reduces, as indicated by the dashed red lines. At 40 kHz, the same trend is seen, indicated by the dashed black line, with a further reduction in the amplitude. The HRAGIM can thus provide the desired wave bending effect for airborne ultrasound, however with a considerable loss in transmission.

The fabrication of such a minute design is made possible by a high resolution, large area projection microstereolithography system capable of fabricating microscale, high-aspect ratio features over a wide area. In contrast to other methods such as fused deposition modeling and UV projection waveguide systems, this approach is ideal for samples with high structural complexity and with a feature size ranging from microns to centimeters. A three-dimensional computer model is first made of the metasurface, which is then sliced into 2D patterns. These patterns are projected via a UV digital micromirror device (DMD), focused onto the surface of a photosensitive monomer, which cures under UV exposure. The cured layer, in the shape of the 2D slice pattern, is then lowered to resupply liquid resin on its surface. The pattern projection is repeated to form the subsequent layer. To expand the scalability of the printed metasurface, the projection system moves with an optical scanning system to project patterns on multiple areas of the liquid surface, producing large scale metasurfaces with a microscale resolution. This process is illustrated in Fig. 2(a). The fabricated sample has a width of 3.43 mm, is made up of 6 periods, and is thus 61.12 mm in length, as seen in Fig. 2(b). The material used for fabrication was a custom formulated UV-cured 1,6-hexanediol diacrylate polymer, with a low dosage of photoabsorber (Sudan 1; CAS 842–07-9), which has a Young's modulus of 512 MPa and a density of 1.1 g/cm3. The smallest features in this sample are the walls of the unit cells that are 86 μm in thickness. However, it must be noted that if required, the resolution offered by this fabrication technique can be further improved by utilizing higher focusing power of the reduction lenses and lowering the sensitivity of the photocurable polymer used.38 The fabricated minuscule prototype is then experimentally validated in a 2D waveguide of height 6 mm, to ensure that only the fundamental mode can propagate inside [this is, hence, also the depth of the third dimension of the sample as seen in Fig. 2(b)]. The waveguide shown in Fig. 2(c) is made with laser-cut acrylic plates to confine the transmitted ultrasonic wave in a quasitwo-dimensional space. The entire experimental setup is essentially a scaled-down version of the scan stage used for audible sound in previous works.7,8,12,39 It would thus be ideal to use a speaker array to generate a plane wave with a Gaussian amplitude profile to impinge on the metasurface. In this experimental setup, for the sake of convenience, two single frequency (40 kHz) Murata speakers (MA40S4S) were employed as the source—fitted with tapered horn-type waveguide adapters, to generate a quasiplane wave in the far field. Perfectly matched layers (PMLs) were used in the simulations to minimize reflections from the boundaries, which were imitated in experiments by using absorbing foam on the two sides. A two-dimensional linear scan stage was programmed to map the field on the transmissive side of the metasurface. This is done by translating a glass tube of diameter 1 mm over a rectangular region, with a step size of 2 mm. It should be noted that the scattering due to the tube is negligible as it is much smaller than the wavelength of the propagating sound. At every step, the glass tube guides the sound to a microphone (Murata MA40S4R). An LM358-based operational amplifier was used as the preamplification system, from where the signal is transmitted to an NI PCI-6251 data-acquisition board. The collected time-domain signals were then Fourier transformed to the frequency domain to generate the complex field pattern at 40 kHz. Figure 3 shows the results from the measurement in comparison to two sets of full wave simulations—with and without dissipation. The case without dissipation shown in Fig. 3(a) reaffirms the scalability of the GIM: similar wavefront bending is seen for both 3.43 kHz (top) and 40 kHz (bottom). However, when the thermoviscous effects are incorporated, the effect of loss is more significant in the case of 40 kHz, as discussed previously. The normalized wave field in the case with dissipation in Fig. 3(a) is then compared with our measurement results shown in Fig. 3(b). A reasonable agreement is seen between the numerical and experimental results, for both acoustic pressure and intensity fields. Minor deviations can be attributed to fabrication defects, measurement errors, and acoustic structure-interactions25 that are not taken into account in the numerical simulations.

FIG. 2.

(a) Schematic illustration of the fabrication process. (b) The ultrasound metasurface prototype with a thickness of 3.43 mm and a length of 61.70 mm. The inset shows the zoom of the sample, where the shunted Helmholtz resonators and straight channels can be seen. (c) The scaled down experimental set-up of the two-dimensional linear scan stage. Two ultrasonic speakers are connected to tapered horn-type wave guide adapters, which are shown on the right (red), to ensure the incident wavefront. The tube (black) on the left maps the scan area and guides the sound to a receiver.

FIG. 2.

(a) Schematic illustration of the fabrication process. (b) The ultrasound metasurface prototype with a thickness of 3.43 mm and a length of 61.70 mm. The inset shows the zoom of the sample, where the shunted Helmholtz resonators and straight channels can be seen. (c) The scaled down experimental set-up of the two-dimensional linear scan stage. Two ultrasonic speakers are connected to tapered horn-type wave guide adapters, which are shown on the right (red), to ensure the incident wavefront. The tube (black) on the left maps the scan area and guides the sound to a receiver.

Close modal
FIG. 3.

(a) Simulated acoustic fields (Pascal), which show wavefront bending at both 3.4 kHz (top) and 40 kHz (bottom), for cases with (right) and without (left) dissipation. (b) Comparison of the normalized acoustic pressure and intensity fields between simulation and measurement.

FIG. 3.

(a) Simulated acoustic fields (Pascal), which show wavefront bending at both 3.4 kHz (top) and 40 kHz (bottom), for cases with (right) and without (left) dissipation. (b) Comparison of the normalized acoustic pressure and intensity fields between simulation and measurement.

Close modal

In conclusion, we have designed, fabricated, and experimentally characterized a minuscule HRAGIM that operates at 40 kHz. The scalability of such a design is demonstrated, and the role of thermoviscosity is discussed. It is clearly seen that dissipation has a greater effect on transmission at higher frequencies. However, its effect on phase is relatively negligible, and wavefront modulation can therefore be realized. Although the presence of thermoviscous dissipation can act as a limitation for transmissive applications such as bending and focusing, it can be fruitful to engineer compact devices for tunable asymmetric transmission7,26 and sound absorption,27 among others.11,40 It is hoped that this study will bring about exotic possibilities to the research in acoustic metasurfaces, especially in miniaturized acoustic devices. Such a metasurface based design and their realization through emerging additive manufacturing techniques can hence be readily scaled down to operate at much higher frequencies, to find applications as compact acoustic devices for sensing, levitation, noncontact ultrasonic imaging, and therapeutic ultrasound.

See the supplementary material for diffractive efficiencies of HRAGIMs scaled to operate at other frequencies, relevant simulation results, and the functionality of the 40 kHz HRAGIM at nearby frequencies.

Chen Shen, Yangbo Xie, and Steven Cummer were supported by a Multidisciplinary University Research Initiative grant from the Office of Naval Research (Grant No. N00014-13-1-0631) and an Emerging Frontiers in Research and Innovation grant from the U.S. National Science Foundation (Grant No. 16-41084). Huachen Cui and Xiaoyu Zheng were supported by a funding from the Air Force Office of Scientific Research (No. FA9550-18-1-0299).

1.
B.
Assouar
,
B.
Liang
,
Y.
Wu
,
Y.
Li
,
J. C.
Cheng
, and
Y.
Jing
,
Nat. Rev. Mater.
3
,
460
(
2018
).
2.
H.
Ge
,
M.
Yang
,
C.
Ma
,
M. H.
Lu
,
Y. F.
Chen
,
N.
Fang
, and
P.
Sheng
,
Natl. Sci. Rev.
5
,
159
(
2018
).
3.
S. A.
Cummer
,
J.
Christensen
, and
A.
Alù
,
Nat. Rev. Mater.
1
,
1
(
2016
).
4.
Y.
Zhu
,
X.
Fan
,
B.
Liang
,
J.
Cheng
, and
Y.
Jing
,
Phys. Rev. X
7
,
1
(
2017
).
5.
C.
Shen
,
J.
Xu
,
N. X.
Fang
, and
Y.
Jing
,
Phys. Rev. X
4
,
1
(
2014
).
6.
G.
Memoli
,
M.
Caleap
,
M.
Asakawa
,
D. R.
Sahoo
,
B. W.
Drinkwater
, and
S.
Subramanian
,
Nat. Commun.
8
,
14608
(
2017
).
7.
Y.
Li
,
C.
Shen
,
Y.
Xie
,
J.
Li
,
W.
Wang
,
S. A.
Cummer
, and
Y.
Jing
,
Phys. Rev. Lett.
119
,
035501
(
2017
).
8.
J.
Li
,
C.
Shen
,
A.
Díaz-Rubio
,
S. A.
Tretyakov
, and
S. A.
Cummer
,
Nat. Commun.
9
,
1
(
2018
).
9.
Z.
Tian
,
C.
Shen
,
J.
Li
,
E.
Reit
,
Y.
Gu
,
H.
Fu
,
S. A.
Cummer
, and
T. J.
Huang
,
Adv. Funct. Mater.
29
,
1808489
(
2019
).
10.
C.
Shen
,
Y.
Xie
,
J.
Li
,
S. A.
Cummer
, and
Y.
Jing
,
J. Appl. Phys.
123
,
124501
(
2018
).
11.
T.
Liu
,
S.
Liang
,
F.
Chen
, and
J.
Zhu
,
J. Appl. Phys.
123
,
091702
(
2018
).
12.
Y.
Xie
,
W.
Wang
,
H.
Chen
,
A.
Konneker
,
B. I.
Popa
, and
S. A.
Cummer
,
Nat. Commun.
5
,
1
(
2014
).
13.
W.
Wang
,
Y.
Xie
,
B. I.
Popa
, and
S. A.
Cummer
,
J. Appl. Phys.
120
,
195103
(
2016
).
14.
Y.
Li
,
S.
Qi
, and
M. B.
Assouar
,
New J. Phys.
18
,
043024
(
2016
).
15.
Y.
Li
,
X.
Jiang
,
B.
Liang
,
J. C.
Cheng
, and
L.
Zhang
,
Phys. Rev. Appl.
4
,
1
(
2015
).
16.
Y.
Xie
,
A.
Konneker
,
B. I.
Popa
, and
S. A.
Cummer
,
Appl. Phys. Lett.
103
,
1
(
2013
).
17.
Y.
Li
,
B.
Liang
,
X.
Tao
,
X. F.
Zhu
,
X. Y.
Zou
, and
J. C.
Cheng
,
Appl. Phys. Lett.
101
,
233508
(
2012
).
18.
K.
Song
,
J.
Kim
,
S.
Hur
,
J.-H.
Kwak
,
S.-H.
Lee
, and
T.
Kim
,
Sci. Rep.
6
,
32300
(
2016
).
19.
J.
Qian
,
J.-P.
Xia
,
H.-X.
Sun
,
S.-Q.
Yuan
,
Y.
Ge
, and
X.-Z.
Yu
,
J. Appl. Phys.
122
,
244501
(
2017
).
20.
C.
Shen
,
A.
Díaz-Rubio
,
J.
Li
, and
S. A.
Cummer
,
Appl. Phys. Lett.
112
,
183503
(
2018
).
21.
G. Y.
Song
,
Q.
Cheng
,
T. J.
Cui
, and
Y.
Jing
,
Phys. Rev. Mater.
2
,
065201
(
2018
).
22.
Y.
Xie
,
C.
Shen
,
W.
Wang
,
J.
Li
,
D.
Suo
,
B.-I.
Popa
,
Y.
Jing
, and
S. A.
Cummer
,
Sci. Rep.
6
,
35437
(
2016
).
23.
M.
Molerón
,
M.
Serra-Garcia
, and
C.
Daraio
,
New J. Phys.
18
,
033003
(
2016
).
24.
X.
Jiang
,
Y.
Li
, and
L.
Zhang
,
J. Acoust. Soc. Am.
141
,
EL363
(
2017
).
25.
N. J.
Gerard
,
Y.
Li
, and
Y.
Jing
,
J. Appl. Phys.
123
,
124905
(
2018
).
26.
F.
Ju
,
Y.
Tian
,
Y.
Cheng
, and
X.
Liu
,
Appl. Phys. Lett.
113
,
121901
(
2018
).
27.
C.
Shen
and
S. A.
Cummer
,
Phys. Rev. Appl.
9
,
54009
(
2018
).
28.
L.
Zigoneanu
,
B. I.
Popa
, and
S. A.
Cummer
,
Phys. Rev. B
84
,
1
(
2011
).
29.
B.
Liu
,
J.
Zhao
,
X.
Xu
,
W.
Zhao
, and
Y.
Jiang
,
Sci. Rep.
7
,
1
(
2017
).
30.
Y.
Xie
,
Y.
Fu
,
Z.
Jia
,
J.
Li
,
C.
Shen
,
Y.
Xu
,
H.
Chen
, and
S. A.
Cummer
,
Sci. Rep.
8
,
1
(
2018
).
31.
A.
Marzo
,
S. A.
Seah
,
B. W.
Drinkwater
,
D. R.
Sahoo
,
B.
Long
, and
S.
Subramanian
,
Nat. Commun.
6
,
1
(
2015
).
32.
G. T.
Clement
,
H.
Nomura
,
H.
Adachi
, and
T.
Kamakura
,
Phys. Med. Biol.
58
,
6263
(
2013
).
33.
X.
Zheng
,
W.
Smith
,
J.
Jackson
,
B.
Moran
,
H.
Cui
,
D.
Chen
,
J.
Ye
,
N.
Fang
,
N.
Rodriguez
,
T.
Weisgraber
, and
C. M.
Spadaccini
,
Nat. Mater.
15
,
1100
(
2016
).
34.
X.
Zheng
,
H.
Lee
,
T. H.
Weisgraber
,
M.
Shusteff
,
J.
Deotte
,
E. B.
Duoss
,
J.
Kuntz
,
M. M.
Biener
,
Q.
Ge
,
J. A.
Jackson
,
S. O.
Kucheyev
,
N. X.
Fang
, and
C. M.
Spadaccini
,
Science
344
,
1373
(
2014
).
35.
N.
Yu
,
P.
Genevet
,
M. A.
Kats
,
F.
Aieta
,
J.-P.
Tetienne
,
F.
Capasso
, and
Z.
Gaburro
,
Science
334
,
333
(
2011
).
36.
Y.
Li
,
B.
Liang
,
Z.-M.
Gu
,
X.-Y.
Zou
, and
J.-C.
Cheng
,
Sci. Rep.
3
,
2546
(
2013
).
37.
L. E.
Kinsler
,
A. R.
Frey
,
A. B.
Coppens
, and
J. V.
Sanders
,
Fundamentals of Acoustics
, 4th ed. (
John Wiley & Sons
,
2012
).
38.
X.
Zheng
,
J.
Deotte
,
M. P.
Alonso
,
G. R.
Farquar
,
T. H.
Weisgraber
,
S.
Gemberling
,
H.
Lee
,
N.
Fang
, and
C. M.
Spadaccini
,
Rev. Sci. Instrum.
83
,
125001
(
2012
).
39.
C.
Shen
,
Y.
Xie
,
J.
Li
,
S. A.
Cummer
, and
Y.
Jing
,
Appl. Phys. Lett.
108
,
2
(
2016
).
40.
Y.
Zhu
,
J.
Hu
,
X.
Fan
,
J.
Yang
,
B.
Liang
,
X.
Zhu
, and
J.
Cheng
,
Nat. Commun.
9
,
1
(
2018
).

Supplementary Material