## 2 Simulation Method and Models  
A typical approach to study the optical response of the plasmonic molecule is a semianalytical multiple multipole program [4; 12; 13] based on boundary discretization. This method was good for understanding the variation of the electromagnetic field with gradient structure modification, and was valid when \(ka<<1\). Spectra of nanoparticle clusters are analyzed based on their structural symmetry using symmetry-adapted linear combinations of the individual nanoparticle plasmon modes [4]. In this paper, we get the full vertical three-dimensional electromagnetic field (E, H) of the nanostructures by calculating the Maxwell's equations using the finite element method software COMSOL Multiphysics (V5.2, Radio Frequency module, Stockholm, Sweden). Scattering spectra and energy diagrams of the graphene-based trimeric plasmonic molecules are obtained by employing the frequency domain and the eigenfrequency solvers of the commercial software. Then the evolution of the collective mode of the trimer is investigated with the help of irreducible representations of the relevant point group, which is similar to other studies [4; 12]. More significantly, the tunable properties of the plasmonic molecule are demonstrated by modifying the chemical potential of the graphene without changing the geometric parameters.  
The geometry of the proposed structure is schematically shown in Figure 1a, where three graphene nanodisks are arranged in a triangle and lie on the SiO\({}_{2}\) substrate. The nanocluster is illuminated by transverse magnetic (TM) polarization light with electric field E along the longitudinal direction. The extinction cross-section \(\sigma_{ext}\), composed of scattering cross-section \(\sigma_{sc}\) and absorption-cross section \(\sigma_{abs}\), is expressed as  
\[\sigma_{ext}=\sigma_{sc}+\sigma_{abs}, \tag{1}\]with  
\[\sigma_{sc}=\frac{1}{I_{0}}\iint(\overrightarrow{n}\cdot\vec{S_{sc}})dS, \tag{2}\]  
where \(\overrightarrow{n}\) indicates the normal vector pointing outwards from the plasmonic molecule, and \(\vec{S_{sc}}\) and \(I_{0}\) are the scattered intensity (Poynting vector) and the incident electromagnetic energy intensity, respectively. The integral is taken over the closed surface of trimer,  
\[\sigma_{abs}=\frac{1}{I_{0}}\iint QdV, \tag{3}\]  
where \(Q\) is the power loss density in the graphene molecule and the integral is taken over the volume.  
The graphene nanodisk is treated as a thin film with an equivalent thickness of \(\Delta\) in our study. Thus, it is characterized by an equivalent permittivity of \(\varepsilon_{g,\mathrm{eq}}\) (\(\Delta\)) [40; 41], which is a function of both \(\Delta\) and the surface conductivity of graphene given by \(\varepsilon_{g,\mathrm{eq}}=-i\sigma_{g}\eta_{0}/(k_{0}\Delta)\). Here, \(\eta_{0}\) (\(\approx\)377 \(\Omega\)) is the impedance of air, \(k_{0}=2\pi/\lambda\) represents the wave number of the free space, and \(\sigma_{g}\) is the contributions from the interband electron-electron transition \(\sigma_{inter}\) and the intraband electron-photon scattering \(\sigma_{intra}\). According to the Kubo formula [34; 42; 43], graphene's surface conductivity is written as \(\sigma_{total}=\sigma_{inter}+\sigma_{intra}\),  
\[\sigma_{inter}=i\frac{e^{2}}{4\pi\hbar}\ln[\frac{2|\mu_{c}|-\hbar(\omega-i \tau^{-1})}{2|\mu_{c}|+\hbar(\omega-i\tau^{-1})}], \tag{4}\]  
\[\sigma_{intra}=i\frac{e^{2}k_{B}T}{\pi\hbar^{2}(\omega+i\tau^{-1})}[\frac{\mu _{c}}{k_{B}T}+2\ln(\exp(-\frac{\mu_{c}}{k_{B}T})+1)], \tag{5}\]  
where \(\mu_{c}\) is chemical potential, which is tunable by either chemical doping or external gate voltage. \(\tau\) is the momentum relaxation time of electrons, and it is set to 0.5 ps throughout this paper. Twenty cores of Xeon CPU and 512 Gb memory were employed to run the calculations. It took 47 h to complete each run.