Influence of Surface Roughness on Strong Light-Matter Interaction of a Quantum Emitter-Metallic Nanoparticle System

We investigate the quantum optical properties of strong light-matter interaction between a quantum emitter and a metallic nanoparticle beyond idealized structures with a smooth surface. Based on the local coupling strength and macroscopic Green’s function, we derived an exact quantum optics approach to obtain the field enhancement and light-emission spectrum of a quantum emitter. Numerical simulations show that the surface roughness has a greater effect on the near-field than on the far-field, and slightly increases the vacuum Rabi splitting on average. Further, we verified that the near-field enhancement is mainly determined by the surface features of hot-spot area.

Many studies have focused on realizing strong coupling between a quantum emitter (QE) and nanophotonic environments 7,[32][33][34][35][36] . Nanoparticle-on-mirror (NPoM) 7,37,38 and coupled spheres 39 plasmonic nanostructures with sub-nanometer separation can show strong plasmon-exciton interaction. However, in practice, all fabricated MNPs suffer from nanoscale size inaccuracies and surface roughness that make the distance of QE and nanostructures uncertain [40][41][42][43][44] . Previous studies have shown the differences between the optical response properties of realistic and idealized MNPs [45][46][47][48][49][50] . Specifically, the surface roughness has a slight impact on the far-field response like the scattering cross section 41,46,47,49,[51][52][53] , whereas the near-field enhancement strongly depends on surface features 43,48,49,54 . Therefore, it is reasonable to expect that the imperfect shape of MNPs will modify the light-matter interaction, especially when QEs are very close to MNPs. However, to the best of our knowledge, previous studies have only investigated weak coupling or the far-field response, and the role of surface roughness in the strong coupling regime are still missing in both theory and experiment.
In this work, we show how surface roughness modifies the strong plasmon-QE interaction. For this purpose, a QE is placed just a few nanometers away from an MNP with rough surface. To capture the essential features of surface roughness, we consider a simple model in which a QE coupled to a 20-nm-radius MNP 55 . For detecting the optical response, a far-field spectrum detector is placed 1000 nm away from the MNP, as shown in Fig. 1(a). The distance between the QE and the center of the idealized MNP is 22 nm, and the dielectric function of the MNP is characterized by the Drude mode ω ω ωγ . We choose typical parameters for silver with ε ∞ = 6, ω m = 7.90 eV and γ m = 21 meV 39 . The dipole moment of QE is d = 24D, and the coupled system is embedded in vacuum. Here, the MNP size and the distance between the MNP and the QE (>1 nm) are in which the Drude model very accurately describes the optical response of a metal 56  Δh (see Methods section for more details). In this study, we use the parameters σ = 10 nm and Δh = 1 nm; however, the mean radius of rough MNP is still 20 nm, which agrees with the radius of an idealized MNP. Figure 1(b) shows one realization of rough MNPs with σ = 10 nm and Δh = 1 nm.

Results
Theory. To describe quantum light-matter interaction, we derive a general formalism for obtaining the spontaneous emission properties in a fully quantized approach for both emitters and optical field. We consider the interaction of a two-level QE and MNPs; therefore, we adopt a self-consistent procedure of electromagnetic field in absorbing medium based on the dyadic Green function [57][58][59] . Under this framework, the electric field operator is given by where ω is the angular frequency, c is the speed of light. ε 0 is the vacuum permittivity, and ε and ε I are respectively the permittivity and its imaginary part of the MNP. ˆ † f and f are respectively the creation and annihilation operators of the total electric-field represented by a continuum of harmonic oscillators. The dyadic Green's function G (r, r′, ω) in Equation (1) denotes the field with frequency ω at r responding to a point source located at r′. The general definition of G (r, r′, ω) is the solution of the wave equation ∇ × ∇ × G (r, r′, ω) − (ω/c) 2 ε (r, ω) G (r, r′, ω) = Iδ (r − r′), where I is the unit tensor.
By using Equation (1), the system Hamiltonian under the rotating wave approximation can be written as 0ˆˆˆ † † † where σ = † e g and σ = − g e are the raising and lowering operators of two-level QE, respectively. The QE is located at position r a , with the atom transition frequency ω a and transition dipole moment d = dμ, where μ is unit orientation vector. |e〉 and |g〉 are the upper and lower states of the QE, respectively.
The temporal evolution of system at time t can be obtained by solving the following wave function e g 3 0 where e, 0 represents the QE is in excited state and its surrounding is unexcited, and ω g 1 r , ( , ) indicates that the QE has returned to the ground state and the surrounding optical field is excited. Here, we suppose that QE is initially in the upper state and the optical field is in the vacuum, that is, C e (0) = 1 and C g (r, ω, 0) = 0. Therefore, we obtain the evolution equations of amplitude C e (t) and C g (r, ω, t) as a a 0  where  denotes the principal value. C e (ω) is the local spontaneous emission spectrum (LSES) of QE that corresponds to polarization spectrum 55,62 or dipole spectrum 63 . Γ(r a , ω) and Δ(r a , ω) are the local coupling strength (LCS) 60,61,64,65 and level shift 2,55,66,67 , respectively. The LCS is usually expressed in terms of the local density of states (LDOS) which is defined as 68 a a a 2 ρ (r, ω) is also called the projected LDOS because G(r a , r a , ω) is projected into μ direction, which is, the orientation of the transition dipole moment of QE 55 . In vacuum, the LDOS is ρ 0 = ω 2 /(π 2 c 3 ). Then, LCS can be rewritten as Γ(r a , ω) = πωd 2 /(3ħε 0 )ρ(r a , ω).
LDOS is an essential quantity for describing the emission characteristics of QE from weak to strong coupling regime in micro-and nanostructures 69,70 . Therefore, it is necessary to reveal the influence of surface roughness on LDOS. Moreover, it is widely accepted that the strong coupling is inextricably linked to high field enhancement; therefore, we use the definition of the enhancement of LDOS is the imaginary part of Green's tensor in homogeneous medium. Here ε b is the dielectric constant of homogeneous medium.
The far-field detected spectrum (FFDS) of the emitted light by QE is defined through a double-time integration over the first order quantum correlation function . Therefore, it is straightforward to obtain the analytical expression of FFDS.
is the propagation function from QE to detector. The propagation function describes the essential role of light propagation in FFDS through the Green's propagator G (r d , r a , ω) for two space-point. Equation (10) clearly shows that FFDS contains information about both the local dynamics and the propagation process.
For investigating the influence of surface roughness on Green's propagator G (r d , r a , ω), we define a factor, in a manner similar to Equation (9), called far-field enhancement (FFE).
At the end of this section, we further explain our theoretical method. First, without making Markovian approximation, the LSES and FFDS are valid in both the weak and the strong coupling regime. Second, we can obtain the reversible dynamics of QE in the strong coupling regime 71-73 through the Fourier transformation of LSES C e (ω). Finally, our method can also be extended to multiple QEs system to find out the underlying cooperative effect under strong light-matter interaction.
Based on the above theoretical method, in the following section, we study the differences in optical response between MNPs with smooth surface and rough surface. For evaluating the LCS and propagation function in Equation (10), we choose the MNPBEM toolbox 74-76 based on the boundary element method (BEM) as the numerical simulation platform. The BEM approach only requires discretizing the MNP surface, making it extremely suitable and efficient for plasmonics studies 47,49 . Hereafter, we use a dipole emitter to represent the two-level QE.  Fig. 2(c) is longer than that in Fig. 2(a), indicating that LDOS is more sensitive to surface features when the dipole is perpendicular to the MNP. The same conclusion can also be drawn from far-field enhancement, as indicated by the longer error bar in Fig. 2(d) than in Fig. 2(b). A comparison of the standard error of averaged LDOS and that of far-field enhancement clearly shows that the surface roughness produces a comparatively small effect in the far-field and a dramatic one in the near-field. Furthermore, we find that the LDOS of rough MNPs reduces by ~10% and ~12% on average for x-and z-oriented dipoles compared to idealized MNPs, respectively. Because the surface roughness can shift the resonant peak of LDOS, the averaged LDOS of rough MNPs is reduced and broadened slightly compared to that of idealized MNPs. However, the peak frequency shows no visible change between idealized and rough MNPs in terms of both LDOS and far-field enhancement.

Near-and Far-Field Enhancement.
Local Surface Features Analysis. The above results indicate that the height fluctuation of the surface greatly impacts the near-field enhancement. However, it is difficult to directly determine what type of rough surface can increase the strength of the light-matter interaction. We know that a tip on surface should enhance LDOS owing to the lightning rod effect, whereas the magnitude of enhancement and existence of other effects remain unknown.
In this section, we investigate the optical properties of the MNP with a bump. The bump is placed on the MNP surface at various positions and different heights to change the bump-dipole distance. Considering the dipole orientation and symmetry of the sphere, the bump is sited along both x-and y-directions when using an x-oriented dipole but only along the y-direction when using a z-oriented dipole. Figure 3((A-C), left-hand side) shows a schematic diagram of a bump sited along the x-direction at various positions. Figure 2(b) and (d) indicate that the surface roughness has almost no effect on the far-field, therefore, we focus on the effect of the bump on near-field enhancement. Figure 3(a,b) shows LDOS for MNPs with a small bump for an x-oriented dipole, and Fig. 3(d,e) shows that for a z-oriented dipole. The bump height is 0.5 nm in Fig. 3(a) and (d) and 1.0 nm in Fig. 3(b) and (c). The numerical results clearly show that all LDOSs with bumps are larger than that of the idealized MNP (black dashed line in figures). For the same bump height and direction, the shorter bump-dipole distance, the larger LDOS; for a fixed bump position, LDOS for a 1-nm high bump is larger than that for a 0.5 nm high bump. Therefore, the numerical results confirm that the bump enhances the LDOS, and it shows larger enhancement with shorter bump-dipole distance. In our configuration, the largest the enhancement of LDOS is produced by a bump with 1 nm height and θ = 0.4 rad, whose LDOS is nearly two times that of an idealized MNP. When the bump is far from the dipole, for example, with θ = 0.7 rad (green line in Fig. 3), LDOS of MNPs with a bump converges with that of idealized MNPs. We also note that for a bump with the same height and position, the LDOS for a z-oriented dipole is slightly larger than that for an x-oriented dipole. This agrees with the result in the previous section that the impact of surface roughness is greater when the dipole is perpendicular to MNP. Figure 3(a) and (b) show the peak of LDOS red shifts when a bump is located along the x-direction for an x-oriented dipole; the red-shifting is larger if the bump is closer to the dipole. However, the peak is blue-shifted for a bump located along the y-direction. The blue-shifting changes nonmonotonously with θ and has a maximum of ~0.013 eV at θ = 0.55 rad (blue line in Fig. 3). Similar behavior is also found when using a z-oriented dipole, as shown in Fig. 3(d) and (e). In this case, the maximum blue-shifting is ~0.011 eV. In addition, the bump direction is nonsignificant because the electric field distribution is circularly symmetric (right panel of Fig. 4).
It is worth to point out that a bump on the surface broadens LDOS owing to the appearance of a new plasmon mode at lower energy, and the secondary peak shows larger red-shifting with increasing bump height. This can be seen from a comparison of the red line in Fig. 3(d) with that in Fig. 3(e). For a 1-nm-high bump located along the x-direction with an x-oriented dipole, the secondary peak combined with the main peak of LDOS produces a broad range of enhancement of up to nearly 10 5 (red dashed line in Fig. 3(b)). This feature may enable further exploring the possibility of expanding and controlling the range of enhancement of LDOS. This new plasmonic mode appears because of the bump breaking the surface symmetry. To support this claim, we recover the surface symmetry by replacing the single bump by a ring bump (D, left-hand side of Fig. 3). The LDOS of an MNP with a ring bump shows no obvious broadening compared to that of an idealized MNP for various θ, as shown in Fig.  3

(c) and (f).
Briefly, our calculations show that a bump close to the dipole can increase the LDOS magnitude, shift the resonant frequency slightly, and broaden the LDOS.  The results shown in Figs 2 and 3 indicates that a bump generally takes effect when it is close enough to the dipole. Therefore, the LDOS of rough MNPs may primarily be determined by the surface features of the hot-spot area because of the short distance to the dipole. Here, the hot-spot area is defined as the surface in the range of θ π   0 / 8, where θ is the angle of zenith in spherical coordinates. We choose this range to guarantee that the localized surface plasmon energy is predominately located in this region, as shown in Fig. 4. The origin of the spherical coordinates is the MNP center, and the dipole location is (22 nm, 0, 0).
To verify the above hypothesis, we need to compare the difference in LDOS of two MNPs whose hot-spot areas are identical whereas other surface areas are significantly different. Therefore, we create a partially smooth-rough MNP (smooth-rough MNP, model (III) in Fig. 5(a)) by smoothing the rough surface of the hot-spot area for comparison with an idealized MNP (smooth MNP, model (I)). For contrast, we also generate an idealized MNP with a rough hot-spot area (rough-smooth MNP, model (IV)) and compare it to a completely rough MNP (model (II)). Then, we use these four types of MNPs to form a simulation group; in each group, we use identical realization of the rough surface. Additionally, to avoid the discontinuity of the junction of two areas, we use the interpolation method. Figure 5(a) shows an example of a simulation group.
We generate four simulation groups, G 1 − G 4 , to perform numerical simulations. Figure 5(b) shows the corresponding results. According to the hypothesis, the LDOS of smooth MNP and smooth-rough MNP and that of rough MNP and rough-smooth MNP should share similar features, respectively. The numerical results show the predicted phenomenon: F i of MNPs with the same surface features as the hot-spot area within a group only shows a minor difference. However, a new plasmonic mode may only appear for a completely rough MNP but disappear for a partially rough MNP, like the rough-smooth MNP. However, in most cases, the surface features of hot-spot area are crucial.

Near-and Far-Field Spectra.
In this section, we study the influence of surface roughness on the emission spectra. The LSES and FFDS for rough MNPs are calculated with a specific transition frequency ω a and averaged. We only focus on the influence of surface roughness on the peak frequency, as indicated by the triangle with the horizontal error bar in the emission spectra. Figure 6(a) and (d) show the averaged LSES of rough MNPs for x-and z-oriented dipoles, respectively. For an idealized MNP, of the splitting of LSES are ~92 meV for an x-oriented dipole and ~138 meV for a z-oriented dipole. Numerical simulation results show that the averaged splitting of rough MNPs is larger in both cases, being ~95 meV for an x-oriented dipole and ~143 meV for a z-oriented dipole. The splitting enlarges because both peaks in the averaged LSES shift outward: the peak at lower energy is red-shifted whereas that at higher energy is blue-shifted. We note that the peaks in the LSES also remain in the FFDS; however, light emitted to the far-field adds an additional peak at ~2.77 eV, as shown in Fig. 6(b) and (e). This additional peak is caused by the propagation function  ω r ( , r , ) r d a , whose features are seen in Fig. 2(b) and (d). The error bars of the LSES and FFDS clearly indicate that the peaks in the FFDS are less affected by the surface roughness than the peaks in the LSES, confirming the former conclusion that the surface roughness has less effect on the far-field. Moreover, the deviations of the frequency of the peaks in FFDS between idealized and rough MNPs are too small to observe.
The anticrossing behavior of LSES is the footprint of strong coupling. To examine the anticrossing behavior of LSES, we locate the maximum of peaks in the LSES with various transition frequencies, and we show the results in Fig. 6(c) and (f). We observe a small increase in the vacuum Rabi splitting after introducing the surface roughness for both dipole orientations; however, the increase for the z-oriented dipole is larger. It should be noted that in the anticrossing curves of rough MNPs, greater deviation from idealized MNP emerges at band energy of ~2.96 eV Finally, we investigate how a single bump affects the features of LSES and FFDS. The models are the same as those shown in Fig. 3 and Fig. 7 shows the results of the MNP with a 0.5-nm-high bump. The bump enlarges splitting in both the LSES and the FFDS. Therefore, the existence of a bump is beneficial for achieving strong coupling, especially for a z-oriented dipole and small θ. For an x-oriented dipole, when the bump is located along the x-direction, the splitting in the LSPS and FFDS is larger than the bump located along the y-direction while other conditions remain the same. It is generally more difficult to observe the Rabi splitting in FFDS because of the broadening and asymmetry of the peaks. However, we note that when the bump is located along the y-direction, the splitting is more symmetric compared to that in an idealized MNP, as shown in Fig. 7(e). This situation occurs because the bump enhances the peak at higher energy, which is difficult to observe for an idealized MNP. Furthermore, the bump can narrow the peaks and make the splitting in FFDS easy to observe with a z-oriented dipole. This can be seen by comparing the gray dashed line and solid lines in Fig. 7(f). In addition, Fig. 7 also shows that the lower-energy peak of splitting shows larger shifting than a high-energy peak as θ decreases from 0.7 rad to 0.4 rad. This is because the lower-energy peak is ~2.9 eV and is close to the peak of LDOS, where the variance is larger; therefore, it is susceptible to the surface roughness. By contrast, the peak located at ~2.77 eV in the FFDS shows negligible change. It should be noted that the FFDS for an idealized MNP with a z-oriented dipole has a peak near 2.85 eV, as indicated by the black arrows in Fig. 7(f). This peak is reduced by the bump even with θ = 0.7 rad; it is a feature of the bump that can be observed in the FFDS.

Discussion
Our comparative study of the field enhancement and emission spectra between idealized MNP and rough MNP clearly demonstrates the influence of surface roughness on strong light-matter interaction. In determining the far-field response, an idealized model is adequate even the QE is close to MNP. By contrast, the difference becomes substantial referred to near-field enhancement. The magnitude of the field and its characteristics, such as resonant frequency and plasmonic mode, may change dramatically. In particular, three important features are noticeable. First, the surface roughness can slightly increase the vacuum Rabi splitting on average. Second, the near-field enhancement is strongly dependent on the local surface features of the hot-spot area; however, the random fluctuation of the surface may create a new plasmonic mode and make the LDOS and LSES unpredictable. Finally, the FFDS originates from LSES but is weighted by the propagation function ω r r ( , , ) . The splitting in LSES is in the region of the high-order LSPR mode as it dominates the dynamics of QE that strongly couples with MNP; however, the peaks of the propagation function are predominately located at the dipole LSPR. This mismatch reduces the impact of surface roughness on the far-field and makes the FFDS steady.
In summary, we have studied the influence of surface roughness on near-and far-field responses by simulating plenty of rough MNPs coupled with QE. Both near-and far-field enhancement of QE can be investigated using the BEM method. The LDOS may be greatly influenced by the surface roughness; however, on average, the LDOS and anticrossing are robust to the surface roughness of MNPs. Though the surface roughness of MNPs has obvious effect on LSES, its effect on FFDS is negligible. The results of the current study are obtained at a QE-MNP distance of 2 nm; however, there is no significant difference and limitations in applying it to a ultranarrow gap in composite systems like dimers and NPoM, except for the consideration of nonlocal damping [77][78][79][80][81][82][83] . Our research is of practical importance, and it may be inspiring and helpful in designing experiments for strong plasmon-QE interactions.

Methods
Generation of Gaussian Random Rough Surface. We generate the Gaussian random fluctuation by attaching the Gaussian stochastic phase factor to the Fourier coefficients in the spectrum density, which is defined as where σ 2 is the variance of the height fluctuation, l x and l y are the correlation length along x-and y-directions, respectively, and Δφ is the Gaussian stochastic variable. Then, the rough surface whose height obeys the Gaussian distribution can be obtained by the inverse Fourier transformation x y 1 where −1  denotes the inverse Fourier transformation, Δh is the maximum height variation of the surface roughness. Finally, h(x, y) is projected on the surface of an idealized MNP to obtain a rough MNP.
Data availability. The datasets generated and analysed during the current study are available from the corresponding author on reasonable request. Figure 7. LSES of idealized MNP with a 0.5-nm-high bump located along the x-direction (a) and y-direction (b) for x-oriented dipole and along the y-direction for z-oriented dipole (c). Curves A, B, and C correspond to models A, B, and C in Fig. 3, respectively, whose θ values are 0.4 rad (red line), 0.55 rad (blue line), and 0.7 rad (green line). (d-f) are the corresponding FFDS for models (a-c). Transition frequencies are indicated by black dots and hollow dots on the curves for MNPs with a bump and idealized MNP, respectively.