In this paper, we present the results of numerical simulations of the optical spectra of a three-sphere photonic molecule. The configuration of the system was continuously modified from linear to triangular, in-plane with the fundamental mode excited in one of the spheres and perpendicular to it. We found the relative insensitivity of the spectra to the in-plane deviation from the linear arrangement up to about 110°. For larger angles, the spectra show significant modification consisting of the major spectral peaks splitting and shifting. On the contrary, the spectra are quite sensitive to out-of-plane molecule deviation, even at small angles. Thus, the spectra of photonic molecules can be modified by changing the mutual positions of the constituent resonators, which can be useful in reconfigurable photonic devices.

## I. INTRODUCTION

Optical whispering gallery modes (WGM) are resonances of axisymmetrical optical resonators possessing especially long lifetimes *τ* (or corresponding narrow spectral linewidths *γ* ∼ 1/*τ* in scattering experiments). The electromagnetic field of WGMs is confined to a small spatial region in the vicinity of the surface of the resonator resulting in significant enhancement of the energy density of these modes. WGMs of spherical resonators are characterized by their angular momentum (order) *l*, polar number *m* describing the component of the angular momentum along the polar axis of a corresponding spherical system of coordinates, and radial number *s* specifying the distribution of the field in the radial direction. Scattering resonances of spherical particles are known as Mie resonances.^{1} WGMs correspond to Mie resonances of large orders *l* ≫ 1 and small (order of unity) radial numbers *s* and are characterized by high Q-factors defined as *Q*_{l,s} = *ω*_{l,s}/*γ*_{l,s}, where *ω*_{l,s} and *γ*_{l,s} are frequency and the linewidth of the resonance of order *l* and radial number *s*. In spherical resonators, these parameters do not depend on the polar number *m* resulting in 2*l* + 1 degeneracy of the resonators. The Q-factors of WGMs can reach values of 10^{7} − 10^{9}, which in combination with small mode volumes make WGMs attractive for applications as various types of sensors, filters, nonlinear optical devices, etc. (see reviews in Refs. 2 and 3).

The electromagnetic field of WGMs falls off with the distance from the surface of the resonator in a quasi-exponential manner, similar to the behavior of the wave functions of atom-bound electrons. Even though in the case of WGMs such an evanescent behavior persists only up to a certain distance from the resonator’s center, after which it crosses over to spherical wave-like oscillations, this analogy with atomic wave functions allowed for the system of several WGM resonators positioned in the close proximity of each other to be called photonic molecules.^{4–10} The most research so far has been done on photonic molecules formed by several two-dimensional disk resonators.^{8–10} In a more difficult case of coupled spherical resonators, only diatomic photonic molecules^{4,7,11} and the linear chains of spherical resonators^{6,12} were studied. The reasons for this disparity are quite clear—three-dimensional resonators have a much larger number of modes, which all interact with each other when resonators form a photonic molecule, increasing thereby the computational complexity of the problem. In a diatomic molecule or linear chain, the remaining axial symmetry allows one to avoid having to deal with interactions between modes with different polar numbers *m* by choosing the polar axis along the axis of symmetry. This reduces the number of coupled equations that need to be solved and makes the problem more tractable. Also, from the experimental point of view, it is much easier to assemble multi-resonator planar structures of the disk resonators then build three-dimensional photonic molecules. However, recent experiments^{13} with levitating droplets demonstrated the feasibility of placing three or more almost perfectly spherical droplet resonators in any kind of geometrical configuration making an investigation of the modes of such photonic molecules more timely.

## II. MODEL

In this work, we study normal modes of the system of three identical spherical resonators arranged in configurations ranging from an equilateral triangle formed by the centers of the resonators almost touching each other to a linear chain. Designating positions of the centers of the resonators by their position vectors *r*_{i}, we can present the field in the region outside of the resonators as a linear combination of vector spherical harmonics (VSH) $Nm,l(3)r$ (TM polarization) and $Mm,l(3)r$ (TE polarization),^{1} where upper index three indicates the radial dependence in the form of outgoing spherical waves,

Upper index in the expansion coefficients $al,m(i)$ and $bl,m(i)$ numerates different spheres and runs from 1 to 3. The internal field in each sphere can be presented as

where upper index 1 in the VSHs indicates that the radial dependence is given by the Bessel functions of the first kind, which remain finite at the centers of the spheres. To facilitate applications of the Maxwell boundary conditions at the surface of each sphere, one can use the so-called translation theorem^{14} expressing VSH centered on one sphere in terms of VSHs centered at different spheres. This procedure results in the following system of equations for expansion coefficients of the scattered field:

where $\alpha l(TE)$ and $\alpha l(TM)$ are single sphere Mie coefficients defined as the ratios of the expansion coefficients of the scattered and incident fields, and $Al,ml1,m1rj\u2212ri$, $Bl,ml1,m1rj\u2212ri$ are so-called translation coefficients appearing after application of the translation theorem, which describe radiative coupling between spheres. For convenience of the readers, we provide explicit formulas for the translation coefficients in the Appendix to this article following Ref. 1. We also introduced a TE-polarized external incident field characterized by expansion coefficients $cm,l(inc)$ assuming that it only excites the first sphere (*i* = 1). Once scattering field coefficients $al,m(i)$ and $bl,m(i)$ are found, the expansion coefficients of the internal fields can be immediately computed as

where *n*_{d}, *n*_{w} are indices of refraction of the resonators and surrounding media, *x* = *kR*_{d} is the size parameter, *R*_{d} is the radius of the sphere, *k* is the vacuum wave number of electromagnetic field, and $\u2032$ signifies differentiation with respect to the entire argument of the respective spherical Bessel or Hankel functions.

The complexity of the system of equations defined by Eqs. (3) and (4) is determined by the translation coefficients $Al,ml1,m1rj\u2212ri$ and $Bl,ml1,m1rj\u2212ri$. For instance, when all spheres are arranged along a straight line, these coefficients might be made diagonal in polar indices^{1} $Al,ml1,m1=Al,ml1\delta m,m1,Bl,ml1,m1=Bl,ml1\delta m,m1$ by choosing the polar axis of the spherical coordinate system used to define VSHs along this line. In this coordinate system, Eqs. (3) and (4) with different values of *m* decouple significantly simplifying the numerical analysis of the problem.^{4,6,11,12}

In a non-linear three-resonator photonic molecules, the advantages of this coordinate system are lost, so we introduce a coordinate system with the polar axis and *XY* plane defined in such a way that a mode excited in the sphere 1 in the absence of its interaction with two other spheres would have been a fundamental WGM characterized by *m* = *l*. The *Y*-axis of our system is chosen to always pass through the centers of spheres 1 and 2, while the position of the sphere 3 continuously changes from a linear three-resonator chain (the center of sphere 3 is on the *Y* axis) to an equilateral triangle (the center of sphere 3 is on the *X* axis) and can also move outside of *XY* plane (see Fig. 1). Translation coefficients in Eqs. (3) and (4) depend on the relative positions of the centers of the resonators *r*_{1,2} = *r*_{1} − *r*_{2}, *r*_{1,3} = *r*_{1} − *r*_{3}, and *r*_{2,3} = *r*_{2} − *r*_{3}, which in the coordinate system shown in Fig. 1 are given by their spherical coordinates as

where angles *φ*_{1,3} and *φ*_{2,3} are related to each other via equation *φ*_{2,3} − 2*φ*_{1,3} + *π*/2 = 0 with *φ*_{1,3} changing between *π*/6 for the equilateral configuration and *π*/2 for the linear chain, and angles *ϑ*_{2,3} = 2*ϑ*_{1,3} describe deviation of the center of the 3rd sphere from the aligned equatorial planes of spheres 1 and 2.

Experiments of Ref. 13 were conducted with droplet resonators having radius *R*_{d} ≈ 50 *µ*m and WGMs excited in the spectral range of around *λ* ≈ 700 nm. Taking into account the relative refractive index of the droplets in water *n*_{d} = 1.516, one can estimate the order of the excited WGMs to be around *L* ≈ 2*πn*_{d}*R*_{d}/*λ* ≈ 600. This means that even for a single *l* = *L*, one has to deal with 2*L* + 1 = 1201 coupled equations. Taking into account additional equations with *l* both greater and smaller than *L* will further increase the computational complexity of the problem. The situation can be simplified, however, if one notices that single sphere Mie coefficients $\alpha l(TE,TM)$ for higher order modes have very strong frequency dependence and are very small when calculated at frequencies deviating from their own resonance values by more than the resonance width. This observation indicates that when working in the vicinity of a given frequency one can neglect contributions from all modes whose own resonances are spectrally separated from the frequency range of interest. The extreme case of such an approximation is the so-called resonance approximation in which we keep only coefficients with *l* = *L* and neglect cross-polarization coupling assuming that the incident light is tuned to excite only a particular polarization. For this approximation to work, of course, the free spectral range of a single resonator, which can be estimated as △*λ*_{FSR} ∼ *λ*^{2}/(2*πn*_{d}*R*_{d}), must be much larger than the spectral range *δλ*_{L} occupied by the resonances characterized by a given *L* and different polar numbers *m*. The latter can be estimated using, for instance, a bi-sphere photonic molecule as $\delta \lambda L\u223c(\lambda 2\Gamma L/c)AL,0L(2Rd,0,0)$, where Γ_{L} is the radiative decay rate of a WGM with *l* = *L* and *c* is the vacuum speed of light and the translation coefficients are computed in the coordinate system with its polar axis passing through the centers of the spheres. In this case, polar number *m* remains a conserving quantity and the spectral range occupied by the WGMs split by the interaction is determined by modes with *m* = 0, for which the interaction is the strongest (see, e.g., Ref. 11). If one introduces parameter Γ_{x} = (2*πR*_{d}/*c*)Γ, one can relate parameters △*λ*_{FSR} and *δλ*_{L} as

For experiments of Ref. 13, we can estimate this ratio to be about 0.08, which justifies using the resonance approximation in our calculations. Thus, neglecting all inter-modal and cross-polarization coupling and leaving only expansion coefficients corresponding to the TE resonant term with *l* = *L*, we are left with a much smaller system of equations

Noting that $AL,mL,m1rj\u2212ri=AL,mL,m1ri\u2212rj$ and introducing (2*L* + 1) × (2*L* + 1) matrices

with indices running from − *L* to *L*, we can rewrite the system of Eq. (7) for a three-resonator molecule in the form of three matrix equations

Column vectors *B*^{(i)} (*i* = 1, 2, 3) in Eq. (9) contain 2*L* + 1 expansion coefficients *b*^{(i)} for each of the resonators.

## III. RESULTS

Calculations of translation coefficients with *L* as large as 400 present additional computational challenges due to precision problems while computing Bessel functions of such a high order. While these difficulties are not insurmountable, the purpose of this work is to gain general qualitative understanding of the properties of such structures rather than to try to predict quantitative results of the experiments similar to the one of Ref. 13. Thus, we considered in our simulations WGMs with *L* = 50, *n*_{d} = 1.32, *n*_{w} = 1 and assumed that the excitation of a single sphere would have only excited a fundamental mode with *M* = *L*. As long as results are presented in terms of dimensionless size parameter *x*, the results are independent of *R*_{d} and the resonance wavelength.

We begin by computing spectra of a linear chain over a broad spectral interval containing single sphere resonances of both TE and TM polarizations with *L* ranging between 50 and 53. The results of the computations are presented in Fig. 2. The main feature seen on this graph is the appearance of the repeating patterns of resonances grouped around the main peak whose frequency coincides with the frequency of a single sphere resonance. Each groups of peaks are well separated from each other confirming our assertion that the resonance approximation shall work well in the case under consideration. The spectral pattern consisting of a main peak corresponding to a single sphere resonance surrounded by multiple secondary peaks was first observed in Ref. 12. It was shown in that paper that in a linear chain of spherical resonators with odd numbers of the spheres there always exists a collective mode of the chain with a frequency exactly equal to that of a single sphere resonance. The distribution of the field among the resonators at this frequency according to Ref. 12 is such that only 1st and 3rd spheres are illuminated while the middle sphere remains dark. Distribution of the intensities shown in Fig. 2 reproduces this behavior.

Now we focus on a spectral interval around the single sphere resonance with *L* = 50 and study how the spectral pattern changes when the linear arrangement of the spheres becomes distorted. Assuming that the equatorial planes of all spheres remain aligned (*ϑ*_{2,3} = *ϑ*_{1,3} = 0) when the angle between vectors *r*_{1,2} and *r*_{2,3} changes from zero to 2*π*/3, we generate spectra shown in Fig. 3. The most remarkable feature of these spectra is their surprising stability with respect to the changing alignment. Essentially, all the way through the values of the angles between *r*_{1,2} and *r*_{2,3} close to 110°, the changes in the spectral pattern are insignificant. The only modifications seen in the graphs are slight movement of the peaks at the end of the spectral range toward the central peak and narrowing of the latter. The stability of the spectrum can be explained by two properties of the translation coefficients responsible for coupling between the resonators. First, the translation coefficients rapidly decrease with increasing distance between the centers of the spheres allowing for introduction of the nearest-neighbor approximation used, e.g., in Ref. 12. Second, a fundamental mode (concentrated in the equatorial plane of a resonator and characterized by polar number *M* = *L*) of one sphere interacts most strongly with a fundamental mode of another sphere because the field overlap decreases as polar number *M* changes away from *L* (the field of the modes with *M* < *L* can be thought of as localized around planes inclined with respect to the equatorial plane). Thus, as long as the centers of all spheres remain in the *XY* plane and the interaction between the first and the third spheres can be neglected, the change in the alignment of the centers of the spheres does not affect the position of the main peak—the main spectral feature—while generating only small changes in other spectral features. The most important of these secondary effects is the shift of the peaks at the end of the spectral region toward the center resulting in narrowing of the entire spectral rage corresponding to a given resonance.

Once the angle between vectors *r*_{1,2} and *r*_{2,3} increases beyond 110°, the spectrum starts changing more significantly. It is interesting that these changes take place over very small variations of the angle beyond the critical value, which also can be explained by a sharp dependence of the translation coefficients on the distance between the centers of the spheres. The main reason for the observed spectral modifications is the increasing role of the interaction between first and third spheres as the distance between their centers decreases. This interaction is responsible for the splitting of the center peak into two narrower peaks as well as for the enhancement and splitting of the secondary peaks at the end of the spectral interval associated with *L* = 50 modes. The spectral changes are also accompanied by changes in the distribution of intensities of the field within and between the spheres.

The misalignment of the equatorial planes of the resonators has a more profound effect on the spectra of the photonic molecules. With increasing angle *ϑ*_{2,3}, the height of the central peak decreases, while all other resonances become much more prominent (see Fig. 4). This effect can also be understood by thinking of the field profile in the resonators as linear combinations of VSH defined in a coordinate system with the polar axis perpendicular to the equatorial planes of the first two spheres. When the equators of the all spheres are aligned the strongest interaction is between the fundamental VSH, which results in the strongest resonance at the unperturbed single-sphere frequency. As the third sphere starts deviating from the *XY* plane, VSHs with different polar numbers *m* composing the field in the sphere enjoy better overlaps with VSHs in other spheres resulting in shifting spectral pattern. This trend of diminishing the weight of the single sphere resonance and enhancing the resonances at the ends of the spectral range is observed throughout the entire range of *ϑ*_{2,3} from 0 to 90°. The details of the spectra are quite sensitive to the position of the third sphere: changes of the angle as small as several degrees result in significant movements of the peaks, their overlap and splitting. Apparently, such sensitivity can be also attributed to the strong dependence of the overlap between different vector spherical components comprising the field of different spheres on their mutual positions. In the supplementary material, we present videos showing spectral changes with changing misalignment of the equatorial planes of the individual resonators. The sensitivity of the spectra to the alignment manifests itself via apparent jumps in the video.

## IV. CONCLUSION

In this paper, we presented the results of numerical simulations of the optical spectra of a three-sphere photonic molecule for different configurations of the system. Among the most interesting results of our simulations is the relative insensitivity of the spectra to the deviation of the arrangement from the linear up to angles between lines connecting centers of the spheres of about 110° as long as equatorial planes of the spheres are coincident. The most important change in this situation is the narrowing of the spectral region corresponding to excited mode. We found that misalignment of the equatorial planes has a more profound effect: as the third sphere moves out of the common equatorial plane, we observe the significant modification of the spectra, which consists of spectral peaks merges and splits. The results of our computations show that the spectra of photonic molecules can be modified by changing the mutual positions of the constituent resonators, which can be useful in reconfigurable photonic devices.

## SUPPLEMENTARY MATERIAL

The videos showing the evolution of the spectra of the structure with changing angles between the second and third resonators are shown in the supplementary material.

## ACKNOWLEDGMENTS

This research was supported by the United States–Israel Binational Science Foundation (BSF) (Grant Nos. 2016670 and 2020683), the U.S. National Science Foundation (NSF) (Grant No. 1711801), ICore: the Israeli Excellence Center “Circle of Light” (Grant No. 1802/12), and the Israel Science Foundation (Grant Nos. 537/20 and 1572/15).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Vladimir Shuvayev**: Data curation (equal); Formal analysis (equal); Software (equal); Visualization (equal); Writing – review & editing (equal). **Stanislav Kreps**: Data curation (equal); Validation (equal); Visualization (equal). **Tal Carmon**: Conceptualization (equal); Funding acquisition (equal); Validation (equal); Writing – review & editing (equal). **Lev Deych**: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.

### APPENDIX: TRANSLATION COEFFICIENTS

Translation coefficients are defined as in Ref. 1,

where

and

are Wigner 3*j* symbols.