text:
# Bose-Einstein Condensates in Optical Lattices: Band-Gap Structure and Solitons

## I Introduction

Trapping of neutral atoms in optical lattices Salomon et al. (1987) provides an effective and powerful means for controlling the properties of Bose-Einstein condensates (BECs), as confirmed by a number of the recent experiments (see, e.g., Refs. Anderson and Kasevich (1998); Orzel et al. (2001); Burger et al. (2001); Cataliotti et al. (2001); Christiani et al. (2002); Morsch et al. (2002); Denschlag et al. (2002)). Coupling between the BEC droplets localized in weakly interacting minima of a standing light wave creates a reconfigurable matter-wave structure that can be easily manipulated by varying the strength of the periodic potential and its well spacing.

Many properties of such arrays of the BEC droplets can be understood, within the framework of the mean-field approach, by employing the tight-binding approximation borrowed from solid state physics Javanainen (1999); Trombettoni and Smerzi (2001a); Abdullaev et al. (2001). Being based on the assumption of strong localization of the condensate wave functions in the individual potential wells of the lattice, the tight-binding approximation and the resulting discrete lattice models provide a limited description of the BEC dynamics in an optical lattice. An alternative reduction of the continuous system to the discrete vector lattice involves employing a complete set of on-site Wannier states Alfimov et al. (2002a), and is more general than a tight-binding approximation. However, it is to be clarified whether such a reduction can provide predictions on the dynamics and stability of a BEC in a periodic potential more readily than the full mean-field model. A more adequate analysis of both periodic and localized modes of the continuous Gross-Pitaevskii (GP) equation with a periodic lattice potential is still missing. Some of the recent studies of the GP model Zobay et al. (1999); Bronski et al. (2001); Konotop and Salerno (2002); Alfimov et al. (2002b) point towards interesting features of the BEC dynamics that should be observed beyond the applicability of discrete models.

The main purpose of this paper is to analyze the structure and stability of the stationary states of a condensate in an optical lattice, the so-called nonlinear modes Kivshar et al. (2001), in the framework of a continuous mean-field model. This approach allows us to consider more general types of stationary states and BEC dynamics without many of the restrictions of the discrete models. We solve the continuous one-dimensional GP equation numerically and describe families of stationary states, which possess different spatial structure and stability properties.

First, we find spatially extended stationary nonlinear periodic modes of the condensate in a lattice. These periodic modes possess different symmetry, and provide a generalization of both the Bloch modes in the theory of linear periodic systems and of the extended states of a condensate in a finite optical lattice with a confining potential Steel and Zhang (1998). We demonstrate that the matter-wave spectrum for such periodic BEC states has a band-gap structure, which survives in the presence of the mean-field nonlinearity in the GP equation.

We also find novel types of condensate states which cannot exist in the linear limit. These are spatially localized modes of the condensate which exist in the band gaps of the matter-wave spectrum and, therefore, they should be called gap solitons. The existence of the matter-wave gap solitons was first suggested in the framework of the coupled-mode theory Zobay et al. (1999), in analogy to optical gap solitons in Bragg gratings. Remarkably, gap solitons - spatially localized atomic wavepackets - may exist even in a condensate with repulsive interatomic interactions Zobay et al. (1999). Their experimental demonstration, however, remains a challenge. Here, we present a comprehensive theory of the structure and linear stability analysis of different families of gap solitons with respect to variations in the lattice parameters and effective nonlinearity defined by the number of atoms in the condensate. We discuss the key scenarios of instability of spatially localized modes in optical lattices and the physical mechanism for creating different multiple states of the BEC droplets.

We would like to emphasize that many of the effects described here in the framework of the GP equation with a periodic potential, such as the mode stability and the existence of the localized modes in different band gaps, cannot be obtained in the framework of the tight-binding approximation since the coherent and spectral properties of these two models are different. This latter issue is similar to the BEC dynamics in a double-well potential where a rigorous coupled-mode wave theory should be based on the analysis of nonlinear modes of the entire structure, rather than condensate wave functions of isolated wells Ostrovskaya et al. (2000).

The paper is organized as follows. In Sec. II we present our model and reduce it to the one-dimensional GP equation, under the assumption of a cigar-shaped condensate. Section III includes the analysis of both linear and nonlinear matter-wave spectra in an optical lattice. We revisit the structure of the band-gap spectrum for the Bloch-like modes of the noninteracting condensate, and extend this concept to the nonlinear case, demonstrating the existence of the nonlinear Bloch-like modes and the nonlinearity-dependent band-gap structure of the matter-wave spectrum in optical lattices. The study of the BEC gap solitons in different band gaps is summarized in Sec. IV, where we also present our results on the stability of gap solitons, for both repulsive and attractive atomic interactions. Section V concludes with a summary and further perspectives.

## II Model

We consider the macroscopic dynamics of a BEC loaded onto an elongated optical lattice, similar to those created in a number of experiments (see, e.g., Refs. Salomon et al. (1987); Anderson and Kasevich (1998); Orzel et al. (2001)). The condensate dynamics is described by the three-dimensional GP equation,

$$i\mathrm{}\frac{\mathrm{\Psi }}{T}=\left[\frac{\mathrm{}^2}{2m}^2+V(X,Y,Z)+\mathrm{\Gamma }|\mathrm{\Psi }|^2\right]\mathrm{\Psi },$$

where $`\mathrm{\Psi }(X,Y,Z;T)`$ is the wave function of the cigar-shaped condensate, $`(Y,Z)`$ are the directions of strong transverse confinement, and $`X`$ is the direction of the lattice. The combined potential of the optical lattice and magnetic trap $`V(X,Y,Z)`$ can be written as, $`V(X,R^2)=E_0\mathrm{sin}^2(\pi X/d)+(1/2)m(\omega _{x}^{}{}_{}{}^{2}X^2+\omega _{}^2R^2)`$, where $`R^2=Y^2+Z^2`$, $`E_0`$ is the well depth of the optical lattice, $`d`$ is the characteristic lattice constant, and $`\omega _i`$ are trapping frequencies in the corresponding directions. The parameter $`\mathrm{\Gamma }=4\pi \mathrm{}^2a_s/m`$ characterizes the s-wave scattering of atoms in the condensate which introduces an effective nonlinearity in the mean-field equation; it is positive for repulsive interactions and negative for attractive interactions.

Measuring time in units of $`\omega _{}^1`$, the spatial variables in units of the transverse harmonic oscillator length, $`a_0=(\mathrm{}/m\omega _{})^{1/2}`$, the wave function amplitude in units of $`a_0^{3/2}`$, and the potential in units of $`\mathrm{}\omega _{}`$, we obtain the following dimensionless GP equation:

$$i\frac{\mathrm{\Psi }}{t}=\left[\frac{1}{2}^2+V(x,r^2)+g_{3\mathrm{D}}|\mathrm{\Psi }|^2\right]\mathrm{\Psi },$$
(1)
where $`t=T\omega _{}`$, $`(x,r)=(X,R)/a_0`$, and $`g_{3\mathrm{D}}=4\pi (a_s/a_0)`$ . The potential now takes the form $`V(x,r^2)=V_0\mathrm{sin}^2(Kx)+(1/2)(\mathrm{\Omega }^2x^2+r^2)`$,where $`V_0=E_0/(\mathrm{}\omega _{})`$, and $`K=\pi a_0/d`$. The ratio of the confinement strengths for the magnetic trap $`\mathrm{\Omega }=\omega _x/\omega _{}`$ varies from $`10^1`$ Greiner et al. (2001); Burger et al. (2001); Cataliotti et al. (2001); Orzel et al. (2001) to $`1/\sqrt{2}`$ Christiani et al. (2002); Morsch et al. (2002); Denschlag et al. (2002).

To simplify our analysis, we consider a cigar-shaped condensate in a strongly elongated trap ($`\mathrm{\Omega }10^1`$) and a quasi-one-dimensional optical lattice in the direction of a weak confinement. Under these assumptions, the model can be reduced to a one-dimensional GP equation by assuming separable solutions of the form, $`\mathrm{\Psi }(r,x;t)=\mathrm{\Phi }(r)\psi (x;t)`$. Due to the tight confinement of the condensate in the transverse direction, we assume that the spatial structure of the function $`\mathrm{\Phi }(r)`$ is well described by a solution of the two-dimensional radially symmetric quantum harmonic oscillator problem, $`_{}^2\mathrm{\Phi }+\mathrm{\Phi }r^2\mathrm{\Phi }=0`$. Taking the ground state solution, $`\mathrm{\Phi }(r)=C\mathrm{exp}(r^2/2)`$, where $`C^2=1/\pi `$ is found from the normalization condition, $`_{\mathrm{}}^{\mathrm{}}|\mathrm{\Phi }|^2𝑑y𝑑z=1`$, and integrating Eq. (1) over the transverse coordinate, we derive an effective one-dimensional GP equation in the dimensionless variables,

$$i\frac{\psi }{t}=\frac{1}{2}\frac{^2\psi }{x^2}+V(x)\psi +\sigma g_{1\mathrm{D}}|\psi |^2\psi ,$$
(2)
where $`g_{1\mathrm{D}}=2|a_s|/a_0`$, and the coefficient $`\sigma =\text{sgn}(a_s)=\pm 1`$ characterizes the type of the s-wave interactions. Due to the weak magnetic confinement in $`x`$ direction, it is possible to omit the term $`\mathrm{\Omega }^2`$ from the effective potential, thus assuming that the 1D lattice potential , $`V(x)=V_0\mathrm{sin}^2(Kx)`$, is uniform and is not affected by the presence of the additional confining potential. Furthermore, it is convenient to reduce the parameter space of the model i.e. effectively set $`g_{1\mathrm{D}}=1`$) by rescaling the wave function as $`\psi \psi \sqrt{g_{1\mathrm{D}}}`$. In physical terms, it means that the peak density of the condensate is used to characterize the effective mean-field nonlinearity.

The model (2) can describe the BEC in a one-dimensional optical lattice over a wide range of experimental parameters, such as lattice depth, well spacing, and the effective nonlinearity. Indeed, Eq. (2) is invariant with respect to the following scaling transformation: $`\{K,V_0,x,t,\psi \}\{K/\eta ,V_0/\eta ^2,x\eta ,t\eta ^2,\psi /\eta \}`$, where $`\eta `$ is a dimensionless constant. However, a significant point of difference between the lattice in our model and those in experiments Salomon et al. (1987); Anderson and Kasevich (1998); Orzel et al. (2001) is the absence of an additional confining potential (e.g., gravity or magnetic trapping).

We determine the experimentally relevant range of dimensionless parameters for our model by taking the case of $`{}_{}{}^{87}\text{Rb}`$ atoms, which sets $`m=1.44\times 10^{25}`$ $`\mathrm{kg}`$ and $`a_s=5.3`$ $`\mathrm{nm}`$. The characteristic lattice constant, $`d`$, is determined by the angle between the intersecting laser beams forming the lattice and, in the current experiments, can be varied in the range $`0.41.6`$ $`\mu \mathrm{m}`$ Christiani et al. (2002); Morsch et al. (2002). The lattice depth, $`E_0`$, scales linearly with the light intensity, and varies between zero and $`E_0^{\mathrm{max}}20E_{\mathrm{rec}}`$, where $`E_{\mathrm{rec}}=\mathrm{}^2\pi ^2/(2md^2)`$ is the lattice recoil energy Christiani et al. (2002); Morsch et al. (2002). Taking the transverse frequency as $`\omega _{}2\pi \times 10^2`$ $`\mathrm{Hz}`$ Orzel et al. (2001); Burger et al. (2001), we obtain the following regions of the model parameters: $`2K8`$ and $`4.5\times 10^1V_0^{\mathrm{max}}7.0\times 10^2`$, respectively. In what follows, for the sake of computational convenience, we will make use of the scaling properties of the model and use two different sets of the parameters: $`K=0.1,V_00.1`$ ($`\eta =5`$), and $`K=0.4,V_01.0`$ ($`\eta =20`$).

We note that the number of atoms in the condensate wavefunction, $`𝒩`$, is now given by the expression $`𝒩=N/g_{1\mathrm{D}}`$, where

$$N=_{\mathrm{}}^{\mathrm{}}|\psi (x,t)|^2𝑑x$$
(3)
is an integral of motion for the model (2).

## III Band-Gap Spectrum of Matter Waves

### III.1 Linear Bloch modes

Stationary states of the condensate in an infinite periodic potential of an optical lattice (‘nonlinear modes’) are described by solutions of Eq. (2) of the form: $`\psi (x,t)=\varphi (x)\mathrm{exp}(i\mu t)`$, where $`\mu `$ is the corresponding chemical potential. The steady-state wave function $`\varphi (x)`$ obeys the time-independent GP equation:

$$\frac{1}{2}\frac{d^2\varphi }{dx^2}\left[V_0\mathrm{sin}^2(Kx)\varphi \mu \varphi \right]\sigma |\varphi |^2\varphi =0,$$
(4)
where the lattice depth $`V_0`$ and period $`K`$ are the experimentally controlled parameters, and they both can be varied in our model.

In the case of noninteracting condensate, which formally corresponds to $`\sigma =0`$, Eq. (4) is linear in $`\varphi `$ and, as has been shown in Berg-Sørensen and Mølmer (1998), the Bloch-function formalism applies. In this case (4) can be written as the Mathieu’s equation:

$$\frac{d^2\varphi }{dy^2}[2q\mathrm{cos}(2y)p]\varphi =0,$$
(5)
where $`y=Kx`$, $`q=V_0/(2K^2)`$, and $`p=2(q+\mu /K^2)`$. By employing the Floquet-Bloch theorem, the condensate wave function can be presented as a superposition of Bloch waves,

$$\varphi (y)=b_1P_1(y)e^{i\lambda y}+b_2P_2(y)e^{i\lambda y},$$
(6)
where $`P_{1,2}(y)`$ are periodic functions with the least period $`\pi `$, $`b_{1,2}`$ are constants, and $`\lambda `$ is the characteristic of the Floquet exponent. The structure of all solutions to the Mathieu’s equation can be determined by Floquet theory (this text-book analysis can be found, e.g., in Ref. Jordan and Smith (1977)). According to this theory, the spectrum of Eq. (5) consists of domains (‘bands’) of the $`(p,q)`$ space in which there exist only amplitude-bounded, oscillatory solutions. The bands are separated by the regions (‘gaps’) where unbounded solutions exist. The band edges correspond to exactly periodic solutions.

Figure 1 presents the band-gap diagram for the extended solutions of Eq. (5) which describe noninteracting condensed atoms in an optical lattice. The results are presented for the parameter domain $`(\mu ,V_0)`$ relevant to our problem. The shaded areas correspond to the regions where the Floquet exponent has a real part (which is an odd or even integer) and, therefore, spatially oscillating (but not in general periodic) solutions to the linear model equation (4) exist. If $`\lambda `$ has a nonvanishing imaginary part, the solutions do not oscillate (decay or grow exponentially with $`x`$); such solutions belong to the gaps separating the bands of oscillatory solutions. The first two gaps where unbounded solutions exist, numbered $`\mathrm{I}`$ and $`\mathrm{II}`$ in Fig. 1 are bordered by the band edges representing the linear Bloch waves with periodicity of $`2\pi /K`$ and $`\pi /K`$, respectively.

It is in these gaps of the linear spectrum that the nonlinear localization of matter waves in the form of gap solitons occurs, and it is this localization that could enable the creation of solitons in the BEC with repulsive interatomic interaction. The chemical potentials corresponding to such localized waves lie in the band gaps of the spectrum of the noninteracting condensate (marked areas in Fig. 1).

The dashed line in Fig. 1 marks the value of the lattice depth for which the the study of the gap modes presented below is conducted. We note that this value is greater than the threshold value $`V^{}=\mu `$ that roughly marks the transition between the quasi-unbounded ($`V_0<V^{}`$) and strongly bounded ($`V_0>V^{}`$) condensate wavefunctions (see also discussion in Ref. Drese and Holthaus (1997)). The tight-binding or ”fragmented” condensate regime corresponds to $`V_0V^{}`$, where the bands ”collapse” to discrete levels of the Wannier states Berg-Sørensen and Mølmer (1998).

The knowledge of a complete band-gap spectrum of the matter-waves in an optical lattice provides us with important clues on the existence and stability regions of different types of BEC localized states that may be excited in an optical lattice. The frequently employed nearest-neighbor tight-binding approximation based on a discrete model Trombettoni and Smerzi (2001a, b); Abdullaev et al. (2001) is inferior to the analysis of the complete continuous model in that it describes only one (the first) band surrounded by two semi-infinite gaps. The following sections aim to provide us with understanding of the details of the formation and stability of ”gap modes” - the nonlinear localized states of the condensate in a lattice. Despite the recent advances in the treatment of the continuous GP models Bronski et al. (2001); Wu et al. (2002); Alfimov et al. (2002b), such an analysis is still missing 111At the time of final preparations of this publication a paper containing stability analysis of fundamental matter-wave solitons in the first gap has appeared on Los Alamos preprint archive: K.M. Hilligsøe et al., arXiv:cond-mat/0208136. The results of our analysis are consistent with the findings of this work.

### III.2 Nonlinear Bloch modes

The effect of interatomic interaction, described by the cubic nonlinearity in the mean-field GP equation, results in an effective shift of the chemical potential for the extended periodic states, i.e. Bloch waves of BEC, on the band edges. This leads to the corresponding shifts of the band-gap edges, and the gaps can even close up for large densities of the repulsive condensates. The nonlinearity-induced shift of the bands can be described by the multi-scale perturbation theory developed for the Mathieu’s equation with cubic nonlinearity (4) and a shallow potential ($`V_01`$) (see, e.g. Ref. El-Dib (2001)). Here, we study the extended periodic stationary states of the interacting condensates in the form of nonlinear Bloch waves, by solving Eq. (4) numerically with standard relaxation techniques.

In correspondence to the linear odd and even Bloch waves that satisfy Eq. (5), the nonlinear Bloch waves display similar spatial structures, in which the droplets of the neighboring wells are out of phase or in phase with each other (see Fig. 2). The families corresponding to the two distinct stationary states for different condensate densities are shown in Fig. 2 for the repulsive condensate ($`\sigma =+1`$). The departure of the nonlinear Bloch states from the linear limit (recovered at $`N0`$) indicates the degree to which the mean-field nonlinearity affects the structure of the band gaps in the matter-wave spectrum. However large this departure is, the spatial periodicity and symmetry of the linear Bloch waves corresponding to different bands (see, e.g., Ref. Abramowitz and Stegun (1965)) is preserved in the nonlinear regime. The ”in-phase” and ”out-of-phase” extended periodic states of the condensate corresponding to the lower and upper edges of the lowest band (see Fig. 2, bottom panel) represent a continuous-wave analogue to the so-called ”unstaggered” and ”staggered” modes of the discrete lattice model, respectively Abdullaev et al. (2001). However, none of the higher-order band modes can be described by the discrete model, in the framework of the nearest-neighbor tight-binding approximation, and its validity is therefore restricted to the lowest band only.

## IV Matter-Wave Gap Solitons

### IV.1 Spatially localized modes

We have numerically found different families of spatially localized BEC states in an optical lattice, for both repulsive and attractive interaction. Such states can be regarded as strongly localized soliton-like modes created in band gaps of the optical-lattice spectrum, and therefore they should be compared with gap solitons found for other periodic systems beyond the coupled-mode theory (see, e.g., Refs. Sukhorukov and Kivshar (2001, 2002) and references therein). Several families of the lowest-order gap modes and their spatial structure are presented in Figs. 3 to 7 for $`\sigma =+1`$. The gap solitons exist in all band gaps, excluding the semi-infinite gap of the spectrum below the first band, which is analogous to the total internal reflection gap in optical Bragg gratings.

One of our main results is that each of the band gaps includes a branch of the lower-order fundamental gap modes possessing a single central peak. These fundamental localized states, or the BEC droplets, are strongly localized near a minimum of a potential well, as shown for the first two gaps by the branches A and H in Fig. 3 and the corresponding profiles in Figs. 4, 5, and 6. Higher-order modes can be thought of as bound states of several fundamental modes, including the continuous counterparts of ”twisted modes” of the discrete model Sukhorukov and Kivshar (2002) and multiple states of the localized modes in nonlinear lattices Kivshar et al. (1998).

In the first gap (see Figs. 4 and 5 ), we show only the lowest-order bound states of two (states B, D, F, and G) or three (states C and E) in-phase (states D, E, and F) or out-of-phase (states B, C, and G) fundamental modes. These can be ”even”, i.e. centered between potential wells (states B, C, D, and E) or ”odd”, i.e. centered on the potential minima (states F and G). For low condensate densities, the number of atoms in the localized states that correspond to the bound states of two or three fundamental modes are, correspondingly, two or three times larger than that in a fundamental state (see Fig. 3). The families of the higher-order formed of the same number of fundamental states are degenerate in $`N`$ for low densities; this degeneracy lifts for larger densities, i.e. for larger mean-field nonlinearities. It is interesting to note that, along with the fundamental and its bound states, there exist higher-order localized modes and their corresponding bound states that do not have analogues in the discrete lattice models. A representative mode (K) and its out-of-phase bound state (L) are shown in Fig. 6.

The spatially localized modes exist only in the gaps of the Bloch-wave spectrum. Closer to the band edges, the mode structure changes dramatically. In particular, all the localized modes of BEC with repulsive interaction ( $`\sigma =+1`$) develop extended ”staggered” oscillating tails near the upper gap edges, as demonstrated for the fundamental mode and two of its bound states in the second gap (Fig. 6, right column at $`\mu =0.09`$). The spatial symmetry and periodicity of these tails corresponds to the Bloch-wave structure at the corresponding band edge (see Fig. 6, second column). Inside the bands, these localized states (two examples are shown in Fig. 7) exist as the modes on an extended, oscillating background with the structure of the corresponding Bloch wave \[cf. Figs. 7 and 2(c)\]. These states, localized on a background, are analogous to the ”antidark” optical gap solitons described in Ref. Conti and Trillo (2001), and can, in principle, be observed experimentally.

For the attractive interatomic interaction ($`\sigma =1`$), additionally to the gap modes described above, novel types of localized modes and their bound states are found in the semi-infinite gap, i.e. below the lowest edge of the first band of the matter-wave spectrum (see Figs. 1 and 2). The modes of this type, shown in Fig. 8 are supported by the attractive interaction, and they are similar to conventional bright solitons with exponentially decaying tails that have no nodes and are weakly modulated by the lattice potential. For any fixed value of the lattice amplitude, $`V_0`$, with growing chemical potential, the localized modes approach the lower edge of the first band of the linear spectrum and develop oscillating tails with the structure defined by the structure of the lowest-order Bloch waves \[cf. Fig. 2(a)\].

### IV.2 Stability of localized modes

Stability of extended and localized modes in nonlinear systems is a very important issue, since only dynamically stable modes are likely to be generated and observed in experiments. To determine the stability properties of the localized modes, we consider small perturbations to a solution of the GP equation (2) in the form:

$$\psi (x,t)=\varphi (x)e^{i\mu t}+\epsilon \left[u(x)e^{i\beta t}+w^{}(x)e^{i\beta ^{}t}\right]e^{i\mu t},$$
(7)
where $`\epsilon 1`$, and $`\varphi (x)`$ is the steady state localized nonlinear mode. We linearize Eq. (2) around the localized solution and obtain, to the first order in $`\epsilon `$, the linear eigenvalue problem for the perturbation modes:

$$\left[\begin{array}{cc}\widehat{}_0& \sigma \varphi ^2\\ \sigma \varphi ^2& \widehat{}_0\end{array}\right]\left(\begin{array}{c}u\\ w\end{array}\right)=\beta \left(\begin{array}{c}u\\ w\end{array}\right),$$
(8)
where

$$\widehat{}_0\frac{1}{2}\frac{d^2}{dx^2}V(x)2\sigma |\varphi |^2+\mu .$$

In this presentation, the modes describing the development of instability have either purely imaginary or complex eigenvalues $`\beta `$; in this latter case the instability is called an oscillatory instability. We solve the eigenvalue problem (8) for the condensate localized states found numerically in the previous section by standard finite-difference methods or as a matrix eigenvalue problem with Fourier or polynomial interpolants for the differential operators.

It is well known from the theory of optical gap solitons in nonlinear Bragg gratings, that the linear stability analysis of gap solitons is not trivial and presents significant computational difficulties Schöllmann et al. (1999); Barashenkov and Zemlyanaya (2000); Schöllmann (2000); Schöllmann and Mayer (2000). We have confirmed that the same difficulties arise in the stability analysis of matter-wave gap solitons. The reason for this is the occurence of weak oscillatory instabilities of gap solitons, which are associated with poorly localized perturbation eigenfunctions. Detection of such instabilities for a finite continuous system, whether by finite-differences or pseudo-spectral methods, is highly sensitive to the boundary conditions, the size of the spatial domain, and thus to the periodicity of the system. In addition, the discretization of the continuous system of the spatial extent $`L`$ with $`M`$ grid points has been shown to produce spurious unstable eigenvalues. In general, the convergence of the perturbation eigenvalues as $`\{L,M/L\}\mathrm{}`$ is exponentially slow Schöllmann et al. (1999); Barashenkov and Zemlyanaya (2000). The effective and accurate detection of oscillatory instabilities is possible for $`L\mathrm{}`$, and the infinite domain boundary conditions can be simulated, for example, by expanding the wavefunctions in the basis of Bloch waves (or generalized plane waves in optical case Schöllmann (2000); Schöllmann and Mayer (2000)). The accuracy of this method is limited by the number of Bloch functions in the expansion. Alternatively, the Evans function method used, e.g., in Sukhorukov and Kivshar (2001, 2002) is capable of treating exact boundary conditions at the edges of the spatial domain.

One should bear in mind, however, that the stability results obtained for the infinite domain may not be confirmed by the dynamical simulations of the time-dependent GP equation due to their essentially finite spatial domains. (For instance, optical gap solitons found to be stable on an infinite domain can be unstable on a finite domain Barashenkov and Zemlyanaya (2000)).The problem is not only a numerical one since the physical extent of the optical lattices is always limited by the size of the external trapping potential. Below, we present stability analysis for a large but finite lattice without an external potential. The virtue of this analysis is in its qualitative predictions with regard to the types and scenarios of the instability that is exhibited by the localized modes of the BEC in a 1D lattice. Quantitatively accurate, experimentally relevant calculations of both localized states and their stability should ideally treat a finite lattice with the parabolic trapping potential.

A typical structure of the eigenvalue spectrum is presented in Fig. 9. Since the perturbation eigenfunctions are supported by the periodic potential of the lattice, the full spectrum of the eigenvalue problem (8), i.e. the values of $`\mu \pm \beta `$, has the band-gap structure identical to that presented in Fig. 1. In addition, for any fixed $`\mu `$, there exist modes localized inside the gaps of the continuous spectrum, the so-called ”internal modes”. These linear gap modes exist due to the distributed ”defect” in the lattice potential introduced by the nonlinear localized states. As the chemical potential grows, the internal modes of the localized states move towards the upper edges of the continuum bands (see Fig. 9) and generate an oscillatory instability through the resonance mechanism previously described in a different physical context Sukhorukov and Kivshar (2001). Because different ”internal” perturbation modes with different symmetries can cross into the band and resonate with the linear Bloch waves, the instability can trigger various types of spatial dynamics, including the symmetry breaking instability.

Applying the stability analysis to the modes identified in Sec. IV A, we found that all of the higher-order gap modes identified for our model exhibit oscillatory instability in the regions of the existence domain near the upper band edges. For the parameters used in our calculations the strong exponential instability associated with a purely imaginary $`\beta `$ does not occur has not been detected for gap modes. A typical evolution of the internal perturbation mode in a gap and birth of the oscillatory instability for large $`\mu `$ is shown in Fig. 10 (top), for the even, out-of-phase gap mode. The time scale for development of such weak instabilities is large, which means that oscillatory unstable modes can be dynamically stable on the scales comparable with the typical lifetime of the condensate.

For the lowest-order fundamental gap modes, our analysis has shown the existence of internal modes which, if excited, may be responsible for long-lived oscillations. It also produced very narrow regions of extremely weak oscillatory instability ($`\mathrm{Im}\beta 10^410^5`$) in the close vicinity of the band edges. The accuracy of these calculations is low since the gap mode is poorly localized in this region. As the computed values $`\mathrm{Im}\beta 0`$ with increased accuracy of the method, it is most likely that this oscillatory instability is unphysical and the fundamental gap solitons are linearly stable.

In contrast, some of the higher-order localized modes of the attractive condensate ($`\sigma =1`$), in the semi-infinite gap below the first band, have been found to experience the strong (exponentially growing) linear instability, characterized by a pair of purely imaginary eigenvalues of the perturbation modes $`\beta `$. These modes originate from neutral modes of the excitation spectrum (with $`\mathrm{Re}\beta =0`$). An example of such instability development via bifurcation of the neutral mode into the imaginary plane is presented in Fig. 11 for an odd bound state of fundamental bright solitons. Remarkably, the instability growth rate drops as the mode moves away from the lowest band edge. This is due to the fact that the condensate density grows and, ultimately, strong localization of large-amplitude BEC droplets occurs at the individual wells. In this limit the lattice effects are weak, and the property of a localized state is described by the free-space GP equation, for which the (inverted for our choice of $`\mu `$) Vakhitov-Kolokolov stability criterion states that the fundamental bright solitons of attractive BEC are stable as long as $`N/\mu >0`$, i.e. in the entire semi-infinite gap (see Fig. 8, top).

In order to check the results of the stability analysis of the stationary solutions for the BEC modes, we have employed direct numerical simulations of the time-dependent equation (2), using a split-step Fourier method which preserves the normalization of the condensate wave function with a high accuracy. To confirm the stable dynamics of gap solitons, the results were obtained by running the simulations for 8000 dimensionless time units which is equivalent to a time of 10.6 seconds. Barrett et al. Barrett et al. (2001) who created the first BEC directly in an optical trap, found that for a condensate of $`3.5\times 10^4`$ atoms the $`1/e`$ lifetime was about 3.5 seconds. However, Stamper-Kurn et al. Stamper-Kurn et al. (1998) observed that in a purely optical trap the $`1/e`$ lifetime at low densities was longer than $`10`$ $`\mathrm{s}`$, although at that timescale significant loss occurs as the population goes from $`10^7`$ atoms to $`10^5`$ atoms. The results of our dynamical simulations showing the different instability developments for oscillatory and linearly unstable higher-order states of repulsive and attractive condensates are presented in Figs. 10(b) and 11(b), respectively.

## V Conclusions

We have analyzed the stationary states and stability of the condensate in an optical lattice, in the framework of the GP equation with a periodic potential. We have confirmed some of the earlier predictions for the condensate structure based on the effective tight-binding model, as well as revealed a number of novel features not described by the effective discrete models. In particular, we have extended the concept of the Bloch waves to the case of a nonlinear model and described the families of the nonlinear periodic states which generalize the corresponding Bloch-wave modes of linear periodic systems. We have demonstrated that the spectrum of such periodic modes possesses a band-gap structure, which depends on the number of atoms, i.e. the effective nonlinearity in the mean-field theory.

In the gaps of the matter-wave spectrum, we have numerically found the families of novel spatially localized modes of the condensate in an optical lattice trapped by a few lattice minima. Such modes should be compared with optical gap solitons, and they can exist even for the condensates with repulsive interaction. We have analyzed the structure and, more importantly, the linear stability of the gap solitons which belong to different band gaps, for both repulsive and attractive interatomic interactions. We believe that these novel types of spatially localized, coherent excitations of the condensate can potentially be observed in experiments.

Our studies extend the analysis of the BEC in an optical lattice beyond the applicability limits of the frequently used discrete models based on the tight-binding approximation, as well as they call for the future extensive analysis of the condensate stability in optical lattices described by the three-dimensional GP equation with a periodic potential. In particular, similar to the condensates in a parabolic trap, we expect that optical lattices can support dark solitons, excited on the different types of the nonlinear Bloch waves. The stability of such dark solitons should differ dramatically from the stability of the conventional dark solitons. These studies are now in progress.

## Acknowledgements

We are indebted to N. Robins for help with numerical simulations and to A.A. Sukhorukov for critical reading of this manuscript. We thank C.J. Williams and P.V. Elyutin for useful discussions and suggestions. This work was partially supported by the Australian Research Council and the Australian Partnership for Advanced Computing.


format: text

text:
# Polychromatic solitons in a quadratic medium

## I Introduction

Wave mixing at different carrier frequencies, of which generation of higher harmonics is a well-known example, has had a long history of investigation Franken61 ; Armstrong62 . In an optical medium whose symmetry group lacks a center of inversion, the lowest-order nonlinear response is quadratic, which gives rise to three-photon interactions. The resultant three-wave mixing is resonant when the condition $`\omega _3=\omega _1+\omega _2`$, imposed on frequencies of the three waves, is met. In the particular case when $`\omega _1=\omega _2`$, the process reduces to the second harmonic generation (SHG).

Solitons in quadratic \[$`\chi ^{(2)}`$\] media has been a major area of research recently (see Ref. Etrich01  for a review). The first experimental observation of $`\chi ^{(2)}`$ solitons was reported for type-I SHG, which involves exactly one component of each harmonic, in the $`(2+1)`$-dimensional (bulk) geometry Torruellas95 . The observation of solitons in the $`(1+1)`$-dimensional geometry (planar nonlinear waveguide) followed shortly afterwards Schiek96 . Much theoretical work has been performed for the solitons in both type-I and type-II SHG, the latter case involving two components of the fundamental harmonic, corresponding to different polarizations, and a single component of the second harmonic Stegeman ; Buryak95 ; Tran95 ; Malomed96 ; Buryak97 ; JenaTypeII . The SHG process in an isotropic medium, with two polarizations at both harmonics, has been considered too Boardman .

More complex cases of multi-resonance wave mixing in quadratically nonlinear media have not received much attention because of serious difficulties with their experimental realization using the birefringence-based wave-vector matching schemes. However, the recent rapid development of the quasi-phase-matching (QPM) technique has changed the situation. The technique, originally proposed long ago Armstrong62 , relied upon periodic structures (usually, periodically poled ones) with an alternating sign of the quadratic nonlinearity. Recently, the QPM technique has been extended from periodic to Fibonacci-series-based structures Fradkin99 , and further to fully quasiperiodic ones Fradkin02  (see also Ref. Johansen ). This makes it possible to essentially relax conditions on the wave-vector difference for the matching to take place. Results presented in Ref. Fradkin02  show, both theoretically and experimentally, not only high values of effective $`\chi ^{(2)}`$ coefficients, but also that one can attain the wave-vector matching simultaneously for several sets of waves involved in the nonlinear interactions. In particular, the development of the QPM technique has been an incentive to study multi-resonance systems Kivshar99a ; Kivshar99b ; Towers00 ; Sukhorukov01 , a possible application of which may be design of soliton-based logic gates. Indeed, a soliton may naturally be used as a bit of information, and the interactions of solitons can potentially support logic operations.

In this work, we develop a model of nonlinear mixing between two fundamental-harmonic waves with different frequencies in a quadratic medium. Via the $`\chi ^{(2)}`$ nonlinearity, they generate three different wave components of the second harmonic. Note that interactions between waves with different frequencies in optical media with cubic ($`\chi ^{(3)}`$) nonlinearity is a well-known topic, which has extremely important applications to the WDM multichannel format of data transmission in fiber-optic telecommunication links Agrawal ; however, it appears that the issue has not yet been studied for $`\chi ^{(2)}`$ media, in which the mixing may be realized in both spatial and temporal domains.

We formulate the model in section II, and produce its stationary soliton solutions in section III. These may be both fully polychromatic solitons (PCSs), including all the five field components, and special (simple) solutions of “A” and “B” types, which amount to ordinary two-component SHG solitons in mutually orthogonal (non-intersecting) subsets of the five-wave system. In section IV, we demonstrate, by means of direct simulations, that the solitons of all these types are stable in a broad region in the system’s parameter space. In the same section, we simulate collisions between the simple solitons. When they overlap, the nonlinear interaction generates a component which was absent in both of them before the collision, which makes the final result of the collision strongly inelastic: the former simple solitons develop extra components and become polychromatic ones. Additionally, one or two extra solitons are generated by the collision. We demonstrate that the collision between the initial simple solitons may be a basis for an all-optical XOR logic gate. The paper is concluded by section V.

## II The model

We consider the interaction of five waves in a diffractive dielectric medium with the quadratic nonlinear susceptibility. The carrier frequencies of the waves satisfy resonant conditions, $`\omega _3=2\omega _1,\omega _4=\omega _1+\omega _2,\omega _5=2\omega _2`$, so that $`\omega _1`$ and $`\omega _2`$ may be classified as carriers of the fundamental-harmonic components, while $`\omega _3`$, $`\omega _4`$ and $`\omega _5`$ represent three components of the second-harmonic wave group, generated by the two fundamental-harmonic components via the $`\chi ^{(2)}`$ nonlinearity.

Assuming, as usual, that the wave envelopes $`E_1,E_2`$ and $`E_3,E_4,E_5`$ of these components are slowly varying ones, a system of five equations, coupled parametrically through the components of the $`\chi ^{(2)}`$ nonlinear susceptibility tensor, can be derived from the Maxwell’s equations to govern the evolution of the waves in the spatial domain (the derivation follows the well-known procedure worked out for the usual type-I and type-II $`\chi ^{(2)}`$ interactions, see a detailed account given in the review Etrich01 ):

$`2ik_1{\displaystyle \frac{E_1}{z}}+{\displaystyle \frac{^2E_1}{x^2}}+2\stackrel{~}{\chi }_1E_1^{}E_3e^{i\mathrm{\Delta }k_{3,1,1}}+\sqrt{2}\stackrel{~}{\chi }_2E_2^{}E_4e^{i\mathrm{\Delta }k_{412}}`$ $`=`$ $`0,`$
$`2ik_2{\displaystyle \frac{E_2}{z}}+{\displaystyle \frac{^2E_2}{x^2}}+\sqrt{2}\stackrel{~}{\chi }_2E_1^{}E_4e^{i\mathrm{\Delta }k_{412}}+2\stackrel{~}{\chi }_3E_2^{}E_5e^{i\mathrm{\Delta }k_{522}}`$ $`=`$ $`0,`$
$`2ik_3{\displaystyle \frac{E_3}{z}}+{\displaystyle \frac{^2E_3}{x^2}}+\stackrel{~}{\chi }_1E_1^2e^{i\mathrm{\Delta }k_{311}}`$ $`=`$ $`0,`$ (1)
$`2ik_4{\displaystyle \frac{E_4}{z}}+{\displaystyle \frac{^2E_4}{x^2}}+\sqrt{2}\stackrel{~}{\chi }_2E_1E_2e^{i\mathrm{\Delta }k_{412}}`$ $`=`$ $`0,`$
$`2ik_5{\displaystyle \frac{E_5}{z}}+{\displaystyle \frac{^2E_5}{x^2}}+\stackrel{~}{\chi }_3E_2^2e^{i\mathrm{\Delta }k_{522}}`$ $`=`$ $`0,`$
where $`k_1`$ through $`k_5`$ are the corresponding carrier wave numbers, $`z`$ and $`x`$ are the propagation and transverse coordinates, $`\mathrm{\Delta }k_{lmn}=k_lk_mk_n`$ is the wave-vector mismatch, and $`\stackrel{~}{\chi }_n8\pi \omega _n^2/c^2\chi ^{(2)}`$, $`\chi ^{(2)}`$ being an element of the quadratic susceptibility tensor. The extra factors of $`2`$ in front of the coefficients $`\stackrel{~}{\chi }_1`$ and $`\stackrel{~}{\chi }_3`$ in the first two equations reflects the fact the equations may be derived from a Lagrangian, and the factor $`\sqrt{2}`$ in front of all the terms containing $`\stackrel{~}{\chi }_2`$ is added by definition, to simplify subsequent rescalings.

Equations (1) are normalized by measuring $`x`$ and $`z`$, respectively, in units of the input beam size $`r_0`$ and diffraction length $`z_d=r_0^2k_4`$ at the frequency $`\omega _1+\omega _2`$. Introducing dimensionless fields

$`E_1`$ $`=`$ $`(u_1/\sqrt{2}\gamma )\mathrm{exp}\left(i\beta _1z\right),E_2=(u_2/\sqrt{2}\gamma )\mathrm{exp}\left(i\beta _2z\right),E_3=(u_3/2\gamma )\mathrm{exp}\left(2i\beta _1zi\mathrm{\Delta }k_{311}z\right),`$ (2)
$`E_4`$ $`=`$ $`(u_4/\sqrt{2}\gamma )\mathrm{exp}\left(i(\beta _1+\beta _2)zi\mathrm{\Delta }k_{412}z\right),E_5=(u_5/2\gamma )\mathrm{exp}\left(2i\beta _2zi\mathrm{\Delta }k_{522}\right),`$ (3)
and $`\gamma \stackrel{~}{\chi }_1r_0^2`$, a normalized system of equations is obtained:

$`2{\displaystyle \frac{k_1}{k_4}}\left(i{\displaystyle \frac{u_1}{z}}\beta _1u_1\right)+{\displaystyle \frac{^2u_1}{x^2}}+\chi _1u_1^{}u_3+\chi _2u_2^{}u_4`$ $`=`$ $`0,`$
$`2{\displaystyle \frac{k_2}{k_4}}\left(i{\displaystyle \frac{u_2}{z}}\beta _2u_2\right)+{\displaystyle \frac{^2u_2}{x^2}}+\chi _2u_1^{}u_4+\chi _3u_2^{}u_5`$ $`=`$ $`0,`$
$`2{\displaystyle \frac{k_3}{k_4}}\left(i{\displaystyle \frac{u_3}{z}}\left(2\beta _1\mathrm{\Delta }k_{311}\right)u_3\right)+{\displaystyle \frac{^2u_3}{x^2}}+{\displaystyle \frac{1}{2}}\chi _1u_1^2`$ $`=`$ $`0,`$ (4)
$`2{\displaystyle \frac{k_4}{k_4}}\left(i{\displaystyle \frac{u_4}{z}}\left(\beta _1+\beta _2\mathrm{\Delta }k_{412}\right)u_4\right)+{\displaystyle \frac{^2u_4}{x^2}}+\chi _2u_1u_2`$ $`=`$ $`0,`$
$`2{\displaystyle \frac{k_5}{k_4}}\left(i{\displaystyle \frac{u_5}{z}}\left(2\beta _2\mathrm{\Delta }k_{522}\right)u_5\right)+{\displaystyle \frac{^2u_5}{x^2}}+{\displaystyle \frac{1}{2}}\chi _3u_2^2`$ $`=`$ $`0.`$
where $`\chi _n\stackrel{~}{\chi }_n/\stackrel{~}{\chi }_1`$ ($`n=1,2,3`$, i.e., $`\chi _11`$), $`\beta _{1,2}`$ being two phase-velocity shifts. To reduce the number of parameters in Eqs. (4), the fields and coordinates can be rescaled further. Defining $`u_n=\beta _1U_n`$, $`z=Z/\beta _1`$ , and $`x=X/\sqrt{|\beta _1|}`$, we obtain

$`i{\displaystyle \frac{U_1}{Z}}+{\displaystyle \frac{^2U_1}{X^2}}\alpha _1U_1+\chi _1U_1^{}U_3+\chi _2U_2^{}U_4`$ $`=`$ $`0,`$
$`i{\displaystyle \frac{U_2}{Z}}+{\displaystyle \frac{^2U_2}{X^2}}\alpha _2U_2+\chi _2U_1^{}U_4+\chi _3U_2^{}U_5`$ $`=`$ $`0,`$
$`2i{\displaystyle \frac{U_3}{Z}}+{\displaystyle \frac{^2U_3}{X^2}}\alpha _3U_3+{\displaystyle \frac{1}{2}}\chi _1U_1^2`$ $`=`$ $`0,`$ (5)
$`2i{\displaystyle \frac{U_4}{Z}}+{\displaystyle \frac{^2U_4}{X^2}}\alpha _4U_4+\chi _2U_1U_2`$ $`=`$ $`0,`$
$`2i{\displaystyle \frac{U_5}{Z}}+{\displaystyle \frac{^2U_5}{X^2}}\alpha _5U_5+{\displaystyle \frac{1}{2}}\chi _3U_2^2`$ $`=`$ $`0.`$
where $`\alpha _1=1`$, $`\alpha _2=\beta _2/\beta _1`$, $`\alpha _3=(4\beta _12\mathrm{\Delta }k_{311})/\beta _1`$, $`\alpha _4=(2\beta _1+2\beta _22\mathrm{\Delta }k_{412})/\beta _1`$, and $`\alpha _5=(4\beta _22\mathrm{\Delta }k_{522})/\beta _1`$. We assume that $`\left|k_4k_2k_1\right|k_n`$, and everywhere, except for the phase-mismatch parameters introduced above, one may set $`k_3=k_4=k_5`$, i.e., the three components of the second-harmonic wave have similar wave numbers.

Equations (5) assume that the three parametric interactions (“vertices”), which couple, respectively, the wave triplets (1,1,3), (2,2,5), and (1,2,4), are nearly phase matched (as it was mentioned above, a possible way of achieving this may be provided by the QPM technique based on quasiperiodic structures Fradkin02 ). It is straightforward to see that, like the model describing the type-II SHG Tran95 ; Malomed96 ; Buryak97 , Eqs. (5) have two Manley-Rowe invariants, namely, the total power,

$$Q=_{\mathrm{}}^{\mathrm{}}|U_1|^2+|U_2|^2+4\left(|U_3|^2+|U_4|^2+|U_5|^2\right)dx,$$
(6)
and the power imbalance,

$$Q_{\mathrm{imb}}=_{\mathrm{}}^{\mathrm{}}|U_1|^2|U_2|^2+4\left(|U_3|^2|U_5|^2\right)dx.$$
(7)
The present model does not include walkoff terms (group-velocity mismatch). While a detailed discussion of the walkoff is beyond the scope of this work, it is relevant to mention that QPM and similar techniques, such as tandem structures tandem , make it possible to suppress the walkoff tandem ; Stolzenberger . We also note that QPM can give rise to an effective cubic nonlinearity Clausen , which may compete with the underlying quadratic interactions Bang ; Johansen . For this reason, cubic terms should sometimes be added to a dynamical model, but this is not an issue for immediate consideration in the present context. As concerns the physical realization of the system, estimates using typical values of the relevant physical parameters in such quadratically nonlinear materials as lithium niobate and KTP Fradkin02 ; Johansen  suggest that a necessary (quasi) period of the QPM structure is $`10`$ $`\mu `$m, and the power and transverse size of the soliton beam are expected to be $`20`$ $`\mu `$m and $`100`$ mW, respectively.

If still more resonances are allowed, other essential wave components may appear, for instance, those corresponding to the combinational frequencies $`\omega _{12}2\omega _1\omega _2`$ and $`\omega _{21}2\omega _2\omega _1`$ (obviously, they belong to the fundamental-harmonic wave set). Denoting the corresponding wave numbers as $`k_{12}`$ and $`k_{21}`$, we see that these new components will indeed be essential if the wave-number triplets (12,2,3) and/or (21,1,5) are nearly phase-matched. The accordingly modified system will include seven components (four pertaining to the fundamental harmonic, and three to the second harmonic) and five vertices. However, such a situation seems much more exotic (five simultaneous resonances) than the more generic possibility of three simultaneous resonances underlying the model considered in the present work.

## III Stationary soliton solutions

Particular exact solutions of Eqs. (5) for stationary solitons can be sought for in the form

$$U_n=A_n\mathrm{sech}^2\left(\lambda X\right),n=1,\mathrm{},5,$$
(8)
where $`A_n`$ are amplitudes, and $`\lambda `$ is the inverse width of the soliton. Inserting this (8) into Eqs. (5), a solution can be obtained if $`\alpha _1=\alpha _2=\alpha _3=\alpha _4=\alpha _5\alpha `$ and $`\lambda =\sqrt{\alpha }/2`$. The amplitudes are found to be

$`A_1`$ $`=`$ $`\pm \sqrt{\xi }A_2,A_2=\pm {\displaystyle \frac{6\lambda ^2}{\chi _2^2\xi +\chi _3^2/2}},A_3={\displaystyle \frac{\chi _1A_1^2}{12\lambda ^2}},`$
$`A_4`$ $`=`$ $`{\displaystyle \frac{\chi _2A_1A_2}{6\lambda ^2}},A_5={\displaystyle \frac{\chi _3A_2^2}{6\lambda ^2}},`$ (9)
where $`\xi (\chi _3^2/2\chi _2^2)/(\chi _1^2/2\chi _2^2)`$.

In addition to the exact solution based on Eqs. (8) and (9), other particular solutions can be found setting $`U_{2,4,5}=0`$ or, alternatively, $`U_{1,3,4}=0`$. Equations (5) then reduce to the well-known type-I SHG model Stegeman ; Buryak95 ; Etrich01 ,

$`i{\displaystyle \frac{V}{Z}}+{\displaystyle \frac{^2V}{X^2}}V+V^{}W`$ $`=`$ $`0,`$ (10)
$`2i{\displaystyle \frac{W}{Z}}+{\displaystyle \frac{^2W}{X^2}}\rho W+{\displaystyle \frac{1}{2}}V^2`$ $`=`$ $`0,`$ (11)
where $`V`$ and $`W`$ are the fields at the fundamental and second harmonics. In this case, there is a well-known special exact solution corresponding to $`\rho =1`$ in Eq. (11) Karamzin74 ,

$$V=\left(3/\sqrt{2}\right)\mathrm{sech}^2\left(X/2\right),W=\left(3/2\right)\mathrm{sech}^2\left(X/2\right).$$

When $`\rho 1`$, a family of stationary soliton solutions to Eqs. (10) and (11) can be found numerically Buryak95  (or approximated analytically by means of the variational method Vika ). We will refer to the general solution for the case when $`U_1=V`$, $`U_3=W`$, and $`U_{2,4,5}=0`$ as an “A” type soliton, while the opposite “B” type is defined as the one with $`U_2=V`$, $`U_5=W`$, and $`U_{1,3,4}=0`$. Both these types will also be called “simple” solitons.

General solutions for the polychromatic (five-wave) soliton can be found from the stationary version (the one with $`U_n/Z=0`$) of Eqs. (5) by means of the standard numerical methods for two-point boundary-value problems. In Figs. 1 and 2, comparison is made between the general five-wave solitons and the particular solutions generated by Eqs. (10) and (11) at equal values of all the parameters. In Fig. 1, one can see that the A soliton is much larger in amplitude, while its width is not widely different from that of the polychromatic one. As a consequence of this, the power $`Q`$ \[see Eq. (6)\] of the A-type soliton shown in Fig. 1 is $`Q=36`$, while for PCS in the same figure, $`Q=8`$ \[the other invariant is $`Q_{\mathrm{imb}}=0`$ in the case considered, see Eq. (7)\].

This drastic difference in the powers can be understood: in the case of the full PCS, one has two nonlinear terms in the first two equations (5), rather than one term in the case of the simple solitons; therefore, the amplitude necessary to compensate the spreading out of the beam due to the diffraction term is, roughly, twice as small in comparison with the simple solitons, or, eventually, the power is $`4`$ times as small. To check whether the power of the PCS is indeed essentially lower than that of simple solitons in the general case, in Fig. 2 we show a typical example of the change of the two powers with the variation of $`\alpha _2`$ or, equivalently, $`\rho `$ in Eq. (11). As is seen, PCS indeed persistently maintains a lower power $`Q`$ than its A-type counterpart. On varying the different parameters $`\alpha _n,\chi _n`$ the value of $`Q`$ of the polychromatic soliton may be increased (or decreased) in value from that portrayed in Fig. 2 but it only exceeds that of the simple solitons when $`\chi _2<0.7`$. In the limit of $`\chi _20`$ Eqs. (5) decouple and the polychromatic soliton tends to a “double” simple soliton i.e. both “A” and “B” type solitons existing in one envelope. The double simple soliton by definition has twice the power of single simple soliton. Provided the nonlinearity coefficient $`\chi _2`$ is large enough ($`\chi _20.7`$) it is possible to conclude that PCSs in a quadratically nonlinear medium may be produced from waves with different frequencies at a much lower net input power than the ordinary SHG solitons, i.e., it may be essentially easier to generate PCSs in the experiment and use them in potential applications. Of course, these results are meaningful provided that these solitons are stable.

## IV Stability and interactions of the solitons

### IV.1 The stability of the polychromatic solitons

To test the stability of PCSs, we solved the full system of Eqs. (5), employing the numerical split-step fast-Fourier-transform method. The simulations were performed with a computational grid of 2048 points, the transverse and propagation step-sizes being, respectively, $`h_x0.02`$ and $`h_z=0.01`$. The integration domain had the transverse size $`100`$ ($`50<x<+50`$), which is by far larger than any $`x`$ size relevant to the solitons, see Figs. 4, 5, and 9 below. Absorbing layers were placed at the edges of the computational domain to prevent reentering of radiation. Specially monitoring interaction of the radiation waves with the absorbers, we have verified that, in all the cases considered, no reflection took place indeed.

Evolution of initial configurations close to the stationary solutions was simulated for a variety of parameters. To impose initial perturbations, values of the initial amplitude and width of the pulses were altered against the exact stationary solutions. From the results of the simulations, we have concluded that PCSs survive, clearly remaining stable, as long as the simulations were run, the maximum simulation length being 300 diffraction lengths of the soliton.

To further test robustness of PCSs, we ran numerical experiments in which the solitons were successfully generated from initial Gaussian pulses, launched in the pump fields $`U_1`$ and $`U_2`$ (with no initial second-harmonic components, which corresponds to the generation of SHG solitons in real experiments Stegeman\_experiment ). The results demonstrate that PCSs not only are stable against small perturbations, but also play the role of strong attractors in the system. An example of PCS generation from the Gaussian beams is illustrated by Fig. 3. Typically, there is a period of strong relaxation of the beams, where the amplitudes fluctuate and excess energy is radiated away. The second-harmonic components $`U_{3,4,5}`$ are generated, and the soliton arranges itself to a quasi-stationary form, which then propagates, in a stable fashion, over a distance in excess of 100 diffraction lengths (as long as the simulations were run). As a result of many runs of systematic simulations, PCSs have been found to be attractors in a broad range in the system’s parameter space. We concentrated mainly on the parameter space $`\alpha _n=\mathrm{1..5}`$ and $`\chi _n=\mathrm{1..2}`$ where the PCS definitely had a lower value $`Q`$ than the simple solitons. In connection to this, it is relevant to note that the above-mentioned families of the simple A and B solitons are also stable in the ordinary SHG model, except for a small region of their existence domain Buryak95 ; Etrich01 .

### IV.2 Collisions between orthogonal simple solitons

An interesting possibility is to consider collisions of the mutually orthogonal simple solitons of the above-mentioned A and B types. The overlapping between the colliding solitons will give rise to the generation of the field $`U_4`$, which is absent in both A and B solitons in their pure form, and the issue is how the generation of this field will affect the interaction between the solitons. Equations (11) have the property of the Galilean invariance, so “moving” solitons (in fact, the solitary beams propagating at an angle in the planar waveguide) can be constructed by the transformation

$`V(X,Z)`$ $`=`$ $`V_0\left(XCZ\right)\mathrm{exp}\left(iCX/2iC^2Z/4\right),`$ (12)
$`W(X,Z)`$ $`=`$ $`W_0\left(XCZ\right)\mathrm{exp}\left(iCXiC^2Z/2\right),`$
where $`C`$ is the “velocity” (actually, a slope) of the “moving” soliton.

We collided the A and B solitons with an initial separation between their centers $`x_0=20`$, varying their velocities $`\pm C`$. A representative set of values of the parameters for which the results are displayed below are $`\alpha _{1,\mathrm{},5}=1`$, $`\chi _1=\chi _3=1`$, and $`\chi _2=2`$; many simulations run for other values gave quite similar results.

In the case of moderate collision velocities (and zero phase difference between the solitons), the generation of the field $`U_4`$ in the course of the collision gives rise to a third polychromatic soliton with the zero velocity. Trajectories of the initial A and B solitons alter, and they appear from the collisions as PCSs too, i.e., the collision adds to them components which were initially absent. The three post-collision solitons have roughly equal powers, with a mirror symmetry in the power distribution amongst individual components, see Fig. 4. At higher collision velocities (again, for the zero phase difference), four PCSs are generated by the collision, all having a nonzero velocity (see Fig. 5). As the collision velocity increases, less and less interaction takes place, until the initial solitons pass through one another unchanged (elastically). In all these cases, 90$`\%`$ of the initial power is typically converted into the resulting PCSs, i.e., radiation loss is not a dominant actor in the collision dynamics.

Results produced by many runs of the simulations for the collision of the in-phase (zero-phase-difference) solitons are summarized in Figs. 6 and 7. The constant initial separation $`x_0`$ between the solitons in the simulations means that the increasing velocity corresponds to an increasing incidence angle $`\theta _\mathrm{i}`$. As may be naturally expected, the collision alters the trajectories of the solitons most in the case of low velocities. At higher velocities, the outer solitons keep essentially the same velocity after the collision as they had before it, while the two additionally generated inner solitons are significantly slower.

The most significant inelastic effects occur when the solitons collide with the velocities $`\pm 0.4`$. In this case, the exit trajectory is altered the most (see Fig. 6), and, as per Fig. 7, the largest part of the net power is transferred into the newly generated harmonic components of the outer solitons.

The interaction between the A and B solitons produces nontrivial results also in the case when they have zero initial velocities, but their tails overlap at the initial position. For the initial separation $`x_0=5`$, the result of the interaction with $`C=0`$ is displayed in Fig. 8, in the form of the distribution of the fields $`U_1`$ and $`U_2`$ produced by the interaction. Again, three PCSs appear, the central one with the zero velocity, and two outer solitons with nonzero velocities. In the course of the evolution, a phase difference develops across the component fields. In fact, it is this phase difference which repulses the outer solitons, lending them a nonzero velocity.

It is well known from theoretical Stegeman\_theory  and experimental Stegeman\_experiment  studies of collisions between SHG solitons of the usual type that the result strongly depends on their relative phase at the collision point (provided that the colliding solitons have nearly identical amplitudes): they attract each other and therefore interact strongly in the case of the zero phase difference, and repel each other if the phase difference between the fundamental-harmonic components is $`\pi `$. However, in the present system the effect of the phase difference on collisions between the simple solitons of the orthogonal types is much weaker. Simulations performed with various phase differences (including $`\pi /2`$) between the two orthogonal simple solitons does not display any change in the post-collision dynamics of the beams when compared to the zero phase difference case. Also it may be seen from the structure of the underlying system (5) that phase differences between the “A” and “B” solitons will have little effect. Indeed, these equations are exactly invariant against the transformation $`u_1u_1\mathrm{exp}\left(i\varphi _1\right)`$, $`u_3u_3\mathrm{exp}\left(2i\varphi _1\right)`$, $`u_2u_2\mathrm{exp}\left(i\varphi _2\right)`$, $`u_5u_5\mathrm{exp}\left(2i\varphi _2\right)`$ and $`u_4u_4\mathrm{exp}\left(i\varphi _1+i\varphi _2\right)`$ with an arbitrary phase shift $`\varphi _{1,2}`$.

The interactions between the simple solitons in this system may find application as a basis for an all-optical logic gate. Indeed, consider two solitons of the “A” and/or “B” types, propagating in such a way that they will collide. If the solitons belong to the same type (i.e., the configuration is AA or BB), then they attract or repel one another, depending on their relative phase Stegeman\_experiment ; Stegeman\_theory , or, if the relative velocity is high enough, they simply pass through one another. But if they belong to the opposite types (an AB configuration), at least three PCSs are produced by the collision even at high velocities. This is a behavior which is expected from an exclusive OR gate, alias XOR. The actual outcome can be easily established by checking the $`U_4`$ content in the output (by means of an appropriate optical filter). An advantage of such a design of the XOR gate is that any output beam will be a stationary soliton (even in the case of a strongly inelastic collision), which makes it convenient for further manipulations (cascadability).

## V Conclusion

We have introduced the simplest five-component model of polychromatic solitons (PCSs) in a quadratically nonlinear optical medium. The existence and stability of both polychromatic and two types of simple solitons have been demonstrated. An essential peculiarity of PCSs is that their power is much smaller, at the same values of the control parameters, than the power of the usual two-component (simple) solitons. We have also performed systematic simulations of collisions between mutually orthogonal simple solitons, concluding that the collisions are strongly inelastic (including the interaction between two solitons with the zero initial velocity), giving rise to transformation of the simple solitons into polychromatic ones, and generation of one or two additional PCSs. A value of the relative velocity at which the inelastic effects are strongest has been found. We have also shown that the collision may serve as a basis to design a simple all-optical logic gate of the XOR type.

###### Acknowledgements.

We appreciate valuable discussions with A.V. Buryak. This work was supported, in a part, by the US-Israel Binational Science Foundation under the grant No. 1999459, and by a matching grant from the Tel Aviv University.


text:
# Thermal conductivity of anisotropic and frustrated spin-1/2 chains

(August 14, 2002; revised: September 18, 2002)

## Abstract

We analyze the thermal conductivity of anisotropic and frustrated spin-$`1/2`$ chains using analytical and numerical techniques. This includes mean-field theory based on the Jordan-Wigner transformation, bosonization, and exact diagonalization of systems with $`N18`$ sites. We present results for the temperature dependence of the zero-frequency weight of the conductivity for several values of the anisotropy $`\mathrm{\Delta }`$. In the gapless regime, we show that the mean-field theory compares well to known results and that the low-temperature limit is correctly described by bosonization. In the antiferromagnetic and ferromagnetic gapped regime, we analyze the temperature dependence of the thermal conductivity numerically. The convergence of the finite-size data is remarkably good in the ferromagnetic case. Finally, we apply our numerical method and mean-field theory to the frustrated chain where we find a good agreement of these two approaches on finite systems. Our numerical data do not yield evidence for a diverging thermal conductivity in the thermodynamic limit in case of the antiferromagnetic gapped regime of the frustrated chain.

Introduction -  Transport properties of low-dimensional spin systems have attracted recently interest both from the experimental and theoretical side. A particular motivation comes from the observation that magnetic excitations of one-dimensional spin systems significantly contribute to the thermal conductivity which is manifest in many experiments on materials such as the spin-ladder systemsolo00 ; hess01 ; kudo01  $`\text{(Sr,La,Ca)}_{14}\text{Cu}_{24}\text{O}_{41}`$ and the spin-chain compounds $`\text{SrCu}\text{O}_2`$ and $`\text{Sr}_2\text{Cu}\text{O}_3`$solo01 . Assuming elementary excitations to carry the thermal current and using a relaxation time ansatz for their kinetic equation one finds extremely large mean-free paths being, for example, of the order of 1000Å in $`\text{La}_5\text{Ca}_9\text{Cu}_{24}\text{O}_{41}`$hess01 . Although the magnitude of the mean-free path is currently an issue of intense discussion, the question arises whether heat transport in low-dimensional spin systems is ballistic, i.e., whether intrinsic scattering of magnetic excitations is ineffective to render the thermal conductivity finite. From the theoretical point of view this issue is related to the value of the so-called (thermal) Drude weightzotos97  $`D_{\mathrm{th}}`$ which is the zero-frequency weight of the thermal conductivity $`\kappa `$. A nonzero value of $`D_{\mathrm{th}}`$ corresponds to a diverging thermal conductivity. This scenario is trivially realized if the energy-current operator is a conserved quantity, which is the case for the spin-$`1/2`$ Heisenberg chainniemeyer71 ; zotos97 . For a number of other models like the frustrated chain, the dimerized chains or the spin ladder the energy-current operator is not conserved and the question of nonzero $`D_{\mathrm{th}}`$ is a challenging topic.
In this paper, we establish various numerical and analytical techniques to analyze the thermal Drude weight and to compute the temperature dependence of $`D_{\mathrm{th}}(T)`$. We study the model Hamiltonian $`H=_lh_l`$ with the local energy-density given by

$$h_l=J\{(S_l^+S_{l+1}^{}+\text{H.c.})/2+\mathrm{\Delta }S_l^zS_{l+1}^z+\alpha \stackrel{}{S}_l\stackrel{}{S}_{l+2}\}.$$
(1)
The $`XXZ`$ model ($`\alpha =0`$) is integrable whereas for nonzero frustration the model becomes nonintegrable. Recently, Klümper and Sakai kluemper02  obtained $`D_{\mathrm{th}}(T)`$ for $`\alpha =0`$ and $`0\mathrm{\Delta }1`$ by using the Bethe ansatzkluemper  which allows us to test our approaches in this regime. Exact diagonalization of finite systems up to $`N=14`$ sites has been applied by Alvarez and Gros gros02  to investigate $`D_{\mathrm{th}}(T)`$ for the isotropic Heisenberg chain (i.e., $`\mathrm{\Delta }=1`$), the frustrated chain, and the spin ladder. Our numerical analysis goes beyond this by allowing $`\mathrm{\Delta }1`$ and extension to larger systems with $`N18`$ sites.
Thermal conductivity -  The thermal conductivity $`\kappa `$ is defined by $`j=\kappa T`$ and is given by the following expressionmahan :

$$\kappa (\omega )=\beta _0^{\mathrm{}}𝑑te^{i\omega t}_0^\beta 𝑑\tau j(ti\tau )j.$$
(2)
$`j`$ is the energy-current operatorcurrent  and $`\beta =1/T`$ is the inverse temperature. The current operator satisfies the equation of continuity: $`_th_l=i[H,h_l]=(j_{l+1}j_l)`$. For exchange interaction of arbitrary range, i.e., $`[h_{l\pm m},h_l]0`$ for $`mm_0`$, this implies

$$j_l=i\underset{m,n=0}{\overset{m_01}{}}[h_{lm1},h_{l+n}].$$
(3)
In our case we have $`m_0=2`$ (see Eq. (1)) leading to

$$j=\underset{l}{}j_l=i\underset{l}{}[h_{l2}+h_{l1},h_l+h_{l+1}].$$
(4)
Note that the current operator derived from Eq. (4) includes the proper limiting form for $`\alpha 1`$ where one recovers the current operators of two decoupled chains.
The conductivity $`\kappa (\omega )`$ may be decomposed according to $`\mathrm{Re}\kappa (\omega )=D_{\mathrm{th}}(T)\delta (\omega )+\kappa _{\mathrm{reg}}(\omega )`$ into a singular part at zero frequency and a regular part $`\kappa _{\mathrm{reg}}(\omega )`$. The quantity of interest is the thermal Drude weight $`D_{\mathrm{th}}(T)`$ which can be computed via (see, e.g., Ref. zotos97, )

$$D_{\mathrm{th}}(T)=\frac{\pi \beta ^2}{ZN}\underset{\genfrac{}{}{0pt}{}{m,n}{E_m=E_n}}{}e^{\beta E_m}|m|j|n|^2.$$
(5)
$`Z`$ is the partition function and $`N`$ the number of lattice sites. Note that we exclusively use periodic boundary conditions. If $`j`$ is a conserved quantity, the conductivity reduces to $`\kappa (\omega )=D_{\mathrm{th}}(T)\delta (\omega )`$ and expression (5) simplifies to $`D_{\mathrm{th}}=\pi \beta ^2j^2/N`$.

Mean-field theory -  Using a Jordan-Wigner transformationmahan  the spin operators are mapped to spinless fermionic operators $`c_l^{()}`$. In the case of $`\mathrm{\Delta }=0`$ and zero frustration $`\alpha `$ the corresponding Hamiltonian is diagonal in momentum space and reads $`H=_kϵ_kc_k^{}c_k`$ with a tight-binding dispersion $`ϵ_k=J\mathrm{cos}(k)`$. A nonzero value of $`\mathrm{\Delta }`$ or $`\alpha `$ leads to a four-fermion interaction term that can be treated approximately by Hartree-Fock (for details, see, e.g., Ref. brenig97, ), resulting in a renormalization of $`ϵ_k`$ to $`\stackrel{~}{ϵ}_k=J(1+2A(\mathrm{\Delta }2\alpha ))\mathrm{cos}(k)`$. The parameter $`A=\frac{1}{\pi }_0^\pi 𝑑k\mathrm{cos}(k)f(\stackrel{~}{ϵ}_k)`$ has to be determined self-consistently where $`f(ϵ)=1/(\mathrm{exp}(\beta ϵ)+1)`$ is the Fermi function. Using Eq. (5) the thermal conductivity $`\kappa (\omega )=D_{\mathrm{th}}(T)\delta (\omega )`$ can be computed directly. Here we focus on the case of $`0\mathrm{\Delta }1`$ and $`\alpha <\alpha _{\mathrm{crit}}`$hase93 , which is known to exhibit gapless spinonlike excitations. Results for $`D_{\mathrm{th}}(T)`$ of the $`XXZ`$ model are shown in Fig. 1 (thick dashed lines). For comparison Bethe-ansatz resultskluemper02  are included in the figure (solid lines). In the case of $`\mathrm{\Delta }=0`$, no approximations are necessary in the Jordan-Wigner approach and consequently, Jordan-Wigner and Bethe-ansatz results are identical. For $`\mathrm{\Delta }>0`$, the main observation is that the mean-field theory produces qualitatively the right picture of the temperature dependence of $`D_{\mathrm{th}}(T)`$. Both the slope of $`D_{\mathrm{th}}T`$ at low temperatures and the position of the maximum are well predicted. Deviations at high temperatures are due to the neglect of many-particle excitations in the mean-field approximation.
Bosonization -  In the continuum limit the physics of the anisotropic spin-$`1/2`$ chain at low energies is described by the Luttinger-liquid Hamiltonianschulz

$$H=\frac{1}{2}𝑑x\left(vK(_x\mathrm{\Theta })^2+\frac{v}{K}(_x\varphi )^2\right),$$
(6)
where $`\varphi `$ is a bosonic field in $`1+1`$ dimensions and $`\mathrm{\Theta }`$ is the dual field $`_x\mathrm{\Theta }=\frac{1}{K}_\tau \varphi `$. $`K`$ is the Luttinger parameter and $`v=(J\pi /2)\frac{\mathrm{sin}\gamma }{\gamma }`$ is the spinon velocity where the anisotropy is parametrized via $`\mathrm{\Delta }=\mathrm{cos}(\gamma )`$ here. The local current operator $`j(x)`$ is again given by the equation of continuity: $`_xj(x)=_th(x)`$. We obtain

$$j=v^2𝑑x_x\varphi (x)_x\mathrm{\Theta }(x).$$
(7)
The Drude weight follows from $`D_{\mathrm{th}}=\pi \beta ^2j^2/N`$. Thus we have to evaluate the two-point function $`j(x,\tau )j(0,0)`$, $`\tau `$ being the time variable. The computation is similar to the procedure for the susceptibility in Ref. eggert94, . We change to coordinates $`z=v\tau +ix`$ and $`\overline{z}=v\tau ix`$. By decomposing $`\varphi (z,\overline{z})=\phi (z)+\overline{\phi }(\overline{z})`$ into its chiral parts and using the respective two-point functions such as $`\phi (z)\phi (w)=(K/4\pi )\mathrm{ln}(zw)`$ we obtain

$$j(x,\tau )j(0,0)=2\frac{v^2}{(4\pi )^2}\left(\frac{1}{z^4}+\frac{1}{\overline{z}^4}\right).$$
(8)
Before performing the space integration the imaginary time direction is compactified by mapping the plane ($`z`$) into the strip ($`\zeta `$) using $`z(\zeta )=\mathrm{exp}(2\pi \zeta /\beta )`$ leading to the replacement $`v\tau \pm ix(v\beta /\pi )\mathrm{sin}(\pi \frac{v\tau \pm ix}{v\beta })`$ in Eq. (8). After the change of variables $`u=\mathrm{tan}(\pi \tau /\beta );w=i\mathrm{tan}(i\pi x/(v\beta ))`$ we finally find

$$D_{\mathrm{th}}(T)=\frac{\pi ^2}{3}vT.$$
(9)
This coincides with Klümper’s and Sakai’s analytic expressionkluemper02  for the low-temperature limit of the $`XXZ`$ model if the velocity $`v`$ is equal to $`v=(J\pi /2)\frac{\mathrm{sin}\gamma }{\gamma }`$. However, the result is more generally valid for models with the continuum limit given by the Luttinger-liquid Hamiltonian.

Exact diagonalization (ED) -  In this part we present our results for $`D_{\mathrm{th}}(T)`$ obtained by exact diagonalization for finite systems with $`N18`$. We start with the discussion of different values of the anisotropy $`\mathrm{\Delta }`$ at zero frustration. Figure 2 shows $`D_{\mathrm{th}}`$ for the isotropic case $`\mathrm{\Delta }=1`$, for a gapped, antiferromagnetic system ($`\mathrm{\Delta }=10`$) and in the ferromagnetic regime ($`\mathrm{\Delta }=1,2`$). While we show in Fig. 2(a) that we reproduce the results by Alvarez and Grosgros02  for system sizes of $`N14`$, our analysis extends this case to $`N18`$. This is due to exploiting both conservation of total $`S^z`$ and momentum $`k`$ in the exact diagonalization. By comparing with the curve obtained from the Bethe ansatzkluemper02  (solid line in Fig. 2(a)) it can be seen that for a system of $`N=18`$ sites the thermodynamic limit is reached for temperatures around $`T0.3J`$ for $`\mathrm{\Delta }=1`$. At low temperatures $`D_{\mathrm{th}}(T)`$ is exponentially suppressed due to the finite-size gap in the case of an even number of sites and divergent for an odd number. The latter is due to the degeneracy of the ground state in case of odd-numbered systems.
For the gapped, antiferromagnetic case we choose $`\mathrm{\Delta }=10`$, shown in Fig. 2 (b), where the finite-size effects of the two-spinon gap are small. The data for $`D_{\mathrm{th}}`$ are convergent for $`T3J`$, but substantial finite-size effects are still present in the vicinity of the maximum, i.e., at small temperatures compared to the two-spinon gap $`8.055126J`$cloizeaux66 . At low temperatures the thermal Drude weight $`D_{\mathrm{th}}(T)`$ is expected to be exponentially suppressed in the thermodynamic limit. In the inset of Fig. 2 (b) $`D_{\mathrm{th}}(T)`$ is plotted logarithmically versus $`1/T`$. If one fits $`D_{\mathrm{th}}(T)\text{exp}(\delta /T)`$ to the numerical data at low temperaturesT\_range  one finds $`\delta =8.056J`$ for $`N=18`$ sites with similar values of $`\delta `$ found for other $`N`$. This compares well to the two-spinon gapcloizeaux66 . Hence we conclude that mainly the elementary excitations contribute to the thermal conductivity at low temperatures.
In the ferromagnetic regime ($`\mathrm{\Delta }1`$) (results are shown for $`\mathrm{\Delta }=1`$ and $`\mathrm{\Delta }=2`$ in Fig. 2(c)) our main observation is that convergence with $`N`$ is very good at all temperatures. For example, we find for $`\mathrm{\Delta }=2`$ that the relative difference between the finite-size data for $`N=16`$ sites and $`N=17`$ is negligibly small, namely $`|D_{\mathrm{th}}^{N=17}(T)D_{\mathrm{th}}^{N=16}(T)|/D_{\mathrm{th}}^{N=17}(T)<0.008`$ for $`T>0.05J`$. If one extracts $`\delta `$ from a fit of $`D_{\mathrm{th}}(T)\mathrm{exp}(\delta /T)`$ to the numerical dataT\_range  for $`\mathrm{\Delta }=2`$, we find $`\delta 0.97J`$ which coincides with the one-triplet gap $`(\mathrm{\Delta }+1)J`$ that can easily be obtained from a spin-wave computation. The fast convergence is even more remarkable for the case $`\mathrm{\Delta }=1`$ where the numerical data are consistent with $`D_{\mathrm{th}}(T)T`$ at low temperatures. In addition, there is no qualitative difference between even- and odd-numbered systems due to the ferromagnetic nature of the interaction for $`\mathrm{\Delta }1`$.
Frustrated chain -  Now we turn to the case of nonzero frustration. Since $`[H,j]0`$ here, care has to be taken about off-diagonal matrix-elements of $`j`$ if degeneracies occur. However, since we use classification by momentum $`k`$ and $`S^z`$ degeneracies are lifted and do not play a crucial role. In Fig. 3 we show $`D_{\mathrm{th}}(T)`$ obtained numerically from Eqs. (4) and (5) for even system sizes with $`N18`$ and $`\alpha =0.35,\mathrm{\Delta }=1`$. A central result of this paper is that, while we observe a finite Drude weight at temperatures $`T>0`$ and all system sizes investigated, we still find a substantial reduction of the Drude weight with increasing system size at high $`T`$. This is in sharp contrast to the $`XXZ`$ model, where finite-size effects are small at high temperatures (see e.g. Fig. 2 (a)). These observations clearly point to a vanishing of the Drude weight in the thermodynamic limit for $`\alpha =0.35`$. However, the question of dissipationless thermal transport at arbitrary $`\alpha >0`$ remains to be studied in more detail.

Finally, we compare our mean-field approach with numerical results on finite systems and with nonzero frustration. In Fig. 4 we present the thermal Drude weight of systems with $`N=16`$ sites for $`\alpha =0,0.05,0.15`$ and $`\mathrm{\Delta }=1`$, i.e., in the gapless regime. As is obvious from this figure, there is a good agreement between ED and the mean-field approach regarding the temperature dependence of the thermal Drude weight. The general features of $`D_{\mathrm{th}}(T)`$ are a reduction of the absolute value of $`D_{\mathrm{th}}`$ on increasing $`\alpha `$, a shift of the position of the maximum to lower temperatures and thus a crossing of the curves for different $`\alpha `$ at low temperatures which are present in both the ED and mean-field results. Deviations at high temperatures are again understandable due to the neglect of many-particle excitations in the effective one-particle picture.
Conclusion -  We performed a detailed analysis of the thermal Drude weight for anisotropic and frustrated spin-$`1/2`$-Heisenberg chains by using mean-field theory, bosonization and ED. In the case of the $`XXZ`$ model we demonstrated the applicability of these techniques for computing the temperature dependence of the thermal Drude weight. Using ED we obtained results on finite systems of $`N18`$ sites for arbitrary values of the anisotropy $`\mathrm{\Delta }`$. In the ferromagnetic regime ($`\mathrm{\Delta }1`$) of the $`XXZ`$ chain the numerical data converge to the thermodynamic limit at arbitrary temperature for moderately small system sizes ($`N18`$). The analytical results compare well with the Bethe ansatzkluemper02  in the gapless regime of the $`XXZ`$ model and to our numerics in the case of the frustrated chain on finite systems. Our numerical data at $`\alpha =0.35`$ mark a clear difference between the integrable $`XXZ`$ case and the nonintegrable one at $`\alpha =0.35`$: while in the former case the Drude weight remains finite in the thermodynamic limit we have clear indications for a vanishing Drude weight at high temperatures in the latter case. Extended analysis of these findings for frustrated and dimerized chains and spin ladders will be the subject of a forthcoming paper.
Acknowledgments -  This work was supported by the DFG, Schwerpunktprogramm 1073, and by a DAAD-ANTORCHAS exchange program. We acknowledge helpful discussions with J.V. Alvarez, B. Büchner, C. Gros, and C. Hess.


text:
# Scaling at the onset of chaos in a network of logistic maps with two types of global coupling

## Abstract

We study a network of logistic maps with two types of global coupling, inertial and dissipative. Features of the clusterization process are revealed and compared, which are associated with presence of each type of couplings. For the parameter region near the onset of chaos we outline scaling properties in dynamical behavior of the network and illustrate them on the Kaneko phase diagrams.

*Institute of Radio-Engineering and Electronics of RAS,
Zelenaya 38, Saratov, 410019, Russia*

*Saratov State University,
Astrakhanskaya 83, Saratov, 410026, Russia*

Recently a notable interest is attracted to systems of globally coupled elements \[1-4\]. In particular, such models are discussed in context of description of neural networks, multi-modal lasers, electronic and chemical systems. They are regarded as perspective for nontrivial application in information processing to reproduce, at least in some respects, analogous processes in biological systems.

One particular class of networks, rather simple but possessing many nontrivial peculiarities in the dynamical behavior, was introduced by Kaneko . It is a set of logistic maps, each coupled to each other with equal strength. Studies of Kaneko and other authors revealed a reach phenomenology of dynamics of this model. In particular, a phenomenon called clusterization has been discovered – this is a spontaneous formation of groups of elements in process of dynamics of the network with instantaneous states exactly coinciding for elements relating to each group (cluster). A possibility of coexistence of attractors with different number of clusters at the same parameter values of the network gives an opportunity for information storage; control of them may be used for information processing.

In the original work of Kaneko and in many subsequent researches, a special concrete type of coupling was assumed: it tends to equalize instantaneous states of the interacting elements and is called the dissipative coupling . However, it is not the only possible type of interaction of the elements. As found from the renormalization group (RG) analysis of the onset of chaos in coupled period-doubling maps and in coupled map lattices, in general case a weak coupling introduced by smooth functions of state variables is a combination of two relevant types of coupling – inertial and dissipative \[5-8\]. Both types of coupling possess distinct properties in respect to the RG transformation and manifest different scaling regularities at the onset of chaos in coupled logistic maps. A necessity of consideration of globally coupled models with two types of coupling was mentioned e.g. in Refs.\[ 9, 10\], but neither convenient models, nor detailed studies of features of dynamics were reported there.

In the present paper we introduce a model of a network of logistic maps with dissipative and inertial global coupling. Some features of dynamics and structure of phase diagrams will be discussed and compared for the networks with these types of coupling. In particular, we outline properties of universality and scaling in dynamics of the globally coupled maps near the onset of chaos. Due to richness of dynamical phenomena, this region may be especially interesting for applications in the information processing.

Let us consider a network of coupled logistic maps governed by the following equations:

$$x_{n+1}(i)=(1\epsilon _2)f(x_n(i))\epsilon _1x_n(i)+\frac{\epsilon _1}{N}\underset{j=1}{\overset{N}{}}x_n(j)+\frac{\epsilon _2}{N}\underset{j=1}{\overset{N}{}}f(x_n(j)),$$
(1)
where $`f(x)=1\lambda x^2`$ is a nonlinear function corresponding to the logistic map, $`n`$ designates the discrete time, index $`i`$ enumerates elements of the network, $`N`$ is the total number of the elements, $`\epsilon `$ and $`\epsilon `$ are parameters of coupling.

Note that two last terms in Eq. (1) do not depend on the index $`i`$, in other words, they are the same for all elements of the network. Hence, they may be regarded as two mean fields responsible for two types of global coupling:

$`F_n^{(\text{1})}=N^1\underset{j=1}{\overset{N}{}}x_n(j)`$ and $`F_n^{(2)}=N^1\underset{j=1}{\overset{N}{}}f(x_n(j))`$. (2)

As follows from previous works on coupled logistic maps \[5-10\], the term with $`f`$($`x`$) in Eq.(1) and the mean field $`F^{(2)}`$ represent the pure dissipative coupling. The linear term and the term $`F^{(\text{1})}`$ correspond to a combination of inertial and dissipative coupling, but the dissipative contribution is rather small. Quantitatively, the coefficients $`\epsilon _1`$ and $`\epsilon _2`$ are expressed via the coefficients of inertial and dissipative coupling $`\epsilon _I`$ and $`\epsilon _D`$ as follows :

$$\epsilon _I=\epsilon _1,\epsilon _D=\epsilon _20.088\epsilon _1.$$
(2)
In the system with both types of coupling we expect to observe the same phenomenon of clusterization discovered by Kaneko for the dissipative case. Indeed, the main reason –coincidence of conditions of motion for elements with equal instantaneous states is valid irrespectively to the number of mean fields.

To analyze dynamics of the globally coupled network in dependence on the parameters let us follow the approach of Kaneko based on the concept of phases.

Been given a point in the parameter space ($`\lambda \epsilon `$, $`\epsilon `$) we consider an ensemble of identical independent networks, each with its own, randomly chosen initial conditions. Performing a sufficiently large number of iterations, we estimate a number of clusters been formed in each representative of the ensemble. Accounting statistics of the final states, we classify the phase as coherent (one-cluster state dominates in the ensemble), ordered (states with few clusters are of the highest probability), partially ordered (clusters with small and large numbers of elements appear with comparable probabilities), turbulent (number of clusters is of order of the total number of elements in the network). On the phase diagrams of Fig.1 the respective phases are indicated by symbols C, O, PO, and T. Panel (a) relates to the case of inertial coupling, this is cross-section of the parameter space by plane ($`\lambda `$ $`\epsilon _1=\epsilon _I,\epsilon _2=0.088\epsilon _I`$). Panel (b) corresponds to dissipative coupling and reproduce the diagram computed by Kaneko. It represents a cross-section of the parameter space by plane ($`\lambda `$ $`\epsilon _1=0_I,\epsilon _2=\epsilon _D`$) and is presented here for comparison. (Here we restrict ourselves by frames of a rough approach to the definition of phases, as in the original work of Kaneko, because it is sufficient to reach our aims and to demonstrate scaling. However, it has been shown recently that detection of presence or absence of clusterization, in a sense of exact coincidence of instantaneous states of elements, especially in chaotic domain, require much more subtle numerical technique and analysis .)

Observe that in the domain of subcritical dynamics of individual elements the network with inertial coupling demonstrates more reach dynamics than in the case of dissipative coupling. Presence of phases O and PO indicates that formation of multiple clusters and multistability takes place here, in domain of regular dynamics. Perhaps, this implies that the networks with inertial coupling should be considered as more perspective candidates for applications in information processing than those with dissipative coupling.

In Refs.\[5-10\] it was stated that coupled maps with inertial and dissipative coupling near the onset of chaos demonstrate a property of scaling. The network with global coupling may be thought as a set of elements coupled pairwise. Hence, analogous scaling regularities must be valid for the global coupling, and the statement may be formulated as follows.

Let us suppose that at the parameters $`\lambda \epsilon `$, $`\epsilon `$ we detect some phase for the ensemble of systems (1) with initial conditions randomly chosen from interval $`\left|x\right|<C`$. Then, for ensemble with random initial conditions from interval $`\left|x\right|<C\mathrm{/}\left|\alpha \right|`$ at the point $`\lambda +(\lambda \lambda )\mathrm{/}\delta `$, $`\epsilon _I\mathrm{/}\alpha `$, $`\epsilon _D\mathrm{/}2`$, the same kind of phase will be observed, but with doubled time scale of motion. Here $`\lambda `$ is the parameter value of the period-doubling accumulation in an individual element of the network (the logistic map), $`\alpha `$–2.5029…and $`\delta `$4.6692…are the Feigenbaum universal constants.

Straightly speaking, the above statement is true only asymptotically: smaller the neighborhood of the critical point ($`\lambda `$, 0, 0), more accurate the scaling property holds. But actually it works very well, even in sufficiently large scales. It may be seen from comparison of the phase diagrams of Fig 1 (a) and (c), (b) and (d). The bottom pictures present results obtained with rescaling of the parameters and initial conditions in accordance with the formulated rules.

In, conclusion, we have introduced a model network of logistic maps with two types of global coupling. We argued that the Kaneko classification of phases suggested for the network of dissipatively coupled maps remains valid for our generalized model. In contrast to the dissipative coupling, the inertial one ensures reach dynamics (multi-stability, states of multiple clusters – the phases O and PO) not only in supercritical (chaotic) domain of the individual elements, but also in subcritical domain and at the border of chaos, where the dynamics is more regular. This circumstance indicates that the networks with inertial coupling are more interesting for possible applications in information processing than those with dissipative coupling.

For the first time we have demonstrated scaling regularities in phase diagrams of the models with global coupling. Note significance of this observation. As follows from RG argumentation developed earlier for coupled maps and for coupled map lattices, the scaling property appears as an attribute of the universality class. It means that analogous dynamical regimes and the same regularities will be intrinsic not only to our logistic map model, but for globally coupled networks composed of any other period-doubling elements that relate to the Feigenbaum universality class (say, of forced nonlinear dissipative oscillators) . So, the results are expected to be common for globally coupled networks of different physical nature (electronics, biology, economics, etc.).

The work was supported, in part, by RFBR (grant No 00-02-17509) and by Ministry for Industry, Science and Technology of Russian Federation (State Contract No 40.020.1.1.1168).

FIGURE CAPTION

Fig.1. Kaneko’s phase diagrams in cross-section of the parameter space of the globally coupled network of logistic maps by planes of pure inertial coupling (a,c) and pure dissipative coupling (b,d). The bottom pictures are built with the use of scaling rules formulated in the main text. Observe similarity of the diagrams (a) and (c), (b) and (d). The symbols designate Kaneko’s phases: C – coherent, O – ordered, PO – partially ordered, T – turbulent.