text:
# Coherent dimer formation near Feshbach resonances in Bose-Einstein condensates

## Abstract

The results of a recent experiment with 85Rb Bose-Einstein condensates are analyzed within the mean-field approximation including dissipation due to three-body recombination. The intensity of the dissipative term is chosen from the three-body theory for large positive scattering lengths. The remaining number of condensed atoms in the experiment, obtained with applied magnetic field pulses, were used to adjust the intensity of the dissipative term. We found that the three-body recombination parameter depends on the pulse rise time; i.e., for longer rise times the values found become consistent with the three-body theory, while for shorter pulses this coefficient is found to be much larger. We interpret this finding as an indication of a coherent formation of dimers.

Bose-Einstein condensation of dilute atomic gases has been observed for the first time in 87Rb atoms Cornell . Since then, it has been condensed atoms of 23Na Ketterle , 7Li Hulet , 1Fried , metastable 4He Cohen , 85Rb Cornish , 41Modugno  and, more recently, 133Cs Tinoweber . In Ref. Tiesinga , it was theoretically predicted that Feshbach resonances could vary the scattering length of atoms in systems of dilute alkali gases over a wide range of values. A Feshbach resonance occurs when the quasimolecular bound state energy is tuned to the energy of two colliding atoms, applying an external magnetic field to the system. The phenomenon was first realized in a Bose-Einstein condensate in Ref. Inouye . This opened the possibility of exploring new regimes of BEC, allowing changes of the two-body scattering length from negative to positive and also from zero to infinity. This technique was used for the condensation and collapse control of 85Rb atoms in the hyperfine state ($`F=`$ 2, $`m_F=2`$Cornish ; Roberts ; Donley . It was also demonstrated in Ref. Stenger  that strongly enhanced inelastic three-body collisions occurs near Feshbach resonances. In a more recent experiment, in Ref. Claussen , it was explored the region of very large scattering lengths (up to $`4000`$ Bohr radius).

The scattering length $`a`$ has been observed to vary as a function of the magnetic field $`B`$, according to the theoretical prediction Moerdijk , as

$$a=a_b\times \left(1\frac{\mathrm{\Delta }}{BB_r}\right),$$
(1)
where $`a`$ is the scattering length, $`a_b`$ is the background scattering length, $`B_r`$ is the resonance magnetic field of 85Rb, and $`\mathrm{\Delta }`$ is the resonance width.

In the case of 85Rb, one has resonance width $`\mathrm{\Delta }`$ 11.0 G, resonance field $`B_r`$ 154.9 G and background scattering length $`a_b450a_0`$ Claussen , where $`a_0`$ is the Bohr radius. So, given the experimental functional dependence $`B=B(t)`$, one can determine $`a=a(t)`$ and, consequently, the dynamics of the system in such physical conditions. Further, mainly in strong interaction regime, a Bose-Einstein condensate shows inelastic loss processes that cause its depletion. The dominant process of losses has been verified to be the three-body recombination Claussen ; Burt , with a time dependence concerning a simple constant rate equation. All earlier observations in BEC experiments were consistent with a description of mean field including such a loss process. But, in experiments with 85Rb realized in the strong interaction regime Claussen , this picture indicates model breakdown. Bose-Einstein condensates initially stable were submitted to magnetic field pulses carefully controled in the vicinity of 85Rb Feshbach resonance, aiming to test the strongly interacting regime for diluteness parameter $`\chi =na^3`$ varying from $`\chi =`$ 0.01 to $`\chi =`$ 0.5. The loss of atoms from BEC occurred in impressively short time scales (up to two hundreds of $`\mu `$s) and disagrees with previous theoretical predictions Cornish . Such experiments reveal higher loss of atoms in shorter magnetic field pulses applied on BEC and, previously, one knew that as longer is the time spent near a Feshbach resonance as higher is the loss of atoms from BEC Cornish  (consistent with a mean field approach where the inelastic loss term has a constant dissipative rate). According to Ref. Claussen , the results could indicate the existence of a new physics, that cannot be described by the Gross-Pitaevskii (GP) formalism. Motivated by this observed discrepancy, we investigate the dynamics of 85Rb Bose-Einstein condensates submitted to such external conditions and time scales, when we vary the s-wave two-body scattering length in the region of strong interaction. So, we consider a generalized mean field approach that includes the time dependence for both mean field coupling and three-body recombination parameter. Previous results considered a mean field approach with constant mean rate for the three-body recombination and constant value for the mean field coupling Kagan ; Saito . In the present work, our first task is to reproduce the experimental data with time dependent parameters or at least verify possible limitations of the time dependent mean field approach. Next, we consider the magnitude of the recombination rate as described in literature Fedichev ; Nielsen1 ; Esry ; Bedaque ; Tomio ; Nielsen2 ; aform1 ; aform2 , but with time functional dependence.

For describing a BEC in the framework of the mean field approximation, we remind that such an approach is valid for very dilute systems when the average inter particle distances $`d`$ are much larger than $`|a|`$ and the particle wavelengths are much larger than $`d`$ Yukalov ; Dalfovo . Besides, it is important to pay attention to the time scales present in the system: the physical conditions must not change fast enough in order to allow the replacement of a true interatomic potential by the contact interaction. It is possible in principle for the rate of change to be larger than $`\mathrm{}/ma^2`$ for extremely very short changes in the interactions Legget , but it is reasonable to assume valid this time dependent approach at least for longer pulses.

We begin our description from an effective Lagrangian of the nonconservative system as in Abdullaev , in which one describes the dynamics of a trapped Bose-Einstein condensate in spherical symmetry with such a GPE generalization, in which one also considered losses from BEC by three-body recombination. This Lagrangian leads to the following equation of the system:

$`i\mathrm{}{\displaystyle \frac{\mathrm{\Psi }}{t}}=\left[{\displaystyle \frac{\mathrm{}^2}{2m}}^2+{\displaystyle \frac{m\omega ^2r^2}{2}}+U_0|\mathrm{\Psi }|^2\mathrm{i}\mathrm{}{\displaystyle \frac{K_3}{4}}|\mathrm{\Psi }|^4\right]\mathrm{\Psi },`$ (2)
where $`\omega `$ is the geometric mean trapping frequency, $`U_0U_0(t)=4\pi \mathrm{}^2a(t)/m`$ and $`K_3K_3(t)`$ is the recombination loss parameter. The wave-function $`\mathrm{\Psi }=\mathrm{\Psi }(\stackrel{}{r},t)`$ is normalized to the number of atoms $`N`$. The three-body recombination rate $`K_3`$ is here introduced for describing atomic losses from the condensate when three atoms scatter to form a molecular bound state (dimer) and a third atom; so, the kinetic energy of the final state particles allows them to escape from the trap. Other nonconservative processes as amplification from thermal cloud and dipolar relaxation are neglected, since the latter has a much smaller effect than three-body recombination Claussen ; Burt  and in JILA experiments Claussen  the thermal cloud is negligible (only 1000 atoms in a sample of 17500).

A theoretical prediction of $`K_3`$ is a hard task since it is sensitive to the detailed behavior of the interaction potential aform1 ; aform2 . However, such a calculation becomes simpler if we consider that $`a`$ is the only important length scale (reasonable in the weakly bound s-wave state limit) and this has been considered in many works Fedichev ; Nielsen1 ; Esry ; Bedaque ; Tomio ; Nielsen2 ; aform1 ; aform2 . Following Ref. Fedichev , the recombination rate is written as

$$K_3(t)\kappa \frac{\mathrm{}}{m}[a(t)]^4,$$
(3)
where $`\kappa `$ should correspond to the universal value $`\kappa =`$ 3.9. But in Ref. Nielsen1  and in Ref. Esry , it was found $`0\kappa 65`$ and in Ref. Bedaque , $`0\kappa 67.9`$, that is, $`\kappa `$ is not universal (in general $`K_3`$ depends on a three-atom scaleTomio ).

In Fig. 1, we schematically give two characteristic pulses employed in JILA, in which $`B_0`$, $`B_h`$, $`B_m`$, $`t_r`$ and $`t_h`$, correspond to the initial field, the hold field, the final field, the rise time and the hold time of the pulse, respectively. For describing the dynamics of condensates subjected to (1), we consider in our calculations the same experimental parameters and conditions used in Claussen : as $`a`$ is known to be a function of the magnetic field $`B`$ by means of (1), we only use time dependent shapes of experimental magnetic field pulses employed in JILA (triangular and trapezoidal pulses). The hold field is $`B_h=`$ 156.7 G, and the end field in which one measures the remaining number of particles $`N_r`$ of the system is $`B_m=`$ 164.5 G at $`t`$ = 700 $`\mu `$s. Further, we put initial field $`B_0166`$ G, corresponding to a harmonic oscillator state of the system ($`a0`$), applied to an initial sample of $`N_0`$=16500 condensed atoms of 85Rb. Further, in our approach we have used spherical symmetry with mean geometric frequency $`\omega =(\omega _r^2\omega _z)^{1/3}`$, for simulating the cylindrical geometry of JILA (radial: $`\omega _r=2\pi \times 17.5`$ Hz and axial: $`\omega _z=2\pi \times 6.8`$ Hz). Our time dependent calculations started with a Gaussian shape wave-function which we numerically evolve by means of Eq. (2), using the Crank-Nicolson algorithm, as in Refs. caos ; pre . We analyzed the loss of condensed atoms like in Claussen , by considering hold times from $`t_{hold}=`$ 0 (triangular) or units of $`\mu `$s (shorter trapezoidal pulses) to nearly hundreds of $`\mu `$s (longer trapezoidal pulses). The behavior of the scattering length, as a function of total time of the pulse, in the region of strong interaction atom-atom, follows similarly the behavior of the employed field pulses given in Fig. 1. In triangular pulses, there is a sharp peak in the resonance region; and a plateau with maximum value for $`a`$ and $`\xi `$ (minimum value of the field), when $`B=B_h`$ is kept constant during the time $`t_h`$.

In the upper frame of Fig. 2, we used symmetric rise and fall times ($`t_r=12.5\mu `$s, $`t_r25\mu `$s and $`t_r75\mu `$s), to determine the remaining fraction of atoms $`N_r/N_0`$ in the 85Rb BEC as function of the hold time, by adjusting our curve with the first point of Fig. 2 in Ref. Claussen . For $`t_r=12.5\mu `$s, we found $`\kappa 1800`$, very far away of the values described earlier in the literature Fedichev ; Nielsen1 ; Esry ; Bedaque ; Tomio ; Nielsen2 ; aform1 ; aform2 . The results show the same exponential decay and good concordance with experimental data, mainly for short $`t_h`$ (circles in upper frame of Fig. 2). However, as one can realize, for longer hold times, experimental data point out a higher dissipation when compared with our simulations. For other short rise times in this frame or longer rise times (lower frame), we have similar behavior, but we have to decrease $`\kappa `$ for better adjusting to experimental data. So, the comparison with experiments show that $`\kappa `$ depends significantly on values of $`t_r`$ and $`t_h`$. As we know Stoof , the mean field approach should make more adequately if we were in a slower process. So, we also tried to verify if our results would give a lower value of $`\kappa `$ (inside interval described in literature) if we calibrated our calculation with the last point of the longest pulse of JILA Claussen  ($`t_r`$ 252.6 $`\mu `$s). Really, we found $`\kappa 100`$ for this case (closer to values described in literature Fedichev ; Nielsen1 ; Esry ; Bedaque ; Tomio ; Nielsen2 ; aform1 ; aform2 ) and the results reproduce the experimental data for longer rise time but they do not make very well for shorter rise time, as we can observe in lower frame of Fig. 2. The very large value of $`\kappa `$ leads us to conclude that there is a coherent formation of dimers, that occurs up to nearly $`t_r`$ 100$`\mu s`$, as one can realize observing the lower frame of Fig. 2 for shorter rise times. It is physically sensible that, for shorter pulses, the presence of coherence in the formation of dimers would be more plausible than for longer pulses. So, the present calculation included the variation of $`\kappa `$ with the time parameters of the pulse and the amplification of the three-body recombination rate can be associated with the coherent formation of dimers in the inelastic collisions. Together with the significant loss from the coherent formation of dimers, one should also observe a burst of atoms carrying the excess of energy, which for 85Rb with the maximum value of $`a=4000a_0`$, would be above 70nK. Indeed, in the JILA experiment, it was seen a significant number loss from the BEC for pulses lasting only few tens of microseconds, which were accompanied by a burst of few thousand relatively hot ($``$ 150nK) atoms that remained in the trap Claussen . In our description, for each hot atom one dimer is also formed. Therefore, the burst of atoms should be accompanied by a burst of weakly bound dimers. Using the observed temperature and momentum conservation of the recombination process, we predict that the dimers are also relatively hot ($``$ 75nK).

We note that, if we consider a constant value of $`\kappa `$, in all the cases (for any choices of $`t_h`$) we obtain a decreasing behavior of $`N_r`$, as we increase the rise time $`t_r`$. This is in contrast with the experimental results. So, we consider to adjust the values of $`\kappa `$ that approximately better describe the experimental data for each given rise time $`t_r`$; i.e., $`\kappa `$ is taken as a function of the rise time ($`\kappa =\kappa (t_r)`$). Our results are shown in Table I. There is an obvious uncertainty in the given numbers of Table I, that are related to our approximate theoretical fitting and experimental data fluctuations. However, based on such results, the behavior of $`\kappa (t_r)`$ can be approximately described by a linear or exponential decreasing function. Here we consider the following simple functional time-rise dependences of $`\kappa `$:

$`\kappa (t_r)`$ $`=`$ $`2300\mathrm{exp}(0.01\times \omega t_r)\text{or}`$ (4)
$`\kappa (t_r)`$ $`=`$ $`19007\omega t_r.`$ (5)
With such decreasing functional time-rise dependence of $`\kappa `$, we have verified, as shown in Fig. 3, that we have the qualitative behavior observed in the experimental results given in Ref. Claussen  of the remaining number of atoms $`N_r`$ versus $`t_r`$. We have considered several values of the hold time $`t_h`$. For very small hold time, we can also verify the same experimental results that presents a minimum of $`N_r`$ as a function of $`t_r`$. With a linear functional $`t_r`$ dependence of $`\kappa `$, the same qualitative behavior can also be reproduced; which differs quantitatively from the exponential behavior that we show. So, we conclude that the time dependent mean field approach can describe all the experimental data, if the three-body recombination coefficient depends on the rise time in this short time scale. Such dependence on the rise time can be explored, once the uncertainties in the experimental results are reduced and by considering a better fitting of data (improving the values given in Table I). In our interpretation, the higher values of $`\kappa `$ for smaller values of $`t_r`$ are indicating the coherent formation of another species (dimers) in the condensate.

In summary, we report in this work indications based on our calculations that a recent experiment realized in JILA Claussen  is evidencing the coherent formation of dimers from inelastic collisions in 85Rb Bose-Einstein condensates (BEC). We have solved numerically the nonconservative Gross-Pitaevskii equation in spherical symmetry, for condensed systems with very large repulsive two-body interaction, varying in time, due to application of magnetic field pulses, according to Eq. (1). According to the theoretical predictions of three-body recombination rates, we used a dissipation parameter proportional to the quartic power of the scattering length $`a(t)`$ and the observed experimental pulse shapes to calculate the time evolution of the remaining number of atoms in the condensate. We studied this observable as a function of the hold and rise times. The experimental results of $`N_r`$ versus $`t_h`$ can be described in the mean-field approach only with very large $`\kappa `$, when the rise time is small. This indicates a coherent formation of dimers in BEC. For longer pulses, when the coherent dimer formation tends to disappear, we found that $`\kappa `$ approaches the maximum value given by the theoretical predictions for large scattering lengths, $`\kappa 70`$ Nielsen1 ; Bedaque . We parameterize this surprising behavior of the three-body recombination rate considering $`\kappa =\kappa (t_r)`$ and so it was possible to describe the property of a lower dissipation for longer pulses. Finally, it is natural to see a burst of relatively hot atoms and dimers (carrying the excess of dimer binding energy) accompanying a significant loss from the condensate for short pulses when the coherent dimer formation occurs.

We thank partial support from Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) of Brazil.


format: text

text:
# Limits on the WIMP-Nucleon Coupling Coefficients from Dark Matter Search Experiment with NaF Bolometer

## 1 Introduction

Several lines of evidence indicate that the universe contains a large amount of nonbaryonic dark matter. The most plausible candidates for this nonbaryonic dark matter are weakly interacting massive particles (WIMPs) provided by the supersymmetry in the form of the lightest supersymmetric particles (LSPs). LSPs can be directly detected through their elastic scattering off target nuclei in a detector . Many experiments have been performed for the direct detection of WIMPs. Among them, the DAMA experiment has reported an annual modulation signal from the four years’ measurement with nine 10 kg NaI(Tl) detectors . The parameters reported by the DAMA experiment are $`M_{\mathrm{WIMP}}=(52_8^{+10})`$ GeV for the WIMP mass and $`\sigma _{\mathrm{WIMP}\mathrm{p}}^{\mathrm{SI}}=(7.2_{0.9}^{+0.4})\times 10^6`$ pb for the WIMP-proton spin-independently (SI) coupled cross section. On the other hand, the CDMS  and EDELWEISS  experiments have excluded most of the parameter space claimed by the DAMA experiment under the assumptions of standard WIMP interactions and a halo model. Since the cross section for the WIMP-nucleus elastic scattering is given by the sum of SI and spin-dependent (SD) terms, the dark matter search experiments should be performed both via SI and SD interactions. In an SI and SD mixed framework , $`\sigma _{\mathrm{WIMP}\mathrm{p}}^{\mathrm{SD}}`$ could be as high as 0.6 pb for 30–80 GeV WIMPs in the DAMA allowed region under the assumptions of certain parameter sets. Since the limits for SD WIMP-proton cross section, $`\sigma _{\mathrm{WIMP}\mathrm{p}}^{\mathrm{SD}}`$, obtained by the conventional method  are highly WIMP model dependent, both $`\sigma _{\mathrm{WIMP}\mathrm{p}}^{\mathrm{SD}}`$ and SD WIMP-neutron cross section, $`\sigma _{\mathrm{WIMP}\mathrm{n}}^{\mathrm{SD}}`$, limits or the limits in terms of the fundamental WIMP-nucleon coupling coefficients $`a_\mathrm{p}`$ and $`a_\mathrm{n}`$ should be derived .

We, the Tokyo group, reported our previous results with the lithium fluoride (LiF) bolometer at Kamioka Observatory (2700 m.w.e.) , and set the new limits in the $`a_\mathrm{p}`$$`a_\mathrm{n}`$ plane. We compared our results with those of the UKDMC experiments , and tightened the UKDMC’s limits that are close to the DAMA’s SD allowed region in Ref. .

Recently, we developed a sodium fluoride (NaF) bolometer array, and performed a dark matter search measurement from 2002 through 2003. Since $`{}_{}{}^{23}\mathrm{Na}`$ has a larger expectation value of the neutron spins than $`{}_{}{}^{7}\mathrm{Li}`$ as shown in Table 3, we were able to tighten the limits from the previous experiment  in the $`a_\mathrm{p}`$$`a_\mathrm{n}`$ plane.

In this paper, we present the results from a dark matter search experiment performed at Kamioka Observatory using the NaF bolometer array.

## 2 Experimental set-up

The present data were collected at the Kamioka Observatory under a 2700 m.w.e. rock overburden located at Mozumi Mine of the Kamioka Mining and Smelting Co. in Kamioka-cho, Gifu, Japan. Details of the underground laboratory and the experimental set-up are given in Ref. .

The detector is an array of eight NaF bolometers with a total mass of 176 g. Each bolometer consists of $`2\times 2\times 2\mathrm{cm}^3`$ NaF crystal and a neutron transmutation doped (NTD) germanium thermistor glued to it. Each NTD thermistor is biased through a 300 M$`\mathrm{\Omega }`$ load resistor, and the voltage change across the thermistor is fed to the voltage sensitive amplifier with a gain of 10 placed at the 4 K cold stage. The signal from the cold stage amplifier is in turn amplified by the voltage amplifier with a gain of $`4.86\times 10^2`$ placed at the room temperature, and is digitized with a sampling rate of 1 kHz. All the data are recorded continuously without any online trigger for a complete offline analysis.

We refer to each bolometer as from D1 to D8. Five of the NaF crystals are purchased from OHYO KOKEN KOGYO Co. Ltd., two of them are from RARE METALLIC Co. Ltd., and the other is from MATERIAL TECHNOLOGY Co. Ltd. We used the NTD thermistors, which are referred to as NTD3 in Ref. , for the eight bolometers. The concentration of radioactive contaminations in the NaF crystals was previously checked with a low-background germanium spectrometer placed at Kamioka Observatory (2700 m.w.e.), and was found to be less than 0.1 ppb for U, 1 ppb for Th, and 2 ppm for K. In order to eliminate remaining grains of the polishing powder, surfaces of the NaF crystals were etched with ultra-pure water in a clean bench of class 100.

The bolometer array is encapsulated in an inner lead shield with a thickness of 2 cm. In the previous measurement, this inner shield was made of over 200 years old lead. We used a new inner shield with a low-radioactivity lead older than 400 years for this experiment.

## 3 Measurement and analysis

From 23rd December 2002 through 24th January 2003, a dark matter search measurement was performed. In this work, D6 and D7 showed poor signal to noise ratios due to the incomplete thermal conduction between the crystal and the thermistor and due to the high electric noises, respectively. We therefore analyzed the data obtained with the other six bolometers.

We performed numerical pulse shaping to all the data using the optimal filter, or the Wiener filter . The energy of each event was determined by the peak height of a pulse after shaping. The shaped waveform $`U(\omega )`$ in the Fourier space is given by

$`U(\omega )={\displaystyle \frac{S^{}(\omega )C(\omega )}{|N(\omega )|^2}},`$ (1)
where $`S(\omega )`$, $`C(\omega )`$ and $`N(\omega )`$ are Fourier transformations of the ideal signal waveform, the measured waveform and the noise, respectively, and the asterisk denotes complex conjugation. Averaged signal of the gamma-ray source of $`{}_{}{}^{137}\mathrm{Cs}`$ was used as the ideal signal waveform, while the noises were obtained by gathering data while no pulse, whose height above a sufficiently low level, exists. A good lineality of this shaping method was confirmed using simulated pulses which were obtained by adding the ideal signals of known energy to the raw data. Since the signals shaped from the real events have the symmetrical waveform, the asymmetrical waveform from any electric noise events or any events that were thought to have occurred within the thermistors were eliminated by the waveform analysis. The event selection efficiency was estimated by analyzing the simulated pulses in the same manner as the real data because the deficiency is mainly due to the threshold smeared by the baseline fluctuation caused by the noises. Energy resolution of each detector was also estimated by a similar method, and listed in Table 1. Energy threshold of each detector was set at the energy in which the efficiency was equal to 0.8, and the events above this threshold were selected as real bolometer signals. Finally, those events which have associated signals in other bolometers within $`\pm 10`$ ms were discarded because the expected rate of the coincident WIMP events is extremely low.

Calibrations were carried out using the $`{}_{}{}^{60}\mathrm{Co}`$ and $`{}_{}{}^{137}\mathrm{Cs}`$ gamma-ray sources by the same manner as described in Ref. . Obtained gains for six bolometers are listed in Table 1. Gains and linearities were checked every four or five days, and the gains were found to be stable within $`\pm 5\%`$.

Fig. 1 shows the spectra obtained with the six bolometers together with one of the spectra (D6) obtained in the previous measurement with the LiF bolometer . Live times of the six bolometers are listed in Table 1. The total exposure is 3.38 $`\mathrm{kg}\mathrm{days}`$. From this measurement with the NaF bolometer, we obtained the same background levels as the previous measurement.

From the obtained spectra, we derived the limits on the SD WIMP-proton elastic scattering cross section, $`\sigma _{\mathrm{WIMP}\mathrm{p}}^{\mathrm{SD}}`$, in the same manner as used in Refs. . The limits are shown in Fig. 2. The best limit is 27 pb for $`M_{\mathrm{WIMP}}=20\mathrm{GeV}`$. Here we used the astrophysical and nuclear parameters listed in Table 2. We used the $`\lambda ^2J(J+1)`$ values calculated by assuming the odd group shell model , where $`\lambda `$ is $`\mathrm{Land}\stackrel{´}{\mathrm{e}}`$ factor and $`J`$ is the total spin of the nucleus.

As discussed previously, these $`\sigma _{\mathrm{WIMP}\mathrm{p}}^{\mathrm{SD}}`$ limits are fraught with potentially significant WIMP model dependence. We then derived the model-independent limits on $`\sigma _{\mathrm{WIMP}\mathrm{p}}^{\mathrm{SD}}`$ and $`\sigma _{\mathrm{WIMP}\mathrm{n}}^{\mathrm{SD}}`$ in the same procedure as described in Ref. . The SD WIMP-nucleus cross section, $`\sigma _{\mathrm{WIMP}\mathrm{N}}^{\mathrm{SD}}`$, is given by the following equation:

$`\sigma _{\mathrm{WIMP}\mathrm{N}}^{\mathrm{SD}}`$ $`=`$ $`4G_\mathrm{F}^2\mu _{\mathrm{WIMP}\mathrm{N}}^2\times {\displaystyle \frac{8}{\pi }}(a_\mathrm{p}S_{\mathrm{p}(\mathrm{N})}+a_\mathrm{n}S_{\mathrm{n}(\mathrm{N})})^2{\displaystyle \frac{J+1}{J}},`$ (2)
$`=`$ $`4G_\mathrm{F}^2\mu _{\mathrm{WIMP}\mathrm{N}}^2\times \left(\sqrt{C_{\mathrm{p}(\mathrm{N})}^{\mathrm{SD}}}\pm \sqrt{C_{\mathrm{n}(\mathrm{N})}^{\mathrm{SD}}}\right)^2,`$ (3)
where $`G_\mathrm{F}`$ is the Fermi coupling constant, $`\mu _{\mathrm{WIMP}\mathrm{N}}`$ is the WIMP-nucleus reduced mass, $`a_\mathrm{p}`$ and $`a_\mathrm{n}`$ are the WIMP-proton and WIMP-neutron coupling coefficients, respectively, $`S_{\mathrm{p}(\mathrm{N})}`$ and $`S_{\mathrm{n}(\mathrm{N})}`$ are the expectation values of the proton and neutron spins within the nucleus N, respectively, $`C_{\mathrm{p}(\mathrm{N})}^{\mathrm{SD}}`$ and $`C_{\mathrm{n}(\mathrm{N})}^{\mathrm{SD}}`$ are the proton and neutron contributions to the total enhancement factor of nucleus N, respectively. The sign in the square of Eq. (3) is determined by the relative sign of $`a_\mathrm{p}S_{\mathrm{p}(\mathrm{N})}`$ and $`a_\mathrm{n}S_{\mathrm{n}(\mathrm{N})}`$. $`S_{\mathrm{p}(\mathrm{N})}`$ and $`S_{\mathrm{n}(\mathrm{N})}`$ calculated with recent shell models are shown in Table 3. From the WIMP-nucleus cross section limit, $`\sigma _{\mathrm{WIMP}\mathrm{N}}^{\mathrm{SD}\mathrm{lim}}`$, which is directly set by the experiments, the respective limits on WIMP-proton and WIMP-neutron cross section, $`\sigma _{\mathrm{WIMP}\mathrm{p}(\mathrm{N})}^{\mathrm{SD}\mathrm{lim}}`$ or $`\sigma _{\mathrm{WIMP}\mathrm{n}(\mathrm{N})}^{\mathrm{SD}\mathrm{lim}}`$, are given by assuming that all events are due to WIMP-proton and WIMP-neutron elastic scatterings in the nucleus, such as:

$$\sigma _{\mathrm{WIMP}\mathrm{p}(\mathrm{N})}^{\mathrm{SD}\mathrm{lim}}=\sigma _{\mathrm{WIMP}\mathrm{N}}^{\mathrm{SD}\mathrm{lim}}\frac{\mu _{\mathrm{WIMP}\mathrm{p}}^2}{\mu _{\mathrm{WIMP}\mathrm{N}}^2}\left(\frac{C_{\mathrm{p}(\mathrm{N})}^{\mathrm{SD}}}{C_\mathrm{p}^{\mathrm{SD}}}\right)^1,$$
(4)
$$\sigma _{\mathrm{WIMP}\mathrm{n}(\mathrm{N})}^{\mathrm{SD}\mathrm{lim}}=\sigma _{\mathrm{WIMP}\mathrm{N}}^{\mathrm{SD}\mathrm{lim}}\frac{\mu _{\mathrm{WIMP}\mathrm{n}}^2}{\mu _{\mathrm{WIMP}\mathrm{N}}^2}\left(\frac{C_{\mathrm{n}(\mathrm{N})}^{\mathrm{SD}}}{C_\mathrm{n}^{\mathrm{SD}}}\right)^1,$$
(5)
where $`\mu _{\mathrm{WIMP}\mathrm{p}}`$ and $`\mu _{\mathrm{WIMP}\mathrm{n}}`$ are the WIMP-proton and WIMP-neutron reduced mass, respectively, and $`C_\mathrm{p}^{\mathrm{SD}}`$ and $`C_\mathrm{n}^{\mathrm{SD}}`$ are the enhancement factors of proton and neutron themselves, respectively. Ratios of the enhancement factors calculated by the shell models are shown in Table 3.

The obtained model-independent limits on $`\sigma _{\mathrm{WIMP}\mathrm{p}}^{\mathrm{SD}}`$ and $`\sigma _{\mathrm{WIMP}\mathrm{n}}^{\mathrm{SD}}`$ are shown in Fig. 3. We then derived the limits in the $`a_\mathrm{p}`$$`a_\mathrm{n}`$ plane from them. According to Ref. , the allowed region in the $`a_\mathrm{p}`$$`a_\mathrm{n}`$ plane is defined as:

$$\left(\frac{a_\mathrm{p}}{\sqrt{\sigma _{\mathrm{WIMP}\mathrm{p}(\mathrm{N})}^{\mathrm{SD}\mathrm{lim}}}}\pm \frac{a_\mathrm{n}}{\sqrt{\sigma _{\mathrm{WIMP}\mathrm{n}(\mathrm{N})}^{\mathrm{SD}\mathrm{lim}}}}\right)^2<\frac{\pi }{24G_\mathrm{F}^2\mu _{\mathrm{WIMP}\mathrm{p}}^2},$$
(6)
where the relative sign inside the square is determined by the sign of $`S_{\mathrm{p}(\mathrm{N})}/S_{\mathrm{n}(\mathrm{N})}`$, and the small proton-neutron mass difference is ignored. Obtained limits in the $`a_\mathrm{p}`$$`a_\mathrm{n}`$ plane for WIMPs with mass of 10 GeV, 50 GeV and 100 GeV are shown in Fig. 4. The limits on the respective coupling coefficients are, for instance, $`|a_\mathrm{p}|<16`$ and $`|a_\mathrm{n}|<68`$ for $`M_{\mathrm{WIMP}}=50\mathrm{GeV}`$.

In Fig. 4, we compared our results with those of our previous results  and those of UKDMC experiment . Though our $`\sigma _{\mathrm{WIMP}\mathrm{p}}^{\mathrm{SD}}`$ limits derived from this measurement are almost the same as those derived from the previous results with the LiF bolometer, we improved the limits in the $`a_\mathrm{p}`$$`a_\mathrm{n}`$ plane significantly. This is because $`{}_{}{}^{N}23`$ has a larger expectation value of the neutron spins than $`{}_{}{}^{L}7`$. It is clearly seen that our results obtained with LiF or NaF detectors set limits complementary to the results obtained with NaI detectors of UKDMC experiment because the spin factor of $`{}_{}{}^{F}19`$ are large and the sign of $`S_{\mathrm{p}(\mathrm{N})}/S_{\mathrm{n}(\mathrm{N})}`$ for $`{}_{}{}^{F}19`$ is opposite to those for $`{}_{}{}^{L}7`$, $`{}_{}{}^{N}23`$, and $`{}_{}{}^{I}127`$. As described in Ref. , we did not compare our results with the SD annual modulation allowed regions shown in Fig. 2  because the DAMA group had not shown the results from $`{}_{}{}^{N}23`$ and $`{}_{}{}^{I}127`$ independently.

The remaining background in the energy region between 10 and 90 keV restricting the current sensitivity is supposed to be due to the uranium and thorium contamination in a crystal holder made of the OFHC copper. A use of the crystal holder with a sufficiently radio-pure material such as polytetrafluoroethylene (PTFE) and active shields between the crystals and their holder would help to reduce the background from the holder and the other outside materials.

## 4 Conclusion

We have performed the dark matter search experiment at Kamioka observatory with the NaF bolometer. With a total exposure of 3.38 $`\mathrm{kg}\mathrm{days}`$, we have confirmed and tightened the limits in the $`a_\mathrm{p}`$$`a_\mathrm{n}`$ plane obtained by the previous results with the LiF bolometer. For the WIMPs with a mass of 50 GeV, we have successfully excluded more than half of the parameter region in the $`a_\mathrm{p}`$$`a_\mathrm{n}`$ plane allowed by the previous results with the LiF bolometer. Our results set limits complementary to the results obtained by NaI detectors of UKDMC experiment. Since the limits from the UKDMC experiments are close to the DAMA’s SD allowed region in Ref. , we should also be able to compare our limits in the $`a_\mathrm{p}`$$`a_\mathrm{n}`$ plane if DAMA publishes relevant data to calculate these parameters. It is worth noting that this result has been obtained by a single detector with both $`{}_{}{}^{23}\mathrm{Na}`$, which has relatively large $`C_{\mathrm{n}(\mathrm{N})}^{\mathrm{SD}}/C_\mathrm{n}^{\mathrm{SD}}`$, and $`{}_{}{}^{19}\mathrm{F}`$, which has the largest $`C_{\mathrm{p}(\mathrm{N})}^{\mathrm{SD}}/C_\mathrm{p}^{\mathrm{SD}}`$, in the detector material.

## Acknowledgements

We would like to thank all the staff of Kamioka Observatory of the Institute for Cosmic Ray Research, the University of Tokyo, to whose hospitality we owe a great deal in using the facilities of the Observatory. This research is supported by the Grant-in-Aid for COE Research by the Japanese Ministry of Education, Culture, Sports, Science and Technology.


text:
# Scattering of dipole-mode vector solitons: Theory and experiment

## I Introduction

An understanding of the interaction of simple physical objects leading to the formation of more complex objects is an ultimate goal of fundamental research in many fields of physics. Recent progress in generating spatial optical solitons— the self-trapped states of light with particle-like properties— in various nonlinear bulk media allows to study the truly two-dimensional self-trapping of light and different types of interaction of multi-dimensional solitary waves, including the formation of more complicated localized states book .

Spatial optical solitons have attracted considerable attention as possible building blocks of all-optical switching devices where light is used to guide and manipulate light itself book . The robust nature of spatial optical solitons displayed in the propagation and interaction science  allows us to draw an analogy with atomic physics, treating spatial solitons as “atoms of light”. Furthermore, when several light beams generated by coherent sources are combined to produce vector or composite solitons, this process can be viewed as the formation of composite states or “molecules of light”.

Recently, the existence of robust “molecules of light” in the form of dipole-mode vector solitons was predicted theoretically us  and also verified experimentally exper . The dipole-mode solitons (or ‘dipoles’, for simplicity) originate from the trapping of a dipole-mode optical beam by an effective waveguide created by a mutually incoherent fundamental beam of nearly radial symmetry. The first observation of this novel type of optical vector soliton was reported in Ref. exper , where the dipoles were generated using two different methods: the phase imprinting technique and a symmetry-breaking instability of a vortex-mode composite soliton, another type of fundamental (radially symmetric) composite soliton created by the incoherent coupling of two optical beams moti ; scat . It is worth to mention that vector solitons can be also created in certain processes involving coherent interaction of waves such as second harmonic generation. In this particular situation, the two constituent beams forming a soliton molecule are fundamental and its second harmonic, respectively clusters2HG .

The concept of vector solitons as ‘molecules of light’ can be compared with photonic microcavity structures, micrometer-size “photonic quantum dots” that confine photons in such a way that they act like electrons in an atom photonic . When two of these “photonic atoms” are coupled together, they produce a “photonic molecule” whose optical modes bear a strong resemblance to the electronic states in a diatomic molecule like hydrogen photonic2 . The self-trapped states of light we study here can be viewed as somewhat similar photonic structures where, the photonic trap and the beam it guides are both made of light and create self-trapped photonic atoms and molecules.

In this paper we present a comprehensive study of the scattering properties of the dipole-mode vector solitons and analyze, in particular, the interaction between these objects and other self-trapped structures such as scalar optical solitons and other dipoles. We describe a number of interesting effects observed in numerical simulations of such interactions, for both isotropic saturable and anisotropic nonlocal nonlinear models. These include the absorption of a soliton by a dipole and the replacement of the soliton with a dipole component, the field momentum redistribution that can be viewed as the transformation of a linear momentum into an angular momentum with the subsequent dipole spiraling, etc. Additionally, we verify experimentally some of our analytical predictions by studying the generation and scattering of the composite spatial solitons in photorefractive nonlinear crystals. The versatility of the phenomena described here makes dipole-mode vector solitons of great importance not only because of the fundamental interest in nonlinear physics but also because of potential promising applications in all optical switching and integrated optics.

The paper is organized as follows. Section II studies the scattering of scalar solitons and dipole-mode vector solitons in the framework of an isotropic saturable nonlinear medium. In this section, we also present the most important examples of the dipole-soliton interactions. Section III includes the studies of the formation and interaction of the dipole-mode solitons in an anisotropic nonlocal nonlinear model which is used for describing the nonlocal anisotropic nonlinearities of photorefractive media. In Sec. IV, we summarize the results of our experimental studies of the interaction of the dipole solitons in photorefractive nonlinear crystals. Finally, Sec. V concludes the paper.

## II Soliton scattering in a saturable isotropic medium

### II.1 Model and solitons

We consider here the propagation of two light beams interacting incoherently in a saturable nonlinear medium. In the steady state regimes in the paraxial approximation, the mutual beam interaction can be described by a system of two coupled nonlinear Schrödinger (NLS) equations us ; exper ; moti ; spiral ,

$`i{\displaystyle \frac{u}{z}}={\displaystyle \frac{1}{2}}\mathrm{}_{}u+F(I)u,`$ (1a)
$`i{\displaystyle \frac{v}{z}}={\displaystyle \frac{1}{2}}\mathrm{}_{}v+F(I)v,`$ (1b)
where $`u(𝐫_{},z)`$ and $`v(𝐫_{},z)`$ are the dimensionless envelopes of the beams which are self-trapped in the cross-section plane $`𝐫_{}=(x,y)`$ and propagate along the direction $`z`$. The function $`F(I)=I(1+sI)^1`$ characterizes a saturable nonlinearity of the medium, where $`s`$ is a dimensionless saturation parameter $`(0<s<1)`$ and $`I=|u|^2+|v|^2`$ is the total beam intensity.

We would like to mention that we consider here only the stationary propagation of light excluding any nonlinearity-mediated temporal effects. In fact, the experimental observations indicate that this is a common situation in many nonlinear systems involving free (no feedback) propagation of optical beams in both fast (such as atomic) and slow (photorefractive, thermal) nonlinear media. Typically, the only dynamics which may occur in such cases is related to switching effects and dies out within the time scale determined by the temporal response of the medium. However, the temporal in noninstantaneous media such as photorefractive crystals may lead to a number of novel transiting effects (see, e.g., Refs. temp ; temp0 ; temp1 ; temp2 ).

Equations (1) describe different types of spatially localized composite solutions. The dipole-mode vector soliton (or “a molecule of light”) is a stationary state which is composed of a nodeless beam in the $`v`$ component and a dipole beam (or a pair of out-of-phase solitons) in the $`u`$ component. Solitons in the $`u`$ component have opposite phases and thus they repel each other, but the role of the complimentary beam $`v`$ is to stabilize the structure making it robust. A numerical analysis of the linearized equations (1) shows no signs of linear instability of this composite structure us , as was also recently confirmed by means of the asymptotic analytical theory yang . Moreover, it was shown us  that such robust dipole-mode vector solitons exist for a wide range of the beam powers $`P_u=|u|^2𝑑𝐫_{}`$ and $`P_v=|v|^2𝐫_{}`$. Since we are interested in showing stability far from the regime in which one beam is dominant, all numerical experiments are performed using as initial conditions stationary states in which $`P_uP_v`$.

We are interested in the dynamics of the dipole soliton under the action of finite external perturbations introduced by its collision with other objects. The word “finite” emphasizes the fact that we can no longer make use of linearized equations and that we must deal with the full system (1). This fact, combined with the complex structure of the dipole which lacks radial symmetry, makes analytical predictions on the dipole dynamics very difficult. Nevertheless, as will be shown below, one may extract some general rules on which qualitative predictions may be based.

The idea is that the dipole can be seen as a bound state of a soliton beam (in $`v`$) plus a pair of vortices with opposite charges (in $`u`$) and, therefore, many of the effects observed in the composite beam collisions and described below can be understood once the mutual interaction of these simpler objects is studied.

One of the components of the dipole is a soliton beam (to be referred to as soliton hereafter). Spatial solitons are stable localized states which have no nodes and which are the states of minimum energy of the system for a fixed power. When two of these solitons are in different, mutually incoherent beams (say, one in $`u`$ and the other one, in $`v`$), they interact incoherently and attract each other. Thus, during an incoherent interaction two solitons may either become bound or scatter. In the former case, we have an example of what we call a molecule of light, which is typically referred to as a “vector soliton”. However, when two solitons are derived from the same beam they interact coherently and the outcome of their mutual interaction depends on their phase difference. When this quantity is small or zero, solitons experience mutual attraction, whereas if their mutual phases differ by $`\pi `$, they repel each other.

Another nonlinear structure that should be mentioned in this context is a vortex-mode composite soliton introduced in moti  which in our model (1) is stable only in the vicinity of the bifurcation point yang . Thus, the vortices may only be stabilized by co-propagating with a very large soliton beam (e.g., when a vortex in the linear beam $`u`$ is guided by an effective waveguide created in the component $`v`$). Otherwise, a composite state of a vortex plus a soliton constitutes an unstable molecule of light.

A dipole can be seen as a pair of vortices as described above or, alternatively, as a bound state of two solitons with a phase difference of $`\pi `$. While in principle, these solitons should repel each other, the system is stabilized due to the interaction with a soliton-induced waveguide created by the other, mutually incoherent, component.

### II.2 Numerical results for the soliton collisions

#### II.2.1 Soliton-dipole scattering

The first type of numerical simulations we present here consists in shooting a scalar soliton against a dipole-mode vector soliton. All the simulations discussed here have been performed using a split-step operator technique using FFT, with grid sizes of up to 512$`\times `$512 points covering a rectangular domain of 68$`\times `$34 adimensional units. The initial data are always a combination of stationary states. For instance, when a soliton is launched against a dipole, we start with

$`u(𝐱,0)`$ $`=`$ $`u_{dipole}(𝐱)+u_{soliton}(𝐱𝐝)e^{i𝐩_\mathrm{𝟎}𝐱},`$ (2)
$`v(𝐱,0)`$ $`=`$ $`v_{dipole}(𝐱).`$ (3)
Here $`𝐝=(d_x,d_y)`$, $`d_xd_y`$, $`d_y`$ is the impact parameter, and $`𝐩_\mathrm{𝟎}`$ is proportional to the initial (linear) momentum of the incoming scalar soliton. The initial data $`u_{dipole}`$, $`u_{soliton}`$, and $`v_{dipole}`$ are obtained numerically by a suitable minimization procedure outlined in Ref. us .

The result is an inelastic collision in which the soliton becomes deflected and the dipole gains both linear and angular momenta. The whole process is depicted in Fig.1. Soliton scattering occurs when the incident beam has medium to large linear momentum or when it has an appropriate initial phase. For instance, in Fig.1 the incident soliton has sign $`()`$ and it crashes against the part of the dipole with $`(+)`$ sign. A conservation law forces the dipole to rotate and the soliton becomes deflected, sometimes as much as by a 90 degree angle.

When the linear momentum of the incident soliton is large, it moves too fast to suffer a destructive influence from the dipole. In Fig.2 we plot the exchange of the linear momentum between the soliton and dipole as a function of the impact parameter. The effective interaction is clearly attractive: the soliton coming from below ($`d_y<0`$) feels the drag of the dipole above it and gets deflected upwards ($`p_y>0`$), while the dipole moves downwards.

The second family of numerical experiments is performed with solitons which are slow and, as is usual in scattering processes, the effects of the interaction process may be more drastic. For some impact parameters the soliton gets too close to the lobe of the dipole with the smallest phase difference and fuses with it with some emission of radiation and a subsequent rotation of the dipole. This is well reflected in Fig.3 (radiation is not seen).

#### II.2.2 Dipole-dipole collisions

The third family of numerical simulations corresponds to shooting dipoles against each other. These collisions, which resemble interaction of atomic molecules provide a rich source of phenomena depending on the mutual orientation of the dipoles and on the initial energy. Figure 4 summarizes the main results observed. There we see three cases (a-c) in which the dipole solitons are preserved. The figure shows an in-phase collision with weak interaction \[Fig.4(a)\], an out-of-phase collision with repulsion \[Fig.4(b)\], and an example of the collision with nonzero impact parameter in which two vortex states are created and they decay into a pair of spiralling solitons \[Fig.4(c)\].

The last case, Fig.4 (d), shows an interesting inelastic process when two dipoles fuse into a more complex state which then decays creating a new dipole and a pair of simple solitons. All these processes may be understood in terms of the phase of the lobes of each dipole as described above.

## III Soliton scattering in anisotropic nonlocal media

### III.1 Composite solitons

The dipole-mode vector solitons considered so far were restricted to those realized in isotropic nonlinear media. However, up to now the majority of experimental observations of dipole-mode and multi-pole vector solitons have been performed in photorefractive nonlinear crystals which are known to exhibit anisotropy in their nonlinear response anisotropy . In effect, even circularly symmetric optical beams induce strongly asymmetric refractive index changes which significantly affect the formation of spatial solitons as well as their interaction.

In this section we employ the commonly accepted model for the photorefractive nonlinearity that takes into account its most important properties wieslaw  to investigate some of the previously discussed examples of interactions of vector solitons.

The interactions we consider here involve dipole-mode vector solitons and scalar solitons. The dipole-mode vector solitons consist of two mutually incoherent optical beams with the envelopes $`u`$ and $`v`$, propagating in a bulk anisotropic nonlocal nonlinear medium such as a biased photorefractive crystal. When the characteristic spatial scales are larger than the photorefractive Debye length and the diffusion field may be neglected, the steady-state propagation along the $`z`$ axis of a photorefractive crystal with an externally applied electric field along the $`x`$ axis is described by the equations:

$$\begin{array}{c}i\frac{u}{z}+\frac{1}{2}^2u=\frac{\gamma }{2}\frac{\phi }{x}u,\hfill \\ i\frac{v}{z}+\frac{1}{2}^2v=\frac{\gamma }{2}\frac{\phi }{x}v,\hfill \\ ^2\phi +\phi \mathrm{ln}(1+I)=E_0x_0\frac{}{x}\mathrm{ln}(1+I),\hfill \end{array}$$
(4)
where $`\gamma `$ and $`E_0`$ are the normalized nonlinearity coefficient and external field, respectively, $`I|u|^2+|v|^2`$ is the total intensity, $`=\widehat{x}(/x)+\widehat{y}(/y)`$, and $`\phi `$ is the dimensionless electrostatic potential induced by the light with the boundary condition $`\phi (\stackrel{}{r}\mathrm{})0`$. The propagation coordinate $`z`$ is measured in units of the diffraction length, and the transverse coordinates are normalized by the characteristic beam size $`x_0`$. The above system of equations was solved numerically by applying concurrently the finite difference and split-step Fast Fourier methods to the electrostatic potential and propagation equations, respectively and using a 256x256 grid size computational window.

### III.2 Soliton-dipole scattering

The first two of the numerical simulations we present here consist of colliding a scalar soliton into one lobe of a dipole-mode vector soliton. The scalar soliton is coherent with the dipole component of the vector soliton and can thus be considered part of the same beam in the theoretical model (4) above. The initial data for all cases presented in this section reflect the experimental conditions in that they are not exact solutions to the model. The fundamental component of the vector soliton and the scalar beam are Gaussian beams and the dipole component is the first derivative of a Gaussian beam along $`y`$ coordinate.

The first example of the collision between the dipole-mode vector soliton and the scalar soliton is shown in Fig.5. Here the scalar soliton is in phase with the lobe of the dipole with which it collides. This leads to a strong attraction between the scalar beam and the lobe and eventual absorption of the former followed by, as in the isotropic case, a rotation of the whole dipole-mode vector soliton. In a photorefractive crystal the complete rotation of a dipole soliton is prohibited by the anisotropy of the nonlinear refractive index change optlet . Hence, unlike the isotropic case (see Fig.1), the vector soliton exhibits only angular oscillations about the vertical axis. As Fig.5 clearly shows the collision leads to dipole soliton experiencing also a lateral shift which is due to transfer of a linear momentum from the scalar soliton to the vector soliton. As our numerical simulations show the outcome of this collision can be more dramatic if the intensity of the scalar beam is increased. In such a case the collision may result in the break up of the vector soliton such that the two out-of-phase lobes of the dipole are no longer trapped. Such an effect occurs when the intensity of the scalar beam is comparable to that of the fundamental component of the vector soliton.

In the next example shown in Fig.6 the scalar soliton and the lobe of the dipole it interacts with are out-of-phase. Both solitons propagate initially along parallel trajectories. Because of the phase relation the interaction is now repulsive leading to rotation of the dipole. Again the rotation is limited by the anisotropic refractive index distribution. On further propagation, the vector and scalar solitons are clearly repelled and the vector soliton will reorientate to the stable vertical direction once the effect due to presence of the scalar beam is negligible. Our simulations show that by increasing the angle between initial trajectories of both solitons one can induce even stronger rotation of the dipole but this constitutes so drastic perturbation to the dipole that it often breaks-up so that the dipole lobes are no longer trapped.

### III.3 Dipole-dipole scattering

The next few examples of the numerical simulations involve the collisions of two dipole-mode vector solitons. Unless stated otherwise, the dipole components of the vector solitons are coherent with each other while the fundamental components are mutually incoherent.

In the example depicted in Fig.7 the two identically oriented vector solitons collide centrally propagating along the direction of the applied electric field. The mutual interaction is now attractive (phases of dipoles coincide) and since the intersection angle (or transverse velocities) of the solitons is relatively small the solitons fuse forming a new dipole-mode vector soliton. This fused soliton undergoes internal oscillations as its parameters are quite different from the exact soliton solution. Further simulations show that after emitting some radiation the structure reaches a steady state.

In Fig.8 the phase of one of the dipoles has been shifted by $`\pi `$ and hence each lobe now collides with an out of phase lobe of the other dipole. This leads to the mutual repulsion of the dipole-mode vector solitons, which the weak attraction of the incoherent Gaussian beams is unable to counteract. Further simulations show that the structure of each soliton is well preserved throughout the interaction and therefore this collision could be considered as an almost elastic collision. Note very close similarity between this result and simulations of dipoles interaction in isotropic medium shown in Fig.4.

If the initial trajectories of both solitons are chosen such that the collision is non-central then both dipoles usually undergo spatial rotation as shown in Fig.9. This rotation is initiated by the mutual repulsion of the out-of-phase lobes in respective dipole components. Again, the rotation is limited by the anisotropic induced refractive index change as the dipoles force their way past each other. As the solitons pass the zone of interaction they undergo oscillatory rotation about the vertical axis which is again due to the orientation perturbation caused by the collision as can be seen in (d),(e) and (f).

Fig.10 shows another example of the noncentral collision of the two dipoles for the case when the initial angle between soliton trajectories is increased to 0.21 degree. As the impact parameter (i.e. lateral mismatch of initial trajectories) is smaller than in the preceding case, the repulsive interaction of dipole components is weaker. This, combined with the larger mutual velocity of the solitons (larger intersection angle) allows the solitons to pass through one another with a weak, but still visible, intermediate tilt of their axes. The final state corresponds to the perturbed dipole states with an excited internal oscillation which is reflected in an unequal energy distribution between both lobes in each dipole.

Finally, Fig.11 shows an example of two dipole-mode vector solitons colliding along their long (y) axes. In this particular example the fundamental components of both dipole solitons are coherent and in-phase. The initial orientation of the dipoles was such that the two directly interacting lobes were out-of-phase hence their interaction was repulsive. On the other hand the interaction of fundamental components is strongly attractive. Since the mutual velocity of the solitons is less (smaller intersection angle) the fundamental components fuse in collision trapping two out-of-phase lobes. In this way new dipole mode vector solitons was formed with its dipole component being constituted by two out-of-phase lobes coming from two different solitons. Notice that because of symmetry of the problem this new vector soliton does not exhibit any transversal motion. The remaining the outmost lobes are no longer bound and separate as fundamental solitary beams as can be seen in Fig.11(e).

## IV Experimental results

The generation of an isolated dipole-mode vector soliton was reported earlier in Ref. exper ; anton ; motzek . Fig.12 shows schematically experimental set-up used in those and current studies.

The dipole-mode soliton can be created using a few different processes: (i) phase imprinting, when one of the beam components is sent through a phase mask in order to imprint the required phase structure, and (ii) symmetry-breaking instability of a vortex-mode composite soliton or (iii) superposition of two oppositely charged vortex beams. In this way, one obtains a dipole-like structure with a phase jump along its transverse direction that is perpendicular to the optical axis of the crystal \[see Fig.13(a)\]. This dipole-like beam is then combined with the second, node-less beam and the resulting composite beam is focused on the input face of the photorefractive SBN crystal (the crystal has the same parameters as in exper ), biased with a DC field of 1-2.5 kV applied along its optical axis. The exit and input facets of the crystal are viewed by CCD cameras and images stored in the computer. To control the degree of saturation of the nonlinearity, we illuminated the crystal with a wide beam derived from a white light source. In our experiments the ratio of the peak intensity of the soliton beam to intensity of white light illumination was always of the order of unity. To ensure that both beams forming a composite structure are mutually incoherent, one of them is reflected by a mirror mounted on a piezo-electric transducer (PZT). When the transducer is driven by a 1kHz signal it introduces a phase modulation into the beam. As the photorefractive crystal is slow, it cannot follow fast changes of the relative phase of both beams making them effectively incoherent. Propagating in a self-focusing saturable medium, such a composite input beam creates a dipole-mode vector soliton, as shown in Fig.13(b-d). As discussed above, the fundamental component creates an effective asymmetric waveguide that guides a dipole-like mode in the form of two out-of-phase solitary beams that otherwise would repel and fly apart.

To observe the soliton-dipole interaction effects, we launched a scalar soliton beam against the dipole soliton. The input state is shown in Fig.13(e), where the dipole-mode soliton is presented by its two-lobe $`u`$ component only. The relative phase of the soliton and the upper lobe of the dipole is close to $`\pi `$. Therefore, when the soliton interacts with a dipole, it gets deflected (repulsive interaction) and transforms a part of its linear momentum into an angular momentum of the dipole that starts rotating clockwise which is clearly visible in Fig.13 (f-h). When the scalar soliton is removed, the vector soliton rotates back until it realigns to its original stable orientation orientation which is vertical in the experimental situation discussed here. Fig.14 illustrates the dynamics of this process showing the intensity of both, dipole (top row) and fundamental (bottom row) components. It should be noted that the above experimental observations are in accordance with the results of numerical modelling. In particular, the intensity distributions shown in Fig.13(e-h) closely resemble theoretical results displayed in Fig.1(a) and Fig.6(c) Also, the observed dynamical behavior of the dipole agrees well with our (not shown here) numerical simulations.

In the subsequent graph (Fig.15) we show the interaction of two closely spaced dipole mode vector solitons. These solitons were propagating in parallel. Graphs Fig.15(a-b) shows the intensity of dipole and fundamental components at the input face of the crystal. Fig.15(c) depicts the dipole components of both solitons seen at the exit of the crystal. This picture was obtained by superimposing two graphs corresponding to different solitons, hence it displays propagation without interaction. The dipole-soliton interaction is shown in the bottom row of this figure. Graphs in Fig.15(d-f) show the total intensity (d) and the intensity of dipole (e) and fundamental components of the interacting dipole-mode solitons. Because of $`\pi `$ phase difference between nearest lobes the interaction leads to a repulsion of the corresponding lobes and the clockwise rotation of both vector solitons. This behavior is analogous to that found in our numerical simulations shown in Fig.9(b-c). In both, theory and experiment, the mutual rotation of the interacting dipoles is evident.

## V Conclusions

We have studied interactions of dipole-mode composite solitons with different nonlinear localized structures such as scalar solitons and other dipole-mode solitons. Our studies demonstrate that, apart from the robustness of the dipole-mode solitons against strong perturbations such as the interaction with other solitons (which is a consequence of their linear stability predicted earlier), in many cases the dynamics of the dipole-mode composite solitons can be understood qualitatively as that of the bound states of simpler solitons. Such dynamics involves two major degrees of freedom of the composite soliton, namely, the rotation of the structure as a whole and the relative oscillation of the lobes of the dipole, which both can be excited in collisions. In some cases, the excitation of the dipole soliton is so strong that the dynamics of the composite structure resembles that of a pair of spiraling beams analyzed earlier in Ref. spiral . This is only one of many interesting phenomena observed in our numerical simulations of dipole-mode soliton scattering which also include the excitation of rotational motion by collision with a scalar soliton, annihilation or strong deflection of the incident soliton, etc. Even richer effects are observed when two dipole-mode solitons collide with a nonzero impact parameter. It is remarkable that the rich dynamics observed here may be understood qualitatively in terms of the structure of the colliding objects and the relative phases of the dipole components. To make our results more realistic providing a comparison with experiment, we have extended our analysis and have studied the anisotropic nonlocal model that is more relevant for describing photorefractive nonlinearities. Finally, we have verified some of our theoretical predictions experimentally employing the self-trapping effect in nonlinear photorefractive crystals.

A rich variety of the scattering effects described theoretically and verified experimentally might make the dipole-mode solitons attractive candidates for realistic applications in the field of integrated optics.

## Acknowledgements

This work was partially supported by the Planning and Performance Fund of the Australian National University, the Australian Research Council, the grant BFM2000-0521 (Ministerio de Ciencia y Tecnolog a) and the grant PAC02-002 (Consejeria de Ciencia y Tecnologia de la Junta de Comunidades de Castilla-La Mancha), Australian-German Joint Reserach Co-operation Scheme, Australian-German project exchange program by DAAD, and the DFG-Graduiertenkolleg ”Nichtlineare kontinuierliche Systeme”.

Figure Captions

* Fig.1. Soliton-dipole scattering. (a) Snapshots of the intensity profile of each of the beams taken at few points along their propagation distance. (b) Three-dimensional plot of the total intensity $`|u|^2+|v|^2`$, which shows the dipole rotation induced in the collision. (c) Same as (b), but with the $`u`$ and $`v`$ components shown separately.
* Fig.2 (a, b) Components of the linear momentum of the incident soliton (solid line) and dipole (marked by circles) after an inelastic collision with a large incident momentum, $`𝐩_𝐮u^{}ud𝐫_{}`$, as a function of the impact parameter $`d_y`$, which shows the crucial role of the dipole asymmetry. Total $`p_y`$ does not vanish because of radiation (not seen in the figure).
* Fig.3 Absorption of a soliton by a dipole. (a) Intensity profile of each of the beams at various value of their propagation distance –the darker the more intense; (b) Three-dimensional plot of the total beam intensity; (c) Same as in (b) with the $`u`$ and $`v`$ components separated.
* Fig.4 Collisions of two dipoles with the zero \[(a), (b), and (d)\] and nonzero \[(c)\] impact parameter, and different orientation of the dipole prior the collision \[cf. (a) and (d)\].
* Fig.5 Soliton absorption in the anisotropic nonlocal model. The top row shows the dipole components of the vector solitons while the bottom row shows the Gaussian components when a colliding soliton is absorbed by the dipole. Relative intensities: Gaussian=1.4, Dipole=1.1, Coherent beam=1.0. Frames are taken at $`z=`$0.0, 1.1, 2.6, 4.1, 5.6 diffraction lengths, respectively. Collision angle=$`0.24^o`$. The applied field is in the horizontal direction.
* Fig.6 Soliton-dipole interaction in the anisotropic model. The top row shows (a) initial dipole components, (b) stable orientation of the dipole components with no interacting scalar soliton at $`z=`$2.3, and (c) rotated dipole when in the presence of a scalar soliton also at $`z=`$2.3. The scalar soliton is of out of phase with the dipole lobe closest to it. The bottom row shows the corresponding fundamental components. The applied field is in the horizontal direction. Relative intensities: fundamental component =1.2, dipole=1.1, scalar beam=1.0.
* Fig.7 Dipole fusion in the anisotropic model. The top row shows the dipole components of the vector solitons while the bottom row shows the Gaussian components. Relative intensities: fundamental=1.0, dipoles=1.8. Frames are taken at $`z=`$0.0, 3.4, 4.8, 5.9, 7.0, respectively. Collision angle=$`0.13^o`$.
* Fig.8 Elastic collision of two dipoles in the anisotropic model. The top row shows the dipole components of the vector solitons, while the bottom row shows the Gaussian components. Relative intensities: Gaussian=1.0, Dipole=1.8. Frames are taken at $`z=`$0.0, 1.5, 3.0, 4.5, 6.0, respectively. Collision angle=$`0.21^o`$.
* Fig.9 Dipole rotation in the anisotropic model. The top row shows the dipole components of the vector solitons while the bottom row shows the Gaussian components. Relative peak intensities: Gaussian=1.0, Dipole=1.2. Frames are taken at $`z=`$0.0, 2.3, 3.4, 4.5, 5.6, 6.7, respectively. Collision angle=$`0.19^o`$.
* Fig.10 The same as in Fig.9 but for the collision angle =$`0.21^o`$. Frames are taken at $`z=`$0.0, 2.3, 3.0, 3.7, 4.4, respectively.
* Fig.11 Dipole collision in the anisotropic model. The top row shows the dipole components of the vector solitons while the bottom row shows the fundamental components. Relative intensities: Gaussian=0.65, Dipole=1.1. Frames are taken at $`z=`$0.0, 1.5, 3.7, 5.2, 6.7, respectively. Collision angle=$`0.17^o`$.
* Fig.12 Experimental setup. PM-phase mask, F-filter, $`\lambda /2`$ \- half-waveplate, P-polarizer, PBS-polarizing beam splitter, PZTM-mirror mounted on piezo-electric transducer, O- microscope objective, F-filter, V-DC biasing field applied to the crystal, CCD - camera.
* Fig.13 Experimental results. Top row: the formation of the dipole soliton. (a) Initial intensity of the dipole component; (b) total intensity of the vector soliton after 10mm propagation in a biased SBN crystal; (c,d) the dipole and fundamental components of the vector soliton after propagation. Bottom row: interaction of dipole and scalar solitons. (e) Initial intensity of the dipole and soliton beams; (f) the same after interaction in a biased SBN crystal; (g-h) dipole and fundamental components of the rotated vector soliton after the interaction. Voltage: 1.3 kV. Intensities: dipole- 1.5$`\mu W`$, the fundamental component - 0.9 $`\mu W`$, scalar soliton - 1.6$`\mu W`$.
* Fig.14 Temporal evolution of the dipole component of the tilted vector soliton after the scalar soliton is blocked. (a) initial intensity of the scalar soliton and vector solitons; (b) dipole component immediately after the scalar beam is blocked; (c-g) temporal evolution of the dipole component. The time step between frames is 0.5sec.
* Fig.15 Experimental results for the interaction of two dipoles. Top row: (a-b) intensity distribution of the two dipoles and two fundamental components; (c)superimposed images of the dipole components of two independently propagating vector solitons (no interaction); Bottom row: soliton interaction. (d) total intensity distribution of two interacting solitons; (e-f) their dipole and fundamental components. Voltage: V=1.1kV. Intensities of both dipoles: 1.5, 1.6 $`\mu W`$ Intensities of the fundamental components: 1.6, 1.7$`\mu W`$.


text:
# References

On some erroneous comments on the literature of neutrino oscillations in the website ‘Neutrino Unbound’ of C.Giunti

J.H.Field

Département de Physique Nucléaire et Corpusculaire, Université de Genève

24, quai Ernest-Ansermet CH-1211Genève 4.

E-mail: john.field@cern.ch

## Abstract

A number of misleading or incorrect comments by C.Giunti on seven arXiv preprints that I have written on the theory of neutrino oscillations are discussed. The essential new features of my approach are also briefly reviewed

PACS 03.65.Ta, 14.60.Pq, 13.20.Cz
Keywords ; Quantum Mechanics, Neutrino Oscillations.

In providing and maintaining the website ‘Neutrino Unbound’  C.Giuni is providing an enormous service to the neutrino physics community. In the comments on the literature in the ‘Theory of Neutrino Oscillations’ section, however, there is a mixture of personal scientific judgements on a limited sub-set of subjectively ‘important’ issues as well as a systematic omission of any comment on certain key contradictory arguments presented in the literature. There are also factual errors in the comments. In an e-mail exchange I had with Giunti on the subject of his webpage he claimed that, in any case, only his personal opinion was expressed. However language is sufficiently rich to discriminate clearly between a fact that one knows (e.g. that such-and-such a subject was discussed in such-and-such a paper leading to such-and-such a conclusion) and one’s personal judgement of the correctness, or not, of some scientific theory. In the following I try my best to make this discrimination.

In the papers I have written recently on the subject of the basic theory of neutrino oscillations  I make different initial physical assumptions (the most important one of which I rediscuss at the end of this note) to those made by essentially all authors on the subject from Pontecorvo onwards. It is then not surprising that I obtain different results. These different basic assumptions that underlie my calculations are not pointed out in any of Giunti’s comments on my papers. Instead it is wrongly claimed in several different places that there are ‘mistakes’ in my calculations. Giunti also claims in his comment on one of his own papers that the correctness of the standard formula for the oscillation phase has been ‘shown’, when, indeed, it has not. Only different initial assumptions are made to those in my papers and different results are obtained.

I have shown in a recent paper  that a clear and simple experimental test exists, comparing oscillations of neutrinos originating either from pion or kaon decays, to discriminate unambigously between the ‘standard’ formula for the oscillation phase  and my predictions.

In the following, I copy verbatim some of the comments on the literature from ‘Neutrino Unbound’ and add my remarks to them.

The description of neutrino and muon oscillations by interfering amplitudes of classical space-time paths, Field, J. H., hep-ph/0110064. Comment: Again the factor of 2 mistake in the phase! This error has been explained in \[Carlo Giunti, Chung W. Kim, Found. Phys. Lett. 14 (2001) 213-229, hep-ph/0011074 \]. \[C.G.\].

The factor 2 difference in the contribution of neutrino propagation to the oscillation phase results from the constraint of space-time geometry on the neutrino trajectories. Since the different mass eigenstates are detected at a common space-time point, travel the same distance and have different velocities, they are necessarily produced at different times in the interfering amplitudes. This assumption was not discussed in the paper of Giunti and Kim cited above. The statement ‘This error has been explained in…’ is therefore false. In fact, it is simply assumed in the paper of Giunti and Kim (in Eqn(22)) that the neutrinos have equal velocities in space-time (i.e. a common production time) but different kinematical velocities $`p/E`$. As explained in much detail in Reference  these manifestly incompatible assumptions (equal space-time velocities but different kinematical velocities) underestimate by a factor of 2 the contribution of neutrino propagation to the oscillation phase and result in the neglect of an important contribution to the oscillation phase from the propagator of the source particle. It is also shown in Reference  that the same initial assumptions, when made in the standard Gaussian wave-packet treatment, leads to the same underestimation of the oscillation phase.

A covariant path amplitude description of flavour oscillations: The Gribov-Pontecorvo phase for neutrino vacuum propagation is right, Field, J. H., hep-ph/0211199. Comment: Same as \[Field, J. H., hep-ph/0110064\] and \[Field, J. H., hep-ph/0110066\]. Mistake explained in \[Carlo Giunti, Chung W. Kim, Found. Phys. Lett. 14 (2001) 213-229, hep-ph/0011074\] and Giunti, C., Physica Scripta 67 (2003) 29-33, hep-ph/0202063\]. The attribution of the ”factor of two” mistake to Gribov and Pontecorvo is an historical aberration. The claim that the ”factor of two” discrepancy in the Gribov and Pontecorvo paper \[Gribov, V. N., Pontecorvo, B., Phys. Lett. B28 (1969) 493.\] was unnoticed before \[ Field, J. H., hep-ph/0211199\] is pure fantasy. The fact is that nobody speculated about it. \[C.G.\].

There is no ‘mistake’, unless the assumption that the space-time and kinematical velocities are the same is one. In any case, nothing is ‘explained’ in the references cited. The key assumption, that is the basis of my calculations, is not even discussed in the first of them. In the second one, it is discussed, and it is claimed that an analagous analysis of the optical Young double slit experiment excludes the possiblity of different production times. However the path amplitude analysis presented is incorrect as it does not take into account the contribution to the interference phase of the propagator of the source. All this is explained in Reference .

Maybe lots of people knew about the factor 2 difference in Gribov and Pontecorvo, but, as I am not a clairvoyant, there is no way I could know this. All I said was that, as far as I know, the difference was never discussed in the literature.

The phase of neutrino oscillations, Giunti, C., Physica Scripta 67 (2003) 29-33, hep-ph/0202063. Comment: It is shown that the standard phase of neutrino oscillations is correct, refuting the claims of a factor of two correction presented in \[ S. De Leo, G. Ducati, P. Rotelli, Mod. Phys. Lett. A15 (2000) 2057-2068, hep-ph/9906460, Field, J. H., hep-ph/0110064., Field, J. H., hep-ph/0110066.\]. The wave packet treatment of neutrino oscillations presented in \[C. Giunti, C. W. Kim, U. W. Lee, Phys. Rev. D44 (1991) 3635-3640, C. Giunti, C. W. Kim, Phys. Rev. D58 (1998) 017301, hep-ph/9711363.\] is improved taking into account explicitly the finite coherence time of the detection process.

The argument given in this paper to refute the results obtained in my papers, as cited above, is based on an incorrect path-amplitude analysis of the Young double slit experiment in optics. This analysis is discussed in Reference  . In Giunti’s comment on this paper (see above) there is no mention that his analysis has been shown in it to be incorrect, and that the correct one gives the same result as the classical wave theory of light. In any case, it is not possible ‘That it is shown that the standard phase of neutrino oscillations is correct’. One can only discuss the correctness, or not, of the mathematical derivation given some initial set of physical assumptions. To prove the correctness of the result, that of the initial asumptions must also be demonstrated. In all papers, authored or co-authored by Giunti, this initial assumption is that the different mass eigenstates have different kinematical velocities but the same space-time velocity. This assumption has never been proved. My personal opinion is that the assumption is incorrect, but I cannot prove this. The experiment I suggested in Reference  can, however, resolve, unambigously, this question.

Coherence in Neutrino Interactions, C. Giunti, hep-ph/0302045. From the abstract: The claim in \[Field, J. H., hep-ph/0301231.\] is refuted in a pedagogical way.

As shown in Reference  the modification proposed in this paper to the Standard Model amplitude for pion decay avoids the experimental constraint of Reference  demonstrating the incoherence of the production of different mass eigenstates in pion decay, but also forbids the possibility of neutrino oscillations following pion decay! The derivation of the standard oscillation phase requires, instead, coherent production of the different mass eigenstates. Since neutrino oscillations following pion decay apparently exist, the predictions of this paper are experimentally excluded, so that it cannot, in any way, ‘refute’ the results of Reference . Giunti’s only comment on Reference  is a quotation from the Bible, the meaning and relevance of which is unclear to me.

The crucial assumption that distinguishes my treatment of the quantum mechanics of neutrino oscillations from all of the papers cited in the section ‘Theory of Neutrino Oscillations’ of ‘Neutrino Unbound’ (with the exception of Reference ) is that mentioned above in connection with Giunti’s comments on my paper hep-ph/0110064. (now superseded by hep-ph/0303151). This is the condition that the velocities of the neutrino mass eigenstates are described consistently both kinematically, and in space-time. That is to say that the neutrinos propagate over macroscopic distances as classical particles. Such an assumption is fundamental to Feynman’s path amplitude formulation of quantum mechanics. To each such classical trajectory corresponds a probability amplitude. That the Feynman path amplitude description gives a correct quantitative description of nature has been demonstrated in many atomic physics experiments involving space-time interference effects, some of which are discussed in Reference .

The universal assumption that is made when the ‘standard’ oscillation phase is derived is that the different neutrino eigenstates are produced at a unique time since they are assumed to be different components of a coherent ‘lepton flavour eigenstate’. As first pointed out by Shrock , in the Standard Model the different mass eigenstates are actually produced incoherently in different physical processes, and so not necessarily at the same time, which is the case when a coherent state is produced. It is shown in Reference  that the production of such a coherent state in pion decay is excluded by the measured $`e/\mu `$ branching ratio.

The other important point where my approach differs from many of the other treatments in the literature is in the importance accorded to the introduction of ad hoc Gaussian ‘wave-packets’ in the description of neutrino oscillations. An in-depth criticicism of such approaches is presented in Reference . In a later paper on the same subject by Giunti  not only were no counter arguments to my criticisms given, but my paper was not even cited! As discussed in Reference  there are indeed momentum wave packets in the exact path amplitude description of pion decay due to the unobserved smearing of the physical mass of the recoiling muon. The effect on the formula for the oscillation phase is, however, completely negligible.

Another subject that has caused much ink to flow and contradictory positions to be adopted is the role of different assumptions concerning the calculation of the kinematical velocities of the neutrinos: exact energy-momentum conservation, equal momenta or equal energies. It is pointed out in Reference  that the same result, at O($`m^2`$), is obtained whatever assumption is made, so that unless O($`m^4`$) corrections to the oscillation phase are of interest, this debate is an irrelevant one.