## Abstract

We introduce a system with one or two amplified nonlinear sites (‘hot spots’, HSs) embedded into a two-dimensional linear lossy lattice. The system describes an array of evanescently coupled optical or plasmonic waveguides, with gain applied to selected HS cores. The subject of the analysis is discrete solitons pinned to the HSs. The shape of the localized modes is found in quasi-analytical and numerical forms, using a truncated lattice for the analytical consideration. Stability eigenvalues are computed numerically, and the results are supplemented by direct numerical simulations. In the case of self-focusing nonlinearity, the modes pinned to a single HS are stable and unstable when the nonlinearity includes the cubic loss and gain, respectively. If the nonlinearity is self-defocusing, the *unsaturated* cubic gain acting at the HS supports *stable* modes in a small parametric area, whereas weak cubic loss gives rise to a *bistability* of the discrete solitons. Symmetric and antisymmetric modes pinned to a symmetric set of two HSs are also considered.

## 1. Introduction

Modes of fundamental significance to nonlinear optics [1,2] and plasmonics [3–12] are dissipative spatial solitons that result from the simultaneous balance among diffraction, self-focusing nonlinearity, loss and compensating gain. Stability is a crucially important issue in the theoretical analysis of dissipative solitons. An obvious necessary condition for the stability of localized modes is the stability of the zero background around them. The basic complex Ginzburg–Landau (CGL) equation, which includes the bandwidth-limited linear gain and nonlinear loss acting on a single field, is unable to produce stable dissipative solitons, because the action of the linear gain on the zero background makes it unstable. On the other hand, dissipative solitons can be fully stabilized in systems of linearly coupled CGL equations [13,14] that model dual-core waveguides, with the linear gain and loss acting in different cores [11,15–19], including the -symmetric version of the system that features the balance between the gain and loss [20,21]. Stable solitons can also be generated by a single CGL equation with cubic gain ‘sandwiched’ between linear and quintic losses, which may be realized in optics as a combination of linear amplification and power-dependent absorption [22–32].

Another method for creating stable localized modes makes use of linear gain applied at a ‘hot spot’ (HS, i.e. a localized amplifying region embedded into an ordinary lossy waveguide [33–35] or a Bragg grating [36]). Models with multiple HSs [37–41], and similar extended amplifying structures [42,43], have also been studied. HSs can be built by implanting gain-producing dopants into a narrow segment of the waveguide [44], or, alternatively, by focusing an external pump beam at the designated position of the HS in a uniformly doped waveguide. In addition to models with the localized direct (phase-insensitive) gain, systems including the localization of parametric gain were also developed [45].

Dissipative solitons can be stably pinned to the HS owing to the balance between the local gain and uniform loss in the bulk waveguide. For narrow HSs, modelled by the delta-functional distribution of the gain, solutions for the pinned dissipative solitons are available in an analytical form [33,37,41]. Furthermore, models with mutually balanced gain and loss applied in the form of *δ*-functions at separated points [46], or at a single location, in the form of a gain–loss dipole described by the derivative of the *δ*-function [47], make it possible to find -symmetric solitons pinned to these points. Other one- and two-dimensional HS-pinned modes, including stable vortices fed by the gain confined to an annular-shaped area [48–52], can be found numerically [34,35,38–40].

While dissipative solitons in uniform media are always unstable against the blowup in the absence of the higher-order (quintic) nonlinear losses [53,54], it is worth noting a counterintuitive result [35] demonstrating that *stable* dissipative localized modes in uniform linearly lossy media may be supported by *unsaturated* localized cubic gain alone. Stable dissipative solitons were also predicted in a setting that combines the uniformly distributed linear gain in the *D*-dimensional space and nonlinear loss growing from the centre to periphery faster than *r*^{D}, where *r* is the radial coordinate [55].

The class of models with the localized gain includes lattice systems. In [56], the one-dimensional model was introduced for a linear lossy lattice with a single or two amplified (active) sites embedded into it, which represent the HSs in the discrete setting. It was assumed that the nonlinearity was carried solely by the same active sites. This system, which may be considered as a variety of discrete CGL equations [57–67], can be implemented in the experiment using arrays of optical waveguides [68] or arrayed plasmonic waveguides [69–71]. In particular, it suggests possibilities for selective excitation of particular core(s) in the arrayed waveguides, if the system is uniformly doped, but only the selected cores are pumped by an external laser beam. In [56], discrete solitons pinned to the HS in the lattice system were found in an analytical form, similar to the soliton solutions available in the discrete linear Schrödinger equation with embedded nonlinear elements [72–74], and the stability of the localized models was investigated by means of numerical methods.

This work aims to extend the analysis for the two-dimensional lattice system, with one or two active nonlinear sites embedded into the linear lossy bulk lattice. The experimental realization of such a two-dimensional medium is also possible in nonlinear optics, using waveguiding arrays permanently written in bulk silica [75,76]. An essential distinction from the one-dimensional counterpart [56] mentioned above is that two-dimensional localized lattice modes cannot be found analytically, even when only a single nonlinear site is embedded into the linear matrix. Nevertheless, we demonstrate that semi-analytical solutions can be obtained for truncated (finite-size) lattices.

The paper is organized as follows. The discrete two-dimensional CGL equation is introduced in §2. Section 3 presents semi-analytical results for the truncated lattices. Results of the linear-stability analysis for the HS-pinned lattice solitons against small perturbations are reported in §4. In §5, we extend the stability analysis, considering the crucially important issue of the onset of the zero-solution instability. A brief consideration of the double-HS system is presented in §6. The paper is concluded by §7.

## 2. The model

As said above, we consider the two-dimensional generalization of the one-dimensional lattice model introduced in [56]:
2.1
where *m*,*n*=0,±1,±2,… are discrete coordinates on the lattice, *δ*_{m,0} and *δ*_{n,0} are the Kronecker symbols, and the coefficient of the linear coupling between adjacent cores is scaled to unity. Further, *γ*>0 is the linear loss in the bulk lattice, *Γ*_{1}>0 and *Γ*_{2} account for the linear gain and linear potential applied at the HS site (*m*=*n*=0), whereas *B* and *E* account for the Kerr nonlinearity and nonlinear loss/gain (for *E*>0/*E*<0) acting at the HS.

In optics, the two-dimensional discrete equation, as well as its one-dimensional counterpart, can be derived by means of well-known methods [57–68,77]. In the application to two-dimensional arrays of plasmonic waveguides, which can be built as a set of metallic nanowires embedded into a bulk dielectric [69–71], equation (2.1) can be derived in the adiabatic approximation, with the exciton field eliminated in favour of the photonic one. It is relevant to mention that the well-known *staggering transformation* [68], , where the asterisk stands for the complex conjugate, simultaneously reverses the signs of *Γ*_{2} and *B*, thus rendering the self-focusing and defocusing signs of the nonlinearity, which correspond to *B*>0 and *B*<0, respectively, mutually convertible. This circumstance is essential, in particular, for modelling arrays of plasmonic waveguides, where the intrinsic excitonic nonlinearity is always self-repulsive. Below, we fix the signs of *Γ*_{2} and *B* by setting *Γ*_{2}>0, i.e. the corresponding linear potential at the HS is *attractive*, whereas *B* may be positive (self-focusing), negative (self-defocusing) or zero. Unless *B*=0, this coefficient is normalized to be *B*=±1. These two cases are considered separately below, along with the case of *B*=0, when the nonlinearity is represented solely by the cubic dissipation localized at the HS.

## 3. The analysis for the truncated lattice

We seek solutions to equation (2.1) for stationary localized modes with real propagation constant *k* as
3.1
Outside the HS site, *m*=*n*=0, equation (2.1) gives rise to the linear stationary equation
3.2
which, unlike its one-dimensional counterpart (cf. [56]), does not admit exact analytical solutions. An approximate solution can be constructed on a truncated lattice. The simplest version of the truncation is shown in figure 1, where four independent amplitudes obeying equation (3.2) are defined: *U*_{0} (at the centre), *U*_{1} (in the first rhombic layer surrounding the centre) and *U*_{2,3} (two independent amplitudes in the second rhombic layer). The amplitudes in all other layers are neglected: *U*_{4,5,…}=0. At (*m*,*n*)≠(0,0), the so truncated equation (3.2) yields
3.3
where central amplitude *U*_{0} is defined to be real, and the complex coefficient is
3.4

The remaining nonlinear equation (2.1) at the HS site, *m*=*n*=0, reduces to a complex equation relating real peak intensity and propagation constant *k*:
3.5
which can be solved numerically. Subsequently, the amplitudes in the surrounding layers are obtained from equation (3.3).

The accuracy of the truncated-lattice analysis can be improved by including more rhombic layers in the calculation. Although the respective versions of equations (3.3) and (3.5) become more complicated, they are still easy to solve numerically. Table 1 shows numerically found *U*_{0} and *k* for the different truncation settings at two different sets of parameters. Naturally, convergence of the solutions for the localized modes is observed with the increase of the number of rhombic layers.

## 4. The linear-stability analysis

The stability of the pinned modes found as outlined above was studied by means of the linearization procedure [78]. To this end, perturbed solutions were taken as
4.1
where *V* _{m,n}(*z*)=*X*_{m,n}(*z*)+i*Y* _{m,n}(*z*) is a complex-valued perturbation with infinitesimal amplitude *ϵ*. Substituting this, along with *U*_{m,n} split into the real and imaginary parts as per equation (3.2), into equation (2.1) results in a linear system that governs the evolution of perturbations *X*_{m,n} and *Y* _{m,n}:
4.2
The eigenvalue problem is obtained by substituting and in equations (4.2), the underlying stationary mode *u*_{m,n}(*z*) being linearly stable if all eigenvalues λ have Re(λ)≤0. For the use in the stability analysis, modal amplitudes *U*_{m,n} were found numerically by solving the full system obtained from equation (2.1) by the substitution of expression (3.1), rather than from the truncated-lattice approximation, although the difference is very small. The numerical solution was carried out with periodic boundary conditions, imposed onto a domain of a size which was much larger than the width of any soliton considered below.

Figure 2 shows a typical stable mode and its stability-eigenvalue spectrum. The mode is peaked at the HS and symmetric about it (the shape of this mode and its propagation constant are very close to their counterparts produced by the truncated-lattice approximation).

### (a) The self-focusing nonlinearity: *B*=+1

To develop a systematic analysis, we first address the stability of pinned modes in the self-focusing regime (*B*=1). Figure 3 shows the solution branches as functions of the localized linear gain *Γ*_{1} at different values of the cubic dissipation *E*. Stable solution branches can be found only at *E*>0, i.e. in the presence of the cubic loss. For instance, at *E*=0.1, the pinned modes with peak amplitudes *U*_{0,0}>1.277 are linearly stable, and unstable at *U*_{0,0}<1.277. Naturally, the stable and unstable modes belong to portions of the branches where the amplitude is, respectively, a growing or decreasing function of *Γ*_{1}. The transition between the curves which contain the unstable segment and those which do not contain it happens, in the case displayed in figure 3, at *E*=*E*_{cr}≈0.18.

For the parameters considered here, stable and unstable solutions exist at *Γ*_{1}≥0.7675, whereas at *Γ*_{1}<0.7675 any initial excitation decays into the zero solution, as there is not enough energy input to support non-trivial modes. Figure 4 shows some typical examples of such stable and unstable solutions. The stability is greatly enhanced by increasing the magnitude of the cubic loss, *E*. At *E*=1, all the pinned modes are stable, even at very large values of *Γ*_{1}. On the other hand, in the absence of the cubic dissipation (*E*=0), or in the presence of the cubic gain (*E*<0), all the modes are unstable. These unstable modes have essentially the same profile as those in figure 4, and therefore they are not shown here.

Figure 3 shows that all the solution branches emerge from the critical value of the linear gain, *Γ*_{1}≈1.152, which is explained below. Furthermore, the unstable branch corresponding to *E*=0 approaches a vertical asymptote at *Γ*_{1}=0.5, near which the peak amplitude increases drastically. Finally, figure 5 shows typical evolution of peak amplitudes of the unstable modes corresponding to different values of *E*, as obtained from full simulations of equation (2.1). In the case of *E*>0, the cubic loss stabilizes the system, and therefore the unstable mode evolves into a stable one existing at the same *Γ*_{1} (figure 3). However, at *E*≤0 (the cubic gain, or zero loss), unstable modes quickly blow up, as there are no stable branches that might serve as attractors in this case.

### (b) The self-defocusing nonlinearity: *B*=−1

Figure 6 shows solution branches in the self-defocusing regime, with *B*=−1. In the presence of a small cubic loss, the solution branch exhibits a *bistability* for a certain range of values of linear gain *Γ*_{1}. For instance, at *E*=0.1, stable solutions with different amplitudes coexist in the region of 1.152≤*Γ*_{1}≤1.319. The two stable branches are connected by an unstable one, whose modal profile and stability spectrum are similar to those found in the self-focusing case, see figure 4*a*,*b*, and therefore they are not displayed here. An example of the bistability is shown in figure 7. Generally, the mode with the smaller peak amplitude is broader. When *E* increases, the unstable branch eventually gets stabilized by the strong cubic loss and disappears from the bifurcation diagram, which happens, in the case shown in figure 7, at *E*=*E*_{cr}≈0.58. At *E*>*E*_{cr}, all solutions are stable, and the bistability does not take place.

A remarkable finding is that self-defocusing gives rise to stable modes even in the absence of the cubic loss, and in the presence of weak cubic gain (*E*≤0). For small values of the cubic gain, bistability similar to that presented in figure 8 is observed. When the cubic gain grows (i.e. *E* is getting more negative), the stable branch with the larger peak amplitude disappears. All the solutions become unstable at still larger strengths of the cubic gain. Figure 9 shows typical examples of the evolution of the peak amplitudes of unstable modes corresponding to different values of *E*, as produced by simulations of equation (2.1). In the case of *E*>0, the initial unstable mode evolves into the closest stable one available at the same value of *Γ*_{1}. As concerns unstable modes found at *E*<0, they blow up rather than evolve into stable modes, if the latter ones exist in the case under consideration.

### (c) The case of *B*=0

Finally, we address the stability of the pinned modes in the case of *B*=0, when the nonlinearity is represented solely by the cubic loss or gain. Figure 10 shows the stability of solution branches in this situation. In particular, the solutions are always unstable in the presence of the cubic gain (*E*<0), and are always stable in the presence of the cubic loss (*E*>0). If the cubic dissipation vanishes as well, i.e. *E*=0, equation (2.1) becomes linear, admitting a single stable solution at a specific value of the HS gain, which is *Γ*_{1}=1.152 for the chosen parameters, that provides for the compensation of the background loss. This solution obviously has an arbitrary amplitude, which explains why the corresponding branch is vertical in figure 10. Further, figure 11 shows examples of stable and unstable modes and their (in)stability spectra. Unstable modes in the case of *B*=0 evolve into the zero solution rather than blowing up, in spite of the action of the local cubic gain, in that case.

## 5. Onset of instability of the zero solution

Here, we again address the linearized version of equation (2.1), setting *B*=*E*=0, and study the stability of the zero solution around the HS, which is an obvious condition necessary for the background stability of pinned modes in the nonlinear system. The onset of the background instability exactly corresponds to the existence of the above-mentioned pinned-mode solution to the linearized equation, which implies the equilibrium between the bulk loss (∼*γ*) and local gain (∼*Γ*_{1}), in the presence of the HS's attractive potential (∼*Γ*_{2}). Unlike the only counterpart of the present model [56], for the two-dimensional lattice, this solution of the linear equation cannot be found in an analytical form. However, the truncated-lattice approximation presented in §3 provides an efficient way to study the onset of the background instability. In particular, the onset is determined by the critical value of *Γ*_{1} that makes equation (3.5) with *B*=*E*=0 solvable. These values, found for different truncated configurations, are summarized in table 2. As the number of rhombic layers included in the calculations increases, the critical value of *Γ*_{1} slowly converges to the above-mentioned numerically found value, *Γ*_{1}≈1.152, from which the solution branches emanate in figures 3, 6 and 10.

## 6. Dual hot spots: symmetric and antisymmetric modes

Finally, we briefly consider pinned modes in the lattice with two identical HSs. The truncated-lattice approximation was first used to derive simplified equations, such as equation (3.5), to describe this configuration. To illustrate the situation, we place two HSs at *n*=0 and *m*=±1. Pinned modes supported by the dual HS can be naturally classified as symmetric, antisymmetric and asymmetric (provided that the latter species exists) [72–74]. While amplitudes of the solutions at the two HSs may be assumed real without loss of generality, the symmetric and antisymmetric modes are those with *U*_{−1,0}=*U*_{1,0} and *U*_{−1,0}=−*U*_{1,0}, respectively. The simplest truncated-lattice approximations, with a single rhombic layer surrounding the two HSs, are shown in figure 12. For the symmetric mode (figure 12*a*), a calculation similar to that presented in §3 yields
6.1
with *K* given by equation (3.4), whereas the antisymmetric mode (figure 12*b*) gives
6.2

Table 3 presents numerically found values of the propagation constant, *k*, and HS amplitude, *U*_{1}, for both the symmetric and antisymmetric modes. Figure 13 shows that stable pinned modes, both symmetric and antisymmetric, can be supported by the dual HS. Peak amplitudes of these modes (see the figure caption) are in good agreement with the prediction of the truncated-lattice approximation.

## 7. Conclusion

We have introduced the two-dimensional discrete dynamical system based on the bulk linear lossy lattice, into which one or two nonlinear sites with linear gain (HSs) are embedded. The system can be implemented in the form of an array of optical or plasmonic waveguides, admitting, in particular, selective excitation of individual cores, by the local application of the external pump to the uniformly doped array. The analysis of localized modes pinned to the HSs was developed both semi-analytically (using truncated lattices) and numerically. It has been found that, in the case of the self-focusing nonlinearity at the single HS, the pinned modes are stable (unstable) when the nonlinearity acting at the HS contains the cubic loss (gain). On the other hand, it is worthwhile to note that, in the case of the self-defocusing nonlinearity, the HS with an unsaturated cubic gain supports stable modes in a rather narrow parameter region, and bistability occurs in the case of weak cubic loss. Symmetric and antisymmetric modes pinned to dual HSs were also discussed briefly.

A challenging issue for the subsequent analysis is the possibility of the existence of asymmetric modes supported by the symmetric dual HSs. On the other hand, a natural extension may be to embed one or two HSs into a nonlinear lossy lattice.

## Funding statement

Partial financial support was provided by the Hong Kong Research Grants Council through General Research Fund contract no. HKU 711713E.

## Footnotes

One contribution of 19 to a Theme Issue ‘Localized structures in dissipative media: from optics to plant ecology’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.