Definition of Polaritonic Fluctuations in Natural Hyperbolic Media:   Applications to Hexagonal Boron Nitride, Bismuth Selenide and the Spontaneous   Emission Sum Rule

Sean Molesky; Zubin Jacob

arXiv:1702.01862·physics.optics·March 27, 2019

Definition of Polaritonic Fluctuations in Natural Hyperbolic Media: Applications to Hexagonal Boron Nitride, Bismuth Selenide and the Spontaneous Emission Sum Rule

Sean Molesky, Zubin Jacob

PDF

TL;DR

This paper develops a framework to quantify electromagnetic fluctuations in hyperbolic media like hexagonal boron nitride and bismuth selenide, revealing significant field enhancements and implications for light-matter interactions.

Contribution

It introduces a novel approach linking hyperbolic response traits to field coupling, enabling fluctuation calculations without nonlocality assumptions.

Findings

01

Electromagnetic field fluctuations can be over 120 times larger than in vacuum.

02

The optical sum rule is satisfied but does not limit emission enhancement.

03

Broad spectral regions exhibit significantly increased electromagnetic fluctuations.

Abstract

The discovery of photonic hyperbolic dispersion surfaces in certain van der Waals bonded solids, such as hexagonal boron nitride and bismuth selenide (a topological insulator), offers intriguing possibilities for creating strongly modified light-matter interactions. However, open problems exist in quantifying electromagnetic field fluctuations in these media, complicating typical approaches for modeling photonic characteristics. Here, we address this issue by linking the identifying traits of hyperbolic response to a coupling between longitudinal and transverse fields that can not occur in isotropic media. This description allows us to calculate the influence of hyperbolic response on electromagnetic fluctuations without explicitly imposing a characteristic size (model of nonlocality), leading to formally bounded expressions so long as material absorption is included. We then apply this…

Equations151

⟨ E (r, ω) \otimes E^{*} (r^{'}, ω) ⟩ = \frac{ω Θ ( ω , T )}{π} Im {\overset{ˇ}{G} (r, r^{'}, ω)},

⟨ E (r, ω) \otimes E^{*} (r^{'}, ω) ⟩ = \frac{ω Θ ( ω , T )}{π} Im {\overset{ˇ}{G} (r, r^{'}, ω)},

F (ω) = Tr [Im {\overset{ˇ}{G} (r, r^{'}, ω)}] = \int d V_{k} Tr [Im {e^{i k \cdot (r - r^{'})} \overset{ˇ}{G} (k, ω)}],

F (ω) = Tr [Im {\overset{ˇ}{G} (r, r^{'}, ω)}] = \int d V_{k} Tr [Im {e^{i k \cdot (r - r^{'})} \overset{ˇ}{G} (k, ω)}],

w_{_{L}}

w_{_{L}}

w_{_{T}}

w = w_{_{L}} + w_{_{T}},

w = w_{_{L}} + w_{_{T}},

k \times w_{_{L}}

k \times w_{_{L}}

k \cdot w_{_{T}}

k \cdot E_{_{L}} (k, ω)

k \cdot E_{_{L}} (k, ω)

B_{_{L}} (k, ω)

k \times E_{_{T}} (k, ω)

i c^{2} k \times B_{_{T}} (k, ω)

iω E_{_{L}} (k, ω)

k \cdot D_{_{L}} (k, ω)

k \cdot D_{_{L}} (k, ω)

\overset{ˉ}{B}_{_{L}} (k, ω)

k \times \overset{ˉ}{E}_{_{T}} (k, ω)

i c^{2} k \times \overset{ˉ}{B}_{_{T}} (k, ω)

iω D_{_{L}} (k, ω)

E_{_{L}} (k, ω)

E_{_{L}} (k, ω)

B_{_{L}} (k, ω)

\overset{ˉ}{E}_{_{L}} (k, ω)

\overset{ˉ}{E}_{_{L}} (k, ω)

\overset{ˉ}{B}_{_{L}} (k, ω)

ϵ (k, ω) = 0

ϵ (k, ω) = 0

\overset{ρ}{ˉ} (k, ω) = i ϵ_{_{0}} k \cdot \overset{ˉ}{E}_{_{L}} (k, ω),

\overset{ρ}{ˉ} (k, ω) = i ϵ_{_{0}} k \cdot \overset{ˉ}{E}_{_{L}} (k, ω),

\overset{ˉ}{E}_{_{T}} (k, ω)

\overset{ˉ}{E}_{_{T}} (k, ω)

\overset{ˉ}{E}_{_{L}} (k, ω)

\overset{ˉ}{B}_{_{T}} (k, ω)

ω (I_{3} k^{2} - k \otimes k - k_{_{0}}^{2} \overset{ϵ}{ˇ} (k, ω)) \overset{ˉ}{A}_{_{T}} (k, ω) + k_{_{0}}^{2} \overset{ϵ}{ˇ} (k, ω) k \overset{ˉ}{V} (k, ω) = 0,

ω (I_{3} k^{2} - k \otimes k - k_{_{0}}^{2} \overset{ϵ}{ˇ} (k, ω)) \overset{ˉ}{A}_{_{T}} (k, ω) + k_{_{0}}^{2} \overset{ϵ}{ˇ} (k, ω) k \overset{ˉ}{V} (k, ω) = 0,

k \overset{ϵ}{ˇ} (k, ω) (ω \overset{ˉ}{A}_{_{T}} (k, ω) - k \overset{ˉ}{V} (k, ω)) = 0,

k \overset{ϵ}{ˇ} (k, ω) (ω \overset{ˉ}{A}_{_{T}} (k, ω) - k \overset{ˉ}{V} (k, ω)) = 0,

k = ϵ_{_{P}} (k, ω) k_{_{0}},

k = ϵ_{_{P}} (k, ω) k_{_{0}},

k \overset{ϵ}{ˇ} (k, ω) k = 0,

k \overset{ϵ}{ˇ} (k, ω) k = 0,

ϵ_{_{U}} (k, θ, ω) = 0,

ϵ_{_{U}} (k, θ, ω) = 0,

ϵ_{_{U}} (k, θ, ω) = s (θ)^{2} ϵ_{_{P}} (k, ω) + c (θ)^{2} ϵ_{_{A}} (k, ω) .

ϵ_{_{U}} (k, θ, ω) = s (θ)^{2} ϵ_{_{P}} (k, ω) + c (θ)^{2} ϵ_{_{A}} (k, ω) .

\overset{ˉ}{V} (k, ω) = ω k \overset{ϵ}{ˇ} (k, ω) \overset{ˉ}{A}_{_{T}} (k, ω) / (k \overset{ϵ}{ˇ} (k, ω) k)

\overset{ˉ}{V} (k, ω) = ω k \overset{ϵ}{ˇ} (k, ω) \overset{ˉ}{A}_{_{T}} (k, ω) / (k \overset{ϵ}{ˇ} (k, ω) k)

(\overset{ˇ}{I} \frac{k ^{2}}{k _{_{0}}^{2}} - \frac{k \otimes k}{k _{_{0}}^{2}} - \overset{ϵ}{ˇ} (k, ω) + \frac{ϵ ˇ ( k , ω ) k \otimes ϵ ˇ ( k , ω ) k}{k ϵ ˇ ( k , ω ) k}) \overset{ˉ}{A}_{_{T}} =

(\overset{ˇ}{I} \frac{k ^{2}}{k _{_{0}}^{2}} - \frac{k \otimes k}{k _{_{0}}^{2}} - \overset{ϵ}{ˇ} (k, ω) + \frac{ϵ ˇ ( k , ω ) k \otimes ϵ ˇ ( k , ω ) k}{k ϵ ˇ ( k , ω ) k}) \overset{ˉ}{A}_{_{T}} =

0.

k = ϵ_{_{E}} (k, θ, ω) k_{_{0}},

k = ϵ_{_{E}} (k, θ, ω) k_{_{0}},

ϵ_{_{E}} (k, θ, ω) = \frac{ϵ _{_{A}} ( k , ω ) ϵ _{_{P}} ( k , ω )}{ϵ _{_{U}} ( k , θ , ω )} .

ϵ_{_{E}} (k, θ, ω) = \frac{ϵ _{_{A}} ( k , ω ) ϵ _{_{P}} ( k , ω )}{ϵ _{_{U}} ( k , θ , ω )} .

\overset{ˉ}{E}_{_{L}} (k, ω) = \hat{k} ϵ_{_{H}} (k, θ, ω) \overset{ˉ}{E}_{_{T}} (k, ω),

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Definition of Polaritonic Fluctuations in Natural Hyperbolic Media: Applications to Hexagonal Boron Nitride, Bismuth Selenide and the Spontaneous Emission Sum Rule

Sean Molesky

Department of Electrical Engineering, Princeton University, Princeton, NJ 08544, USA

[email protected]

Zubin Jacob

School of Electrical and Computer Engineering, Brick Nanotechnology Center, Purdue University, West Lafayette, IN 47907, USA

Department of Electrical and Computer Engineering, University of Alberta, Edmonton, AB T6G 2V4, CAN

Abstract

The discovery of photonic hyperbolic dispersion surfaces in certain van der Waals bonded solids, such as hexagonal boron nitride and bismuth selenide (a topological insulator), offers intriguing possibilities for creating strongly modified light-matter interactions. However, open problems exist in quantifying electromagnetic field fluctuations in these media, complicating typical approaches for modeling photonic characteristics. Here, we address this issue by linking the identifying traits of hyperbolic response to a coupling between longitudinal and transverse fields that can not occur in isotropic media. This description allows us to calculate the influence of hyperbolic response on electromagnetic fluctuations without explicitly imposing a characteristic size (model of nonlocality), leading to formally bounded expressions so long as material absorption is included. We then apply this framework to two exemplary areas: the optical sum rule for modified spontaneous emission enhancement in a general uniaxial medium, and thermal electromagnetic field fluctuations in hexagonal boron nitride and bismuth selenide. We find that while the sum rule is satisfied, it does not constrain the enhancement of light-matter interactions in either case. We also show that both hexagonal boron nitride and bismuth selenide possess broad spectral regions where the magnitude of electromagnetic field fluctuations are over 120 times larger, and over 800 times larger along specific angular directions, than they are in vacuum.

In 1987, Yablonovitch conceived the photonic crystalYablonovitch (1987); John (1987) as a means of physically forbidding electromagnetic field fluctuations over a finite bandwidth to improve the performance of semiconductor lasers and solar cells. Variations of this idea, tracing back as far as Purcell’s pioneering work on nuclear magnetic resonancePurcell (1946), appear ubiquitously in contemporary optics. Manipulation of a system’s field fluctuations characteristics (the photonic density of states, two-point correlations, ect.) provides a means of controlling an extensive list of phenomena, including field enhancementLi et al. (2017); Khurgin and Sun (2017), sub-wavelength confinementArcari et al. (2014); Men et al. (2014); Goban et al. (2015), thermal propertiesGuo et al. (2012a); Yan et al. (2017); Jin et al. (2017), spontaneousHoang et al. (2015); Pelton (2015); Pick et al. (2017) and coherent emissionNowack et al. (2007); Wu et al. (2015a), and atom-atom interactionsDong et al. (2017); Cortes and Jacob (2017). This strong relation between fluctuational properties and system response has naturally led to reciprocal forms of the problem originally proposed by YablonovitchMunday et al. (2012); Belacel et al. (2013); Raman et al. (2013); Sauvan et al. (2013); Lin et al. (2016). Namely, to what extent can electromagnetic field fluctuations inside a system be enhanced over a given spectral bandwidth? Remarkably, realizable complements of the photonic bandgap, media in principal offering unbounded enhancement, do exist. Hyperbolic (indefiniteSmith and Schurig (2003)) media are widely stated to posses a broadband photonic dispersion singularityJacob et al. (2010, 2012) leading to electromagnetic field fluctuations of magnitude bounded only by second order effects.

Yet, while this picture of singular (indefinite) field fluctuations has been highly successful for interpreting functional applications of hyperbolic media spanning the domains of imagingJacob et al. (2006); Liu et al. (2007); Guo et al. (2012b), nanophotonicsYao et al. (2011); Ginzburg et al. (2013); Milton et al. (2013); Galfsky et al. (2015); Kim et al. (2015); Iorsh et al. (2015), and certain semi-classicalJacob and Shalaev (2011); Newman et al. (2013); Shekhar and Jacob (2014); Slobozhanyuk et al. (2015) and thermal interactionsMolesky et al. (2013); Jacob (2014); Gu et al. (2014); Li et al. (2015); Dyachenko et al. (2016); Dai et al. (2017), there are a number of situations of interest where it can not be easily applied. In particular, an unbounded magnitude of electromagnetic fluctuations presents difficulties for connecting semi-classical and quantum optical processes with existing results in macroscopic quantum electrodynamicsMilonni (1995); Barnett and Loudon (1996); Tip (1997) when the fields inside a hyperbolic media must be described directly. (As opposed to situations where boundary conditions can be used to reformulate all quantities of interest in terms of external fields.) With the growing realization that many natural materials intrinsically exhibit high quality hyperbolic responseEsslinger et al. (2014); Dai et al. (2014); Caldwell et al. (2014); Yoxall et al. (2015); Wu et al. (2015b); Korzeb et al. (2015); Nemilentsau et al. (2016), there are at least two motivations for reexamining this issue. First, there is interest in exploring whether hyperbolic response could be useful in proposed technologies relying on quantum optical effectsCortes et al. (2012); Jacob (2012); Kapitanova et al. (2014); Ferrari et al. (2015), which would benefit from a simplified description. Second, as all hyperbolic media transition smoothly between frequency windows of hyperbolic and (normal) elliptic or isotropic response, the lack of a single, transparently consistent, formulation to simultaneously treat both settings is bothersome.

Significant steps have been made towards this goal since the start of the decade. In 2011, Maslovski and SilveirinhaMaslovski and Silveirinha (2011) (wire metamaterial) and Poddubny et al.Poddubny et al. (2011) (general hyperbolic media) showed that the observable effects of the fluctuation singularity are smoothed once the emitter (atom) is given a finite size. Like the more familiar longitudinal fluctuation divergences encountered in isotropic media, it was found that considering field interactions with a spatially finite charge distribution leads to bounds proportional to $1/a^{3}$ and $1/a$ specific to the hyperbolic regime, with $a$ denoting the characteristic size. (For later reference, it is worthwhile to note that the particular form of the distribution was observed to have a weak influence on these terms.) This was shortly followed by the work of Potemkin et al.Potemkin et al. (2012), formally rederiving the functional structure of the Green function in real spaceChen (1983); Felsen and Marcuvitz (1994), in principle setting hyperbolic and isotropic media on equal footing. The singular terms found in the expressions of this work ( $\delta$ distribution terms) corrected earlier simplified guesses at the Green function singularity for uniaxial media given by WeiglhoferWeiglhofer (1989). Nearly concurrently, studies by Poddubny et al.Poddubny et al. (2012) and Yan et al.Yan et al. (2012) explored the influence of microscopic models. Considering a cubic lattice of hyperbolically polarizable dipoles, and the interlaced metal dielectric multilayer structures frequently used for creating hyperbolic metamaterials respectively, both investigations found that nonlocal response strongly altered results based on homogeneous approximations. Of particular relevance to this article, the later study showed that nonlocality in the polarization response produces another mechanism bounding fluctuations in hyperbolic media, leading to a proportionality of $1/\beta^{3}$ , with $\beta$ representing the length scale of nonlocality. (This form matches what we find for the purely longitudinal contribution, but misses the polaritonic characteristics. Additional discussion of these developments are given in the review by Poddubny et al.Poddubny et al. (2013). Note that the language used in these studies, which is conventional, denotes what we will refer to as the magnitude of fluctuations or fluctuation density, as the Purcell factor or photonic density of states. It is our belief, at least in the context of this article, that this alternate terminology will help avoid potential confusion.)

Nevertheless, while the above findings do in fact contain all the ingredients required to characterize a general uniaxial medium, they remain unsatisfactory in practice. Most importantly, there is no directly calculable quantity that qualifies the singular nature of the fluctuation density in the hyperbolic case, nor any obvious connection with the widely successful intuitive understanding of these media. This information is unequivocally contained in the real space Green function as it has been previously presented, but it is not easily accessible. (Recovering the results we will present here from real space Green function is not simple, and requires conscientious treatment of the formally singular parts of this expression.) The two finite-size normalization approaches face similar difficulties. While technically accurate, fluctuations in hyperbolic media appear to become inextricably tangled with the characteristic length scale. This is both a clear break from isotropic media, where (as we will review) the transverse fluctuation density is directly determined from the permittivity parameters, and a cumbersome impediment for performing calculations in natural hyperbolic media, as in most cases only local permittivity data is readily available.

Here, we confront these shortcomings, consider in the plasma physics community as early as 1962Arbel and Felsen (1962); Fisher and Gould (1969), through the essential coupling between transverse and longitudinal electric fields that arises in any anisotropic medium. This routinely overlooked connection allows us to produce a characteristic regularized field fluctuation density bounded by material absorption, even in the hyperbolic case; and to analytically link the polariton excitationsIshii et al. (2013) known to occur in hyperbolic media with the near-field optical and thermal properties that they exhibit. We then apply this framework to study two areas of interest. The first of these is the optical sum rule for modified spontaneous emission enhancement. This general result states that the frequency integrated transverse fluctuation density is a constant of any photonic system. (For instance, in photonic crystals the suppression of fluctuations in the band gap is compensated by a corresponding enhancement at the band edge van Hove singularities.) Based on previously reported descriptions of the uniaxial Green function, it is not clear why this well established results continues to hold in the presence of hyperbolic response. The second application is a calculation of the thermal fluctuation density (total electromagnetic energy density) for hexagonal boron nitride and bismuth selenide. In these investigations, we find that while the sum rule is valid, it does not extend to the polariton modes that define hyperbolic response. We also observe that both hexagonal boron nitride and bismuth selenide have broad spectral regions where the fluctuation density is over 120 times larger, along specific directions over 800 times larger, than it is in vacuum.

From the fluctuation dissipation theoremKubo (1966); Rytov et al. (1988),

[TABLE]

and our stated objective amounts to regularizing the fluctuation density (FD) :

[TABLE]

for hyperbolic media. (Where $\Theta\left(\omega,T\right)$ is the energy of a harmonic oscillator at frequency $\omega$ and temperature $T$ , $\check{\mathcal{G}}\left({\bf r},{\bf r}^{\prime},\omega\right)$ the dyadic Green function of the medium, $\check{\mathcal{G}}\left({\bf k},\omega\right)$ it’s Fourier transform, $\int dV_{\bf k}$ an integral over reciprocal space and $\text{Im}\left\{...\right\}$ the imaginary part.) The text is organized into five sections. The first three cover our theoretical work leading to equations (27), (34), (39) and (40), along with brief reviews of relevant background information. The last two explore our chosen applications.

I Polariton Excitations in Anisotropic Media

We begin by decomposing the Maxwell equations in terms of a chosen direction in reciprocal space vector ${\bf k}$ . Letting

[TABLE]

be the projection of a vector ${\bf w}$ along $\hat{\bf k}$ , and onto the plane perpendicular to $\hat{\bf k}$ respectively, any vector ${\bf w}$ can be represented as

[TABLE]

where ${\bf w}_{{}_{L}}$ and ${\bf w}_{{}_{T}}$ are referred to as the longitudinal and transverse components. From these definitions

[TABLE]

so that the Maxwell equations in vacuum separate to become

[TABLE]

with $\epsilon_{{}_{0}}$ and $\mu_{{}_{0}}$ denoting the permittivity and permeability of vacuum, $\rho\left({\bf k},\omega\right)$ the charge density, ${\bf j}\left({\bf k},\omega\right)$ the current density, ${\bf B}\left({\bf k},\omega\right)$ the magnetic field, and ${\bf E}\left({\bf k},\omega\right)$ the electric field. Assuming that the relative permeability is negligibly different than vacuum, $\check{\mu}\left({\bf k},\omega\right)=\check{I}$ , as we will throughout, macroscopic averaging of (5a)-(5d) produces

[TABLE]

where ${\bf D}\left({\bf k},\omega\right)=\check{\epsilon}\left(k,\omega\right){\bf E}\left({\bf k},\omega\right)$ is the electric displacement field, and the $f$ subscript is introduced as a shorthand to mark that the quantity is free, i.e. separate from the microscopic densities that have been averaged over in producing (6a)-(6d) from (5a)-(5d). The overline $\bar{X}$ serves as a similar reminder that the electric and magnetic fields appearing in (6a)-(6d) are spatially averaged, and not equivalent to the identically named fields in (5a)-(5d).

I.1 Microscopic densities

The microscopic equations, (5a) and (5b), show that the longitudinal electric and magnetic fields are entirely determined by their respective charge densities,

[TABLE]

so that they vanish in the absence of charge. (Since the most important properties related to these quantities are the Coulomb self-energy and charge momentumCohen-Tannoudji et al. (1997), we interchangeably refer to them as Coulombic.) A homogeneous solution to (6a)-(6d) is hence purely transverse, and from (6c), determined by ${\bf E}_{{}_{T}}\left({\bf k},\omega\right)$ . (We refer to all such quantities that define homogeneous solutions as normal variables.)

I.2 Isotropic media

A treatment of isotropic media follows largely by analogy. From the relation between ${\bf D}\left({\bf k},\omega\right)$ and $\bar{\bf E}\left({\bf k},\omega\right)$ , defined by the scalar relative permittivity $\check{\epsilon}\left(k,\omega\right)=\epsilon\left(k,\omega\right)I_{3}$ , the longitudinal electric and magnetic fields are

[TABLE]

Therefore, so long as $\epsilon\left(k,\omega\right)\neq 0$ , the normal variables of an isotropic medium are again transverse and equivalent to $\bar{\bf E}_{{}_{T}}\left({\bf k},\omega\right)$ .

The caveat to this congruence is the appearance of the polarization condition

[TABLE]

in equation (6a). When this condition is met, the displacement field may be zero even if the longitudinal electric field is not. As the remaining macroscopic equations do not depend on the longitudinal electric field, these Coulombic modes evolve independent of the transverse electromagnetic solutionsFerrell (1958); Berreman (1963); Newman et al. (2015). Averaging (5a) directly

[TABLE]

showing that each Coulombic solution is a mechanical macroscopic oscillation of the microscopic charge density, mediated by the electric field.

The appearance of these Coulombic solutions make (6a)-(6d) fundamentally different than a scaled vacuum. The fact that $\epsilon\left({\bf k},\omega\right)$ exists because of the presence of charges is inescapable, even after macroscopic averaging. However, the effects resulting from the two solution types can usually be treated independently as they tend to exhibit markedly different behaviour.

I.3 Anisotropic (uniaxial) media

For anisotropic media the relative permittivity tensor, $\check{\epsilon}\left(k,\omega\right)$ , can not simplified and a more careful analysis is required. Rewriting (6a)-(6d) in the Coulomb gauge,

[TABLE]

a homogeneous solution requires both

[TABLE]

and

[TABLE]

with $k_{{}_{0}}=\omega/c$ , $\bar{{\bf A}}\left({\bf k},\omega\right)$ standing for the electromagnetic vector potential, and $\bar{V}\left({\bf k},\omega\right)$ the scalar potential. In order to satisfy (13) there are three distinct possibilities.

(O): If in addition to being perpendicular to ${\bf k}$ , $\bar{{\bf A}}_{{}_{T}}\left({\bf k},\omega\right)$ is constrained to directions perpendicular to ${\bf k}\;\check{\epsilon}\left(k,\omega\right)$ then $\bar{V}\left({\bf k},\omega\right)=0$ . As in the microscopic picture, the normal variables are then transverse. Evaluating (12), simplifying to a uniaxial response as we will throughout, produces the s-polarized, or ordinary, wave condition

[TABLE]

with $\bar{{\bf A}}_{{}_{T}}\left({\bf k},\omega\right)$ confined to the direction $\hat{\bf s}=\left[-s\left(\phi\right),c\left(\phi\right),0\right]$ relative to the unit direction in reciprocal space $\hat{\bf k}=\left[s\left(\theta\right)c\left(\phi\right),s\left(\theta\right)c\left(\phi\right),c\left(\theta\right)\right]$ . (Our labeling convention for uniaxial media is shown in Fig.1.)

(C): If ${\bf k}\;\check{\epsilon}\left(k,\omega\right){\bf k}=0$ then $\bar{V}\left({\bf k},\omega\right)$ can be non-zero independent of the value of $\bar{{\bf A}}_{{}_{T}}\left({\bf k},\omega\right)$ . These purely longitudinal solutions generalize the Coulombic modes of an isotropic media, (9), with the updated criterion

[TABLE]

accounting for the reduced symmetry of $\check{\epsilon}\left(k,\omega\right)$ . For uniaxial anisotropy, (15) reduces to

[TABLE]

with

[TABLE]

We will refer to this projection as the uniaxial permittivity of the medium.

(AP): If $\bar{{\bf A}}_{{}_{T}}\left({\bf k},\omega\right)$ is not perpendicular to ${\bf k}\;\check{\epsilon}\left(k,\omega\right)$ , then (13) forces the equality

[TABLE]

and (12) becomes

[TABLE]

Given the directional restrictions on $\bar{{\bf A}}_{{}_{T}}$ , satisfaction of this constraint requires $k$ to be a solution of the p-polarized, or extraordinary, wave condition

[TABLE]

with $\bar{{\bf A}}_{{}_{T}}\left({\bf k},\omega\right)$ confined to the direction $\hat{\bf p}=$ $\left[-c\left(\theta\right)c\left(\phi\right),-c\left(\theta\right)s\left(\phi\right),s\left(\theta\right)\right]$ , and $\epsilon_{{}_{E}}\left(\bf k,\theta,\omega\right)$ defined as the extraordinary permittivity

[TABLE]

(The similarity between (21) and the excitation condition of a surface plasmon polaritonFord and Weber (1984) is not coincidental.)

Substitution into $\bar{V}\left({\bf k},\omega\right)=\omega{\bf k}\;\check{\epsilon}\left(k,\omega\right)\bar{{\bf A}}_{{}_{T}}\left({\bf k},\omega\right)/$ $\left({\bf k}\;\check{\epsilon}\left(k,\omega\right){\bf k}\right)$ shows that for the extraordinary family of solutions

[TABLE]

where $\bar{{E}}_{{}_{T}}\left({\bf k},\omega\right)$ is the undetermined scalar magnitude of the transverse component of the electric field,

[TABLE]

is the relative degree of polarization anisotropy between the optical axis and plane, and

[TABLE]

defines the hyperbolic permittivity.

(22) contains the essential physics that will guide the rest of the manuscript: in an anisotropic media the normal variables are a mixture of transverse and longitudinal fields. Averaging (5a) as before,

[TABLE]

making it is apparent that for any extraordinary solution the electromagnetic (transverse) oscillation is accompanied by a Coulombic charge oscillation. In light of this fundamental coupling, we refer to these solutions as anisotropic polaritons (AP). (Notice that the same analysis can be applied to the magnetic field and relative permeability tensor $\check{\mu}\left(k,\omega\right)$ .) From B., it is clear that in isotropic media such excitations are impossible. The global direction of the electric field for a Coulombic mode is uniquely fixed by the propagation direction of the charge oscillation. This means that longitudinal electric fields can not couple to magnetic fields, and hence can not be electromagnetic. (24) and (25) show that the Coulombic part of an AP type mode grows proportionally with the degree of anisotropy of the medium, (23), and is resonant with zeros of the uniaxial permittivity, (17). Accordingly, these properties also characterize the distinguishing features of hyperbolic response.

II Normal Variable Decomposition of the Anisotropic Green Function

Substituting (6c) into (6d), the electric field inside an anisotropic medium obeys the equation

[TABLE]

The dyadic Green function of a uniaxial medium is the formal inverse of this relation,

[TABLE]

where all reciprocal vectors have been normalized by $k_{{}_{0}}$ , and, recalling our previous definitions,

[TABLE]

In isotropic media, $\epsilon_{{}_{H}}\left(k,\omega\right)$ reduces to zero while $\epsilon_{{}_{U}}\left(k,\omega\right)$ and $\epsilon_{{}_{E}}\left(k,\omega\right)$ become the isotropic permittivity $\epsilon\left(k,\omega\right)$ so that (27) simplifies to

[TABLE]

(For the remainder of the article, the U and I superscripts mark that the results applies specifically to either uniaxial or isotropic media.)

As an operator, the Green function (27) determines the electric field generated by a point current source as a modal expansion of the three homogeneous solution families.

(O): Ordinary electromagnetic

[TABLE]

(C): Coulombic

[TABLE]

(AP): Anisotropic polariton

[TABLE]

Taking the trace of the imaginary part of the uniaxial Green function, the FD in reciprocal space is then

[TABLE]

(The trace of the coupling matrices $\hat{\bf p}\otimes\hat{\bf k}$ and $\hat{\bf k}\otimes\hat{\bf p}$ are zero.) The poles of this function show that the first term again represents ordinary (O) excitations, while the final two terms form the combined contributions of the anisotropic polariton (AP) and purely longitudinal (C) families, which couple due the presence of shared longitudinal fields. Regardless of this mixing, the influence of these two types of excitations are distinguished by their poles. The transverse $\;\hat{\bf p}\otimes\hat{\bf p}$ term has $|k^{2}-\epsilon_{{}_{E}}\left(k,\omega\right)|^{2}$ (AP) type poles, but not $|\epsilon_{{}_{U}}\left(k,\omega\right)|^{2}$ (C) type poles.

III Model-Independent Characteristic

To this point, all permittivity factors have been written as functions of the magnitudes $k$ and $\omega$ . From the inverse Fourier transformation, these dependencies correspond to spatial and temporal nonlocality. For most problems in macroscopic electromagnetic, the spatial dependence is dropped, producing

[TABLE]

For the FD, this approximation raises difficulties. Ignoring nonlocality removes all spatial dependence of the Coulombic solutions, making the associated longitudinal field scale invariant. There is then no suppression of arbitrarily high momentum contributions, and the Coulombic piece of the FD diverges. (This issue is of course not unique to macroscopic electromagnetics, and is commonly dealt with in quantum field theories by the introduction of regulators, parameter renormalization and counter termsSchwartz (2014).)

Yet, while problematic for properties like the excitation lifetime of embedded emittersBarnett et al. (1999), for isotropic media the local approximation is nevertheless largely valid; and many relevant geometry independent characteristics can be determined from the homogeneous Green function without any treatment of nonlocality. This follows from a careful consideration of the length scales involved. For optical and infrared energies, the wavelength of light is almost always large compared to the scale of nonlocality (e.g. the scale of the material lattice), and the difference between $\epsilon\left(k,\omega\right)$ and $\epsilon\left(\omega\right)$ is usually insubstantial for any externally excited electric field. (A fact that also makes accurate probing of the permittivity response at optical and infrared frequencies above $k/k_{{}_{0}}\approx 5$ Shekhar et al. (2017) difficult.) This approximate indistinguishability of $\epsilon\left(k,\omega\right)$ and $\epsilon\left(\omega\right)$ is equivalent to the assumption that all Coulombic poles can be moved to arbitrarily large reciprocal vector without tangibly altering the Green function for all pertinent reciprocal vectors. Such a change has has two primary effects. First, the propagation length of Coulombic solutions becomes vanishingly small. Second, interaction of an external field with Coulombic solutions becomes possible only in the presence of a physical discontinuity. Combining these two features, longitudinal fields can only exist at the interface of two media (surface charge densities), and therefore, only influence surface effects. Volume characteristics, related to the propagation of an external field once it enters the mediumCohen-Tannoudji et al. (1997), are conversely contained only in the electromagnetic modes, and hence require no knowledge of the Coulombic component.

The mixed fields appearing in AP type modes would seem to preclude a similar separation of domains in anisotropic media. The asymptotic behavior of (34) shows that the final term of

[TABLE]

diverges as $k^{3}$ in the limit of local permittivity response. Since this term contains contributions from the AP type modes, it is implausible that a meaningful characterization of hyperbolic response can be captured without it. In turn, this would mean that the FD of a hyperbolic media (or anisotropic media generally) is only describable once a microscopic model (or other physically motivated normalization) is specified.

To work around this apparent difficulty, we begin by expanding all absolute values and imaginary parts of (34), treating $k$ as a real variable, and extending the resulting expression over the entire complex $k$ plane. Jordan’s lemma and the Cauchy integral theorem then imply that if $k^{2}\mathcal{F}^{{}_{U}}\left({\bf k},\omega\right)$ tends to a constant value as $|k|\to\infty$ the value of $\mathcal{F}^{{}_{U}}\left(\omega\right)$ equals the residues of $\mathcal{F}^{{}_{U}}\left({\bf k},\omega\right)$ . (Here, we are taking infinite semi-circle contours in the upper and lower half spaces depending on the value of ${\bf k}\cdot{\bf r}$ , with the convention that the square root function is cut along the positive real axis.) Assuming nonlocality leads to $\text{Im}\left\{\epsilon_{{}_{A}}\left(k,\omega\right)\right\}$ and $\text{Im}\left\{\epsilon_{{}_{P}}\left(k,\omega\right)\right\}$ having scaling at least $\propto 1/k^{2}$ , so that there is no contribution from the path at arbitrarily large $k$ , these residues split into two distinct classes.

(i) Poles of the form $k_{{}_{S}}=\pm\sqrt{\epsilon_{{}_{X}}\left(k_{{}_{S}},\omega\right)}$ resulting from the wave equations that have a second order dependence on spatial nonlocality, i.e. poles that tend to $k_{{}_{S}}=\pm\sqrt{\epsilon_{{}_{X}}\left(\omega\right)}$ in the limit of local response. (Here $k_{{}_{S}}$ and $\epsilon_{{}_{X}}$ are place holder labels that could apply to either ordinary or AP modes.)

(ii) Poles that have a first order dependence on the exact characteristics of nonlocality, and thus tend to towards arbitrarily large values in the limit of local response.

Following the same reasoning as isotropic media, volume characteristics, for low to moderate $k$ fields, must depend only on the first class pole. Therefore, we can conclude that the sum of these residues is the correct generalization of the transverse FD of an isotropic material to an anisotropic setting.

The validity of the above argument rests on the asymptotic $k$ scaling of the imaginary part of permittivity response being stronger than $\propto 1/k^{2}$ , but this is a generally valid assumption for any medium. As $\epsilon\left(k,\omega\right)-1=\chi\left(k,\omega\right)$ is a response function, it is analytic for all but an finite set of points in the complex $k$ plane for a given $\omega$ Abrikosov et al. (1975). Correspondingly, there is a convergent Laurent series expansion in complex $k$ such that

[TABLE]

for $M<|k|<\infty$ , where $M$ is magnitude of the largest $\bf k$ pole of $\chi\left(k,\omega\right)$ , and $\left\{c_{n}\left(\omega\right)\right\}$ are the frequency dependent coefficients of the expansion. From the definition of the Fourier transform, $\chi\left(k,\omega\right)$ is related to the real space susceptibility as

[TABLE]

and the required convergence of this expression for all values of $r$ is only guaranteed if $\left(n\geq 1\right)\Rightarrow c_{n}\left(\omega\right)=0$ . (Although we have assumed inversion symmetry, $\epsilon\left({\bf k},\omega\right)=\epsilon\left(k,\omega\right)$ , throughout the manuscript all arguments and forms we present hold generally.)

With this result in hand, the residues determined by the wave equation poles sum to give

[TABLE]

for O type modes, and

[TABLE]

for AP type modes, with $k_{{}_{O}}$ and $k_{{}_{E}}\left(\theta\right)$ standing for the solutions to (14) and (20) nearest those of the local approximation in the complex plane. In (40), the first term results from the transverse field, and the second from the longitudinal field.

Convincingly, (39) and (40) are also precisely the result obtained by considering the normal variables found in the second section. Directly, (27) identifies the transverse part of the AP excitations with the second term of (34). Since there is no question as to the convergence of this term, the existence of (40) (the longitudinal part) must follow as a requirement of the Maxwell equations. (Either approach to (40) is independent of the specific form of nonlocality considered.) Further, the permittivity dependence of the Coulombic piece of this term

[TABLE]

is the same as that observed in calculating the power radiated by a dipole in a losses hyperbolic mediumFelsen and Marcuvitz (1994); Zhukovsky et al. (2013), and is essentially an angular version of the FD associated with a surface plasmon polariton excitationFord and Weber (1984). (Note that in this contribution the usual distinctions between one and two-sheeted hyperboloidsCortes et al. (2012) - metallic response along the optical axis, or in the optical plane - are nearly absent.)

Computation of the FD resulting from the second class of poles, on the other hand, does requires a specific model of nonlocalityHorsley and Philbin (2014); Pollard et al. (2009) (just as in isotropic media). To focus our discussion we will not pursue these details. Still, there are some general characteristics worth noting. Considering the longitudinal component of (34), once $k^{2}$ surpasses $|Re\left\{\epsilon_{{}_{E}}\left(k,\theta,\omega\right\}\right)|$ the prefactors $\left(|k^{2}-\epsilon_{{}_{P}}\left(k,\omega\right)|/|k^{2}-\epsilon_{{}_{E}}\left(k,\theta,\omega\right)|\right)^{2}$ and $\left(|k^{2}-\epsilon_{{}_{A}}\left(k,\omega\right)|/|k^{2}-\epsilon_{{}_{E}}\left(k,\theta,\omega\right)|\right)^{2}$ will both quickly approach unity. (Assuming nonlocality does not drastically increase the peak magnitude of the polarization response for real $k$ .) Once this condition is achieved, the final term of (34) is increasingly well approximated as

[TABLE]

This is again the exact result found by considering the normal variables of the Coulombic solutions independently, and a straightforward extension of the Coulombic FD encountered in isotropic media,

[TABLE]

For $k$ where the above approximation is valid, the residues from these (C) type poles can safely be attributed to pure Coulombic modes.

Running contrary to this discussion, it should be emphasized that in cases where strong interactions with Coulombic solutions are expected (for example embedded emitter) that the Green function forms provided previously, in particular (34), are well suited to computation. They are valid in any general uniaxial (or isotropic) medium, and can be normalized using any specific model of non-locality or standard quantum field theory approach.

IV The Sum Rule for Modified Spontaneous Emission Enhancement in Hyperbolic Media

The sum rule for modified spontaneous emission enhancement, formulated by Barnett and LoudonBarnett and Loudon (1996, 1998), states that it is not possible to alter the total relative rate of spontaneous emission into purely electromagnetic (transverse) modes. If the properties of a medium enhance the relative rate of spontaneous emission in one spectral range, they must equally suppress this relative rate in another. Mathematically, this is written as

[TABLE]

where

[TABLE]

is the relative rate of spontaneous emission of a single level emitter of frequency $\omega$ , with transition dipole moment ${\bf d}$ , at position ${\bf r}$ in a medium described by $\check{G}\left({\bf r},{\bf r}^{\prime},\omega\right)$ , $\pi_{{}_{T}}$ is the transverse projection operator, and

[TABLE]

is the rate of spontaneous emission in vacuum.

As anisotropy does not change the frequency dependence of the Green function, an argument for the general validity of the sum rule, given by ScheelScheel (2008), theoretically guarantees that the transverse part of the uniaxial Green function

[TABLE]

must satisfy (44) so long as the permittivities considered satisfy the Kramers-Kronig relations. An illustrative example of this result, assuming local Lorentzian polarization responses

[TABLE]

for $\epsilon_{{}_{A}}\left(\omega\right)$ and $\epsilon_{{}_{P}}\left(\omega\right)$ with $\omega_{\rho}=\left\{500_{{}_{A}},\;700_{{}_{P}}\right\}\;cm^{-1}$ , $\omega_{{}_{0}}=\left\{600_{{}_{A}}\;,1000_{{}_{P}}\right\}\;cm^{-1}$ and $\gamma=\left\{5_{{}_{A}}\;,10_{{}_{P}}\right\}\;cm^{-1}$ is provided in Fig.2. From the graph, we observe that the regions of hyperbolic response are essentially featureless, and, as in isotropic media, enhancement follows the magnitude of polarization. In fact, Fig.2(D) shows that the orientationally averaged enhancement of this transverse part is nearly equivalent to considering the planar and axial permittivities separately and summing the result. (Precisely, replacing (47) with the sum of $(2/3)\;\mathcal{F}^{I}_{{}_{T}}\left(\omega\right)$ with $\epsilon\left(k,\omega\right)=\epsilon_{{}_{P}}\left(\omega\right)$ and $(1/3)\;\mathcal{F}^{I}_{{}_{T}}\left(\omega\right)$ with $\epsilon\left(k,\omega\right)=\epsilon_{{}_{A}}\left(\omega\right)$ , produces a result that is not significantly different form that of the true hyperbolic medium.)

These results for the transverse enhancement of the spontaneous emission, which hold to arbitrarily low absorption ( $\gamma$ ), follow from the normal variable picture. Any property of a linear macroscopic medium should be consistent with some arrangement of dipoles in a vacuum. Since scattering from a dipole does not introduce new electromagnetic modes, there is no way that the sum can be modified. From (40) the resonant effects of hyperbolic response for AP type excitations occur in the Coulombic fieldCortes and Jacob (2017). In taking the strictly electromagnetic (transverse) part of (40) these features are ignored.

Orientationally averaged spontaneous emission enhancement resulting from the Coulombic portion of the AP FD

[TABLE]

is plotted in Fig.3. (For numerical convenience in the remainder of this article the FD will be taken to be vacuum normalized by the prefactor appearing in (49).) Comparing with Fig.2, it is clear that this enhancement does not obey the sum rule: it is an additional positive contribution that grows arbitrarily large as material absorption is decreased. By itself, this fact is not particularly unusual. In a general isotropic medium the absorption of energy into matter is not limited by the number of electromagnetic modes (44), and so neither is the enhancement contribution of Coulombic modes. Yet, there are key distinctions between these two cases.

(i) The AP enhancement of the FD does not diverge in the limit of local permittivity response (nonlocality is a second order effect). This is not the case for purely Coulombic enhancementHorsley and Philbin (2014).

(ii) The AP enhancement of the FD is not simply related to the magnitude of the polarization, or total density of charge carriers, as has been shown for Coulombic enhancement in isotropic mediaScheel (2008). Instead, it depends principally on the magnitude of anisotropy, $|\epsilon_{\Delta}\left(\omega,\theta\right)|^{2}$ and material absorption, $\left(Re\left\{\epsilon_{{}_{U}}\left(\omega,\theta\right)\right\}|\epsilon_{{}_{U}}\left(\omega,\theta\right)|^{2}\right)^{-1}$ .

(iii) The AP type enhancement of the FD shows a unique scaling with material absorption, which is stronger than the scaling exhibited by either the transverse O type (47) or longitudinal C type (42) enhancement, Fig.4.

These final two points are salient for potential applications involving embedded quantum emitters. Although nothing can be stated unequivocally without knowledge of the nonlocal response features, the scaling trends indicate that polaritonic solutions become dominate if either the material loss, $Im\left\{\epsilon\right\}$ , is small, or the anisotropy is large.

V Thermal Fluctuations in Hexagonal Boron Nitride and Bismuth Selenide

Like the degree of relative spontaneous emission enhancement, the thermal energy density in the electric and magnetic fields is similarly set by the FD through the relation

[TABLE]

with

[TABLE]

Using (1), (39), (40), and (6c), the electromagnetic energy density in the O and AP type modes of a uniaxial medium is then

[TABLE]

with

[TABLE]

and

[TABLE]

denoting the relative electric and magnetic contributions, and $\mathcal{U}_{BB}\left(\omega,T\right)$ the energy density of an ideal blackbody.

The results of this expression (ignoring pure Coulombic contributions) for hexagonal boron nitride and bismuth selenide, normalized by $\mathcal{U}_{BB}\left(\omega,T\right)/2$ for direct comparison with the FD, are plotted in Fig.5. Four aspects of this figure warrant attention.

(i) The solid green curves, denoting the contribution of AP type modes, confirm that in real media either a high degree of anisotropy $\epsilon_{{}_{\Delta}}\left(\theta,\omega\right)$ or low material absorption $\epsilon_{{}_{U}}\left(\theta,\omega\right)$ may lead to a large polaritonic FD. Bisumth selenide exhibits substantial material absorption, yet nevertheless, large enhancement results from the extreme difference between the axial and planar permittivity components. Hexagonal boron nitride possess much less anisotropy, but has lower material absorption, leading to a nearly identical FD.

(ii) The energy density of the magnetic field is often roughly an order of magnitude larger than that of the electric field (total lines).

(iii) (51) equates the green lines with the FD contributions of AP (bold line) and O (dashed line) type modes. Both hexagonal boron nitride and bismuth selenide show broad spectral regions where the FD is over 120 times larger than vacuum.

(iv) Moving to Fig.6, the AP component of the FD has extreme angular dispersion, concentrated along the critical angles determined by $Re\left\{\epsilon_{{}_{U}}\left(\theta,\omega\right)\right\}=0$ . Along the cone set by this angle, the polaritonic FD is over 800 times larger than vacuum FD in both bismuth selenide and hexagonal boron nitride.

Given the degree of enhancement observed, one may question whether the use of local response in (53) and (54) is accurate. Without experimental evidence, this is open; but we are inclined to believe that the approximation does hold. In either material, the largest absolute value of the permittivities $\epsilon_{{}_{A}}\left(\omega\right)$ , $\epsilon_{{}_{P}}\left(\omega\right)$ , and $\epsilon_{{}_{E}}\left(\omega\right)$ is $\approx 400$ , and the smallest $\approx 0.1$ . Based on these bounds, in the local approximation, all wave equation poles occur below $20\;k_{o}$ . Taking the largest lattice spacing present in either material, $\approx 3\;nm$ , this upper limit of $k$ still corresponds to less than $1\%$ of the Brillouin zone for wavelengths longer than $6\;\mu m$ . (The smallest wavelength of hyperbolic response considered occurs in hexagonal boron nitride at $6.25\;\mu m$ .) As such a small change will only minimally modify the probed bandstructure around the dominant optical phonon featuresSerrano et al. (2007); Zhang et al. (2010); Wu et al. (2015b), substantial variation of the permittivity response should not be expected.

It is interesting to compare the electromagnetic energy plotted in Fig.5 with the full near-field energy density above a half space of hyperbolic mediaGuo et al. (2012a); Guo and Jacob (2014). For this purpose, equation (16) of Guo et al.Guo and Jacob (2014) is plotted in Fig.7. (Other examples of calculated thermal properties that can be compared with volume FDs appear in a recent work by Liu et al.Liu and Xuan (2016).) At the nearest observation points considered, the near-field calculation produces values larger than the associated bulk value, signaling that enhancement at these distances is driven primarily by the excitation of surface charges. All the same, recalling that in Fig.5 no C type or surface polariton modes are included, the two figures show good agreement. The additional peaks seen in Fig.7 match the surface polariton condition $Re\left\{\epsilon\left(k,\omega\right)\right\}=-1$ . This observation indicates that (39) and (40) are the correct measures of local electromagnetic fluctuation characteristics in a hyperbolic medium, akin to the index of an isotropic medium. Further support of this claim is seen in experimentally reported confinement factors for hexagonal boron nitride resonatorsCaldwell et al. (2014), which are within $30~{}\%$ of the FD we have found for hexagonal boron nitride at these wavelengths.

VI Summary

In summary, we have shown that material absorption analytically quantifies electromagnetic fluctuations in hyperbolic (generally uniaxial) media in a manner entirely analogous to the isotropic media. From this result, we have studied the sum rule for modified spontaneous emission enhancement, and have found that it does not apply to the key polaritonic features of a hyperbolic response. We have also investigated the density of electromagnetic fluctuations (electromagnetic thermal energy density) inside both hexagonal boron nitride and bismuth selenide. We have found that both media have broad spectral regions where this quantity is over 120 times (along specific angular directions 800 times) larger than it is in vacuum. Our results unify the computation of electric field fluctuations in uniaxial and isotropic settings, and should prove useful for testing the potential of hyperbolic systems in emerging optical technologies.

Acknowledgments

This work has been supported by funding from the National Science and Engineering Research Council of Canada, Alberta Innovates Technology Futures, and the National Science Foundation. We thank Cris Cortes and Prashant Shekhar for valuable discussions in preparing the article, as well as Chensheng Gong and Todd Van Mechelen for their contributions to the initial investigation of the ideas we have pursued here.

Bibliography91

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Yablonovitch (1987) E. Yablonovitch, Phys. Rev. Lett. 58 , 2059 (1987).
2John (1987) S. John, Phys. Rev. Lett. 58 , 2486 (1987).
3Purcell (1946) E. M. Purcell, Phys. Rev. 69 , 681 (1946).
4Li et al. (2017) K. Li, N. J. Hogan, M. J. Kale, N. J. Halas, P. Nordlander, and P. Christopher, Nano Lett. 17 , 3710 (2017).
5Khurgin and Sun (2017) J. B. Khurgin and G. Sun, in Quan. Plas. (Springer, 2017) pp. 303–322.
6Arcari et al. (2014) M. Arcari, I. Söllner, A. Javadi, S. L. Hansen, S. Mahmoodian, J. Liu, H. Thyrrestrup, E. H. Lee, J. D. Song, S. Stobbe, et al. , Phys. Rev. Lett. 113 , 093603 (2014).
7Men et al. (2014) H. Men, K. Y. Lee, R. M. Freund, J. Peraire, and S. G. Johnson, Opt. Express 22 , 22632 (2014).
8Goban et al. (2015) A. Goban, C.-L. Hung, J. Hood, S.-P. Yu, J. Muniz, O. Painter, and H. Kimble, Phys. Rev. Lett. 115 , 063601 (2015).