## Abstract

Theory and simulations employing coupled-cavity rate equations in an active photonic lattice show excitation of coherent low frequency collective oscillations in the photon and carrier densities, analogous to photonic sound. For parameters just below the lattice stability threshold, long range wave propagation results from external excitation of a few cavities. Above threshold long range coherent oscillations are self-excited without external stirring.

©2004 Optical Society of America

## 1. Introduction

Closely packed microlaser arrays, such as guided mode VCSEL cavities and active defect superlattices, interact through their evanescent fields, Fig. 1. Fringe field interference during stimulated emission introduces active cavity coupling, whereby photons confined in one cavity induce emission in neighboring cavities. Thusly coupled cavities lead to *active photonic lattice* behavior[1] that is similar in certain aspects to solid crystal, whereby the radiation envelopes play the role of atomic wavefunctions.

Specifically a weak interaction among eigenmodes that are laterally confined at each site, implying either paraxial *k*
_{⊥}≡*k*
_{x,y}≪*k _{z}* or laterally evanescent

*k*

_{x,y}→

*iκ*

_{x.y}modes, generates propagating pass-bands

*ω*=

*ω*(

_{o}*k*)+Ω(

_{z}*K*

_{⊥}) where

*ω*is the unperturbed eigenfrequency and Ω is the slow frequency of the modulation envelope with lattice wavevector

_{o}*K*

_{⊥}. (By contrast, interference among radiating sites

*k*

_{⊥}~

*k*creates photonic bandgaps inside a continuum[2].) The situation may arise by the clustering of microlaser cavities characterized by index- or gain-guided stand-alone modes. Similar mechanism applies to fringe interactions among active defect-localized modes[3] in a photonic lattice; in that case the pass-band lies inside the lattice band-gap. One may also consider a passive dielectric photonic crystal and load each site with active material emitting at a bang-gap frequency; an evanescent mode is excited in each site and the lattice turns into a coupled defect lattice.

_{z}Active coupling differs from passive dielectric interference[2] in that it involves stimulated cross-cavity emission and absorption via photon “tunneling”. It is therefore accompanied by a periodic modulation of the carrier density over the cavity lattice. The slow coupled oscillations in the photon and carrier number constitute a “photonic sound” excited in the lattice. The role of the traditional collision-driven pressure is taken over here by stimulated photon absorption and emission. Radiation pressure is negligible because the electron recoil during emission is very small. The observed group velocity, of the order *µm*/*ps*, is comparable to the sound velocity. Near the stability boundaries of the dispersion relation the decay constant of these waves is small and long range propagation effects appear.

The coupled equations for an active photonic lattice, the spontaneous lattice relaxation into Bloch steady-states, and the stability analysis for the lattice oscillations around steady-state has been presented earlier[1]. In this paper we focus on the numerical observation of photonic sound waves over active photonic lattices. In the stable parameter regime these waves are excited by a periodic external driving of selected array sites around the steady-state values. Sizable 1-D and 2-D arrays (up to 1×73 and 21×21) are simulated with either periodic or finite boundaries. For parameters well inside the stable regime the excited waves are decaying over few sites; they become long range waves reaching the system boundaries for parameters near the stability boundary. Above the stability boundary we observe self-excitation of non-linear waves without any external stirring and under uniform in space and time applied laser biasing. They appear as photonic “convection cells”, superficially resembling convection patterns in uniformly heated layers. Chaotic lattice patterns emerge without stirring at even higher bias.

We consider 1-D or 2-D slab array, Fig. 1, with a periodic *xy* (lateral)plane arrangement of laser cavities centered at lattice points **R**
* _{ij}*=

*i*

**b**+

_{x}*j*

**b**

_{y}. Each stand-alone cavity is characterized by profiles of complex gain coefficient

*g*, mirror losses

*µ*, dielectric constant

*∊*and passive absorption

*α*of the general form

*ν*(

**r**)=

*ν*+

_{o}*δνχ*(

_{ν}**r**) where

*ν*stands for

*g*,

*µ*,

*∊*,

*α*respectively and

*ν*denotes a uniform background. The eigenmodes of the uncoupled cavities are

_{o}*U e*

^{ikz-i}^{ω}

^{t}where

*U*(

**r**) is the transverse mode profile and

*z*the axial direction. The material property variations over the lattice have the general form

*ν*(

**r**)=

*ν*+

_{o}*δν*∑

_{ij}*χν*(

**r**-

**R**); the ij-th site experiences the other sites as a perturbation of the form Δ

_{ij}*ν*(

_{ij}**r**)=

*δν*∑

_{i′≠i, j′≠j}

*χν*(

**r**-

**R**

**). For weakly interacting cavities the collective eigenmodes are written in the tight-binding approximation as**

_{i}^{′}**j**^{′}*ε _{ij}* being the ij-cavity complex amplitude. We will focus on the evolution of the discrete lattice envelope

*ε*(

_{ij}*t*), factoring out the fact dependence

*e*-

^{ikz}*. By substituting the form (1) inside the Maxwell-Bloch equations, averaging over axial length*

^{iωt}*z*, and taking the integral projections with

*U*around each lattice site, the coupled lattice equations retaining near neighbor interactions

*i*

^{′}=

*i*±1,

*j*

^{′}=

*j*±1 for the carrier density

*N*, the photon density ℵ and the phase difference Φ

_{ij;i′j′}≡

*φ*

*-*

_{ij}*φ*

*i*

^{′}

*j*

^{′}among cavities are[1]

$$-{\upsilon}_{g}{g}_{r}\mathit{ln}{\hat{N}}_{\mathit{ij}}2{\Upsilon}_{g}\sum _{i\text{'},j\text{'}}\sqrt{{\aleph}_{\mathit{ij}}{\aleph}_{i\prime j\prime}}\mathit{cos}{\Phi}_{\mathit{ij};i\prime j\prime}$$

$$+\sum _{i\prime ,j\prime}{\upsilon}_{g}\phantom{\rule{.2em}{0ex}}2\phantom{\rule{.2em}{0ex}}{Z}_{\mathit{ij};i\prime j\prime}\sqrt{{\aleph}_{\mathit{ij}}{\aleph}_{i\prime j\prime}}\mathit{cos}{\Psi}_{\mathit{ij};i\prime j\prime}$$

$$+\sum _{i\prime ,j\prime}\frac{{\upsilon}_{g}}{2}2\left\{{Z}_{ij;i\prime j\prime}\sqrt{\frac{{\aleph}_{i\prime j\prime}}{{\aleph}_{\mathit{ij}}}}\mathit{sin}{\Psi}_{\mathit{ij};i\prime j\prime}-{Z}_{\mathit{ij}-1;i\prime j\prime -1}\sqrt{\frac{{\aleph}_{i\prime j\prime -1}}{{\aleph}_{i,j-1}}}\mathit{sin}{\Psi}_{i,j-1;i\prime ,j\prime -1}\right\}$$

Here *g _{r}*,

*gi*are the real and imaginary parts of the gain coefficient, ℵ is given in terms of the original electric field amplitude

*ε*=

*Ee*as ℵ=

^{iφ}*∊εε**/8

*π.ħω*and

*N*̂

*≡*

_{ij}*N*/

_{ij}*N*with

_{tr}*N*the transparency density. The cavity coupling strength enters via $Z\equiv \sqrt{{Z}^{2}+{\Xi}^{2}}$ and the combined phase Ψ

_{tr}_{ij;i′j′}=Φ

_{ij;i′j′}+

*ϑ*where

_{ij}*ϑ*=

*tan*

^{-1}Ξ/

*Z*and

*Z*and Ξ given by

involve, respectively, the lump-summed contributions from the imaginary dielectric (real gain *g _{r}*, absorption

*δα*, mirror losses

*δµ*) and the real dielectric (imaginary gain

*g*and refractive

_{i}*δ∊*) response. The coefficients ϒ

*, ϒ*

_{g}*, ϒ*

_{µ}*and ϒ*

_{α}_{∊}express thintercavity overlap strengths among cavity eigenmode profiles

*U*(

**r**-

**R**). These modes are lattice orthonormal ∫

_{ij}**dr**

^{2}*U*(

**r**-

**R**)

_{ij}*U*(

**r**-

**R**

_{i′j′})=∫

**dr**

^{2}*U*(

**r**)

*U*(

**r**±Δ

*R*

_{i-i′,j-j′)}=

*δ*

_{i-i′,j-j′}thus the near-neighbor Δ

**R**

_{i-i′,j-j′}=±

*b*coupling strengths

*Y*≡∫

_{a}**dr**

*U*(

**r**)

*χa*(

**r**)

*U*(

**r**±

**b**) stem entirely from the variations in the gain, mirror reflectivity

*δµ*, absorption

*δα*, and dielectric constant

*δ∊*from their respective “floor” values. The different strengths of

*Yν*,

*ν*=

*g*,

*µ*,

*∊*,

*α*reflect the differences in the spatial profiles

*χν*(

**r**); for identical profiles

*χ*all coupling terms assume the same value ϒ. Finally the term

*S*

_{ij;i′j′}=Λ

_{g}*ζg*̂

_{r}lnN_{i′j′}+Λ

_{µ}*µ*-Λ

_{∊}*δα*̂ lump-sums the phase-independent contributions Λ

*≡∫*

_{ν}**dr**

^{2}*U*(

**r**±

**b**)

*χν*(

**r**)

*U*(

**r**±

**b**) from inter-cavity interactions. Notice that the carrier density evolution involves only real gain-coupling terms. For small coupling strengths ϒ

*≪1, and for around steady-state variations, the carrier and photon densities scale as*

_{ν}*N*≃

*N*+𝒪(ϒ), ℵ≃

_{o}*ℵ*+𝒪(ϒ), where the subscript

_{o}*means the steady-state values for an isolated-cavity; for linear oscillations we may approximate*

_{o}*Z*

_{ij;i′j′}≃

*Z*and Ξ

_{o}_{ij;i′j′}Ξ

*with ℵ*

_{o}*inside Eqs. (5)–(6).*

_{o}The collective array behavior has been investigated via numerical integration of the coupled lattice Eqs. (2)–(4), starting with random initial conditions. Under lattice-periodic *J*
* _{ij}*(

**r**)=

*J*(

_{o}χJ*r*), constant in time drive currents we observe a spontaneous relaxation to steady-states

*dN*/

*dt*=

*dℵ*/

*dt*=0. For periodic arrays the lattice envelope assumes a Bloch wave form, where all amplitudes are constant

*N*=

_{ij}*N*,

_{o}*ℵ*=

_{ij}*ℵ*and the phase-difference among all neighbor sites is fixed Φ

_{o}^{x,y}

_{ij;i′j′}≡

*φ*-

_{ij}*φ*

_{i′j′}=

**K**

*·*

_{mn}*b*

_{x,y}and proportional to an inverse lattice vector

*K*=(

_{mn}*m*,

*n*)2

*π*/

*L*. In general

Any value Φ yields a steady state for the rate Eqs. (2)–(4) but configurations maximizing the collective gain are favored; thus Φ is either zero or *π* for *Z* positive or negative. The steady-state values *N _{o}*,

*ℵ*are obtained analytically from the zeros of Eqs. (2)–(4).

_{o}We proceed to examine the collective evolution of small perturbations about given steady state values Φ* _{ij}*=Φ

*,*

_{o}*ε*=

_{ij}*ε*,

_{o}*N*=

_{ij}*N*. Collective modes over a periodic lattice conform with the discretized version of a Bloch wave sampled over lattice sites

_{o}*δA* being the lattice envelope of either carrier or photon density. Similar variations in the intercavity phase shift, of amplitude *δ*Φ* _{o}*, are caused by out-of-phase electric field changes

*δε*relative to the steady state field

*ε*. Replacing independent phase variations by a lattice-periodic envelope phase manifests the transition to collective behavior with long range coherence (in general situations far from equilibrium, Φ

_{o}*in Eqs. (2)–(4) are independent variables). The description resembles the dynamic behavior of a periodic, spring-coupled mass system. Due to the periodicity, linearization of Eqs. (2)–(4) about any site leads to identical stability equations. Expanding around equilibrium and using Eq. (8) yields the following stability matrix for perturbations about steady-state*

_{ij}The elements *D _{XY}*≡

*∂*(

*X*Ẋ)/

*∂Y*with

*X*,

*Y*being either of

*N*,

*ε*, Φ are found from the rhs of Eqs. (2)–(4) at steady-state - the linearization is performed using complex E notation and the results are then recast in terms of ℵ∝|

*ε*|

^{2}. The roots of the cubic characteristic equation

yield the dispersion relation for lattice perturbation wavenumber *κ*, affected by the perturbed steady state period **K**. It has been shown that the real (non-oscillating) root is negative (stable) λ_{3}<0 over the entire parameter space of interest. Lattice stability is thus determined by the Γ_{1,2} sign of the complex root pair λ_{1,2}. Linear stability depends primarily on the ratio *g _{i}*/

*g*between the imaginary and real gain coefficients[1]. Unconditional stability exists for

_{r}*g*/

_{i}*g*≤1, yielding a dispersion of stable, time decaying collective oscillations. For

_{r}*g*/

_{i}*g*>1 the stability is conditional and depends on the coupling strength

_{r}*Ƶ*. At larger coupling strengths we observe transition to limit cycles, and eventually chaotic oscillations. Finally, for

*g*/

_{i}*g*≪1, the complex roots are approximated by the eigenvalues of the upper 2×2 submatrix.

_{r}Dispersion Eq. (10) represents complex frequency *ω* for a real wavenumber *κ*. It applies to initial value problems, regarding the time evolution of periodic lattice perturbations. Equivalently, complex wavenumber solutions *κ*→*κ*+*iξ* exist for real excited frequencies *ω*. That applies to boundary value problems, related to the excitation of selected lattice sites with real frequency *ω _{o}*, and subsequent wave propagation. Near the boundary of the stability regime, where Γ≪Ω, imposing ℑ

*ω*=

*ℑ*(λ)=0 to the analytic continuation of Eq. (10) with complex κ yields the decay constant ξ,

where *V _{gr}* is the lattice group velocity. Since ξ is proportional to the decay rate Γ, when the lattice parameters are well inside the stability region the decay constant is fairly large,

*ξb*

_{x,y}~1, and driven waves decay over few lattice sites. Figure 2 shows various snapshots of the photon density envelope

*ℵ*obtained during a simulation of 21×21 finite VCSEL arrays. Three central sites

_{ij}*i*=10,

*j*=10-12 are driven by a sinusoidal drive current of amplitude

*I*/

_{dr}*I*=1.1 superimposed on the steady-state drive

_{th}*I*/

_{o}*I*=3.1 for the entire array. The coupling strength in all cases is

_{th}*Z*=0.0066 and yields in-phase

*K*

_{x,y}

*b*

_{x,y}=0 steady-state. In Fig. 2a (46ns) the ratio

*g*/

_{i}*g*=0.5 lies well inside the stability regime and yields evanescent waves within two lattice sites, regardless of the drive frequency

_{r}*ω*. Increasing the ratio

_{o}*g*/

_{i}*g*=1.5 near the stability boundary excites small

_{r}*ξb*

_{x,y}≪1, long range waves propagating over the entire array, Figs. 2b (44ns) and 2c (37ns). The wave periods

*κ*

_{x,y}

*b*

_{x,y}=

*π*/2 and

*κ*

_{x,y}

*b*

_{x,y}=

*π*respectively match the drive current frequencies

*ω*=Ω(

_{o}*κ*

_{x,y};

*K*=0) from Eq. (10). The wave decay is in large part due to the 1/

*r*dependence of the “point”-source radiation, due to energy conservation.

Figure 3 plots the dynamic in time evolution of a 1×49 array excited at *j*=25, *g _{i}*/

*g*=1.5 and same drive currents as before. The coupling strength

_{r}*Z*=0.0133 corresponds to a site separation-to-radiation waist ratio

*b*/

*ω*≃2. The excited wave pattern of

*κb*=

*π*/2 (four sites per wavelength) matches the applied drive frequency

*ω*=Ω(

_{o}*κ*). The amplitude decay rate in 1D is due solely to the decay constant γ value. For typical site separation

*b*=4

*µm*, the group velocity estimated from the Mach cone is

*υ*≃2.4×10

_{gr}^{4}

*cm*/

*s*, of the order of the sound velocity in solids.

When a lattice parameter crosses the stability boundary there is no steady-state relaxation. Instead a limit cycle emerges where the lattice undergoes a periodic sequence of coherent non-uniform patterns, without external stirring. Figure 4 shows two snap-shots of the electric field envelope for a finite 21×21 array with *g _{i}*/

*g*=4 and

_{r}*Ƶ*=0.00155. Surface of section plots

*N*(

_{ij}*t*) vs.

_{n}*ℵ*(

_{ij}*t*) at a given site, recorded at time increments

_{n}*t*, exhibit periodic non-linear cycles. The trajectories for all lattice sites fall into the

_{n}*same*curve, each site separated by a phase shift around the cycle relative to its neighbors[1]. The transition from a fixed point

*N*=

_{ij}*N*to a limit cycle is an example of a Hopf bifurcation. For even higher coupling strengths

_{o}*Ƶ*we observe transitions to intermittent and finally fully spatio/temporary chaotic states, under constant and uniform external biasing.

## References and links

**1. **S. Riyopoulos, “Coherent phase locking, collective oscillations and stability in coupled VCSEL arrays,” Phys. Rev. A **66**, 53820 (2002). [CrossRef]

**2. **J. D. Joannopoulos, R. D. Meade, and J. N. Winn in *Photonic Crystals* (Princeton University, Princeton N.J., 1995).

**3. **O. Painter, J. Vuckovic, and A. Scherer. “Defect Modes of a two-dimensional photonic crystal in an optically thin dielectric slab,” J. Opt. Soc. Am. B **16**, 275 (1999). [CrossRef]

**4. **A. Yariv, Y. Xu, R. K. Lee, and A. Scherer, “Coupled resonator optical waveguide: a proposal and analysis,” Opt. Lett. **24**, 711–713 (1999). [CrossRef]