Anapole-assisted giant electric field enhancement for surface-enhanced coherent anti-Stokes Raman spectroscopy

The coherent anti-Stokes Raman spectroscopy (CARS) techniques are recognized for their ability to detect and identify vibrational coherent processes down to the single-molecular levels. Plasmonic oligomers supporting full-range Fano-like line profiles in their scattering spectrum are one of the most promising class of substrates in the context of surface-enhanced (SE) CARS application. In this work, an engineered assembly of metallic disk-shaped nanoparticles providing two Fano-like resonance modes is presented as a highly-efficient design of SECARS substrate. We show that the scattering dips corresponding to the double-Fano spectral line shapes are originated from the mutual interaction of electric and toroidal dipole moments, leading to the so-called non-trivial first- and second-order anapole states. The anapole modes, especially the higher-order ones, can result in huge near-field enhancement due to their light-trapping capability into the so-called “hot spots”. In addition, independent spectral tunability of the second Fano line shape is exhibited by modulating the gap distance of the corner particles. This feature is closely related to the electric current loop associated with the corner particles in the second-order anapole state and provides a simple design procedure of an optimum SECARS substrate, where the electric field hot spots corresponding to three involved wavelengths, i.e., anti-Stokes, pump, and Stokes, are localized at the same spatial position. These findings yield valuable insight into the plasmonic substrate design for SECARS applications as well as for other nonlinear optical processes, such as four-wave mixing and multi-photon surface spectroscopy.

www.nature.com/scientificreports/ the dark modes and hybridization of plasmons of different multipolar symmetry. The narrower spectral linewidth compared to the standard plasmon resonances as well as the large local field enhancement provided by the dark modes, make the plasmonic Fano-based structures an ideal platform for diverse applications such as refractive index chemical and biological sensing 24,[28][29][30] , surface-enhanced spectroscopy 31,32 , low-threshold nano-lasers 33 , and novel on-chip photonic device designs 34,35 . Metallic nanoparticle clusters in the form of the trimer 36 , quadrumer 30 , pentamer 24 , heptamer 37 , and higherorder 25 assemblies are one of the most promising ways to generate Fano-like resonance line shapes. The structural and chemical characteristics of the clusters such as size, shape, thickness, and material of the particles present in the cluster play an important role to obtain the desired Fano response. In addition, inter-particle gap distances between the adjacent particles define the coupling strength between LSPR modes of the elementary nano-scatterers.
In terms of the analysis of such atypical line profiles in oligomer plasmonic structures, several different techniques have been conducted over the recent years. Plasmon hybridization theory 38 has been extensively exploited to intuitively describe the interactions in plasmonic systems, where the structure under consideration is subdivided into two or more subsystems with known properties. The modes of each subsystem are added constructively or destructively to form a bonding (symmetric) and an anti-bonding (anti-symmetric) mode 39 , leading to the formation of Fano-like resonance line shape. In addition, multiple studies have been devoted to theoretically investigate and interpret the origin of these anomalous resonances in metallic oligomers, such as the coupled mass oscillator model [40][41][42] ,the coupled-mode formalism 43 , the circuit model of Fano resonance 44 , and the rigorous quantum electro-dynamic formulation 45,46 . These simplistic methods, however, are inadequate to provide a complete scenario of Fano-like spectral profiles in plasmonic oligomers, forcing the use of full-wave theories 45,46 .
Regarding the application, plasmonic Fano-based oligomer structures have been of particular interest to act as a substrate for single molecular detection as well as vibrational imaging. In this context, Halas et al. 47 demonstrated that a plasmonic quadrumer substrate supporting a Fano-like resonance mode can be applied to acquire single-molecule detection sensitivity using the coherent anti-Stokes Raman spectroscopy (CARS) 48,49 . CARS is a specific type of four-wave mixing 50 (FWM) non-linear optical process where two laser beams and a Stokes field coherently interact through the third-order susceptibility of the identified material to get access vibrational properties. To achieve an efficient Raman signal enhancement in a plasmonic SECARS substrate depends on simultaneous enhancement of electric field at the pump, Stokes, and anti-Stokes frequencies as well as having hot spots with the same spatial positions 51 . This indicates individual contribution of each of the three wavelengths in SECARS enhancement factor (EF). He et al. 51 presented a plasmonic asymmetric trimer assembly supporting double-Fano line shapes which brings all of the three CARS frequencies hot spots to the same spatial location. The achieved EF for SECARS signal in this structure is around 2 × 10 11 for the gap distance of g = 12 nm . Recently, Zhang et al. 52 theoretically investigated the plasmonic crisscross dimer exhibiting an electric dipole bright mode and magnetic dipole dark modes, which can be matched with the wavelength of pump, Stokes and anti-Stokes lights. Their proposed design achieves a significant amplification of the CARS signal reaching near a factor of 10 13 for the 8-nm gap size.
Extending the aforementioned idea, we propose and theoretically analyze a metallic octamer structure representing double-Fano spectral line profile, which is constructed by an assembly of thin-film disk-shaped nanoparticles. Compared to the similar proposals based on multi-disks CARS substrates 53 , our design enables efficient generation of a single hot spot governing the CARS emission through a symmetrically centered nanogap. Our theoretical study provides an in-depth analysis of scattering dips in plasmonic oligomers, and their corresponding electric current distribution gives insight to design an optimum substrate for SECARS applications. To begin with, the Cartesian decomposition method in the finite element based commercial software (COMSOL Multiphysics) is exploited to calculate the contribution of the five leading multipoles in the total scattering cross-section (SCS). Decomposition results reveals that the two characteristic dips originate from a destructive interference between electric and toroidal 54-56 dipole moments, indicating the first-and second-order anapole states 57 excited in the system. The effects of the structure size as well as the illumination angle on the evolution of Fano line shapes is then explored. These analyses provide a simple design procedure for an optimum multi-resonance SECARS substrate. In this regard, independent spectral tunability of the first-order anapole state is exhibited by modulating the gap distance of the corner particles. Benefiting from such an independent tunability behavior, three excited resonance modes with large field enhancements can be achieved whose hot spots are spatially overlapping and their wavelengths spectrally match with the wavelengths of pump, Stokes, and anti-Stokes beams. When the inter-disk gap spacing is reduced to g = 8 nm , the near-fields are further enhanced due to the improved energy confinement, and the maximum SECARS EF is approaching to 10 15 . Therefore, it is expected that the proposed oligomer design could be useful for the applications based on biosensing and nonlinear surface spectroscopy.

Results and discussion
Double-Fano resonant structure. Let us first propose a simple plasmonic oligomer which serves as a basis for our engineered SECARS substrate. The structure presented here is a finite triangular lattice of diskshaped gold nanoparticles as shown in Fig. 1. The geometry has lattice constant a = 180 nm and disk radius R = 84 nm which is equivalent to the gap distance of g = 12 nm . The height of the disk-based structure is always kept constant at h = 20 nm , which has little influence on the scattering spectrum of the oligomer. Here, similar to Ref. 58 , we have focused on relatively thin gold films that can be realized benefiting from advances in lithography and self-assembly techniques 59 . In terms of the shape of the particles, nano-disks have several advantages such as easy fabrication, polarization-independent property, and relatively simple control of parameters. Note that although the near field enhancement considerably increases with the reduction of the gap size, the gap dis-tance is not considered sub-10-nm to avoid the proximity effects that deteriorates the field enhancement due to nonlocal response 60 . This undesired effect is due to the significant localization of electromagnetic fields caused by the coupling between the adjacent nanoparticles 61 . The dispersive relative permittivity of gold is taken from the experimental data in Ref. 62 . For simplicity, the surrounding dielectric environment is assumed to be a vacuum (ε d = 1) in all simulations. According to the point made in Ref. 51 , introduction of the dielectric substrates or probe molecules in practice only shifts the scattering dips to longer wavelengths accompanied with a slight increase in linewidth because of dielectric screening. Apparently, deposing the particles on top of a metallic substrate results in the improvement of the performance due to the contribution of image modes 63 .
The total SCS spectrum of our proposed structure is shown in Fig. 2. This scattering spectrum exhibits two characteristic dips at the wavelengths of 1 = 705 nm and 2 = 890 nm representing excitation of the Fanolike resonance modes. To get some insight into the observed Fano line shapes, Cartesian multipolar analysis is implemented through the following paragraphs.
Theoretical investigation of asymmetric double-Fano resonance mode. Generally, when an incident plane wave interacts with a nanoparticle, the re-emitted beam at an observation point located in the far-field zone can be described in terms of radiation from different multipoles. To identify the contribution of each multipole, we exploit decomposition of the total induced fields inside the meta-molecule under-study into approximate Cartesian multipole basis including contributions of the toroidal dipole moments. In the framework of approximate Cartesian decomposition, the contribution of each multipole moment in the total radiation of any particle with arbitrary shape can be calculated by starting from the induced electric currents and then simplifying the special functions with the small argument asymptotic terms in the long-wavelength regime 64 (See Methods).
The multipolar SCS is computed to elucidate the nature of the different optical resonances sustained by the structure, particularly those regard the Fano-like resonance line shapes. Figure 2 includes contribution of the five leading order multipole moments in the total SCS spectrum and their summation obtained by Eq. (17) (See "Methods"). Here, contribution of the Cartesian electric and magnetic dipoles are represented by P and MD, electric and magnetic quadrupoles are denoted with EQ and MQ, respectively, and the toroidal dipoles are labeled by TD. As can be seen in Fig. 2a, the sum of Cartesian multipole contributions is in a good agreement with directly calculated total SCS (See "Computational methods") throughout most parts, except for the region with shorter wavelength. The small discrepancy between these spectra can be explained by the long-wavelength approximations made in Cartesian multipole decomposition 64,65 . To further illustrate, we have employed the multipole expansion method in spherical coordinate, showing the results for the first four modes as well as their summation (See Fig. S1 in Supplementary Information). We see in Fig.2a that the most contributive multipoles to the meta-molecule response are the electric dipole as well as closely matched toroidal dipole and magnetic quadrupole responses. Furthermore, the total interference of the Cartesian electric dipole P and the toroidal dipole TD moments (represented as ED by the green curve) reveals the partial scattering cancellation with Fano-like spectral windows. Therefore, it can be concluded that the observed Fano-like resonances are closely attributed to the anapole excitation, i.e. the destructive interference between the electric and toroidal dipoles due to their identical radiation pattern. Furthermore, anapole states naturally lead to the formation of localized fields with strong enhancements 66 , which opens up new possibilities in spectroscopic and biochemical sensing and quantum nanophotonic applications 67 . Here, the observed scattering dips are referred as the first-and secondorder anapole states [68][69][70][71] , labeled with AM 1 and AM 2 in Fig.2, respectively. To explore origin of the cancellation of P and TD moments observed in ED spectrum, considering dominant components in the direction of applied electric field (x-axis), the phases of Cartesian electric φ P and toroidal dipole moments φ TD and their differences �φ = φ P − φ TD are shown in Fig. 2b. It is visible that in our spectral range of interest, the phase difference is www.nature.com/scientificreports/ very close to π radiant, indicating that P and TD moments are almost out of phase. This fulfills the cancellation observed in Fig. 2a.
To further illustrate excitation of the so-called anapole modes, the surface charge density distributions on the particles at the position of characteristic dips are shown in Fig.2c,d, where the observation plane is considered in the mid-height of the structure. According to the dipolar moment orientation of the particles (dark arrows in figure), the AM 1 mode of our proposed oligomer shows dual displacement electric current loops within the structure. These moments are aligned head-to-tail and form a pair of out-of-plane counter-oriented magnetic dipoles, which corresponds to the superposition of a magnetic quadrupole moment and a head-to-tailed loop of magnetic dipole moments perpendicular to the meta-molecule surface. The latter dipolar configuration implies the excitation of a so-called dynamic toroidal dipole moment 69,72,73 directed parallel to the x-axis. On the other hand, the AM 2 mode possesses two pairs of such poloidal current loops and result in four field zeroes along the y-direction, indicating a clear combination of the AM 1 mode and an accompanied standing wave character. This phenomenon can be explained by the formation of hybrid Mie-Fabry-Perot modes 74 or the superposition of several internal modes 75 . Further, Fig.2e,f represents the spatial distributions of near-field amplitude |E| |E inc | of the oligomer at two characteristic anapole wavelengths, where the near field enhancements provided by the hot spots of AM 1 and AM 2 modes are 44.7 and 91.5, respectively.
To go further into the aforementioned mode interaction and to elucidate the constitution influence, in what follows, we investigate the effect of the nanoparticles size as well as the illumination conditions on overall spectra. We expect to observe that changing the value of these parameters influences both the direct excitation as well as the mutual coupling between resonant modes.
Dependence of double-Fano resonance mode on the oligomer's size. Here, we illustrate how the oligomer size affect evolution of the double-Fano spectral line shapes, a fact that critically determines the plasmonic structure design for SECARS applications. To this end, the size of oligomer under-investigation is scaled down by decreasing radius of the individual nanoparticles while keeping the radius to gap ratio fixed and equal to 7. The remaining optical and geometrical parameters are the same as Fig.1. From Fig.3a, it can be clearly observed that the Fano line profile with lower resonance wavelength ( AM 2 state) gradually disappears with decreasing the radius Figure 2. (a) The simulated SCS spectra for our proposed oligomer as well as contributions of the five dominant multipole moments to the total radiated field. "Sum Scat" states for the SCS as the sum of the multipole contributions;"Total Scat" states for the total SCS calculated directly in COMSOL (See "Computational methods"). The labels P, MD, EQ, MQ, and TD denote electric dipole, magnetic dipole, electric quadrupole, magnetic quadrupole, and toroidal dipole modes, respectively . Also, AM 1 and AM 2 represent the so-called first-and second-order anapole states with characteristic wavelength, respectively . (b) The phases of the electric dipole φ P and toroidal dipole φ TD moments as well as their difference �φ versus wavelength. The surface charge density distribution is shown at (c) AM 2 wavelength and (d) AM 1 wavelength. Distribution of electric field amplitude |E| |E inc | in the mid-height plane of the oligomer at (e) the second-and (f) first-order anapole modes. The SCS derived from ED results as the coherent superposition of P and TD contributions, symbolically represented as P+TD. www.nature.com/scientificreports/ R and the next one (representing AM 1 state) turns into a Lorentzian-shaped resonance band. Such transitions in the resonance line shapes can be quantitatively interpreted as follows: The multipole expansion analysis starts from the induced polarization current density J(r) in the nanoparticle, which given the convention e jωt for the harmonic electromagnetic fields, is calculated by 76 where ω is the frequency of the incident light in vacuum, E(r) is the induced electric field, and ε p = ε 0 ε r,p and ε d = ε 0 ε r,d are the dielectric permittivities of the nanoparticles and background medium, respectively, ε 0 being the vacuum permittivity. The relative permittivities ε r,p and ε r,d are equal to the square of the corresponding refractive indices n p , n d .
Let recall the recently obtained exact expression for the dipolar vector of electric field in the Cartesian coordinate basis 64 : where k d = ω/v d and v d = c/n d are respectively the wave vector and speed of light in the surrounding environment, and r = rr is the radius vector of a volume element inside the scattering medium (volume V). Making the small argument approximation to the spherical Bessel functions with terms up to fourth order (k d r) 4 (See Eq. (7) in "Methods"), and then grouping the contributions with the same power of (k d r) , we obtain where the contributions coming from l -th order spherical Bessel function of the first kind j¯l(·) are indicated. Here, the (k d r) 0 term, P α in Eq. (2), is the well-known small source approximation of the electric dipole moment of a current density distribution J α , and the (k d r) 2 term, T α in Eq. (2), represents the so-called toroidal dipole moment. In addition, (k d r) 4 term shown by R T α in Eq. (2) is referred to as the mean-square radii of the toroidal dipole moment 77 . Notably, the T α integral contains contributions from the two integrals in Eq. (1): It is the sum of the terms of order (k d r) 2 coming from both j 0 (k d r) in the first integral of Eq. (1) and j 2 (k d r) in the second one. By further decreasing the particle's size parameter k d r (e.g. having around 30 nm in radius disks), the highest www.nature.com/scientificreports/ order terms of j 0 (k d r) become negligible compared to the first one, leading to exclusion of term corresponding to l = 0 in T α formula. Accordingly, the toroidal dipole expression tends to a negligible value and causes the mutual coupling required for the excitation of a non-radiating anapole state, namely P = jk d T , to be destroyed. Following this discussion, Fig.3b,c clearly illustrate that the P multipolar contribution shows a transition to Lorentzian line profile (Fig.3b) and the TD becomes negligible (Fig.3c). This analysis not only gives insight regarding the design of SECARS substrate, but also serves as an additional verification making the anapole nature of the resultant scattering dips manifest.
Dependence of double-Fano resonance mode on the excitation angle. To further demonstrate the tailoring overall spectra characteristic, the angle of incident illumination is also varied from 0 to 90 and the spectra evolution is investigated. This study provides insights into how increasing excitation angle will affect the spectral response of our designed nanostructure. Considering the initial oligomer design (Fig.1), it is found that the observed anapole states progressively vanish from the total scattering spectra as the incident angle becomes more grazing (Fig.4a). The reason of this behavior can be clarified through identifying the multipolar contributions in terms of different incident angles, as shown in Fig.4b-e). It can be concluded that such evolution of scattering spectra is mainly a consequence of changes in the direct coupling between the incident plane wave and the different multipoles sustained by the structure. So, the amplitude of MD and EQ multipoles increases by the incident angle as a result of an increase of their interplay with the illuminated plane wave. On the other hand, P and TD multipoles amplitude is weakened and their spectral line profile noticeably changes with increasing the incident angle as well. These performances can be qualitatively interpreted as a consequence of a change in their direct coupling with incident plane wave combined with changes in their mutual interactions, respectively. It is pointed out that spectrum evolution of the MQ is close to that of TD given the existed correlation between these two multipoles along with an electric octopule moment 54 . In order to accomplish the goal of our study, one may take into account dependence of the dielectric function of gold nano-disks which can be affected by surface scattering of nano-scale metallic particles. The surface scattering can modify dispersion of gold nanoparticles which are much smaller than the incoming wavelength. As the thickness of the metal increases, discrepancy between the dielectric function of bulk and nano-scale materials tends to be negligible 78,79 . For all of the results discussed above, we disregarded modification of gold nano-disks dielectric function affected by surface scattering mechanism. To carry out validation of our approach, it should be noticed that our proposed oligomer is illuminated through an s-polarized plane wave whose electric field is parallel to the surface of nano-disks having notably large diameters (168 nm in diameter). Therefore, excitation of free electrons will be dominant along the transverse direction and, as a result, height of the nano-disks will not considerably affect scattering performance of our plasmonic structure. To clearly illustrate, in Fig.5 the scattering cross section for different heights of nano-disks are calculated in both normal and oblique incidence (θ = 60 • ) . It can be observed that by increasing the thickness of gold nano-disks, apart from a slightly higher scattered signal of the thicker disks together with minor spectral shifts of the modal resonances potentially used in SECARS, no significant changes are experienced in the total scattering spectra. www.nature.com/scientificreports/

Spectral manipulation of anapole states for SECARS applications.
In this section, we demonstrate that the optical response of our proposed plasmonic oligomer can be manipulated such that the incident beams at the pump ω p and Stokes ω s frequencies interact coherently with medium to generate the emission of the CARS photons at the blue-shifted anti-Stokes frequency ω As , i.e., ω As = 2ω s − ω p 51 . Note that in the CARS process, one of a four-wave mixing processes involving three laser fields interacting with the sample, not only the frequencies of the incident and scattered light, but also the wavevectors of each beam including phases should be considered. This fact demonstrates to be critical in the first experimental implementation into microscopy using noncollinear excitation of pump and Stokes visible dye laser beams 80 . In the plane-wave approximation for a nonabsorbing material, the CARS signal intensity depends linearly on the term sinc 2 � � k · � l/2 , where l is the thickness of sample, k the phase mismatch � k = � k p − 2 � k s + � k As and � k i i = p, s, As are the wavevectors for each one of the interacting (pump, Stokes and anti-Stokes) lights. As a consequence, generation of the CARS signal needs the phase-matching condition � � k · � l ≪ π to be fulfilled, which reveals the coherent nature of the CARS process. Furthermore, for a very thin film the phase-matching condition is usually satisfied both in forward and backward directions relative to the propagation direction of the excitation beams 81 . When using SECARS substrates, the effective thickness where light interaction is enhanced typically lies under a subwavelength scale, as it occurs in our proposal, thus enabling the phase-matching condition to be fully satisfied. Nevertheless, such condition is further constrained to excitation under normal incidence when using a coplanar array of metallic octamer structures 82 .
The observed non-radiating anapole states are interesting for CARS applications due to the several reasons: (1)Given the lack of scattering and radiation loss in a dipole channel, a highly concentrated electromagnetic near-field in the center of the structure is generated boosting nonlinear effects, (2)The higher-order anapole states provide stronger energy localization and narrower resonance windows, an exotic feature that is advantageous for nonlinear higher-harmonic generators, (3)The anapole wavelengths can be easily measured through identifying the far-field scattered spectrum of white light from the structure. As a consequence, our observed anapole states are nearly ideally suited for the SECARS process.
Revisiting the surface charge density distributions corresponding to AM 1 and AM 2 anapole modes shows that the corner particles of the oligomer under study only contribute to the formation of external current loops (not the central ones) in case of AM 2 mode whereas these particles are essential elements to produce current loops corresponding to AM 1 mode. Since electric field hot spots of two anapole states are mostly concentrated in the gap region of the central dimmer as seen in Fig.2d-e), it can be concluded that the central current loops play the main role in characteristics of the observed Fano-like line profiles. It has been shown previously that the depth of Fano-like spectral windows can be modulated by varying the LSPR coupling strength between adjacent nanoparticles. Furthermore, the spectral position of scattering dips can be tuned through controlling the radius of associated current loops 83,84 . Following this idea, it is expected that modulating the gap distance of the corner particles here will change the coupling strength between the nanoparticles as well as the size of main current loops in case of AM 1 mode and the size of external current loops in case of AM 2 anapole mode. Therefore, adjusting the gap separation between corner particles will result variations of both resonance depth and spectral position for the AM 1 mode. Concerning the AM 2 mode, on the other hand, such an adjustment will only affect depth of the scattering dip without displacing its spectral position. The reason of this performance can be inferred from the size of the main (internal) current loops of the AM 2 mode remaining unchanged. To examine the above qualitative analysis, we vary the gap distances of the corner particles on the top and bottom sides by considering different values of deviations d from the initial value g = 12 nm . We numerically calculate the total scattering spectra of the resultant system. The simulated responses are shown in Fig.6 where S = 2R + g + d , demonstrating the agreement with the aforementioned analysis of the anapole modes. www.nature.com/scientificreports/ Interestingly, the obtained results exhibit independent spectral tunability of the first-order anapole state ( AM 1 mode) with respect to the variations of gap distance of the corner particles on the top and bottom sides. This feature is the design key to simply adjust spectral response of the first-order anapole mode independent from the second one for SECARS applications.
It should be most advantageous for SECARS enhancements to tune the two anapole states at the pump and Stokes wavelengths while the superradiant shoulder (at the blue side of the both anapole states) at the CARS wavelength. This strategy results minimizing the radiation losses at both excitation wavelengths while maximizing the far-field coupling at the wavelength of anti-Stokes emission 85 . The reduction of the scattered light at the wavelengths associated with the anapole states ensures efficient energy coupling into the nanostructure, thus generating highly localized, intense field enhancements. In addition, the highly-radiative shoulder provides an enhanced propagation of CARS signal to the far-field region. It can be seen the SECARS condition is fulfilled for the deviation value of d = 40 nm as shown in Fig.7. The corresponding wavelengths are As = 653 nm , p = 720 nm , s = 803 nm , where the Stokes frequency and the pump laser excitation frequency matches to the first-and second-order anapole states respectively and the high-frequency shoulder overlaps with the anti-Stokes mode. These three involved frequencies in scattering spectrum obviously satisfy ω As = 2ω s − ω p .
Efficiency of the SECARS substrate is identified by maximizing the total electromagnetic field EF at the anti-Stokes frequency through the relation 86 : where g p , g s , and g As are the localized electric field enhancements corresponding to the pump, Stokes, and anti-Stokes frequencies, respectively. Thus, an appropriate plasmonic substrate staisfying Raman dipole-forbidden band distribution condition of the molecule will present strong localized hot spots depending on three characteristic frequencies. Spatial overlap of the hot spots at three resonant modes that form the "mixed frequency coherent mode" is important for SECARS application and it will be useful to reach a strong EF. The spatial electric field distributions at three involved wavelengths of our designed substrate are shown in the Fig.7b-d. The fields are monitored at the middle height of the oligomer (z = 10 nm) and they confirm high simultaneous electric field enhancement with the hot spots at the same spatial locations (gap of the central dimer), both necessary conditions to yield an efficient SECARS substrate. The corresponding electric field EFs are g s = 73.3 , g p = 40.6 , and g As = 26.7 leading to the total EF of G SECARS = 1.05 × 10 13 .
Through adjusting geometrical parameters of our proposed oligomer, the wavelengths of anti-Stokes, pump, and Stokes modes can be freely tuned, providing a set of different wavelengths to be used in SECARS applications. In Table 1, we have summarized the value of SECARS EF for different radius of nanoparticles where the radius to gap ratio is considered to be fixed (R/g = 7) . For the design parameters of R = 77 nm and g = 11 nm , anti-Stokes, pump and Stokes wavelengths are obtained as As = 634 nm , p = 695 nm , s = 769 nm , respectively. Here, the pump frequency falls into the visible region which is of great interest in most SECARS applications. Using this design, one can utilize a short wavelength laser to pump the substrate and materials. In addition, it should be noted that gap spacing between nanoparticles plays a crucial role in the resonance modes and hot spot field enhancement. Shrinking the gap distance to g = 8 nm while the remaining parameters being the same as those in Fig. 2 yields a As Figure 6. Simulated total SCS of our proposed oligomer considering various gap distances for corner particles denoted by the deviation S = a + d . The other geometrical parameters are chosen as the same as those in Fig.1.
The results show the independent spectral tunability of the first-order anapole mode ( AM 1 mode) from the second one. It can be seen that the SECARS condition is fulfilled for the deviation value of d = 40 nm. www.nature.com/scientificreports/ S3 in Supplementary Information) which is almost two orders of magnitude larger than those having the same gap size reported in the literature 52 , along with a red-shift of the characteristics wavelengths. Let us conclude our discussion by highlighting that our proposal might potentially optimize its performance through a fullwave Maxwell inverse design 87 . For instance, tapered ultranarrow gaps enable field enhancements exceeding 10 2 in Ref. 88 , suggesting that surface engineering of the central dimer in our nanostructure thus shrinking the interparticle gap domain while maintaining the gap distance fixed might be favorable for the generation of higher CARS signals. In addition, by optimizing the oligomer's geometrical parameters, we can tune the plasmonic substrate geometry so that its multiple resonances matches to the frequencies of the pump, Stokes, and anti-Stokes lights in the visible as well as near infrared spectral range.

Conclusion
In summary, we have demonstrated an engineered assembly of metallic disk-shaped nanoparticles to be a highlyefficient free-standing substrate for SECARS application. To this end, we first propose a simple plasmonic oligomer structure supporting double-Fano resonance line shapes, which serves as a basis for our proposed SECARS substrate. Expanding the total field scattered by the whole structure into Cartesian multipolar basis reveals that the two observed characteristic dips corresponds to the so-called non-trivial first-and second-order anapole states. Distribution of electric current loops associated with these anapole modes suggests independent spectral tunability of the first-order anapole state through modulating gap distance of the corner particles. This feature along with giant electric field enhancements provided by anapole states offers a simple design procedure in the  Figure 7. (a) Simulated total SCS spectra of our designed SECARS substrate with deviation value of d = 40 nm where the vertical dashed green, blue, and red lines respectively denote the As = 653 nm , p = 720 nm , and s = 803 nm satisfying the SECARS condition as ω As = 2ω p − ω s . Spatial distributions of the electric field intensity corresponding to the three involved (b) anti-Stokes, (c) pump, and (d) Stokes signals, when monitored at the middle height of the oligomer (z = 10 nm). www.nature.com/scientificreports/ context of SECARS applications. Notably, the electric field "hot spots" corresponding to three involved wavelengths, i.e., anti-Stokes, pump, and Stokes, are localized at the same spatial position, which can be exploited to achieve a huge CARS signal field enhancement. The estimated total enhancement factor of SECARS process for our proposed design is up to 10 15 providing the highest enhancement ever-obtained for the inter-disk spacing of g = 8 nm.

Methods
Multipole expansion method. In standard electrodynamic textbooks 89,90 , the spherical multipoles result from the multipolar decomposition of electromagnetic fields in terms of the spherical vector harmonics. However, the complexity in interpreting the spherical basis is caused the spherical multipoles not to be analyzed directly. To improve this situation, recently exact expressions for the localized charge current density multipoles in the Cartesian basis were developed 64,91 . It is noted that the obtained expression set implicitly comprises toroidal multipole moments.
Exact Cartesian multipole decomposition. In the Cartesian coordinate system with the origin coinciding with the mass center of the nanoparticle, the ordinary ED moment of the nanoparticle is expressed by Eq. (1) in the main text, and the other leading multipoles up to MQ are described as where d M α , Q E αβ , Q M αβ are respectively the magnetic dipole (MD), electric quadrupole (EQ), and magnetic quadrupole (MQ) tensors, δ is the Dirac delta function, and subscripts α, β = x, y, z . Commonly, a long-wavelength approximation is utilized making the above expressions sufficiently straightforward. The so-called toroidal dipole moments are derived explicitly by this kind of formulation. , and Q T αβ denote respectively magnetic dipole, electric quadrupole, magnetic quadrupole, and toroidal quadrupole moments of Cartesian multipole decomposition. In addition, R 2 m and R 2,Q m αβ represent the mean-square radii corrections for the magnetic dipole and magnetic quadrupole moments, respectively 77 .
The partial scattering cross-sections corresponding to each multipole moment are given by assuming |E in | as the electric field amplitude of the incident plane wave. It is convenient to define also the magnitudes C P sca and C T sca as which correspond to effective cross-sections of the Cartesian electric and toroidal dipoles, respectively. These two dipoles may interfere constructively or destructively, depending on their relative phase difference, and their combined contribution to the electric dipole scattering cross-section is described in Eq. (11). Considering the other higher order terms are neglected in the investigated nanoantenna, the total scattering cross-section C total sca from the above multipoles can be approximated as It is well worth mentioning that in spite of popularity of the approximate Cartesian multipoles (known as Cartesian multipole decomposition elements), they cannot exactly reconstruct the electrodynamic radiation fields and the related scattering phenomena in case of large particles as well as high refractive-index particles.
Computational methods. Numerical simulations of the structure discussed in this study were carried out using the finite element method (FEM) with the implement of COMSOL Multiphysics, where Perfectly Matched layers (PML) were utilized to avoid spurious reflections at the surrounding boundaries of plasmonic nano- www.nature.com/scientificreports/ systems. The induced electric field E(r) inside the structure was calculated in the framework of the scatteredfield formulation. To directly determine the total scattering cross section, the scattered flux in a closed surface surrounding the oligomer was calculated using the Poynting theorem. The surface charge density distributions are obtained by considering the skin effect and applying Gauss' law during FEM calculations 92 . In all cases, a linearly-polarized light with the electric field polarization along the gap axis of the central dimer (x-direction) was used as the excitation source, where the propagation direction coincides with the z-axis. Here, k in , E in , and H in are the wave vector, electric field, and magnetic field of the incident light, respectively. To reach sufficiently high simulation accuracy the particles need to be finely meshed. Here the maximum mesh size is 8 nm. The parameter sweep module in this study is utilized for the spectrum simulation covering a range from 550 nm to 1200 nm.