Thickness Insensitive Nanocavities for 2D Heterostructures using Photonic Molecules

Two-dimensional (2D) heterostructures integrated into nanophotonic cavities have emerged as a promising approach towards novel photonic and opto-electronic devices. However, the thickness of the 2D heterostructure has a strong influence on the resonance frequency of the nanocavity. For a single cavity, the resonance frequency shifts approximately linearly with the thickness. Here, we propose to use the inherent non-linearity of the mode coupling to render the cavity mode insensitive to the thickness of the 2D heterostructure. Based on the coupled mode theory, we reveal that this goal can be achieved using either a homoatomic molecule with a filtered coupling or heteroatomic molecules. We perform numerical simulations to further demonstrate the robustness of the eigenfrequency in the proposed photonic molecules. Our results render nanophotonic structures insensitive to the thickness of 2D materials, thus owing appealing potential in energy- or detuning-sensitive applications such as cavity quantum electrodynamics.


INTRODUCTION
Photonic crystals have spatially periodic dielectric constants that facilitate the engineering of photonic bandstructure via geometry [1]. By introducing defects into the periodic structure, photonic confinements can be generated for cavities, waveguides and directional couplers [2,3]. A single-mode photonic crystal cavity with a mode volume of the order of ∼ (λ/n) 3 can be used to enhance the strength of light-matter coupling towards the quantum limit [4], whilst allowing individual components to be combined into integrated photonic circuits.
The conventional fabrication methods used to realize photonic crystals is to etch a periodic array of air holes in a high-index semiconductor slab such as Si and GaAs [5][6][7]. By improving the design [8] and fabrication methods [9], such conventional photonic crystal cavities have been optimized towards the theoretical limit. For example, Qfactors have been pushed into ultra-high (> 10 5 ) regime, and the experimentally measured performance can reach the theoretical calculations [10]. Such a cavity is ideal to study and control light-matter interactions by directly embedding quantum emitters such as quantum dots [11][12][13] into the dielectric structure. This type of cavityquantum-dot system has been widely applied in quantum devices [14,15], and extended to study multi-modal lightmatter interactions such as multi-photon processes [16,17], electromagnetic controllability [18][19][20] and phononic effects [21,22].
However, conventional photonic crystal cavities are less suitable for integration with 2D atomically thin semiconductors. This is since it is non-trivial to embed 2D semiconductors into the dielectric slab during growth, their properties are strongly influenced by the dielectric disorder, and encapsulation with an ultra-flat, inert insulator is essential for high-performance 2D devices [23][24][25][26][27]. Hexagonal boron nitride (hBN) is a layered 2D insula-tor and has proven to not only be an ideal substrate for other 2D materials, but also host a variety of defect quantum emitters having potential interests for nanophotonics. The integration of hBN in photonic crystals raises new challenges in the nanofabrication [28,29]. Meanwhile, in contrast to conventional cavity materials such as Si and GaAs whose thickness can be precisely controlled during growth, the thickness of 2D materials is difficult to control and even accurately identify. The indeterminate thickness of 2D heterostructures will result in shifts in the frequency of 2D-material cavity from the designed value [28,30]. Recently, we have proposed a hybrid nanobeam cavity to solve the problem of etching of 2D heterostructures [31,32], but the problem of 2D-material thickness fluctuation has not been fully investigated.
Here, we report several approaches to realize a robust resonance frequency with 2D heterostructures using photonic molecules. Such a photonic molecule consists of two coupled cavities and usually host coupled modes with spatially symmetric (S) and antisymmetric (AS) character arising from the mode coupling [33][34][35][36]. For a single cavity, the resonance frequency generally shifts linearly with the thickness of the 2D heterostructure. Therefore, we exploit photonic mode coupling to achieve nonlinear shifts of the eigenfrequencies of coupled cavities. In this case, the linear slope of frequency shift is suppressed to zero at the local inflection points, enabling a resonance frequency robust to the thickness of the 2D materials. Based on the coupled mode theory, we find that such nonlinear eigenfrequencies can be achieved either using a homoatomic molecule with frequency filtered coupling or using heteroatomic molecules around the diabolical point. These predictions are further demonstrated quantitatively using 2D and 3D finite-difference time-domain (FDTD) calculations.

COUPLED MODE THEORY
For a cubic cavity, the wave vector of the cavity mode is proportional to the inverse of the spatial extent, and thereby, the resonance frequency monotonically decreases with the increasing cavity size. This is generally the same to a single nanocavity based on conventional photonic crystal structures. We consider the bare frequency f 1 of a single cavity that varies linearly with the hBN thickness t hBN . Thus, the aim is to construct a nanophotonic molecule with eigenfrequencies f ± that vary non-linearly and nonmonotonically with t hBN . As such, the first order derivative (linear slope) of f ± to t hBN is zero at the minima or maxima points, which means the resonance frequency is robust within a range of hBN thickness.
We consider two individual cavities with the bare frequency f 1 (2) and the photon loss rate γ 1 (2) . The molecule is described using the coupled mode theory [37] according to 1 2π where a 1(2) is the field amplitude of two cavities and g is the coupling strength between them. Then the resonance frequencies f ± are calculated by solving the eigenstates of Eq. 1 as where ∆f = f 1 −f 2 and ∆γ = γ 1 −γ 2 . For brevity, in this work, we define all parameters in frequency units, i.e., the coupling strength is hg when converted to the energy unit where h is the Planck constant. For the results of coupled mode theory based on Eq. 1 and 2, we omit the unit THz for brevity.
By solving Eq. 1 and 2, we obtain several cases having nonlinear and nonmonotonic f ± , as presented in Fig. 1. We set cavity C1 (denoted by red) the frequency to f 1 = −t hBN and the loss rate to γ 1 = 2, and take it as the one that we want to suppress the sensitivity of the frequency to the thickness. As schematically shown in Fig. 1(a), we first consider a homoatomic molecule, meaning that the cavity C2 (denoted by blue) is the same as C1. In this case, the two eigenfrequencies are f ± = f 1 ± g, with a splitting of 2g. Therefore, to achieve a nonlinear f ± we need a nonlinear g. One possible way is to introduce an optical filter [38,39], which allows only the photons within certain frequencies to exchange between the two cavities. Then it is expected that the coupling factor g is frequency dependent, e.g., following a Gaussian filter as where σ g is the bandwidth of the filter and g 0 is the coupling strength at f 0 . In this case, we obtain that when g 0 = √ eσ g , eigenfrequencies f + at t = f 0 − σ g and f − at t = f 0 + σ g have a zero value for both first-order (linear slope) and second-order derivative which means robust eigenfrequency regimes. Fig. 1(b) presents the eigenfrequencies with g 0 = 1.648, f 0 = 0 and σ g = 1 for example. The purple line is g(f ), the light red line is the bare frequencies of single cavity f 1 , and the solid red and blue lines are the eigenfrequencies f ± . As shown, there are regions with robust f + or f − , denoted by the gray dashed arrows.
The homoatomic molecule depicted schematically in Fig. 1(a) requires a filtered photonic mode coupling. In contrast, direct coupling in the free space between two cavities is much easier to achieve from a technological perspective. In this case, the coupling strength g is determined by the separation of the two cavities, whilst the frequency dependence is not significant. As discussed above, if g is nearly constant, the eigenfrequencies f ± of a homoatomic molecule will still vary linearly. Therefore, we construct heteroatomic molecules to achieve nonlinear f ± in the free space, as schematically shown in Fig.  1(c). The cavity C1 (red) is identical to that used in the homoatomic case which has f 1 = −t hBN and γ 1 = 2. We set the cavity C2 (blue) as a lower Q-factor with the larger loss rate γ 2 = 20. Eq. 1 and 2 indicate that, when the coupling strength exceeds the decay rate g > ∆γ/4, the molecule operates in the strong coupling regime with a splitting. In contrast, in the weak coupling regime g < ∆γ/4, it operates without an anticrossing [40]. Moreover, the eigenfrequencies f ± exhibit strong non-linearity and non-monotonicity around the diabolical point g = ∆γ/4, the degenerate state at the threshold between weak and strong coupling [41,42]. Fig. 1(d) shows the calculated eigenfrequencies with f 2 = −2t hBN , meaning that the low-Q cavity C2 shifts more rapidly than the high-Q cavity C1 with the hBN thickness. The light lines (red and blue) are the bare frequencies f 1 (2) , whilst the dotted lines (red and blue) are the eigenfrequencies f ± at the diabolical point g = 4.5. As shown, the high-Q mode (red dotted line) exhibits strong non-linearity and non-monotonicity. In contrast, the solid lines (red and blue) are the case of g = 4.24 which is slightly below the diabolical point. Since the system is in a weak coupling regime, the superposition of the two cavities is not significant, and thereby, cavity C1 dominates in the spatial profile of the high-Q mode. Moreover, in this case, the first-order (linear slope) and second-order derivative of the high-Q mode frequency (red line) to t hBN are both zero at the resonance t hBN = 0, thus generating robust eigenfrequency regimes denoted by the gray dashed arrows in the figure. Fig. 1(e) presents the results with f 2 = 0, which means the low-Q cavity C2 is not an hBN cavity. The light lines (red and blue) are the bare frequencies f 1 (2) , whereas the dotted lines (red and blue) represent the situation at the diabolical point g = 4.5. For f 2 = 0, we cannot use a g which is slightly below the diabolical point, because the high-Q mode (red) blue shifts at t hBN < 0 whilst redshifts at t hBN > 0. As a result, the indeterminacy (slope) around the resonance is enhanced rather than suppressed. In contrast, the robust eigenfrequencies are achieved by g slightly above the diabolical point, as shown by the case of g = 4.6 (solid red-blue lines) in Fig. 1(e). In this case, the eigenfrequencies exhibits strong non-monotonicities, which provides the minima and maxima points with zero linear slopes. Meanwhile, the mixture of red and blue color around the resonance t hBN = 0 denotes that the eigenstates are the superposition of modes from both two cavities, which means their spatial profile still involves cavity C1, the one we want to suppress the indeterminacy. Therefore, the case of g = 4.6 in Fig. 1(e) is also valid to achieve robust resonance frequencies.

2D-MATERIAL CAVITIES AND MODE COUPLINGS
The excitonic properties of 2D semiconductors are degraded when placed into an inhomogeneous environment such as dielectric disorder [25]. In this case, the exciton linewidth is broadened, and some novel single-photon emitters will be suppressed [43,44]. Therefore, 2D heterostructures with full hBN encapsulation is necessary to achieve the pristine excitonic and optical qualities [23,24]. Nowadays, the most common method used to build 2D heterostructures and integrate them into a nanocavity is dry viscoelastic stamping followed by transfer on top of pre-fabricated photonic cavities [45][46][47][48]. However, the attachment limits the cavity-matter overlap and decreases the cavity Q-factor. Another way is to directly fabricate nanophotonic structures using hBN [28,29], with the problem that hBN cannot readily be etched [49]. To solve these problems, we recently proposed a hybrid nanobeam cavity design [31,32] that integrates hBN into nanocavities whilst avoiding etching of nanostructures into hBN layers. As presented in Fig.  2(a), all nanoscale trenches are etched in Si 3 N 4 using mature fabrication technologies, whilst the hBN layer is not perforated. Such nanocavities have been experimentally demonstrated to simultaneously achieve pristine excitonic quality, high Q-factor and large cavitymatter overlap, providing an ideal platform to explore the fundamental physics of 2D materials and their lightmatter interaction in nanophotonic systems [50][51][52].
As presented in Fig. 2(b), the cavity photonic confine- ment is formed by locally chirping the photonic crystal periodicity around the center point. The distance between nanoscale trenches (periodicity) follows a Gaussian profile where a i is the ith distance schematically shown in Fig.  2(a), a is the lattice constant, and A, σ defines the Gaussian profile. We fix the nanoscale trenches with the length h x = 120 nm, depth h z = 150 nm and lattice constant a = 270 nm, whilst the Si 3 N 4 has a thickness of 200 nm. The designed hBN thickness t hBN is 50 nm. At the bottom of Fig. 2(a), we present typical modes for two cavities with A = 0.1, σ = 4 and A = −0.1, σ = 4 calculated using 2D FDTD. For A = 0.1, σ = 4, the fundamental mode M1 is at 407 THz with a Q-factor of 4.8 × 10 6 and the second mode M2 is at 393 THz with a Q-factor of 1.8 × 10 4 . For A = −0.1, σ = 4, the fundamental mode M1 is at 424 THz with a Q-factor of 5.9 × 10 3 and the second mode M2 is at 440 THz with a Q-factor of 820.
Since the refractive index of hBN is not well established [53,54], we use a fixed value of 2.0 for both hBN and Si 3 N 4 , while the mesh size used was 10 nm, providing a precision better than 0.01%. We observe that for A > 0, the frequency of mode M2 is lower than M1, while it is higher when A < 0. This different mode frequency sequence arises from the differ-ent photonic bandgap as presented schematically in Fig.  2(b). For A < 0, the periodicity at the center is increased. The larger periodicity means a longer wavelength (lower frequency), thus the bandgap confinement is concave at the center. Thereby, the higher order mode (larger mode volume) has a higher frequency. Vice versa, for A > 0 the bandgap confinement is convex, and thereby, the higher order mode has a lower frequency. The cavity with A < 0 has a lower Q-factor, and the field mainly distributes in the air gaps. Therefore, the convex bandgap confinement is usually used in 1D nanobeam photonic crystal cavities [47,48].
The photonic molecule can be constructed by coupling two cavities within one nanobeam or in two separated nanobeams. For the case within one nanobeam, we define the periodicity profile as where the center point of the two cavities are at xc 1 (2) with the Gaussian profile A 1(2) and σ 1 (2) . The coupling strength is controlled by the distance between them ∆xc = xc 1 −xc 2 , as shown in Fig. 2(c). For the case with two separated nanobeams, we set xc 1 = xc 2 = 0, and the coupling strength is controlled by the distance between two nanobeams ∆z as shown in Fig. 2(d). In Fig.  2(c)(d) we present the distance-dependent spectra calculated using 2D FDTD, in the case of molecules with one or two nanobeams, respectively. A homoatomic profile A 1(2) = 0.2 and σ 1(2) = 2 are used. The symmetric (S) and antisymmetric(AS) modes have opposite frequency sequences in the two cases since the bandgap confinement in the x direction along the nanobeam is convex for A > 0 whilst concave in the z direction. The S and AS modes of the photonic molecule can be generally described by the coupled mode theory in Eq. 1 and 2, but it is not strict. For example, it predicts the homoatomic molecule with f 1 = f 2 to have a symmetric splitting with eigenfrequencies f ± = f 1 ± g and the loss rate of S and AS modes remains γ 1 . However, in the FDTD calculation results, e.g., the case in the upper panel of Fig. 2(c) which corresponds to ∆xc = 9, the Qfactor is 1.2×10 5 for AS mode while 5.3×10 5 for S mode. This is due to the limitations of the coupled mode theory. Firstly, in Eq. 1 and 2 we use a real g to describe the coupling, but in reality, this term could also involve an imaginary part which results in a non-Hermitian system [36,55]. The imaginary coupling results in different loss rates of S and AS modes. Secondly, the coupled mode theory is valid for individual cavities that couple to each other through the evanescent field. When the distance between the two cavities is large, i.e., ∆xc > 8 in Fig.  2(c) and ∆z > 500 nm in (d), we observe a nearly symmetric splitting which indicates that the two cavities are approximately individual. In contrast, the splitting becomes asymmetric as the distance decreases, because the two confinements merge into one and the two cavities are no longer individual. The comparison between individual and non-individual cases is schematically presented in the inset of Fig. 2(c). Nonetheless, the calculation results in Fig. 2(c)(d) reveal that for both types of the photonic molecule, when the coupling strength g is smaller than 2.5 (corresponding to ∆xc = 9 or ∆z = 550 nm), the eigenfrequencies split symmetrically thus can be well described by the coupled mode theory.
Before exploring the nonlinear eigenfrequencies for photonic molecules, we first study the hBN thickness dependence of the bare frequency of a single cavity. Usually, apparent optical contrast is used to quickly identify the thickness of hBN flakes [56], but the error of this method remains ± several layers, corresponding to a thickness of ∼ 2 nm. In Fig. 3(a) we present the hBN thickness (t hBN ) dependence of cavity modes calculated using 2D FDTD. For several different values of A and σ, the cavity mode exhibits a similar frequency shift slope of -0.6 THz/nm. This means that in the 2D space, it is difficult to have two cavities with remarkably different t hBN dependent frequency shifts. Therefore, to achieve the bare frequencies as predicted in Fig. 1(d) of the two cavities, we extend the system to 3D space. The dry etching of hBN normally has an etching angle that can be optimized by the recipe [49], but it is difficult to make it completely vertical. Due to the non-vertical characteristic, the hBN thickness t hBN also impacts the nanobeam width d ′ y . In this work, we set the etching angle to be 45 • . As schematically shown in the upper panel of Fig. 3(b), for the same designed width d y (resist beam width), if t hBN increases 5 nm, the actual nanobeam width d ′ y will increase by 10 nm. The effect of increasing d ′ y from 100 to 110 nm is obviously different from that of 400 to 410 nm. Therefore, we use a high-Q thick cavity with the designed nanobeam width d y of 400 nm, and a low-Q thin cavity with d y of 100 nm, to achieve remarkably different frequency shift slopes. The frequencies and Q-factors of the two cavities calculated using 3D FDTD are presented in Fig. 3(b). As shown, the high-Q cavity with A = 0.25, σ = 6, d y = 400 nm has the frequency shift slope of −0.8 THz/nm, whilst the low-Q cavity with A = 0.15, σ = 1, d y = 100 nm has the slope of −1.7 THz/nm. The 2D FDTD results in Fig. 2(c)(d) reveal that the coupling strength of the photonic molecule can be controlled by the distance in x and z direction. Similarly, the coupling can also be controlled by the distance in y direction when we extend the system to 3D space [31]. The 2D FDTD is a valid approximation of the 3D system, when the two cavities have the same nanobeam width d y in the y direction. Based on the results in Fig. 3, we obtain that different d y in 3D space is only necessary to achieve remarkably different shifts of bare frequencies, which corresponds to Fig. 1(d). In contrast, uniform d y is enough to achieve the bare frequencies in Fig. 1(b)(e). Therefore, we use 2D FDTD to demonstrate the robust eigenfrequencies for the cases of Fig. 1(b)(e) and 3D FDTD for the case of Fig. 1(d). From simple to complex, we next first show the case of Fig. 1(e) in Fig. 4, then the case of Fig. 1(d) in Fig. 5 and finally discuss the case of Fig. 1(b) in Fig. 6.  Fig. 1(e).

HETEROATOMIC MOLECULES
We firstly discuss the heteroatomic molecule comprising an hBN/Si 3 N 4 hybrid cavity and a second cavity without hBN, corresponding to the case predicted in Fig.  1(e). Here we approximate the photonic molecule in the xz 2D space as shown in Fig. 4(a), and control the coupling g between the two cavities by varying the distance ∆z. This approximation is of course not strict, but the key point, i.e., varying the distance between the two cavities to achieve a photonic molecule around the diabolical point, is the same for the real 3D structure and the 2D approximation.
We use A = 0.2, σ = 4 for the high-Q cavity C1 and A = −0.3, σ = 1 for the low-Q cavity C2. Since the coupled mode theory only applies to relatively large distance ∆z > 500 nm (Fig. 2(d)), we first calculate the hBN thickness t hBN dependent spectra at ∆z = 600 nm and present the result in Fig. 4(b). The four sharp peaks shifting significantly with t hBN are from the high-Q cavity C1, and the two broad peaks exhibiting slight shifts are from the low-Q cavity C2. As shown, the mode M2 of cavity C1 (C1M2) exhibits a clear anticrossing to the mode M1 of cavity C2 (C2M1). The mode C1M2 has a bare frequency of 415 THz and a Q-factor of 2.5 × 10 4 at the designed thickness t hBN = 50 nm, and the mode C2M1 has a bare frequency of 441 THz and a Q-factor of 350. Increasing the distance ∆z results in a reduction of the coupling strength g between these two modes. As presented in Fig. 4(c), the splitting at resonance decreases with increasing ∆z, and the system finally reaches the  Fig. 3(b). (b) Calculated spectra with the distance ∆y = 580 and 680 nm. Inset shows the frequency of high-Q mode at ∆y = 680 nm, which agrees well to the prediction in Fig. 1(d).
state slightly above the diabolical point as predicted in Fig. 1(e). The inset in Fig. 4(c) shows that the resonance frequency of AS mode varies by 0.14 THz within the indeterminacy 2 nm of t hBN (17 − 19 nm), one order of magnitude smaller than the case in a single cavity (1.2 THz). As discussed in the context of Fig. 3(b), the frequency shift slope of a single cavity can be controlled by varying the nanobeam width d y in 3D space. Therefore, we construct a photonic molecule based on these two cavities. Hereby, we use A = 0.25, σ = 4, d y = 400 nm for the high-Q cavity C1 and A = 0.15, σ = 1, d y = 100 nm for the low-Q cavity C2. The mode coupling is controlled by the distance ∆y as schematically shown in Fig. 5(a). In Fig. 5(b) we present the spectra with different ∆y calculated using 3D FDTD. The resonance point of the two cavities is at t hBN = 61.5 nm, consistent with the bare frequencies of single cavities in Fig. 3(b). When ∆y = 580 nm, we observe a significant anticrossing indicating that the system is in the strong coupling regime. In contrast, when ∆y = 680 nm, the coupling strength reduces and the anticrossing is suppressed. Around the resonance t hBN = 61.5 nm, the t hBN dependent frequency shift of the high-Q mode is suppressed, as presented in the inset in Fig. 5(b). The 3D calculation results in Fig.  5 agree well with the state slightly below the diabolical point predicted in Fig. 1(d).

HOMOATOMIC MOLECULES
In Fig. 1(b) we predict that the homoatomic molecule with a filtered mode coupling g(f ) is another method to achieve nonlinear resonance frequency and improve the robustness to t hBN . The optical filter is usually realized by integrating material with nonlinear absorption allowing only the photons at certain frequencies to pass. However, the simulation of such complex materials in the FDTD calculation is non-trivial and we qualitatively explore the homoatomic molecule using a frequency selective coupler. As discussed in the context of Fig. 3, 2D FDTD is a valid approximation to describe the case of Fig. 1(b), and thereby, we propose the nanophotonic device defined in the xz plane as presented in Fig. 6(a). The two nanobeam cavities are located on each side, and an array of the 2D hexagonal photonic crystal lattice is in the middle. Only the photons within the photonic band of the hexagonal photonic crystal are transmitted and exchanged between the two cavities. Therefore, the hexagonal photonic crystal serves as a frequency-selective filter for the mode coupling.
We set A = 0.2, σ = 4 for the two nanobeam cavities. The lattice constant of the hexagonal photonic crystal is set as a h = 270 nm, and we vary the radius of the air holes r to control the photonic band. The calculated photonic bands are presented in the upper panels in Fig.  6(b), and the dashed lines denote the lower edge of the bandgap. This value is 392, 405 and 426 THz for the radius r of 80, 90 and 100 nm, exhibiting a positive dependence as expected. The calculated t hBN dependent spectra are presented in the bottom panels in Fig. 6(b). Multiple peaks are observed because the nanophotonic structure in Fig. 6(a) is very complex. Here we mainly focus on the peak with the largest intensity and its splitting, which are denoted by the black arrows on top of the spectra. The splitting occurs when the hBN thick-ness t hBN exceeds a threshold. Moreover, this threshold in the cases of r = 80, 90 and 100 nm corresponds to the frequency of 392, 402 and 415 THz, respectively. These values are consistent with the band of the hexagonal photonic crystal denoted in the upper panels in Fig. 6(b), indicating that the hexagonal photonic crystal is a valid filter for this photonic molecule.
We note that the nanophotonic structure in Fig. 6(a) is a simple example to qualitatively explore the filtered coupling. We use the distance of 900 nm because larger distance will suppress the coupling into the weak regime as shown in Fig. 2(d). As such, the hexagonal photonic crystal has only three rows of air holes, which is not enough for a highly efficient filter. These limitations result in the relatively weak non-linearity, i.e., the t hBN dependent shift of the splitting peak is reduced but not enough to obtain a robust (zero linear slope) resonance frequency. A potential optimization approach is introducing ∆xc between the two nanobeam cavities. In this case, we can control (through ∆xc) the wave vector direction of the photons that exchange between the two cavities, and match the wave vector to some points in the Γ − M − K − Γ photonic band. Thereby, the hexagonal photonic crystal will not only filter the frequency but also the wave vector, and thereby, the efficiency will be improved.

CONCLUSION
In summary, we have demonstrated that the nonlinearity of the mode coupling in the photonic molecules results in coupled mode eigenfrequencies that are robust to the thickness of 2D heterostructure. Although atomically thin semiconductors have well-defined emission energies compared to the random energy of traditional quantum dots, the indeterminate thickness of their hBN encapsulation leads to indeterminate cavity resonance frequencies. For a single cavity, the frequency shift of 0.6 THz/nm is obtained in the 2D FDTD calculations, and a value of 1.7 THz/nm is obtained in the 3D calculations. By comparison, the linewidth of the free excitons in the atomically thin semiconductors such as monolayer transition metal dichalcogenides at low temperatures is ∼5 meV (1.2 THz) and that of the localized excitons is 0.1 − 0.2 meV (∼50 GHz). As a result, it is difficult to achieve the light-matter resonance directly in a single cavity, and complex external methods have to be used to control the detuning [31,52]. Our results reveal that photonic molecules are capable of exhibiting thickness-independent resonance frequencies, thus having great potential in the 2D-material nanophotonics, especially the detuning sensitive applications such as lasing and exciton-photon polaritons.

RESEARCH FUNDING
S. Y. and P. J. acknowledge support from the National Science Fund for Distinguished Young Scholars (52225507), the National Key Research and Development Program of China (No. 2021YFF0700402), and the Fundamental Research Funds for the Central Universities. J. F., C. Q. and P. J. gratefully acknowledge the German Science Foundation (DFG) for financial support via grants SPP-2244, as well as the clusters of excellence MCQST (EXS-2111) and e-conversion (EXS-2089). C. Q. gratefully acknowledges the Alexander v. Humboldt foundation for financial support in the framework of their fellowship programme. P. J. acknowledges support from the China Scholarship Council (202006280067).

AUTHOR CONTRIBUTIONS
C. Q., J. F. and S. Y. conceived the project. P. J. and C. Q. performed the calculations. All authors discussed the results and wrote the manuscript.

CONFLICT OF INTEREST
Authors state no conflict of interest.

DATA AVAILABILITY
The datasets generated during and/or analyzed during the current study are available from the corresponding authors on reasonable request.