text:
# Interaction of a soliton with a local defect in a fiber Bragg grating

## 1 Introduction

Optical fibers carrying a Bragg grating (BG) written on their cores have attracted a great deal of attention for various reasons . Fibers equipped with homogeneous or chirped BG are a basis for the design of many elements used in optical telecommunications, such as filters , add-drop WDM multiplexers , dispersion compensators , and others. They also find important applications as sensors, since the optical properties of the fiber-based BG are strongly affected by local stress, see, e.g., a recent work  and reviews . Besides that, BG-carrying fibers are interesting in their own right, as a unique physical medium supporting various optical excitations. In particular, strong dispersion induced by BG opens a way to observe slow light in this medium, which is a topic of great current interest .

The BG-induced dispersion in combination with the Kerr nonlinearity of the fiber give rise to localized pulses, which are usually called gap solitons (or, more generically, BG solitons ). After being predicted analytically  (see also the review ), BG solitons have been observed in experiments . They have attracted a lot of interest as their dynamical properties are very different, in many respects, from the properties of usual solitons in optical fibers (except for the regime in which the amplitude of the BG soliton is small, then it is virtually equivalent to a usual nonlinear-Schrödinger soliton ). In particular, quest for very slow solitons, that are expected to exist in fiber BGs , is a challenge .

The subject of this work is interaction of a soliton with a local defect in the BG-carrying fiber. A local defect in this medium may trap light via four-wave mixing without the formation of a soliton per se , but the interaction of solitons with defects is especially interesting . First, it is a fundamental physical problem – in particular, because the defect may help to stop the soliton and thus to create a self-sustained pulse of standing light (leaving aside the problem of the fiber-loss compensation, which is a subject of another work ). Second, narrow soliton pulses interacting with local inhomogeneities may increase the sensitivity when the grating is used as a sensor. Moreover, the trapping of soliton in a fiber BG may be used to design optical-memory devices. Further potential applications include optical delay lines and oscillators, see the concluding section of this paper.

The objective of the present work is a detailed study of the interaction of BG solitons with a local point-like defect. In particular, it will be shown that a moving soliton can be captured by the defect. We will also study the stability of trapped solitons, concluding that they may be both unstable and stable, depending on the values of the system’s parameters. It is necessary to mention that the interaction between a BG soliton and local defect of the grating was considered in a very recent paper , which was chiefly focused on a detailed analysis of the case of finite-size defects. We, instead, consider point-like defects, which makes it possible to find an *exact* solution for trapped solitons, and an analytical potential of the interaction between a moving soliton and the defect. The study of the soliton-defect interaction reveals features that were not reported in Ref.  – most notably, reversal of the sign of the interaction from attraction to repulsion in some cases. Another interesting property which we present here is the fact that a stable stationary trapped soliton is, as a matter of fact, an attractor: up to the accuracy of numerical simulations, it is found to exists at a single value of its energy, but initial solitons from a broad area in the system’s parameter space relax to this stable state.

The paper is organized as follows: in section 2, the model is formulated, and analytical results are presented in section 3. Section 4 contains numerical results for the stability of trapped solitons, which show the existence of the attractor. In section 5, simulations of collisions of moving solitons with the defect are presented, the most remarkable feature being the above-mentioned reversal of the sign of the interaction; this, in particular, gives rise to two different regions in the parameter space, where the moving soliton is captured or repelled by the defect. The paper is concluded by section 6.

## 2 The model

We approximate a local disruption (or enhancement) of the resonant Bragg grating due to the presence of a defect by a $`\delta `$-function perturbation added to the standard nonlinear BG model for the right- and left-traveling waves $`u(x,t)`$ and $`v(x,t)`$. Taking into regard that the local perturbation may also affect the fiber’s refractive index, cf. Ref. , we arrive at the following model, written in the usual notation :

$`i{\displaystyle \frac{u}{t}}+i{\displaystyle \frac{u}{x}}+v+({\displaystyle \frac{1}{2}}|u|^2+|v|^2)u`$ $`=`$ $`\kappa \delta (x)v+\lambda \delta (x)u,`$ (1)
$`i{\displaystyle \frac{v}{t}}i{\displaystyle \frac{v}{x}}+u+({\displaystyle \frac{1}{2}}|v|^2+|u|^2)v`$ $`=`$ $`\kappa \delta (x)u+\lambda \delta (x)v,`$ (2)
where the ordinary ratio between the self-phase-modulation (SPM) and cross-phase-modulation (XPM) coefficients, $`1:2`$, is adopted, the local defect is placed at the point $`x=0`$, and the real constants $`\kappa >0`$ and $`\lambda `$ measure, respectively, the local suppression of the grating ($`\kappa <0`$ corresponds to a locally enhanced Bragg reflection) and local perturbation of the refractive index.

It should be mentioned that the region with a locally changed refractive index may also give rise to partial reflection, independently from that induced by BG. The corresponding perturbation is described by the cross-coupling linear terms proportional to the $`\delta `$-function, which are similar to those (proportional to $`\kappa `$) that are already present in Eqs. (1) and (2); however, the coefficients in front of these additional terms may be complex, due to the phase shift of the reflected wave. Thus, in the most general case, the terms $`\kappa \delta (x)v`$ and $`\kappa \delta (x)u`$ with real $`\kappa `$ in Eqs. (1) and (2) can be replaced, respectively, by $`\left(\kappa _1+i\kappa _2\right)\delta (x)v`$ and $`\left(\kappa _1i\kappa _2\right)\delta (x)u`$, with two real parameters, $`\kappa _1`$ and $`\kappa _2`$. However, in this work we confine the consideration to the case of real $`\kappa `$. As for the parameter $`\lambda `$, it must be strictly real in the conservative model (its imaginary part accounts for a local gain, if it is present in the model ).

If $`\kappa `$ is positive and small, the corresponding perturbation terms in Eqs. (1) and (2) imply that the Bragg reflection is disrupted in a region whose size is much smaller than a characteristic length necessary for the complete reflection, which is, typically, $`\mathrm{\hspace{0.17em}1}`$ mm , and which also determines a characteristic intrinsic spatial scale of the soliton. If, on the other hand, $`\kappa 1`$, this implies that the size of the disruption region is comparable to the reflection length. In the latter case, an exact form of the spatially distributed perturbation should be, strictly speaking, taken into regard (see Ref. ). Nevertheless, the use of the $`\delta `$-function approximation may be justified in this case too, as the fundamental features of the interaction of a soliton with the inhomogeneity may be weakly sensitive to the exact shape of the inhomogeneity. The latter property is well known to be true in many other applications of the soliton perturbation theory . In accordance with this, in the numerical part of the work we will assume that $`\kappa `$ takes values between $`0`$ and $`1`$. On the other hand, in the case $`\kappa <0`$, $`\left|\kappa \right|`$ may take any value, as in this case the $`\delta `$-function terms in Eqs. (1) and (2) have the literal meaning of a local perturbation in the form of the enhanced Bragg reflection.

It is relevant to mention that this difference in the interpretation of positive and negative values of $`\kappa `$ is very similar to the well-known difference between the positive and negative values of a perturbation parameter that determines the local suppression or enhancement (corresponding, respectively, to the so-called micro-resistor or micro-short) of the tunneling supercurrent in the standard model of a long Josephson junction, which is based on the sine-Gordon equation with a $`\delta `$-functional perturbation . Indeed, a micro-resistor may only make the existing critical density of the supercurrent smaller, while a micro-short may arbitrarily increase the local value of the critical density. As concerns the refractive-index perturbation parameter $`\lambda `$ in Eqs. (1) and (2), it may take arbitrary values, as it literally accounts for a local perturbation.

Notice that Eqs. (1) and (2) can be derived from the Hamiltonian,

$$H=_{\mathrm{}}^+\mathrm{}\left[\frac{1}{2}i\left(u^{}\frac{u}{x}+v^{}\frac{v}{x}+u\frac{u^{}}{x}v\frac{v^{}}{x}\right)\frac{1}{4}\left(\left|u\right|^4+\left|v\right|^4\right)\left|u\right|^2\left|v\right|^2\left(u^{}v+uv^{}\right)\right]+H_{\mathrm{int}}.$$

In this expression, the asterisk stands for the complex conjugation, the overall sign of $`H`$ is chosen so that the SPM and XPM terms correspond to self-focusing, and the perturbation part of the Hamiltonian, which accounts for the interaction of the soliton with the local defect, is

$$H_{\mathrm{int}}=\left[\kappa \left(u^{}v+uv^{}\right)+\lambda \left(\left|u\right|^2+\left|v\right|^2\right)\right]|_{x=0}.$$
(3)
## 3 Analytical results

We start by seeking for stationary solutions to Eqs. (1) and (2) in the form

$$u=U(x)\mathrm{exp}\left(i\omega t\right),v=V(x)\mathrm{exp}\left(i\omega t\right),$$
(4)
which leads to equations

$`\omega U+iU^{}+V+(\sigma |U|^2+|V|^2)U`$ $`=`$ $`\kappa \delta (x)V+\lambda \delta (x)U,`$ (5)
$`\omega ViV^{}+U+(\sigma |V|^2+|U|^2)V`$ $`=`$ $`\kappa \delta (x)U+\lambda \delta (x)V,`$ (6)
the prime standing for $`d/dx`$. A family of exact solutions to the unperturbed equations (5) and (6), with $`\kappa =\lambda =0`$, is well known  (these solutions are frequently called gap solitons):

$$U(x)=\sqrt{\frac{2}{3}}(\mathrm{sin}\theta )\mathrm{sech}\left(\left(x\xi \right)\mathrm{sin}\theta \frac{i}{2}\theta \right),V=U^{},$$
(7)
where $`\xi `$ is the coordinate of the soliton’s center, and $`\theta `$ is a real parameter taking values

$$0<\theta <\pi ,$$
(8)
which also determines the frequency in Eqs. (4),

$$\omega =\mathrm{cos}\theta .$$
(9)
### 3.1 An exact solution for the pinned soliton

An *exact* solution for a stationary soliton trapped by the local defect can be found in the form \[cf. Eq. (7)\]

$$U=\sqrt{\frac{2}{3}}(\mathrm{sin}\theta )\mathrm{sech}\left[\left(x+a\mathrm{sgn}x\right)\mathrm{sin}\theta \frac{i}{2}\theta \right],V=U^{},$$
(10)
where an expression for the parameter $`a`$ is obtained by integrating Eqs. (5) and (6) in an infinitesimal vicinity of the point $`x=0`$:

$$\mathrm{tanh}\left(a\mathrm{sin}\theta \right)=\sqrt{\frac{\kappa \lambda }{\kappa +\lambda }}\frac{\mathrm{tanh}\left(\sqrt{\kappa ^2\lambda ^2}/2\right)}{\mathrm{tan}\left(\theta /2\right)}.$$
(11)
An example of this exact solution, with $`\kappa =0.4`$ and $`\theta =0.5\pi `$, is shown in Fig. 1.

In Eq. (11), we assume that $`\kappa ^2>\lambda ^2`$ (the opposite case will be considered separately below), and $`\sqrt{\kappa ^2\lambda ^2}`$ is realized in such a way that it has the same sign as $`\kappa `$ itself. In particular, in the case $`\lambda =0`$, Eq. (11) takes the form

$$\mathrm{tanh}\left(a\mathrm{sin}\theta \right)=\frac{\mathrm{tanh}\left(\kappa /2\right)}{\mathrm{tan}\left(\theta /2\right)}.$$
(12)
The intensity of the field at the central point of the trapped soliton is

$$\left|U\left(x=0\right)\right|^2=\frac{2}{3}\frac{\mathrm{sin}^2\theta }{\mathrm{sinh}^2\left(a\mathrm{sin}^2\theta \right)+\mathrm{cos}^2\theta }.$$
(13)
Note that $`\left|U(x)\right|`$ is continuous in the solution (10) across the point $`x=0`$, while the phase of $`U(x)`$ has a discontinuity at this point.

If $`\kappa <0`$, Eqs. (11) and (12) yield $`a<0`$, hence, as it follows from Eq. (10), in this case $`\left|U(x)\right|`$ has two symmetric local maxima at $`x=a`$, and a local minimum at the central point $`x=0`$. On the contrary to this, if $`\kappa >0`$, we have $`a>0`$, which implies that the solution has a single maximum at $`x=0`$. These features can be readily understood by dint of the interaction Hamiltonian (3). Indeed, together with the relation $`V=U^{}`$, see Eq. (7), the expression (3) yields the following value of the interaction Hamiltonian for the trapped soliton:

$$\left(H_{\mathrm{int}}\right)_{\mathrm{trapped}}=2\left(\kappa \lambda \right)\left|U\left(x=0\right)\right|^2.$$
(14)
With regard to the constraint $`\kappa ^2>\lambda ^2`$ adopted above, we conclude that this contribution increases the net Hamiltonian if $`\kappa `$ is negative, and decreases it in the opposite case. Therefore, following the general principle stating that the ground state of a system corresponds to the minimum of $`H`$, the system tries to make $`\left|U\left(x=0\right)\right|^2`$ smaller if $`\kappa <0`$, and larger if $`\kappa >0`$.

The result for the case $`\kappa ^2<\lambda ^2`$ can be obtained from Eq. (11) by analytical continuation \[or directly from Eqs. (5) and (6)\]:

$$\mathrm{tanh}\left(a\mathrm{sin}\theta \right)=\sqrt{\frac{\lambda \kappa }{\lambda +\kappa }}\frac{\mathrm{tan}\left(\sqrt{\lambda ^2\kappa ^2}/2\right)}{\mathrm{tan}\left(\theta /2\right)},$$
(15)
where $`\mathrm{sgn}\left(\sqrt{\lambda ^2\kappa ^2}\right)\mathrm{sgn}\lambda `$. In the particular case $`\kappa =0`$, Eq. (15) simplifies to

$$\mathrm{tanh}\left(a\mathrm{sin}\theta \right)=\frac{\mathrm{tan}\left(\lambda /2\right)}{\mathrm{tan}\left(\theta /2\right)},$$
(16)
cf. Eq. (12). Equations (15) and (16) yield $`a>0`$, i.e., the single maximum of $`\left|U(x)\right|`$ at $`x=0`$, if $`\lambda `$ is negative, and $`a<0`$, i.e., a local minimum at $`x=0`$ and two maxima at $`x=a`$, if $`\lambda `$ is positive. Taking into regard the present constraint $`\kappa ^2<\lambda ^2`$, it is easy to check that these conclusions are again in accordance with the principle of minimization of the interaction Hamiltonian (14).

In the special case $`\kappa =\lambda `$, both Eqs. (11) and (15) amount to

$$\mathrm{tanh}\left(a\mathrm{sin}\theta \right)=\frac{\kappa }{\mathrm{tan}\left(\theta /2\right)}\frac{\lambda }{\mathrm{tan}\left(\theta /2\right)}.$$
(17)
Lastly, in the case $`\kappa =\lambda `$, both expressions (11) and (15) yield $`a=0`$, which is easy to understand: due to the relation $`V=U^{}`$, the unperturbed soliton (7) remains an exact solution to Eqs. (5) and (6) if $`\kappa =\lambda `$.

Obviously, the expressions (11), (12), and (15) through (17) make sense if they yield $`\left|\mathrm{tanh}\left(a\mathrm{sin}\theta \right)\right|<1`$, otherwise the solution for the trapped soliton does not exist. The latter restriction means that the exact solution exists not in the whole interval (8) in which the unperturbed solution (7) was found, but rather in a narrower band,

$$\theta _{\mathrm{min}}<\theta <\pi ,$$
(18)
where, in both cases $`\kappa ^2{}_{>}{}^{<}\lambda _{}^{2}`$,

$$\theta _{\mathrm{min}}=2\mathrm{tan}^1\left(\left|\sqrt{\frac{\kappa \lambda }{\kappa +\lambda }}\mathrm{tanh}\left(\sqrt{\kappa ^2\lambda ^2}/2\right)\right|\right).$$
(19)
This expression takes the form $`\theta _{\mathrm{min}}=\lambda `$ in the case $`\kappa =0`$, and in the case $`\lambda =\kappa `$ it is $`\theta _{\mathrm{min}}=2\mathrm{tan}^1\left(\left|\kappa \right|\right)`$.

The existence limit for the trapped solitons, as given by Eqs. (18) and (19), was also verified by numerical methods, solving Eqs. (5) and (6) with some approximations for $`\delta (x)`$, see details below. The numerical solutions (not shown here, as they are not especially interesting by themselves) demonstrate good agreement with the analytically predicted existence limit.

### 3.2 A remark about the soliton stability

Some (incomplete) conclusions about the stability of the trapped solitons can be made on the basis of the known Vakhitov-Kolokolov (VK) criterion , which states that a necessary condition of the stability of the soliton is $`dE/d\omega <0`$, where $`\omega `$ is the soliton’s frequency, see Eq. (9), and $`E_{\mathrm{}}^+\mathrm{}\left[\left|u(x)\right|^2+\left|v(x)\right|^2\right]𝑑x`$ is the energy (norm) of the solution. In the general case, the expression for the energy of the exact solution (10) is cumbersome, but it takes a simple form in the case $`\lambda =0`$, on which the numerical simulations will be focused below:

$$E=\frac{8}{3}\left(\theta \left[\frac{\pi }{2}\mathrm{sin}^1\left(\mathrm{sech}\kappa \right)\right]\mathrm{sgn}\kappa \right).$$
(20)
The calculation of the derivative $`dE/d\omega `$ by means of Eqs. (20) and (9) shows that all the trapped-soliton solutions satisfy the VK criterion. However, this result is not sufficient for claiming stability of the solutions. Indeed, in the unperturbed case ($`\kappa =0`$) the VK criterion per se predicts that all the solitons (7) are stable. However, studies of the stability, performed originally within the framework of the variational approximation , and then by means of accurate numerical methods , had demonstrated that, in fact, the interval (8) in which the solitons exist is divided into two parts: they are stable if

$$0<\theta <\theta _{\mathrm{cr}},\theta _{\mathrm{cr}}1.011\left(\pi /2\right),$$
(21)
and unstable if $`\theta _{\mathrm{cr}}<\theta <\pi `$. Notice that $`\theta _{\mathrm{cr}}`$ is just slightly larger than $`\pi /2`$, which is the value at which the gap soliton (7) has the largest amplitude and smallest width. This example clearly suggests that direct investigation of the solitons’ stability is necessary in the present model too, which will be done by means of numerical methods in section 4.

### 3.3 The potential of the soliton-defect interaction

Besides the exact findings presented above, more general analytical results can be obtained in an approximate form, assuming that $`\kappa `$ and $`\lambda `$ are small parameters. In the lowest (adiabatic) approximation of the perturbation theory , one can find an effective potential of the interaction between the soliton, treated as a quasi-particle, and the local defect. To this end, the unperturbed wave form, based on Eqs. (4) and (7), is substituted into the expression (14) for the interaction Hamiltonian, which yields

$$U_{\mathrm{int}}(\xi )=\frac{8}{3}\frac{\left(\lambda \kappa \right)\mathrm{sin}^2\theta }{\mathrm{cosh}\left(2\xi \mathrm{sin}\theta \right)+\mathrm{cos}\theta }$$
(22)
(recall $`\xi `$ is the coordinate of the soliton’s center). In particular, in the case $`\lambda =0`$ that will be considered below, the potential (22) clearly implies attraction in the case $`\kappa >0`$, and repulsion in the opposite case. The latter circumstance strongly suggests that the trapped soliton cannot be stable in the case $`\kappa <0`$, $`\kappa ^2>\lambda ^2`$ (nor in the case $`\lambda >0`$, $`\lambda ^2>\kappa ^2`$).

## 4 Numerical simulations of the stability of the pinned soliton

### 4.1 The approximation for the delta-function

To simulate the stability of the trapped solitons and the interaction of a free soliton with the local defect, we have to adopt an approximation to represent the $`\delta `$-functions in Eqs. (1),(2), (5) and (6). In this work, we use a numerical scheme in which the coordinate $`x`$ is represented by $`501`$ grid points $`x_j`$, $`j=250,\mathrm{},1,0,+1,\mathrm{}+250`$. As an approximation to the $`\delta `$-function, we have defined the following function on a set of $`2N+1`$ grid points in the central part of the integration domain, located symmetrically around zero:

$$\stackrel{~}{\delta }\left(x_{n\left(N+1\right)}\right)\begin{array}{cc}A\mathrm{cos}\left(\frac{n\left(N+1\right)}{2N+1}\pi \right)\hfill & \mathrm{for}n=1,\mathrm{},2N+1,\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}$$
(23)
The normalization factor is defined so as to maintain the canonical normalization of the $`\delta `$-function, $`_{\mathrm{}}^+\mathrm{}\stackrel{~}{\delta }(x)𝑑x_j\stackrel{~}{\delta }\left(x_j\right)\mathrm{\Delta }x=1`$, which yields

$$A=\left[\mathrm{\Delta }x\underset{n=1}{\overset{2N+1}{}}\mathrm{cos}\left(\frac{n\left(N+1\right)}{2N+1}\pi \right)\right]^1,$$
(24)
$`\mathrm{\Delta }x`$ being the spacing of the grid (in fact, it was $`0.04)`$.

In most cases presented below, we use $`N=2`$ \[then Eq. (24) with $`\mathrm{\Delta }x=0.04`$ yields $`A=\left(\left(1+\sqrt{5}\right)\mathrm{\Delta }x\right)^1\mathrm{7.\hspace{0.17em}726}`$\], which makes the $`\delta `$-function quite narrow indeed. Besides that, we also used another approximation for the $`\delta `$-function, in which a simple rectangular regularization was used instead of (23). In some cases, sharp edges of the rectangular approximation produced problems with the convergence, but in those cases when results converged, they were virtually the same as generated by the approximation (23).

It should be noted that, in the case $`\kappa A>1`$ \[recall $`\kappa `$ is the perturbation parameter in Eqs. (1) and (2)\], the local value of the coefficient $`\kappa \stackrel{~}{\delta }(x)`$ may exceed $`1`$, which is, strictly speaking, an unphysical situation, as it was discussed above. However, we consider such a situation as modeling a physical one with a broader defect and smaller amplitude, when the integral strength of the perturbation, $`\kappa _{\mathrm{}}^+\mathrm{}\stackrel{~}{\delta }(x)𝑑x`$, takes the same value, while the coefficient $`\kappa \stackrel{~}{\delta }(x)`$ does not exceed $`1`$. Comparing results produced in simulations by the formally unphysical perturbation, and by its physical counterpart, we have verified that this assumption is, to a large extent, correct, with some minor differences observed only when the defect strength $`\kappa `$ is small. This issue will be further discussed below.

### 4.2 The stationary soliton

In simulations of the stability of the trapped solitons and the interaction of a free soliton with the local defect (the interaction is the subject of the next section), we concentrated on the case $`\lambda =0`$, which is the most interesting one, as it focuses on the effects induced by the defect in the BG proper. We will mainly assume $`\kappa >0`$ (although the opposite case will also be briefly considered), which corresponds to the most natural situation with the local perturbation suppressing the local Bragg reflectivity (rather than enhancing it). This is also justified because, as conjectured above and shown below, all the trapped solitons are unstable in the case $`\kappa <0`$.

Our first objective is to analyze the stability of a stationary soliton trapped by the local defect. To this end, instead of using the exact solution given by Eqs. (10) and (11) through (17), we directly constructed trapped-soliton solutions to Eqs. (5) and (6), with the regularized $`\delta `$-function taken as per Eq. (23), by means of the Newton-Raphson method. An example of the thus found stationary soliton is shown in Fig. 1, along with the exact solution (10), for $`\kappa =0.4`$ and $`\theta =\pi /2`$. After obtaining the stationary solutions in the numerical form, it was found that they are always close to the analytical solution for the ideal $`\delta `$-function given by Eqs. (10) and (11), so that the value of the parameter $`\theta `$ corresponding to each numerically found solution could be easily identified.

### 4.3 Stability simulations

The stability of the stationary solutions was directly simulated by means of the split-step method applied to Eqs. (1) and (2), employing the fast Fourier transform. First of all, in the case $`\kappa <0`$ the numerical solution confirms that, as it follows from Eq. (10), $`\left|U(x)\right|`$ has two symmetric local maxima at $`x=a`$, and a local minimum at the central point $`x=0`$. Dynamical simulations confirm that *all* the trapped solitons with $`\kappa <0`$ are unstable, as it was assumed above on the basis of the perturbation theory. The simulations further show that, as a result of the instability, two local maxima at $`x=a`$ start to separate and eventually move away from the defect, so that the trapped soliton splits into two free separating ones. Since all the trapped solitons are unstable for $`\kappa <0`$, the investigation of this case was not carried out further.

Typical results illustrating the stability of the trapped solitons in the case $`\kappa >0`$ are displayed in Fig. 2, which shows the evolution of the initial pulse in the form of the trapped soliton (10), with its center placed at $`x=0`$, the value of the perturbation parameter being $`\kappa =0.08`$ (note that this case is a “strictly physical” one, as it corresponds to $`\kappa A<1`$, see above). In the series of panels displayed in Fig. 2, $`\theta `$ is given values $`\theta _{\mathrm{in}}=0.4\pi ,\mathrm{\hspace{0.17em}0.5}\pi ,\mathrm{\hspace{0.17em}0.7}\pi ,0.9\pi `$. The left part of each panel shows a “side view” of the evolution of the soliton (in terms of values of $`\left|U\right|`$), the horizontal axis being time. The right parts of the panels (sometimes, these are insets) show a “top view” of the evolution, in terms of contour plots representing the solution. The evolution of solitons in the whole region $`\kappa _{}^<\mathrm{\hspace{0.17em}\hspace{0.17em}0.2}`$ (note that it comprises both the “strictly physical” and formally unphysical situations) is found to be essentially the same as shown in Fig. 2 for $`\kappa =0.08`$.

As it is seen from the plot in Fig. 2(b), and from the results for other values of $`\kappa `$, only a single value of $`\theta _{\mathrm{in}}`$ directly gives rise to a stable soliton (it is single up to the resolution provided by the accumulated data, i.e., we cannot rule out that it may be a narrow interval, rather than the single value). It is also found, more importantly, that in *all* the cases when the initial pulse relaxes into a stable soliton, the corresponding final value $`\theta _{\mathrm{stab}}`$ of $`\theta `$ is also a single one, being equal or very close to $`\pi /2`$, up to the accuracy of the numerical data. It is noted that $`\theta =\pi /2`$ yields the maximum amplitude and minimum width of the soliton as per Eq. (10), and it belongs to the stability region (21) of the unperturbed solitons, although being very close to the instability border. We also stress that the value $`\theta =\pi /2`$ always belongs, in the case $`\lambda =0`$, to the existence band (18) of the pinned-soliton solutions.

The exact solution (10) for the trapped soliton with $`\theta =\pi /2`$ takes (provided that $`\lambda =0`$) a very simple form: as it follows from Eq. (12), in this case $`a=\kappa /2`$. Note also that, if $`\kappa `$ is large, the latter result implies that the maximum value of the soliton’s field is very small: according to Eq. (13), $`|U(x=0|^2=(4/3)(\mathrm{cosh}\kappa 1)^1`$. In fact, this means that, in the limit of large $`\kappa `$, the trapped soliton becomes a quasi-linear defect-supported mode, cf. Ref. .

The fact that there is a (presumably) single stable trapped soliton in the present model implies that the model gives rise to an effective  attractor. Although attractors are usually considered as objects specific to dissipative systems, they may occur in conservative nonlinear-wave models due to the possibility of radiation losses.

If $`\theta _{\mathrm{in}}`$ is essentially smaller than $`\theta _{\mathrm{stab}}`$, the soliton cannot relax to the stable state and may only decay into radiation, see Fig. 2(a). In fact, the expression (20) for the energy of the trapped soliton imposes an absolute limit on the values of $`\theta _{\mathrm{in}}^{(0)}`$ of the initial unperturbed soliton (the one corresponding to $`\kappa =0`$) which can self-trap into the stable state. Indeed, the energy $`E_{\mathrm{in}}^{(0)}`$ of the initial unperturbed soliton is given by Eq. (20) with $`\kappa =0`$, i.e., $`E_{\mathrm{in}}^{(0)}=(8/3)\theta _{\mathrm{in}}^{(0)}`$, and, obviously, the energy $`E_{\mathrm{fin}}`$ of the established state can only be smaller than $`E_{\mathrm{in}}`$. Taking for $`E_{\mathrm{fin}}`$ the expression (20) with $`\theta =\theta _{\mathrm{stab}}\pi /2`$, we find $`E_{\mathrm{fin}}=(8/3)\mathrm{sin}^1\left(\mathrm{sech}\kappa \right)`$. Thus, the limitation $`E_{\mathrm{fin}}<E_{\mathrm{in}}^{(0)}`$ takes a final form

$$\theta _{\mathrm{in}}^{(0)}>\left(\theta _{\mathrm{in}}^{(0)}\right)_{\mathrm{min}}\mathrm{sin}^1\left(\mathrm{sech}\kappa \right).$$
(25)
Obviously, $`\left(\theta _{\mathrm{in}}^{(0)}\right)_{\mathrm{min}}`$ monotonically decreases with $`\kappa `$, from $`\pi /2`$ at $`\kappa =0`$ to $`0.2244\pi `$ when $`\kappa =1`$ and further to zero as $`\kappa \mathrm{}`$.

In order to investigate a possible dependence of $`\theta _{\mathrm{stab}}`$ on the defect’s strength $`\kappa `$, simulations were run for a range of values of $`\theta `$ slightly smaller than $`\pi /2`$, with $`\kappa `$ varying from $`0.1\pi `$ to $`0.9\pi `$. Similar to what is depicted in Fig. 2(a) and described above, solitons with $`\theta _{\mathrm{in}}<\pi /2`$ decay, but at a slower and slower rate, as $`\theta _{\mathrm{in}}`$ is getting closer to $`\pi /2`$. For the simulation time $`t=300\pi `$ and $`\kappa =0.1`$, the decrease of the amplitude of $`|u|`$ and $`|v|`$ is larger than $`1\%`$ if $`\theta _{\mathrm{in}}<0.49\pi `$. Of course, the decay stops as $`\kappa 0`$. However, the decay also becomes slower for larger values of $`\kappa `$ (which can be understood too: in the limit of $`\kappa +\mathrm{}`$, the trapped soliton goes over into the above-mentioned linear defect mode, which is stable). For instance, at $`\kappa =0.9`$, the loss of the amplitude observed at $`t=300\pi `$ is less that $`1\%`$ if $`0.47\pi <\theta _{\mathrm{in}}<\pi /2`$. So, there may be two possibilities: either $`\theta _{\mathrm{stab}}`$ is getting slightly smaller than $`\pi /2`$ when $`\kappa `$ is getting larger, or $`\theta _{\mathrm{stab}}`$ remains equal to $`\pi /2`$, but the solitons are less unstable, so the decay rate for $`\theta _{\mathrm{in}}<\theta _{\mathrm{stab}}`$ is smaller for larger values of $`\kappa `$. Because of the very slow decay or settling rate of $`|u|`$ and $`|v|`$ for $`\theta _{\mathrm{in}}`$ close to $`\pi /2`$, we cannot pinpoint a single value of $`\theta _{\mathrm{stab}}`$ value; however, a very conservative estimate is that, at all the values of $`\kappa `$ considered, $`\theta _{\mathrm{stab}}`$ is confined to the interval $`0.47\pi <\theta _{\mathrm{stab}}<0.53\pi `$.

Figure 2(c) shows that if the initial value $`\theta _{\mathrm{in}}`$ is somewhat larger than $`\theta _{\mathrm{stab}}`$, the initial soliton evolves so that $`\theta `$ relaxes to $`\theta _{\mathrm{stab}}`$, shedding off excess energy in the form of radiation waves. If $`\theta _{\mathrm{in}}`$ is still larger for small values of $`\kappa `$ (actually, if $`\theta _{\mathrm{in}}{}_{}{}^{>}\mathrm{\hspace{0.17em}\hspace{0.17em}0.8}\pi `$ at $`\kappa =0.08`$, see Fig. 3), the soliton does not relax into the stable state, but rather decays into radiation, sometimes generating a small-amplitude residual free soliton which is flung away from the point where the defect is located \[see Fig. 2(d)\]; spontaneous symmetry breaking, apparent in the latter case, is quite possible due to the instability. However, for $`\kappa _{}^>\mathrm{\hspace{0.17em}\hspace{0.17em}0.2}`$, the decay of the soliton at large values of $`\theta _{\mathrm{in}}`$ is no longer observed.

Figure 3 summarizes the results in the form of a stability region in the plane ($`\kappa ,\theta _{\mathrm{in}}`$), in which the initial soliton evolves into the stable one. It is quite natural that the stability region expands with the increase of $`\kappa `$, as it is easier for a stronger defect to pull the field and thus capture the soliton. Besides that, Eq. (20) shows that the energy of the trapped soliton decreases with $`\kappa `$, which also makes it easier to create a soliton from a given initial pulse.

Lastly, one may notice that, if the value $`\theta _{\mathrm{min}}`$ limiting the existence region of the trapped solitons, see Eqs. (18) and (19), exceeds $`\theta _{\mathrm{stab}}\pi /2`$, the stable trapped soliton cannot exist. However, in the case under consideration, $`\lambda =0`$ and $`\kappa >0`$, Eq. (19) yields $`\theta _{\mathrm{min}}=2\mathrm{tan}^1\left(\mathrm{tanh}(\kappa /2)\right)`$. Since $`\mathrm{tanh}(\kappa /2)<1`$ for all value of $`\kappa `$, it is obvious that $`\theta _{\mathrm{min}}`$ is always smaller than $`\pi /2`$ (for example, $`\theta _{\mathrm{min}}=0.276\pi `$ when $`\kappa =1`$), hence the condition $`\theta _{\mathrm{stab}}>\theta _{\mathrm{min}}`$ never excludes $`\theta _{\mathrm{stab}}`$ from the existence region.

## 5 Trapping and reflection of a moving soliton by the localized defect

The interaction of a moving soliton with the defect is a problem of obvious interest, both in its own right, and for such applications mentioned above as sensing by means of solitons and the soliton-based dynamical optical memory. To simulate the interactions, solutions for free moving solitons were first generated far from the localized defect in a numerical form by means of the Newton-Raphson method, as stationary solutions in the moving reference frame (these solutions are also available in the analytical form , which shows that the soliton’s velocity $`c`$ may take any value from $`1<c<+1`$). Then, full simulations of Eqs. (1) and (2) were run by means of the split-step method combined with the fast Fourier transform, launching the moving soliton to collide it with the localized defect (it was again set $`\lambda =0`$ and $`\kappa >0`$).

First, the parameter plane ($`c,\kappa `$) was explored, with $`c`$ taking values $`0.05,\mathrm{\hspace{0.17em}0.075},\mathrm{\hspace{0.17em}0.1},\mathrm{\hspace{0.17em}0.2},\mathrm{},0.7`$ and $`\kappa `$ ranging from $`0.07`$ to $`0.9`$. In all the simulations, the soliton was launched with $`\theta =\pi /2`$, which is equal (or very close) to the single value $`\theta _{\mathrm{stab}}`$ which, as it was demonstrated above, is selected by the localized defect as the single value at which the trapped soliton is stable. Figure 4 displays several different generic examples of the interaction of the moving soliton with the defect. In each panel, the lower part shows the evolution of the field $`\left|u(x,t)\right|`$ by means of contour plots, and the upper part shows the waveforms $`\left|u(x)\right|`$ and $`\left|v(x)\right|`$ at the end of the simulation, the location of the defect being marked by a short vertical line.

If the defect is weak, the soliton passes through it, decreasing its velocity due to radiation loss. However, a stronger defect captures the soliton into an oscillatory state, see Fig. 4(a). The possibility of the capture of the soliton (through the radiation loss) is in agreement with the predicted attractive character of the interaction between the soliton and the defect, see Eq. (22).

Simulations demonstrate that, quite naturally, the minimum value of the defect’s strength $`\kappa `$ which is necessary for the capture of a given free soliton increases with its velocity. Figure 5 shows a border between the trapping (B) and transmission (A) regions in the plane ($`c,\kappa `$), for fixed $`\theta =\pi /2`$.

However, the soliton-defect interaction is more complicated. When the defect strength $`\kappa `$ is large, and when the soliton is energetic enough \[either $`c`$ or $`\theta `$ is large – region (C) in Fig. 5\], its energy splits into three parts as a result of the collision: one part is trapped by the defect, and two others are scattered away, in both directions. The energy which is not trapped (either bouncing back or the passing through the defect) self-traps into a secondary soliton. Three examples of this mode of the interaction are shown in Figs. 4(c,d,e), for fixed values $`c=0.5`$ and for $`\theta =0.5\pi `$, and increasing values of the defect’s strength, $`\kappa =0.3,\mathrm{\hspace{0.17em}0.6},`$and $`0.8`$. It can be seen that when $`\kappa `$ becomes larger, shares of the energy which are trapped and reflected back become larger. In regions A and B in Fig. 5, some energy is also scattered in the forward or backward directions; however, in these regions scattered energy is small, and it completely disperses into radiation, without giving rise to any secondary soliton.

Another characteristic set of numerical data can be displayed for a fixed small value of the soliton’s velocity, $`c=0.075`$, while the parameter $`\theta `$ of the moving soliton and the defect’s strength $`\kappa `$ are varied. In this representation, the simulations demonstrate that the defect behaves in a qualitatively different way for small and large values of $`\theta `$. Namely, for $`\theta _{}^<\mathrm{\hspace{0.17em}\hspace{0.17em}0.6}\pi `$, the soliton interacts with the defect attractively, in accordance with the prediction of the perturbation theory, while for $`\theta _{}^>\mathrm{\hspace{0.17em}\hspace{0.17em}0.6}\pi `$ the interaction becomes *repulsive*. The change of the sign of the interaction between the soliton and defect was inferred from the observation of outcomes of the collision.

In the case $`\theta _{}^<\mathrm{\hspace{0.17em}\hspace{0.17em}0.6}\pi `$, the soliton gets captured by the defect, which clearly testifies that the interaction is attractive. The capture is accompanied by emission of a conspicuous jet of radiation. Passage of the soliton through the attractive defect is also possible (this happens, for instance, if the initial velocity of the soliton increases).

In a case when the defect strength is $`\kappa =0.5`$, and the soliton is heavier, having $`\theta =0.7\pi `$, the soliton bounces back from the defect, which is a characteristic feature of repulsion. The reversal of the interaction from the attraction to repulsion, which is observed for the solitons with $`\theta _{}^>\mathrm{\hspace{0.17em}\hspace{0.17em}0.6}\pi `$, is beyond the applicability limits of the simple perturbation theory that gave rise to the effective interaction potential (22). However, it agrees with the findings of the numerical stability analysis reported in the previous section, where it has been shown that solitons with large initial values of $`\theta _{\mathrm{in}}`$ are unstable if trapped by the defect, which may be realized as a result of repulsion.

Figure 7 summarizes the results obtained for the interaction of the moving soliton with the local defect in the present case (fixed velocity). The border between attraction (A) and repulsion (B) regions is shown. To the left of the regions (A) and (B), there is a thin strip (C) where the defect is weak, $`\kappa _{}^<\mathrm{\hspace{0.17em}0.1}`$, and the soliton is heavy enough, $`\theta _{}^>\mathrm{\hspace{0.17em}0.45}`$. In this case, the soliton is neither captured nor reflected by the defect, but rather passes through it with significant radiation losses, see Fig. 6. By extrapolating the border between the regions (A) and (B), it can be inferred that the region (C) is made up of two sub-regions: an upper one, where the defect is effectively repulsive, and a lower one, where the defect is attractive. However, in either case the “heavy” soliton passes through the weak defect.

Lastly, in the case when the defect is very strong, and the soliton is very heavy, the interaction may combine attractive and repulsive features. In this case, the energy of the incoming soliton is split into three parts: one is trapped at the defect, while the rest of the energy is scattered in both the forward and backward directions. The energies scattered away are large enough for self-trapping into secondary solitons. This happens in the region (D) of Fig. 7, and is similar to what happens in the region (C) of Fig. 5, where the defect was very strong, and the soliton was very fast. As the latter outcome of the collision is very similar to what was shown in Figs. 4(c,d,e), separate figures for this case are not included.

As well as in Fig. 5, the region (D) in Fig. 7 is distinguished as a separate one to emphasize that *significant* amounts of energy are scattered away from the defect. In the other regions, some energy is scattered too, but its amount is much smaller.

Concerning the results presented in this section, it is relevant to mention that moving unperturbed solitons also have their intrinsic stability limit, which, however, very weakly depends on the velocity : as $`c`$ increases from $`0`$ to $`1`$, the critical value $`\theta _{\mathrm{cr}}`$, above which the solitons are unstable, increases from $`1.011\left(\pi /2\right)`$ \[see Eq. (21)\] to $`1.017\left(\pi /2\right)`$. This implies that, strictly speaking, the incident soliton is unstable in all the cases shown in Fig. 6; however, its instability is very weak, therefore it does not manifest itself in any way before the collision is completed (for the same reason, formally unstable solitons may seem fairly stable in an experiment).

Finally, we return to the discussion of the validity of Eqs. (23) and (24) as an approximation for the defect when $`\kappa A>1`$, which corresponds to the formally unphysical situation. We stated above that such a situation is an appropriate model for a physical one, with a broader defect and smaller amplitude $`A`$, so that the product $`\kappa \stackrel{~}{\delta }(x)`$ does not exceed $`1`$, provided that both the “physical” and “unphysical” perturbations are subject to the identical normalization, $`_{\mathrm{}}^+\mathrm{}\stackrel{~}{\delta }(x)𝑑x=1`$.

As a possibility to verify this statement explicitly, we took the approximation (23) with $`N`$ $`=19`$ \[in the latter case, the normalization condition (24) yields $`A=1`$\], and performed simulations of the collision of the soliton with this broad defect in the whole $`(\kappa ,\theta )`$ parameter plane, keeping the soliton’s velocity constant, $`c=0.075`$. The results were then compared to those displayed in Figure 7, which had been obtained for the narrow defect (recall $`N=2`$ in the case of Fig. 7). It was found that, to a large extent, the two sets of the results are similar indeed. The agreement is very good for larger values of $`\kappa `$. Some disagreement occurs in the region (C), where $`\kappa `$ is small. In Fig. 7, the vertical part of the border of the region (C) is located at $`\kappa 0.1`$, with its tip at $`\theta 0.45`$. For the broad defect with $`N=19`$, it was found that this border moves out to $`\kappa 0.25`$, and the tip of the region moves to $`\theta 0.55`$, so that the region (C) is wider but shorter. This means that the conclusion whether the solitons can pass through the defect depends not solely on its integral strength $`\kappa _{\mathrm{}}^+\mathrm{}\stackrel{~}{\delta }(x)𝑑x`$, but also on its amplitude $`\kappa A`$. However, since this affects a small part of the $`(\kappa ,\theta )`$ parameter plane, the above conjecture concerning the actual equivalence of the narrow (formally unphysical) and broad approximations for the $`\delta `$-function is quite acceptable.

## 6 Conclusion

In this work, we have studied in detail the dynamics of the interaction of gap solitons with a localized defect, which is realized as a local suppression of the Bragg grating that supports the gap solitons. A family of exact solutions for trapped solitons in the model with the $`\delta `$ \- functional defect, and an approximate expression for the effective potential of the soliton-defect interaction, were found in the analytical form. Direct simulations have demonstrated that this conservative model gives rise to an effective attractor: up to the accuracy provided by numerical data, the trapped soliton is stable at a single value of its intrinsic parameter $`\theta `$ – in fact, at $`\theta =\pi /2`$, that corresponds to the soliton with the largest amplitude and smallest width. Depending on the strength $`\kappa `$ of the defect, initial solitons from a broad domain in the parametric plane ($`\kappa ,\theta `$) relax into the stable state; outside the attraction domain, they decay into radiation (sometimes, generating a residual free soliton).

Simulations of the interaction between moving solitons and the defect produce a number of different outcomes. If the soliton’s velocity and $`\theta `$ are not too large, the soliton is captured by the defect, in accordance with the fact the effective potential of the interaction between them is attractive. When the soliton is “heavy”, with $`\theta _{}^>\mathrm{\hspace{0.17em}0.6}`$, the character of the interaction reverses from attraction to repulsion. In this case, the soliton bounces back. However, if the defect is weak, regardless of whether the interaction is attractive or repulsive, heavy solitons can pass through the defect. If it is strong, and the soliton is heavy or fast enough, the collision may split the soliton in three parts. Some energy is trapped by the defect to form a pinned soliton, while large amounts of energy scattered away may generate secondary solitons (transmitted and reflected ones) of significant amplitude.

These results suggest various applications. Besides the optical-memory design and sensing (defect detection), which were mentioned in the Introduction, the delay of a soliton passing through the repulsive defect can be used to construct an optical delay line. Besides that, an optical resonator may be implemented by means of a soliton trapped in a cavity formed by two localized repulsive defects, with a locally enhanced Bragg reflection. In such setups, the local defects can be switched on and off in a controllable way by means of electrostriction.

## Acknowledgement

B.A. Malomed appreciates hospitality of the Department of Electronic Engineering at the City University of Hong Kong.

## Figure Captions

Fig. 1. The real (even) and imaginary (odd) parts of the field $`U(x)`$ for the trapped soliton, with $`\kappa =0.4`$ and $`\theta =0.5\pi `$. The dashed curves depict the exact analytical solution (10) obtained for the ideal $`\delta `$-function. The solid curves show the numerical solution, with the $`\delta `$ \- function approximated as per Eq. (23) with $`N=2`$. The two pairs of curves are nearly indiscernable. It was also checked that using an approximation closer to the ideal $`\delta `$ \- function brings the numerically found waveform still closer to the exact solution.

Fig. 2. Typical results illustrating the stability and instability of the trapped solitons in the case $`\lambda =0`$ and $`\kappa =0.08`$. The left part of each panel shows a “side view” of the evolution of $`|U|`$, starting with the initial stationary trapped soliton configuration. Right parts of the panels (or insets) show the evolution in terms of contour plots. $`T`$ is the total simulation time. (a) $`\theta _{\mathrm{in}}=0.4\pi `$ is smaller than $`\theta _{\mathrm{stab}}`$. The soliton decays into radiation. (b) $`\theta _{\mathrm{in}}=0.5\pi \theta _{\mathrm{stab}}`$. The field directly self-traps into a stable soliton. (c) $`\theta =0.7\pi >\theta _{\mathrm{stab}}`$. The field evolves into a stable soliton, shedding off excess energy in the form of radiation. (d) $`\theta =0.9\pi `$. This value is much larger than $`\theta _{\mathrm{stab}}`$, and the pulse decays into radiation, generating a small residual soliton thrown away from the defect.

Fig. 3. The region in the plane ($`\kappa ,\theta _{\mathrm{in}}`$) between the two borders gives rise to the stable soliton trapped by the local defect with $`\lambda =0`$. In fact, in the region of relatively large values of $`\kappa `$, the lower border becomes “fuzzy”: as it is explained in the text, no systematic deviation from $`\theta =\pi /2`$ was found, but a deviation within a margin of $`\mathrm{\Delta }\theta =0.03\pi `$ is not ruled out, as the instability of unstable solitons becomes very slow, impeding to monitor their evolution up to a definite result. Initial solitons taken outside the stability region decay (in the region above the upper stability border, a residual small-amplitude soliton can be found, being flung away from the localized defect).

Fig. 4. Collision between a moving soliton with the fixed value $`\theta =\pi /2`$ and the defect. The lower and upper panels show, respectively, the evolution of the field $`|u(x,t)|`$, and the waveforms $`|u(x)|`$ (solid line) and $`|v(x)|`$ (dashed line) at the end of the simulation. (a) The defect with $`\kappa =0.2`$ captures the soliton with the velocity $`c=0.2`$, which is accompanied by conspicuous emission of radiation. (b) In the case of a stronger defect, with $`\kappa =0.6`$, less energy passes through, while larger shares of the energy are trapped at the defect and reflected back. (c) The case of the strongest defect, with $`\kappa =0.8`$. In this case, most energy is reflected back, while the amounts of energy trapped and transmitted through the defect are smaller.

Fig. 5. Regions in the parametric plane ($`c,\kappa `$), with fixed $`\theta =\pi /2`$, in which the moving soliton, respectively, passes the defect, is captured by it, or is split into three parts.

Fig. 6. The collision of the moving soliton with the defect. The parameters are $`c=0.075`$, $`\kappa =0.1,\theta =0.9\pi `$. The upper and lower panels have the same meaning as in Fig. 4. In this case, the soliton overcomes the repulsive barrier and passes through the defect, with a considerable loss of energy.

Fig. 7. A summary of results obtained for the interaction of the moving soliton and local defect, in the case when the soliton’s velocity is kept constant at $`c=0.075`$. In the region A, the defect is attractive for the solitons with small $`\theta `$, and the solitons is captured by it. In the region B, the defect is repulsive for the solitons with larger $`\theta `$; the soliton bounces back. In the region C, the defect is weak. The solitons passes through it, regardless of whether the interaction is attractive or repulsive. In the region D, the defect is strong, and, simultaneously, the soliton is heavy. The interaction combines repulsive and attractive features, and the soliton splits into three parts. One part is trapped, the remaining energy being scattered in forward and backward directions.


format: text

text:
# Topological Effects in Matrix Models representing Lattice Gauge Theories at Large N Talk delivered at TH-2002, International Conference on Theoretical Physics, Paris, UNESCO, July 22-27, 2002.

## 1 Introduction

The techniques producing numerical values for physical quantities affected by the strong interactions have continuously evolved over the last twenty years. They have become quite sophisticated and, coupled with major increases in affordable computing power, generate steady, incremental progress . However, if the sole activity in lattice field theory were that of adding digits of higher accuracy to some basic quantities one might get the impression that this subject, while useful, is rather dull. Luckily, this is not the case: there are lattice activities which do take risks and, at times, produce a breakthrough (e.g. lattice chiral fermions ), or, show that lattice field theory is good for other things than QCD (e.g. Higgs triviality bounds ). There are people who venture into supersymmetric territory , tackle “insoluble problems” like QCD with a chemical potential , use lattice techniques to check scenarios for explaining the baryon number asymmetry , and are inventive in many other ways. In turn, lattice field theory constructions enter mainstream particle physics  and even quantum gravity , providing novel “UV-completions”.

In this talk I shall present a project whose long term objective is to solve 4D planar QCD. The project is in its initial stages and, unlike more typical lattice field theory projects, is far from being guaranteed to succeed. This project is a collaborative effort with J. Kiskis and R. Narayanan . At present, the focus is narrowed down to the calculation of the $`Q^2`$-dependent meson-meson correlation functions in 4D. In 2D, where the planar limit is solvable (as shown by ’t Hooft ), this correlation function is of fundamental importance. Although ’t Hooft’s solution was mainly analytical, to actually get numbers also a numerical procedure is required, even in two dimensions. There is little chance for a similar analytical breakthrough in 4D. But, we hope that a numerical approach - this time containing a stochastic element - will be possible using old ideas of large $`N_c`$ reduction and new lattice fermions.

’t Hooft’s planar limit  in 4D requires three separate limits, one more than in ordinary QCD. One needs to take an ultraviolet cutoff, $`\mathrm{\Lambda }`$, (needed to define the theory beyond its perturbative series) to infinity, an infrared cutoff (typically the volume of an Euclidean torus), V, to infinity and the number of colors, $`N_c`$, at $`g^2N_c=\lambda `$ fixed, also to infinity ($`\lambda `$ is to be traded for a physical scale). The natural order of limits is as follows: first $`V\mathrm{}`$, next $`\mathrm{\Lambda }\mathrm{}`$ and last $`N_c\mathrm{}`$. One knows how to exchange the first two limits without altering the result even when the quarks are massless. Less is known about taking $`N_c\mathrm{}`$ earlier than last, but we shall do that. Otherwise, we are relegated to solving not only QCD first, but also QCD-like theories for different numbers of colors, and then the relevance of getting the planar limit numerically becomes mostly academic. The result  is nevertheless important for our project because such investigations tell us, quantitatively, how far $`N_c=\mathrm{}`$ is from $`N_c=3`$. For us, “solving” planar QCD means that a shortcut to the $`N_c=\mathrm{}`$ limit was used, and the method finds a numerical solution at $`N_c=\mathrm{}`$ with considerable more ease than would be possible in real QCD. This ancient hope sounds more realistic today than it did five years ago also thanks to progress made by string theorists . To be sure, in general, the interchange of the limits $`N_c\mathrm{}`$ with $`\mathrm{\Lambda }\mathrm{}`$ is non-trivial . For example, the theory may contain a set of degrees of freedom whose number does not increase with $`N_c`$, but whose impact on the ultraviolet behavior at any finite $`N_c`$ is large. Nor is the interchange of the $`V\mathrm{}`$ and $`N_c\mathrm{}`$ limits a trivial matter in a theory with massless elementary fields .

## 2 Large $`N_c`$ reduction

We start on a lattice with a “bare” coupling $`\beta =\frac{4}{g^2}`$ and at a finite volume $`V`$. There is very good evidence that the limit $`N_c\mathrm{}`$, with $`b\frac{\beta }{N_c}`$ fixed, exists. It exists for the entire range of positive values of $`b`$. At small values of $`b`$ (strong coupling) large $`N_c`$ reduction  tells us that the $`V\mathrm{}`$ limit is trivial: there is no $`V`$-dependence left at leading order in $`N_c`$; one can set $`V=1`$ and lattice momentum space is generated from the eigenvalues of the gauge field operators. Complications arise when one tries to extend this to weak coupling (large $`b`$) , but there are some options available . This is not dissimilar to the situation in models that admit a weakly coupled dual string description at large $`\lambda `$.

In their work, Eguchi and Kawai  proved that the Schwinger-Dyson equations for single-trace operators in the reduced and full models were the same at infinite $`N_c`$. Later  it was realized that this only established reduction for small enough $`b`$, and that to extend it to large $`b`$ an additional ingredient, for example, “quenching” , was needed.

In modern language, large $`N_c`$ reduction works by “field orbifolding” . This terminology comes from continuum field theories (actually these theories are the low energy sector of some string theories and the term originates from a method of constructing special backgrounds for the latter ), but the dependence on continuous space-time plays little role in the basic combinatorial identity. In our lattice application there is no explicit continuum. Field orbifolding can produce pairs of theories that admit ’t Hooft expansions and, at least at the level of bare diagrams, are simply related by possible rescalings of the coupling constant. The fields in one member of the pair are obtained by subjecting the fields in the other member to a set of constraints. These constraints reduce the symmetry group from $`U(N_c)\times 𝐆`$ to $`_iU(N_c^i)\times 𝐆^{}`$. In the context of large $`N_c`$ reduction one starts with a theory that has $`V=1`$, a single site lattice, and one ends up with a theory with a larger lattice (here chosen as a four dimensional symmetrical torus), $`V=L^4`$. The original theory has a $`U(N_c)\times U(1)^4`$ and the new theory has a symmetry $`\left[_{i=1,L^4}U(N_c/L^4)\right]\times U(1)^4\times Z(L)^4`$. (The $`N_c^i=N_c/L^4`$ and are assumed to be integers.) The two theories have the same $`b`$ and, in the second, the factor $`_{i=1,L^4}U(N_c/L^4)`$ is a lattice gauge symmetry group. The original theory is the reduced model and the main assumption of our project is the hope that the factor $`L^4`$ isn’t really needed to get to the planar limit.

Let me first illustrate the main trick in a simpler model, where the field integration measure, excepting the action term, is flat and unrestricted. This makes it easy to implement the constraint on the original fields by a “decoration” of the original planar diagrams. The “decoration” is a discrete gauge-field theory whose gauge link variables need to be summed over, living on the planar diagram of the original theory.

The model is known as the “Weingarten model” , and it defines a generating functional for random surfaces on a hypercubic lattice - a rudimentary string theory. The partition function of its reduced version  is given by

$$Z_{\mathrm{DW}}^N(b)=\underset{\mu =1}{\overset{d}{}}\underset{1i,jN}{}\left[_{\mathrm{}}^{\mathrm{}}𝑑A_{ij}^\mu 𝑑A_{ij}^\mu \right]e^{\frac{b}{N}_{\mu \nu }trA_\mu A_\nu A_\mu ^{}A_\nu ^{}_\mu trA_\mu A_\mu ^{}}$$
(1)
Although the integral does not converge, the formula is meaningful as a compact definition of an infinite series in $`b`$. Similarly, the large $`N`$ limit of $`F(b)\frac{1}{N^2}\mathrm{log}(Z_{\mathrm{DW}}^N(b))`$ is also meaningful as a series in $`b`$. It is given there by the sum of all planar $`|\mathrm{\Phi }|^4`$-type Feynman diagrams with a factor of $`b`$ attached to each vertex.

The reduced model has a symmetry $`SU(N)\times U(1)^d`$, where the $`SU(N)`$ acts by conjugation simultaneously on all of the $`A_\mu `$ matrices and the $`U(1)`$ factors independently affect the phases of each $`A_\mu `$ matrix. There also is a symmetry group made out of $`2^dd!`$ elements consisting of permutations of the $`\mu `$ indices and conjugations. This “crystallographic group” plays no role in our particular application here, and hence will be ignored (it also was ignored a few paragraphs earlier). The constrained, or, equivalently, space-time-extended model, is defined on a hypercubic toroidal lattice containing $`L^d`$ sites. Its partition functions is:

$`Z_{\mathrm{CDW}}^N^{}(b)={\displaystyle \underset{\mu =1}{\overset{d}{}}}{\displaystyle \underset{x_\mu =0}{\overset{L1}{}}}{\displaystyle \underset{1i,jN}{}}\left[{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑A_{ij}^\mu (x)𝑑A_{ij}^\mu (x)\right]`$
$`e^{\frac{b^{}}{N^{}}_{\mu \nu }_xtrA_\mu (x)A_\nu (x+\mu )A_\mu ^{}(x+\nu )A_\nu ^{}(x)_\mu _xtrA_\mu (x)A_\mu ^{}(x)}`$ (2)
The symmetry of the extended model is $`\left[_xU(N^{})\right]\times U(1)^d\times Z(L)^d`$. The first factor is a local lattice gauge symmetry, the second is a global symmetry and the third a space-time symmetry. Again we ignore the crystallographic group. Actually, the phases associated with the $`U(1)`$’s go only from zero to $`\frac{2\pi }{L}`$ \- larger phases can be reduced by lattice gauge transformations. Again the integral does not converge, but the formula is meaningful as a compact definition of an infinite series in $`b^{}`$. Similarly, the large $`N^{}`$ limit of $`F^{}(b^{})\frac{1}{L^d(N^{})^2}\mathrm{log}(Z_{\mathrm{CDW}}^N^{}(b^{}))`$ is meaningful as a series in $`b^{}`$. It is again given by the sum of all planar $`|\mathrm{\Phi }|^4`$-type Feynman diagrams with a factor of $`b^{}`$ attached to each vertex, but, with propagators that are less trivial than before.

The extended model can be viewed as a constrained version of the reduced model if we pick $`N=N^{}L^4`$, $`b=b^{}`$, write the $`1iN`$ indices as double indices $`(i^{},x)`$ with $`1i^{}N^{}`$, $`0x_\mu L1`$ and restrict the $`[A_\mu ]_{(i^{},x),(j^{},y)}`$ to be zero when $`yx+\widehat{\mu }`$. When $`y=x+\widehat{\mu }`$ there are no additional restrictions. If all $`A_\mu `$’s in the reduced model of eq. (1) obey this constraint, the portion of the trace referring to the $`x,y`$ components of the indices can be written out explicitly and one obtains the extended model of equation (2). An important observation is that the constraint can be viewed as enforcing invariance under a subgroup of the symmetry group of the reduced model:

$`A_\mu =𝒯_\nu (A_\mu )=e^{\frac{2\pi ı}{L}\delta _{\mu ,\nu }}g_\nu ^1A_\mu g_\nu ,\nu =1,..,d`$ (3)
$`\left[g_\nu \right]_{(i^{},x),(j^{},y)}=\delta _{i^{},j^{}}\delta _{x,y}e^{\frac{2\pi ı}{L}y_\nu }`$ (4)
The set of operations $`𝒯_\nu `$ generates a $`Z(L)^d`$ subgroup of the symmetry group of the reduced model. The elements of the group are labeled by $`[s]`$ where $`[s]=(s_1,s_2,\mathrm{},s_d)`$ and $`s_\mu =0,1,\mathrm{},L1`$.

$$𝒯^{[s]}(A_\mu )=e^{\frac{2\pi ı}{L}s_\mu }\left(g^{[s]}\right)^{}A_\mu g^{[s]},g^{[s]}=\underset{\nu =1}{\overset{d}{}}g_\nu ^{s_\nu }$$
(5)
Note that the action of $`𝒯^{[s]}`$ on $`A_\mu `$ differs from its action on $`A_\nu `$ for $`\mu \nu `$ by an overall phase factor. This $`Z(L)^d`$ intersects both original symmetry group factors, $`SU(N)`$ and $`U(1)^d`$. Invariance under the action of the four generators $`𝒯_\nu `$ is equivalent to invariance under the entire $`Z(L)^d`$ group.

Since we integrate freely on the entries of the matrices $`A_\mu `$, we can redefine the constrained model by keeping full matrices in the quadratic term (where the constrained and unconstrained entries of the $`A_\mu `$-matrices do not mix) and only matrices satisfying the constraint in the quartic term. The constraint can be put in effect by first expanding in the exponential of the quartic factor, and then replacing each occurrence of a constrained $`A_\mu `$, independently, by $`𝒯^{[s]}(A_\mu )`$ with an unconstrained $`A_\mu `$, and averaging each of the many $`𝒯^{[s]}`$ factors over the group $`Z(L)^d`$. Only the parts of $`A_\mu `$ that are invariant, making up the matrix $`A_{\mu \mathrm{constr}}`$, survive this group averaging, while the other entries are set to zero. We now leave this group averaging to the end of the calculation, and, instead, do the Gaussian integrals first, producing Feynman diagrams. The propagators are now given by

$$<𝒯^{[t]}(A_\mu )_{ij}𝒯^{[s]}(A_\nu ^{})_{kl}>=e^{\frac{2\pi ı}{L}(t_\mu s_\mu )}g_{il}^{[st]}g_{kj}^{[ts]}\delta _{\mu \nu }$$
(6)
$`𝒯^{[t]}`$ and $`𝒯^{[s]}`$ appear only in this propagator. Moreover, as we see, there is only dependence on $`[ts]`$. Hence, the independent averages over $`[t]`$ and $`[s]`$ can be replaced by a single average over the combination $`[ts]`$. It is only this extra group averaging procedure that distinguishes the two models.

The sum over the $`i,j,\mathrm{}`$ indices produces a factor of $`N`$ for each index loop in the double line notation for propagators in the unconstrained model. In the extended model, because of the constraints, these index sums produce a trace of the product of the various $`g^{[s]}`$ factors encountered round the index loop. The $`g^{[s]}`$ matrices also form a group, $`Z(L)^d`$, but this group is distinct from the $`Z(L)^d`$ made out of the $`𝒯^{[s]}`$ operations. Thus, for every index loop we need to calculate the trace of a $`g^{[t]}`$ that is the parallel transporter round the loop.

From their definition we see that the $`g^{[s]}`$ obey:

$$trg^{[s]}=N\delta _{[s],[0]}$$
(7)
where the delta function forces the group element to be unity. This has the important consequence that the phase factors that appear in the definition of the $`𝒯^{[s]}`$ operations cancel out any time the contribution from the index loops is non-vanishing. (Note that the phase factors are associated with propagators, not index lines.) Thus, we can ignore the phase factors and take the propagators as given by

$$<𝒯^{[t]}(A_\mu )_{ij}𝒯^{[s]}(A_\nu ^{})_{kl}>g_{il}^{[u]}g_{kj}^{[u]}$$
(8)
Group averaging amounts to averaging over all $`L^d`$ values of the $`[u]`$ labels.

More generally, this induces us to think about a situation in which we assign to each propagator in the planar diagram a group element $`gG`$ where $`G`$ is some finite group. $`g`$ is represented by an $`N\times N`$ matrix $`D_{ij}(g)`$ which replaces one of the two index lines making up the propagator. Before the “decoration”, this index line was a Kronecker delta. The other index line of the same propagator is decorated by $`D_{kl}(g^1)`$. The group element associated with each propagator is independently averaged over all of $`G`$ and the decoration introduces a weighting factor given by the product over all index loops (faces) of traces of products of $`D(g)`$ matrices round the loop. The summations over $`G`$ are now easy to perform - a standard exercise in lattice strong coupling expansion, this time defined on a random two dimensional lattice .

For any $`hG`$, $`trD(h)`$ is a class function and hence expandable as:

$$trD(h)=\underset{r}{}c_rd_r\chi _r(h)c_r=\underset{hG}{}\chi _r(h^1)\frac{1}{d_r}trD(h)$$
(9)
The averaging over the $`g`$’s associated to the propagators produce a total factor given by

$$\underset{r}{}c_r^fd_r^2=\underset{r}{}c_r^{2+v}d_r^2$$
(10)
where $`f`$ is the number of faces and $`d_r`$ the dimension of the $`r`$-irreducible representation of $`G`$. In the above equation we make use of the planarity of the diagram. If $`e`$ is the number of edges in the diagram and $`v`$ the number of vertices, we have:

$$fe+v=2e=2vf=2+v$$
(11)
In order to be able to absorb the decoration factor into a rescaling of $`b^{}`$ we need $`c_r=K`$ for all $`r`$. From (9), $`K=\frac{N}{d(G)}`$, where $`d(G)`$ is the number of elements in $`G`$. $`d(r)`$ is the dimensionality of the irreducible representation $`r`$, and $`d(G)=_rd_r^2`$. The contribution of the diagram is $`\left(\frac{b^{}}{N^{}}\right)^v_rc_r^fd_r^2`$. $`K`$ is an integer and we see that the representation of $`G`$ furnished by $`D(h)`$ contains $`K`$ copies of the regular representation. This is a consequence of the requirement to have an overall $`b^v`$ factor. It is easy to check that in our case, where $`G=Z(L)^d`$ (with elements $`g^{[s]}`$), this condition on the representation holds.

Turning now to our special case, we set $`d(G)=L^d`$, $`N^{}=\frac{N}{d(G)}`$ and $`b^{}=b`$. The weight of the above diagram becomes $`N^2b^v\frac{1}{d(G)}`$. From the definition of the free energy of the reduced and extended model we now conclude that

$$F^{}(b)=\frac{d(G)}{N^2}\mathrm{log}(Z_{\mathrm{CDW}}^N^{}(b))=F(b)$$
(12)
The equality of the free energies is the main result of reduction. Relations between other correlation functions can be obtained by adding more terms to the action and taking derivatives with respect to the new, infinitesimal, couplings. The above result was obtained a long time ago  by direct inspection of the series and by mapping the Schwinger-Dyson equations of single-trace operators from one model to the other.

In its extended version, the Weingarten model will turn into a lattice gauge theory if we constrain the $`A_\mu (x)`$ matrices to be simple-unitary and change the integration measure to the $`SU(N)`$ group Haar measure. The relation between the reduced and extended versions holds even after the unitarity constraint has been imposed on both models. This was proven a long time ago, first by the Schwinger-Dyson equation method  and, subsequently, by inspection of lattice strong coupling diagrams  and also using stochastic quantization . A proof by field orbifolding would require a generalization of the analysis above, but since the result is known to hold, there is little doubt that such a generalization exists. Field orbifolding admits generalizations involving nonabelian global symmetries, and it would be interesting to see whether previously unknown results in planar lattice gauge theory could be obtained from this approach.

## 3 Prototype model and its extended variants

The unconstrained, reduced, gauge model consists of 4 matrices $`U_\mu SU(N_c)`$, 4 unitary matrices $`V_\mu SU(N_f)`$ and two rectangular Grassmann matrices $`\overline{\mathrm{\Psi }},\mathrm{\Psi }`$, where one index is in the fundamental representation of $`SU(N_c)`$ (color), and the other index in the fundamental representation of $`SU(N_f)`$ (flavor). The fermion matrices also have a Dirac index, which we often suppress. The model is defined by:

$`S=\beta ^u{\displaystyle \underset{\mu >\nu }{}}C_{\mu \nu }^uC_{\mu \nu }^u+\beta ^v{\displaystyle \underset{\mu >\nu }{}}C_{\mu \nu }^vC_{\mu \nu }^v+Tr\left[\overline{\mathrm{\Psi }}D(U,V)\left(\mathrm{\Psi }\right)\right]`$
$`C_{\mu \nu }^u[U_\mu ,U_\nu ],C_{\mu \nu }^v[V_\mu ,V_\nu ]`$ (13)
The linear operator $`D(U,V)`$ is a reduced version of a lattice Dirac operator, to be introduced later on.

The constrained associated model lives on an $`L^4`$ lattice . The constraints mean that the matrices $`U_\mu ,V_\mu ,\overline{\mathrm{\Psi }},\mathrm{\Psi }`$ are expressible in terms of $`L^4`$ smaller matrices $`u_\mu (x),v_\mu (x),\overline{\psi }(x),\psi (x)`$. The dimensions are $`n_c`$ instead of $`N_c`$ and $`n_f`$ instead of $`N_f`$ where $`n_{c,f}=N_{c,f}/L^4`$. ($`n_{c,f},N_{c,f},L`$ are all integers.)

$`\left[U_\mu \right]_{(ix),(jy)}=[u_\mu (x)]_{ij}\delta _{y,x+\mu },\left[V_\mu \right]_{(ax),(by)}=[v_\mu (x)]_{ab}\delta _{y,x+\mu }`$
$`\left[\mathrm{\Psi }\right]_{(ix),(ay)}=[\psi (x)]_{ia}\delta _{xy},\left[\overline{\mathrm{\Psi }}\right]_{(ay),(ix)}=[\overline{\psi }(x)]_{ai}\delta _{yx}`$ (14)
For small lattice gauge coupling the constrained and unconstrained models will be equivalent in the limit $`N_{c,f}\mathrm{}`$ with $`\frac{N_f}{N_c}`$ held fixed and a simple $`D(U,V)`$. But, for larger coupling the equivalence breaks down and we shall adopt here one of the proposals for keeping the equivalence up even at large lattice couplings, namely quenching .

## 4 Quenching

Quenching amounts to freezing the eigenvalues of $`U_\mu ,V_\mu `$ to sets that are uniformly distributed round the circle. Other than that the unitary matrices are free. One also takes the rescaled coupling, $`\frac{\beta _f}{N_c}`$, to infinity so as to freeze out the flavor gauge group . The $`U_\mu `$ eigenvalues play the role of lattice momenta in boson propagator lines, while the $`V_\mu `$ eigenvalues enter only into the fermionic lines. The operator $`D(U,V)`$ breaks up into $`N_f`$ blocks of size $`4N_c\times 4N_c`$ each. This is the crucial point of our approach: The numerical price of handling fermions has a chance to simplify tremendously if indeed $`N_c`$ does not need to be anywhere nearly as big as the parameter $`n_cL^4`$ typical to numerical QCD simulations and, nevertheless, the planar limit stays well approximated.

## 5 Overlap fermions

We now define $`D`$ for one of the flavor blocks. So long we are focusing only on fermions we can absorb the flavor phase into the color matrices $`U_\mu `$ (whose eigenvalues are quenched as before). Define the unitary operators $`X_\mu `$:

$$X_\mu =\frac{1\gamma _\mu }{2}U_\mu +\frac{1+\gamma _\mu }{2}U_\mu ^{}$$
(15)
The Wilson-Dirac matrix $`D_W`$ is a given by

$$D_W=m+\underset{\mu }{}(1X_\mu )$$
(16)
We choose the parameter $`m`$ to be about $`1`$ (it must be between $`2`$ and $`0`$). It is easy to see that $`H_W\gamma _5D_W`$ is hermitian.

The overlap construction says that the chiral determinant line bundle is represented on the lattice by the inner product of two states, $`+|`$. Here $`+|`$ is the eigen-bra of a system of noninteracting second quantized fermions in a single body potential given by $`H_W`$ with the mass sign switched, while $`|`$ is the eigen-ket of a system of noninteracting second quantized fermions in a single body potential given by $`H_W`$ with the mass chosen as above. One can change the lattice definition without affecting the continuum limit so that the mass in $`H_W`$ is taken to be $`+\mathrm{}`$ for the $`+|`$-state. In that case the $`H_W`$ for the $`+|`$-state can be replaced by $`\gamma _5`$ .

The overlap $`+|`$ can be written in terms of single particle wave-functions as the determinant of a matrix $`A`$ defined below . Choose the chiral representation for the $`\gamma _\mu `$ matrices, so that

$$\gamma _5=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$$
(17)
Let $`\mathrm{\Phi }`$ be a simple unitary matrix of eigenvectors of $`H_W`$, ordered so that

$$\mathrm{sign}(H_W)\mathrm{\Phi }=\mathrm{\Phi }\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),\mathrm{\Phi }=\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)$$
(18)
Then

$$+|=detA$$
(19)
For a vector-like theory we need $`|det(A)|^2`$, which, unlike $`detA`$, is a function, not a section of a bundle. Simple linear algebra tells us that

$$|det(A)|^2=det(A)det(D)$$
(20)
However,

$$\gamma _5\mathrm{sign}(H_W)\mathrm{\Phi }=\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)$$
(21)
Hence,

$$det\left[\frac{1+\gamma _5\mathrm{sign}(H_W)}{2}\mathrm{\Phi }\right]=det\left[\frac{1+\gamma _5\mathrm{sign}(H_W)}{2}\right]=|+||^2$$
(22)
We discover that chiral symmetry can be preserved on the lattice if we use the overlap Dirac operator :

$$D_o=\frac{1+\gamma _5\mathrm{sign}H_W}{2}$$
(23)
The operator $`D`$ is constructed out of $`D_o`$, with the addition of a small mass term for the quarks.

## 6 Topology and Anomalies

Anomalies are expressed by particular Berry-phase factors innate to the $`|`$ state of the overlap . They obstruct smooth gauge invariant sections in cases where anomalies do not cancel.

Topology is expressed by $`H_W`$ having unequal numbers of positive and negative eigenstates (above, where $`A`$ is taken as a square matrix, we tacitly assumed these two numbers to be equal but this need not be so, in which case what we called “$`detA`$” is identically zero). One can easily exhibit explicit sets of $`U_\mu `$ matrices for which there is non-zero topology.

Problems arise if $`H_W`$ has an eigenvalue that is exactly zero. But this will not happen as the continuum limit is approached because one can prove the following inequality :

$$[U_\mu ,U_\nu ]ϵ\frac{\psi |H_W^2|\psi }{\psi |\psi }1|1+m|\delta (ϵ)$$
(24)
where $`\delta (ϵ)`$ goes to zero as $`ϵ`$ goes to zero.

At infinite $`N_c`$, for large enough lattice couplings and a standard lattice action, the commutators indeed become bounded , as required.

## 7 Some history

The first treatment of fermions in reduced models was carried out by Gross and Kitazawa  who invented the momentum “force-feeding” prescription that was derived above as a remnant of the frozen-out flavor gauge sector . These authors also invented a version of reduction that used hermitian matrices, $`A_\mu `$, rather than the unitary link matrices one has on a lattice . The Gross-Kitazawa model can be obtained from the lattice model by approximating $`U_\mu 1+ıA_\mu `$. However, as also noted by Gross and Kitazawa, this approximation does away with anomalies and topology. For many years we did not known how to cure this problem on the lattice, where the link matrices are unitary. The problem was the same as the notorious problem of lattice chirality. As we have seen, recent progress on the chiral fermion problem  makes it now possible to preserve all the necessary fermionic properties also in reduced lattice models removing the main obstacle of principle to making nontrivial use of reduction in numerical simulations.

## 8 Acknowledgments

I am grateful to Laurent Baulieu and to Daniel Iagolnitzer for the invitation to TH-2002. I wish to acknowledge partial support at the Institute for Advanced Study from a grant in aid from the Monell Foundation as well as partial support by the DOE under grant DE-FG02-01ER41165 at Rutgers University.