text:
# POSSIBLE MEASUREMENTS OF GPDs AT COMPASSTalk given at the workshop “Future physics at COMPASS”, Cern (2002)

## 1 Goal of an experiment with the high energy muon beam

The Generalized Parton Distributions (GPDs) provide a unified description of the nucleon. As explained by M. Diehl , they interpolate between the parton distributions and the hadronic form factors. Experimentally the GPDs can be accessed in exclusive measurements such as Hard Exclusive Meson ($`\rho ,\pi \mathrm{}`$) Production (HEMP) and Deeply Virtual Compton Scattering (DVCS). The latter reaction is the simplest from the theoretical point of view but also the most difficult experimentally because one has to select perfectly the final state (one lepton, one proton and one photon) among all the possible reactions. In pratice we can by now start the investigation of Meson Production and we foresee an upgrade of the COMPASS setup for DVCS measurement.

Experiments have already been undertaken at very high energy with the HERA collider to study mainly the gluon GPDs at very small $`x_{Bj}`$ ($`10^2`$). Larger values of $`x_{Bj}`$ have been investigated in fixed target experiments at JLab (6 GeV, with plans for an upgrade at 11 GeV) and HERMES (at 27 GeV). The experimental program using COMPASS at CERN (at 100 and/or 190 GeV) will enlarge the kinematical domain to a large range of $`Q^2`$ and $`x_{Bj}`$ (1.5 $`Q^2`$ 7 GeV2 and 0.03 $`x_{Bj}`$ 0.25) (see Fig. 1). A large range in $`Q^2`$ is required to control the factorisation in a hard, pertubatively calculable amplitude and a soft amplitude which is parametrized by the generalized parton distributions $`H,E,\stackrel{~}{H},\stackrel{~}{E}`$. The GPDs depend on three kinematical variables: $`x`$ and $`\xi `$ parameterize the longitudinal momentum fractions of the partons, while $`t`$ relates to the transverse momentum transfer.

Since the theoretical proof of factorization assumes that the transfer $`t`$ is finite (that is $`t/Q^20`$) , we consider in the following $`|t|`$ smaller than 1 GeV2. Another condition of factorization concerns the helicity of the virtual photon. In case of Hard Exclusive Meson Production it is mandatory to impose that the virtual photon be longitudinal in order to select the perturbative gluon exchange. Experimentally we should consider Rosenbluth separation for $`\pi ^0`$ production, while for $`\rho ^0`$ production we can select longitudinal $`\rho ^0`$s through the angular distribution of the decay products and assume the s-channel helecity conservation. Hard Exclusive Meson Productions seem more complex to analyze as they contain non perturbative information on both the target and the produced meson. Nevertheless they offer the possibility to disentangle different GPDs (vector meson production depends on $`H`$ and $`E`$ only; pseudo-scalar production depends on $`\stackrel{~}{H}`$ and $`\stackrel{~}{E}`$ only) and to separate contributions from different flavors. Forward differential longitudinal $`\rho _L^0`$ electroproduction cross section measurements which provide the largest counting rates, have already been undertaken and are presented in Fig. 2 as a function of c.m. energy W for three values of $`Q^2`$ (5.6, 9 and 27 GeV2). The theoretical curve is an incoherent sum of the quark and gluon contributions . No measurement has been done at $`x_{Bj}`$ larger than 0.05. Thanks to the expertise of the NMC collaboration for these absolute measurements we will explore a larger domain in $`x_{Bj}`$, $`Q^2`$ and $`t`$ with the muon beam available at CERN.

Deeply virtual Compton scattering is accessed by photon lepto-production: $`lpl^{}p^{}\gamma `$. In this reaction, the final photon can be emitted either by the leptons (Bethe-Heitler process) or by the proton (genuine DVCS process). If the lepton energy is large enough (see Fig. 3 with $`E_\mu `$ = 190 GeV, $`Q^2`$= 4 GeV2, $`x_{Bj}`$=0.1), the DVCS contribution dominates over the BH contribution so that the cross section is essentially the square of the DVCS amplitude which, at leading order, has the form:

$$T^{DVCS}_1^{+1}\frac{H(x,\xi ,t)}{x\xi +iϵ}𝑑x\mathrm{}𝒫_1^{+1}\frac{H(x,\xi ,t)}{x\xi }𝑑x\mathrm{}i\pi H(\xi ,\xi ,t)\mathrm{}.$$

(where $`\xi x_{Bj}/2`$ and $`t`$ are fixed by the experiment). At smaller lepton energy (see Fig. 3 with $`E_\mu `$ = 100 GeV and same values of $`Q^2`$ and $`x_{Bj}`$ as above), the interference between BH and DVCS becomes large and offers a unique opportunity to study Compton scattering amplitude including its phase. A careful analysis of the dependence of the cross section on the azimuthal angle $`\varphi `$ between the leptonic and hadronic planes and on $`Q^2`$ allows one to disentangle higher-twist effects and to select the real or imaginary parts of the DVCS amplitude (see the details in the previous text of M. Diehl ). If a longitudinally polarized lepton beam and an unpolarized target are used, the angular analysis and the $`Q^2`$ dependence of the cross section difference $`\sigma (e^{})\sigma (e^{})`$ allow one to select the imaginary part of the DVCS amplitude and thus the GPDs at the specific values $`x=\xi `$. This study is being investigated at HERMES  and JLab . If two muon beams of opposite charge and polarization are used, the angular analysis and the $`Q^2`$ dependence of the sum of cross sections $`\sigma (\mu ^+)+\sigma (\mu ^{})`$ allow also one to select the imaginary part of the DVCS amplitude. Moreover the same method applied to the difference of cross sections $`\sigma (\mu ^+)\sigma (\mu ^{})`$ allows one to select the real part of the DVCS amplitude which, for a given $`\xi `$, is sensitive to the complete dependence on $`x`$ of the GPDs. The deconvolution (over $`x`$) of this formula to extract the GPDs is not yet clearly solved, but comparison to model predictions can easily be made. It is clear that the muon beam of high energy at CERN can offer many possibilities in order to investigate the many-faceted problem of the GPDs knowledge.

Figure 4 shows the azimuthal distribution of the charge asymmetry which can be measured at COMPASS and the strong sensitivity to two different models . The first one is based on a simple parametrization of the GPDs: $`H^f(x,\xi ,t)=H^f(x,\xi ,0)F_1^f(t)/2`$ where $`F_1^f(t)`$ represents the elastic Dirac form factor for the quark flavor $`f`$ in the nucleon. The second one  relies on the fact that the GPDs measure the contribution of quarks with longitudinal momentum fraction $`x`$ to the corresponding form factor as is suggested by the sum rule:

$$_1^{+1}H^f(x,\xi ,t)𝑑x=F_1^f(t).$$

As one can associate the Fourier transform of form factors with charge distributions in position space, one can expect that the GPDs contain information about the distribution of partons in transverse position space. In fact it has been demonstrated that, when $`t`$ is purely transverse which amounts to $`\xi =0`$, then $`H(x,0,t)`$ is the Fourier transform of the probability density to find a quark with momentum fraction $`x`$ at a given distance from the center of momentum in the transverse plane. Qualitatively one expects that quarks with a large $`x`$ come essentially from the small valence “core" of the nucleon, while the small $`x`$ region should receive contributions from the much wider meson “cloud". Therefore one expects a gradual increase of the $`t`$-dependence of $`H(x,0,t)`$ as one goes from larger to smaller values of $`x`$. This suggests the parametrization: $`H(x,0,t)=q(x)e^{t<b_{}^2>}=q(x)/x^{\alpha t}`$ where $`<b_{}^2>=\alpha ln1/x`$ represents the increase of the nucleon transverse size with energy. The domain of small $`x_{Bj}`$ reached at COMPASS is related to the observation of see quarks or meson “cloud" or also gluons and it provides a large sensitivity to this three-dimensional picture of partons inside a hadron.

## 2 General requirements for COMPASS

The highest luminosity reachable at COMPASS is required to investigate these exclusive measurements. The experiment will use 100-190 GeV/c muons from the M2 beam line. Limits on radio-protection in the experimental hall imply that the maximum flux of muon to be expected is of $`210^8`$ muons per SPS spill (5.2s spill duration, repetition each 16.8s). Under these circumstances, we can reach a luminosity of $`=510^{32}`$ cm-2s-1 with the present polarized 6LiD or NH3 target of 1.2 meter long, and only $`=1.310^{32}`$ cm-2s-1 with a future liquid hydrogen target of 2.5 meter long.

In order to get useful cross sections with positive and negative muon beams, it is necessary to perform a precise absolute luminosity measurement. This has already been achieved by the NMC Collaboration within a 1% accuracy . The integrated muon flux was measured continuously by two methods: either by sampling the beam with a random trigger (provided by the $`\alpha `$ emitter Am241) or by sampling the counts recorded in 2 scintillators hodoscope planes used to determine incident beam tracks. The beam tracks were recorded off-line, in the same way as the scattered muon tracks to determine exactly the integrated usable muon flux.

Moreover $`\mu ^+`$ and $`\mu ^{}`$ beams of 100 GeV energy, with the same and as large as possible intensity as well as exactly opposite polarization (to a few %) are required. The muons are provided by pion and kaon decay and are naturally polarized. The pions and kaons come from the collision of the SPS 400 GeV proton beam on a Be primary target. A solution was proposed by Lau Gatignon . It consists in: 1) selecting 110 GeV pion beams from the collision and 100 GeV muon beams after the decay section in order to maximize the muon flux; 2) keeping constant the collimator settings which define the pion and muon momentum spreads (both the collimator settings in the hadron decay section and the scrapper settings in the muon cleaning section) in order to fix the $`\mu ^+`$ and $`\mu ^{}`$ polarizations at exactly the opposite value ($`P_{\mu ^+}=0.8`$ and $`P_\mu ^{}=+0.8`$); 3) fixing $`N_\mu ^{}`$ to $`210^8`$ $`\mu `$ per SPS spill with the longest 500mm Be primary target; 4) using a shorter target to find $`N_{\mu ^+}`$ close to $`210^8`$ $`\mu `$ per SPS spill.

This paragraph presents the experimental procedure to select the exclusive HEMP or DVCS channel and the difference equipments that are required. They are mostly part of the existing high resolution COMPASS spectrometer: muon detection which insures a good resolution in $`x_{Bj}`$ and $`Q^2`$, meson detection and identification in RICH or photon detection in calorimeters of good energy and position resolutions to allow two photons separation. The COMPASS spectrometer intercepts only forward outgoing particles (until 10 degrees) and the photon or meson detection limits the experiment to small $`x_{Bj}`$ values ($`x_{Bj}0.15`$). At these high energies the complete final state, including the low energy recoiling proton, needs to be detected because missing mass techniques are not efficient due to the experimental resolutions (the resolution in missing mass which is required is $`(m_p+m_\pi )^2m_p^2=0.25`$ GeV2 and the experimental resolution which can be achieved is larger than 1 GeV2). Consequently the high resolution COMPASS spectrometer needs to be completed by a recoil detector to measure precisely the proton momentum and exclude other reactions under high luminosity conditions. In the next section we will try to by-pass the necessity of a recoil detector to investigate the cleanest channel: $`\mu p\mu p\rho ^0`$ where $`\rho ^0`$s are identified through their decay in two charged pions accurately measured in the forward COMPASS spectrometer.

## 3 A pragmatic solution with the present setup

With the present COMPASS setup we can undertake Hard Exclusive $`\rho ^0`$ Production (with the largest cross section) as we can benefit of the good expertise for this measurement in the previous NMC and SMC experiments and so take the advantage to produce data in this field as soon as possible. At the same time we can also investigate Hard Exclusive $`\pi ^0`$ Production to see the limit of this setup.

In this context the constraints and limits of these experiments are the following ones:
1) With the present polarized 6LiD or NH3 target the production occurs on quasi-free nucleons in the nucleus or in the coherent scattering on the nucleus.
2) Luminosity determination is needed, but can be realized indirectly by measuring at the same time unpolarized Deep Inelastic Scattering and using known or calculated structure functions $`F_2`$ and $`R`$ (the “EMC” nuclear effects have to be taken into account).
3) The selection of longitudinal $`\rho ^0`$ will be made by the angular distribution of the decay product. No Rosenbluth separation is envisageable for $`\pi ^0`$ production.
4) The absence of a recoil detector prevents the complete exclusivity of the channel.

Precise simulations of exclusive $`\rho ^0`$ and $`\pi ^0`$ production have been performed  and were already presented at the COMPASS meeting in Munich in 2000. The selection of exclusive events can be summarized as follows:
1) Deep inelastic events are selected by cuts on variables depending on the scattered muon kinematics:
$`2Q^2Q_{max}^2`$ and $`35(20)\nu 170(90)`$ GeV for $`\rho ^0`$
$`1Q^2Q_{max}^2`$ and $`15(10)\nu 170(85)`$ GeV for $`\pi ^0`$
The values outside(inside) the brackets correspond to the beam energy of 190(100) GeV.
2) $`\rho ^0`$ and $`\pi ^0`$ are identified through decays: $`\rho ^0\pi ^+\pi ^{}`$ and $`\pi ^0\gamma \gamma `$. Only two hadrons of opposite charge associated with the vertex defined by the incident and scattered muons are required for $`\rho ^0`$ production and only two photons with the incident and scattered muons are demanded for $`\pi ^0`$ production. For $`\rho ^0`$, hadrons have to be identified as pions. It is then required that each pion decay is emitted in the laboratory at an angle smaller than 180 mrad (the acceptance limit of the forward spectrometer) and its momentum is above 2 GeV. For $`\pi ^0`$ it is demanded that each decay photon has energy above 2 GeV and enters either electromagnetic calorimeter ECAL1 or ECAL2. In addition the separation of two photons at the entrance of a calorimeter should be larger than 4 cm.
3) To isolate at best the exclusive $`\rho ^0`$ events, a cut on the inelasticity I 111$`I=\frac{M_X^2M_p^2}{W^2}`$ where $`W^2=(p+q)^2`$ is the total energy squared in the virtual photon-proton system and $`M_{X^2}=(p+qv)^2`$ is the missing mass squared of the undetected recoiling system ($`p`$, $`q`$, $`v`$ are the four-momenta of the target proton, virtual photon and meson respectively.)  is used. In Fig. 5 the inelasticity distribution is shown for the SMC $`\rho ^0`$ sample  for the events with the invariant mass in the central part of the $`\rho ^0`$ invariant mass peak. For the inelasticity distribution the peak at I=0 is the signal of exclusive $`\rho ^0`$ production. Non-exclusive events, where in addition to detected fast hadrons, slow undetected hadrons are produced, appear at $`I0`$. However, due to the finite resolution, they are not resolved from the exclusive $`\rho ^0`$ peak. For the cut $`0.05I0.05`$ defining the exclusive sample the amount of the residual non-exclusive background for the SMC experiment was up to about 10% at large $`Q^2`$.
4) Finally for the $`\rho ^0`$ channel a cut on the invariant mass of the two pions can be applied in order to reduce the non-resonant background. For the SMC sample the invariant mass distribution after selections, including the cut on inelasticity, is shown in Fig. 5. Although the shape of the mass spectrum varies with $`Q^2`$, a mass cut, i.e. $`0.62m_{\pi ^+\pi ^{}}0.92`$ GeV2, allows a selection of $`\rho ^0`$ events with the relatively low amount of non-resonant background.
It is clear that the good resolution on charged particles associated to $`\rho ^0`$ decay allows the criteria 3 and 4 and thus provides a good signature of the exclusive $`\rho ^0`$ channel. The selection of the $`\pi ^0`$ channel depends strongly on the quality of the electromagnetic calorimeters, but constraints cannot be so nicely determined.

The simulation uses event generators based on a traditional parametrization on NMC data for $`\rho ^0`$ production and on 2 models of GPDs for $`\pi ^0`$ production. Secondary interactions of the decay charged pions, absorption of the decay photons in the target, kinematical smearing based on the experimental resolution, trigger acceptance, acceptance for pions and photons and track reconstruction efficiency are considered. A global $`\rho ^0`$ selection efficiency which takes into account secondary interactions, three tracks efficiency, two pion acceptance, a cut on $`I`$ and $`M_{\pi ^+\pi ^{}}`$ and muon trigger acceptance is evaluated to 0.21(0.04) with the medium $`Q^2`$ trigger (which consists of the “middle trigger" and the “ladder trigger" of the muon trigger hodoscopes) and to 0.36(0.30) with the full $`Q^2`$ range trigger (which includes, in addition, the newly implemented “outer trigger"). The values outside(inside) the brackets correspond to the beam energy of 190(100) GeV. A global $`\pi ^0`$ selection efficiency which takes into account secondary interaction, two photon acceptance and muon trigger acceptance is found close to 0.30.

Total cross sections integrated over the $`Q^2`$ and $`\nu `$ acceptance and expected counting rates (with the medium $`Q^2`$ trigger) for a period of 150 days (1 year) assuming an overall SPS and COMPASS efficiency of 25% are presented in table 1. The two limits for $`\pi ^0`$ production correspond to the two models. About two third of the produced $`\rho ^0`$ are longitudinally polarized.

The background to the reaction $`\mu N\mu \rho ^0N`$ has been studied. It is due to the events with at least two slow undetected particles which are outside the acceptance of the spectrometer or for which tracks are not reconstructed due to inefficiency and which pass all selections for exclusive $`\rho ^0`$ events. A first possible source of this background is $`\rho ^0`$ production with diffractive dissociation of the target $`\mu N\mu \rho ^0N^{}`$ with the subsequent decay of the excited state $`N^{}`$ in $`N+k\pi `$. In the lower part of Fig. 6 a schematic drawing of the inelasticity distribution for these events is presented. Although the smearing effects are not taken into account it is clear that the events from diffractive dissociation of the target into the lowest masses excitations of the nucleon will contribute to the exclusive $`\rho ^0`$ sample in the inelasticity cut. An estimation has been found close to 20%.
A second source of background will originate from the large number of inclusive deep inelastic events. Simulation with a sample of 5000 events only using the generators LEPTO and JETSET has given an upper limit of contamination of also 20%.

## 4 An ideal solution with a completed setup

As mentioned in the previous section it is clear that only a recoil detector which allows the low energy recoiling proton detection will help to select exclusive channels as HEMP or DVCS. The latter reaction is surely the most delicate because one has to select a final state with one muon, one photon and one low energy proton among many competing reactions listed below:
1) Hard Exclusive $`\pi ^0`$ Production $`\mu p\mu p\pi ^0`$ where $`\pi ^0`$ decays in two photons, for which the photon with higher energy imitates a DVCS photon, and the photon with smaller energy is emitted at large angle outside of the acceptance or its energy is below the photon detection threshold.
2) Diffractive dissociation of the proton $`\mu p\mu \gamma N^{}`$ with the subsequent decay of the excited state $`N^{}`$ in $`N+k\pi `$. (The low energy pions are emitted rather isotropically).
3) Inclusive Deep Inelastic Scattering with, in addition to the reconstructed photon, other particles produced outside the acceptance or for which tracks are not reconstructed due to inefficiency.
Moreover one has to take into account a background which includes beam halo tracks with hadronic contamination, beam pile-up, particles from the secondary interactions and external Bremsstrahlung.

A simulation has been realized in order to define the proper geometry of the detector complementing the present COMPASS setup and to analyze the operational conditions. The goal was to maximize the ratio of DVCS events over DIS events for a sample of events with one muon and one photon in the COMPASS spectrometer acceptance plus only one proton of momentum smaller than 750 MeV/c and angle larger than 40 degrees (it is the typical kinematics of a DVCS event at small $`t`$). The simulation relies on the event generator program PYTHIA 6.1  which includes most of the known processes  such as Deeply Inelastic Scattering and Deeply Meson Production. The experimental parameters such as maximum angle and energy threshold for photon detection and maximum angle for charged particle detection could then be tuned. With photon detection extended up to 24 degrees and above an energy threshold of 50 MeV and with charged particle detection up to 40 degrees, one observes that the number of DVCS events as estimated with models is more than an order of magnitude larger than the number of DIS events over the whole useful $`Q^2`$ range (see Fig. 7).

The COMPASS setup will be instrumented with two electromagnetic calorimeters ECAL1 and ECAL2 . They are mainly constituted of lead-glass blocks called GAMS. They are cells of 38.4 $`\times `$ 38.4 $`\times `$ 450 mm3. Typical characteristics of such calorimeter are:
\- energy resolution: $`\sigma P_\gamma /P_\gamma =0.055/\sqrt{P_\gamma }`$ \+ 0.015
\- position resolution: $`\sigma _x=6.0/\sqrt{P_\gamma }+0.5`$ in mm
\- high rate capability: 90% of signal within 50ns gate with no dead time
\- effective light yield: about 1 photoelectron per MeV; hence low energy photons of down 20 MeV can be reconstructed.
The separation of the overlapping electromagnetic showers in the cellular GAMS calorimeter is carefully studied in the Ref. . The result of the study shows that at 10 GeV one can reach a 100% level of the separation efficiency for a minimum distance between 2 photon tracks at the entrance of the calorimeter of D = 4 cm. The last value is slightly shifted to D = 5 cm at 40 GeV.
This excellent performance of the calorimeters will provide a key role in the perfect separation between DVCS events and Hard $`\pi ^0`$ events.

One possible solution to complement the present COMPASS setup is presented in Fig. 8. It consists of one recoil detector described below, an extended calorimetry from 10 to 24 degrees, and a veto for charged forward particles until 40 degrees. This calorimeter has to work in a crowdy environment and in a magnetic fringe field of SM1 and therefore it has to be studied further.

At the present time our studies have focused on the possibility to design and successfully operate a dedicated recoil detector. One goal is to identify and measure the protons momenta between a minimum value and 750 MeV/c. A solution consists in a large time of flight setup between a thin segmented cylindrical layer of scintillator counters, about 3m long and surrounding the 2.5 meter long target, and a thick layer at about 1m distance from the first layer. The thickness of the first layer has to be as small as possible in order to detect protons of minimum momentum. With an hydrogen target of 3cm diameter, target wall thickness of about 3mm of equivalent scintillator and a first layer of 4mm, a minimum momentum of 270 MeV/c is reached. All the counters are read at both sides by photomultiplier counters to determine time and position with very accurate resolutions (300ps and 1.8cm). The consequent resolution in momentum varies from 3 to 10%. The resolution in $`t`$ is twice this value, thus it is very desirable to further study all the parameters which can be improved. Moreover the exclusion of extra particles have to be studied with kinematical fits depending on the experimental resolution and/or with low energy $`\pi ^0`$ detection. This detector has to work in a high rate environment. It has to be as large and hermetic as possible within a reasonable cost. The actual realisation of such performances is under active investigation.

Counting rates, given in table 2, have been estimated assuming 6 months of data taking (1 year) assuming an overall efficiency of 25% and considering the present COMPASS setup where the photon detection is limited to 10 degrees plus a proton detection from 250 to 750 MeV/c ($`|t|0.64`$ GeV2). This statistics allows for the $`\varphi `$ distribution presented in Fig. 4 for $`Q^2`$ = 2, 4 and 6 GeV2 and $`x_{Bj}`$ = 0.05 and 0.10. Studies devoted to the $`t`$ dependence of the cross section can be investigated but for that it is quite worthwhile to try to improve the $`t`$ resolution.

We have tested the concept of this detector using the already existing muon beam and a simplified setup (one sector of scintillators with reduced length). The muon beam was scattered off a 10 cm long polyethylene target, mostly equivalent in radiation length to the foreseen long liquid hydrogen target. We used three scintillators read-out at both sides, a 4mm thick close to the target (A), a 5 cm thick 80 cm away from the target (B) and an extra scintillator (C) to know if particles go through B or are stopped in B. The rates observed in the scintillator close to the target, using the nominal intensity of $`210^8`$ muons per spill, is of the order of 1 MHz (mainly due to Möller electrons). It demonstrates that the background environment is acceptable for the time of flight system.
The result of the time of flight operation (see Fig. 9(a)) shows a clear proton signal. With the knowledge of the $`\beta `$ velocity and the energy lost in B for stopped particles one can reconstruct their masses. It is done in Fig. 9(b) where one can see pions, protons and deuterons for raw data, corrected data and the target out contribution which is about two order of magnitude smaller.

The position resolution obtained on A and B and the time of flight resolution are better than 2 cm and 300 ps respectively. Extension to long (3 meter) and thin scintillators have to be studied carefully and technology has to be improved to achieve still better resolution. An efficiency study of such a recoil detector is being performed.

## 5 CONCLUSIONS

This study takes advantage of the high energy of the muon beam available at CERN which provides a large $`Q^2`$ and $`x_{Bj}`$ range and encourages us for the following roadmap. Hard Exclusive Meson Production have to be undertaken as soon as possible with the present setup. A large number of $`\rho ^0`$ events (a few 10K) can be produced in 1 year. The $`\rho ^0`$ channel which decays in $`\pi ^+\pi ^{}`$ is the easiest channel to isolate, the $`\pi ^0`$ channel is more difficult but very important to test the calorimetry performances. A complete experiment with both Hard Exclusive Meson Production with a large set of mesons and Deeply Virtual Compton Scattering has to be envisaged in a next step with a completed COMPASS setup. For this purpose one needs a “long" hydrogen target, a recoil detector and an extension of the calorimetry at larger angles.

COMPASS is the unique place which provides $`\mu ^+`$ and $`\mu ^{}`$ of 100 GeV in order to study carefully two scales of observation $`x_{Bj}=0.05\pm 0.02`$ and $`x_{Bj}=0.10\pm 0.03`$ on a large domain of $`Q^2`$ from 2 to 7 GeV2 and to measure the azimuthal distribution of the Beam Charge Asymmetry which seems very promising to test the geometrical interpretation of GPDs.

## ACKNOWLEDGMENTS

We acknowledge useful discussions with Dietrich von Harrach, Fritz Klein, Alain Magnon and all members of the Saclay COMPASS group.


format: text

text:
# Synchronization in coupled map lattices as an interface depinning

## Abstract

We study an SOS model whose dynamics is inspired by recent studies of the synchronization transition in coupled map lattices (CML). The synchronization of CML is thus related with a depinning of interface from a binding wall. Critical behaviour of our SOS model depends on a specific form of binding (i.e., transition rates of the dynamics). For an exponentially decaying binding the depinning belongs to the directed percolation universality class. Other types of depinning, including the one with a line of critical points, are observed for a power-law binding.

Recently synchronization of chaotic systems attracted a lot of interest FUJISAKA . To a large extent, this interest is motivated by numerous experimental realizations of this phenomena including lasers, electronic circuits, or chemical reactions BOCCAL . Synchronization acquires additional features in spatially extended systems, where it can be regarded as a certain nonequilibrium phase transition. There are increasing efforts to understand the properties of this transition. Relatively well understood is the synchronization transition (ST) in certain cellular automata. Since a synchronized state can be regarded as an absorbing state, for cellular automata, which are discrete systems, the phase transition, as expected, belongs to the Directed Percolation (DP) universality class GRASS . However, continuous systems, as e.g., coupled map lattices (CML) KANEKO , need infinite time to reach a synchronized state and the relation with DP does not seem to hold. Indeed, Pikovsky and Kurths argued PIKKURTHS  that for continuous systems this transition should belong to the so-called bounded Kardar-Parisi-Zheng (BKPZ) universality class TU . Recently, precise numerical calculations confirmed their predictions but only for CML with some continuous maps AHLERS . Surprisingly, the ST for discontinuous LIVI ; AHLERS , or continuous but sufficiently steep LIPDROZ  maps was found to belong to the DP universality class. It would be desirable to understand the critical behaviour of the ST in CML and the present paper might be a step in this direction.

Let us briefly describe the setup which is used to study synchronization in CML AHLERS . In the simplest, one-dimensional case, one takes a single-chain CML of size $`L`$ that is composed of $`L`$ diffusively coupled local maps $`f(u_i)`$ KANEKO . The maps are chaotic and act on continuous site variables $`u_i`$ ($`i=1,\mathrm{},L`$), that are typically bounded ($`0<u_i<1`$). Then, one couples such a spatio-temporally chaotic system with its identical copy, which initially has a different set of site variables. It turns out that the evolution of such a coupled system depends on the coupling strength. For weak coupling the two CML’s are desynchronized and essentially independent. However, for sufficiently strong coupling the system gets synchronized and approaches a state where corresponding pairs of site variables in both copies take the same values. To quantify synchronization one can introduce the synchronization error $`w(i,t)=|u_1(i,t)u_2(i,t)|`$ where a lower index denotes a copy of CML, and its spatial average $`w(t)=\frac{1}{L}_{i=1}^Lw(i,t)`$. To relate this problem with BKPZ one argues that in the continuous limit and close to the synchronized state, the evolution equation for $`w(i,t)`$ is given as a Langevin equation with multiplicative Gaussian noise which then, using the Hopf-Cole transformation $`h=\mathrm{ln}(w)`$, is transformed into BKPZ. Let us notice that CML’s are deterministic systems and the noise has only effective meaning, mimicking their chaotic behaviour. The problem of a relation of such deterministic systems with stochastic counterparts is very interesting and recently is drawing some attention EGOLF . In the above representation the desynchronized phase in CML ($`w(t)>0`$) corresponds to the interface pinned relatively close to the wall ($`h_i<\mathrm{}`$). In the synchronized phase ($`w(t)0`$) the interface depins and drifts away ($`h_i\mathrm{}`$). Perfectly synchronized state ($`w(t)=0`$) is reached only after infinitely long time. The above analysis requires the differentiability of the local map and thus is not applicable to discontinuous maps. However, numerical results show that the relation with BKPZ breaks down also for continuous but sufficiently steep maps LIPDROZ . From a theoretical point of view it would be desirable to understand why such properties of the local map affect the nature of the ST and move into the DP universality class. Let us notice, that if the coupled CML system enters the synchronized state, it will remain in this state forever. Such a state can be thus considered as an absorbing state of the dynamics, although it cannot be reached in any finite time. Well-developed techniques are available to study phase transitions in models with absorbing states HAYE .

It is clear that the problem of synchronization in extended systems is related with a number of very interesting problems in nonequilibrium statistical mechanics such as the KPZ model, nonequilibrium wetting PAST  or directed percolation. Recently, some arguments were given that a notoriously difficult particle system, the so-called PCPD model, might be also related with these modelsHIN . Deeper understanding of these problems and their mutual relations would be certainly desirable.

In the present paper we introduce an interfacial, discrete (SOS) model which is inspired by the dynamics of CML with discontinuous maps. The synchronization of CML is thus related with a depinning of interface from a binding wall. Numerical calculations show that the universality class of the depinning transition in our model depends on the choice of binding, which enters the dynamics through certain transition rates. For an exponentially decaying binding, the depinning belongs to the DP universality class. In this case the overall behaviour of the model is very similar to Sneppen’s model of interface propagation in a random environment that is driven by a certain extremal dynamics SNEPPEN ; OLAMI . For a power-law decaying binding the depinning transition is characterized by a different set of critical exponents. In the case of a rapid decay these exponents are very close to those obtained for the bosonic version of PCPD model KOCKEL . In the case of a slow decay, in the entire unbounded phase the interface remains critical, that in itself is an interesting property of a nonequilibrium system.

In our model discrete site variables $`h_i=1,2,\mathrm{}`$ are defined on a one-dimensional lattice of size $`L`$ ($`i=1,\mathrm{},L`$) with periodic boundary conditions ($`h_{L+1}=h_1`$). In an elementary update we select randomly a site $`i`$ and change the variable $`h_i`$, and possibly neighbouring ones, according to the following rule: (i) with probability $`p(h_i)`$ one sets $`h_i=h_{i+1}=h_{i1}=1`$; (ii) with probability $`1p(h_i)`$ the site variable $`h_i`$ increases by unity ($`h_ih_i+1`$). During a unit of time $`L`$ elementary updates are performed. To complete the definition we have to specify the function $`p(h)`$. To allow a drift toward $`h=\mathrm{}`$ the function $`p(h)`$ must decay to zero for $`h\mathrm{}`$. Numerical results that we present below are obtained for two cases: $`p(h)=\mathrm{e}^{\gamma h}`$ (model I) and $`p(h)=a(h+1)^\gamma `$ (model II), where $`\gamma >0`$ and $`a>0`$ are control parameters of the model.

To make a link with synchronization in CML’s the following remarks are in order. Numerical calculations for CML’s with discontinuous maps show that the Lyapunov exponent that governs the evolution of the synchronization error $`w(i,t)`$ is negative in the vicinity of the transition AHLERS . Approximately, the evolution of $`w(i,t)`$ is thus made of consecutive contractions ($`w(i,t)cw(i,t)`$ and $`c<1`$) that, due to discontinuity of the map, are from time to time interrupted by discontinuous changes that might substantially increase the value of $`w(i,t)`$. To notice a link between the dynamics of synchronization error in CML and our model we introduce new variables $`w_i=\mathrm{e}^{h_i}`$ (let us notice a similarity to the inverse Hopf-Cole transformation). Indeed, the increase of $`h_i`$ by unity according to the rule (ii) decreases $`w_i`$ by a factor e that corresponds to the contraction of $`w(i,t)`$. The first rule mimics the discontinuous jumps of $`w(i,t)`$. Since local maps in CML’s are coupled, a jump at a site $`i`$ also affects its neighbours.

Monte Carlo simulations of our model are similar to those of other models with absorbing states HAYE . For some details related to the fact that the model needs an infinite time to reach an absorbing state see e.g., LIPDROZ . We observed that for sufficiently large $`\gamma `$ the interface depins from the $`h=1`$ wall and drifts away. For smaller $`\gamma `$ the model remains in the active phase with the interface relatively close to the wall. To examine the nature of the phase transition we introduced $`w(t)=\frac{1}{L}_{i=1}^Lw(i,t)`$ that in the steady state is denoted as $`w`$ (in the following we refer to this quantity as activity). Of course in the absorbing phase $`w=0`$ and in the active phase $`w>0`$. Upon approaching the critical point $`w`$ typically exhibits a power-law decay $`w(\gamma \gamma _c)^\beta `$, with a characteristic exponent $`\beta `$ and for the critical point located at $`\gamma =\gamma _c`$. Moreover, we studied the time dependence of $`w(t)`$. One expects that at criticality this quantity has a power law decay $`w(t)t^\mathrm{\Theta }`$, where $`\mathrm{\Theta }`$ is another characteristic exponent. We also used the so-called spreading technique SPREADING . First, we set $`h_i=\mathrm{}`$ for every but one site ($`i_0`$) that was set to unity. Then we monitored the subsequent evolution of the model (actually, of interest are only sites with positive $`w(i,t)`$ i.e., with finite $`h_i`$) measuring the average activity of the system $`w(t)`$, the survival probability $`P(t)`$, and the averaged spread square $`R^2(t)=\frac{1}{w(t)}_iw(i,t)(ii_0)^2`$. One expects that at criticality: $`w(t)t^\eta `$, $`P(t)t^\delta `$, and $`R^2(t)t^z`$, that at the same time defines the critical exponents $`\eta `$, $`\delta `$, and $`z`$.

First, let us describe the results obtained for model I ($`p(h)=\mathrm{e}^{\gamma h}`$). From the steady state measurements of $`w`$ (Fig. 1) we estimate $`\gamma _c=0.791(1)`$, and $`\beta =0.28(1)`$. Such a location of the critical point is confirmed from time-dependent simulations (Fig. 2a) and the spreading method (Fig. 2b). From these data we estimate $`\mathrm{\Theta }=0.160(2)`$, $`\eta =0.317(5)`$, $`\delta =0.16(1)`$, and $`z=1.26(1)`$. The results for $`P(t)`$ and $`R^2(t)`$ are not presented. Obtained values of critical exponents clearly show that model I belongs to the DP universality class for which HAYE : $`\beta =0.2765`$, $`\eta =0.3137`$, $`\mathrm{\Theta }=\delta =0.1595`$ and $`z=1.265`$.

To have a more complete insight into the critical behaviour of our model we calculate the average interfacial width $`W(t)=[_{i=1}^L(h(i,t)h(i,t))]^2^{1/2}`$. At criticality this quantity typically behaves as $`W(t)t^\beta ^{}`$ where $`\beta ^{}`$ is the growth exponent. For model I simulations show almost linear increase of $`W(t)`$ with time (inset in Fig. 3) and we estimate $`\beta ^{}=1.0(1)`$. Such a value of the growth exponent indicates very strong fluctuations, which is confirmed through a visual inspection of the interface profile (Fig. 3).

In principle, from direct calculations for CML’s we can obtain an interface profile as a logarithm of the synchronization error. Our preliminary calculations for CML model examined by Ahlers and Pikovsky AHLERS  show that large fluctuations, similar to those in model I are clearly seen. However, to quantify these fluctuations and calculate e.g., the width $`W(t)`$, simulations on a longer time scale are required. Since for CML dynamics is defined in terms of $`w`$’s rather than $`h`$’s, simulations for longer time, which must probe states of extremely small synchronization error, severely suffer from the finite accuracy of numerical computations. On the other hand our SOS model easily allows us to examine such a long-time regime.

It is already known that some SOS-like models belong to the DP universality class. In particular, Alon et al. DPWETT  introduced a model of interface roughening with edge evaporation. In this model the growth of the interface width is logarithmic in time and is thus much different than in our model. As far as the critical behaviour is concerned, our model seems to be more closely related with the Sneppen’s model of interface spreading in environment with quenched disorder SNEPPEN . His model is driven by extremal dynamics and has the growth exponent $`\beta ^{}`$ close to unity. Later it was established that this model is actually equivalent to the DP model sitting exactly at the critical point OLAMI . Important ingredients of Sneppen’s model are quenched disorder and extremal dynamics. The fact that our model, that misses these features, exhibits essentially the same behaviour is in our opinion quite interesting and worth further examination.

Using the same procedure we studied the model II ($`P(h)=a(h+1)^\gamma `$). First we kept $`a=1`$ fixed and varied only the parameter $`\gamma `$. Some of our numerical results are shown in Fig. 4-Fig. 5. Using these data we estimate: $`\gamma _c=1.74925(5)`$, $`\beta =0.36(1)`$, $`\mathrm{\Theta }=0.185(3)`$, $`\eta =0.03(1)`$, $`\delta =0.445(5)`$, and $`z=1.19(1)`$. All exponents considerably differ from DP exponents. To support these estimations, let us notice that the hyperscaling relation $`\frac{z}{2}=\mathrm{\Theta }+\delta +\eta `$ is satisfied by the above values. It is rather surprising for us to observe that our values of exponents $`\beta `$, $`\mathrm{\Theta }`$ and $`z`$ are in a very good agreement with recent estimation for the so-called bosonic version of a Pair Contact Process with Diffusion (PCPD) KOCKEL . It would be interesting to examine whether this is only a numerical coincidence or there is a deeper relation between these two problems.

We also studied model II for fixed $`\gamma `$ and varying the amplitude $`a`$. Monte Carlo simulations show that for $`\gamma =2`$ the critical behaviour of the model seems to be the same as for $`\gamma =\gamma _c=1.74925`$ and $`a=1`$. It suggests that for a certain range of $`\gamma `$ the depinning transition belongs to the same universality class. However, a different behaviour was observed for $`\gamma =1`$. Indeed, in this case our simulations suggest (Fig. 6) that the depinning transition is characterized by the exponents $`\beta 3.0(3)`$ and $`\mathrm{\Theta }=0.62(5)`$. What is also interesting, the entire unbounded phase $`a<a_c0.07`$ seems to be critical with a power-law decaying order parameter $`wt^x`$, where $`x0.8(5)`$, and $`x`$ might slightly vary with $`a`$. Such a novel behaviour is also worth further studies.

One of the future problems would be to check whether there are some other types of critical behaviour for this kind of SOS models. Clearly, an important ingredient that determines the critical behaviour is the form of the function $`p(h)`$. It would be interesting to explore some other functions (e.g., $`\mathrm{e}^{\gamma h^2}`$ or $`[\mathrm{ln}(h)]^\gamma `$) for a possibly new behaviour. Although we were not able to recover the BKPZ universality class within our approach, there is still a possibility that for a certain choice of $`p(h)`$ such a behaviour might appear. It would be also interesting to check whether a new critical behaviour found in model II has a counterpart in a synchronization transition in CML.

This work was partially supported by the Swiss National Science Foundation and the project OFES 00-0578 ”COSYC OF SENS”.


text:
# Stable vortex and dipole vector solitons in a saturable nonlinear medium

## I Introduction

Spatial solitons have been a subject of many studies since their first theoretical prediction . The previous research on spatial solitons was driven by their promising applications in all-optical devices in which the light guides and steers the light itself . Early works studied optical materials with the Kerr (cubic) nonlinearity which exhibit stable fundamental (single-hump) solitons in one spatial dimension  and collapse in two and three spatial dimensions . Later works focused on optical materials with saturable nonlinear response such as photorefractive crystals (see  and references therein). The nonlinearity saturation suppresses the collapse of fundamental solitons in two and three dimensions , which opens the door for their experimental observation in multi-dimensional optical beams. The instability of higher-order (multi-hump) solitons is not suppressed by the nonlinearity saturation however . Internal modes of fundamental solitons in saturable optical materials have also been reported . These modes are responsible for long-lived shape oscillations.

Recently, incoherent coupling of spatial solitons in photorefractive crystals was proposed and experimentally demonstrated . The mutual trapping of incoherent optical beams leads to many novel spatial solitons such as vortex and dipole vector solitons . The incoherently coupled spatial solitons are described by a system of coupled nonlinear Schrödinger (NLS) equations. A similar system of equations also describes temporal solitons in birefringent optical fibers and wavelength-division-multiplexed systems . Additionally, vortex vector solitons are known in the Bose-Einstein condensation guided by a magnetic trap .

Vortex and dipole vector solitons in saturable optical materials are interesting for both physical and mathematical reasons. Physically, these spatial solitons are novel nonlinear objects. They bifurcate from a coupled state where a fundamental soliton in one component guides a small higher-order mode in the other component. Far from the bifurcation threshold, both components strongly trap each other and form a fully coupled vector soliton. Mathematically, existence and stability of the vortex and dipole vector solitons are challenging problems due to their complexity. The existence of vortex and dipole solitons was established numerically  and with a heuristic variational method . However, analytical expressions for radial-angular dependences of vortex and dipole solitons have not been found. On the stability of vortex solitons, the numerical results of  suggest that vortex solitons with small vortex components are stable and observable for large propagation distances, while vortex solitons with large vortex components are unstable. Numerical results of  show however that all vortex solitons are linearly unstable, and the instability leads to the breakup of vortex solitons into rotating dipole solitons. The discrepancy between  and  raises an open question: are vortex vector solitons with small vortex components really stable or not? On the other hand, dipole solitons were found numerically to be always stable in .

In this paper, we clarify the issues of existence and stability of vortex and dipole vector solitons in a saturable nonlinear medium such as photorefractive crystals. First, we address the existence and uniqueness of vortex and dipole solitons with the perturbation technique . We derive perturbation series expansions for the vortex and dipole solitons near the bifurcation point where the vortex and dipole components are small. The analytical formulas for these solitons are in excellent agreement with our numerical results. We also prove that those vector solitons are unique up to phase-, translation-, and rotation-invariances. Next, we study the linear stability of vortex and dipole solitons with both the spectral analysis and numerical methods. We show that dipole solitons are linearly stable in the entire existence domain, while the vortex solitons with large vortex components are linearly unstable, in agreement with . However, we prove that vortex solitons with small vortex components are linearly stable, confirming the results of , not . Lastly, we study the nonlinear evolution of linearly unstable vortex solitons. We show that an unstable vortex soliton breaks up into a rotating dipole soliton only when the vortex component is below a certain threshold. Above this threshold, an unstable vortex soliton breaks up into two fundamental vector solitons instead.

## II Existence and uniqueness of vortex and dipole solitons

The mathematical model for two incoherently-coupled laser beams in a photorefractive crystal is well known (see, e.g., ). After variable rescalings, the model can be written as a system of coupled equations:

$$i\frac{E_1}{z}+\mathrm{\Delta }E_1+\frac{E_1(|E_1|^2+|E_2|^2)}{1+s(|E_1|^2+|E_2|^2)}=0,$$
(1)
$$i\frac{E_2}{z}+\mathrm{\Delta }E_2+\frac{E_2(|E_1|^2+|E_2|^2)}{1+s(|E_1|^2+|E_2|^2)}=0,$$
(2)
where $`\mathrm{\Delta }`$ is the two-dimensional Laplacian, and $`s`$ is the saturation parameter.

Vector solitons in this model take the form:

$$E_1=u(x,y)e^{iz},E_2=w(x,y)e^{i\lambda z},$$
(3)
where the frequency of the $`E_1`$ wave has been normalized to one, and the frequency of the $`E_2`$ wave is $`\lambda `$. The amplitude functions $`u(x,y)`$ and $`w(x,y)`$ satisfy the nonlinear boundary-value problem with zero boundary conditions on the $`(x,y)`$ plane:

$$\mathrm{\Delta }uu+\frac{u(|u|^2+|w|^2)}{1+s(|u|^2+|w|^2)}=0,$$
(4)
$$\mathrm{\Delta }w\lambda w+\frac{w(|u|^2+|w|^2)}{1+s(|u|^2+|w|^2)}=0.$$
(5)
The system (4)–(5) may have several types of vector solitons localized in two dimensions. The fundamental vector solitons take the form:

$$u=c_u\mathrm{\Phi }(r),w=c_w\mathrm{\Phi }(r),$$
(6)
where $`\mathrm{\Phi }(r)`$ is a real-valued, single-hump function, $`r=\sqrt{x^2+y^2}`$, and $`c_u`$ and $`c_w`$ are arbitrary complex parameters constrained by the relation $`|c_u|^2+|c_w|^2=1`$. These solitons exists at $`\lambda =1`$ where the system (4)–(5) reduces to a scalar equation for $`\mathrm{\Phi }(r)`$, see Eq. (10) below. The vortex vector solitons take a general form:

$$u=\mathrm{\Phi }_u(r)e^{in\theta },w=\mathrm{\Phi }_w(r)e^{im\theta },$$
(7)
where $`n`$ and $`m`$ are topological charges of vortices in the $`u`$ and $`w`$ components, $`\mathrm{\Phi }_u(r)`$ and $`\mathrm{\Phi }_w(r)`$ are real-valued functions, and $`(r,\theta )`$ are the polar coordinates on the $`(x,y)`$ plane. The simplest vortex soliton has $`n=0`$ and $`m=1`$ . The multi-pole vector solitons take yet another general form:

$$u=U(r,\theta ),w=W(r,\theta ),$$
(8)
where $`U(r,\theta )`$ and $`W(r,\theta )`$ are real-valued functions and may have multi-hump profiles on the $`(x,y)`$ plane. The multi-pole solitons with a single hump for $`u(x,y)`$ and multiple humps for $`w(x,y)`$ were approximated in the variational approach  by the ansatz:

$$u=U(r),w=W(r)\mathrm{cos}m\theta ,$$
(9)
where the number of humps in the $`w`$ component is $`2m`$. The simplest multi-pole soliton is a dipole soliton which has $`m=1`$.

In this paper, we study the simplest vortex and multi-pole vector solitons. We will refer to them simply as the vortex soliton and dipole soliton hereafter. The existence and uniqueness of the vortex and dipole solitons can be studied by perturbation methods. The perturbation series expansions are derived in the neighborhood of the bifurcation value $`\lambda =\lambda _0(s)`$, where the $`w`$ component is small. With the perturbation arguments, we show that the bifurcation value $`\lambda _0(s)`$ is the same for both branches of vortex and dipole solitons, and these solitons bifurcate uniquely from $`\lambda =\lambda _0(s)`$.

Setting $`w=0`$ in (4), we find the nonlinear boundary-value problem for the scalar soliton $`u=u_0(r)`$:

$$u_0^{\prime \prime }+\frac{1}{r}u_0^{}u_0+\frac{u_0^3}{1+su_0^2}=0,$$
(10)
where $`u_0(r)`$ is a real-valued function. We take $`u_0(r)`$ to be the fundamental soliton, i.e., $`u_0(r)>0`$ for finite $`r0`$. If $`w`$ is small in (5), we get a linear eigenvalue problem for the first-order correction $`w=w_1(x,y)`$:

$$\mathrm{\Delta }w_1\lambda w_1+\frac{u_0^2}{1+su_0^2}w_1=0,$$
(11)
where $`w_1(x,y)`$ is a complex-valued function in general and $`\lambda `$ is the eigenvalue. The linear equation (11) supports localized solutions of the form $`\varphi (r)e^{\pm im\theta }`$ at a set of discrete values of $`\lambda `$. Since we study the simplest vortex and dipole solitons, we set $`m=1`$, and require $`\varphi (r)`$ to be a non-negative function for $`r0`$. The corresponding eigenvalue $`\lambda _0(s)`$ and eigenfunction $`\varphi (r)`$ satisfy the following reduced equation:

$$\varphi ^{\prime \prime }+\frac{1}{r}\varphi ^{}\left(\lambda _0+\frac{1}{r^2}\right)\varphi +\frac{u_0^2}{1+su_0^2}\varphi =0.$$
(12)
The eigenvalue $`\lambda `$ is unique once $`s`$ is fixed. We normalize the eigenfunction $`\varphi (r)`$ such that it has a maximum value one. Numerically, we compute $`\lambda _0(s)`$ and $`\varphi (r)`$ by the shooting method. Fig. 1a shows the dependence of $`\lambda _0`$ versus saturation parameter $`s`$. Fig. 1b shows the scalar soliton $`u_0(r)`$ and the normalized eigenfunction $`\varphi (r)`$ at $`s=0.5`$, where $`\lambda _0=0.2622`$.

When $`\varphi (r)`$ and $`\lambda _0(s)`$ are known, a general solution for $`w_1(r,\theta )`$ can be written in the form:

$$w_1(r,\theta )=\varphi (r)(\mathrm{cos}\theta +ip\mathrm{sin}\theta ),$$
(13)
where $`p`$ is an arbitrary real parameter. In this general solution, we have removed arbitrary rotations and translations on the $`(x,y)`$ plane as well as an arbitrary phase shift in the $`w`$ component. We note that the solution (13) is identical to the variational ansatz in .

Below, we use the perturbative method and show that there are only two continuations of the solution (13): for $`p=\pm 1`$ and $`p=0`$. When $`p=\pm 1`$, the perturbation series expansion recovers the vortex soliton (7) with $`n=0`$ and $`m=\pm 1`$. When $`p=0`$, the perturbation expansion recovers the dipole soliton (8) with a single hump for $`u`$ and a double hump for $`w`$. For given values of $`s`$ and $`\lambda `$, the two solutions are unique up to phase-, translation- and rotation-invariances. At other values of $`p`$, the solution with $`w`$’s leading-order term as (13) can not exist.

The perturbation series expansions for vector solitons in system (4) and (5) take the form:

$$u=u_0(r)+ϵ^2u_2(r,\theta )+ϵ^4u_4(r,\theta )+\mathrm{O}(ϵ^6),$$
(14)
$$w=ϵw_1(r,\theta )+ϵ^3w_3(r,\theta )+ϵ^5w_5(r,\theta )+\mathrm{O}(ϵ^7),$$
(15)
and

$$\lambda =\lambda _0(s)+ϵ^2\lambda _2(s)+ϵ^4\lambda _4(s)+\mathrm{O}(ϵ^6),$$
(16)
where $`ϵ`$ is a small parameter, $`u_0(r)`$ is the scalar fundamental soliton solving (10), $`w_1(r,\theta )`$ is the first-order correction in the form (13), and the cut-off frequency $`\lambda _0`$ is the eigenvalue of (12). The objective of the perturbation analysis is to uniquely determine the coefficients $`\lambda _2,\lambda _4,\mathrm{}`$ as well as expressions for functions $`u_2,u_4,w_3,w_5`$ and so on. Once the coefficients $`\lambda _0,\lambda _2,\mathrm{}`$ have been obtained, we can compute $`ϵ`$ from $`\lambda `$ in the expansion (16). Once $`ϵ`$ is found, together with functions $`u_0,u_2,w_1,w_3,\mathrm{}`$, we can approximate the vector soliton by the expansions (14) and (15). Below, we will carry out the perturbative calculations to order $`ϵ^3`$.

Substituting the perturbation series (14), (15) and (16) into the original equations (4) and (5), at order $`ϵ^2`$, we get an inhomogeneous equation for $`u_2`$:

$$\mathrm{\Delta }u_2u_2+\frac{u_0^2(2+su_0^2)}{(1+su_0^2)^2}u_2+\frac{u_0^2}{(1+su_0^2)^2}\overline{u}_2=\frac{u_0|w_1|^2}{(1+su_0^2)^2}.$$
(17)
where $`\overline{u}`$ is the complex conjugate of $`u`$. The linearized operator in the left-hand-side of (17) has a non-empty null space spanned by three linearly independent localized eigenfunctions:

$$u_{2h}^{(1)}=iu_0(r),u_{2h}^{(2)}=u_0^{}(r)\mathrm{cos}\theta ,u_{2h}^{(3)}=u_0^{}(r)\mathrm{sin}\theta .$$
(18)
These eigenfunctions correspond to the phase and translational invariances of solitons in the scalar $`u`$ equation. The right-hand-side of (17) is orthogonal to $`u_{2h}^{(1)}`$ because it is real-valued. It is also orthogonal to $`u_{2h}^{(2)}`$ and $`u_{2h}^{(3)}`$ because it has different angular dependence of $`1`$ and $`\mathrm{cos}2\theta `$ (rather than $`\mathrm{cos}\theta `$ and $`\mathrm{sin}\theta `$). Therefore, up to phase-, translation-, and rotation-shifts, a localized solution to (17) is constructed uniquely in the form:

$$u_2=\frac{1}{2}(1+p^2)u_{20}(r)+\frac{1}{2}(1p^2)u_{22}(r)\mathrm{cos}2\theta ,$$
(19)
where functions $`u_{20}(r)`$ and $`u_{22}(r)`$ satisfy the equations

$$u_{20}^{\prime \prime }+\frac{1}{r}u_{20}^{}u_{20}+\frac{u_0^2(3+su_0^2)}{(1+su_0^2)^2}u_{20}=\frac{u_0\varphi ^2}{(1+su_0^2)^2},$$
(20)
$$u_{22}^{\prime \prime }+\frac{1}{r}u_{22}^{}\left(1+\frac{4}{r^2}\right)u_{22}+\frac{u_0^2(3+su_0^2)}{(1+su_0^2)^2}u_{22}=\frac{u_0\varphi ^2}{(1+su_0^2)^2}.$$
(21)
We do not know exact analytical expressions for $`u_{20}(r)`$ and $`u_{22}(r)`$ but can compute them numerically.

At order $`ϵ^3`$, we get the equation for $`w_3`$ as

$$\mathrm{\Delta }w_3\lambda _0w_3+\frac{u_0^2}{1+su_0^2}w_3=\left\{\lambda _2\frac{|w_1|^2+2u_0u_2}{(1+su_0^2)^2}\right\}w_1.$$
(22)
We denote

$$h_1(r)=\frac{\varphi ^2+2u_0u_{20}}{2(1+su_0^2)^2},h_2(r)=\frac{\varphi ^2+2u_0u_{22}}{2(1+su_0^2)^2},$$
(23)
and rewrite Eq. (22) in the form:

$`\mathrm{\Delta }w_3\lambda _0w_3+{\displaystyle \frac{u_0^2}{1+su_0^2}}w_3=\left[\lambda _2(1+p^2)h_1{\displaystyle \frac{1}{2}}(1p^2)h_2\right]\varphi \mathrm{cos}\theta `$ (24)
$`+ip\left[\lambda _2(1+p^2)h_1+{\displaystyle \frac{1}{2}}(1p^2)h_2\right]\varphi \mathrm{sin}\theta {\displaystyle \frac{1}{2}}(1p^2)h_2\varphi \mathrm{cos}3\theta {\displaystyle \frac{1}{2}}ip(1p^2)h_2\varphi \mathrm{sin}3\theta .`$ (25)
The homogeneous part in (24) supports two linearly independent localized solutions $`\varphi (r)\mathrm{cos}\theta `$ and $`\varphi (r)\mathrm{sin}\theta `$. As a result, a localized solution of the non-homogeneous equation (24) exists if and only if the following solvability conditions are satisfied:

$$_0^{\mathrm{}}r\varphi ^2\left\{\lambda _2(1+p^2)h_1(r)\frac{1}{2}(1p^2)h_2(r)\right\}𝑑r=0,$$
(26)
$$p_0^{\mathrm{}}r\varphi ^2\left\{\lambda _2(1+p^2)h_1(r)+\frac{1}{2}(1p^2)h_2(r)\right\}𝑑r=0.$$
(27)
We will show below that these solvability conditions define only two perturbation series solutions for vector solitons. These solutions correspond to the choice $`p=\pm 1`$ or $`p=0`$, which produce vortex and dipole solitons respectively.

### A Vortex solitons

If $`p0`$, we eliminate parameter $`\lambda _2`$ from the system (26)–(27) and find the solvability condition in the form:

$$(1p^2)_0^{\mathrm{}}r\varphi ^2h_2(r)𝑑r=0.$$
(28)
The integral in Eq. (28) only depends on the parameter $`s`$, not on $`p`$. We have checked numerically that this integral never vanishes for any $`s`$. Thus, the condition (28) is satisfied only when $`p=\pm 1`$. In this case, it follows from (13), (19) and (24) that $`w_1=\varphi (r)e^{\pm i\theta }`$, $`u_2=u_{20}(r)`$, and $`w_3=f(r)e^{\pm i\theta }`$. We can continue the perturbation series expansions (14)–(16) to higher orders and find that all $`u_{2n}(n0)`$ corrections are only functions of $`r`$, and all $`w_{2n+1}(n0)`$ corrections have the form $`g(r)e^{\pm i\theta }`$. Thus, the perturbation series solution gives a vortex vector soliton (7) with $`n=0`$ and $`m=\pm 1`$.

When $`p=\pm 1`$, we find from Eq. (26) that the coefficient $`\lambda _2`$ is

$$\lambda _2=\lambda _{2v}(s)\frac{2_0^{\mathrm{}}r\varphi ^2h_1𝑑r}{_0^{\mathrm{}}r\varphi ^2𝑑r}.$$
(29)
The functional dependence of $`\lambda _{2v}`$ versus $`s`$ is computed from this formula and plotted in Fig. 1a (dashed line) alongside the cut-off frequency $`\lambda _0(s)`$. Since the (non-negative) function $`\varphi (r)`$ is normalized to have maximum one, the perturbation parameter $`ϵ`$ determines the amplitude (maximum) of the vortex component $`w`$ with error at the order of $`ϵ^3`$, see (13) and (15). The dependence of $`ϵ`$ versus $`\lambda `$ and $`s`$ can be obtained from the perturbation series (16) with error of order $`ϵ^2`$:

$$ϵ=\sqrt{\frac{\lambda \lambda _0(s)}{\lambda _{2v}(s)}}.$$
(30)
We compare the analytical formula (30) with numerical results for $`s=0.5`$, where $`\lambda _0=0.2622`$ and $`\lambda _{2v}=0.1010`$. A dashed line in Fig. 2(a) shows the amplitude of the vortex component $`w`$ computed from (30). Numerically, vortex solitons are computed from the original system (4)–(5) by the shooting method. The amplitudes of the $`u`$ and $`w`$ components are also shown in Fig. 2. In Fig. 2(b), a profile of $`u(x,y)`$ and $`w(x,y)`$ components for $`s=0.5`$ and $`\lambda =0.4`$ is shown. We can see from Fig. 2a that the agreement between the analytical predictions and numerical values on the amplitudes of the $`w`$-component is good over a wide range of $`\lambda `$ values. Also, the numerically-obtained amplitude of the $`u`$-component depends linearly on $`\lambda `$, which is in agreement with the perturbation series (14) up to error of order $`ϵ^4`$, since $`ϵ^2(\lambda \lambda _0)`$.

### B Dipole solitons

If $`p=0`$, the condition (27) is satisfied, while the condition (26) gives the correction term $`\lambda _2`$ in the form:

$$\lambda _2=\lambda _{2d}(s)\frac{_0^{\mathrm{}}r\varphi ^2(2h_1+h_2)𝑑r}{2_0^{\mathrm{}}r\varphi ^2𝑑r}.$$
(31)
The functional dependence of $`\lambda _{2d}`$ versus $`s`$ is numerically computed and plotted in Fig. 1(a) (solid line). When $`\lambda _2`$ is given by (31), a localized solution $`w_3(r,\theta )`$ of Eq. (24) exists, and this solution is real-valued. We can further show that the perturbation series expansions (14)–(16) can be successfully continued to higher orders of $`ϵ`$, and a dipole-soliton solution (8) can be obtained. This solution is real-valued and has the symmetries

$`u(x,y)=u(x,y),u(x,y)=u(x,y),w(x,y)=w(x,y),w(x,y)=w(x,y).`$ (32)
Similar to the vortex soliton case, the perturbation parameter $`ϵ`$ here gives the amplitude (maximum) of the dipole component $`w`$ with accuracy of O($`ϵ^3`$). The formula for $`ϵ`$ is still (30), but the $`\lambda _2`$ value is now given by Eq. (31) instead of (29). The comparison between the analytical result (30)–(31) and numerical results for dipole solitons is shown in Fig. 3a for $`s=0.5`$. In this case, $`\lambda _0=0.2622`$ and $`\lambda _{2d}=0.1174`$. We see again that the agreement between numerical and analytical results is very good over a wide range of $`\lambda `$ values. In Fig. 3b, the profiles of $`u(x,y)`$ and $`w(x,y)`$ components of a dipole soliton, computed with numerical iteration methods, are displayed.

We note that in view of Eqs. (29) and (31), the integral in Eq. (28) is actually

$$_0^{\mathrm{}}r\varphi ^2h_2(r)𝑑r=(2\lambda _{2d}\lambda _{2v})_0^{\mathrm{}}r\varphi ^2𝑑r.$$
(33)
Inspection of the $`\lambda _{2d}(s)`$ and $`\lambda _{2v}(s)`$ curves in Fig. 1a immediately confirms that the integral (33) is always positive. Thus Eq. (28) holds only when $`p=\pm 1`$.

To summarize, we have shown that there are only two vector solitons of the system (4)–(5) that bifurcate from the cut-off frequency $`\lambda =\lambda _0(s)`$. They are either a vortex soliton (7) or a real-valued dipole soliton (8). The solutions are determined in terms of perturbation series expansions up to order $`ϵ^3`$. Both solitons are unique up to phase, position and rotation invariances. The analytical results are confirmed by numerical calculations. Computations of the perturbation series expansions prove the existence and uniqueness of vortex and dipole vector solitons observed numerically in .

## III Linear stability of vortex and dipole solitons

In this section, we study the linear stability of vortex and dipole vector solitons by spectral analysis, supplemented by numerical computations. Since the linearization operators differ for vortex and dipole solitons, we shall treat the two cases separately.

### A Vortex solitons

To study the linear stability of the vortex solitons (7) with $`n=0`$ and $`m=1`$, we linearize the system (1)–(2) with the perturbation in the form:

$$E_1=e^{iz}\left[\mathrm{\Phi }_u(r)+u_+(r)e^{in\theta +\sigma z}+\overline{u}_{}(r)e^{in\theta +\overline{\sigma }z}\right],$$
(34)
$$E_2=e^{i\lambda z+i\theta }\left[\mathrm{\Phi }_w(r)+w_+(r)e^{in\theta +\sigma z}+\overline{w}_{}(r)e^{in\theta +\overline{\sigma }z}\right].$$
(35)
The linearization problem can be written in the form:

$`i\sigma u_+`$ $`=`$ $`u_+^{\prime \prime }{\displaystyle \frac{1}{r}}u_+^{}+\left(1+{\displaystyle \frac{n^2}{r^2}}\right)u_+Vu_+V_{uu}(u_++u_{})V_{uw}(w_++w_{}),`$ (36)
$`i\sigma u_{}`$ $`=`$ $`u_{}^{\prime \prime }{\displaystyle \frac{1}{r}}u_{}^{}+\left(1+{\displaystyle \frac{n^2}{r^2}}\right)u_{}Vu_{}V_{uu}(u_++u_{})V_{uw}(w_++w_{}),`$ (37)
$`i\sigma w_+`$ $`=`$ $`w_+^{\prime \prime }{\displaystyle \frac{1}{r}}w_+^{}+\left(\lambda +{\displaystyle \frac{(n1)^2}{r^2}}\right)w_+Vw_+V_{uw}(u_++u_{})V_{ww}(w_++w_{}),`$ (38)
$`i\sigma w_{}`$ $`=`$ $`w_{}^{\prime \prime }{\displaystyle \frac{1}{r}}w_{}^{}+\left(\lambda +{\displaystyle \frac{(n+1)^2}{r^2}}\right)w_{}Vw_{}V_{uw}(u_++u_{})V_{ww}(w_++w_{}),`$ (39)
where

$$V=\frac{\mathrm{\Phi }_u^2+\mathrm{\Phi }_w^2}{1+s(\mathrm{\Phi }_u^2+\mathrm{\Phi }_w^2)},V_{uu}=\frac{\mathrm{\Phi }_u^2}{(1+s(\mathrm{\Phi }_u^2+\mathrm{\Phi }_w^2))^2},V_{uw}=\frac{\mathrm{\Phi }_u\mathrm{\Phi }_w}{(1+s(\mathrm{\Phi }_u^2+\mathrm{\Phi }_w^2))^2},V_{ww}=\frac{\mathrm{\Phi }_w^2}{(1+s(\mathrm{\Phi }_u^2+\mathrm{\Phi }_w^2))^2}.$$

The linearized problem can be formulated in the Hamiltonian form :

$$\pm i\sigma u_\pm =\frac{\delta h}{\delta \overline{u}_\pm },\pm i\sigma w_\pm =\frac{\delta h}{\delta \overline{w}_\pm },$$
(40)
where $`h`$ is the energy quadratic form associated with an eigenvector $`𝐮=(u_+,u_{},w_+,w_{})^T`$ and a linearized self-adjoint operator $``$ of the right-hand-sides of the system (36)–(39):

$$h=𝐮,𝐮=i\sigma _0^{\mathrm{}}r𝑑r\left(|u_+|^2|u_{}|^2+|w_+|^2|w_{}|^2\right).$$
(41)
The eigenvalue $`\sigma `$ is defined by the spectrum of the linearized problem (36)–(39) when $`𝐮(r)`$ is localized as $`r\mathrm{}`$ such that the integral (41) makes sense. The eigenvalues could be isolated or embedded into continuous spectrum of the system (36)–(39). The vortex soliton is linearly unstable if there exists an eigenvalue $`\sigma `$ for some $`n`$ such that $`\mathrm{Re}(\sigma )>0`$. We note that if $`(\sigma ,n,u_+,u_{},w_+,w_{})`$ is a solution of the linear system (36)–(39), so are $`(\overline{\sigma },n,\overline{u}_{},\overline{u}_+,\overline{w}_{},\overline{w}_+)`$, $`(\sigma ,n,u_{},u_+,w_{},w_+)`$ and $`(\overline{\sigma },n,\overline{u}_+,\overline{u}_{},\overline{w}_+,\overline{w}_{})`$. Thus, complex unstable eigenvalues $`\sigma `$ always come in quartets, while real and imaginary eigenvalues $`\sigma `$ always come in pairs. We also note that eigenmodes $`(\sigma ,n,u_+,u_{},w_+,w_{})`$ and $`(\overline{\sigma },n,\overline{u}_{},\overline{u}_+,\overline{w}_{},\overline{w}_+)`$ give the identical perturbation in (34) and (35), so do $`(\sigma ,n,u_{},u_+,w_{},w_+)`$ and $`(\overline{\sigma },n,\overline{u}_+,\overline{u}_{},\overline{w}_+,\overline{w}_{})`$.

According to the stability theory of solitary waves in the system of coupled NLS equations , only eigenvalues $`\sigma `$ with negative or zero values of the energy quadratic form (41) may bifurcate to the domain $`\mathrm{Re}(\lambda )>0`$, leading to instabilities. We shall apply this theory and study the spectrum of the linearization problem (36)–(39) at the cut-off frequency $`\lambda =\lambda _0(s)`$, when $`\mathrm{\Phi }_u(r)=u_0(r)`$ and $`\mathrm{\Phi }_w(r)=0`$. We will show that the vortex soliton is linearly stable in the neighborhood of the cut-off frequency $`\lambda _0(s)\lambda <\lambda _c(s)`$, where $`\lambda _c(s)1`$ is the instability threshold. In the limit $`\lambda \lambda _0(s)`$, the linearization problem decomposes into two linear problems:

$$\pm i\sigma u_\pm =u_\pm ^{\prime \prime }\frac{1}{r}u_\pm ^{}+\left(1+\frac{n^2}{r^2}\right)u_\pm \frac{u_0^2}{1+su_0^2}u_\pm \frac{u_0^2(u_++u_{})}{(1+su_0^2)^2},$$
(42)
and

$$\pm i\sigma w_\pm =w_\pm ^{\prime \prime }\frac{1}{r}w_\pm ^{}+(\lambda _0+\frac{(n1)^2}{r^2})w_\pm \frac{u_0^2}{1+su_0^2}w_\pm .$$
(43)
The first linear problem (42) is the stability problem of a scalar fundamental soliton $`u=u_0(r)`$ in a saturable medium. The linear stability of such solitons has been well established (see  for instance), thus unstable eigenvalues $`\sigma `$ do not exist in (42). The continuous spectrum of the system (42) is located at $`\mathrm{Re}(\sigma )=0`$ and $`|\mathrm{Im}(\sigma )|1`$. The continuous spectrum is irrelevant for stability of solitary waves in the system of coupled NLS equations, if no embedded eigenvalues with negative energy quadratic form (41) exist . The discrete spectrum of (42) consists of isolated eigenvalues $`\sigma `$ such that $`\mathrm{Re}(\sigma )=0`$ and $`|\mathrm{Im}(\sigma )|<1`$, including the zero eigenvalue at $`n=0`$ and $`n=\pm 1`$ with three eigenfunctions (18) and three generalized eigenfunctions. Additional eigenvalues for internal modes exist in (42) for $`\mathrm{Re}(\sigma )=0`$ and $`0|\mathrm{Im}(\sigma )|<1`$. These modes have been determined numerically in . Using those numerical results, we have found that the energy quadratic form (41) is positive for all internal modes of the system (42). For instance, only one internal mode with $`n=0`$ exists and has positive value of $`h`$ for $`s=0.5`$ (which corresponds to $`\omega =0.5`$ in , see Fig. 3). We have also checked that no embedded eigenvalues with $`|\mathrm{Im}(\sigma )|1`$ exist in the problem (42) for $`s=0.5`$.

The second linear problem (43) is uncoupled for $`w_+`$ and $`w_{}`$. Since the operator in the right-hand-side of Eq. (43) is self-adjoint, the spectrum of $`\sigma `$ is purely imaginary, i.e. $`\mathrm{Re}(\sigma )=0`$. The continuous spectrum of (43) is located at $`|\mathrm{Im}(\sigma )|>\lambda _0`$. Its discrete spectrum consists of isolated eigenvalues with $`\mathrm{Im}(\sigma )>\lambda _0`$ for $`w_+`$ and $`\mathrm{Im}(\sigma )<\lambda _0`$ for $`w_{}`$ and can be embedded into continuous spectrum. The discrete spectrum of (43) includes two zero eigenvalues at $`n=0`$ with eigenfunctions $`w_\pm =\varphi (r)`$, two zero eigenvalues at $`n=\pm 2`$ with eigenfunctions $`w_\pm =\varphi (r)`$, and two non-zero eigenvalues $`\sigma =\pm i(1\lambda _0)`$ at $`n=\pm 1`$ with eigenfunctions $`w_\pm =u_0(r)`$. The zero eigenvalues at $`n=0`$ are induced by the phase invariance of the $`E_2`$ equation, while the zero eigenvalues at $`n=\pm 2`$ are induced by the symmetry of the uncoupled problem (43). The eigenvalues $`\sigma =\pm i(1\lambda _0)`$ with $`n=\pm 1`$ are induced by the arbitrary polarizations in the fundamental vector soliton (6) at $`\lambda =\lambda _0(s)`$. The latter eigenvalues result in negative values of the energy quadratic form (41):

$$h=(1\lambda _0)_0^{\mathrm{}}u_0^2(r)r𝑑r<0.$$
(44)
We have checked numerically that (43) has no other discrete eigenvalues for $`s=0.5`$.

Applying the stability theory of solitary waves , we count eigenvalues of the problems (42) and (43) that produce negative and zero values of the energy quadratic form $`h`$. At $`\lambda =\lambda _0(s)`$, only two eigenvalues $`\sigma =\pm i(1\lambda _0)`$ give the negative energy (44). Several zero eigenvalues give zero energy at $`n=0,\pm 1,\pm 2`$. However, zero eigenvalues at $`n=0`$ and $`n=\pm 1`$ are preserved at $`\lambda >\lambda _0(s)`$ due to translation, rotation, and complex-phase symmetries of the system (1)–(2). Only two zero eigenvalues of the problem (43) at $`n=\pm 2`$ are not preserved by the symmetry and they can move out of zero for $`\lambda >\lambda _0(s)`$. We shall now consider the shift of these negative-energy and zero-energy eigenvalues for $`\lambda `$ near the cut-off frequency $`\lambda _0(s)`$.

We show first that the negative-energy eigenvalues $`\sigma =\pm i(1\lambda _0)`$ never bifurcate off the imaginary axis for $`\lambda >\lambda _0(s)`$ regardless whether they are embedded or not. Indeed, at any $`\lambda `$ value, the linearization problem (36)–(39) has the exact discrete eigenmode

$$u_+=0,u_{}=\mathrm{\Phi }_w(r),w_+=\mathrm{\Phi }_u(r),w_{}=0$$
(45)
for $`n=1`$ and $`\sigma =i(1\lambda )`$, and

$$u_+=\mathrm{\Phi }_w(r),u_{}=0,w_+=0,w_{}=\mathrm{\Phi }_u(r)$$
(46)
for $`n=1`$ and $`\sigma =i(1\lambda )`$. This result contrasts what happens in the system of coupled NLS equations with Kerr nonlinearities, where the negative-energy discrete eigenvalues, which are embedded in the continuous spectrum, bifurcate to the complex plane and lead to the instability .

We note that the exact eigenmodes (45) and (46) generate an approximate solution of the system (1)–(2):

$$E_1=\mathrm{\Phi }_u(r)e^{iz}\gamma \mathrm{\Phi }_w(r)e^{i\theta +i\lambda z}+\mathrm{O}(\gamma ^2),E_2=\mathrm{\Phi }_w(r)e^{i\theta +i\lambda z}+\gamma \mathrm{\Phi }_u(r)e^{iz}+\mathrm{O}(\gamma ^2),$$
(47)
where $`\gamma `$ is an arbitrary small parameter. Solution (47) is nothing but the original vortex vector soliton under a small rotation in the $`(E_1,E_2)`$ plane. This rotation leaves the original system (1) and (2) invariant and is one of the symmetries of the present problem. Thus, although solution (47) appears like internal oscillations of vortex solitons, this oscillation does not create any energy radiation and is fundamentally different from internal oscillations discussed in .

We show next that the zero eigenvalues at $`n=\pm 2`$ move to the imaginary axis (as conjugate pairs) as $`\lambda >\lambda _0(s)`$ and do not create any instability. We will use the perturbation series expansions and will present calculations only for the case $`n=2`$ (the case $`n=2`$ is similar). When $`\lambda `$ is close to $`\lambda _0(s)`$, we construct an approximate solution to the linearization problem (36)–(39) at $`n=2`$ in the form:

$$u_\pm =ϵu_\pm ^{(1)}(r)+\mathrm{O}(ϵ^3),w_+=\varphi (r)+ϵ^2w_+^{(2)}(r)+\mathrm{O}(ϵ^4),w_{}=ϵ^2w_{}^{(2)}(r)+\mathrm{O}(ϵ^4),\sigma =ϵ^2\sigma _2+\mathrm{O}(ϵ^4),$$
(48)
where $`ϵ`$ is the same small parameter as in expansions (14)–(16). Substituting (48) into the system (36)–(39), we find an exact solution at order $`ϵ`$: $`u_+^{(1)}=u_{}^{(1)}=u_{22}(r)`$, where $`u_{22}(r)`$ solves the equation (21). At order $`ϵ^2`$, we need to solve the non-homogeneous equation for $`w_+^{(2)}(r)`$:

$`w_+^{(2)\prime \prime }+{\displaystyle \frac{1}{r}}w_+^{(2)\prime \prime }\left(\lambda _0+{\displaystyle \frac{1}{r^2}}\right)w_+^{(2)}+{\displaystyle \frac{u_0^2}{1+su_0^2}}w_+^{(2)}=\left(\lambda _2i\sigma _2\right)\varphi {\displaystyle \frac{2\varphi (u_0u_{20}+u_0u_{22}+\varphi ^2)}{(1+su_0^2)^2}}.`$ (49)
The solvability condition for this equation can be simplified by virtue of Eq. (29), and we find that the eigenvalue coefficient $`\sigma _2`$ is given as

$$\sigma _2=2i\frac{_0^{\mathrm{}}r\varphi ^2h_2𝑑r}{_0^{\mathrm{}}r\varphi ^2𝑑r}.$$
(50)
Utilizing Eq. (33), we see that

$$\sigma _2=2i\left[2\lambda _{2d}(s)\lambda _{2v}(s)\right],$$
(51)
whose imaginary part is positive from Fig. 1a. Thus, the energy quadratic form of the bifurcated eigenmode (48) and (51) (up to the order $`ϵ^2`$) is negative:

$$h=ϵ^2\mathrm{Im}(\sigma _2)_0^{\mathrm{}}r\varphi ^2𝑑r<0.$$
(52)
The analytical eigenvalue formula (48) and (51) at $`n=2`$ is plotted in Fig. 4 versus $`\lambda `$ for $`s=0.5`$ (dash-dotted line). Numerically, we have determined these eigenvalues for $`s=0.5`$ and various values of $`\lambda `$, and the results are plotted in Fig. 4 (solid line) as well. When $`\lambda `$ is close to $`\lambda _0`$, the analytical formula agrees well with the numerical values.

We have shown above that the two zero eigenvalues of the system (43) at $`n=\pm 2`$ move to the imaginary axis when $`\lambda >\lambda _0(s)`$, while the two non-zero negative-energy eigenvalues at $`n=\pm 1`$ remain on the imaginary axis. Thus, we conclude that vortex solitons are linearly stable near the cut-off frequency $`\lambda =\lambda _0(s)`$, i.e. vortex solitons with small vortex components are linearly stable. This result confirms the conclusions of  and does not support conclusions of  where all vortex vector solitons were claimed to be linearly unstable.

Unstable eigenvalues of vortex solitons may appear far away from the cut-off frequency $`\lambda _0(s)`$. Indeed, the two imaginary eigenvalues $`\sigma `$ for $`n=\pm 2`$ that bifurcate from zero eigenvalues at $`\lambda >\lambda _0(s)`$ have negative energy (52). When these eigenvalues collide with eigenvalues of positive energy or with continuous spectrum, the oscillatory instability may arise . We confirm this scenario and compute unstable eigenvalues $`\sigma `$ of the linear system (36)–(39) with the numerical shooting method. The unstable eigenvalues are found exactly at $`n=\pm 2`$ and are shown in Fig. 4 for $`s=0.5`$. The unstable eigenvalues appear when $`\lambda >\lambda _c0.402`$, where $`\lambda _c`$ denotes the frequency for onset of instability. These results agree with Fig. 3 of Ref. , where the unstable eigenvalues were found from time-integration of the linearized equations under random-noise initial perturbations. Fig. 5 displays our numerical solutions for the vortex soliton and the corresponding unstable eigenfunction at $`s=0.5`$, $`\lambda =0.5`$ and $`n=2`$. Thus, for the case $`s=0.5`$, unstable eigenvalues exist at $`\lambda >\lambda _c0.402`$, while vortex vector solitons exist at $`\lambda >\lambda _00.2622`$, see Fig. 4. In the interval $`\lambda _0<\lambda <\lambda _c`$, i.e., $`0.2622<\lambda <0.402`$ for $`s=0.5`$, unstable eigenvalues do not exist and the vortex solitons are linearly stable.

We conclude this analysis with two remarks. First, it follows from Fig. 4 for $`s=0.5`$ that the eigenvalues $`\sigma `$ at $`n=\pm 2`$ merge into the continuous spectrum at $`\lambda 0.396`$, while unstable eigenvalues appear at $`\lambda =\lambda _c0.402`$. Our numerical results are inconclusive as to what happens in the narrow interval $`0.396<\lambda <0.402`$. This problem is left open for future studies. And second, when $`\lambda `$ is further away from the cut-off frequency $`\lambda _0(s)`$, the vector vortex solution bifurcates into scalar vortex solutions with $`u=0`$ and $`w=\mathrm{\Phi }_w(r)e^{i\theta }`$, see . The scalar vortex soliton has additional unstable eigenmodes at $`|n|2`$ which have smaller growth rates (see ). We do not study this bifurcation where the family of vector vortex solitons terminates, nor the number of unstable eigenvalues of vector vortex solitons near this bifurcation.

### B Dipole soliton

To study the linear stability of the dipole solitons (8), we linearize the system (1)–(2) with the perturbation:

$$E_1=e^{iz}\left[U(x,y)+\left(u_r(x,y)+u_i(x,y)\right)e^{\sigma z}+\left(\overline{u}_r(x,y)\overline{u}_i(x,y)\right)e^{\overline{\sigma }z}\right],$$
(53)
$$E_2=e^{i\lambda z}\left[W(x,y)+\left(w_r(x,y)+w_i(x,y)\right)e^{\sigma z}+\left(\overline{w}_r(x,y)\overline{w}_i(x,y)\right)e^{\overline{\sigma }z}\right].$$
(54)
Here $`u_r,u_i,w_r`$ and $`w_i`$ are complex functions and are very small. The linearization problem is then written in the form:

$`i\sigma u_i`$ $`=`$ $`\mathrm{\Delta }u_r+u_r\left(V+2V_{uu}\right)u_r2V_{uw}w_r,`$ (55)
$`i\sigma u_r`$ $`=`$ $`\mathrm{\Delta }u_i+u_iVu_i,`$ (56)
$`i\sigma w_i`$ $`=`$ $`\mathrm{\Delta }w_r+\lambda w_r\left(V+2V_{ww}\right)w_r2V_{uw}u_r,`$ (57)
$`i\sigma w_r`$ $`=`$ $`\mathrm{\Delta }w_i+\lambda w_iVw_i,`$ (58)
where

$$V=\frac{U^2+W^2}{1+s(U^2+W^2)},V_{uu}=\frac{U^2}{(1+s(U^2+W^2))^2},V_{uw}=\frac{UW}{(1+s(U^2+W^2))^2},V_{ww}=\frac{W^2}{(1+s(U^2+W^2))^2}.$$

The linearized problem can be formulated in the same Hamiltonian form (40) with the energy quadratic form :

$$h=i\sigma _{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\left(\overline{u}_ru_i+\overline{u}_iu_r+\overline{w}_rw_i+\overline{w}_iw_r\right)𝑑x𝑑y.$$
(59)
At the cut-off frequency $`\lambda =\lambda _0(s)`$, the same analysis as for the vortex solitons shows existence of a pair of eigenvalues $`\sigma =\pm i(1\lambda _0)`$ with negative values of $`h`$ and a number of zero eigenvalues with zero values of $`h`$. We show again that the eigenvalues $`\sigma =\pm i(1\lambda _0)`$ with negative energy never bifurcate into complex domain for $`\lambda >\lambda _0(s)`$. Indeed, for any $`\lambda `$ value, the linearization problem (55)–(58) has the exact solution

$$u_r=W(x,y),u_i=W(x,y),w_r=U(x,y),w_i=U(x,y)$$
(60)
at $`\sigma =i(1\lambda )`$, and

$$u_r=W(x,y),u_i=W(x,y),w_r=U(x,y),w_i=U(x,y)$$
(61)
at $`\sigma =i(1\lambda )`$.

We study the zero eigenvalues of the system (55)–(58) with perturbation series expansions for $`\lambda >\lambda _0(s)`$:

$`\sigma =ϵ^2\sigma _2+\mathrm{O}(ϵ^4),u_r=ϵu_r^{(1)}(r,\theta )+\mathrm{O}(ϵ^3),u_i=\mathrm{O}(ϵ^3),`$ (62)
$`w_r=w_r^{(0)}(r,\theta )+ϵ^2w_r^{(2)}(r,\theta )+\mathrm{O}(ϵ^4),w_i=w_i^{(0)}(r,\theta )+ϵ^2w_i^{(2)}(r,\theta )+\mathrm{O}(ϵ^4).`$ (63)
Here $`ϵ`$ is the same small parameter as in expansions (14)–(16), while the functions $`w_{r,i}^{(0)}(r,\theta )`$ are linear combinations of the eigenfunctions of the null space of the problem (55)–(58) at $`\lambda =\lambda _0(s)`$:

$$w_r^{(0)}=c_1\varphi (r)\mathrm{cos}\theta +c_2\varphi (r)\mathrm{sin}\theta ,w_i^{(0)}=d_1\varphi (r)\mathrm{cos}\theta +d_2\varphi (r)\mathrm{sin}\theta ,$$
(64)
where $`c_1`$, $`c_2`$, $`d_1`$, and $`d_2`$ are constants. Substituting (63) into the system (55)–(58), we find an exact solution at order $`ϵ`$:

$$u_r^{(1)}=c_1u_{20}(r)+\left(c_1\mathrm{cos}2\theta +c_2\mathrm{sin}2\theta \right)u_{22}(r).$$
(65)
where $`u_{20}(r)`$ and $`u_{22}(r)`$ solve the problems (20) and (21). At order $`ϵ^2`$, four solvability conditions are needed for solving the non-homogeneous equations for $`w_r^{(2)}(r,\theta )`$ and $`w_i^{(2)}(r,\theta )`$. Using Eq. (31), we transform the four solvability conditions to the form:

$`\sigma _2c_1=0,\sigma _2d_2=0,i\sigma _2c_2{\displaystyle _0^{\mathrm{}}}r\varphi ^2𝑑r=d_2{\displaystyle _0^{\mathrm{}}}r\varphi ^2h_2𝑑r,i\sigma _2d_1{\displaystyle _0^{\mathrm{}}}r\varphi ^2𝑑r=c_1{\displaystyle _0^{\mathrm{}}}r\varphi ^2\left(2h_1+h_2\right)𝑑r.`$ (66)
If $`\sigma _2=0`$, then $`c_1=d_2=0`$, while $`c_2`$, $`d_1`$ are arbitrary constants. Thus, the zero eigenvalue persists in the system (55)–(58) for $`\lambda >\lambda _0(s)`$ with two eigenfunctions $`w_r=\varphi (r)\mathrm{sin}\theta `$ and $`w_i=\varphi (r)\mathrm{cos}\theta `$. The two eigenfunctions are related to the symmetries of the system (1)–(2) with respect to rotation in $`\theta `$ and shift of the complex phase. If $`\sigma _20`$ however, the system (66) has only the trivial solution: $`c_1=c_2=d_1=d_2=0`$. Therefore, the other two zero eigenvalues do not bifurcate to the imaginary axis but simply disappear for $`\lambda >\lambda _0(s)`$.

We have analytically proved above that the dipole solitons are linearly stable in the neighborhood of the cut-off frequency $`\lambda _0(s)`$. Moreover, contrary to vortex solitons, there are only two eigenvalues of negative energy for $`\lambda >\lambda _0(s)`$, and they remain on the imaginary axis for all $`\lambda `$ values \[see (60) and (61\]. Thus, we conjecture that the dipole solitons are linearly stable in the whole domain of their existence. This conjecture is in agreement with the numerical work of . We again confirm this result by numerical simulations of the system (1)–(2) linearized around the dipole soliton (8). For $`s=0.5`$, we have simulated the linearized system for several values of $`\lambda `$ between $`\lambda =0.3`$ and $`\lambda =0.85`$. We did not find any instability in the linearized system. Since $`\lambda =0.3`$ is close to the cut-off frequency $`\lambda _0=0.2622`$ and $`\lambda =0.85`$ is close to the end frequency $`\lambda =1`$, we conclude that dipole solitons are indeed linearly stable in the whole existence interval.

## IV Nonlinear evolution of perturbed vortex solitons

Here we study the nonlinear evolution of perturbed vortex solitons. The unstable vortex soliton under small random-noise perturbations was found in  to break up into a rotating dipole vector soliton. We will show below that such a breakup scenario holds only when the vortex component of the vortex soliton is below certain threshold. Above that threshold, unstable vortex solitons break up into two rotating fundamental vector solitons instead. We will also show that the vortex solitons with small vortex components are not only linearly but also nonlinearly stable.

We consider first the nonlinear evolution of linearly stable vortex solitons. For this purpose, we have simulated the system (1)–(2) starting with a linearly-stable vortex soliton under various types of small initial perturbations such as random-noise and deterministic ones. We have found that the vortex solitons are also nonlinearly stable for all small perturbations. To demonstrate, we select $`s=0.5`$ and $`\lambda =0.38`$, where the vortex soliton has been shown to be linearly stable (see Fig. 4). As initial perturbations, we choose

$$E_1(r,\theta ,0)=(1+\alpha )\mathrm{\Phi }_u(r,\lambda ),E_2(r,\theta ,0)=(1+\alpha )\mathrm{\Phi }_w(r,\lambda )e^{i\theta },$$
(67)
where $`\alpha `$ is a small perturbation parameter, that measures amplification of the vortex soliton by a factor $`1+\alpha `$. The simulation result with $`\alpha =0.05`$ is shown in Fig. 6. This figure shows that the perturbed vortex soliton persists the nonlinear evolution and exhibits little change of shape even after 300 diffraction lengths. This clearly confirms the linear and nonlinear stability of the vortex soliton with $`s=0.5`$ and $`\lambda =0.38`$. Other perturbations to this soliton give similar evolution results.

We study next the nonlinear evolution of linearly unstable vortex solitons. We have shown in Sec. III A that these solitons possess two unstable eigenmodes $`(\sigma ,n,u_+,u_{},w_+,w_{})`$ and $`(\overline{\sigma },n,\overline{u}_{},\overline{u}_+,\overline{w}_{},\overline{w}_+)`$ with $`n=2`$. The $`\sigma `$ versus $`\lambda `$ graph is shown in Fig. 4 for $`s=0.5`$, while unstable eigenfunctions $`(u_+,u_{},w_+,w_{})`$ for $`s=0.5`$ and $`\lambda =0.5`$ are displayed in Fig. 5. However, we recognize that these two unstable eigenmodes are equivalent in view of Eqs. (34) and (35). Thus, any small initial perturbation to the vortex soliton is projected onto this unstable eigenmode which grows exponentially, while the rest of the initial perturbation disperses away. For convenience, we choose the initial perturbation to be exactly this unstable eigenmode, i.e.,

$$E_1(r,\theta ,0)=\mathrm{\Phi }_u(r)+\alpha (u_+(r)e^{2i\theta }+\overline{u}_{}(r)e^{2i\theta }),$$
(68)
$$E_2(r,\theta ,0)=e^{i\theta }\left[\mathrm{\Phi }_w(r)+\alpha (w_+(r)e^{2i\theta }+\overline{w}_{}(r)e^{2i\theta })\right],$$
(69)
where $`\alpha `$ is a small perturbation parameter. The advantage of this special perturbation is that it shortens the distance for the breakup of the vortex soliton and reduces the radiation noise in the nonlinear evolution of the perturbed solution.

We have discovered two breakup scenarios of the unstable vortex soliton with the initial perturbation (68)–(69). We confirm that unstable vortex solitons with relatively small vortex components indeed break up into a rotating dipole soliton, in agreement with . However, when the vortex component increases above a certain threshold, an unstable vortex soliton breaks up into two rotating fundamental vector solitons rather than one dipole soliton. For example, when $`s=0.5`$ and $`\alpha =0.05`$, the vortex soliton breaks up into a dipole soliton when $`0.402<\lambda _{}^<\mathrm{\hspace{0.25em}0.45}`$, and into two fundamental vector solitons when $`\lambda >0.45`$. Indeed, when $`\lambda =0.45`$ (where the vortex component is relatively small), the time evolution of the perturbed vortex soliton is plotted in Fig. 7. It is seen that this soliton breaks up into a rotating dipole soliton. But when $`\lambda =0.5`$ (where the vortex component is bigger), the time evolution is shown in Fig. 8. Here two rotating fundamental vector solitons are formed after the breakup of the unstable vortex soliton. We have also found that these breakup scenarios are insensitive to the type of initial perturbation imposed, because we have simulated the evolutions with different values of $`\alpha `$ in (68) and (69) as well as with other forms of initial perturbations such as random noise, but the breakup scenarios do not change. To check the numerical accuracy of our simulations, we have used more grid points and wider $`(x,y)`$ intervals and obtained identical results. Furthermore, our results conserve energies of the $`E_1`$ and $`E_2`$ components very well.

Intuitively, it is not difficult to understand the above two breakup scenarios of unstable vortex solitons. When the vortex component of the vortex soliton is small, the instability (with $`n=2`$) breaks up the vortex ($`E_2`$) component into two weak humps, while it does not significantly affect the single-hump shape of the fundamental ($`E_1`$) component since $`E_1`$’s initial amplitude is much higher. During the subsequent evolution, the two humps of the $`E_2`$ component are too weak to break the $`E_1`$ component into two pieces, thus the solution relaxes into a dipole soliton instead of two fundamental solitons. However, when the vortex component of the vortex soliton is sufficiently large, the fundamental component becomes small (see Fig. 2 and ). In this case, instability breaks up both the vortex and fundamental components into two pieces, and two fundamental solitons are formed then.

## V Summary and discussion

To summarize, we have studied both analytically and numerically the existence, uniqueness, and stability of vortex and dipole vector solitons in saturable optical materials in (2+1) dimensions. We have shown that the analytical expressions for vortex and dipole vector solitons can be constructed with perturbation series expansions near the cut-off frequency $`\lambda =\lambda _0(s)`$. We have also shown that only two vector solitons bifurcate from the same cut-off frequency, which are vortex and dipole solitons. Furthermore, we have proved that both vortex and dipole solitons are linearly stable when the vortex and dipole components are small. As the vortex and dipole components increase, the family of vortex vector solitons becomes linearly unstable, while that of dipole vector solitons remains linearly stable in the entire existence domain. We have also shown that unstable vortex solitons break up into a rotating dipole soliton only when the vortex component is relatively small. When the vortex component crosses a certain threshold, the vortex soliton breaks up into two rotating fundamental vector solitons instead. We expect that our results are significant not only for studies of spatial vector solitons in a saturable nonlinear medium but also for studies of Bose-Einstein condensation.

In this paper, we have studied only the simplest vortex and dipole vector solitons which bifurcate from the fundamental $`u`$ and small $`w`$ components. One natural question to ask is the existence and stability of other vortex and multi-pole vector solitons. The perturbation series expansion method developed in this paper is powerful for a systematic study of general vortex and multi-pole vector solitons near their bifurcation points. But this problem lies outside the scope of the present article. We note, however, that vortex solitons (7) with $`|n|>0`$ and $`|m|>0`$ exist, and they are expected to be linearly always unstable because each component has non-zero charge and is linearly unstable by itself . This expectation is consistent with our preliminary numerical simulations on vortex solitons with charges such as $`n=1`$ and $`m=1`$.

Recently, three-component vortex and dipole vector solitons in a saturable medium have been investigated . The authors found that those solitons are linearly unstable provided their total topological charge is non-zero. In view of our results in this paper, that conclusion needs modification. We plan to study that system carefully in the near future.

## Acknowledgments

The authors appreciate helpful discussions with Yu. Kivshar, Z. Musslimani, B. Sandstede and D. Skryabin. The work of D.P. was supported in part by the NSERC grant No. 5-36694 and CFI grant No. 5-26773. The work of J.Y. was supported in part by the Air Force Office of Scientific Research under contract F49620-99-1-0174, and by the National Science Foundation under grant DMS-9971712.