Tailoring the lineshapes of coupled plasmonic systems based on a theory derived from first principles

Coupled photonic systems exhibit intriguing optical responses attracting intensive attention, but available theoretical tools either cannot reveal the underlying physics or are empirical in nature. Here, we derive a rigorous theoretical framework from first principles (i.e., Maxwell’s equations), with all parameters directly computable via wave function integrations, to study coupled photonic systems containing multiple resonators. Benchmark calculations against Mie theory reveal the physical meanings of the parameters defined in our theory and their mutual relations. After testing our theory numerically and experimentally on a realistic plasmonic system, we show how to utilize it to freely tailor the lineshape of a coupled system, involving two plasmonic resonators exhibiting arbitrary radiative losses, particularly how to create a completely “dark” mode with vanishing radiative loss (e.g., a bound state in continuum). All theoretical predictions are quantitatively verified by our experiments at near-infrared frequencies. Our results not only help understand the profound physics in such coupled photonic systems, but also offer a powerful tool for fast designing functional devices to meet diversified application requests.


Introduction
Recently, photonic systems consisting of multiple plasmonic/dielectric resonators coupled in different ways have attracted much attention [1][2][3][4] . Compared to simple systems containing only one type of resonators, coupled systems exhibit more fascinating near-field (NF) properties (e.g., local field enhancement) and far-field (FF) responses manifested by unusual lineshapes, such as Fano resonance [5][6][7] and Rabi oscillations 8,9 , dictated ultimately by how the involved resonators are coupled together. Couplings have offered more opportunities for controlling the NF and FF light environments of such complex photonic systems as desired, making them particularly useful in applications, such as nanolasing 10,11 , fluorescence enhancement [12][13][14] , and information transport [15][16][17] .
Despite great advances on the experimental side, the theoretical understandings of such systems are far from satisfactory, which also hinders the rapid design of appropriate systems with desired NF and FF responses. For example, full-wave simulations require huge computing costs and reveal very little physics. Meanwhile, although many models (e.g., coupled-mode theory (CMT) 18-21 , Fano's formula 22,23 , or effective circuit models 24,25 ) were proposed to analyze the underlying physics, they typically require model parameters fitted from simulation results, and thus cannot predict unknown phenomena before having studied the systems numerically. As an early attempt, a photonic tight-binding method (TBM) 26 , with all involved parameters computable without fitting procedures, was proposed to successfully predict the resonance peak positions of a coupled system. Unfortunately, the TBM provides no information on the entire optical responses (e.g., the lineshapes), which are usually more desired for practical applications. The intrinsic difficulties are that these systems are open in nature, in which different resonators can couple not only with each other via NFs but also, more importantly, with external free space via FF interactions (Fig. 1a). To establish a complete theory to predict the entire optical properties of arbitrarily coupled photonic systems, one needs to rigorously consider both NF and FF interactions on the same foot. While several semi-analytical approaches have recently appeared, they have their own limitations and are not generic enough to study arbitrarily coupled systems in a formal way [27][28][29] .
In this paper, we derive a formal theoretical framework from first principles (i.e., Maxwell's equations), with all involved parameters directly computable without fitting procedures, to predict the optical lineshapes of arbitrarily coupled photonic systems. The obtained equations resemble the empirical CMT but are derived from first principles, and thus have unambiguous physical meanings, as clearly revealed by benchmark calculations against rigorous Mie theory on a model system. After validating our theory through comparison with experimental/numerical results on a realistic plasmonic metasurface, we present how to employ it to tailor the lineshape of a coupled plasmonic system as desired by varying the interresonator coupling. In particular, we show that it is possible to generate a completely "dark" optical mode with vanishing radiative loss (i.e., a bound state in continuum (BIC) 30,31 ) in such systems, although the constituent resonators exhibit moderate radiative losses. All theoretical predictions are quantitatively verified by experimental results on a series of metasurfaces containing plasmonic resonators coupled in different ways.

Establishment of the formal theory
We start by establishing a formal theory applicable to generic coupled open systems. As shown in Fig. 1a, we consider the scatterings of a system consisting of M arbitrary resonators located at different positions in a host medium under certain external illumination. Such an open system can be schematically described by the model depicted in Fig. 1b, where the region containing resonators is connected to the external continuum via N ports with well-defined properties. Formally, we need to solve the following Schrödinger-like equation: where Ψðr; ωÞ is the total wave function, and b V m is the Hamiltonian of the whole system with b H h describing the host medium and b V m the potential contributed by the mth resonator.
To expand the unknown function Ψðr; ωÞ appropriately, we need a complete set of basis wave functions that are orthogonal to each other and normalizable in certain ways. In the same spirit as the TBM 26,32 , here, we define a set of wave functions fψ LEM m ðr; ωÞ; m ¼ 1; :::; Mg, which are the (approximate) solutions of the Hamiltonian b , describing the subsystem containing only the mth resonator. For simplicity, here, we assume that each resonator supports only one mode, and the extensions to more general cases (e.g., resonators exhibiting multiple or degenerate modes) are straightforward. Different from the systems treated by the TBM, which are closed 26 , and thus have well-defined localized eigenfunctions, here, the open systems under study only support leaky eigenmodes (LEM), as explained subsequently.
Suppose that the resonators exhibit high quality (Q) factors; we can use the following approach to obtain ψ LEM m ðr; ωÞ. Shining the subsystem with external illumination, we can solve b H m Ψ m ¼ ωΨ m to obtain Ψ m analytically or numerically, and then obtain the response spectrum of the system. We then identify the resonance frequency ω m of the mth resonator from the maximum of  Fig. 1 Schematics of the system under study and our theory. a Photonic system containing multiple arbitrary resonators coupled together under external illumination. The inset shows a typical optical lineshape of such a system. b Schematics of our theory: under certain external illumination, the total scattered field of the coupled system is a linear combination of leaky eigenmodes (LEM, ψ LEM m ) of different resonators, each containing a near-field part ψ NF m and a far-field tail ψ FF m the response spectrum. Choosing a "background" representing the system at a frequency far from all resonances, we can calculate the background wave function Ψ B by shining the "background" medium with the same external illumination. We finally obtain the desired LEM wave function through ψ LEM m ¼ Ψ m À Ψ B for the mth resonator. We note that fψ LEM m g are quite different from the quasinormal-mode (QNM) functions defined in refs. 27,28,33 . While fψ LEM m gare wave functions of the systems under external illumination at real frequencies ω m , QNM functions are eigenfunctions of the systems without external illumination corresponding to complex eigenfrequencies. Moreover, LEM functions do not diverge at infinity, whereas QNM functions inevitably diverge 34,35 . Therefore, LEM functions are particularly suitable for the lineshape problems studied here, which require external illumination. Examples of how to obtain ψ LEM where ψ NF m and ψ FF m represent the NF and FF parts of the wave function, respectively. Technically, for any given system with well-defined external ports, we can always project ψ LEM m onto the port modes on reference planes of all external ports and then construct ψ FF m by these port modes, which are assumed to fill the entire space. With ψ FF m known, we then obtain ψ NF m numerically based on Eq. (2). The NF functions ψ NF m have good properties to help us perform further analyses. In the vicinity of the scatterer, under the high-Q approximation where the FF part of the wave function is significantly weaker than the NF part, ψ NF m can be approximately viewed as the eigenfunction of Meanwhile, ψ NF m can be normalized since it is well localized around the mth resonator. Moreover, considering that these wave functions are spatially well separated, we find that they approximately satisfy the following orthonormal condition: where the integrals are performed over the entire space. We note that one needs to multiply ψ LEM m by the same normalization constant that is used to normalize ψ NF m , since these two functions are connected by Eq. (2). Equation (4) indicates that fψ NF m ; m ¼ 1; :::; Mg form a set of orthogonal bases to expand the total wave functions in the NF region. Note that the approximation Eq. (4) is widely used in the TBM for treating electrons in solids 26 .
We now identify the FF eigenbases of the system. In the FF region, eigenmodes are just a set of propagating modes k ± q E n o allowed by the system, where +(−) denotes the incoming (outgoing) propagation direction, q labels the mode channel, and k is the wavevector satisfying certain dispersion relation kq (ω). These wave functions satisfy the following orthogonal condition: where the integrals are performed on the reference plane of a particular external port. In principle, extending our theory to study cases with continuum scattering ports 36,37 is also possible, although one needs to compute all parameters related to these scattering channels.
We are now ready to represent Ψ as a linear combination of these basis functions. We have Ψ = Ψ B + Ψ sca , where Ψ sca is contributed by the scatterings of all resonators. In the same spirit as the TBM, Ψ sca can be approximately written as a sum of scattered fields Ψ sca m associated with each individual scatterer. At first glance, one may expect that Ψ sca m ðr; ωÞ must be ψ LEM m ðr; ω m Þ defined previously. However, ψ LEM m ðr; ω m Þ is the scattered wave at resonance frequency ω m , not at arbitrary frequencies as required in Eq. (1). We can amend ψ LEM m ðr; ω m Þ slightly to obtain the form of ψ LEM m ðr; ωÞ for a frequency ω not far from ω m . The NF part [ψ NF m ðr; ω m Þ] is solely determined by ω m , as it is (approximately) an eigensolution of Eq. (3) for eigenfrequency ω m . Since, we will need to utilize the orthonormal properties of ψ NF m ðrÞ offered by Eq. (3) later, here, we take the original form of ψ NF m ðrÞ in constructing the trial wave functions at general frequencies ω ≠ ω m . Meanwhile, the FF part ψ FF m ðr; ω m Þ contains propagating terms depending on the wavevector k q , which must be modified from k q (ω m ) to k q (ω) according to the dispersion relations. We note, however, that ψ LEM m ðr; ωÞ thus obtained neglects the frequency corrections to the FF radiation amplitudes. In principle, such corrections can be taken into account by considering the NF-FF relation of a given source 38,39 . To obtain a concise analytical form for our theory, here, we neglected such corrections, justified by the high-Q approximation. Later, we show that such an approximation works quite well even though the original modes supported by individual resonators do not exhibit extremely high-Q factors.
We can finally construct the total wave function as where {a n } are a set of unknown coefficients representing the strengths of fields scattered by different resonators under external illumination represented by fs þ q g denoting the excitation amplitudes at different incoming ports, and Ψ q B denotes the background wave function obtained when only the qth port is excited with unit amplitude. Substituting Eq. (6) into Eq. (1), projecting both sides by ψ NF m and utilizing the orthogonal condition Eq. (4), we obtain the following equations to determine{a n }: We next multiply both sides of Ψðr; ωÞ defined in Eq. (6) by each FF outgoing basis k À q D , and then perform the field integrations at the reference planes of all ports. Using the orthonormal conditions Eq. (5), we finally obtain the set of equations: , which describe the strengths of scattered fields measured at different external ports. Here, all parameters in Eqs. (7) and (8) are unambiguously defined and can be calculated via the following integrals: where "V" and "S" denote whether the integrals are performed over the entire volume or at the reference plane of a port. The physical meanings of all involved parameters can be clearly seen from their expressions. For example, t mn and X mn represent the coupling strengths between two resonators due to their NF and FF interactions, respectively. Derivations of Eqs. (7)-(9) can be found in Sec. II of the Supplementary Information. It is helpful to explicitly discuss the conditions imposed on our systems to make the derived theory (e.g., Eqs. (7)-(9)) valid. By re-examining Eq. (7) for the single-scatterer case, we find that Im(Γ m ), if it exists, can shift the resonance frequency ω m , and thus, a large Im(Γ m ) implies that ψ NF m is not reasonably chosen. Therefore, the first criterion is ImðΓ m Þ ! 0, which determines the accuracy of our theory at resonance. Meanwhile, we also require ReðΓ m Þ << ω m , which is responsible for the correctness of our theory in describing the entire lineshape. The second criterion can be easily satisfied by a moderate Q value (e.g., Q > 5), as long as the frequency dispersion of the material is not significant and high-order modes are all far from the mode under study. The first criterion, however, requires the resonators to be deep subwavelength in size so that ψ FF m and ψ NF m can exhibit a π/2 phase difference inside the whole region occupied by the resonator 38 , leading to a negligible Im(Γ m ). For plasmonic resonances, such a deep-subwavelength condition is easily satisfied. However, for dielectric resonances, such a condition can only be satisfied in systems with a very high refraction index (n), which pushes the Q factors to even higher values (see Sec. IX in the Supplementary Information for more details).
We note that Eq. (9) is derived for lossless systems, and thus, Γ m must only contain radiation damping. In realistic systems, we also need to consider another parameter Γ a m , representing the damping due to absorption (i.e., replacing Γ m by Γ m þ Γ a m in Eq. (8)). This parameter can be computed using Γ a m represents the Hamiltonian of the realistic lossy systems, while b H 0 m describes the same system with material losses omitted 40 . Equations (7)-(9) are the core results of this paper, which have clear and profound physical meanings. While Eq. (7) describes the dynamics of each mode under certain excitations, Eq. (8) describes the measurable scattering spectra. We note that Eqs. (7) and (8) resemble the two equations in CMT 18,19 , but our theory is different and possesses the following merits. In the empirical CMT, the key parameters defined are usually obtained by fitting with numerical simulations, while the remaining parameters can be derived by energy-conservation and time-reversal arguments 41 . In contrast, here, in our theory, all parameters can be unambiguously evaluated by Eq. (9), and therefore, one can use it to predict the lineshapes of coupled systems before performing numerical simulations on them. Moreover, the empirical CMT cannot explicitly consider the NF couplings between resonators 18 , while in our approach, NF couplings t mn can be unambiguously determined (see Eq. (9)) and explicitly included in determining the lineshape (Eq. (8)).
Although single-resonator parameters (ω res and Γ res ) can be analytically obtained for certain high-symmetry structures for that analytical formulas of scattering coefficients are available 42 , such an approach is not general enough to deal with arbitrary coupled systems without analytical expressions of scattering coefficients and cannot be used to study the couplings between different resonators.

Applications to photonic systems and benchmark tests
We now apply the developed formal theory to photonic systems, described generally by an inhomogeneous permittivity function εðr; ωÞ, in which at each local pointr, the permittivity is εðωÞ ¼ ε 1 ½1 þ ω 2 p =ðω 2 0 À ω 2 þ iωΓ e Þ, where ε ∞ , ω 0 , ω p , and Γ e are all position-and frequencyindependent parameters, describing the local properties of constituent materials. The governing equations (i.e., Maxwell's equations in the frequency domain) can be formally rewritten as Eq. (1) 40 , where the Hamiltonian is given by and the wave function is defined as ΨðrÞ ¼HẼP À V Þ T , withẼ,H, andP denoting the electric, magnetic, and polarization fields, respectively, andṼ ¼ dP=dt describing the polarization current. Consider the lossless case first (i.e., Γ e = 0). The inner product between two wave functions is defined as 26,40 Meanwhile, in the FF region occupied by air, the inner product between two-port modes can be defined as wherec is the light speed in the host medium. This ensures that different port modes are orthogonal and that each mode carries a unit of energy flux 43 . With Eqs. (10)- (12) and supposing that fψ NF m ; ψ FF m g are obtained, one can substitute them into Eq. (9) to compute all parameters (see Sec. III in the Supplementary Information) and then substitute them into Eqs. (7) and (8) to determine the lineshape.
For photonic resonators with regular shapes, fψ NF m ; ψ FF m g can be obtained analytically. For arbitrary resonators, we need to numerically obtain the required wave functions. We emphasize that, however, such numerical calculations are only needed once. Once fψ NF m ; ψ FF m g are obtained, we can predict the lineshapes of the coupled systems without having to perform simulations on them.
We first choose an analytically solvable system-a single gold sphere illuminated by an x-polarized plane wave-to test our theory against Mie theory. As shown in Fig. 2a, consider a sphere located at the origin with radius r m = 0.036λ p and Drude permittivity εðωÞ ¼ ε 0 ½1 À ω 2 p =ω 2 , with ω p and λ p denoting the plasmon resonance frequency and wavelength. Such a problem can be analytically solved by Mie theory 44,45 , yielding an analytical form of Ψ sca ðr; ωÞ. When the scatterer is much smaller than the wavelength of incident light, the electric dipole channel dominates in the frequency range plotted 46 , and thus, we can obtain ω res ¼ ½1 À 8π 2 ðr m =λ p Þ 2 =15 þ :::ω p = ffiffi ffi 3 p and the analytical forms of ψ LEM ; ψ FF , and ψ NF , as well as k ± j i (see Sec. IV in the Supplementary Information). Figure 2a depicts the field distributions of ψ LEM ; ψ FF , and ψ NF : ψ NF exhibits a clear electric dipole resonance feature, and ψ FF represents the FF radiation of an electric dipole located at the origin.
Substituting all wave functions into Eq. (9), we find κ ¼ d ¼ 2:92 10 À2 ffiffiffiffiffi ffi ω p p i and Γ ¼ 4:28 10 À4 ω p . Since there is only one scatterer and one port in the system, we neglect all subscripts without causing confusion. Substituting these parameters into Eqs. (7) and (8), we obtain the scattering spectrum of the nanosphere, defined as σðωÞ ¼ 3π 1 À R j j 2 =ð2η 0 k 2 0 Þ, with η 0 ¼ ffiffiffiffiffiffiffiffiffiffiffi μ 0 =ε 0 p being the vacuum impedance and R s À =s þ , representing the scattering coefficient. The spectrum thus calculated is depicted in Fig. 2b as a solid line, well matching the Mie theory (squares) and FEM calculation (circles) results.
Under the electric dipole approximation, we further simplify the analytical expressions of all involved parameters (see Sec. IV in the Supplementary Information) as withp ¼ R sphereP ðrÞdr (P is the polarization field inside the sphere; see the inset in Fig. 2a), representing the effective dipole moment of the nanosphere. Equation (13) reveals a few important physics difficult to obtain from numerical calculations. First, κ and d, defined as two distinct field integrations (Eq. (9)), surprisingly generate identical results (see Eq. (13)), which is consistent with the time-reversal symmetry argument 19 . Second, Γ takes an expression identical to that derived for a dipole emitter based on Poynting's theorem (see Eq. (8.74) in ref. 38 ), revealing the clear physical meaning of the radiation damping. Finally, Eq. (13) uncovers the relation 2Γ ¼ p j j 2 ω 4 res =ð6πε 0 c 3 Þ ¼ d j j 2 verified by numerical calculations (see Fig. 2c), which ensures energy conservation consistent with Poynting's theorem 38 . We note that these relations were derived by energy-conservation and timereversal arguments in the empirical CMT. Here, they are directly and rigorously demonstrated in our theory simply because our theory is established based on Maxwell's equations, which already satisfy energy-conservation and time-reversal symmetry.
After studying coupled electric dipole resonators to justify our theory against analytical formulas derived in prior literature 47,48 (see Sec. V in the Supplementary Information for details), we implement our theory to study arbitrary photonic coupled systems. As shown in Fig. 3a, the system we consider is a periodic metasurface with unit cells arranged in a hexagonal lattice (with periodicity 550 nm), each containing two different types of nanoparticles (bar and C-shaped resonator) coupled together. All nanoparticles are made of silver and are placed on a semi-infinite dielectric substrate (n = 1.55). Following the general strategy established above, we first perform lossless FEM simulations to study the scattering properties of two model systems, each containing resonators of a particular type arranged in the same hexagonal lattice (see Fig. 3a). Due to the periodic arrangements with deep-subwavelength spacing, only the zero-order transmission/reflection channels survive in the FF. From the calculated reflection spectra (circles) shown in Fig. 3b, c, we identify the resonance frequencies {ω m , m = 1, 2} of the two resonators (see dashed lines in Fig. 3b, c). We then follow the general strategy described in the last section to determine the needed NF and FF wave functions fψ FF m ; ψ NF m ; m ¼ 1; 2g. Substituting these singleresonator properties into Eq. (9), we obtain all needed parameters (see Sec. VI in the Supplementary Information for details) and, in turn, the desired transmission/reflection spectra. The reflectance spectra calculated by our theory are plotted in Fig. 3b, c as black lines, in perfect agreement with FEM simulations (circles) of realistic structures. This is remarkable since we did not perform any fitting procedures in obtaining these spectra. The lineshape of the coupled system predicted by our theory is further confirmed by our experiments. We fabricated three samples according to the designs using the standard electron-beam lithography (EBL) method (see left panel in Fig. 3b-d for their scanning electron microscopy (SEM) images) and experimentally characterized their reflection  Supplementary  Information). The excellent agreement among the FEM, the experimental and our theoretical results unambiguously justify our theory.

Implementations of the theory in lineshape tailoring
We now apply our theory to "design" the lineshape of a photonic system. Figure 3d shows that the interresonator coupling can dramatically change the lineshape of a coupled system, essentially determined by the two "dressed" modes with frequencies and bandwidths fω ± ;Γ ± g. Therefore, we must first understand the properties of the dressed modes fω ± ;Γ ± g.
Consider a two-mode two-port system with two resonators placed on the same plane illuminated by a normally incident wave. Assuming Γ a 1 ¼ Γ a 2 ¼ Γ a for simplicity, we can explicitly rewrite Eq. (7) as Diagonalizing the matrix containing t by an orthogonal transformation M, we obtain the following equation describing the amplitudes of two collective modesã ± : and ΔΓ ¼ Γ 1 À Γ 2 , and a þãÀ ð Þ T ¼ M a 1 a 2 ð Þ T . Since an orthogonal transformation does not change the trace of a matrix, it is sufficient to study Δω ¼ω þ Àω À and ΔΓ ¼Γ þ ÀΓ À , which are determined by t, Δω, and ΔΓ via Here, and in what follows, we have scaled all involved physical quantities (i.e., Δω, ΔΓ, Δω, ΔΓ, and t) by ffiffiffiffiffiffiffiffiffi Γ 1 Γ 2 p to make them dimensionless. Equation (16) shows that even for two resonators with fixed properties, one can still use the interresonator coupling t to change the properties of the "dressed" modes and, in turn, "design" the final lineshape of the coupled system.
The left and right panels in Fig. 4a depict, respectively, how Δω and ΔΓ vary with Δω and t, with ΔΓ set at two different values. We find that while Δω exhibits circular equal-value lines on the Δω~t plane independent of ΔΓ, ΔΓ, exhibits fascinating behavior on the Δω~t plane depending sensitively on ΔΓ. In particular, on each Δω~t phase plane with a fixed ΔΓ, we always find two special lines, defined as ΔΓ ¼ 0 (red lines) and ΔΓ ¼ ± ðΓ 1 þ Γ 2 Þ (green lines), to separate the whole space into four subregions with distinct properties. Physically, while the  Fig. 3 Benchmark test of our theory on a realistic system. a Schematic of the coupled plasmonic system under study. Here, the geometrical parameters are p = 530, d = 30, w = 240, l = 420, R = 110, and a = 85, all in units of nm. b-d Reflectance spectra of periodic metasurfaces containing b bar resonators only, c C resonators only, and d the two resonators coupled together, obtained by our theory (solid lines), FEM simulations (circles), and measurements (triangles). White dashed lines and gray areas denote the frequencies and widths of the resonant modes. The right panels of c and d are SEM images of the fabricated samples with scale bars (white lines) of 500 nm condition ΔΓ ¼ 0 implies that the two dressed modes have identical bandwidths (i.e.,Γ þ ¼Γ À ), the other condition ΔΓ ¼ ± ðΓ 1 þ Γ 2 Þ means that one dressed mode exhibits vanishing radiative damping. Interestingly, these two phase boundary lines rotate as ΔΓ changes, as shown in Fig. 4b.
To illustrate the key features of the four subregions, we purposely choose eight points from a circle on the Δω~t plane with ΔΓ = 2 (see Fig. 4a) and illustrate in Fig. 4c how the reflection spectra of the corresponding systems evolve. Consistent with our expectations, the spectra of systems 1 and 5 only exhibit one peak, as the other mode is completely dark, while the spectra of systems 3 and 7 exhibit two peaks with equal bandwidths. In between these special points, the spectra gradually evolve. Notably, the radiation damping (bandwidths) of the two "dressed" modes can vary continuously from 0 to Γ 1 + Γ 2 , while moving on the circle (see Fig. 4d).
The physics is very clear: now that the dressed modes are appropriate linear combinations of two original modes, their radiation damping must also be linear combinations of that of the two original modes. Therefore, varying Δω and t can dramatically modify the relative portions of the two original modes in constructing the dressed modes and, in turn, efficiently control the radiation damping of the dressed modes. In principle, one can realize any desired lineshapes based on our phase diagram by choosing certain original modes and "tuning" the coupling t. Of particular interest is the appearance of a purely dark mode with infinitely long lifetime, which shares the same physical origin as the BIC and has many interesting applications [49][50][51] .
We now experimentally verify our predictions on lineshape tailoring based on coupled systems constructed by the two resonators studied in Fig. 3a. Since t is solely determined by the overlap between the ψ NF m of two resonators (see Eq. (9)), we understand that changing the resonators' relative configuration can dramatically modify t. Indeed, as we rotate the C-shaped resonator with respect to the bar resonator, we find that t drastically changes (see solid line in Fig. 5a). In particular, increasing the relative angle θ between two resonators can drive t to change from a positive value to a negative value, passing through 0 at a particular angle. Such an intriguing t~θ relation can be simply explained by an effective model for plasmonic coupling established previously 47,48 . Choosing six points on the t~θ curve, as shown in Fig.  5a, we employ our theory to study the optical lineshapes of their corresponding realistic systems. Since the two original modes have fixed properties, these six systems with different t are located on a straight line in the phase diagram passing through two phase boundaries (see Fig.  5b). Their reflection spectra, computed by our theory, are depicted in Fig. 5c as solid lines, exhibiting the expected behaviors. In particular, the spectrum of the third system only exhibits one peak, while that of the fifth system contains two equal-bandwidth peaks, consistent with the phase diagram shown in Fig. 5b. Once again, we emphasize that all spectra are calculated with our theory directly and without any fitting procedures.
We then perform both experiments and simulations to verify the above theoretical predictions. We fabricate samples according to the designs using the standard EBL method, with the right panel in Fig. 5c showing SEM images of the fabricated samples. Illuminating these samples with normally incident light withẼ kŷ, we measure their transmission/reflection spectra and depict the reflection spectra as solid triangles in Fig. 5c. We also perform FEM simulations to calculate their reflection spectra (open circles in Fig. 5c). Both the experimental and simulation results are in excellent agreement with the spectra obtained by our theory (solid lines in Fig. 5c). In particular, the measured/simulated spectra of sample 3 exhibit clear BIC features, while those of sample 5 contain two peaks with equal bandwidths. We also employ our theory to predict the transmission spectra of these systems, which are in excellent agreement with the measured and simulated results (see Sec. VIII in the Supplementary  Information).
The solid line in Fig. 5d depicts how varying t significantly modulates the radiative Q factor of the lowfrequency dressed mode, as predicted by our theory. That the Q factor diverges at a specific point signifies the appearance of a BIC. The symbols are the Q factors of six realistic samples obtained by analyzing their measured reflection spectra. Excellent agreement is noted between the experimental and analytical results. At the frequency where the BIC appears, the radiations from the two individual resonators exactly cancel each other, leading to vanishing of the total radiative loss (see Sec. VI in the Supplementary Information).

Discussion
In summary, we have derived a formal theoretical framework directly from Maxwell's equations to study the optical responses of arbitrarily coupled photonic systems, in which all involved parameters are unambiguously computable without any fitting procedures. After testing it against both Mie theory and numerical simulations on different systems, we illustrate how to employ it to design the lineshape of a coupled system by modulating the couplings between resonators. In particular, we show that one can always choose a specific coupling between two arbitrary resonators to make one of the "dressed" modes in the coupled system completely dark, creating a BIC. All predictions are quantitatively verified by our experiments and simulations at near-infrared wavelengths. In addition to revealing the profound physics underlying the coupling-induced phenomena, our theory also offers a powerful tool to design optical devices with wellcontrolled NF and FF properties, and can be extended to study coupled systems for other types of waves.

Simulations
We employed FEM simulations using the commercial software COMSOL Multiphysics. The permittivity of Ag was described by the Drude model εðωÞ ¼ ε 1 À ω 2 p ωðω þ iΓ e Þ , with ε 1 ¼ 5ε 0 , ω 0 ¼ 0 THz, and ω p ¼ 2π 2176:2 THz. The effective damping rate was set as Γ e ¼ 2π 38:3 THz for the bar structure and Γ e ¼ 2π 27:3 THz for the Cshaped resonator, obtained by fitting with our experimental results. The SiO 2 spacer was considered a lossless dielectric with permittivity ε = 2.42. Additional losses caused by surface roughness and grain boundary effects in thin films, as well as dielectric losses were effectively considered in the fitting parameter Γ e .

Fabrication
All our meta-devices were fabricated following standard EBL and lift off processes. First, the positive resist was successively spin coated on a silica substrate, and exposed with EBL (JEOL 6300) with an acceleration voltage of 100 kV. After exposure, the samples were developed in the solution of isopropanol alcohol and methyl isobutyl ketone. Then, 3 nm Cr and 30 nm Au/Ag were deposited using electron-beam evaporation. Finally, the top patterns were formed after a lift of process. All samples had dimensions of 80 µm × 80 µm.

Optical characterizations
We used a homemade macroscopic spectrometer equipped with a broadband supercontinuum white light source and a fiber-coupled grating spectrometer (Ideaoptics NIR2500) to characterize the optical properties of fabricated samples (see more details in Sec. VII of Supplementary Information).