Optimization of Anisotropic Photonic Density of States for Raman Cooling

TL;DR

This paper explores how anisotropic photonic density of states can be optimized to enhance Raman cooling efficiency in optical refrigeration, addressing the challenge of imbalance in Stokes and anti-Stokes scattering.

Contribution

It introduces a framework for optimizing Raman cooling using anisotropic photonic DoS and demonstrates its importance for systems with lower symmetry.

Findings

Optimization of Raman cooling figure of merit considering crystal orientations

Anisotropic photonic DoS significantly improves cooling efficiency

Applicable to various laser cooling methods with anisotropic media

Abstract

Optical refrigeration of solids holds tremendous promise for applications in thermal management. It can be achieved through multiple mechanisms including inelastic anti-Stokes Brillouin and Raman scattering. However, engineering of these mechanisms remains relatively unexplored. The major challenge lies in the natural unfavorable imbalance in transition rates for Stokes and anti-Stokes scattering. We consider the influence of anisotropic photonic density of states on Raman scattering and derive expressions for cooling in such systems. We demonstrate optimization of the Raman cooling figure of merit considering all possible orientations for the material crystal and two example photonic crystals. We find that the anisotropic description of the photonic DoS and the optimization process is necessary to obtain the best Raman cooling efficiency for systems having lower symmetry. This general result applies to a wide array of other laser cooling methods in the presence of anisotropy.

Figures

Equations

$$W_S=\frac{1}{V}\sum _{\mathbf{k}}{n}_i(n_0+1)|{\mathcal{M}}_S(\mathbf{k}){|}^2\delta ({\omega }_i-{\omega }_0-\omega (\mathbf{k}))$$

$$W_{AS}=\frac{1}{V}\sum _{\mathbf{k}}{n}_i{n}_0|{\mathcal{M}}_{AS}(\mathbf{k}){|}^2\delta ({\omega }_i+{\omega }_0-\omega (\mathbf{k}))$$

$$W_S=\iint \frac{k^{2}dkd\Omega }{(2\pi {)}^3}{n}_i(n_0+1)|{\mathcal{M}}_S(\hat{\Omega }){|}^2\delta ({\omega }_S-ck),$$

$$\frac{d{W}_S}{d\Omega }=n_i(n_0+1)|{\mathcal{M}}_S(\hat{\Omega }){|}^2{D}_0({\omega }_S),$$

$$D_0(\omega )=\int \frac{k^{2}dk}{(2\pi {)}^3}\delta (\omega -ck)=\frac{{\omega }^2}{(2\pi c{)}^3}$$

$$W_S=\sum _{\alpha }$$

$$\times \delta ({\omega }_S-{\omega }_{\alpha }(\mathbf{k})),$$

$$\frac{d{W}_S}{d\Omega }$$

$$\times \delta ({\omega }_S-{\omega }_{\alpha }(\mathbf{k}))\delta (\hat{\Omega }-{\hat{\Omega }}_{\alpha }(\mathbf{k}))$$

$$=n_i(n_0+1)|{\mathcal{M}}_S(\hat{\Omega }){|}^{2}D({\omega }_S,\hat{\Omega }),$$

$$D(\omega ,\hat{\Omega })\equiv \sum _{\alpha }{\int }_{\text{FBZ}}\frac{d\mathbf{k}}{(2\pi {)}^3}\delta (\omega -{\omega }_{\alpha }(\mathbf{k}))\delta (\hat{\Omega }-{\hat{\Omega }}_{\alpha }(\mathbf{k})).$$

$$W_S$$

$$=n_i(n_0+1)\int |{\mathcal{M}}_S(\hat{\Omega }){|}^{2}D({\omega }_S,\hat{\Omega })d\Omega ,$$

$$\delta (\hat{\Omega }-{\hat{\Omega }}_{\alpha }(\mathbf{k}))=\delta (\cos \theta -\cos {\theta }_{\alpha }(\mathbf{k}))\cdot \delta (\varphi -{\varphi }_{\alpha }(\mathbf{k}))$$

$$D(\omega ,\hat{\Omega })\propto \sum _{\alpha }{\int }_{\text{FBZ}}\frac{d\mathbf{k}}{(2\pi {)}^3}\delta (\omega -{\omega }_{\alpha }(\mathbf{k}))\delta (\hat{\Omega }\cdot {\hat{\Omega }}_{\alpha }(\mathbf{k})-1).$$

$$D(\omega )=\int D(\omega ,\hat{\Omega })d\Omega$$

$$\frac{d{S}_S}{d\Omega }={\left(\frac{{\omega }_S}{c}\right)}^4\frac{N\hslash }{M{\omega }_0}(n_0+1)|R(\hat{\Omega }){|}^2$$

$$\frac{d{S}_{AS}}{d\Omega }={\left(\frac{{\omega }_{AS}}{c}\right)}^4\frac{N\hslash }{M{\omega }_0}{n}_0|R(\hat{\Omega }){|}^2$$

$$\frac{dS}{d\Omega }=\frac{n}{n_{i}c}\frac{dW}{d\Omega },$$

$$|R{|}^2=|{\mathbf{e}}_i^T\cdot R_{\mathbf{\xi }}\cdot {\mathbf{e}}_s{|}^2,$$

$$R_x=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& d\\ 0& d& 0\end{array}\right),R_y=\left(\begin{array}{ccc}0& 0& d\\ 0& 0& 0\\ d& 0& 0\end{array}\right),R_z=\left(\begin{array}{ccc}0& d& 0\\ d& 0& 0\\ 0& 0& 0\end{array}\right),$$

$$R_{\mathbf{\xi }}^{\prime }=U{R}_{\mathbf{\xi }}{U}^T,$$

$$U_{[011]}=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1/\sqrt{2}& -1/\sqrt{2}\\ 0& 1/\sqrt{2}& 1/\sqrt{2}\end{array}\right),$$

$$U_{[111]}=\left(\begin{array}{ccc}1/\sqrt{2}& -1/\sqrt{2}& 0\\ 1/\sqrt{6}& 1/\sqrt{6}& -2/\sqrt{6}\\ 1/\sqrt{3}& 1/\sqrt{3}& 1/\sqrt{3}\end{array}\right),$$

$$S_S={\left(\frac{{\omega }_S}{c}\right)}^4\frac{N\hslash }{M{\omega }_0}(n_0+1){\int }_{\Omega }\frac{D({\omega }_S,\hat{\Omega })}{D_0({\omega }_S)}|R(\hat{\Omega }){|}^{2}d\Omega ,$$

$$S_{AS}={\left(\frac{{\omega }_{AS}}{c}\right)}^4\frac{N\hslash }{M{\omega }_0}{n}_0{\int }_{\Omega }\frac{D({\omega }_{AS},\hat{\Omega })}{D_0({\omega }_{AS})}|R(\hat{\Omega }){|}^{2}d\Omega .$$

$$\alpha <\frac{{\omega }_0}{{\omega }_{AS}}{S}_{AS}-\frac{{\omega }_0}{{\omega }_S}{S}_S,$$

$$\text{FoM}\equiv \frac{{\omega }_0}{{\omega }_{AS}}{S}_{AS}-\frac{{\omega }_0}{{\omega }_S}{S}_S.$$

Full text

I Introduction

Since the first proposal of the concept of optical cooling by Peter Pringsheim in 1929 Pringsheim (1929), various experimental demonstrations have been produced in atomic gases Hänsch and Schawlow (1975); Dalibard and Cohen-Tannoudji (1985); Aspect et al. (1988) and in solids Epstein et al. (1995); Seletskiy et al. (2010, 2012); Melgaard et al. (2013, 2014); Zhang et al. (2013). Laser cooling of solids in particular, can be achieved by anti-Stokes photon up-conversion processes in which thermal phonons are annihilated Sheik-Bahae and Epstein (2004, 2009); Ding and Khurgin (2012). The most common approach employs anti-Stokes fluorescence produced by specific electronic transitions in certain rare-earth elements (e.g. Yb3+, Er3+, Tm3+)Epstein et al. (1995); Seletskiy et al. (2010, 2012); Melgaard et al. (2013, 2014). An alternative approach uses strong exciton-phonon interaction in a polar material to achieve cooling via photoluminescence Zhang et al. (2013). Furthermore, in order to minimize non-radiative recombination processes, only direct band gap materials are suitable. Even though these techniques for cooling solids have been remarkably successful Melgaard et al. (2013, 2016), the reliance on specific materials restricts their broad applicability.

Inelastic light scattering, which also provides anti-Stokes photon up-conversion and has fewer restrictions on material properties, has separately emerged as an alternative method for laser cooling of solids. The fundamental challenge for cooling using inelastic light scattering (e.g. Brillouin or Raman scattering) lies in the natural imbalance between the Stokes and the anti-Stokes scattering rates. In a homogeneous material the spontaneous scattering rate for the Stokes process is proportional to $n_{0}+1$ while the scattering rate for the anti-Stokes process is proportional to $n_{0}$, where $n_{0}=[\exp(\hbar\omega_{0}/k_{B}T)-1]^{-1}$ is the Bose-Einstein distribution and $\omega_{0}$ is the phonon frequency. For optical phonons in silicon, for example, this ratio of Stokes transition rate to anti-Stokes transition rate is about 10 to 1 at room temperature. To overcome this challenge researchers have attempted to engineer either the electronic properties Khurgin (2006); Ding and Khurgin (2012); Zhang et al. (2013, 2016) or the optical properties Bahl et al. (2012); Chen and Bahl (2015); Chen et al. (2016); Kim and Bahl (2017) of the system. For instance, successful engineering of electronic properties via excitonic resonance has led to the demonstration of Raman cooling of a longitudinal optical (LO) phonon mode in ZnTe nanobelts Zhang et al. (2016), in which only the anti-Stokes scattering is in excitonic resonance.

The engineering of optical properties, on the other hand, focuses on enhancing the available photonic states Purcell (1946); Gaponenko (2002) for only anti-Stokes scattered photons. Using this concept with optical resonators resulted in the first demonstration of Brillouin cooling of an acoustic phonon mode Bahl et al. (2012), and the understanding that the cooling efficiency is shaped by the photonic density of states (DoS) of the optical modes Kim and Bahl (2017).

In addition to dealing with only a single resonant mode of the system, we further expand this method to incorporate the engineered DoS of all the photonic modes of the system in the context of Raman scattering Chen and Bahl (2015). By combining the effect of Raman selection rules and the photonic DoS, we showed that net Raman cooling can be potentially achieved in a system with low absorption. However, the photonic DoS in the calculation was assumed to be isotropic, which as we will see in the present study is insufficient in dealing with certain photonic systems. Furthermore, the fundamental problem of analyzing Raman scattering in a system considering the full effect of anisotropic photonic DoS has not previously been studied. In this study we produce for the first time an explicit expression for calculating anisotropic DoS in a photonic engineered system and analyze the effect of anisotropic DoS on Raman scattering. We apply this new knowledge to the specific problem of Raman laser cooling in solids and perform a full optimization of the cooling process.

The basic concept of using anisotropic photonic DoS to engineer Raman scattering is illustrated in Fig. 1. The scattering efficiency is strongly influenced by structure at two scales. First, the crystal structure at the microscopic, molecular scale determines available vibrational states and the resulting Raman transitions. For a selected incident light polarization, the response at this microscopic scale will determine the Raman scattering intensity in all directions (for example: left panel of Fig. 1b). If we change the material crystal orientation with respect to a fixed incident polarization, the resulting scattering pattern will also change accordingly (left panel of Fig. 1c).

Second, at the mesoscopic scale, the photonic structure (e.g. photonic crystal) determines the available states for light in the material or the anisotropic photonic DoS (for example: middle panel of Fig. 1b) of the system. If the principal axes of the photonic crystal are rotated, i.e. the patterning of air-spheres is changed, the photonic DoS pattern will also be rotated (middle panel of Fig. 1c).

These two effects at the microscale and mesoscale together determine the Raman efficiency and the resulting scattering pattern is given by the product of the anisotropic DoS and the Raman scattering intensity (right panel of Fig. 1b and 1c). By arranging the relative orientation of the photonic crystal and the material crystal we can thus obtain the best resulting Raman scattering efficiency for various Raman applications. In this paper, we aim to exploit this combined effect at two scales and propose how to optimize the relative orientation for maximizing Raman cooling efficiency.

To demonstrate our approach, we use crystalline silicon as our material crystal since its highly symmetric structure and well-known properties simplify the calculation. In addition, we choose two different photonic crystals with different symmetries to emphasis the effect of anisotropy on Raman cooling efficiency.

To standardize coordinate systems and crystal axes, we use $x$, $y$, and $z$ to denote a fixed coordinate system, which is also used to define the polarization direction of the incident light. Miller indices with bold numbers (e.g. $[\textbf{100}]$, $[\textbf{010}]$) are used to denote directions with respect to the photonic crystal, and the usual Miller indices (e.g. $[100]$, $[11\bar{2}]$) denote directions with respect to the material crystal.

The paper is organized as follows. In section II we review the theory of Raman scattering and produce an expression of anisotropic photonic DoS. In section III we discuss the effect of Raman selection rules on the Raman scattering intensity. Finally, in section IV we discuss the application to Raman laser cooling.

II Raman scattering and anisotropic DoS

In this section, we will provide an explicit expression for calculating anisotropic DoS and discuss how it affects the Raman transition rate. Later, we will directly employ the anisotropic DoS in the calculation of the Raman scattering efficiency, which is directly proportional to the Raman transition rate.

In a bulk homogeneous medium, the total transition rate for Raman scattering with pump frequency $\omega_{i}$ and phonon frequency $\omega_{0}$ is given by Fermi’s golden rule Loudon (1963, 1964):

[TABLE]

for the Stokes process and

[TABLE]

for the anti-Stokes process, where $V$ is the volume of the system, $n_{i}$ is the number of pump photons, and $n_{0}=[\exp(\hbar\omega_{0}/k_{B}T)-1]^{-1}$ is the equilibrium phonon occupation number at temperature $T$. The delta functions in the above equations enforce energy conservation, where $\omega_{i}-\omega_{0}=\omega_{S}$ is the Stokes frequency, $\omega_{i}+\omega_{0}=\omega_{AS}$ is the anti-Stokes frequency, and $\omega(\mathbf{k})$ is the photon-dispersion relation of the system. $\mathcal{M}_{S}$ and $\mathcal{M}_{AS}$ are summations over all possible intermediate transition matrix elements for the Stokes and anti-Stokes process, which include both electron-photon interaction and electron-phonon interaction Loudon (1963). They are in general functions of both frequencies and wave vectors of the photons and phonons participating in scattering process as well as the electron state energies for the virtual transitions. The summations in Eqs. (1) and (2) should also be over all possible final states, which would include both final photon and final phonon states labeled by their respective wave vectors $\mathbf{k}$ and $\mathbf{q}$. However, since we consider only first order non-resonant Raman scattering in this paper, the only phonon modes that participate in the scattering process are the zone center Raman active phonon modes Yu and Cardona (1996). We can therefore neglect the dependence of $\mathcal{M}_{S}$ and $\mathcal{M}_{AS}$ on frequency and phonon momentum $\mathbf{q}$, and sum over only all final photon states labeled by $\mathbf{k}$. In most experiments, the pump frequency is much greater than the phonon frequency $\omega_{i}\gg\omega_{0}$, which allows us to approximate $|\mathcal{M}_{S}(\mathbf{k})|\approx|\mathcal{M}_{AS}(\mathbf{k})|$ Yu and Cardona (1996). As we discuss later in section III, the dependence of $\mathcal{M}_{S}$ and $\mathcal{M}_{AS}$ on $\mathbf{k}$ comes only from the direction of the final photon state $\hat{\Omega}(\mathbf{k})$. We can thus rewrite them as $\mathcal{M}_{S}(\hat{\Omega}(\mathbf{k}))$ and $\mathcal{M}_{AS}(\hat{\Omega}(\mathbf{k}))$.

In a bulk homogeneous medium, the dependence of the transition rate on the direction of scattered light only enters via $\mathcal{M}_{S}$ and $\mathcal{M}_{AS}$ since the photon-dispersion relation is isotropic and given by $\omega(\mathbf{k})=c|\mathbf{k}|$. If we convert the sum into integral in Eq. (1), we have

[TABLE]

where a spherical coordinate system is used, and $\mathbf{k}=(k,\hat{\Omega})$. Here $\mathcal{M}_{S}$ is only a function of direction $\hat{\Omega}$ and is independent of the magnitude $k$. Thus, we can write down the transition rate observed over a unit solid angle as

[TABLE]

where

[TABLE]

is the isotropic photonic DoS per unit solid angle and is independent of $\hat{\Omega}$ Gaponenko (2002). In this case, the photonic density of states factor can be taken out of the integration over solid angle in the calculation of total scattering rate, as demonstrated in Chen and Bahl (2015).

In a photonic crystal system with infinite periodicity, the system is amenable to the standard Bloch treatment and only the modes in the first Brillouin zone are significant. Eq. (1) (similarly for Eq. (2)) then becomes

[TABLE]

where we have converted the sum over $\mathbf{k}$ in Eq. (1) into an integral over all possible wave vectors in the first Brillouin zone. Here we also sum over all bands labeled by $\alpha$. Each wave vector $\mathbf{k}$ is now associated with more than one eigenfrequency. It is therefore important to also label the direction of the wave vector corresponding to each band $\hat{\Omega}_{\alpha}(\mathbf{k})$. In this case both the matrix elements $\mathcal{M}_{S}$ and $\mathcal{M}_{AS}$ and the delta-functions depend on the direction of the scattered light. Due to the anisotropic dispersion relation, we cannot separate the integral over $\mathbf{k}$ into integrals over $dk$ and $d\Omega$ as we did in the isotropic case. However, if we focus our attention on the transition rate for light scattered into a differential solid angle $d\Omega$ around the direction of interest $\hat{\Omega}$, we can write

[TABLE]

and similarly for the anti-Stokes process. Here we have defined the anisotropic photonic density of states along a given direction $\hat{\Omega}$ and at a given frequency $\omega$ as

[TABLE]

The second delta function in the above expression selects only states around the direction of interest $\hat{\Omega}$. The total transition rate can now be written as

[TABLE]

and similarly for the anti-Stokes process.

In carrying out the calculation of the anisotropic DoS in Eq. (8), we can rewrite $\hat{\Omega}$ in terms of polar and azimuthal angles and use the relation

[TABLE]

to evaluate the delta function. However, writing $d\Omega=\sin\theta d\theta d\phi$ may result in uneven spacing for $d\Omega$ at different polar angles. Therefore in our calculation we rewrite Eq. (8) as

[TABLE]

This change in argument in the delta function will introduce an additional constant factor to the resulting DoS. This constant factor can be eliminated by demanding that the integral of the anisotropic DoS in Eq. (9) over entire solid angle gives the total density of states

[TABLE]

With the above expression of anisotropic DoS, we can now evaluate how the photonic DoS varies in space for our particular systems. We choose a diamond-structured photonic crystal and a hexagonal close-packed (HCP) photonic crystal for our DoS calculation. Photonic crystal structures that are similar to our HCP crystal have been discussed in Joannopoulos et al. (2011) and fabricated in Qi et al. (2004) by stacking two dimensional layers. However, the structures reported in these references are the $ABCABC\cdots$ stacking, which effectively converts the photonic crystal into the face-centered cubic structure. Here, our structure is the $ABAB\cdots$ stacking and has the symmetry of the HCP structure.

For the diamond photonic crystal (Fig. 2), air spheres of radius $r=0.3a$ are arranged in diamond structure and patterned into crystalline silicon, where $a$ is the lattice constant of the photonic crystal lattice. For the HCP photonic crystal (Fig. 3), we put air cylinders of radius $r=0.45a$ and height $h=0.75c$ at the corresponding atomic sites of an HCP crystal, where $a$ and $c=\sqrt{8/3}a$ are the lattice parameters for the HCP photonic crystal lattice.

The band structure calculation is carried out using the MIT Photonic-Bands package (MPB) Johnson and Joannopoulos (2001). We use a $24\times 24\times 24$ grid to discretize the unit cell of both diamond and HCP photonic crystals. In the DoS calculation for the diamond photonic crystal, we use 1,110 equally-spaced $k$-points in the irreducible Brillouin zone, and a total of 42,931 $k$-points are used over the first Brillouin zone for the anisotropic DoS calculation. For the HCP photonic crystal, we use 994 equally-spaced $k$-points to discretize the irreducible Brillouin zone. The total number of $k$-points in this case is 22,439 over the first Brillouin zone. This choice of number of $k$-points allows at least 45 $k$-points on average to be around a differential solid angle in the calculation $d\Omega=4\pi/10,000$.

As examples of actual experiments, for both photonic crystals we choose two sets of pump frequencies and plot the anisotropic DoS at the anti-Stokes and Stokes frequencies corresponding to the Raman shift of the triply degenerate optical phonon modes of silicon (one LO and two TO modes, 519 cm*-1*). For the diamond photonic crystal, when the system is pumped at $\omega a/2\pi c=0.594$ the Stokes frequency lies in the complete photonic band gap and thus there are no states available (Fig. 2d). The anti-Stokes frequency in this case is at $\omega a/2\pi c=0.639$ and the corresponding anisotropic DoS has the anisotropic shape shown in Fig. 2c. The peaks in the DoS completely reflect the band structure in Fig. 2a. This can be seen more clearly if we shift the pump to $\omega a/2\pi c=0.625$, leading to an anti-Stokes frequency at $\omega a/2\pi c=0.670$ (Fig. 2e) and a Stokes frequency at $\omega a/2\pi c=0.580$ (Fig. 2f). Now the Stokes frequency lies at the upper edge of the photonic band gap, and there are only states along the L direction in the band diagram. Thus, in the anisotropic DoS plot we see only 8 peaks corresponding to the 8 equivalent L points in the first Brillouin zone.

In comparison to the diamond photonic crystal, the DoS of the HCP photonic crystal exhibits more anisotropic behavior, particularly along the $[\textbf{001}]$ direction, as shown in Fig. 3. In this structure, there is a large band gap for the out-of-plane modes but only a narrow in-plane band gap. Fig. 3c and 3d show the anti-Stokes DoS and the Stokes DoS when the system is pumped at $\omega a/2\pi c=0.478$. The Stokes frequency lies around the bottom edge of the band gap, and have a negligible photonic DoS. The anti-Stokes DoS, on the other hand, has large peaks along the M direction, corresponding to the low-group velocity bands around the M point.

III Raman selection rules and and Raman efficiency with anisotropic DoS

The Raman scattering process in a bulk homogeneous medium can be evaluated using Eqs. (1) and (2). The calculation of the matrix element $\mathcal{M}$ in Eqs. (1) and (2) can be carried out from first principles Deinzer and Strauch (2002); Lazzeri and Mauri (2003). An easier approach is to use the semi-empirical formulas for the Raman scattering efficiency per unit length per unit solid angle Aggarwal et al. (2011); Chen and Bahl (2015); Wagner and Cardona (1983) as follows:

[TABLE]

for Stokes scattering, and

[TABLE]

for anti-Stokes scattering. Here $N$ is the number of unit cells per unit volume, $M$ is the atomic mass, and $\hbar$ is the reduced Planck constant. The Raman tensor $R(\hat{\Omega})$ is directly proportional to the summation over the matrix elements over intermediate virtual electronic states $\mathcal{M}$ in the previous section. However, when written in this form, the independent elements in $R(\hat{\Omega})$ can be directly measured in experiments Loudon (1963, 1964); Yu and Cardona (1996). We also note that the Raman scattering efficiency is related to the Raman transition rate by the relation Loudon (1963)

[TABLE]

where $n$ is the refractive index of the medium, and $n_{i}$ is the number of incident photons. Before we apply the anisotropic photonic DoS to the Raman scattering efficiencies, we first examine how the Raman tensor (or the material crystal) affects the Raman scattering intensity.

For incident light propagating in the direction $\mathbf{k}_{i}$ with polarization $\mathbf{e}_{i}$ scattered into the direction $\mathbf{k}_{s}$ (same as $\hat{\Omega}$ in section II) with polarization $\mathbf{e}_{s}$, the total Raman intensity $|R|^{2}$ can be calculated through Yu and Cardona (1996)

[TABLE]

where $R_{\boldsymbol{\xi}}$ is the Raman tensor for a phonon polarized in the $\boldsymbol{\xi}$ direction (the direction of atomic displacement) and the superscript “$T$” denotes the transpose of the vector. The Raman intensity for this scattering geometry will depend on the choice of $\mathbf{e}_{s}$. However, there is always another possible scattered photon state with polarization $\mathbf{e}_{s}^{\prime}$ that is perpendicular to both $\mathbf{k}_{s}$ and $\mathbf{e}_{s}$. The total contribution to the Raman intensity from these two possible photon states is given by the sum $|\mathbf{e}_{i}^{T}\cdot R_{\boldsymbol{\xi}}\cdot\mathbf{e}_{s}|^{2}+|\mathbf{e}_{i}^{T}\cdot R_{\boldsymbol{\xi}}\cdot\mathbf{e}_{s}^{\prime}|^{2}$, which is the square of the projection of the vector $(\mathbf{e}_{i}^{T}\cdot R_{\boldsymbol{\xi}})$ onto a plane perpendicular to $\mathbf{k}_{s}$. This is thus independent of how the two photon states are decomposed (the choice of $\mathbf{e}_{s}$ and $\mathbf{e}_{s}^{\prime}$), and only depends on the direction of $\mathbf{k}_{s}$, as we mentioned in section II.

The form of the Raman tensor $R_{\boldsymbol{\xi}}$ depends on the symmetry of the corresponding phonon mode Loudon (1964). For diamond-structured crystals ($O_{h}$ point group for zone center phonons) and zinc-blende-type crystals ($T_{d}$ point group for zone center phonons), the matrix representation of the Raman tensor for all three Raman-active phonon modes (one LO, two TO) takes the form

[TABLE]

for phonons polarized along $[100]$, $[010]$, and $[001]$, respectively. Here, $d$ is the only independent Raman tensor element for these crystal structures. For a polar material (zinc-blende crystal), $d$ takes different values for LO and TO due to the phonon-polariton effect Yu and Cardona (1996). For a non-polar material such as silicon, where the LO and the two TO modes are triply degenerate, $d$ takes the same value for all these phonon modes. If the phonons are polarized in a more arbitrary direction $\boldsymbol{\xi}$, $R_{\boldsymbol{\xi}}$ will be the corresponding linear combination of the $R_{x}$, $R_{y}$, and $R_{z}$ matrices.

In our calculations for the selected system (crystalline silicon), we fix the incident light propagation to be along the $x$-axis $(\mathbf{k}_{i}=\mathbf{x}$), with a $z$-polarization of the $\mathbf{E}$ field ($\mathbf{e}_{i}=\mathbf{z}$), and allow the scattered wave vector $\mathbf{k}_{s}$ to be in any direction within the entire $4\pi$ solid angle. For each $\mathbf{k}_{s}$, we include the contributions to the Raman intensity from two possible photon polarizations that are perpendicular to the $\mathbf{k}_{s}$ direction. The wave vector $\mathbf{q}$ of the phonon participating in the scattering process for a selected pair of $\mathbf{k}_{i}$ and $\mathbf{k}_{s}$ is obtained from momentum conservation $\mathbf{q}=\mathbf{k}_{i}-\mathbf{k}_{s}$ (for Stokes scattering). The direction of phonon polarization $\boldsymbol{\xi}$ is simply $\mathbf{q}$ for the LO phonon mode, or is two orthogonal unit vectors perpendicular to $\mathbf{q}$ for the two TO phonon modes, respectively. The corresponding Raman scattering intensity from the Raman selection rules is then evaluated through Eq. (13).

When the material crystal is rotated with respect to a fixed incident polarization, the Raman tensor for the rotated crystal is related to the original Raman tensor by a unitary transformation

[TABLE]

where $U$ is the Euler rotation matrix for the given rotation. Here we note that fixing the incident polarization while rotating the material is equivalent to rotating the incident polarization while fixing the material. Throughout this paper we will use the former operation.

In Fig. 4 we show the Raman scattering intensity $|R(\Omega)|^{2}$ from the Raman selection rules described above for crystalline silicon. Note that we have not included the effect of anisotropic DoS and all the values are normalized by $|d|^{2}$. Here we choose three different orientations in which either the $[001]$, $[011]$, or $[111]$ direction coincides with the incident polarization ($\mathbf{e}_{i}=\mathbf{z}$) axis. The rotation matrices (corresponding to Eq. (15)) for the $\mathbf{z}=[011]$ and $\mathbf{z}=[111]$ orientations are given by

[TABLE]

and

[TABLE]

respectively. Silicon wafers are typically manufactured with three different choices of the out-of-plane axis, specifically the $[001]$, $[011]$, or $[111]$ axis of the crystal structure. It is thus natural to choose these out-of-plane axes as the directions of incident polarization for light, assuming an in-plane waveguiding structure. In the next section we will further allow the material crystal to rotate about this axis (the $\mathbf{e}_{i}=\mathbf{z}$ axis), which can be though as rotating the wafer, or chip, in a real experimental setup.

From Fig. 4 we see that when the material crystal is rotated with respect to a fixed incident polarization, the pattern of Raman intensity is also rotated correspondingly. For the $[111]$ orientation, the shape of the total scattering pattern also changes. However, we note that the total Raman intensity integrated over all directions remains unchanged for all cases.

We can now include the effect of anisotropic DoS to the Raman scattering efficiencies. To do so, we first note that Eqs. (11) and (12) are obtained for scattering processes in free space, and the photonic DoS in vacuum $D_{0}(\omega)$ (Eq. (5)) is assumed. Following our treatment in the previous section to replace the vacuum DoS $D_{0}(\omega)$ by the actual anisotropic DoS $D(\omega,\hat{\Omega})$ and integrate over the entire solid angle, we obtain the total Raman scattering efficiency per unit length Chen and Bahl (2015):

[TABLE]

[TABLE]

From Eqs. (16) and (17), we see that the term in the integral has two of degrees of freedom (in $R(\hat{\Omega})$ and $D(\omega,\hat{\Omega})$) and can be engineered to obtain the best scattering efficiency for a cooling process. In the next section we will apply this principle to optimize the Raman laser cooling process.

In practice, one usually collects scattered light within a finite solid angle given by the experimental apparatus. The integral in Eqs. (16) and (17) are then integrated over the solid angle in which scattered light is collected.

IV Optimization of figure of merit for Raman cooling

Now that we know the analytical forms of both the anisotropic photonic DoS and Raman intensity, we can optimize the Raman scattering efficiencies in Eqs. 16 and 17 for laser cooling. We first recall that the minimum requirement to achieve net cooling via Raman scattering from a purely spontaneous process is given by Chen and Bahl (2015)

[TABLE]

where $\alpha$ is the absorption coefficient (per unit length) of the material. We note that in the above cooling condition, it is assumed that the scattered light can escape the system with 100 % efficiency. Without further knowledge of the absorption coefficient, the most import parameter for Raman cooling is the right hand side of Eq. (18), which we define as the figure of merit

[TABLE]

The optimization of a Raman cooling process is therefore reduced to maximizing the FoM.

In the calculation of $S_{S}$ and $S_{AS}$ (Eqs. 16 and 17) for our particular systems, we use the properties of crystalline silicon, $M=28.09$ amu, $N=2.5\times 10^{22}$ cm*-3*, $|d|=1.9\times 10^{-15}$ cm2 and the phonon energy is set to 519 cm*-1* Aggarwal et al. (2011) corresponding to the triply degenerate zone-center phonon modes (one LO and two TO modes).

We use the two photonic crystals discussed in section II as our engineered DoS systems. Here we choose $a=861$ nm for the diamond photonic crystal such that when pumped at around 207 THz (1448 nm exhibits extremely low optical absorption for undoped silicon Green (2008)) the Stokes sideband lies within the photonic band gap.

The choice of lattice constant for the HCP photonic crystal is more crucial than that of the diamond photonic crystal. Without the knowledge of the photonic DoS, one would expect that a good choice of maximizing the FoM would have the Stokes frequency inside the band gap. However, by observing the DoS in Fig. 3a, we see that if the Stokes sideband is moved into the band gap, the anti-Stokes sideband will be at around a local DoS minimum, which significantly reduces the FoM. We therefore choose $a=693$ nm, which places the Stokes frequency at the bottom edge of the photonic band gap and notably gas nonzero DoS (Fig. 3d) and simultaneously place the anti-Stokes at a DoS peak when pumped at 207 THz.

We can now adjust the relative orientation between the photonic crystal and the material crystal to obtain the highest FoM. To explore the optimization space, we examine cases where we set the $[001]$, $[011]$, or $[111]$ orientation of the material crystal as the $z$-axis of the fixed coordinate system, which is also the direction of incident polarization. The material crystal is now allowed to undergo rotation on a $\mathbf{z}$ directed vector by an angle $\theta_{z}$, while keeping incident wave vector $\mathbf{k}_{i}=\mathbf{x}$ and polarization $\mathbf{e}_{i}=\mathbf{z}$ fixed (illustrated in Figs. 6-8). After fixing a particular material orientation, we allow the photonic crystal to take any orientation within the $4\pi$ solid angle. In our calculations, we specify the orientation of the photonic crystal using the polar and the azimuthal angles $\theta$ and $\phi$ of the $[\textbf{001}]$ axis with respect to the fixed coordinate system (not to be confused with $\theta$ and $\phi$ in section II). In both diamond and HCP cases, we divide $\theta$ and $\phi$ into 100 equally-spaced segments. A total of 10,000 different orientations are used for each photonic crystal in the calculation for a chosen material orientation.

Before discussing the effects of anisotropic DoS on Raman laser cooling, we first highlight the importance of anisotropy on the measurement of scattering. We can thus define a directional FoM, where $S_{S}$ and $S_{AS}$ are only integrated over the solid angle of light collection in an actual experiment. This directional FoM characterizes the fractional contribution to the total cooling by the light scattered within a given solid angle. In Figs. 6 and 6 we show the forward directional FoM for the diamond and the HCP photonic crystals respectively as a function of the photonic crystal orientation for each fixed material crystal orientation (specified by a selected crystal axis as the $z$-axis, and a selected $\theta_{z}$). The scattering efficiencies $S_{S}$ and $S_{AS}$ are integrated over a cone of apex angle $\psi=1$ rad, centered around the positive $\mathbf{x}$ direction, which corresponds to a solid angle of $\Omega=2\pi(1-\cos\frac{\psi}{2})\approx$ 0.77 sr in the forward direction. For each 3D surface plot, the direction of a point on the surface represents the direction of the $[\textbf{001}]$ axis of the photonic crystal with respect to the fixed coordinate system, and the distance between a point on the surface to the origin indicates the value of directional FoM for the particular relative orientation. Note that for the $\mathbf{z}=[001]$ and the $\mathbf{z}=[111]$ material orientations, rotating the material crystal about the $z$-axis does not change the FoM since their total Raman scattering patterns are symmetric with respect to the $z$-axis (Fig. 4a and c). Thus only one FoM surface plot is shown in Figs. 6 and 6 for the $\mathbf{z}=[001]$ and the $\mathbf{z}=[111]$ cases.

On the other hand, for changing photonic crystal orientation, we observe that for both the diamond and the HCP photonic crystals the forward directional FoM varies significantly. At certain photonic crystal orientations, almost no forward scattering is observed. The maximum forward directional FoM occurs for the diamond photonic crystal when the $[\textbf{001}]$ axis coincides with the $\mathbf{x}$ axis (and other equivalent orientations). The HCP photonic crystal, on the other hand, gives a maximum forward directional FoM when the $[\textbf{001}]$ axis lies in the $xy$ plane and makes a 69.1∘ angle from the $x$-axis.

The effect of Raman selection rules can be clearly seen when we rotate the material crystal for the $\mathbf{z}=[011]$ case. As we increase $\theta_{z}$ the total Raman scattering pattern in Fig. 4b also rotates about the $z$-axis, resulting in a larger forward directional FoM at $\theta_{z}=\pi/2$ in Figs. 6 and 6. This shows that the combined knowledge of anisotropic photonic DoS and Raman selection rules is crucial here, since in an actual experiment one may not observe any scattered light if the apparatus is designed to collect light at certain direction only.

We can now discuss total Raman laser cooling from scattered light in all directions. In Figs. 8 and 8 we plot the total FoM for the diamond photonic crystal and the HCP photonic crystal respectively, where the scattering efficiencies $S_{S}$ and $S_{AS}$ in Eq. (19) are integrated over entire $4\pi$ solid angle. We see immediately that since the total FoM includes scattered light in all direction, the variation in its value is much less than the directional FoM. Overall, the FoM surfaces for the diamond photonic crystal are close to spherical, showing that the FoM for the diamond photonic crystal is relatively insensitive to changes in both material crystal orientation and photonic crystal orientation. This is due to the fact that the diamond photonic crystal belongs to the space group $O_{h}^{7}$, which has high symmetry along all directions, and the change in anisotropic DoS under rotation does not contribute too much to the integrated Raman efficiency. The largest difference in FoM among all orientations is only about 4 %.

The HCP photonic crystal, on the other hand, belongs to the space group $D_{6h}^{4}$, where properties along the out-of-plane direction behave differently from properties in the in plane directions. Thus, we see a larger difference in the FoM between cases for $[\textbf{001}]=\mathbf{z}$ and $[\textbf{001}]=\mathbf{x}$. The largest FoM for the HCP is obtained when $[111]=\mathbf{z}$ and $[\textbf{001}]=\mathbf{z}$. This can be understood from Fig. 3c and Fig. 4c, where the large in-plane Raman scattering from the selection rules is enhanced by the anti-Stokes DoS peaks along the M directions. The largest FoM among all these orientations is up to 19 % greater than the FoM of the worst orientation.

In these calculations for the selected examples, we use silicon as the material crystal of the system, which, as mentioned in section III, has the same independent Raman tensor element $d$ in Eq. (14) for both LO and TO phonons. It is worth noting that if instead a polar material is used for Raman cooling, $d$ will be different for the LO and TO phonons. In such a case, instead of optimizing the total FoM, one can separately optimize the LO FoM or the TO FoM for the more prominent Raman mode. For example, in the excitonic resonant Raman cooling experiment reported in Zhang et al. (2016) (ZnTe, zinc-blende structure), the dominant phonon mode in the cooling process is the LO mode. A photonic crystal orientation that maximizes the LO FoM could be used to further enhance the cooling efficiency in the experiment.

In Figs. 9 and 10 we plot the FoM for our particular examples as a function of pump frequency, and compare against available experimental data on optical absorption coefficients. Here, the material crystal orientation is denoted by a choice of $z$-axis ($\mathbf{z}=[001]$, $\mathbf{z}=[011]$, or $\mathbf{z}=[111]$) and the angle $\theta_{z}$, while the photonic crystal is denoted by $\theta$ and $\phi$, which are the polar and the azimuthal angles of the $[\textbf{001}]$ axis with respect to the fixed coordinate system. For the diamond photonic crystal, all FoM curves are almost identical and no significant change in FoM is observed. The large photonic band gap allows a wide range of frequencies (207-216 THz) for potential Raman cooling (when Eq. (18) is satisfied). In this case, the results are similar to those calculated under the isotropic DoS assumption Chen and Bahl (2015).

We can also see that overall the FoM for the HCP photonic crystal is smaller than that of the diamond photonic crystal. This is due to the reduced dielectric material to air volume ratio in the HCP structure in order to open a band gap. However, the importance of arranging the orientation of anisotropic DoS with respect to the material crystal can be seen more clearly in Fig. 10 for the HCP case. In particular, without the knowledge of the anisotropic DoS of the system, a natural choice for fabricating such a photonic crystal would be using the most easily available $\mathbf{z}=[001]$ oriented wafer (red curve in Fig. 10a). By simply using a $\mathbf{z}=[111]$ wafer in place of the $\mathbf{z}=[001]$ wafer, we can obtain an improvement of of FoM by 19 % for a 207 THz pump.

V Conclusion

We have provided the first mathematical description of how Raman scattering is influenced by anisotropic photonic DoS and Raman selection rules. Using this, we demonstrated a optimization procedure with anisotropic DoS for laser cooling by going through all possible relative crystal orientations. In the specific cases that we considered, the highly-symmetric diamond photonic crystals showed only small change in the Raman cooling FoM with respect to changes in material and photonic crystal orientations. For such photonic crystals, the isotropic DoS assumption used in Chen and Bahl (2015) is sufficient to predict the range of pump frequencies for potential Raman cooling.

On the other hand, for photonic systems with lower symmetry (e.g. waveguides), the anisotropic-DoS optimization of the FoM is of greater importance. In the specific HCP photonic crystal that we considered in this paper, the difference between optimized and unoptimized FoM can be up to 19 % and is dictated purely by a choice of crystal orientations. For such systems, isotropic DoS alone is not sufficient to predict the best available cooling efficiency of the system. The additional knowledge provided by the anisotropic DoS and the optimization procedure that examines all possible relative orientations give us the capability to identify the best orientations and frequency range for achieving Raman laser cooling.

We note that even though the 19 % enhancement in total FoM may not be significant, the combined knowledge of both anisotropic DoS and Raman selection rules is crucial for both high- and low- symmetry photonic crystals when only directional scattering is considered. In the cases we considered, depending on the relative orientations the maximum directional FoM can be 20 times greater than the minimum (for the HCP case), which is of great importance in actual experiments.

Our formalism not only provides an approach to analyze net Raman cooling, but can also function as a method to improve the efficiency for other laser cooling methods in semiconductors Zhang et al. (2013, 2016) by aligning the best anisotropic DoS at the frequency of emission.

Acknowledgements.

This work was funded by a US Army Research Office Grant (W911NF-15-1-0588), a US Air Force Office for Scientific Research grant (FA9550-14-1-0217), a University of Illinois Campus Research Board grant (RB15183) and a National Science Foundation grant (DMR-1555153). I.G. acknowledges support from the Illinois Scholars Undergraduate Research Program.