# Large Scale Constrained Clustering With Reinforcement Learning

Benedikt Schesch,1,2 Marco Caserta1

###### Abstract

Given a network, allocating resources at clusters level, rather than at each
node, enhances efficiency in resource allocation and usage. In this paper, we
study the problem of finding fully connected disjoint clusters to minimize the
intra-cluster distances and maximize the number of nodes assigned to the
clusters, while also ensuring that no two nodes within a cluster exceed a
threshold distance. While the problem can easily be formulated using a binary
linear model, traditional combinatorial optimization solvers struggle when
dealing with large-scale instances. We propose an approach to solve this
constrained clustering problem via reinforcement learning. Our method involves
training an agent to generate both feasible and (near) optimal solutions. The
agent learns problem-specific heuristics, tailored to the instances
encountered in this task. In the results section, we show that our algorithm
finds near optimal solutions, even for large scale instances.

## Introduction

In Amazon’s operations, resources are traditionally allocated on a site-by-
site basis. Shifting to a cluster-based model enhances resource distribution
efficiency, such as assigning a single technician with rare skills to multiple
sites, reducing travel distances, and capping travel times within clusters.
This requires the development of an approach that identifies the best clusters
in the network, while respecting unique business constraints.

Mixed-integer solvers, traditionally used for these types of constrained
clustering problems, struggle with large-scale instances, due to the NP-hard
nature of these problems. Recent trends (Bengio, Lodi, and Prouvost 2018;
Bello et al. 2016) involve integrating machine learning, particularly
reinforcement learning (RL), into a combinatorial optimization problem. Our
RL-based approach trains agents to develop specific heuristics, enabling them
to handle large-scale instances effectively.

In the sequel, we present a binary linear formulation for the problem; next,
we provide a description of our method; a computational section showcases the
effectiveness of the method and, finally, we conclude with some remarks.

## Problem Formulation

Our formulation is loosely based on the modularity maximization problem.
(Agarwal and Kempe 2008) Given a set of $n$ sites and their inter-site travel
times $d_{ij}$, we minimize the total intra-cluster travel times, while
penalizing unclustered sites. In addition, the maximal distance between any
two sites assigned to the same cluster must be below a threshold $D$. The
decision variable $x_{ij}$ takes value $0$ if sites $i$ and $j$ are assigned
to the same cluster, and $1$ otherwise.

$\displaystyle\min\quad\sum_{i=1}^{n}\sum_{j=i}^{n}\left((1-x_{ij})d_{ij}+x_{ij}D)\right)\mathds{1}_{d_{ij}<D}$
$\displaystyle\textrm{s.t.}\quad\begin{aligned} &x_{ij}=1\text{ if }d_{ij}\geq
D&i,j\in 1\ldots n\\\ &x_{ij}\leq x_{il}+x_{lj}&i,j,l\in 1\ldots n\\\
&x_{ij}=x_{ji}&i,j\in 1\ldots n\end{aligned}$

In the model, the objective function only considers pair of sites whose
distance is below the threshold $D$, and accounts for the total intra-cluster
distances and the penalty for the unclustered sites. The triangle inequality
ensures that each cluster is a fully connected graph, while the symmetry
constraints reflect the undirected nature of the clusters.

## Solution methodology

### RL Environment.

We detail the RL formulation:

State Space. Each state is a triplet, including the distance matrix, the edges
in the current solution, and the remaining available edges.

Action Space. The agent selects from the available edges.

Transition Function. After choosing an edge, the edge joins the current
solution and is removed from the set of available edges. Further edges needed
to satisfy the triangle inequality are added, while edges violating this
inequality are removed from the available set.

Reward. The reward for an action is the difference in the objective function
resulting from the transition.

Figure 1 depicts an example with 4 cities. Blue and red edges are the
available edged and the edges in the current solution, respectively.
Initially, the agent selects edge (2, 3). Consequently, edge (3, 4) is removed
to respect the triangle inequality. In the second step, the agent picks the
edge (1, 3) which is added to the solution. The edge (1, 2) is also added to
respect the triangle inequality. Since there are no available edges anymore, a
solution has been constructed, with objective function equal to 3 + 5 + 6 + 60
= 74, which is an optimal solution.

Figure 1: Example instance with threshold $D=60$.

This approach allows for easy constraints integration by altering the
transition function. It involves removing edges that breach the constraints,
with the advantage of not requiring formal constraints expression, which can
be complex.

Agent. Our agent, designed to predict the value function and the next edge
selection, leverages a graph neural network (GNN) due to the graph-centric
nature of our problem. GNNs, recognized for their efficacy in handling graph-
structured data (Wu et al. 2019), are well-suited here. In our model, nodes
represent sites and edges denote either available or solution-participating
connections. Each edge is characterized by the internode distance and two
binary indicators signifying its availability and inclusion in the solution.
Each node has a constant feature as a bias term.

The initial feature embedding employs a learned linear transformation. The
network architecture has three EGAT convolutions layers (Kamiński et al.
2022). Action readout uses a linear transformation of the final embedding for
each available edge, with all embeddings sized at 8. This process assigns
scores to edges, guiding the agent’s preference. For the critic, we use a
dense neural network, processing the flattened distance and availability
matrices.

Training. We train our agent in our environment using Proximal Policy
Optimization (PPO) (Schulman et al. 2017) with a learning rate of $10^{-5}$.

## Computational Experiments

We test our approach on two experimental setups, using different distributions
over the distance matrices, to verify that the agent learns distribution
specific heuristics.

Cities Environment. To mimic the natural distribution of cities, we allocate
sites among three major cities on a 2D map ($[0,240]\times[0,240]$). Half of
the site are close to the first city, a third to the second, and a sixth to
the smallest city. Site locations are normally distributed around city
centers, with variances from $\mathcal{U}(80,160)$.

We trained and tested our agent on 128 instances, comparing its performance to
that of a random agent. The results for 18-site and 60-site instances are
shown in Tables 1. Our agent consistently outperformed the random agent,
achieving a median optimality gap of $0\%$ and $1.02\%$ in the 18-site and the
60-site instances, respectively, with a substantial improvement in the quality
of clustering solutions.

| 18 sites | 60 sites
---|---|---
|  Random |  Ours |  Random |  Ours
Mean opt. gap | $16.7\%$ | 0.02% | $31.6\%$ | 1.02%
Median opt. gap | $15.8\%$ | 0.00% | $31.0\%$ | 0.43%
Min opt. gap | $0.00\%$ | $0.00\%$ | $14.3\%$ | 0.00%
Max opt. gap | $50.2\%$ | 1.45% | $59.0\%$ | 5.11%
Table 1: Optimality gaps over 128 instances for the city-like design.

General Environment. In this scenario, the distance matrix is generated with
entries following a $\mathcal{U}(0,240)$ distribution. This is a more
challenging problem than the cities environment due to its higher
dimensionality. We focus on the most difficult instances, for which the Xpress
Solver requires over 2000 seconds to terminate with an optimal solution,
comprising a dataset of 69 examples with an average runtime of 2276 seconds.
We also created a second model, with embedding size equal to 64.

|  Random |  Ours-8 |  Ours-64
---|---|---|---
Mean opt. gap | $5.3\%$ | $3.6\%$ | $\mathbf{3.0}\%$
Median opt gap | $5.3\%$ | $3.6\%$ | $\mathbf{3.1}\%$
Min opt gap | $3.9\%$ | $2.7\%$ | $\mathbf{1.8}\%$
Max opt gap | $6.9\%$ | $4.5\%$ | $\mathbf{4.2}\%$
Table 2: Optimality gaps on 69 “difficult” instances with 64 sites.

Table 2 provides the results. Our agent consistently outperforms the random
agent in terms of optimality gap, though it does not always reach optimal
solutions. The larger agent version shows some improvements, indicating that
scaling and further sophistication in the model are beneficial. Despite not
achieving optimal results, our agent provides a fast approximation, especially
in comparison to the running time of an off-the-shelf solver. Advancements in
the graph neural network structure, beyond current EGAT convolutions, could
potentially enhance its effectiveness and help closing the gap to optimality.

## Conclusion

In this study, we proposed a reinforcement learning-based approach to solve
constrained clustering problems, especially useful to deal with very large
instances. The results demonstrate that this approach provides near optimal
solutions to the constrained clustering problem, with a running time orders of
magnitude below the time required by a solver. The next step in our research
is to improve the network architecture, to further reduce the optimality gap
for difficult instances.

## References

  * Agarwal and Kempe (2008) Agarwal, G.; and Kempe, D. 2008.  Modularity-maximizing graph communities via mathematical programming.  _The European Physical Journal B_ , 66(3): 409–418.
  * Bello et al. (2016) Bello, I.; Pham, H.; Le, Q. V.; Norouzi, M.; and Bengio, S. 2016.  Neural Combinatorial Optimization with Reinforcement Learning.  _CoRR_ , abs/1611.09940.
  * Bengio, Lodi, and Prouvost (2018) Bengio, Y.; Lodi, A.; and Prouvost, A. 2018.  Machine Learning for Combinatorial Optimization: a Methodological Tour d’Horizon.  _CoRR_ , abs/1811.06128.
  * Kamiński et al. (2022) Kamiński, K.; Ludwiczak, J.; Jasiński, M.; Bukala, A.; Madaj, R.; Szczepaniak, K.; and Dunin-Horkawicz, S. 2022.  Rossmann-toolbox: a deep learning-based protocol for the prediction and design of cofactor specificity in Rossmann fold proteins.  _Brief Bioinform_ , 23(1).
  * Schulman et al. (2017) Schulman, J.; Wolski, F.; Dhariwal, P.; Radford, A.; and Klimov, O. 2017.  Proximal Policy Optimization Algorithms.  _CoRR_ , abs/1707.06347.
  * Wu et al. (2019) Wu, Z.; Pan, S.; Chen, F.; Long, G.; Zhang, C.; and Yu, P. S. 2019.  A Comprehensive Survey on Graph Neural Networks.  _CoRR_ , abs/1901.00596.


# Discrete Coordinate Scattering Approach to Nano-optical Resonance

Jun Yang  Department of Physics, Peking University, 100871, Beijing, P.R.
China  H. Dong<EMAIL_ADDRESS>Institute of Theoretical Physics, Chinese
Academy of Sciences, 100190, Beijing, P.R. China  C.P. Sun<EMAIL_ADDRESS>Institute of Theoretical Physics, Chinese Academy of Sciences, 100190,
Beijing, P.R. China

###### Abstract

In this letter, we investigate the coherent tunneling process of photons
between a defected circular resonator and a waveguide based on the recently
developed discrete coordinate scattering methods (L. Zhou et al., Phys. Rev.
Lett. 101, 100501 (2008)). We show the detailed microscopic mechanism of the
tunneling and present a simple model for defect coupling in the resonator. The
Finite-Difference Time-Domain(FDTD) numerical results is explored to
illustrate the analysis results.

Recently, resonant optical cavities has attracted lots of interesting from
both experimental and theoretical aspect vahla03review . Devices based on
these micro-cavities has been widely used for processing optical signals
joannopoulos . In most hardware of these signal processing, the elementary
unit consists of a micro-cavity and a side coupling waveguide vahla03review .
This unit has been experimentally realized based on photonic crystal, most
recently ring made with GeN nanowire pdyang06prl . With such kind photonic
setups, photon blockade was experimentally observed in this unit kimble08 .
For this unit, it is well known that the incident wave is totally reflected by
the ring resonator when it is resonant with the modes of the ring. Thus, the
resonator takes the role of frequency-selection little98opt  in this unit.

Such transfer process in the unit is crucial to design more complicated
devices. To investigate this phenomena, one way is to directly solve the
Maxwell equations. However, it is impossible to get an analytic result to
reveal the key physical mechanism for even a simple configuration. Thus, this
is usually done by numerical simulation, such as Finite-Differential-Time-
Domain(FDTD) method taflove , which is used as a validation approach.
Theoretically, coupled mode theory(CMT) is introduced to study the resonator
scattering properties in this unit little98lightware ; little98opt . However,
this is done in the mode space, and therefore not intuitive in real space.

In this letter, we generalize the most recently developed discrete coordinate
scattering methods lan08  to investigate the detailed scattering properties
for the optical signal processing. Here we model the resonator and the
waveguide as array of coupled cavities. And the coupling between resonator and
waveguide is modeled as hopping of the photon among cavities. We present the
method by a qualitatively analysis of an example of a notched-ring resonator.
Rigorous FDTD numerical simulations is also performed to verify the obtained
analytical result.

Figure 1: (Color online) Frequency-selective filter. (a) A notched ring is
evanescently coupled to the waveguide. The notch couples the CW and CCW
traveling modes; (b) Model of the unit by discrete coupled resonator array.
The notch is characterized as point inhomogeneous of the tunneling strength
along the ring.

The elementary units under investigation is schematically illustrated in Fig.
1(a). The defected ring resonator is evanescently coupled to the signal
waveguide. For perfect ring, it supports the clockwise(CW) and counter
clockwise(CCW) traveling wave and there is no interaction between the two
traveling wave modes. However, the scattering of the defected notch causes the
recombination of the CW and CCW modes. In the context of CMT little98opt , the
notched ring is mapped to the coupled double-ring resonator and the notch in
single ring plays the role of the gap, which couples the double rings. Here,
we treat the defected circular resonator as the discrete coupled cavity array
with coupling constant different at the defected area, which is illustrated in
Fig. 1(b). The array consists of $N$ resonators with tight-binding coupling.
The notch in the ring is marked with a dashed line. The Hamiltonian here reads

$H_{c}=\sum_{j}\omega
a_{j}^{\dagger}a_{j}-J\left(a_{j}^{\dagger}a_{j+1}+\text{h.c.}\right)+\delta
J\left(a_{0}^{\dagger}a_{1}+\text{h.c.}\right),$ (1)

where $a_{j}$ is the creation operator for photon with energy $\omega$ the
$j-th$ cavity and $J$ is the tunneling strength between nearby cavities with
defected induced inhomogeneous coupling $\delta J$. By the Fourier
transformation
$a_{j}^{\dagger}=\sum\exp\left(-ijk\right)c_{k}^{\dagger}/\sqrt{N}$, the
Hamiltonian in the momentum space reads $H_{c}=\sum_{k}\left(\omega-2J\cos
k\right)c_{k}^{\dagger}c_{k}+\delta
J\sum_{k,p}\left[\exp\left(ip\right)c_{k}^{\dagger}c_{p}+h.c.\right]/N$. In
the rotation wave approximation, we keep only $p=\pm k$ term and obtain the
Hamiltonian as
$H_{c}=\sum_{k>0}\varepsilon_{k}\left(c_{k}^{\dagger}c_{k}+c_{-k}^{\dagger}c_{-k}\right)+\sum_{k>0}2\delta
J\left(\exp\left(-ik\right)c_{k}^{\dagger}c_{-k}+\text{h.c.}\right)/N$, where
$\varepsilon_{k}=\omega-2\left(J-\delta J/N\right)\cos k$. Physically,
$c_{k}$and $c_{-k}$ stand for the CW and CCW traveling wave, which has been
introduced phenomenally in CMT. The coupling between the two modes is induced
by the defect of the inhomogeneous of tunneling $\delta J$. Diagonalizing the
Hamiltonian with unitary transformation
$c_{k}^{\dagger}=\exp\left(ik/2\right)(\alpha_{k}^{\dagger}+\beta_{k}^{\dagger})\sqrt{2}$
and
$c_{-k}^{\dagger}=\exp\left(-ik/2\right)(\alpha_{k}^{\dagger}-\beta_{k}^{\dagger})\sqrt{2}$
for $k>0$, we obtain

$H_{c}=\sum_{k>0}\varepsilon_{k}^{+}\alpha_{k}^{\dagger}\alpha_{k}+\varepsilon_{k}^{-}\beta_{k}^{\dagger}\beta_{k},$
(2)

where $\varepsilon_{k}^{\pm}=\left(\varepsilon_{k}\pm 2\delta J/N\right)$ and

$\displaystyle\alpha_{k}^{\dagger}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2N}}\sum_{j}2a_{j}^{\dagger}\cos\left[(j-1/2)k\right]$
(3) $\displaystyle\beta_{k}^{\dagger}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2N}}\sum_{j}2ia_{j}^{\dagger}\sin\left[(j-1/2)k\right].$
(4)

are the eigen-modes of the notched ring. And they were related to photonic
molecule in Ref. bayer98prl ; pdyang06prl . Actually, this single defected
notch induced the split of the degenerate resonant peak of CW and CCW
traveling wave. The splitting can be probed by checking the transmission rate
in the waveguide, which will be investigated in the following parts. Actually,
those modes correspond to a symmetric and antisymmetric state with respect to
the notch. The amplitude for the modes $\alpha_{k}$ and $\beta_{k}$ is
illustrated in Fig. 2. In this discrete model of the ring resonator, the
position $j_{\mathrm{notch}}$ of notch should between the resonator $a_{0}$
and $a_{1}$, namely $j_{\mathrm{notch}}=1/2$. Thus for the $\alpha$ mode, the
amplitude of the field at the notch is at maximum. However, there is no field
in the notch for the $\beta_{k}$ mode.

Figure 2: Ez-field amplitude of $\alpha_{k}$(Blue line) and $\beta_{k}$(Red
dashed line) mode on the notched ring. The position of notch is marked with
dot line. For $\alpha_{k}$ mode, the amplitude of Ez-field at the notch is at
maximum. However, the amplitude is null for $\beta_{k}$ mode.

Figure 3: Plot of the reflection rate with the parameters
$\omega=\omega_{1}=8.772\times 10^{14}Hz$, $J=J^{\prime}=3.3\times 10^{14}Hz$,
$g=3\times 10^{13}Hz$ and $\delta J=2.98\times 10^{13}Hz$. (a) Theoretical
calculation of the reflection of the filter; (b) FDTD simulation of the
reflection and transition rate. Two resonator peaks appear at the cavity mode
$\varepsilon_{k}^{\pm}$.

To probe the optical property of the resonator, a transmission line waveguide
is introduced to couple the signal in the microcavity. Physically, the
waveguide can be optical fiber, nanowire and defected lines on photonic
crystal for different systems. Here, we also model the transmission line
waveguide as the discrete coupled resonator, which is illustrated in Fig.
1(b). Each resonator couples to the near-by resonators by means of photon
hopping. Those coupled-resonator array has been proposed to be a new type of
optical waveguideyariv99opl . The Hamiltonian for the transmission line
waveguide here reads

$H_{l}=-J^{\prime}\sum\left(b_{j}^{\dagger}b_{j+1}+\text{h.c.}\right)+\omega^{\prime}b_{j}^{\dagger}b_{j},$
(5)

where $b_{j}^{\dagger}$ is the bosonic creation operator for photon on the
$j-th$ resonator. The tight-binding Hamiltonian here is a good approximation
for modeling the coupled resonator array, which is proved in Ref. yariv99opl .
And in Ref. lan08 , it is showed the model is equivalent to the simple
dielectric waveguide under the linear dispersion area with right going wave
$\phi_{R}(x)$ and left going wave $\phi_{L}(x)$ as showed in Ref. fan98 . The
waveguide couples to the cavity at the $0$-th resonator as the tunneling
between 0-th resonator and the $j_{0}$-th resonator of the ring resonator
$H_{I}=g\left(a_{j_{0}}^{\dagger}b_{0}+a_{j_{0}}b_{0}^{\dagger}\right)$. The
total Hamiltonian for this elementary unit reads $H=H_{c}+H_{l}+H_{I}$. For
the single photon case, the eigenwave stationary scattering function can be
written as

$\left|\Phi\right\rangle=\sum_{j}u_{j}b_{j}^{\dagger}\left|0\right\rangle+\sum_{k>0}\left(v_{k}\alpha_{k}^{\dagger}+w_{k}\beta_{k}^{\dagger}\right)\left|0\right\rangle,$
(6)

where $u_{j}$ is the amplitude of photon in the $j-th$ resonator on the
waveguide and $v_{j}$($w_{j}$) is the amplitude of photon on the resonator
with mode $\alpha^{\dagger}$($\beta^{\dagger}$). Substituting the stationary
scattering function into the time-independent Schrödinger equation
$H\left|\Phi\right\rangle=E\left|\Phi\right\rangle$, one can obtain equations
for the amplitude of the site as

$\displaystyle-J^{\prime}\left(u_{j+1}+u_{j-1}\right)+\omega^{\prime}u_{j}+\frac{g\delta_{j,0}}{\sqrt{2N}}\sum_{k>0}\left(v_{k}+w_{k}\right)e^{-i\left(1/2-j_{0}\right)k}+\left(v_{k}-w_{k}\right)e^{-ik(j_{0}-1/2)}$
$\displaystyle=$ $\displaystyle Eu_{j},$ (7)
$\displaystyle\varepsilon_{k}^{+}v_{k}+\frac{gu_{0}}{\sqrt{2N}}\left(e^{i\left(1/2-j_{0}\right)k}+e^{i(j_{0}-1/2)k}\right)$
$\displaystyle=$ $\displaystyle Ev_{k},$ (8)
$\displaystyle\varepsilon_{k}^{-}w_{k}+\frac{gu_{0}}{\sqrt{2N}}\left(e^{i\left(1/2-j_{0}\right)k}-e^{i(j_{0}-1/2)k}\right)$
$\displaystyle=$ $\displaystyle Ew_{k}.$ (9)

When the incident energy coincides with the energy of $\alpha$ mode, namely
$E=\varepsilon_{k}^{+}$, the amplitude $u_{0}=0$ by Eq. (8) and $w_{k}=0$ by
Eq. (9). Thus only $\alpha$ mode is excited here and the Ez-field at the notch
is at maximum. When the incident energy is resonant with the energy of $\beta$
mode, namely $E=\varepsilon_{k}^{-}$, the amplitude $u_{0}=0$ by Equation(9)
and $v_{k}=0$ by Equation(8). Here, only $\beta$ mode is excited and the
amplitude of the field is zero. This phenomena has been numerically predicted
in Ref. little98opt .

For the scattering process, the photon incidents from the left side of the
line waveguide and is partially reflected by the notched ring resonator array.
The amplitude in the waveguide reads
$u_{j}=\left[\exp\left(iqj\right)+r\exp\left(-iqj\right)\right]\theta\left(-j\right)+s\exp\left(iqj\right)\theta\left(j\right)$
, where $r$($s$) is the reflection(transmission) coefficient of the photon.
The coefficient of the transition here is obtained

$s=\left(1-\frac{Q}{2iJ^{\prime}\sin q}\right)^{-1},$ (10)

where
$Q=2g^{2}/N\sum_{k>0}\cos^{2}(k(j_{0}-1/2))/\left(E-\varepsilon_{k}^{+}\right)+\sin^{2}(k(j_{0}-1/2))/\left(E-\varepsilon_{k}^{-}\right)$.
The reflection rate $1-|s|^{2}$ peaks at the resonance with the internal ring
mode $\varepsilon_{k}^{\pm}$, where the optical signal is totally reflected.
To see this resonant filtering effect, we plot the transmission rate $|s|^{2}$
verse the the incident energy of photon in Fig. 3(a). For comparison, a
rigorous FDTD taflove  simulation is performed with free available code F2P
qiu  for Transverse-Magnetic mode. The simulation is performed with the
following parameter: $n_{0}=1$, $n_{w}=3.2$, $R_{1}=5a$, $R_{2}=4.4a$,
$\theta=54^{\mathrm{o}}$ and $\delta\theta=1.5^{\mathrm{o}}$, where $a=0.25\mu
m$. The transmission and reflection rate is plotted in Fig. 3(b) around the
resonant frequency $\varepsilon_{k}=0.331\times 2\pi c/a$. In Fig. 3(a), the
parameters for this defect ring resonator unit are chosen to fit the numerical
simulation: $\omega=\omega_{1}=8.772\times 10^{14}Hz$, $J=J^{\prime}=3.3\times
10^{14}Hz$, $g=3\times 10^{13}Hz$ and $\delta J=2.98\times 10^{13}Hz$.
Therefore, our model confirms the function of the filter and its mechanism.
However, we should notice a drawback to this approach: the quantitative
analysis will relies on the numerical simulation. At present, the coupling
constant at the notch and between the transmission line and the resonator
should be fitted by numeric simulation.

In conclusion, we have presented a theoretical investigation of the notched
ring filter with discrete coordinate scattering methods. The detailed
mechanism of the coupling between the clockwise and count-clockwise is
revealed to be the inhomogeneous tunneling in the resonator array system. We
have predicted some experimentally accessible results for the coherent
transmission and reflection property and present detailed explanation to field
amplitude distribution at the notch.

This work is supported by NSFC No.10474104, No. 60433050, and No. 10704023,
NFRPCNo. 2006CB921205 and 2005CB724508.

## References

  * (1) K. J. Vahala, Nature 424, 839 (2003).
  * (2) J. D. Joannoupoulos, R. D. Meade and J. N. Winn, Photonic Crystals: Molding the Flow of Light (Princeton U. Press), 1995\.
  * (3) P. J. Pauzauskie, D. J. Sirbuly, and P. Yang, Phys. Rev. Lett. 96, 143903 (2006).
  * (4) B. Dayan, A. S. Parkins, T. Aoki, E. P. Ostby, K. J. Vahala, H. J. Kimble, Science 319, 1062 (2008).
  * (5) B. E. Little, S. T. Chu and H. A. Haus, Opt. Lett. 23, 1570 (1998).
  * (6) A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Methods (Artech House, Boston,2000)
  * (7) B. E. Little, S. T. Chu, H. A. Haus, J. S. Foresi and J.-P. Laine, J. Lightware Technol. 15, 998 (1997).
  * (8) Lan Zhou, Z. R. Gong, Yu-xi Liu, C. P. Sun and F. Nori, Phys. Rev. Lett. 101, 100501 (2008).
  * (9) A. Yariv, Y. Xu, R. K. Lee and A. Scherer, Opt. Lett. 24, 711 (1999).
  * (10) J. T. Shen, S. Fan, Phys. Rev. Lett. 95, 213001 (2005); ibid. 98, 153003 (2007); Opt. Lett. 30, 2001 (2005)..
  * (11) M. Qiu, F2P: Finite-difference time-domain 2D simulator for photonic devices, http://www.imit.kth.se/info/FOFU/PC/F2P/.
  * (12) M. Bayer, T. Gutbrod, J. P. Reithmaier, A. Forchel, T. L. Reinecke, P. A. Knipp, A. A. Dremin, and V. D. Kulakovskii, Phys. Rev. Lett. 81, 2582 (1998).