Three-Dimensional Coupled-Wave Theory in Photonic Crystal

Turn to CWT in PCSEL

Coupled-Wave Theory

Photonic crystal surface-emitting lasers (PCSELs) can realize one-chip, ultra-large-area, coherent semiconductor lasers. However, it is difficult to accurately evaluate the characterization of PCSELs by traditional simulation method such as FDTD or FEM etc. because of ultra-large computation. Thus, the coupled-wave theory (CWT), which is a semi-analytic computational method, is proposed by research group of Susumu Noda from Kyoto University.
The coupled-wave theory in PCSELs traces its history back to 2006 (Ref. 1), and is expanded to three dimension by incorporating the surface emission and high-order coupling effects (Ref. 2). 3D-CWT has been fully explained in the doctoral thesis of Yong LIANG (Ref. 3) and been redefined the coupled matrix by Takuya Inoue (Ref. 4).
In this post, I have changed the symbols of some formulas to avoid conflicts.

Infinite Phonic Crystal

Hermitian and non-Hermitian optical couplings inside PCSELs. a Hermitian couplings between the four fundamental waves $(R_x, S_x, R_y, S_y)$ inside a PCSEL. b Non-Hermitian couplings via radiated waves, where a backside reflector is used for the control of the non-Hermitian coupling coefficient $i/mu$.

Let us start from photonic crystal with a square lattice of air holes which are not so symmetrical in a multilayer dielectric slab. The four basic waves ($R_x$, $S_x$, $R_y$, $S_y$) which play main role in photonic crystal are in the same direction as the primitive vectors of the cell. They represent the amplitudes of waves propagating in their direction.

Theory Derivation

Consider transverse electric (TE) polarization in photonic crystal, because it is normally laser mode in MQW structure.

$$\nabla \times \nabla \times \vec{E}(\vec{r}) = k_0^2 \tilde{n}^2(\vec{r})\vec{E}(\vec{r})$$

Where $\vec{E}$is electric field, $k_0$ is free-space wave vector, $\tilde{n}$ is complex refractive index and $\vec{r}$is position vector. $k_0^2 \tilde{n}^2(\vec{r})\simeq k_0^2n^2(\vec{r})+2\mathrm{i}k_0n_0(z)\tilde{\alpha}(z)$, where $n^2(\vec{r})$ is periodic, $n_0(z)$ is the average refractive index of the material at $z$, and $\tilde{\alpha}(z)$ is gain ($>0$) or loss ($<0$) in each region.
Because photonic crystal is periodic structure, it is obvious the electric field is also periodic. By using Bloch’s theorem, $\vec{E}(\vec{r})=(E_x(\vec{r}),E_y(\vec{r}),0)$ can be expanded as

$$E_j(\vec{r})=\sum_{m,n}E_{j,m,n}(x,y,z)\mathrm{e}^{-\mathrm{i}m\beta_0 x-\mathrm{i}n\beta_0 y},\ j=x,y$$

and $n^2(\vec{r})$ can be expanded as

$$n^2(\vec{r}) = n^2(z)+\sum_{m,n}\xi_{m,n}(z)\mathrm{e}^{-\mathrm{i}m\beta_0 x-\mathrm{i}n\beta_0 y},\ m\neq0,n\neq0$$

Where $\beta_0$ is propagation constant and in photonic crystal $\beta_0=2\pi/a$ ($a$ is lattice constant), and $\xi_{m,n}(z)$ is the high-order Fourier coefficient term. And then, combining the above three formulas, … (See Ref. 3, page 57).
We have retained the spatial and spatial derivative terms of $E_x$ and $E_y$ with respect to $x$ and $y$, so we can use them in a finite system in next section.
So, We now can obtain four basic waves in

$$\begin{aligned} E_{x,1,0}=0,\ E_{y,1,0}=R_x(x,y)\Theta_0(z) \\ E_{x,-1,0}=0,\ E_{y,-1,0}=S_x(x,y)\Theta_0(z) \\ E_{x,0,1}=R_y(x,y)\Theta_0(z),\ E_{y,0,1}=0 \\ E_{x,0,-1}=S_y(x,y)\Theta_0(z),\ E_{y,0,-1}=0 \\ \end{aligned}$$

Where $\Theta_0(z)$ is the field profiles in the z-direction of four basic waves, which is the same as the field profile of the fundamental guided mode for a multilayer structure with no periodic structure. The profile $\Theta_0(z)$ for basic waves can be determined by the waveguide mode profile,

$$\frac{\partial^2\Theta_0}{\partial z}+\left[\epsilon_0(z)k^2-\beta^2\right]\Theta_0 = 0$$

Where $\beta$ is the propagation constant, which satisfies $\beta \simeq \beta_0$ under condition of $\Gamma$ point resonance. $\beta$ and $\Theta (z)$ can be calculated by employing the transfer matrix method. Then, normalize $\Theta(z)$ as $\int_{-\infty}^{\infty}|\Theta(z)|^2 \mathrm{d}z = 1$.
Now, we can calculate four basic waves. Here, take $R_x$ as an example. Combine them and neglected the second spatial derivative terms because of the slow variation of $R_x$ compared to $\mathrm{e}^{\mathrm{i}\beta_0x}$, we can get

$$(\beta-\beta_0)R_x\Theta_0 + 2\mathrm{i}k_0n_0(z)\tilde{\alpha}(z)R_x\Theta_0-2\mathrm{i}\beta_0\frac{\partial R_x}{\partial x}\Theta_0 = -k_0^2 \sum_{m'\neq1,n'\neq0}\xi_{1-m',0-n'}E_{y,m',n'}$$

… (See Ref. 3, page 59). Finally, we can rewritten equations above in matrix form as

$$\left( \begin{array}{c} \delta + i \frac{\alpha}{2} \end{array} \right) \left( \begin{array}{cc} R_x \\ S_x \\ R_y \\ S_y \end{array} \right) = C\left( \begin{array}{cc} R_x \\ S_x \\ R_y \\ S_y \end{array} \right) + \mathrm{i}\left( \begin{array}{cc} \partial{R_x}/\partial x \\ -\partial{S_x}/\partial x \\ \partial{R_y}/\partial y \\ -\partial{S_y}/\partial y \end{array} \right)$$

If we consider an infinite photonic crystal, derivative terms is equal to zero because of periodicity,

$$\left( \begin{array}{c} \delta + i \frac{\alpha}{2} \end{array} \right) \left( \begin{array}{cc} R_x \\ S_x \\ R_y \\ S_y \end{array} \right) = (C_{Hermitian} + C_{non-Hermitian} + C_{non-\Gamma}) \left( \begin{array}{cc} R_x \\ S_x \\ R_y \\ S_y \end{array} \right)$$

$$C_{Hermitian} = \left( \begin{array}{cccc} \kappa_{11} & \kappa_{1D} & \kappa_{2D+} & \kappa_{2D-} \\ \kappa_{1D}^{*} & \kappa_{11} & \kappa_{2D-}^{*} & \kappa_{2D+} \\ \kappa_{2D+} & \kappa_{2D-} & \kappa_{11} & \kappa_{1D} \\ \kappa_{2D-}^{*} & \kappa_{2D+} & \kappa_{1D}^{*} & \kappa_{11} \end{array} \right)$$

$$C_{non-Hermitian} = \left( \begin{array}{cccc} i\mu & i\mu \mathrm{e}^{i\theta_{pc}} & 0 & 0 \\ i\mu \mathrm{e}^{-i\theta_{pc}} & i\mu & 0 & 0 \\ 0 & 0 & i\mu & i\mu \mathrm{e}^{i\theta_{pc}} \\ 0 & 0 & i\mu \mathrm{e}^{-i\theta_{pc}} & i\mu \end{array} \right)$$

$$C_{non-\Gamma} = \left( \begin{array}{cccc} \beta_{k;1,0}-\beta_{0} & 0 & 0 & 0 \\ 0 & \beta_{k;-1,0}-\beta_{0} & 0 & 0 \\ 0 & 0 & \beta_{k;0,1}-\beta_{0} & 0 \\ 0 & 0 & 0 & \beta_{k;0,-1}-\beta_{0} \end{array} \right) \thicksim \left( \begin{array}{cccc} k_{x} & 0 & 0 & 0 \\ 0 & -k_{x} & 0 & 0 \\ 0 & 0 & k_{y} & 0 \\ 0 & 0 & 0 & -k_{y} \end{array} \right)$$

$$\beta_{k;m,n} = \sqrt{(k_{x}+m\beta_{0})^{2}+(k_{y}+n\beta_{0})^{2}}$$

Where

$$\kappa_{i,j} = -\frac{\omega_0^2}{2\beta_0 c^2} \int \xi_{i,j} |\Theta_0 (z)|^2 \mathrm{d}z$$

$$\zeta_{p,q}^{(r,s)} = -\frac{\omega_0^4}{2 \beta_0 c^4}\iint_{PC}\xi_{p,q}\xi_{-r,-s}G(z,z')\Theta_0(z')\Theta_*(z)\mathrm{d}z'\mathrm{d}z$$

$$\chi_{j,p,q}^{(r,s)} = -\frac{\omega_0^2}{2 \beta_0 c^2} \sum_{\sqrt{m^2+n^2}>1} \xi_{p-m,q-n}\varsigma_{j,m,n}^{(r,s)},\quad j=x,y$$

$$\xi_{m,n} = \frac{1}{a^2} \int_{unitcell} n_{pc}^2 (\textbf{r}) e^{i(m\beta_0x+n\beta_0y)} \mathrm{d}x \mathrm{d}y$$

$$\begin{aligned} \left( \begin{array}{c} \int_{PC}E_{x,m,n}(z) \Theta_0^*(z)\mathrm{d}z \\ \int_{PC}E_{y,m,n}(z) \Theta_0^*(z)\mathrm{d}z \end{array} \right) & = \frac{1}{m^2+n^2} \left( \begin{array}{cc} n & m \\ -m & n \end{array} \right) \left( \begin{array}{cc} -m u_{m,n}^{(1,0)} & -mu_{m,n}^{(-1,0)} & nu_{m,n}^{(0,1)} & nu_{m,n}^{(0,-1)} \\ nv_{m,n}^{(1,0)} & nv_{m,n}^{(-1,0)} & mv_{m,n}^{(0,1)} & mv_{m,n}^{(0,-1)} \end{array} \right) \left( \begin{array}{cc} R_x \\ S_x \\ R_y \\ S_y \end{array} \right) \\ & \triangleq \left( \begin{array}{cc} \varsigma_{x,m,n}^{(1,0)} & \varsigma_{x,m,n}^{(-1,0)} & \varsigma_{x,m,n}^{(0,1)} & \varsigma_{x,m,n}^{(0,-1)} \\ \varsigma_{y,m,n}^{(1,0)} & \varsigma_{y,m,n}^{(-1,0)} & \varsigma_{y,m,n}^{(0,1)} & \varsigma_{y,m,n}^{(0,-1)} \end{array} \right) \left( \begin{array}{cc} R_x \\ S_x \\ R_y \\ S_y \end{array} \right) \end{aligned}$$

Where,

$$\omega_0 = k_0 c$$

$$u_{m,n}^{(r,s)} = k_0^2\iint_{PC} \xi_{m-r,n-s}(z)G_{m,n}(z,z')\Theta_0(z')\Theta_0^*(z)\mathrm{d}z'\mathrm{d}z$$

$$v_{m,n}^{(r,s)} = -\int_{PC} \frac{1}{n_0^2(z)}\xi_{m-r,n-s}(z)|\Theta_0(z)|^2\mathrm{d}z$$

and when reflections between individual layers can be neglected,

$$G_{m,n}(z,z') = -\frac{i}{2\beta_z}\mathrm{e}^{-i\beta_z|z-z'|}$$

$$\beta_z = k_0 n_0(z)$$

Besides,

$$\kappa_{1D} = -\kappa_{2,0}-\chi_{y,1,0}^{(-1,0)}+\mathrm{Re}\left(\zeta_{1,0}^{(1,0)}\right)\mathrm{e}^{i\theta_{pc}}$$

$$\kappa_{2D+} = -\chi_{y,1,0}^{(0,1)}$$

$$\kappa_{2D-} = \chi_{y,1,0}^{(0,-1)}$$

$$\kappa_{11} = \chi_{y,1,0}^{(1,0)}+\mathrm{Re}\left(\zeta_{1,0}^{(1,0)}\right)$$

$$\mu=\mathrm{Im}\left(\zeta_{1,0}^{(1,0)}\right)$$

$$\theta_{pc} = \mathrm{\pi}+\mathrm{arg}\left(\zeta_{1,0}^{(-1,0)}\right)-\mathrm{arg}\left(\zeta_{1,0}^{(1,0)}\right)$$

In Brief

$$\left( \begin{array}{c} \delta + i \frac{\alpha}{2} \end{array} \right) \left( \begin{array}{cc} R_x \\ S_x \\ R_y \\ S_y \end{array} \right) = (C_{Hermitian} + C_{non-Hermitian} + C_{non-\Gamma}) \left( \begin{array}{cc} R_x \\ S_x \\ R_y \\ S_y \end{array} \right)$$

$$C_{Hermitian} = \left( \begin{array}{cccc} \kappa_{11} & \kappa_{1D} & \kappa_{2D+} & \kappa_{2D-} \\ \kappa_{1D}^{*} & \kappa_{11} & \kappa_{2D-}^{*} & \kappa_{2D+} \\ \kappa_{2D+} & \kappa_{2D-} & \kappa_{11} & \kappa_{1D} \\ \kappa_{2D-}^{*} & \kappa_{2D+} & \kappa_{1D}^{*} & \kappa_{11} \end{array} \right)$$

$$C_{non-Hermitian} = \left( \begin{array}{cccc} i\mu & i\mu \mathrm{e}^{i\theta_{pc}} & 0 & 0 \\ i\mu \mathrm{e}^{-i\theta_{pc}} & i\mu & 0 & 0 \\ 0 & 0 & i\mu & i\mu \mathrm{e}^{i\theta_{pc}} \\ 0 & 0 & i\mu \mathrm{e}^{-i\theta_{pc}} & i\mu \end{array} \right)$$

$$C_{non-\Gamma} = \left( \begin{array}{cccc} \beta_{k;1,0}-\beta_{0} & 0 & 0 & 0 \\ 0 & \beta_{k;-1,0}-\beta_{0} & 0 & 0 \\ 0 & 0 & \beta_{k;0,1}-\beta_{0} & 0 \\ 0 & 0 & 0 & \beta_{k;0,-1}-\beta_{0} \end{array} \right) \thicksim \left( \begin{array}{cccc} k_{x} & 0 & 0 & 0 \\ 0 & -k_{x} & 0 & 0 \\ 0 & 0 & k_{y} & 0 \\ 0 & 0 & 0 & -k_{y} \end{array} \right)$$

$$\beta_{k;m,n} = \sqrt{(k_{x}+m\beta_{0})^{2}+(k_{y}+n\beta_{0})^{2}}$$

Here, $C_{Hermitian}$, $C_{non-Hermitian}$, $C_{non-\Gamma}$ is only related to photonic crystal materials and structures (see Theory derivation). And by solving matrix eigenvalues we can get $\delta$ and $\alpha$, where $\delta = (\beta^2-\beta_0^2)/2\beta_0 \simeq \beta-\beta_0=n_{eff}(\omega-\omega_0)/c$. $n_{eff}$ is the effective refractive index for the fundamental guided mode and $\alpha$ is the mode loss.

Finite Phonic Crystal

$$\left( \begin{array}{c} \delta + i \frac{\alpha}{2} \end{array} \right) \left( \begin{array}{cc} R_x \\ S_x \\ R_y \\ S_y \end{array} \right) = C\left( \begin{array}{cc} R_x \\ S_x \\ R_y \\ S_y \end{array} \right) + \mathrm{i}\left( \begin{array}{cc} \partial{R_x}/\partial x \\ -\partial{S_x}/\partial x \\ \partial{R_y}/\partial y \\ -\partial{S_y}/\partial y \end{array} \right)$$

The equation above can be discretized using the staggered grid finite-difference method. And rewrite it in difference form

$$\frac{1}{2} L\left( \begin{array}{c} \delta + i \frac{\alpha}{2} \end{array} \right) \left( \begin{array}{cc} R_x^{j+1,k}+R_x^{j,k} \\ S_x^{j+1,k}+S_x^{j,k} \\ R_y^{j,k+1}+R_y^{j,k} \\ S_y^{j,k+1}+S_y^{j,k} \end{array} \right) = \frac{1}{2} LC \left( \begin{array}{cc} R_x^{j+1,k}+R_x^{j,k} \\ S_x^{j+1,k}+S_x^{j,k} \\ R_y^{j,k+1}+R_y^{j,k} \\ S_y^{j,k+1}+S_y^{j,k} \end{array} \right) + \mathrm{i} \frac{L}{a} \left( \begin{array}{cc} R_x^{j+1,k}-R_x^{j,k} \\ -S_x^{j+1,k}+S_x^{j,k} \\ R_y^{j,k+1}-R_y^{j,k} \\ -S_y^{j,k+1}+S_y^{j,k} \end{array} \right)$$

Circles mean the known coupled matrix C, triangle mean the unknown field components are updated using the finite-difference scheme

The boundary condition of a finite photonic crystal (with a side length $L$) can be defined as

$$R_x(0, y) = S_x(L, y) = R_y(x, 0) = S_y(x, L) = 0$$

In Ref. 3, with the lattice constant $a=295$ nm, and the device length L=70 μm, the laser is divided into 14 sections for which the eigenvalues ($αL$ and $δL$) converge well.
Finally, if we define

$$P_{stim}=\alpha\iint_0^L(|R_x|^2+|S_x|^2+|R_y|^2+|S_y|^2)\mathrm{d}x\mathrm{d}y$$

$$P_{edge}=\int_0^L(|R_x|^2+|S_x|^2)|_{x=0,L}\mathrm{d}y+\int_0^L(|R_y|^2+|S_y|^2)|_{y=0,L}\mathrm{d}x$$

$$P_{rad}=2\kappa_{v.i}\iint_0^L(|\xi_{-1,0} R_x+\xi_{1,0} S_x|^2+|\xi_{0,-1} R_y+\xi_{0,1} S_y|^2)\mathrm{d}x\mathrm{d}y$$

Here, $P_{stim}$ describes the total stimulated power inside the laser structure, $P_{edge}$ represents the power escaping from the edges of the laser cavity (i.e. the in-plane loss), and $P_{rad}$ represents the radiation power emitted from the device surface. Here, $\kappa_{v.i}=-\Im \{ \frac{k_0^4}{2\beta_0}\iint_{PC}G(z,z')\Theta_{0}(z')\Theta_{0}^*(z)\mathrm{d}z'\mathrm{d}z \}$ is closely associated with the out-of-plane coupling remains.

Time-Dependent

Consider transverse electric (TE) polarization in photonic crystal again, but this time is time-dependent:

$$\nabla \times \nabla \times \vec{E}(\vec{r}) = -\frac{1}{c^2}\frac{\partial^2}{\partial t^2} [\tilde{n}^2(\vec{r})\vec{E}(\vec{r})]$$

Where $c$ is the speed of light in vacuum. Then, expand this equation according to Bloch’s theorem (See Ref. 5). Finally, we obtain

$$\begin{aligned} \frac{\partial}{\partial t}\left( \begin{array}{cc} R_x \\ S_x \\ R_y \\ S_y \end{array} \right) = & \frac{c}{n_g}\left[-\mathrm{i}\delta+\frac{\Gamma g(N)-\alpha_{in}}{2}\right]\left( \begin{array}{cc} R_x \\ S_x \\ R_y \\ S_y \end{array} \right) - \frac{c}{n_g} \left( \begin{array}{cc} \partial{R_x}/\partial x \\ -\partial{S_x}/\partial x \\ \partial{R_y}/\partial y \\ -\partial{S_y}/\partial y \end{array} \right) \\ & -\gamma\left( \begin{array}{cc} R_x \\ S_x \\ R_y \\ S_y \end{array} \right) + \frac{c}{n_g}C\left( \begin{array}{cc} R_x \\ S_x \\ R_y \\ S_y \end{array} \right)+\left( \begin{array}{cc} f_1 \\ f_2 \\ f_3 \\ f_4 \end{array} \right) \end{aligned}$$

Where, $C$ is a 4 × 4 matrix that represents the direct and indirect cross-coupling of the basic waves. $g(N)$ denotes the optical gain of the active layer, which is positive inside the gain section and negative inside the loss section ($g(N)=\frac{g_{max}(N-N_{tr})}{N+[g_{max}/(-g_0)]N_{tr}}\frac{1}{1+\epsilon U}$). $n_g$ is the group index of the guided mode, $\Gamma$ is the optical confinement factor inside the active layer, $\delta$ is the deviation from the Bragg condition due to the carrier-induced refractive index change inside the active layer ($\delta=\delta_0+\frac{2\pi}{\lambda}\Gamma\Delta n$), $\alpha_{in}$ is the material loss, $\gamma$ is the rate of the refractive index change (whose effect is minor) and $f_i$ ($i$ = 1–4) is random noise due to spontaneous emission. (Ref. 6)

Conclusion

In this post, I organized the algorithm steps and briefly described the proof of 3D-CWT. However, if you need to have a detailed understanding of the derivation of CWT, Ref. 3 and Ref. 5 are sufficiently detailed. If you need the relevant code for 3D-CWT, please feel free to contact me.

Reference

1. Coupled-wave model for square-lattice two-dimensional photonic crystal with transverse-electric-like mode | Applied Physics Letters | AIP Publishing
2. Three-dimensional coupled-wave model for square-lattice photonic crystal lasers with transverse electric polarization: A general approach
3. Three-dimensional coupled-wave theory for photonic-crystal surface-emitting lasers
4. General recipe to realize photonic-crystal surface-emitting lasers with 100-W-to-1-kW single-mode operation | Nature Communications
5. Comprehensive analysis of photonic-crystal surface-emitting lasers via time-dependent three-dimensional coupled-wave theory
6. Photonic-crystal lasers with two-dimensionally arranged gain and loss sections for high-peak-power short-pulse operation