Anisotropic mode excitations and enhanced quantum interference in quantum emitter-metasurface coupled systems

This study proposes a nanophotonic structure that supports quantum interference (QI) between orthogonal decay channels in multilevel quantum emitters within the framework of the quantum master equation. The Green functions of the electric field are obtained by applying boundary conditions in the presence of two-dimensional metasurfaces. We demonstrate distinct in-plane excitation features of the surface plasmon modes (SPMs) with the anisotropic metasurfaces tailored to conductivity components. In particular, we observed that the Purcell factor of transitions with orthogonal polarizations experiences unequal enhancements, owing to the anisotropic propagation of the SPMs. This property depends only on the anisotropy of the metasurfaces; thus, it is easily manipulated. Using this platform and considering experimentally achievable material parameters, we predict a strong interference effect in three-level quantum emitters. In principle, this enables the generation of maximum QI. Our study provides a method for realizing QI systems and has potential applications in highly integrated, tuneable quantum devices.


Introduction
Tailoring photonics for monolithic integration enables the fabrication and control of electronic-photonic systems at the nanoscale [1][2][3]. In the last decade, light-matter interaction has been widely studied in two-dimensional layered structures [4]. Generally, most of these structures support the excitation of surface plasmons that exhibit small propagation losses and large field confinements [5,6]. It has been reported that the hybrid polariton modes of isotropic graphene sheets are electrically [7] or chemically [8] tuneable, which have recently gained significant research attention in application fields such as plasmonic waveguides [9], tuneable terahertz platforms [10], and photodetectors [11]. Motivated by these intriguing features, many studies focus on the manipulation of polariton modes. Among them, two-dimensional anisotropic metasurfaces (TDAMs), such as graphene nanoribbons [12], black phosphorus [13,14], van der Waals materials [15], and structures based on lithography and etching technologies [16,17], are promising candidates for integrated quantum photonics devices because their optical responses can easily be manipulated in the laboratory. Owing to the peculiar optical properties of the above TDAMs, they have found applications in hyperlensing [18], cloaking [19], on-chip light manipulation [20], near-field radiative heat transfer [21], biosensors for molecule and virus detection [22], polariton lasing [23], optoelectronics [24], and electro-optic modulators [25].
In the quantum scheme, the spontaneous emission of individual atoms [26] via vacuum modes can occur due to vacuum fluctuations. Purcell proposed that the lifetime of quantum emitters (QEs) placed in a cavity generally suffers strong modifications [27]. More recently, various types of optical structures have been designed to control the spontaneous decay property of QEs, including metamaterials [28][29][30][31], waveguides or cavities [32][33][34][35][36], photonic crystals [37][38][39], and nanostructures [40][41][42]. This new paradigm unites quantum optics and nanophotonics, thereby yielding a plethora of applications, such as superradiance [43], squeezing in resonance fluorescence [44], indistinguishable single-photon generation [45], and entanglement enhancement in the quantum field [46]. It has been noted that an anisotropic photonic environment is a key requirement to produce a significant quantum coherence and interference effect in three-level QEs [47]. Agarwal theoretically predicted the QI between the orthogonal transition channels of a V-type Zeeman atom in a cavity formed by conducting plates [48], where the near-degenerate atomic transitions become (anti)parallel in the strong QI regime. This mechanism leads to many optical phenomena, such as lasing without inversion [49], protecting quantum entanglement [50], electro magnetically induced transparency [51,52], atomic location [53,54], and spontaneous emission cancellation [55]. Moreover, effective control of the QI effect has been demonstrated by coupling QEs to bulk metamaterials [56], metallic nanostructures [57], and periodic waveguide structures [58]. Compared with the structures mentioned above, TDAMs facilitate an atomic-scale thickness and high tuneability [59,60]; thus, they possess unique capabilities to control light-matter interactions [61,62]. Although still in the early stage, they have been extensively employed in the emerging field of quantum photonics [63]. In a recent study, Jha et al investigated the QI effect supported by metasurfaces made of subwavelength-scale nanoantennas [64]. The results indicated that polarization-dependent phase shifts of different decay channels can be realized by reasonably designing the dimensions of the nanoantennas. Consequently, the metasurface creates an anisotropic quantum vacuum and enables prominent QI in orthogonal levels of QEs placed in the far-field region.
In this study, we propose a new interference mechanism in hybrid exciton-plasmon systems composed of a V-type three-level QE and TDAM, where the formation of an anisotropic photonic environment strongly relies on the excitation feature of the SPM. For real materials, both the dispersion and propagation properties of the SPM [15,65] are easily controlled by manipulating the optical response of the metasurfaces. Here, we use a quantum-mechanical formalism for the electromagnetic field excitations [42]. Within this quantization scheme, the excitation features of the SPM and the spontaneous decay rates of the QE in the proximity of the metasurface can be numerically evaluated using the Green function method. We demonstrate the possibility of realizing anisotropic field excitations in the near-field region according to numerical and simulation results. Unequal transition-rate enhancements in orthogonal channels are illustrated depending on the dipole orientation with respect to the crystallographic axes. Consequently, we show that the anisotropic Purcell effect is highly tuneable by engineering the anisotropy of the metasurfaces. When the isofrequency contour exhibits extremely flattened patterns, the anisotropic Purcell effect becomes particularly prominent. Finally, we investigate the influence of material anisotropy on QI strength in orthogonal decay channels of the three-level QE. We demonstrate the tuneability of the QI effect in the proposed systems and predict the generation of maximum QI under particular conditions. The remainder of this paper is organized as follows. In section 2, we derive the dispersion relation of the SPM and numerically investigate the isofrequency contours. The simulation used to verify the accuracy of the numerical results is also described. In section 3, the Purcell factors with respect to the orthogonal in-plane transitions are investigated following the proper quantization scheme. Finally, we analyse the influence of surface anisotropy on the QI strength. Additionally, we consider real two-dimensional surfaces and discuss the experimental implementation of the proposed system. The conclusions are drawn in section 4. The appendices provide further details of the analyses presented in sections 2 and 3.

Mode excitations of the in-plane SPM for different surface topologies
We begin by considering the case of a linear polarized current dipole located in the upper half space of an infinitesimally thin, free-standing TDAM (with host medium dielectric permittivity and permeability equal to the vacuum values), which has a conductivity tensor that can be defined as In the case of isotropic surfaces, either transverse-magnetic (TM) or transverse-electric (TE) SPMs can be supported [66]; however, in the anisotropic case, both the TM-and TE-polarized SPMs can be supported simultaneously [17]. Moreover, in the near-field region, the spontaneous decay of QEs is primarily determined by the excitation features of the SPM. We inspect the plasmon dispersion relation to understand the behaviour of the SPM. Following the method mentioned in reference [67], the scattered Green functions are derived by introducing the electric Hertzian potential and imposing boundary conditions (see appendix A for more details). By determining the zero(s) of the denominator in the Green functions arising from equations (A.21)-(A.24), we obtain the dispersion relation satisfying where ω is the frequency of the current dipole, ε 1 and μ 1 are vacuum permittivity and permeability, k 1 = ω/c denotes the vacuum wave number, and β 1 = k 2 1 − k 2 represents the vertical component of the wave number in the source region. It is clear that the solution of equation (2) defines the wave front propagation of the SPM supported by the TDAM. However, in the anisotropic case, the wave front and propagation direction do not coincide. Instead, the energy-propagation direction of the SPM is typically defined by the group velocity ∇ k ω(k), where k = k x e x + k y e y (e x and e y are unit vectors along the x and y directions) describes the in-plane wave vector [68]. As demonstrated by equation (2), an analytical expression of the group velocity is not available, owing to the complicated dispersion relation. However, the definition of the group velocity indicates that the energy-propagation direction of the surface modes must be orthogonal to the isofrequency contours.
However, the dispersion relation actually predicts different topologies for the SPM tailored to the conductivity components. The elliptic topology arises when both the diagonal components of the conductivity tensor have the same signs in their imaginary parts, i.e., for Im[σ xx ], Im[σ yy ] > 0 or Im[σ xx ], Im[σ yy ] < 0. However, in the latter case, the surfaces can only support the excitation of TE-polarized SPM, thereby leading to a poor field confinement, which is of little practical interest to the topic of enhanced light-matter interactions. The hyperbolic topology arises when the TDAM becomes capacitive along one optical axis and inductive along the orthogonal one, which requires that the imaginary parts of the conductivity components satisfy the relation Im[σ xx ] · Im[σ yy ] < 0. In this case, both the TEand TM-polarized SPMs can be simultaneously supported by the surfaces. As we will discuss later in this section, the dispersion topologies for the SPM are primarily affected by the optical response of the surfaces. Under particular conditions, manipulation of the surface conductivities results in anisotropic mode excitations, dipole radiations, and QI enhancement.
To unveil the mechanisms behind the light-matter interaction mediated by the TDAM, we start by investigating the isofrequency contours for different surface topologies. Figures 1(a) and (b) illustrate the isofrequency contours in the elliptic dispersion case, for which the metasurfaces support only the TM-polarized SPM because the imaginary part of the conductivities is positive. It is clear that when one of them becomes much larger, for example, σ i xx σ i yy , the k surface exhibits a quasielliptic form that is elongated along the direction of the smaller conductivity component. This property is well illustrated by the red solid and green dashed curves, where the flattening of the k surface depends on the ratio of the imaginary parts between the diagonal components of the conductivity tensor.
It should be noted that the wave number of the supported SPM decreases owing to the increase in the optical response (the absolute values of the imaginary part of the conductivity components) of the metasurfaces. This is primarily due to the decrease in the field confinement of the SPM. Moreover, according to the dispersion relation described by equation (2), the branch points related to the different propagating plasmonic modes in the k x plane with TM and TE polarizations (denoted by the superscripts p and s) are k p y = k 1 1 − 4ε 1 /μ 1 σ 2 yy and k s y = k 1 1 − μ 1 σ 2 xx /4ε 1 , respectively. These expressions clearly imply that with an increase in the conductivity components, the wave number of the TM surface modes gradually decreases to the vacuum value k 1 ; thus, the field becomes poorly confined near the metasurface. Here, only the case of σ i xx < 0, σ i yy > 0 is considered, and the results are easily extended to another hyperbolic case with conductivity components that satisfy σ i xx > 0, σ i yy < 0. In contrast to the previous case, the isofrequency contour of the SPMs forms hyperbolas with larger wave numbers. With the assumption of large in-plane wave number for the excited SPM (k x , k y k 1 ), the solutions of equation (2) can be approximated as k y = ±k x |σ i xx |/ σ i yy , which indicate the plasmons with wave vectors on hyperbola asymptotes. This description implies narrow beams of the SPM carrying energy along the directions normal to the hyperbolas. In agreement with the qualitative analysis presented in figure 1(c), the surface anisotropy significantly modifies the energy flow of the hyperbolic plasmons.
As mentioned previously, hyperbolic metasurfaces support the simultaneous excitation of the TM and TE plasmon modes. The inset of figure 1(d) shows that the phase transition between these two critical plasmon modes occurs due to the manipulation of surface anisotropy. For small conductivities, the isofrequency contour of the TM modes is represented by hyperbolas, where the dispersion of the TE modes exhibits a close loop in the centre (orange dashed-dotted line in the inset). In this case, the TE modes are poorly confined near the metasurfaces. The circumstance changes considerably for large conductivities, where the dispersion of the TE modes exhibits hyperbolas (green dashed line in the figure), and the TM modes become poorly confined. In general, the hyperbolic metasurfaces can support a large amount of the overall local density of states compared with the elliptic case, owing to the dominant excitation of either the TM or TE plasmon modes.
To better understand the field excitation features, we plot the spatial distribution and dispersion relation of the SPM using the finite element method (COMSOL Multiphysics, COMSOL Inc., Sweden [69]). The model sets a linear current dipole placed above the TDAM at a distance of L = 10 −3 λ a ; the operation frequency is fixed at 60 THz. As shown in figure 2(a), in the extremely elliptic case, where the isofrequency contour is elongated along the x optical axis, the energy flow of the SPM can be well guided along particular directions at small angles with respect to the y optical axis. In contrast to the isotropic case, for the metasurfaces considered here, the in-plane anisotropy also induces a dominant difference among the field components. This feature is well illustrated in the insets, which show detailed views of the field amplitudes E x and E y near the current source. As indicated, the field amplitudes at the position of the current source are , and it is clear that our numerical results achieve a high accuracy in describing light-matter interactions. Finally, the FFT of the spatial field distribution is given in figure 2(b). The result is in good agreement with the isofrequency contour obtained via the Green function approach, which is depicted by the red solid line in figure 1(a).
Additional to the material anisotropy, the dipole-metasurface distance also influences the anisotropic excitation of the field amplitudes E x and E y . It can be observed in the inset of figure 2(a) that the ratio of the field components continuously decreases with an increase in distance L, which implies a decline in the coupling strength between the dipole and metasurfaces. Physically, the phenomenon can be understood as follows. For small dipole-metasurface distances, the excitation feature of the dipole field is strongly modified by the optical response of the surfaces. Consequently, prominent discrepancies between the in-plane field amplitudes can be produced, owing to strong excitations of the anisotropic SPM. With a further increase in the distance L, the coupling between the dipole and SPM weakens [70], thereby leading to a sharp decrease in the ratio E o x /E o y . This ratio decreases to unity when the distance is sufficiently large, which indicates an isotropic excitation of the field amplitudes. Notably, for systems that consist of QEs interacting with anisotropic metasurfaces, this behaviour demonstrates that the dipole decouples from the surfaces, and the spontaneous decay only occurs via free vacuum modes.
In figures 2(c) and (d), the spatial distributions and FFT of the excited plasmon field with a hyperbolic dispersion are also investigated. Generally, the unidirectional excitation of the SPM on metasurfaces can be observed, and the radiation energy is channelled along specific directions orthogonal to the hyperbola asymptotes. In figure 2(d), the hyperbolic isofrequency contour can be clearly identified through the FFT of the spatial field, which is in good agreement with the numerical result shown in figure 1(c).

Anisotropic Purcell effect and enhanced quantum interference of three-level QEs
Based on the discussions in the previous section, it is clear that TDAMs support the anisotropic excitation of in-plane field amplitudes. In this section, we mainly focus on the quantum interference phenomena of QEs coupled with the TDAM. We consider the case of a three-level QE interacting with the SPM excited on anisotropic metasurfaces, which is illustrated by figure 3. It is probable that dipoles with different in-plane polarizations will display distinct spontaneous decay behaviours, owing to the anisotropic excitations of the SPM. Here, two excited states of the QE are assumed to be close; they are denoted by the states |1 and |2 with energies ω 1 and ω 2 , respectively, and the energy of the ground state is equal to zero and labelled as |g . In our notation, the dipole moment operator of the QE is denoted by l = l σ + 1 e − +σ + 2 e + + h.c., with l defining the magnitude of the dipole, and e − = (e x − ie y )/ √ 2 and e + = (e x + ie y )/ √ 2 refer to the left-and right-rotating unit vectors of the transitionsσ + 1 = |1 g| andσ + 2 = |2 g|, respectively. Under the rotating-wave approximation, the Hamiltonian of the system can be written aŝ In equation (3),f † (r, ω) andf (r, ω) are the bosonic vector field operators for the elementary excitations of the photonic reservoir in the presence of the TDAM [71], ω 1 ≈ ω 2 = ω a represents the resonance frequency between the two degenerate upper levels and the lower level of the three-level QE, with the energy of the lower level considered zero. The QE is placed at the position r a = (0, 0, 10 −3 λ a ), that is, with a distance of 10 −3 λ a above the metasurface. The electric-field vector operator in the frequency domain is as follows [71]:Ê (+) (r, ω) = iωμ 0 dr Ḡ (r, r , ω) ·Ĵ n (r , ω), whereḠ(r, r , ω) is Green's dyadic of the electromagnetic field evaluated at the frequency ω and positions r and r . The noise current operatorĴ n (r , ω) takes the formĴ n (r , ω) = ω ε 0 Im[ε(r , ω)]/πf (r, ω), which acts as the current source of the field and satisfies the following wave equation: where ε(r, ω) and μ(r, ω) are position-and frequency-dependent permittivity and permeability, respectively. At zero temperature, the creation and annihilation field operators satisfy the following commutation relations: Starting from the Hamiltonian given in equation (3) and the quantization scheme described by equations (4)- (7), we transform into the interaction picture and apply second-order perturbation theory. Under the Born-Markovian approximation, the motion equations for the expectation value of the Pauli operators appear to be: [72] d dt In the above equations, γ 1 and γ 2 are spontaneous decay coefficients for the transitions |1 → |g and |2 → |g , respectively, which are modified by the photonic reservoir in the presence of a TDAM. κ 1 and κ 2 represent the QI between two decay channels. Clearly, nonvanished coefficients κ j will produce a QI effect. In this case, the expectation values of the operators σ 11 , σ 22 , and σ 12 couple with each other during the evolution. The spontaneous decay and QI coefficients can be expressed in forms of the Green functions as follows: and where Im[Ḡ xx(yy) (r a , r a , ω a )] represents the imaginary part of the xx (yy) component of the Green dyadic function, Γ x and Γ y denote the spontaneous decay rates of dipole moments with magnitude m and polarizations along two optical axes of the TDAM. It should be noted that in the derivation of equations (11) and (12), the imaginary part ofḠ xy (r a , r a , ω a ) has been neglected because it vanishes in the structure under consideration (i.e. for metasurfaces with vanished transverse conductivity, see appendix A for more details). The relative strength p is defined to measure the degree of QI, with the form [73]  According to equation (13), it is clear that the QI effect disappears in the case of a three-level QE interacting with s free vacuum or isotropic photonic modes, which can be attributed to the vanishing of κ. As we have previously shown, TDAMs support the anisotropic excitations of the SPM with large field confinement. Moreover, the metasurfaces support the anisotropic excitation of the in-plane field amplitudes. Thus, the production and enhancement of the QI effect is expected if the QE is coupled with the metasurfaces. If the photonic environment near the QE becomes highly anisotropic, such that Γ x(y) Γ y(x) is satisfied, the transition channels are effectively parallel to each other, and maximum QI with the relative strength p ≈ ±1 can be achieved.
In figure 4, the Purcell factors (F x(y) = Γ x(y) /Γ 0 , which demonstrates the enhancement of the spontaneous decay rate Γ x(y) in the presence of TDAM; Γ 0 is the decay rate in the free vacuum) of the in-plane dipoles with orthogonal polarizations and the relative strength of the QI are numerically investigated according to equations (11)-(13) and the Green functions given in the appendix. Notably, for the dipole-metasurface distance under consideration, total coupling between the dipoles and SPM can be realized [70]. As illustrated by figures 4(a) and (b), both the x-and y−polarized dipoles display large Purcell factors in the extremely elliptic case (i.e., σ i yy σ i xx or σ i xx σ i yy is satisfied). This can be well understood by combining the results shown in figures 1 and 2. For example, the case of σ i yy σ i xx actually indicates quasielliptic isofrequency contours elongated along the k x axis (red solid curve in figure 1(a)). Under this circumstance, the energy propagation of the supported SPM primarily distributes in small angles with respect to the y axis ( figure 2(a)). As a result, the overall density of states accounting for the spontaneous decay of the x−polarized dipole becomes dominant, which leads to a strong enhancement in the anisotropic Purcell effect (F x F y ). Additionally, sharp declines in the Purcell factors can be observed with an increase in surface conductivities. This behaviour mainly originates from the decrease in the wave numbers of the supported SPM for large conductivities. For small wave numbers, the plasmon field is poorly confined and unable to couple energy from the incoming waves with large wave numbers (green dotted curve in figure 1(a)). In the extreme case of σ i xx , σ i yy σ 0 (σ 0 = e 2 /4 is the conductance quantum), the TM surface modes are poorly confined near the metasurface; thus, coupling between the QE and SPM is negligible. Instead, the free vacuum modes play a key role in determining the spontaneous decay of the dipoles. Based on the above analysis, maximum QI occurs in the extremely elliptic regimes, which is illustrated by figure 4(c).
To clearly identify the influence of surface topology on the spontaneous decay and QI, we also investigate the Purcell factor and the relative strength of the QI in the hyperbolic regime. As demonstrated in figures 5(a) and (b), the spontaneous decay rate increases considerably compared to the elliptic dispersion case. Combined with the results shown in figure 1, it can be deduced that metasurfaces with a hyperbolic topology support the excitation of the SPM with a larger total density of states, owing to the hyperbola feature of the isofrequency contour.  Another remarkable feature of the dipole decay is that a prominent anisotropic Purcell effect can also be achieved in the hyperbolic-dispersion case. This is owing to the progressive shift of the hyperbolic branches towards different optical axes as one of the conductivity components becomes considerably larger than the other. Under this condition, the unidirectional excitation of the SPM results in the selective enhancement of the dipole transitions. With a simultaneous increase in the conductivity components, the isofrequency contour moves away from the origin in the k−space (green dashed curve in figure 1(d)). Consequently, the metasurface cannot support the coupling of incident waves with small wave numbers, which has a destructive effect on the enhanced Purcell factors. Generally, as indicated by figure 5(c), the launching of highly anisotropic SPMs in the extremely hyperbolic case enables the generation of maximum QI between orthogonal decay channels of the three-level QE.
Finally, we consider the dispersive interaction between three-level QEs and anisotropic black phosphorus metasurfaces to emphasise the experimental implementation of our proposed approach. Heterostructures containing black phosphorus or graphene have been previously identified as good candidates for TDAMs with high tuneabilities [12][13][14]. In these structures, the intraband excitations owing to carrier motions can be well described by the Drude model [74], where the interband transitions at higher frequencies are included by introducing the step-absorption function. Based on the above features, the anisotropic conductivities of multilayer black phosphorus can be approximately defined as follows [75]: where the first and second terms on the right-hand side describe the contributions from the intraband and interband transitions, respectively. The effective mass of the electron along the j direction is denoted by m j , τ is the intraband relaxation rate, and j and I j represent the resonant frequencies and relevant strengths of the interband transitions for the jj component of the conductivity tensor, respectively. The carrier concentration of the TDAM, which is denoted by n c in our notation, is experimentally tuneable via bias voltage or chemical doping [76,77].
Before studying the QI property of the proposed system, it is necessary to clarify the dependence of the surface topology on the carrier concentration n c . It is clearly demonstrated in figure 6(a) that both the elliptic and hyperbolic topologies (with abbreviations Ellip. and Hyper.) can be achieved under different conditions. At low carrier concentrations, the optical response of the TDAM is weak (red curves), and fewer electrons participate in the collective oscillation of plasmons. In this case, a relatively low phase-transition frequency between two dispersion topologies can be observed (red dash-dotted line). With an increase in frequency, the interband transition becomes dominant, thereby leading to a sharp negative dip in the imaginary part of the conductivity σ i yy . At large carrier concentrations (blue curves), the optical response of the TDAM becomes considerably stronger. Therefore, a blue-shift of the hyperbolic dispersion regime can be clearly identified.
In figure 6(b), we present a contour plot to illustrate the issue of how a TDAM modifies the relative strength of the QI. Clearly, the strength of the QI can be effectively controlled by manipulating the carrier concentration of the metasurface. It is shown that the critical frequency accounting for the phase change of the surface topology suffers shifts according to the variation of the carrier concentration (red dashed curve). At lower frequencies , the QI strength gradually increases with an increase in the carrier concentration. This behaviour can be well explained by including the results shown in figure 6(a), where the metasurface supports strong excitations of a highly anisotropic plasmon field in the extremely elliptic regime.
Moreover, when the metasurface works in the hyperbolic regime, the dependence of the QI strength on the carrier concentration exhibits rich features. In general, the QI strength reaches both the positive and negative maximums by crossing zero as the excitation frequency varies. This is owing to the dominant contributions of the interband transition at high frequencies, which even forms a dip in the conductivity near the interband transition frequency y (as shown in figure 6(a)). It may seem somewhat counterintuitive that, according to equation (14), an infinite dip occurs at the frequency y . However, this does not occur. For real materials, the interband scattering also contributes to conductivity, which only allows the appearance of the dip with a finite depth [78]. However, the scattering rate only broadens the transition width and does not affect the main physics. Thus, for simplicity, we have omitted it in the approximated form of the conductivity. According to figure 6(a), the optical property of the metasurface experiences dramatic changes as a response to the variation in carrier concentration, ranging from the case of |σ i yy | |σ i xx | to |σ i xx | |σ i yy | in the hyperbolic regime. As discussed previously, in both cases, the metasurface supports strong enhancements of the anisotropic Purcell effects with respect to F x F y and F y F x , which adequately explains the behaviour of the QI by manipulating the carrier concentration.

Conclusions and outlook
This study investigates the excitation of the anisotropic SPM and its influence on the Purcell effect as well as the QI effect between the orthogonal decay channels of a three-level QE in the presence of a TDAM. The Green functions of the electric field for the system under consideration are derived to numerically study the excitation features of the SPM and its influence on the QI effect. The numerical results predict the anisotropic excitation of the in-plane field amplitudes, which is in good agreement with the simulation results. We have extended this remarkable feature to the generation and enhancement of QI. To provide a deeper insight, the anisotropic Purcell effect modified by metasurfaces is also investigated. We have shown that the in-plane transitions can be unequally enhanced, owing to the anisotropic excitation of the SPM. Furthermore, in extremely anisotropic cases, the metasurfaces even support a strong enhancement of the anisotropic Purcell effect, which is the main reason for the appearance of maximum QI. An important aspect of our study consists in showing that the relative strength of QI in three-level QEs coupled with a TDAM can be effectively controlled by engineering the anisotropy of metasurfaces. In real TDAMs, the in-plane anisotropy is determined by the conductivity tensor. As we have presented in section 3, the experimental implementation of the proposed system also demonstrates the feasibility of realizing maximum QI with a flexible manipulation in the carrier concentration of the materials. Our study can be extended to other coupled systems with similar configurations [79,80]; thus, it provides a promising route in the development of highly integrated photonic devices at the nanoscale [81], which may have potential applications in various quantum technologies.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.

Appendix A. Derivation of the Green functions
Consider multi-layered structures constructed by piecewise-constant materials. In the presence of an electric current, the electric and magnetic fields in region n can be generally expressed as [67] E n (r) = k 2 n + ∇∇· Π n (r), ( A . 1 ) where ε n (μ n ) is the relative permittivity (permeability), and k n = ω √ ε n μ n and Π n (r) are the wave number and electric Hertzian potential in the region n, respectively. In this study we consider ε 2 = ε 1 and μ 2 = μ 1 with the assumption that the permittivity and permeability of the host medium are equal to the vacuum values. We assume that the current source exists in the region 1 with time dependence e iωt , and it is denoted by j 1 . Then, in the single-layered case, the Hertzian potentials in the upper (region 1) and lower (region 2) spaces are In the above equations, Π P 1 (r) is the principle part (free propagating modes induced by the electric current without boundary reflections) of the Hertzian potential due to the existence of the current source in the upper space (region 1), and Π S n (r) denotes the scattered potentials in different spatial regions. The principle part of the corresponding Green's dyadic is denoted byḡ P 1 r, r , whereḡ R 1 r, r represents the reflection part that is responsible for the fields reflected by the layered-material,ḡ T 1 r, r is the transmission part owing to the field penetrations through the metasurface scattered into the lower space (region 2), and Ω is the support of the current. In our notation, the principle part of Green's dyadic can be written as where k = k x e x + k y e y is the in-plane wave vector with a magnitude denoted by |k| = k = k 2 x + k 2 y , and its z−component is β 1 = k 2 1 − k 2 . ρ = xe x + ye y is the in-plane component of the position vector r = (x, y, z), which can be varied when evaluating the electric field at different spatial points (similarly, r = (x , y , z ) and ρ denote the position vectors of the current source).Ī is the unit dyadic, the distance between the field and source points (r and r ) is denoted by R = ρ 2 + (z − z ) 2 , with ρ = (x − x ) 2 + y − y 2 representing its projection on the x-y plane. The scattered part of Green's dyadic can be obtained by applying the following boundary conditions: where J s =σ · E is the electric surface current on the boundary. We introduce the two-dimensional Fourier transform on the vector potential: and substitute the Hertzian potential Π n (r) into the boundary conditions given in equations (A.6) and (A.7). After rearranging, we obtain the following equations: In the above equations, Π nj (r) represents the j component of the Hertzian potential in the region n, and the relation k 1 = k 2 is satisfied for the model under consideration. Different from the isotropic case, wherein the z−directed current only induces the potential Π nz (r), in the presence of anisotropic surfaces, the vertical current also induces the potentials Π nx (r) and Π ny (r). Combined with the property that the in-plane currents induce both the horizontal and vertical potentials (i.e., Π nx (r) and Π nz (r) for the x−directed current, Π ny (r) and Π nz (r) for the y−directed current), according to (A.10) and (A.11), the equations for the potentials Π nx (r) and Π nz (r) can be easily determined: Here, the superscripts R and T in the coefficients denote the reflection and transmission parts of the field, respectively. Following the methods mentioned above, the equations for the other scattered coefficients can also be derived. Finally, the integrals related to the scattered Green's dyadic takes the following form: The scattered coefficient g R jl (k) in equation (A.20) clearly contributes to the electric field in the upper space. In the following, we provide detail expressions of the terms g R xx (k), g R yy (k), g R zx (k), and g R zy (k), which are important for obtaining the results of this study: In equations (A.21)-(A.24) the coefficient η = μ 1 /ε 1 , the denominator of the Green functions, defined as D(k) = 2σ xx k 2  In the above expressions, J 0 (kρ) is the zero-order Bessel function,Ḡ P jj r, r andḠ R jj r, r (j = x, y) represent the principle and reflection parts, respectively, where the latter is modified by the metasurfaces and becomes predominant for strong plasmon field excitations. The inner integral functions are Based on the definitions in equations (11) and (12), it is clear that the Green functions in equations (A.25)-(A.27) play a key role in evaluating the spontaneous decay rates of dipoles with different polarizations. Additionally, we notice that in the presence of a TDAM (the boundary conditions are given in equations (A.10) and (A.11)), for a unit strength current source, the x− and y−component of the excited field can be expressed as follows: where the off-diagonal Green functions take the following forms:

Appendix B. Approximation method adopted in the numerical calculations
As indicated in equations (A.25)-(A.32), a numerical evaluation of the Green functions is time-consuming, owing to the two-dimensional feature of the integrals. However, we are primarily interested in the strong QE-SPM coupling cases; therefore, it is feasible to describe the contribution of the SPM by only keeping the residue part of the integrals. Within this approximation, the integrals can be reduced to one-dimensional; thus, a considerable acceleration in the speed of numerical evaluation can be expected. To verify the validity of this method, we present an example that considers the integral defined in equation (A.28), where the dispersion of the SPM is solely determined by the zeros of the denominator in the reflection coefficient g R xx (k).
For specific values of k y , we assume that the roots of D(k) = 0 can be denoted by k sp x . k sp x are functions of k y , according to the definition of D(k). Notably, there are four solutions (accounting for the modes with different wave numbers and propagating features) of the supported wave modes k sp x for every k y ; however, only one of them is allowed under particular conditions. For example, the relation Im k sp x < 0 must be satisfied to obtain a decaying wave that propagates away from the source for the spatial region x − x > 0. In this case, the residue contribution of the SPM on the integral f xx k y (defined in equation (A.28)) takes the following form: In figure B1 we compare the exact numerical result and the residue contribution of the function f xx k y , where the former is obtained by performing the integral defined in equation (A.28), and the latter is obtained by simply evaluating the residue term in equation (B.1) as functions of the in-plane wave number k y . The QE is located at a distance of z a = 10 −3 λ a above the anisotropic metasurface, and different behaviours of the integrand are studied. Subplot (a) demonstrates the imaginary part of the integral function f xx k y versus k y at x = 0. In this case, only the imaginary part of the Green functions evaluated at the spatial position r a = (0, 0, 10 −3 λ a ) contributes to the spontaneous decay of the QE (see equation (A.25)). We note that, except for a subtle difference in the magnitude at small k y and a slight shift in the phase space, the residue approximation is in good agreement with the exact numerical result.
Additional to the source point, the validity of the approximated method under the influence of the propagating SPM is also studied by evaluating the integral function at x = 0.1λ a . Under this condition, the integrand oscillates, owing to the nonvanished factor e −ik sp x (x−x ) . Similarly, as demonstrated in subfigure (c) of figure B1, it is clear that the approximation method possesses high accuracy in evaluating the integral. Within this approximation, we obtain figures 4-6 by numerically evaluating the Green functions. Notably, the approximation method is valid only when the coupling between the QE and SPM makes a dominant contribution to the spontaneous decay [70]. If the QE is too close to the surface, the spontaneous decay process will primarily occur via material dissipation [82]. In another case, where the QE is far from the metasurface, the branch-cut term [83] of equation (A.28) that describes the radiation continuum into space cannot be neglected when performing the integral.