Local density of optical states in the three-dimensional band gap of a ﬁnite photonic crystal

A three-dimensional (3D) photonic band gap crystal is an ideal tool to completely inhibit the local density of optical states (LDOS) at every position in the crystal throughout the band gap. This notion, however, pertains to ideal inﬁnite crystals, whereas any real crystal device is necessarily ﬁnite. This raises the question as to how the LDOS in the gap depends on the position and orientation inside a ﬁnite-size crystal. Therefore, we employ rigorous numerical calculations using ﬁnite-difference time domain simulations of 3D silicon inverse woodpile crystals ﬁlled with air or with toluene, as previously studied in experiments. We ﬁnd that the LDOS versus position decreases exponentially into the bulk of the crystal. From the dependence on dipole orientation, we infer that the characteristic LDOS decay length (cid:2) ρ is mostly related to far-ﬁeld dipolar radiation effects, whereas the prefactor is mostly related to near-ﬁeld dipolar effects. The LDOS decay length has a remarkably similar magnitude to the Bragg length for directional transport, which suggests that the LDOS in the crystal is dominated by vacuum states that tunnel from the closest interface toward the position of interest. Our work leads to design rules for applications of 3D photonic band gaps in emission control and lighting, quantum information processing, and in photovoltaics.

(Received 26 February 2020; accepted 1 May 2020; published 8 June 2020) A three-dimensional (3D) photonic band gap crystal is an ideal tool to completely inhibit the local density of optical states (LDOS) at every position in the crystal throughout the band gap. This notion, however, pertains to ideal infinite crystals, whereas any real crystal device is necessarily finite. This raises the question as to how the LDOS in the gap depends on the position and orientation inside a finite-size crystal. Therefore, we employ rigorous numerical calculations using finite-difference time domain simulations of 3D silicon inverse woodpile crystals filled with air or with toluene, as previously studied in experiments. We find that the LDOS versus position decreases exponentially into the bulk of the crystal. From the dependence on dipole orientation, we infer that the characteristic LDOS decay length ρ is mostly related to far-field dipolar radiation effects, whereas the prefactor is mostly related to near-field dipolar effects. The LDOS decay length has a remarkably similar magnitude to the Bragg length for directional transport, which suggests that the LDOS in the crystal is dominated by vacuum states that tunnel from the closest interface toward the position of interest. Our work leads to design rules for applications of 3D photonic band gaps in emission control and lighting, quantum information processing, and in photovoltaics.

I. INTRODUCTION
Controlling the properties of matter by means of quantum light lies at the heart of quantum optics and cavity quantum electrodynamics (cQED). A prime example is the control of the radiative rate of elementary emitters such as atoms, ions, molecules, or quantum dots. Such control is essential for myriad applications ranging from miniature lasers and lightemitting diodes [1,2], via single-photon sources for quantum information processing [3], to solar energy harvesting [4]. To explore such new applications, a suitably tailored dielectric environment is required wherein the vacuum fluctuations, which play a central role in spontaneous emission [5,6], are controlled. Much after the early realization by Purcell [7] that an emitter's environment such as a cavity controls the emission rate, spontaneous emission control has become one of the main drivers of the burgeoning field of nanophotonics [8][9][10][11][12]. Following the seminal predictions by Bykov and by Yablonovitch, emission control was first studied on photonic crystals [1,13]. Emission control has also successfully been pursued with many different nanophotonic systems and * mavidis@iesl.forth.gr † Present address: ASML Netherlands B.V., 5504 DR Veldhoven, Netherlands. ‡ w.l.vos@utwente.nl many different quantum emitters, for instance, atoms and dye molecules in Fabry-Pérot microcavities [14,15], quantum dots in pillar microcavities [16,17], ions in whispering-gallerymode microspheres [18][19][20], dye molecules in plasmonic nanocavities and on nanoantennae [21][22][23][24][25], or dye in metamaterials [26,27].
In the weak-coupling approximation in cQED that is also known as the Wigner-Weisskopf approximation [28], spontaneous emission of an excited quantum emitter is precisely described by Fermi's golden rule [29] wherein the radiative decay rate is linearly proportional to the local density of optical states (LDOS). The LDOS counts the available number of electromagnetic modes each weighted by their strength at each point r 0 and the projection of their electric field along the axes x, y, z [30][31][32]. The LDOS depends sensitively on the close environment of the emitter. Interestingly, the LDOS not only controls spontaneous emission and blackbody radiation, but also plays a role in van der Waals and Casimir dispersion forces and in Förster resonant energy transfer between different emitters [33]. Since the LDOS represents the density of vacuum fluctuations, it controls the amount of vacuum noise experienced by a qubit [34].
From theory, it is well known that the LDOS is radically inhibited at frequencies within the three-dimensional (3D) band gap in an infinite 3D photonic crystal [31,[35][36][37][38][39][40][41]. The LDOS vanishes at any position in the unit cell, and thus throughout the whole crystal, as well as for all dipole orientations. Concerning photonic crystal experiments, the first studies were reported on 3D crystals without 3D band gap [42][43][44][45][46][47][48][49][50][51][52][53][54], or on band gap crystals with low-efficiency emitters [55,56], and in parallel there was a theoretical study on the anomalous Lamb shift [57]. Leistikow et al. studied efficient quantum dots in inverse woodpile photonic band gap crystals and observed exponential time-resolved decay, typical of weak coupling [58]. A 10× inhibited spontaneous emission rate was observed inside the band gap. Since the emission was averaged over many emitters, it was inferred that a single quantum dot at the center of the crystal would be up to 160× inhibited. To date, however, these results have not been interpreted by theory or numerical calculations.
It is obvious that experimental studies and devices employ finite photonic band gap crystals as energy can radiate from the boundaries of the finite crystal. Consequently, states from the infinite surrounding vacuum tunnel into the crystal, 1 leading to a nonzero LDOS and DOS inside the band gap [59,60]. Therefore, it is natural to wonder how the LDOS in the gap depends on the position and orientation of the emitter inside the crystal. For two-dimensional (2D) photonic crystals, Asatryan et al. found in numerical calculations that the LDOS decreases exponentially from the surface into the crystal [61]. Hermann and Hess found a strong position and orientation dependence of spontaneous emission within the unit cell of an inverse opal and saw that the inhibition in the band gap is on the order of two magnitudes, even for relatively small crystals [62]. Kole reported an exponentially growing inhibition at the center of a spherical inverse opal photonic band gap crystal [63]. Leistikow et al. proposed that the LDOS averaged over a unit cell decreases exponentially with position for frequencies inside the 3D band gap, with a characteristic length scale, the so-called LDOS decay length, but no prediction was offered for the dependence within the unit cell [58].
Thus, it appears that calculations of the 3D LDOS in a 3D photonic band gap crystal are scarce in the literature, due to their extensive computational cost and complexity. Therefore, we systematically investigate in this work the position and orientation dependent inhibition of the LDOS in the band gap of 3D inverse woodpile crystals with finite support. Despite recent progress on analytical approaches in nanophotonics [64,65], there are to date no known analytic solutions for realistic 3D crystals; hence we have embarked on a numerical study to address the questions above. We study the role of the position and interpret the computational results by an analytical expression for the expected behavior of the LDOS. We also study the role of the dipole orientation, and compare it to theoretically known behavior [66]. Since we decided to investigate the experimental results of Leistikow et al. [58], we have chosen to study the inverse woodpile crystal structure that was originally proposed by Ho et al. [67]. In our study, we find remarkable physical features, namely, that (1) the LDOS decreases exponentially with position in the crystal, (2) the magnitude of the exponential length scale, the LDOS decay length ρ , is mostly determined by far-field radiation effects whereas the amplitude prefactor is mostly determined by near-field effects, and (3) the magnitude of the LDOS decay length ρ is remarkably close to the Bragg length-that typifies directional transport [68][69][70]-which implies that the LDOS is strikingly directional.

A. The structure of the finite crystal
The inverse woodpile photonic crystal has a primitive unit cell that is illustrated in Fig. 1. The crystal structure consists of two orthogonal 2D arrays of identical cylindrical pores with radius r p = 0.24a running parallel to the x and z axes [67]. 2 The lattice constants are a (in the y direction) and c (in the x and z directions) in a ratio a/c = √ 2 for the crystal structure to be cubic with a diamond-like symmetry. We discuss the LDOS as a function of the reduced frequencyω that is defined asω ≡ ωa/(2π c 0 ) with c 0 the speed of light in vacuum. The backbone of the crystal has the dielectric constant ε b = 12.1, typical of silicon in the near-infrared and telecom spectral ranges. The cylindrical pores are considered to be either empty (ε p = 1) or filled with a dielectric with ε p = 2.25 that is typical for liquids such as toluene that are used to suspend quantum dot emitters in experiments; see Ref. [58]. In the experimentally relevant spectral range, silicon and toluene are essentially lossless. The finite crystals have an extent of N unit cells along each of the x, y, and z axes with a total volume of V = N 3 unit cells. Figure 2(a) shows the band structure of the infinite crystal with empty pores calculated using the plane wave expansion method [74]. The shaded area in Fig. 2(a) indicates the 3D photonic band gap with a broad relative bandwidth ω/ω mid = 25.0% centered atω mid = 0.585, in good agreement with earlier work [75][76][77]. Figure 2(b) shows the band structure for the crystal filled with toluene. Due to the decreased dielectric contrast the 3D photonic band gap has a reduced relative width ω/ω mid = 6.4%. The band gap is centered at a lower frequency nearω mid = 0.49 due to the increased effective average dielectric constant [78]. It is seen that in both air-and toluene-crystal cases the Y stopgap is larger than the X and the Z. This is sensible since in the y direction one encounters a dielectric contrast of two sets of pores whereas along the x and y directions one encounters dielectric contrast for only one set of pores.

B. Computation of the local density of states
It is well known that the LDOS ρ (i) (ω, r 0 ) at a point r 0 projected along the i axis (i = {x, y, z}) is proportional to the total power P (i) (ω, r 0 ) radiated by an electric point dipole current source J(ω, r 0 ) = −iωp(ω)δ(r − r 0 ) with dipole moment p(ω) that points along the unit vectorê i of the i axis. 3 It is therefore convenient to normalize the total power emitted inside a nanostructured medium to the power P (i) 0 (ω, r 0 ) emitted by the same dipole in a homogeneous isotropic medium with the same dielectric constant ε as where the dipole sits in the nanostructure. The normalized power is equal to the ratio of the LDOS in the nanostructured medium and the LDOS in a homogeneous medium with dielectric constant ε, and reads [79,80] Using Poynting's theorem [81], the power P (i) (ω, r 0 ) radiated by the dipole at position r 0 is equal to the inner product of the dipole moment and the local electric field E(ω, r 0 ) at the position of the dipole, where we use complex notation and consider steady state (time average).
To calculate the power radiated by the dipole inside the finite-size photonic crystals we used the open-source implementation MEEP [82] of the finite-difference time domain (FDTD) method [83]. The finite-size crystal is surrounded by a uniform dielectric buffer with the same dielectric constant as that of the low-ε material in the pores. The computational volume is bounded on all sides by perfectly matched layers of thickness a to emulate infinite space. A dipolar point source is placed at the position of interest r 0 with a Gaussian spectrum with a central frequency equal to the midgap frequenciesω mid = 0.58 andω mid = 0.49 for empty and toluene-filled crystals, respectively. The full width at half maximum (FWHM) of the source spectrum was chosen to be equal to ω = 0.8 to cover all the spectral features of interest.
To assess possible numerical artifacts of our method, we have compared the computed LDOS at the center of a dielectric Mie sphere with analytical results [84], where the details are presented in the Appendix. For the best resolution (smallest grid size) and for a frequency range around the central frequency of the Gaussian pulse, we find convergence up to 3% outside Mie resonances, and about 10% near Mie resonances as shown in Fig. 8. The spatial grid size of = a/30 was used in the photonic crystal calculations, since this gave the best match with the analytic test results for a Mie sphere (see the Appendix), while keeping the computation time within reasonable bounds. The calculations were performed on a workstation with an Intel Core i7 processor with 8 CPU cores at 3.4 GHz clock speed and with 32 GB RAM. To keep the simulations tractable, we studied 3D finite crystals with a volume V = N 3 = 3 3 unit cells. The simulation times were equal to 600(a/c 0 ); the real computation time was around 5000 s in order to achieve sufficient convergence of our calculations.

A. Local density of states versus emitter position
We turn to the dependence of the LDOS ρ (i) (ω, r 0 ) on the position r 0 inside the crystal at frequencies ω inside the 3D photonic band gap. We study the LDOS along trajectories in three different high-symmetry directions, where we make sure that all trajectories are centered. First, we consider the LDOS along the axis of the central pore pointing in the z direction, as shown in Fig. 3. Figure 3(a) illustrates the (x = 0, y = 0, z) positions where the LDOS is probed. This set of probe positions are all in the same embedding medium (either air or toluene), which facilitates the interpretation. Figure 3(b) presents the calculated LDOS for the silicon-air crystal at the midgap frequencyω mid = 0.58, and Fig. 3(c) shows the LDOS for the toluene-filled crystal at the midgap frequency (ω mid = 0.49). The silicon-air data strongly decrease from the crystal surface to the center of the crystal. For x-and z-oriented dipoles, the normalized LDOS tends from about 1 to 5 × 10 −2 , corresponding to a relative inhibition of 20× at the center. For the y-oriented dipole, the normalized LDOS tends from about 0.4 to 5 × 10 −2 , corresponding to a relative inhibition of 8× at the center. In the toluene-filled crystal, see Fig. 3(c), similar trends appear, although with smaller inhibitions of about 2× to 3×, since the refractive-index contrast and thus the photonic strength is less than in the airfilled crystal. Aside, we note that the LDOS near the crystal surface is slightly enhanced (for x-and z-oriented dipoles) or slightly decreased for y dipoles, which we tentatively attribute to surface modes [85] or to the fact that the vacuum modes are reflected by the crystal surface thus leading to interference just outside the surface, similar to the well-known Fresnel interference just outside a mirror [81].
Let us briefly compare to the experiments by Leistikow et al. [58], who studied the emission of quantum dot nanocrystals suspended in toluene that were embedded in silicon inverse woodpile structures. In the corresponding Fig. 3(c), we observe a substantial inhibition of the LDOS, in agreement with the experimental observations. In the current situation, the inhibition at the center of the crystal is less (2 to 8 times) than in the experiments (more than 10 times), which is sensible since in the present case the crystal is smaller (N 3 = 3 3 ) than the ones in the experiments (N 3 = 12 3 ). There are aspects where no definite statements can be made, for instance, since the current results pertain to a single dipole that has a definite orientation, whereas in the experiment an ensemble of quantum dots was studied that sampled many positions throughout the whole crystal (80% of the whole volume) and whose dipole orientations were random.
Since the trend of the LDOS versus z position in Fig. 3(b) is exponential within the domain that is computationally tractable here, we interpret the data with a model consisting of two exponentials: The main characteristic is the LDOS decay length (i) ρ that pertains to dipole orientationê i . In the case of a strong inhibition of the LDOS, as is the case in a broad 3D photonic band gap, (i) ρ will be small, and (i) ρ increases for less inhibition. 4 As discussed below, the LDOS length (i) ρ is connected to far-field radiation effects of the dipole. In Eq. (3) each exponential originates from one of the two opposite (x, y) surfaces of the crystal (at z/a = ±1.5/ √ 2 = ±1.06), hence the plus and minus signs with twice the same characteristic LDOS decay length (i) ρ . And A i is a prefactor that equals half the LDOS at the center of the crystal [since ρ (i) (z = 0)/ρ 0 = 2A i ]. As discussed below, A i appears to be connected to near-field effects of the dipole. In the modeling of the computed LDOS data with Eq. (3), we exclude the two data points near the surface to avoid complications due to surface and edge states and Fresnel interference. The solid curves in Fig. 3(b) and Fig. 3(c) are the fitted curves according to Eq. (3) for both silicon-air and silicon-toluene crystals and each of the three dipole orientations. The exponential model tracks the calculated LDOS data better in the toluene-filled crystal than the air-filled crystal, likely since in the former case the LDOS shows a weaker spatial dependence due to the reduced dielectric contrast; hence deviations are expected to be smaller. The resulting LDOS decay lengths and the prefactors are listed in Table I for both air-filled and toluenefilled crystals. Table I shows that for the air-filled crystal the LDOS decay lengths are consistently smaller than for the toluene-filled crystals for all dipole orientationsê x ,ê y , andê z . The shorter LDOS decay lengths are a direct consequence of the higher dielectric contrast in the air-filled crystal, which results in a broader gap (see Fig. 2) and thus stronger inhibitions. In their study on silicon-toluene crystals, Leistikow et al. [58] inferred the LDOS decay length to be equal to about (i) ρ /a = 1. This is in fair agreement (between 3% and 35% greater) with the results in Table I, which is a gratifying consistency between the experimental and computed results.

B. Model parameters and far field and near field
When considering all parameters in Table I, it is instructive to discuss the role of the dipole orientationê i on both the characteristic LDOS length (i) ρ and the amplitude prefactor A i . Starting with the air-filled crystal, we observe that thê e x -andê z -oriented dipoles have smaller LDOS decay lengths ( (x) ρ = 0.286, (z) ρ = 0.267) than theê y -oriented dipole ( (y) ρ = 0.449). This result can be rationalized by a simple model wherein we consider a dipole to have a far-field radiation pattern typical of a homogeneous medium, namely predominantly in its equatorial plane [81]; see Fig. 4(a). Hence, thê e x dipole radiates predominantly in the yz plane in the crystal. The light that would propagate in this plane notably encounters the Y and the Z high-symmetry directions where the gap is wider (and intermediate directions) as seen in Sec. II A. Hence, we naturally expect a strong inhibition in the yz plane, which agrees qualitatively with the small LDOS decay length for theê x orientation. A similar reasoning holds for theê z dipole, whose equatorial plane is the xy plane in the crystal that again includes the Y stop gap, and thus the LDOS decay length is also small. Conversely, in the case of theê y dipole, the equatorial plane is the xz plane in the crystal. This plane contains the relatively narrower X and Z stop gaps (but not the broad Y gap). Hence, less inhibition is expected than for the other orientations, which agrees well with the observed longer LDOS decay length. Thus, we conclude that arguments based solely on the far-field radiation pattern of the dipole located within the photonic crystal serve to explain the relative strength of the characteristic LDOS length observed for different dipole orientations.
We now turn to the role of the dipole orientation on the prefactor A i . Here, we observe that theê z dipole exhibits the smallest prefactor (A z = 0.012), whereas theê x and thê e y dipoles have almost twice greater and closely similar prefactors (A x = 0.021 and A y = 0.026). To understand this behavior, we recall that in the near-field regime a dipole has the strongest field component E(ω, r 0 ) i in the same direction 235309-5  [81], as illustrated in Fig. 4(b). Let us first consider theê y dipole orientation that has the maximum field in the y direction E(ω, r 0 ) y . In the y direction the E(ω, r 0 ) y field crosses the air-silicon interface within a short distance equal to y = 0.12a = 0.12 × 0.585λ = λ/14. Therefore, this near field experiences a polarization in the high-index material that enhances the near field. The enhanced near field is apparently scattered to far-field radiation (by the interface), which leads to an enhanced LDOS and thus a larger prefactor A y . Conversely, theê x andê z dipoles exhibit near fields in the x and z directions where the field is completely inside the uniform air-filled pore. Thus the concomitant near fields experience no polarization enhancement by the silicon, hence the smaller values of A x and A z prefactors. Based of this reasoning we conclude that the A i prefactors are mainly associated with the near-field distributions of theê i -oriented dipoles. For convenience, the results of our discussion are summarized in Table II.
In the case of the toluene-filled crystal, the behavior of the LDOS amplitude A i is similar to that in the silicon-air crystal. As seen in Table I, the dipole with polarizationê z exhibits the smallest amplitude (A z = 0.101), followed by the dipoles polarized alongê y andê x that have similar amplitudes (0.151 versus 0.176). We thus conclude that the near field has the same impact in this case. The characteristic LDOS length (i) ρ , however, does not exhibit the same pattern as in the silicon-air case. As shown in Table I the strongest inhibition in this case, that is, the smaller LDOS length, is found for theê x -oriented dipole, followed by theê y dipole and theê z dipole. The mismatch between the air and toluene cases is possibly caused by the fact that in the silicon-toluene crystal, the reduced refractive index contrast probably leads to an increase of the directional Bragg length to be larger than the crystal size of 3 × 3 × 3 unit cells considered here. In the toluene case, it is probably not meaningful to interpret the LDOS inside the crystal with band structure features, since the infinite crystal is not reached in the computations. In the silicon-air crystal, the Bragg length is sufficiently small compared to the crystal size that the infinite crystal limit is effectively reached; hence a reasoning invoking the interference associated with the stop gaps in the band structure is meaningful.

C. Comparison between LDOS decay length and Bragg length
To put the LDOS decay length in perspective, we compare it to the well-known Bragg length L B [68][69][70] that describes the exponential decay of a directional incident light beam with a frequency inside a photonic gap. This directional decay is described by a nonzero imaginary part of the wave vector k due to Bragg diffraction. The imaginary part of the wave vector is inversely equal to the Bragg length L B .
For silicon-air inverse woodpile crystals with the same pore radii as here, the Bragg length L B was computed by Devashish et al. [77] by the finite-element method. For x-polarized incident plane waves, they found L (x) B = 0.262a, and L (y) B = 0.428a for y-polarized illumination. These values are similar to the LDOS decay lengths (x) ρ = 0.286a and (y) ρ = 0.449a in Table I for x-and y-oriented dipoles, respectively. Considering the difference between the underlying physics, namely the LDOS versus directional propagation, in other words, the real part of the Green's function [86] versus the imaginary part of the Green's function, it is remarkable for the two different length scales to match so closely.
To further support our interpretation, we consider in the schematic in Fig. 5 a dipole at two different positions inside a finite crystal, where we assume the positions to be on the z axis as in Fig. 3. The dipole emits in many different directions in wave vector space (in wave vector space, a crystal as in Fig. 3 is more extended in the horizontal wave vector direction). Let us first consider anê x -oriented dipole that radiates mostly in the (y, z) plane. Since the dipole has a frequency within the band gap, the radiation in any direction will be exponentially damped, since the wave vector is in every direction complex. Thus, the radiation in the z direction is less damped close to the crystal surface (at position r 0,2 ) than deeper inside at position r 0,1 . Radiation in the perpendicular y direction is equally damped at the different positions, since in this direction the dipole is everywhere at the same distance from FIG. 5. Schematic representation of the position vector r = (r x , r y , r z ) and the complex wave vector k = (k x + ik x , k y + ik y , k z + ik z ) in the crystal under study comprising N = 3 unit cells. Inside the band gap where we calculate the LDOS the imaginary part of k is nonzero. The position vector lies along the axis defined in each of the position dependence case studies in Fig. 3, Fig. 6, and Fig. 7. the crystal-vacuum interface (similar considerations pertain to the x direction). Radiation in an oblique direction with wave vector k will also be increasingly damped when the dipole is located at increasing depth in the crystal. The behavior seen in Fig. 3 suggests that apparently the behavior of the LDOS with z position is mostly determined by the z-directed radiation (which is in its purest sense described by the Bragg length), and hardly by the y-directed or other oblique radiation with wave vector k. Thus, whereas the LDOS usually integrates over a broad spectrum of field modes with wave vectors corresponding to all directions, apparently the radiation tending in the closest vacuum-crystal interface dominates the spectrum. On the other hand, the fact that the LDOS integrates over a broad spectrum, instead of a single wave vector as in the Bragg length, explains perhaps why the LDOS decay length is somewhat larger than the corresponding Bragg length.

D. LDOS along different trajectories
We show in Figs. 6 and 7 the normalized LDOS on the trajectories along the x axis and the diagonal path on the xz plane, respectively. In both of these cases the calculations refer only to the silicon-air crystal. Once again, the positions of the emitters are shown in the schematics of the upper panels of Fig. 6 and Fig. 7, respectively, while the LDOS values are While moving across the air-Si interface, the normalized LDOS does not reveal a smooth and continuous behavior, which makes it impossible to use a simple exponential model such as Eq. (3). Indeed, similar strongly varying behavior across the low-and high-index regions within a unit cell has already been noted in Ref. [59] for the much simpler case of a finite-size Bragg stack (also known as a "one-dimensional photonic crystal").
To highlight this behavior in our 3D crystal, we only draw guides to the eye that mark the trend of the LDOS in each direction. They are marked as black solid curves in Fig. 6(b) and Fig. 7(b). In both cases, it appears that the LDOS reveals abrupt variations while tending across the Si regions, which are highlighted in gray. Interestingly, for the x polarization when moving toward the vacuum-crystal interface from the center, the LDOS shifts down across the silicon region in the LDOS calculated on the x axis [ Fig. 6(b)], whereas it shifts up in the LDOS values calculated on the diagonal of the xz plane [in Fig. 7(b)].

E. Practical consequences
Let us briefly discuss a number of practical implications of our work, namely how to apply 3D photonic band gaps to emission control, quantum information processing, and photovoltaics.
In the field of spontaneous control, since the radiative rate is proportional to the LDOS, controlling the LDOS is a key step [12]. Hence it is clear that a 3D photonic band gap offers an extreme spontaneous emission control. In the field of photovoltaics, it has been realized that an efficient absorber is equivalent to an efficient emitter [87]. Hence a 3D photonic band gap could offer a control means to photovoltaics. In the field of quantum information science, it is relevant to shield qubits from ubiquitous vacuum fluctuations that lead to decoherence of the quantum states [34,88]. One solution to this challenge is to place the qubits (assuming they are dipolar) in a 3D photonic band gap that covers the relevant frequency range of the qubits. Our work provides a design rule, namely where to place a dipolar emitter inside a photonic band gap crystal for a certain emission control, and equivalently, where to place a dipole for a certain absorption control, and again equivalently where to place a qubit for a certain decoherence control.
For instance, if one requires the density of vacuum fluctuations-hence the LDOS-to be shielded by a factor 10×, Fig. 3 shows that this is feasible for dipoles placed anywhere between −0.5a z +0.5a about the center. For dipoles operating at optical frequencies corresponding to 1500 nm in the telecom range, this position freedom corresponds to a relatively large range of about 700 nm. A slight limitation to our study is that we only consider positions in the low-index medium of the photonic crystal nanostructure, although these positions occupy no less than 80 vol % of the whole crystal volume [89]. The results in Fig. 3 also show that a tenfold shielding of the vacuum fluctuations is robust with respect to the orientation of the transition dipole moment of the dipole.
It is exhilarating that a silicon-air crystal has a significant inhibition of the LDOS, in view of the relatively small crystal size of only V = 3 × 3 × 3 unit cells. For dipolar emitters (qubits) operating at optical frequencies corresponding to 1500 nm, this would correspond to a small 3D silicon nanophotonic device with a volume as small as V = 4.2 μm 3 . Such a robustness with respect to the crystal size is due to the small LDOS decay length that is much less than one lattice spacing. In parallel to this paper, an experimental study of the directional stop bands of (necessarily finite) 3D photonic band gap crystals [90] is also reaching the conclusion that relatively small micron-sized crystals are powerful tools to control directional transport.

IV. CONCLUSIONS
In this paper, we have presented a computational study of the inhibition of the LDOS in the 3D photonic band gap of a finite-size 3D photonic crystal. In particular, we focused on crystals with the silicon inverse woodpile structure that were recently studied experimentally. To this end, we considered the LDOS dependence on the emitter's position and orientation inside the crystal. Our calculations showed that except for special cases, it is generally not possible to model the LDOS decrease away from the vacuum-crystal interface with a simple exponential model. However, where the exponential model did work, the LDOS decay length turned out to be surprisingly similar to the Bragg length. As for the impact of crystal size on the LDOS suppression, we found that a crystal only as large as comprising 3 × 3 × 3 unit cells and with good dielectric contrast (silicon-air) already provided more than ten times the inhibition of the LDOS around its center. Therefore, for experiments designed to shield quantum systems from vacuum fluctuations, very small volume devices may well be sufficient to fulfill the design requirements on LDOS suppression.

ACKNOWLEDGMENTS
It is a great pleasure to thank D. Devashish

APPENDIX: NUMERICAL CALCULATION OF THE LDOS
We numerically calculate the LDOS by relating the electric field at the location of the electric point dipole to the power radiated by the dipole; see Eq. (2). The electric field is obtained by placing a point dipole source at point r 0 with a dipole moment parallel to the x, y, or z axes. The transient amplitude of the dipole moment is described by a short Gaussian pulse to generate sufficient bandwidth to cover the entire frequency spectrum of interest. After the initial excitation has vanished, we obtain the electric field component parallel to the dipole moment at r 0 versus time t at every time step and take the Fourier transform to obtain the frequency-resolved field E(ω, r 0 ). This approach has also been used in earlier studies too; see for instance Refs. [62,[91][92][93].
To validate the calculation of the LDOS with the MEEP FDTD code, we compare the results of the FDTD method with the analytical results, namely, the modification of the LDOS at the center of a dielectric sphere [84,94]. We consider a sphere of radius a = 1 (reduced units) made of a dielectric material with real dielectric constant ε b = 12.1. Figure 8 associated with the Mie resonances of the sphere. The numerical results were obtained with a Gaussian pulse centered atω = 0.4 and width ω = 0.9. The numerically computed LDOS using the FDTD method is shown in the same figure as discrete points for various grid resolutions defined as the number of sampling points over a radial distance. Good agreement is found between the analytical and numerically computed LDOS specially at higher resolution (30 grid points).
In Fig. 8(b), we quantify the convergence between analytical and numerical results by showing the relative difference (in percentage) between the numerical and analytical results. We observe that the convergence is better at frequencies outside the resonances-up to 3% near the central frequency 0.4 of the spectrum-than for the frequencies around the resonances, typically up to 10% near the central frequency 0.4. The most extreme differences appear at the upper edge of the spectrum where the precision is limited by the low intensity of the excitation pulse in the computation. This is expected due to the fact that the lifetime of the modes on-resonance is much greater, and hence, more computational time is required for these frequencies. The numerical calculations are in excellent agreement with the analytical results, to within 10% near resonances and 3% for off-resonance frequencies. Therefore, we conclude that the simulated results provide a faithful representation of the physics under study.