## Abstract

We study nonlinear effects in one-dimensional (1D) arrays and two-dimensional (2D) lattices composed of metallic nanoparticles with the nonlinear Kerr-like response and an external driving field. We demonstrate the existence of families of moving solitons in 1D arrays and characterize their properties such as an average drifting velocity. We also analyse the impact of varying external field intensity and frequency on the structure and dynamics of kinks in 2D lattices. In particular, we identify the kinks with positive, negative and zero velocity as well as breathing kinks with a self-oscillating profile.

## 1. Introduction

Within decades of fruitful developments, the field of plasmonics has become a prominent area of research with applications ranging from super-imaging [1] to optical sensors [2], nonlinear devices [3] and photovoltaics [4]. Combining strong surface plasmon resonances, high intrinsic nonlinearities and deep subwavelength scales, plasmonic nanostructures offer a unique playground to develop novel concepts for light manipulation at the nanoscale. The recent advances in the field of nonlinear plasmonics include (but are not limited to) the study of subwavelength solitons in metal–dielectric multilayer structures [5] and arrays of metal nanowires [6], the discovery of plasmonic kinks and oscillons in arrays of metal nanoparticles [7,8] and a detailed analysis of their properties [9,10], as well as the analysis of all-optical switching [11–13], generation of terahertz radiation from silver nanoscale dimers [14,15], and a dimer nanoradar with periodically rotating scattering pattern [16,17].

In this paper, we consider both one-dimensional (1D) arrays and two-dimensional (2D) lattices composed of metal nanoparticles and study theoretically a number of novel nonlinear effects in these subwavelength plasmonic systems assuming the nonlinear Kerr-like response of metal nanoparticles and the presence of an external driving field. We demonstrate generation of stationary, drifting and breathing kinks (or switching waves) as well as long-lived trapped and walking plasmonic solitons.

## 2. Model and governing equations

Figure 1*a*–*d* shows the geometries of our problems: identical spherical silver nanoparticles arranged into a 1D chain and a 2D square lattice. We assume that these arrays are embedded into an SiO_{2} host medium with a permittivity *ε*_{h}=2.1 being driven by an arbitrary electric field at a frequency close to the frequency of the surface plasmon resonance of an individual particle, *ω*_{0}. We assume that the nanoparticle size is much less than the wavelength, and the particle radii, *a*, and the lattice constant, *d*, satisfy the condition *d*=3*a*, so that we can employ the point-dipole approximation [18]. Taking into account intrinsic metallic nonlinearity, we write a particle's dielectric permittivity as follows:
2.1
where , , [19], *χ*^{(3)} is the cubic susceptibility and *E*^{(in)} is the local field inside a particle (hereinafter we assume time dependence).

In general, several different mechanisms contribute to *χ*^{(3)}, including pondermotive forces [20], electron confinement [21–23] and interband thermal modulation [24,25]. This leads to a strong dependency of both real and imaginary parts of *χ*^{(3)} on the type of metal, particle size as well as external pulse duration and frequency [26]. However, when a silver nanoparticle with *a*≥5 nm is driven at *ω*∼*ω*_{0} by a laser pulse whose energy is significantly less than the silver ablation threshold, electron confinement nonlinearity associated with intra-band electron transitions dominates, resulting in a purely real cubic susceptibility *χ*^{(3)}≃3×10^{−9} esu [22,23].

To study nonlinear dynamics of the nanoparticle arrays, we employ the slow varying amplitude approximation which is fully applicable, because each particle acts as a resonantly excited oscillator with slow (in comparison with the light period) inertial response. In the case of a chain of metal nonlinear nanoparticles, the governing equations for the dimensionless envelopes of transversal (with respect to the chain axis) polarization components of the particles can be written as [7–9]
2.2
where
and are the dimensionless slowly varying amplitudes of the particle's dipole moments and external electric field, respectively, , describes thermal and radiation losses of particles (note the last term can be neglected as long as *k*_{0}*a*≪1), , *Ω*=(*ω*−*ω*_{0})/*ω*_{0} and *τ*=*ω*_{0}*t*. Here, we use the same normalization for particle polarization envelopes and external electric field as in [7–9].

To obtain the governing equations for a 2D lattice, one should make the following replacements in equation (2.2): in indexes and in [10]. We point out that this model takes into account all particle interactions through the dipole fields, and it can be applied to both finite and infinite arrays.

## 3. Stability analysis

An essential prerequisite for the existence of plasmonic kinks and solitons is a stable continuous external field background featuring a bistability area. Thus, we first investigate stability of the homogeneous stationary solution of equation (2.2), when d/*dτ*=0, and . To this end, we assume a periodically modulated perturbation of the form , where ‘*’ means complex conjugation, *K* is the modulation wavenumber, *A*_{⊥} is an amplitude of perturbation. Having substituted into equation (2.2), we derive the expression for the instability growth rate:
where . A transition from to has been made via the replacement |*n*−*n*′|=*j* and taking into account symmetry structure of the series. In the case of a 2D lattice, one should change and .

Obviously, the conditions λ_{⊥}≥0 at *K*=0 and λ_{⊥}≥0 at *K*≠0 define the areas of bistability and modulational instability (MI) of the homogeneous stationary solution of equation (2.2) on the plane (*Ω*, ), as shown in figure 2. Importantly, MI occupies only a small region of intensities close to the upper threshold of the bistability area, making possible kink and soliton generation. On the contrary, the longitudinal polarization of the electric field with respect to the chain axis or lattice plane extends MI over the whole upper branch of the bistability zone [7–9]. Hence, longitudinally polarized kinks and solitons are unavailable for observation at a homogeneous background.

It should be noted that a perfect homogeneous transversal excitation of a chain corresponds to just a normal light incidence, whereas in the case of a 2D lattice it is beyond practical realization owing to the wave nature of light. Nevertheless, the significant recent progress in generation of longitudinally polarized light beams [27] manifests the feasibility of a *quasi-homogeneous* 2D transversal excitation.

Interestingly, MI and bistability zones emerge from the same point for a nanoparticle chain, as shown in figure 2*a*. This bicritical bifurcation point, also known as Lifshitz point, is characterized by coexistence of the first- and second-order phase transitions, and its appearance in a nonlinear plasmonic system shows a universal character of nonlinear dynamics for systems of different nature [28–30].

## 4. Kinks

To study kink dynamics, we perform numerical simulations of equation (2.2) for a finite array at initial conditions corresponding to different branches of the homogeneous stationary solution over different parts of the array. To suppress the impact of edge effects on kink dynamics, we consider an array of 200 nanoparticles for a chain and 301×301 nanoparticles for a lattice. Light intensity is taken to be a step-like function of time. Constant results in a constant kink velocity, *v*. That is why a step-like time dependence of allows us to control *v*. The velocity sign is defined as positive if the motion of the switching wave leads to the expansion of the region occupied by the upper branch of bistability; otherwise, the velocity is marked as negative. Hereinafter, we treat kink and soliton velocities as time-averaged values because the discrete nature of nanoparticle arrays leads to small-amplitude harmonic oscillations in the instantaneous value of *v* [31].

The result for a chain is presented in figure 3*a*–*c*, where the kink velocity is switched from negative to positive values, passing through zero. A rigorous analysis manifests that there is no smooth transition between kinks with negative, positive and zero velocities, as shown in figure 3*c*. Moreover, these states possess different structures: kinks with positive and zero velocities are characterized by the one-particle width; whereas the width of kinks with negative velocity extends for 12 particles. We also note that negative-velocity kinks exist for the same intensities as stationary (zero-velocity) kinks do. Therefore, to provoke a domain wall to move in the opposite direction, one should force the dipole next to the domain wall to switch into the lower branch of the bistability curve until the domain wall starts moving. An opposite transition can be realized through only variation in the intensity. We attribute this sharp difference between negative- and positive-velocity kinks to a long-range dipole–dipole interaction inherited in our model. Importantly, the dependency shown in figure 3*c* is typical for any frequency fixed inside a bistability area.

Contrary to their 1D counterparts, 2D kinks with negative and positive velocities are characterized by the same widths of nine nanoparticles. Moreover, they could have vanishing velocity at one value of the intensity only, which is also referred to as Maxwell point [32]. When *Ω*≥0 kink velocity is a monotonically increasing function of the intensity; whereas at *Ω*<0 kinks with negative and positive velocities form a hysteresis loop, as shown in figure 3*d*,*e*. This feature gives a straightforward recipe to generate plasmonic kinks with a self-pulsing profile: if one excites a lattice by, for instance, a longitudinally polarized (with respect to a wavevector) Gaussian beam with a properly chosen maximal value and a waist [27], then the kink velocity will oscillate between positive and negative branches inside a bistability loop. The period and the radii of such self-oscillations are defined by parameters of the Gaussian beam. For the breathing kink, shown in figure 3*f*, they are Δ*τ*≈2400 (which is equal to 0.5 ps) and nanoparticles. Thus, this principle can be used to create an ultrafast all-optical plasmonic modulator whose spatio-temporal action will be encoded by spatial gradient of the external field.

## 5. Walking solitons

Next, we consider solitons in the 1D plasmonic chain driven by a plane wave at normal incidence. Following the standard Newton iteration scheme, we find a soliton family and simultaneously determine its stability, as shown in figure 4*a*. Interestingly, the branch of stable solitons coexist with zero (or close to zero) velocity kinks (compare figures 3*a* and 4*a*). This underlines a common nature of these nonlinear modes since discrete solitons can be interpreted as two tightly bound kinks with opposite polarities. Additionally, this soliton family can be treated as a consequence of the so-called snaking bifurcation considered earlier in Ginzburg–Landau models [31,33].

We note that a long-range coupling formally gives rise to asymmetric solitons, but all these branches are unstable. Additionally, we reveal the generation of stable edge states (see figure 4*a*). Their amplitude is a monotonically increasing function of , which is restricted by a bistability threshold and some critical value of . At higher intensities, the amplitude of the edge state becomes strong enough to be attracted by the upper bistable branch, leading to spontaneous generation of two kinks at the edges of the chain [7].

Within the homogeneous excitation, solitons always stand at rest regardless their width because the effective periodic potential created by the chain requires a finite value of the applied external force to start soliton motion [32]. In order to study soliton mobility, we excite the chain by the tilted light incidence: , where 0≤*Qd*≤*k*_{0}*d*=0.7. Characteristic results are summarized in figure 4*b*–*d*.

Figure 4*b* shows soliton velocity as a function of *Qd*. Interestingly, the threshold value of *Qd*, when the soliton starts to move, depends on , but this dependency is not monotonic. This can be explained by the fact that contributes to both the chain potential and the force acting on the soliton; whereas *Qd* contributes only to the force. Another interesting feature of the velocity dependence on *Qd* is its non-monotonic character: instead of expected continuous growth (because increasing *Qd* corresponds to the rise of the applied force), we observe a slight decrease when *Q* approaches *k*_{0}. We attribute this phenomenon to MI which extends over all bistability region for *Q*∼*k*_{0} [9] and results in formation of small-amplitude background, which inhibits soliton motion.

We note that solitons do not keep their shape while drifting because of long-range coupling, as shown in figure 4*c*. In addition, the chain edges are opaque for solitons. Finally, figure 4*d* shows MI-induced multiple spontaneous soliton generation for the light propagation along the chain axis. In this case, it is difficult to identify individual soliton's velocity because of strong interaction between solitons. Therefore, the green curve does not reach the value *Q*=*k*_{0}.

Next, we analyse 2D solitons in the lattices. Numerical simulations show that 2D solitons can be divided into the classes of trapped and walking localized modes, as shown in figure 5. We also identify numerically the existence of domains on the parameter plane ‘intensity versus frequency’ corresponding to trapped and walking solutions and show them in figure 2*b*. A family of trapped solitons includes both symmetric and asymmetric solutions with deep subwavelength extension of just a few nanoparticles. Importantly, these sorts of solitons always remain at rest regardless of the inhomogeneity of the applied field.

Walking solitons are characterized by a wider localization with approximately 150 nanoparticles, which does not depend on *Ω*, , and the width of the step defined by the initial conditions. We find that walking solitons bifurcate at *Ω*=0 similar to 2D kinks. When *Ω*≥0, these solitons possess a stable stationary profile, drifting with a comparatively small velocity of approximately 10^{−4}*c*, where *c* is the speed of light, under impact of the edge effects and colliding elastically with the lattice boundaries and between each other. However, at *Ω*<0 the walking solitons transform into oscillons [8] which move much faster with a velocity of approximately 10^{−2}*c*, merging into one oscillon and radiating through the lattice boundaries. Finally, we note that a powerful tool for steering of the soliton motion is the phase or amplitude inhomogeneity in the applied field, because the drifting direction is defined by an inhomogeneity gradient.

## 6. Concluding remarks

In our calculations summarized above, we deal with dimensionless variables without specifying the exact geometrical parameters of the system and the applied external field. However, for the discussions of experimental realizations of the presented ideas, we have to convert the parameters into real values. The use of the dipole approximation for modelling silver nanoparticles within classical electrodynamics implies that the particle radii should be in the range between 5 and 20 nm, and the surface plasmon resonance frequency , so that the soliton localization in nanoparticle arrays may reach 0.05λ.

It is also insightful to estimate the maximal pulse duration, because the saturation intensity of the external field –10^{−5}, corresponding to an intensity of approximately 1–10 MW cm^{−2}, can lead to thermal damage. Using the value of the ablation threshold of 3.96 J cm^{−2} obtained for silver particles in a SiO_{2} host matrix in the picosecond regime of illumination [34] and taking into account the amplification of the electric field inside the Ag nanoparticle due to surface plasmon resonance, we evaluate the maximal pulse duration as 5 ns for the peak intensity of 0.75 MW cm^{−2} (). Since this value is much greater than the time needed for the generation of nonlinear dissipative modes, being limited to several picoseconds, we expect that all predicted nonlinear phenomena can be readily observed in experiments.

Thus, we have studied the transversal nonlinear modes in 1D arrays and 2D lattices of nonlinear metal nanoparticles driven by an external laser beam. We have demonstrated the generation of stationary, drifting and breathing kinks as well as trapped and walking plasmon solitons. We believe that our findings may pave a way towards further experimental studies of nonlinear plasmonic nanostructures, which could have important implications for active nanophotonic devices operating beyond the diffraction limit.

## Funding statement

This work was supported by the Ministry of Education and Science of Russia (project 14.124.13.6087-MK), Government of the Russian Federation (grant no. 074-U01) and the Australian National University.

## 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.