In recent years, there have been many interests in studying electromagnetic superscattering because of its special applications^[1–11] in antenna design, optical sensing and detection, bioimaging, optical tagging to spectroscopy, scattering-particle-assisted passive radiation cooling, and photovoltaic devices. To realize superscattering for a subwavelength object, complementary metamaterials based on transformation optics have been proposed initially to drastically enhance the scattering cross section of a particle^[12,13]. This approach requires both electrical and magnetic anisotropic inhomogeneous parameters. In order to relax the requirement, Ruan and Fan proposed a subwavelength object with a plasmonic-dielectric-plasmonic layer structure, and numerically demonstrated that the total scattering cross section of a subwavelength structure could far exceed the single-channel limit when multiple plasmonic resonant modes were aligned^[14]. Later, many theoretical works have been performed to generate superscattering by various structures^[14–23]. Subwavelength multilayer structures have become widely used to overcome the single-channel limit by exciting degenerated confined waves, such as localized surface plasmons in metal^[14] or graphene^[15–17] and multipolar Mie resonances in anisotropic high index dielectrics^[18–20]. Recently, a low or near zero refractive index environment has been proposed as a novel mechanism to enhance the scattering from tiny objects, which has been proved in some works^[24,25].

All of the above-mentioned superscattering phenomena exhibit a single-frequency regime. Multifrequency superscattering may enable improved sensitivity and accuracy for extensive applications ranging from sensing, bioimaging, and optical tagging to spectroscopy^[26,27]. Only very recently was a multifrequency superscatterer first, to the best of our knowledge, proposed by using a subwavelength hyperbolic structure^[28]. Since then, multifrequency superscattering from a deep-subwavelength spoof plasmonic structure has been also discussed^[29]. Moreover, a subwavelength metasurface-based all-metal or metal dielectric structures with low material losses has also been proposed to experimentally achieve multifrequency superscattering by degenerate resonances of confined surface waves^[30,31]. Very recently, a graphene-based dual-band superscatterer at mid-infrared frequencies is proposed, and the phenomena of dual-band superscattering are tunable^[32]. However, the obtained multifrequency superscattering arises mainly from the constructive interference of scattering fields from different channels in the forward direction. In this paper, we have achieved simultaneously, within a graphene-based subwavelength structure, multifrequency superscattering and superscattering shaping with different engineered scattering patterns induced by the interferences of multipoles resonantly overlapped. Due to multiple polaritonic dispersion lines in a graphene-based subwavelength structure, the phenomena of multifrequency superscattering are found to have a high tolerance to some structural variations and realistic graphene loss. In contrast with multilayer plasmonic structures, the performance of the proposed structure can be easily refined for the desired frequency by adjusting the chemical potential and loss of the graphene layer.

In this work, we consider a two-dimensional subwavelength cylindrical structure. The structure consists of six layers of dielectric material, each of which is coated with a graphene shell, as shown in Fig. 1(a). Such a six-layer subwavelength structure is the main focus of this analysis, because six layers is sufficient for multifrequency superscattering and superscattering shaping. The conductivity of graphene can be expressed as the Kubo formula: ${\sigma }_g(\omega ,{\mu }_c,\tau ,T)={\sigma }_{\mathrm{intra}}+{\sigma }_{\mathrm{inter}}$^[17]. The TM-polarized plane wave with the wave vector $k$ along the $x$ direction and magnetic field along the $z$ direction is incident from the left side of the cylinder structure. Based on the transfer matrix approach of multilayered cylindrical structures^[32–34], the incident and scattered magnetic field in the outermost surrounding air region can be expressed as $${\stackrel{\rightharpoonup }{H}}_{\text{inc}}=\stackrel{\rightharpoonup }{z}{H}_0\sum _{m=-\infty }^{\infty }{i}^m{J}_m(k_b\rho )e^{im\varphi },$$(1)$${\stackrel{\rightharpoonup }{H}}_{\mathrm{sca}}=\stackrel{\rightharpoonup }{z}{H}_0\sum _{m=-\infty }^{\infty }{i}^m{s}_m{H}_m^{(1)}(k_b\rho )e^{im\varphi },$$(2)where $H_0$ is the amplitude of the incident field, $\varphi$ and $\rho$ are the azimuthal angle and radial component, respectively, and $J_m$ and $H_m^{(1)}$ denote the Bessel function and the Hankel function of the first kind with order $m$, respectively. $k_b=\omega /c$ is the wavenumber in free space, and $s_m$ is the scattering coefficient with the $m$th angular momentum channel. The magnetic field in the subwavelength cylindrical structure for ${\rho }_{n-1}<\rho <{\rho }_n(n=2,3,\dots ,6)$ and $\rho <{\rho }_1$ can be respectively written as $${\stackrel{\rightharpoonup }{H}}_n=\stackrel{\rightharpoonup }{z}{H}_0\sum _{m=-\infty }^{+\infty }{i}^m[b_{m,n}{J}_m(k_n\rho )+c_{m,n}{H}_m^{(1)}(k_n\rho )]e^{im\varphi },$$(3)$${\stackrel{\rightharpoonup }{H}}_1=\stackrel{\rightharpoonup }{z}{H}_0\sum _{m=-\infty }^{\infty }{i}^m{a}_m{J}_m(k_1\rho )e^{im\varphi },$$(4)where ${\rho }_1$ is the radius of the innermost layer of the dielectric, $k_n=\omega /c\sqrt{{\epsilon }_n}$ is the wavenumber in the $n$th layer of the dielectric medium, and $b_{m,n}$, $c_{m,n}$, and $a_m$ are the magnetic-field coefficients with order $m$. According to Maxwell’s equation, the corresponding electric field can be calculated by $\stackrel{\rightharpoonup }{E}=\frac{i}{\omega \epsilon {\epsilon }_0}\nabla \times \stackrel{\rightharpoonup }{H}$. The electric and magnetic field at the interface $\rho ={\rho }_i(i=1,2,\dots ,6)$ can be restrained by the following boundary conditions for the graphene shells between adjacent layers of dielectric^[35]: $$\stackrel{\rightharpoonup }{n}\times ({\stackrel{\rightharpoonup }{E}}_{n+1,\varphi }-{\stackrel{\rightharpoonup }{E}}_{n,\varphi })=0,$$(5)$$\stackrel{\rightharpoonup }{n}\times ({\stackrel{\rightharpoonup }{H}}_{n+1,z}-{\stackrel{\rightharpoonup }{H}}_{n,z})=\stackrel{\rightharpoonup }{J}.$$(6)

Here, $J=\sigma (\omega )E_{\varphi }$ is the surface current of graphene satisfying Ohm’s law. To derive the scattering coefficient $s_m$, dynamical matrix $D_{m,n}$ and $M_{m,n+1}$ are defined as $$D_{m,n}({\rho }_n)=\left[\begin{array}{cc}{J}_m(k_n{\rho }_n)& H_m^{(1)}(k_n{p}_n)\\ \frac{\sqrt{{\epsilon }_n}{J}_m^{\prime }(k_n{\rho }_n)}{{\epsilon }_n}& \frac{\sqrt{{\epsilon }_n}{H}_m^{\prime (1)}(k_n{\rho }_n)}{{\epsilon }_n}\end{array}\right],$$(7)$$M_{m,n+1}({\rho }_n)=\left[\begin{array}{cc}{J}_m(k_{n+1}{\rho }_n)-i{\sigma }_n{\eta }_1& H_m^{(1)}(k_{n+1}{p}_n)-i{\sigma }_n{\eta }_2\\ \frac{\sqrt{{\epsilon }_{n+1}}{J}_m^{\prime }(k_{n+1}{\rho }_n)}{{\epsilon }_{n+1}}& \frac{\sqrt{{\epsilon }_{n+1}}{H}_m^{\prime (1)}(k_{n+1}{\rho }_n)}{{\epsilon }_{n+1}}\end{array}\right],$$(8)where ${\eta }_1=k_{n+1}{J}_m^{\prime }(k_{n+1}{\rho }_n)/(\omega {\epsilon }_0{\epsilon }_{n+1})$ and ${\eta }_2=k_{n+1}{H}_m^{\prime (1)}·(k_{n+1}{\rho }_n)/(\omega {\epsilon }_0{\epsilon }_{n+1})$. By substituting Eqs. (3)–(6) into Eqs. (7) and (8), the coefficients $b_{m.n}$ and $c_{m,n}$ can be connected to $b_{m,n+1}$ and $c_{m,n+1}$ by the defined dynamical matrix $D_{m,n}$ and $M_{m,n+1}$ as follows: $$D_{m,n}({\rho }_n)\left[\begin{array}{l}{b}_{m,n}\hfill \\ c_{m,n}\hfill \end{array}\right]=M_{m,n+1}({\rho }_n)\left[\begin{array}{l}{b}_{m,n+1}\hfill \\ c_{m,n+1}\hfill \end{array}\right].$$(9)

Then, from the recurrence relation of the coefficients of the adjacent layers, the relation between the coefficients of the innermost layer and the outermost layer can be described by the following formula: $$\left[\begin{array}{cc}{J}_m(k_1{\rho }_1)& 0\\ \frac{\sqrt{{\epsilon }_1{J}_m^{\prime }(k_1{\rho }_1)}}{{\epsilon }_1}& 0\end{array}\right]\left[\begin{array}{c}{a}_m\\ 0\end{array}\right]=T_m\left[\begin{array}{c}1\\ s_m\end{array}\right],$$(10)where $$\begin{gathered} T_m=M_{m,2}({\rho }_1)\{\underset{n=2}{\overset{6}{\mathrm{\Pi }}}[D_{m,n}({\rho }_n){]}^{-1}\cdot M_{m,n+1}({\rho }_n)\} \\ =\left[\begin{array}{cc}{T}_{m,11}& T_{m,12}\\ T_{m,21}& T_{m,22}\end{array}\right]. \end{gathered}$$(11)

Next, the scattering coefficient obtained according to Eqs. (12) and (13) can be expressed as $$s_m=\frac{x{T}_{m,21}-y{T}_{m,11}}{y{T}_{m,12}-x{T}_{m,22}},$$(12)where $x=J_m(k_1{\rho }_1)$, $y=\sqrt{{\epsilon }_1}{J}_m^{\prime }(k_1{\rho }_1)/{\epsilon }_1$. Finally, the normalized scattering cross section (NSCS) can be calculated by $$\mathrm{NSCS}=\sum _{m=-\infty }^{+\infty }{|s_m|}^2.$$(13)

Here, NSCS is normalized by $2\lambda /\pi$, where $\lambda$ is the wavelength of the incident wave. According to the electromagnetic multipole analysis of a cylindrical structure^[20,36,37], for TM-polarized plane waves, $s_0,s_{\pm 1},s_{\pm 2},s_{\pm 3}$, and $s_{\pm 4}$ correspond to the magnetic dipole (MD), electric dipole (ED), electric quadrupole (EQ), magnetic octupole (MO), and electric octupole (EO), respectively. Theoretically, the phenomenon of superscattering can be realized by making the NSCS exceed the single-channel limit. To enlarge the scattering cross section, material and geometrical parameters can be tuned so that different electromagnetic multipoles simultaneously appear at the same frequency. In this system, contributions of ED, EQ, MO, and EO are mainly considered because the contributions of other modes are so small that they could be ignored. Significantly, resonant overlapping of different multipoles can not only realize superscattering but also result in different scattering angular distributions.

To highlight the underlying physical properties, we first consider a lossless subwavelength cylindrical structure by setting $1/\tau \to 0$, where $\tau$ is the relaxation time of charge carriers of graphene. The scattering spectra (NSCS) as functions of the frequency for the structure in Fig. 1(a) are shown in Fig. 1(b). Both the total scattering spectrum and the partial spectra from channels with $m=0$, $m=\pm 1$, $m=\pm 2$, $m=\pm 3$, and $m=\pm 4$ are plotted. The geometrical parameters are ${\rho }_1=1493\mathrm{nm}$ and $d_n={\rho }_n-{\rho }_{n-1}=500\mathrm{nm}(n=2,3,\dots ,6)$. Relative permittivities of dielectrics are ${\epsilon }_1=1.22$, ${\epsilon }_2=1.21$, ${\epsilon }_3=2.06$, ${\epsilon }_4=1.15$, ${\epsilon }_5=1.21$, and ${\epsilon }_6=2.3$. For all graphene shells, the chemical potential is ${\mu }_c=0.72\mathrm{eV}$ at temperature $T=300\mathrm{K}$, and the thickness can be negligible. In the system, resonances from channels with $m=\pm 1$, $m=\pm 2$, $m=\pm 3$, and $m=\pm 4$ are mainly considered because the scattering cross sections from other channels are too small to be ignored. At 3.79715 THz, the total scattering cross section mostly comes from resonant channels with $m=\pm 1$ and $m=\pm 2$, and the scattering cross sections from these four channels almost all reach the signal-channel limit, which makes the total reach four times the signal-channel limit. The total scattering cross section at 9.77464 THz exceeds four times the signal-channel limit, which originates from a combination of the resonance peaks from channels with $m=\pm 2$ and $m=\pm 3$. Scattering cross sections of these peaks almost all reach signal-channel limit, and the contributions from channels with $m=\pm 1$ are negligible. Similarly, the resonance peaks at 12.65616 THz from channels with $m=\pm 1$ and $m=\pm 3$ are all nearly equal to the signal-channel limit, and they are overlapped well, which results in the total scattering cross section reaching four times the signal-channel limit. The resonance peaks with $m=\pm 2$ and $m=\pm 4$ are overlapped at 14.74908 THz, and these peaks are almost all equal to the signal-channel limit. Obviously, the contributions from channels with $m=\pm 1$ cannot be neglected, so the total scattering cross section is greater than four times the signal-channel limit at 14.74908 THz. This implies that multifrequency superscattering can be realized by overlapping electric and magnetic resonances with different orders.

The underlying mechanism of the resonant overlapping of multipoles in the subwavelength cylinder structure can be qualitatively understood from the dispersion of plasmon polaritons supported by a corresponding planar structure in Fig. 2. Starting from resonance modes supported in the corresponding planar structure, the resonances of the subwavelength cylinder structure can be understood using the whispering gallery condition based on the Bohr model^[14,17]. More specifically, the subwavelength cylinder supports the $m$th-order resonance at a frequency $f$, when the propagation constant $\beta$ of all the modes for the corresponding planar structure satisfies the condition: $\beta 2\pi R_{\mathrm{eff}}=2\pi m$, where $R_{\mathrm{eff}}=\sum _{i=1}^6{\rho }_i/6$ is the effective radius of the cylindrical structure. The dispersion curves of the equivalent planar structure are analyzed by the reflection pole method, which suggests that the surface modes of a layered medium are observed in the generalized reflection coefficient. It is performed after obtaining the scattering parameters from the transmission matrix^[32,38]. There are five branches in the dispersion lines due to the mode coupling of the graphene plasmons. Four black vertical solid lines $|m|=1$, $|m|=2$, $|m|=3$, and $|m|=4$ correspond to different modes ED, EQ, MO, and EO, respectively. Obviously, there are plenty of resonant frequencies of the multiple plasmonic modes in the dispersion lines. Here, only resonant frequencies with overlapped multipoles are analyzed. Two resonances at 3.79715 THz and 12.65616 THz from channel $|m|=1$ are observed, and the positions of the corresponding resonance modes are marked by ${\mathrm{ED}}_1$ and ${\mathrm{ED}}_2$. Three resonances at 3.79715 THz, 9.77464 THz, and 14.74908 THz from channel $|m|=2$ are observed, and the positions of the corresponding resonance modes are marked by ${\mathrm{EQ}}_1,{\mathrm{EQ}}_2$, and ${\mathrm{EQ}}_3$. It is observed that a coincidence of ${\mathrm{ED}}_1$ and ${\mathrm{EQ}}_1$ at the frequency of 3.79715 THz interprets the combination of resonances from channels with $m=\pm 1$ and $m=\pm 2$ in Fig. 1(b). Two resonances at 9.77464 THz and 12.65616 THz from channel $|m|=3$ are observed, and the positions of the corresponding resonance modes are marked by ${\mathrm{MO}}_1$ and ${\mathrm{MO}}_2$. It is clear that ${\mathrm{EQ}}_2$ and ${\mathrm{MO}}_1$ appear close to the frequency of 9.77464 THz, which interprets the combination of resonances from channels with $m=\pm 2$ and $m=\pm 3$ in Fig. 1(b). Similarly, resonances from channels with $m=\pm 1$ and $m=\pm 3$ can be overlapped at the frequency of 12.65616 THz because ${\mathrm{ED}}_2$ and ${\mathrm{MO}}_2$ appear close to the corresponding frequency. One resonance at 14.74908 THz from channel $|m|=4$ is observed, and the position of the corresponding resonance mode is marked by ${\mathrm{EO}}_1$. Clearly, ${\mathrm{EQ}}_3$ and ${\mathrm{EO}}_1$ appear close to the frequency of 14.74908 THz, and this means that the resonances from channels with $m=\pm 2$ and $m=\pm 4$ can be overlapped at the corresponding frequency. Therefore, multifrequency superscattering with resonantly overlapped multipoles can be well understood through the dispersion lines.

Figure 3 intuitively shows the performance of the superscattering pattern shaping in the form of far and near fields for the cylinder structure in Fig. 1(a) at the four superscattering resonant frequencies. Far-field angular scattering amplitude $\mathrm{\Gamma }(\varphi )$ can be analyzed by the derived formula^[20]: $$\mathrm{\Gamma }(\varphi )=\sqrt{2/\pi k}|s_0+2\sum _{m=1}^{\infty }{s}_m\cos (m\varphi )|.$$(14)

Theoretically, various scattering pattern shaping can also be realized by resonantly overlapping different electromagnetic multipoles because these multipoles show different phase symmetries and parities^[20,39,40]. Specifically, the scatterings exhibit opposite parities for multipoles of adjacent orders and the same nature (e.g.,  ED and EQ exhibit odd and even parity, respectively) or for multipoles of different natures but the same order (e.g.,  MD and ED exhibit even and odd parity, respectively).

Figure 3(a) shows the superscattering pattern shaping by resonantly overlapping ED and EQ at the resonant frequency of 3.79715 THz. Far-field scattering angular distribution is shown by a blue solid curve, and the scattering near field $|H_z^{\mathrm{sca}}|$ is shown in the background pattern, which is in accord with the blue solid curve. As shown in Fig. 3(a), the scattering is enhanced in the forward direction but suppressed in the backward direction, because the scatterings of ED and EQ exhibit odd and even parities, respectively. The two multipoles are in phase in the forward direction and out of phase in the backward direction. In fact, in the forward direction, the two resonantly overlapped multipoles are always in phase according to the optical theorem^[20,41]. Additionally, the scattering is also enhanced from 90° to 135° and 225° to 270° in the backward half-scattering circle and suppressed from 45° to 90° and 270° to 315° in the forward half-scattering circle. This can be understood analytically by $\mathrm{\Gamma }(\varphi )\propto |\cos \varphi +\cos (2\varphi )|$ according to Eq. (16) with purely overlapping ED and EQ. The ideal scattering angular distribution with solely overlapped ED and EQ is plotted by red crosses, which is consistent with the blue solid curve from the contributions of other multipoles, except ED and EQ are so small that they can be neglected. For the near field, the total magnetic field $H_z$ distribution at 3.79715 THz is shown in Fig. 3(b). The cylinder structure leaves a large “shadow” in the forward direction, and the size of the “shadow” in the forward direction is much larger than the diameter of the cylindrical structure. This arises from the fact that the radiated magnetic fields of resonantly overlapped multipoles interfere destructively with respect to the incident field in the forward direction, as is required by the optical theorem^[39,41]. Clearly, the near field $H_z$ matches well with the scattering field in Fig. 3(a).

Figure 3(c) illustrates the superscattering pattern shaping with resonantly overlapped EQ and MO by showing the scattering-field distribution at the frequency of 9.77464 THz. The far-field scattering angular distribution is plotted by a blue solid curve, which matches well with the scattering near field shown in the background pattern. Like Fig. 3(a), the scattering is enhanced in the forward direction but suppressed in the backward direction due to opposite parities of EQ and MO (even and odd, respectively). The two overlapped multipoles are in phase and out of phase, respectively, in the forward and backward directions. It is observed that the scattering is also suppressed from 30° to 45° and 315° to 330° in the forward half-scattering circle, which can be understood analytically by $\mathrm{\Gamma }(\varphi )\propto |\cos (2\varphi )+\cos (3\varphi )|$ when EQ and MO are only overlapping. The red crosses plot the ideal scattering pattern with the two modes purely overlapped. The slight discrepancy between blue solid curve and red crosses originates from the fact that the contribution of the ED is negligible. For the near field, the total field $H_z$ at 9.77464 THz is shown in Fig. 3(d). A remarkable “shadow” can also be seen in the forward direction, and the “shadow” is much larger than the diameter of the cylindrical structure, which agrees well with the scattering field in Fig. 3(c).

Figure 3(e) illustrates the superscattering pattern shaping with resonantly overlapped ED and MO by showing the scattering-field distribution at the frequency of 12.65616 THz. The scattering near and far fields correspond to the background pattern and the blue solid curve, respectively, and they are well matched. Clearly, the scattering is enhanced symmetrically both in the forward and backward directions due to the same odd parity of the scatterings of ED and MO. Additionally, the interferences of overlapped ED and MO can also suppress the scattering at $\varphi =45°$, $\varphi =90°$, $\varphi =135°$, $\varphi =225°$, $\varphi =270°$, and $\varphi =315°$ due to $\mathrm{\Gamma }(\varphi )\propto |\cos \varphi +\cos (3\varphi )|$. When only considering the overlapped ED and MO, the ideal scattering angular distribution plotted by red crosses matches well with the blue solid curve because other contributions except for ED and MO are too small to be neglected. For the near field, the total magnetic field at the resonant frequency of 12.65616 THz is shown in Fig. 3(f), and the cylinder structure leaves a large “shadow” larger than the structure diameter, not only in the forward direction but also in backward direction, which is consistent with the scattering field in Fig. 3(e).

Figure 3(g) shows the superscattering pattern shaping by resonantly overlapping EQ and EO at the frequency of 14.74908 THz in the form of the scattering far field plotted by the blue solid curve and the scattering near field shown in the bottom pattern. The scattering amplitudes are symmetrically enhanced both in forward and backward directions, which is caused by the same even parity of the EQ and EO. It can be also seen that the scattering is suppressed from 22.5° to 45°, 67.5° to 112.5°, 135° to 157.5°, 202.5° to 225°, 247.5° to 292.5°, and 315° to 337.5° according to $\mathrm{\Gamma }(\varphi )\propto |\cos (2\varphi )+\cos (4\varphi )|$. The ideal scattering pattern with purely overlapped ED and EO plotted by red crosses does not coincide well with the blue solid curve because the contribution of the ED cannot be neglected at the frequency of 14.74908 THz. For the near field, at 14.74908 THz, the remarkable “shadow”, which is much larger than the structure diameter, can be seen in both the forward and backward directions in Fig. 3(h), and the total magnetic field matches well with the scattering field shown in Fig. 3(g).

To facilitate the realization of the experiment and practical application, the effects of the material losses and the structural changes are analyzed in detail in Fig. 4. First, considering the optical losses of graphene, as shown in Fig. 4(a), the superscattering with resonantly overlapped multipoles still exists at corresponding resonant frequencies when changing the relaxation time $\tau$, and the NSCSs of corresponding resonance peaks are slightly reduced due to the addition of graphene’s losses. Then, to illustrate the high tolerance to structural design, the total scattering spectra are plotted in Fig. 4(b) by varying thickness $d_1$ of the innermost layer of the dielectric and the corresponding scattering spectra of $d_1=1400\mathrm{nm}$, $d_1=1500\mathrm{nm}$, and $d_1=1600\mathrm{nm}$ are shown by the green, blue, and red solid curves, respectively. Except for $d_1$, the other parameters are the same as those in the structure of Fig. 1(a). As the thickness $d_1$ increases gradually, though the positions of the resonant frequencies will move to the lower frequencies, superscattering with resonantly overlapped multipoles can still be obtained, and the NSCSs at the corresponding resonant frequencies are nearly unchanged. Compared with metal materials, graphene has the advantage of adjustability. Next, the tunability of graphene shells is illustrated by changing the chemical potential ${\mu }_c$. The corresponding scattering spectra of ${\mu }_c=0.66\mathrm{eV}$, ${\mu }_c=0.73\mathrm{eV}$, and ${\mu }_c=0.8\mathrm{eV}$ are shown by the green, blue, and red solid curves in Fig. 4(c), respectively. When the chemical potential increases, the positions of resonant frequencies will move towards the higher frequencies, and the superscattering phenomenon still exists with the NSCSs of corresponding resonance peaks almost unchanged. From Fig. 4, when considering the material losses and changing the structural parameters, superscattering with resonantly overlapped multipoles can be still obtained, which indicates that superscattering for a graphene-based subwavelength structure has a good tolerance to material losses and variations of structural parameters.

In conclusion, we propose that a graphene-based subwavelength structure can be a promising and versatile platform for the demonstration of simultaneous multifrequency superscattering and superscattering shaping with different engineered scattering patterns. Based on Bohr’s model and dispersion engineering, it is further demonstrated that such efficient superscattering pattern shaping is induced by the interferences of multipoles resonantly overlapped. In contrast with multilayer plasmonic structures, the performance of the proposed structure can be easily refined for the desired frequency by adjusting the chemical potential and loss of the graphene layer. Owing to multiple polaritonic dispersion lines in a graphene-based subwavelength structure, the phenomena of multifrequency superscattering are found to have a high tolerance to some structural variations and realistic graphene loss. In addition, more exotic scattering phenomena are expected in such a platform, such as dynamically tunable multifrequency cloaking and other multifrequency scattering-based devices.