## Abstract

We study numerically the counterpropagating vector solitons in SBN:60 photorefractive crystals. A simple theory is provided for explaining the symmetry-breaking transverse instability of these solitons. Phase diagram is produced that depicts the transition from stable counterpropagating solitons to bidirectional waveguides to unstable optical structures. Numerical simulations are performed that predict novel dynamical beam structures, such as the standing-wave and rotating multipole vector solitonic clusters. For larger coupling strengths and/or thicker crystals the beams form unstable self-trapped optical structures that have no counterparts in the copropagating geometry.

©2004 Optical Society of America

## 1. Introduction

Whenever a new class of solitary waves [1] is discovered, the question of their stability is raised. It pertains to structural stability of a homogeneous wave with regard to symmetry-breaking transverse and modulational instabilities [2] during its evolution. Loosely speaking, transverse means instabilities induced in the plane transverse to the propagation direction, and modulational denote spatial changes to the wave as it evolves. Recently introduced counterpropagating (CP) vector solitons in photorefractive (PR) media [3,4] present no exception.

The formation and interactions of spatial screening solitons [5] have been studied mostly in the copropagation geometry, with few exceptions [3,4,6,7]. Steady-state CP solitons were considered theoretically in one transverse dimension (1D), in Kerr and local PR media. In addition, Ref. [4] introduces a time-dependent model for the generation of CP solitons, and Ref. [7] contains an experimental observation of stable CP solitons in an SBN:60 crystal, in the form of narrow stripes.

Here we present theoretical and numerical study of 2D CP vector solitons in a similar crystal, display symmetry-breaking instabilities of such solitons when the propagation distance and/or the coupling strength are varied, and demonstrate the rich dynamics of different self-trapped beam structures that can form in the crystal. We formulate a simple theory that explains the split-up transverse instability of CP solitons as a first-order phase transition and display the corresponding threshold curve in the parameter plane. We present evidence of a second phase transition when, due to modulational instabilities, the steady-state asymmetric waveguides loose stability to time-dependent periodic and aperiodic optical structures. We perform numerical simulations of CP beams in PR media with time-dependent nonlinearity, and predict novel stable dynamical beam structures that possess no counterparts in the copropagating geometry.

The study of CP beams is performed in a geometry that resembles an actual experimental setup. A Nd:YAG laser beam at 532 nm is split and focused onto the opposite faces of an SBN:60 crystal. The beam components are assumed to be colinear and partially coherent in the medium, which in experiments can be achieved by reflecting one component off the vibrating mirror or by propagating it through a diffuser. The *c* axis of the crystal is placed perpendicular to the propagation direction. To exploit the dominant electro-optic component *r*_{33}
≈200 pm/V of a typical SBN sample, the incident laser beams are linearly polarized, parallel to the c axis. A DC electric field of the order of 1 kV/cm is applied across the crystal, along the *c* axis, and the crystal is illuminated by a uniform white light, to create artificial dark conductivity. Such a geometry is appropriate for the formation of screening spatial solitons.

The slowly varying beam components *F* and *B* counter-propagate in the crystal in the *z* direction, perpendicular to the *c* axis, which is also the *x* axis of the coordinate system. The space charge field generated in the crystal couples to the electro-optic tensor, giving rise to a change in the index of refraction, of the form Δ*n*=-${{n}_{\mathit{0}}}^{\mathit{3}}$
*r*_{eff}
*E*/2, where *n*_{0}
is the unperturbed index, *r*_{eff}
is the effective component of the electro-optic tensor, and *E* is the *x* component of the space charge field. The optical field is given as the sum of the amplitudes [*Fexp* (*ikz*)+*Bexp* (-*ikz*)+*cc*]/2, *k* being the wave vector in the medium, so that the total light intensity (uniform plus beams) is modulated in the *z* direction

where *I*_{0}
=|*F*|^{2}+|*B*|^{2} is the sum of individual beam intensities, *ε* is the degree of coherence of the beams, *m*=2*FB*
^{*}/(1+*I*_{0}
) is the modulation depth, and Δ*ϕ* is the phase difference between forward and backward beams. The intensity is measured in units of the background light intensity. It modulates the space charge field as well, which acquires the form

where **E**
_{0} is the homogeneous part of the space charge field, not to be confused with the external electric field, and **E**
_{1} is the first Fourier component, proportional to *ε*. It is **E**
_{0} that screens the external field, and **E**
_{1} is responsible for the formation of gratings (with the wave-number *2k*) in the index of refraction along the *z* direction.

The propagation of beams in the crystal is governed by the paraxial wave equations

where Δ is the transverse Laplacian, Γ=(*kn*
_{0}
*x*
_{0})
^{2}*r*_{eff}*E*_{e}
is the coupling strength, *E*_{e}
being the external electric field, and *E*_{0}
and *E*_{1}
are the dimensionless *x* components of the space charge field. The equations are put in the dimensionless form using the rescaling (*x*,*y*)→(*x*/*x*_{0}
, *y*/*x*_{0}
), *z*→*z*/*L*_{D}
, (*F*,*B*)→(*F*,*B*) exp (-*i*Γ*z*). Here *x*_{0}
is the typical beam waist and *L*_{D}
=${{\mathit{2}\mathit{\text{kx}}}_{\mathit{0}}}^{\mathit{2}}$
is the diffraction length. Propagation equations can be put in a universal dimensionless form that contains no parameters or coupling constants. All the parameters are then hidden in the scaling quantities and the initial and boundary conditions. We prefer the form given here, with one explicit intensive control parameter Γ. The corresponding extensive control parameter is the crystal length *L*.

In a local, isotropic approximation to the space charge field one assumes the following relaxation-type dynamics of its components [4,8]:

where *τ*=*τ*_{0}
/(*1*+*I*_{0}
) is the relaxation time of the crystal, which also depends on the total intensity. The nonlocal, anisotropic theory of CP solitons [9] suggests that for the geometries of interest here the coherent aniso-description offers similar results to the incoherent iso-description. The gratings induced in the *z* direction affect the propagation of beams little, in contrast to the standard two-wave mixing in PR media, when the gratings wavevector is aligned with the c axis. Therefore, we will consider only incoherent beams, *ε*=0, so that *E*_{1}
is absent. The propagation equations are then coupled only through the dependence of *E*_{0}
on *I*_{0}
.

## 2. Theory of beam displacement

To explain the behavior of beams observed in numerical simulations, we adopt the particle point of view on solitons, whereby they are considered as bundles of focused rays boring an optical path through the crystal [10]. The trajectory of a soliton is represented by the expectation value of its transverse coordinates:

where *A* is the soliton envelope and *I*_{t}
is the total transverse intensity, *I*_{t}
=∫∫|*A*|^{2}
*dxdy*, which plays the role of an effective mass. The motion of the soliton particle is governed by the Hamilton equations for the center of mass

and a similar pair of equations for ‹*y*› and *p*_{y}
, where the conjugate momenta *p*_{x}
and *p*_{y}
are represented by the optical direction cosines, and the potential *V*=Γ*E*_{0}
is the (negative) change in the refractive index. This amounts to viewing an optical soliton as a particle of mass proportional to the intensity, moving in a potential created by the change in the refractive index, which is caused by the soliton itself. Such a point of view is akin to the usual mechanics, except that the role of time is played by *z*, and the “dynamics” is in 2D.

Particle picture is obtained from the ray optics. Quantization of the ray optics leads to the paraxial wave optics. In fact, the construction of rays from the wave optics is equivalent to the construction of classical mechanics from the quantum mechanics. The paraxial wave equation then corresponds to the time-dependent Schrödinger equation,

where **p**=-*i*∇is the transverse gradient and *ρ*=(*x*,*y*) is the transverse position vector. The expectation values are determined using *A* for the wave function. Hence, one can view an optical soliton as a “quantum mechanical” object in the transverse plane, whose wave function is given by the slowly varying envelope of the beam, the potential by the induced change in the refractive index, and the momentum by the transverse gradient. The transition from the wave picture to the particle picture, i.e. the geometrical optics approach to spatial solitons, is well defined, as the size of beams is much larger than the wavelength, diffraction is absent, and incoherent beams are considered.

The transverse displacement of the trajectory along *z* axis obeys the Ehrenfest theorem [11],

The idea is then to express a transverse displacement of the center of mass **d**=(*d*_{x}
,*d*_{y}
) in terms of the expectation values of the gradient of the space charge field. The expectation values are evaluated for the shifted and unshifted states *A*_{d}
=*A* (*ρ*+**d**, *z*) and *A*_{0}
=*A* (*ρ*, *z*), and then subtracted:

In evaluating the expectation values one utilizes the relation between the shifted and unshifted states *A*_{d}
=exp(**d**·∇)*A*_{0}
=exp(*i*
**d·p**)*A*_{0}
. Assuming that **d** is small one obtains

where

and *ê* is the unit vector in the transverse plane. The equation for **d**(*z*) is a harmonic oscillator equation, with the solution *d*(*z*)=*a* sin(*K*
^{1/2}
*z*)+*b*cos(*K*
^{1/2}
*z*), where the constants *a* and *b* are fixed by the boundary conditions. In the case of head-on collision of CP beams, the lowest order steady state mode for the *F* beam has **d**(0)=0 at the entrance face of the crystal and **d**(*L*)≠0 at the exit face. This is similar to a vibrating elastic string whose one end is fixed and the other free. The lowest order mode for the *B* beam is the mirror image. For such a state there exists an obvious threshold condition

To see the form of this threshold condition in the (*L*, Γ) plane one must include Γ. *K* seems to be linear in Γ. However, the integral in Eq. (11) carries another Γ. It comes from the scaling used to write dimensionless propagation equations. Had we used the scaling where no Γ appears, the threshold condition would appear the same, (*L′K′*
^{1/2})_{c}=*π*/2, but the quantities in the primed coordinates would be connected to the unprimed through *L′*=ΓL, and *K*=Γ^{2}
*K′*. Since *K′* does not contain Γ, it can be transferred to the other side of the threshold equation, and the threshold line in the (*L*, Γ) plane acquires the form (Γ*L*)_{c}=*const*. Thus, the theory predicts a simple “*pV*=*const*.” equation of state, with unshifted solitons existing below the threshold curve, and shifted solitons, or bidirectional waveguides, existing above the curve.

A note of caution is required here. The theory developed thus far applies to steady-state situations. No mention is made of the “true” dynamics. Time enters into the picture through the explicit dependence of *K* on *t*. A variety of dynamical behaviors becomes allowed. Analytical analysis becoming prohibitive [12], we resort to numerical analysis, in both 1D and 2D.

## 3. Numerical simulations

The spatial propagation equations and the temporal equations for the space charge field are solved concurrently. The numerical procedure consists in solving Eqs. (4) for the components of the space charge field in time, with the light fields obtained at every step as guided modes of the induced common waveguide. This is achieved by an internal spatial relaxation loop, i.e. nested within the temporal loop, that utilizes a beam-propagation method for the right- and left-propagating components. The spatial loop is iterated until convergence, and the temporal loop is advanced for a time step. The convergence in the temporal loop signifies that steady states are found, however it need not be reached. In that case time-dependent, dynamical states are observed. The procedure is described in Refs. [12, 13].

To capture the transition from a CP soliton to a waveguide clearly, we consider head-on collision of two identical Gaussian beams. In the absence of the other, each beam focuses into a soliton. We are interested in what happens when they are both present, and when the coupling constant Γ and the crystal length *L* are both varied. We note no qualitative difference between the 1D and 2D behavior, and present the 1D case. The situation is depicted in Fig. 1. It is seen that in the plane (*L*, Γ) of control parameters there exists a critical curve below which stable CP solitons exist. At the critical curve a new type of solutions appears, in which the two components do not overlap anymore, but split and cross each other. A few examples are depicted in the insets to Fig. 1. As the beams split, a portion of each beam remains guided by the other, therefore we term these solutions bidirectional waveguides. Both the solitons and the waveguides are steady-state solutions. According to the theory, one can easily define higher order steady-state modes. However, they are not easily detectable, since the system may become dynamically unstable before reaching them. As one moves away from the critical curve, into the region of higher couplings and longer crystals, a new critical curve is approached, where the steady-state waveguides loose stability. The second critical curve is also drawn in Fig. 1, and is similar in shape to the first one. The shape of these curves suggests an inverse power law dependence, in accordance with the theory presented here. In fact, the analytical expression for the fitting curves drawn in Fig. 1 contains a constant term close to Γ
_{th}
=2, and the terms linear and quadratic in 1/*L*. This suggests van der Waals-type corrections to the “equation of state” Γ*L*=*const*., and implies that the simple linear theory presented here is insufficient to account for all the varied behavior of the full nonlinear system.

At and beyond the second critical curve dynamical solutions emerge. The time dependence varies from periodic to aperiodic. A richer dynamical behavior is observed in 2D, as compared to 1D, since there one has a larger phase space at disposal, and can launch beams carrying angular momentum and/or topological defects in their structure. Some 2D examples are presented below.

Common to all simulations are the following data: diffraction length 5.55 mm, transverse scaling length 10 µm, laser wavelength 532 nm, electro-optic coefficient 180 pm/V, bulk refractive index *n*_{0}
=2.35. The value of the external electric field (of the order of 1 kV/cm) is used to fine-tune the coupling strength. The propagation lengths are given in units of *L*_{D}
, and the time in units of *τ*_{0}
. All initial fields are head-on. All the simulations are movie files (MOV), depicting temporal evolution of the optical field in either the transverse (*x,y*), or the longitudinal (*y,z*) plane.

In 2D, below the first critical curve, stable CP solitons are observed. In Fig. 2 we present a stable displaced bidirectional waveguide just above the threshold curve, for Γ=9.2, *L*=1.8, obtained by colliding two Gaussian beams of 15 µm FWHM. As seen, the Gaussians settle fast into a quasi-stable CP soliton, remaining steady until *t*≈160*τ*, at which time they split within a few cycles to new transverse positions. A part of each beam remains at the old position, being guided by the other beam. This is in accordance with theory, which predicts that the center of mass of each beam will become transversally displaced at the exit face of the crystal as soon as the threshold is reached. For this rotationally symmetric geometry the direction in which the beams split is random. For rotationally non-simmetric beam configurations there exists a preferential direction, and this is the direction in which K grows the fastest. According to the threshold condition, the value of the critical crystal length is then the shortest, and, as the threshold is reached, the instability will preferentially grow in that direction.

The enhanced stability of dipole beams, as compared to other beam structures in PR crystals, has been noted [14]. This applies to CP dipoles as well. In an earlier publication [15] we reported symmetry-breaking displacement of CP dipole-mode vector soliton (dipole beam plus the fundamental), similar to the displacement of simple CP solitons. Here we display the robust nature of two CP dipole beams, placed either parallel-parallel or parallel-perpendicular to the external field (Fig. 3). The initial fields are two incoherent dipole pairs, made of anti-phased Gaussian beams. Within a couple of cycles in each geometry the dipole beams reach steady state, in the form of two bean-like deformed dipoles. They remain very stable, in comparing with other beam compositions considered here.

Interest in optical beams carrying angular momentum and topological defects in the structure has risen recently, owing to their useful features for pattern formation. In contrast to dipole-mode vector solitons, beam compositions containing optical vortices as components are noted for instability [16,17]. In the copropagating geometry a vector beam made of a fundamental Gaussian and a vortex beam disintegrates into a deformed two-peak central beam and a dipole beam. The fragmentation of the vortex component can proceed in more than two filaments, depending on the charge and the size of the vortex. An interesting theoretical difference is noted in the copropagation geometry: in the anisotropic modeling of PR nonlinearity, which is closer to the experimental situation, the breakup proceeds much faster than in the isotropic modeling. While the anisotropic vortex breaks within one *L*_{D}
, the isotropic vortex can propagate for tens of *L*_{D}
’s before disintegrating.

As mentioned, in the CP case there is not much difference between the incoherent isotropic modeling and the coherent anisotropic modeling. Hence we observe that vortices generally disintegrate within one *L*_{D}
. In addition, owing to time-dependent nature of the isotropic model employed, we observe interesting dynamical effects. In the next few figures we report different cases of two colliding vortices: a standing-wave structure when the vortices carry the same charge, a stable rotating structure when the charges are the opposite, and various unstable outcomes of vortex-vortex collisions.

Figure 4 depicts standing-wave beams that are the result of the breakup of two identical CP vortices of topological charge +1. From (a) to (c) only one parameter is changed, the crystal thickness is increased from 0.8 to 1.1 to 1.5. After some time the vortices break in each case, and it is seen that the stable output changes from a tripole to a quadrupole. In Fig. 4(b) the beam first breaks into a meta-stable quadrupole, which eventually reverts to the stable tripole configuration. Such multipole clusters – dipoles, tripoles, quadrupoles, etc. – constitute the basic stable breakup modes of colliding vortices.

Figure 5 presents a stable rotating beam structure resulting from the collision of two CP vortices with opposite charges. For the chosen set of parameters the beams continue to rotate indefinitely. They can be considered true rotating propeller solitons [18], as they rotate in time and can have more than two blades. Input vortices are the same as in Fig. 4(a), only the sign of the charge of vortex *B* is reversed. The beams co-rotate, and such a state in time represents a limit-cycle. It can not be accessed by the usual steady-state theory of vector solitons [14, 18].

The initial beam structure of Fig. 5 is used to present spatio-temporally unstable states in Fig. 6(a), only the crystal thickness is increased to 1.5. By increasing *L*, or Γ, unstable situations arise. After a long quasi-periodic regime the system jumps to a disordered state. In our limited time of observation no repetition of the mode structure is observed. Instabilities can also be reached by changing the initial composition of beams. Initial conditions in Fig. 6(b) are the same as in Fig. 6(a), only the initial vortices are made wider. The system reverts to a periodic motion, however not a simple rotation of a stable structure. Finally, another unstable state is displayed in Fig. 6(c), it is still with wide initial vortices, but with a smaller coupling strength and crystal thickness. When these beams are made narrower, as in Fig. 5(a), the system similarly rotates (not shown). In the wider case the beams perform a motion similar to Fig. 6(a). However, it should be noted that the beams in that figure have acquired a higher-order transverse mode structure before disintegrating, namely a ring plus a central peak. In Fig. 6(c) they attain only a simple ring structure. Nonetheless, the dynamics on the ring looks similar, and the final state is a disordered spatiotemporal mode structure.

In summary, we have considered various counterpropagating self-trapped beam structures in an SNB:60 photorefractive crystal. A time-dependent model for the beam propagation and interaction is treated numerically, and various self-focused solutions presented, including pure CP vector solitons. A symmetry-breaking transverse instability of these solitons is noted, and a simple theory provided for explaining such a transition. A corresponding phase diagram is produced that depicts the transition from stable CP solitons to bidirectional waveguides and unstable beam structures. Simulations are performed displaying novel dynamical beam structures, such as the standing-wave and rotating vector solitonic clusters. For larger coupling strengths and thicker crystals the beams form unstable optical structures that have no analogues in the copropagating geometry.

## Acknowledgments

Work at the Institute of Physics, Belgrade is supported by MNTR, under grants OI 1475 and 1478. Work at the Institute of Applied Physics, WWU Münster was in part supported by DFG. MT and PM are supported by FNRS. Work at ULB was also funded in part by the IAP program of the Belgian government.

## References and links

**1. **For an overview see the Special Issue on solitons, Ed. M. Segev, Opt. Phot. News 13, No. 2 (2002).

**2. **Y. S. Kivshar and D. E. Pelinovsky, “Self-focusing and transverse instabilities of solitary waves,” Phys. Rep. **331**, 117–195 (2000). [CrossRef]

**3. **O. Cohen et al., “Spatial vector solitons consisting of counterpropagating fields,” Opt. Lett. **27**, 2013 (2002). [CrossRef]

**4. **M. Belic et al., “Self-trapped bidirectional waveguides in a saturable photorefractive medium,” Phys. Rev. A **68**025601, (2003).

**5. **S. Trillo and W. Torruellas eds., *Spatial Solitons* (Springer, New York, 2001).

**6. **M. Haelterman, A. P. Sheppard, and A. W. Snyder, “Bimodal counterpropagating spatial solitary-waves,” Opt. Commun. **103**, 145 (1993). [CrossRef]

**7. **O. Cohen et al., “Collisions between optical spatial solitons propagating in opposite directions,” Phys. Rev. Lett.89, 133901 (2002); “Holographic solitons,” Opt. Lett.27, 2031 (2002). [CrossRef] [PubMed]

**8. **L. Solymar, D. J. Webb, and A. Grunett-Jepsen, *The physics and applications of photorefractive materials*, (Clarendon Press, Oxford, 1996).

**9. **M. Belić et al., “Anisotropic nonlocal modelling of counterpropagating photorefractive solitons,” Opt. Lett.29 xxxx (2004).

**10. **M. Belić, A. Stepken, and F. Kaiser, “Spatial screening solitons as particles,” Phys. Rev. Lett. **84**, 83 (2000). [CrossRef]

**11. **P. R. Holland, *The quantum theory of motion*, (University Press, Cambridge, 1995).

**12. **O. Sandfuchs, F. Kaiser, and M. R. Belić, “Self-organization and Fourier selection of optical patterns in a photorefractive feedback system,” Phys. Rev. A **64**, 063809 (2001). [CrossRef]

**13. **M. Belić, J. Leonardy, D. Timotijević, and F. Kaiser, “Spatio-temporal effects in double phase conjugation,” J. Opt. Soc. Am. **B12**, 1602 (1995).

**14. **J. J. Garcia-Ripoll et al., “Dipole-mode vector solitons,” Phys. Rev. Lett.85, 83 (2000);W. Krolikowski et al., “Observation of dipole-mode vector solitons,” Phys. Rev. Lett.85, 1424 (2000). [CrossRef]

**15. **K. Motzek et al., “Dynamic counterpropagating vector solitons in self-focusing media,” Phys. Rev. E68 06xxxx (2003). [CrossRef]

**16. **M. R. Belić et al., “Isotropic vs. anisotropic modeling of photorefractive solitons,” Phys. Rev. A **65**, 066609 (2002).

**17. **G. F. Calvo et al., “Two-dimensional soliton-induced space charge field in photorefractive crystals,” Opt. Commun. **227**, 193 (2003). [CrossRef]

**18. **T. Carmon et al., “Rotating propeller solitons,” Phys. Rev. Lett. **87**, 143901 (2001). [CrossRef] [PubMed]