text:
# QUANTUM CHAOS IN THE HEAVY QUARKONIA

## I Introduction

Level statistics, i.e. the statistical theory describing the fluctuation properties of the quantum spectra was originally developed in nuclear physics in early fifties. Pioneering work dates back to Wigner Wig51  who gave the first detailed theoretical analysis of the nuclear spectra from the viewpoint of random matrix theory. Further developments of this subject were carried out by Porter Porter , Mehta Mehta  and Dyson Dyson . In particular, it has been established that the spectra of the complex quantum many-body systems have a behaviour similar to that of the random matrices Brod .

Later the development of the statistical spectroscopy was stimulated by quantum chaology which arose in early eighthies as a quantum theory of classical chaotic systems. It was established that quantized spectra of classical chaotic system spectra have the same statistical properties as those of the random matrices. It is well known, that the classical chaotic dynamics is characterised by an exponential divergence of the neighboring trajectories in phase-space. How far the complexity of the classical motion is reflected in the wave dynamics of the corresponding quantum system (described by a linear wave equation) is an interesting aspect of the study of deterministic non-linear dynamics. Numerical evidence has shown that the statistical properties of the energy spectra of the quantum systems depend crucially on the degree of chaos in the underlying classical dynamics. Namely, the spectrum of the quantum operator, describing classically non-integrable systems, has the same statistical properties as those of random matrices. It has been established that certain distribution functions and fluctuation measures for a wide variety of complex nuclear spectra follow a universal behavior and that moreover they agree very well with the results obtained by ensemble averaging over certain matrix ensembles with randomly distributed matrix elements. In particular, for systems with time-reversal symmetry, good agreement was found with a Gaussian orthogonal ensemble(GOE) Eckh88 ; Gutz90 . It has also been established, that for integrable systems the nearest neighbor spacing distribution is dominated by level clustering and is expected to be Poisson-like. Later, this has been confirmed in a number of studies of model systems like billiard balls, simple rotators etc.Eckh88 ; Bohig ; Selig  as well as in realistic systems such as hydrogen atom in a magnetic field Wint ; Del . Details of the progress of level statistics and random matrix theory are given in a recent review Guhr . Here we deal with the spectral statistics of the quark-anitquark systems(quarkonia). Quarkonium spectrosopy has been the subject of extensive experimental as well as theoretical studies AL -Mukh93 . Theoretical description of quarkonium properties can be done quantum mechanically as well as within the QCD approach. Quantum mechanical description of hadrons is reduced to solving wave equations with Coulomb plus a confining potential. Considerable progress has been made in the calculation of the spectra of quarkonium Quig -Zhu  and baryons Mukh93 . At present a large amount of theoretical and experimental data is available. Recently energy fluctuations and quantum chaos in hadrons and QCD has become a subject of theoretical studies Hal -Gu . In partuclar, it has been found that QCD is governed by quantum chaos in both confinement and deconfinement phases Hal ; Bit . The statistical analysis of the measured meson and baryon spectra shows that there is quantum chaos phenomenon in these systems Pasc . The study of the charmonium spectral statistics and its dependence on color screening has established quantum chaotic behaviour of this system Gu . It was claimed that such a behaviour could be the reason for $`J/\mathrm{\Psi }`$ supression Gu .

In the present work we concentrate our attention on level fluctuations and quantum chaos in heavy quarkonia due to the confining forces. It is well known that many realistic and model confined dynamical systems such as an atom in a magnetic field, particle motion in resonators, quantum dots, billiards exhibit chaotic dynamics. Quarks in the hadrons are also confined systems, and to some extent are equivalent to motion in billiards(domain with hard walls where elastic scattering of the particle will occur ) Bun ; Ber ; Zas . Depending on the shape of the billiard domain the motion can be regular or chaotic Ber . If one uses for quark- antiquark potential Coulomb plus harmonic oscillator potential, the system becomes equivalent to the atom in a magnetic filed, which is also described by a Coulomb plus harmonic oscillator potential whose dynamics is chaotic in classical and quantum cases Del ; Wint . An analysis of the quarkonium would suggest that the dynamics of such a system could be chaotic in classical as well as in quantum cases . Therefore the present paper may be viewed as an attempt to use well-known technique to the quarkonium systems to establish their chaotic behaviour.

Charmonium and bottomonium spectral statistics based on the solution of the Schrödinger equation for this system is considered . Two types of potentials for quark-antiquark interaction are considered: Coulomb plus linear and Coulomb plus harmonic oscillator potentials. Solving the Schrödinger equation with these potentials by the scaling variational method we obtain the energy levels of the quarkonium system. This paper is organized as follows: In the next section we solve the Schrödinger equation analytically using the scaling variational method and find the mass sectra of charmonium and bottomonium. In section 3 we calculate the statistical properties of the spectra of charmonium and bottomonium and in section 4 some conclusions and final remarks are presented.

## II Schrödinger equation for quarkonium

To obatain the quantum spectra the following Schrödinger equation is solved

$$(\frac{1}{2\overline{M}}\mathrm{}\frac{Z}{r}+\lambda r^k)\mathrm{\Psi }=E\mathrm{\Psi },$$
(1)
where $`Z=\frac{4}{3}\alpha _s`$, with $`\alpha _s`$ being the strong coupling constant, $`\lambda `$ is the parameter of the confining potential, $`k=1,2`$ and $`\overline{M}=m_q/2`$ is the reduced mass of quarkonium (for simplicity we consider quarks with the same flavor). To solve this equation we use the scaling variational method. This method Castro  was shown to work well for Coulomb plus confining type potentials. Recently this method was used for the calculation of doubly heavy baryon spectra Mat00  and heavy-light meson spectra Matr02 . For analysis, eq. (1) may be written as

$$(H_0+H^{})\mathrm{\Psi }_{nl}=E_{nl}\mathrm{\Psi }_{nl},$$
(2)
where

$$H_0=\frac{1}{2}\mathrm{}\frac{Z}{r},$$
(3)
and

$$H^{}=\lambda r^k.$$
(4)
According to the scaling variational method we choose the trial wave function of the form

$$\mathrm{\Psi }_{nl}(b,r)=b^{3/2}\varphi _{nl}(br),$$
(5)
where $`\varphi _{nl}(br)`$ are the the eigenfunctions of the Hamiltonian $`H_0`$, i.e.

$$H_0\varphi _{nl}=\epsilon _{nl}\varphi _{nl},$$

satisfying the the normalization condition

$$<\varphi _{nl}|\varphi _{n^{}l^{}}>=\delta _{nn^{}}\delta _{ll^{}}.$$

and the eigenvalues are given as

$$\epsilon _{nl}=\frac{Z^2}{2n^2}.$$

Then the eigenvalues of eq.(1) can be calculated as

$$E_{nl}=ϵ_{nl}(b_0),$$
(6)
where

$$ϵ_{nl}(b)=b^3<\varphi (b,r)|H_0+H^{}|\varphi (b,r)>,$$

and $`b_0`$ is found by minimizing this energy, i.e.

$$\frac{ϵ}{b}=0.$$
(7)
The quantity $`ϵ_{nl}(b)`$ may be calculated as:

$$ϵ_{nl}(b)=n^2(\frac{b^2}{2}Zb\frac{qb^k}{k}),$$

where

$$q=k\lambda n^2<\varphi _{nl}|r^k|\varphi _{nl}>.$$

These matrix elements are calculated analytically Lan . For $`k=1`$ we have

$$q=\lambda n^2\frac{1}{2}(3n^2l(l+1)).$$

where $`n`$ and $`l`$ are the principal and orbital quantum numbers, respectively. Then the extremum condition gives

$$b_0=\frac{Z}{3}+\frac{2^{\frac{1}{3}}Z^2}{3(27q+2Z^3+3^{\frac{3}{2}}(q(27q+4Z^3))^{\frac{1}{2}})^{\frac{1}{3}}}+$$

$$\frac{3(27q+2Z^3+3^{\frac{3}{2}}(q(27q+4Z^3))^{\frac{1}{2}})^{\frac{1}{3}}}{2^{\frac{1}{3}}3}.$$

For $`k=2`$ we have

$$q=5n^2+13l(l+1).$$

The extremum condition (7) leads to an equation for defining $`b_0`$

$$b_0^4Zb_0^3q=0.$$

Solving this equation we get $`b_0`$ that leads to the energy eigenvalues for the case $`k=2`$. Thus the solutions of the Schrödinger equation for Coulomb plus confining potentials are obtained analytically using the scaling variational method.

## III Spectral statistics

Now an analysis of the statistical properties of the spectra is to be considered. One of the main characteristics of the statistical properies of the spectra is the level spacing distribution Eckh88 ; Gutz90 ; Bohig  function. We calculate the nearest-neighbor level-spacing distribution for quarkonium spectra. The nearest neighbor level spacings are defined as $`S_i=\stackrel{~}{E}_{i+1}\stackrel{~}{E}_i`$, where $`\stackrel{~}{E}_i`$ are the energies of the unfolded levels, which are obtained by the following way: The spectrum $`\{E_i\}`$ is separated into smoothed average part and fluctuating parts. Then the number of the levels below $`E`$ is counted and the following staircase function is defined:

$$N(E)=N_{av}(E)+N_{fluct}(E).$$

The unfolded spectrum is finally obtained with the mapping

$$\stackrel{~}{E}_i=N_{av}(E_i).$$

Then the nearest level spacing distribution function $`P(S)`$ is defined as the probability of $`S`$ lying within the infinitesimal interval $`[S,S+dS]`$.

For the quantum systems which are chaotic in the classical limit this distribution function is the same as that of the random matrices Eckh88 ; Guhr . For systems which are regular in the classical limit its behaviour is close to a Poissonian distribution function. This distribution is usually taken to be a Gaussian with a parameter $`d`$:

$$P(H)exp(Tr\{HH^+\}/2d^2),$$

and the random matrix ensemble corresponding to this disribution is called the Gaussian ensemble. For Hamiltonians invariant under rotational and time-reversal transformations the corresponding ensemble of matrices is called Gaussian orthogonal ensemble (GOE). It was established Bohig ; Selig ; Ber85  that GOE describes the statistical fluctuation properties of a quantum system whose classical analog is completely chaotic. The nearest neighbor level spacing distribution for GOE is described by Wigner distribution:

$$P(S)=\frac{1}{2}\pi Sexp(\frac{1}{4}\pi S^2).$$
(8)
The common way to study of the level statistics is to compare the calculated nearest-neighbor level-spacing distribution histogram with the Wigner distribution.

For systems whose classical motion is neither regular nor fully chaotic (mixed dynamics) the level spacing distribution will be intermediate between the Poisson and GOE limits. Several empirical functional forms for the distribution were suggested for this case Guhr . If the Hamiltonian is not time-reversal invariant, irrespective of its behavior under rotations, the Hamiltonian matrices are complex Hermitian and the corresponding ensemble is called Gaussian unitary ensemble (GUE). If the Hamiltonian of the systems is time-reversal invariant but not invariant under rotations, then the corresponding ensemble is called the Gaussian symplectic ensemble. Quarkonium system is time reversal and rotational invariant. Therefore in the case of chaotic motion its nearest-neighbor level spacing distribution is expected to be GOE-type. Now we carry out a detailed analysis of the spectra to establish the nearest-neighbor level spacing distribution function to determine the existence of quantum chaos in our system.

The following values of the $`c`$ and $`b`$ quark masses and potential parameters are choosen in these claculations:
For charmonium:
$`m_c=1.486GeV`$, $`\alpha _s=0.32`$, $`\lambda =0.2GeV^2,0.4GeV^2,0.6GeV^2`$.
For bottomonium :
$`m_b=4.88GeV`$, $`\alpha _s=0.22`$, $`\lambda =0.4GeV^2,0.6GeV^2,0.8GeV^2`$. First 1000 eigenavleus are used to plott each histogram. In our claculations we vary only the principal quantum number $`n`$, keeping the orbital quantum number $`l`$ constant.

Figs. 1 and 2 show the level spacing distribution for charmonium and bottomium, respectively for Coulomb plus linear potential ($`k=1`$) case, and in Fig.s 3 and 4 the level spacing distributions are plotted for Coulomb plus harmonic oscillator potential ($`k=2`$) case, at various values of the confining potential parameter. It is clear that for the smaller values of the parameter $`\lambda `$ this distribution is close to GOE, that means the quantum system is becoming chaotic. It allows us to estimate the values of $`\lambda `$ where a transition from regular to chaotic behaviour will occur. Furthermore we find that for the same values of $`\lambda `$ the bottomonium motion is more regular than that of the charmonium. Comparison of the level spacing distribution for $`c\overline{c}`$ and $`b\overline{b}`$ systems shows us that the transition from regular to chaotic motion for these systems occurs at different values of $`\lambda `$. So, for $`c\overline{c}`$ system (Figs.1 and 2) we have GOE-type distributions at $`\lambda =0.2GeV^2`$ and $`\lambda =0.4GeV^2`$ and Poissonian distribution at $`\lambda =0.6GeV^2`$, while the motion of $`b\overline{b}`$ system is still chaotic at $`\lambda =0.6GeV^2`$ (Figs.3 and 4). At $`\lambda =0.4GeV^2`$ and $`\lambda =0.8GeV^2`$ the level spacing distribution for bottomonium is GOE and Poisson type, respectively. Thus it may be concluded from this comparison that the motion of $`b\overline{b}`$ system is more chaotic than that of the $`c\overline{c}`$ system.

Indeed, if one writes the Schrödinger equations for charmonium and bottomonium in the units $`m_q=\mathrm{}=1`$, the value of the confining potential for bottomonium is much smaller than that of charmonium. The chaotic nature of the motion in the quarkonium potential can be also understood from the following fact: formally the Hamiltonian for Coulomb plus harmonic oscillator potential is equivalent to the Hamiltonian of an atom in an uniform magnetic field, which is a nonintegrable system. It is well known, that the Hamiltonian of an atom in an uniform magentic field can be reduced to the Hamiltonian of coupled nonlinaer oscillators Del , which exhibits chaotic dynamics. With this analogy the change in the confining potential parameters to some extent is equivalent to the change of the shape of billiard domain.

## IV Conclusions

Summarizing, in this work we have applied the technique well-known in nuclear physics, atomic physics and other areas of physics to treat quantum chaos in heavy quarkonia, using for their description the Schrödinger equation with Coulomb plus confining potentials, a non-integrable system. It should be emphasized that this system is quite similar to the case of an atom in a constant magnetic field. It is well known that such a system is non-integralble and exhibits chaotic dynamics in both classical as well as quntum casesDel ; Wint . The analysis of the nearest level spacing distributions of these systems shows, that for some smaller values of the confining potential parameter the motion of the heavy quarkonium is chaotic while for larger values of $`\lambda `$ the level spacing distribution is a Poisson type, that implies that the motion is becoming regular. Previously Gu  the quantum chaos is related to the color screening parameter, while the present work treats the transition from regular to a chaotic motion related to a confining potential parameters. The quarkonia can exhibit regular as well as chaotic dynamics depending on the value of the confining force. Study of the level statistics and quantum chaos in hadrons could play the same role in particle physics as it did in the case of nuclear physics. Therefore quantum chaology of hadrons could become one of the exciting topics in particle physics. More comprehensive analysis of the statistical properties of the heavy quarkonia spectra performed by solving the Dirac equation with Coulomb plus a confining potential is needed. This would establish the role played by relativistic effects in the approach to a level-spacing distribution.

## Acknowledgements

The work of DUM is supported by NATO Science Fellowship of Natural Science and Engineering Research Council of Canada (NSERCC). The work of FCK is supported by NSERCC. The work of URS is spupported by a Grant of the Uzbek Academy of Sciences (contract No 33-02).


format: text

text:
# Adiabatic Compression of Soliton Matter Waves

## 1 Introduction

The discovery of Bose-Einstein condenstates (BEC) in atomic vapors of alkali metals has opened new possibilities to observe matter waves of soliton type. These excitations were first discovered in BEC with repulsive interatomic interactions, in the form of dark solitons , and more recently, in BEC with attractive interatomic interactions, in the form of bright solitons moving on zero backgrounds . Since soliton matter waves are of primary importance for developing concrete applications, it is of interest to devise methods which allow to control them. One possibility is to vary the atomic scattering length by means of external magnetic fields, i.e. by using Feshbach resonances. In Ref., it was shown that an abrupt change in time of the atomic scattering length, can lead to the splitting of the soliton with the generation of new solitons. The fragmentation of the wave obviously decreases the number of atoms contained in the original pulse, this being undesirable for applications such as atom lasers.

In the present paper we suggest to use adiabatic variations of the atomic scattering length, both in time and in space, as an effective way for controlling soliton’s parameters and to induce changes in their shape which could be useful for applications. In contrast to abrupt variations, adiabatic changes make possible to preserve the integrity of the soliton (no splitting occurs), this leading for bright solitons to the compression of the pulse with the increase of the matter density, and for dark solitons to the compression of the hole with a decrease of the background level (matter from background move inside the hole). These phenomena are shown to exist both in presence and in absence of a parabolic confining potential. In the last case we derive analytical equations for soliton parameters in terms of adiabatic invariants and soliton perturbation theory. The results of this analysis is found in good agreement with direct integration of the 1D Gross-Pitaevskii equation (GPE). The effect of a confining parabolic potential on the adiabatic dynamics of the soliton is also investigated by numerical simulations. We find that, for bright solitons, except for the oscillatory motion around the bottom of the trap, the phenomena of pulse compression is practically the same as in absence of the trap (this is particularly true for solitons initially at rest in the bottom of the trap). For dark soliton the phenomenon is only qualitatively preserved, due to the boundaries effects introduced by the trap. The possibility to compress BEC solitons could be an experimental tool to investigate the range of validity of the 1D GPE. Since the quasi one dimensional regime is valid for low densities, it would be indeed interesting to see how far one can compress a soliton in a real experiment by means of adiabatic changes of the scattering length. Effects of adiabatic perturbations on soliton dynamics were also investigated in the context of Josephson junctions  and in nonlinear optics . In contrast with Josephson junctions and optical fibers, which require structural changes or preparation of new samples, the study of adiabatic nonlinear perturbations on BEC solitons appears more natural and easy to perform, since the strength of the nonlinear interaction can be changed by using only external fields. We also remark that soliton dynamics in a quasi one-dimensional BEC under time-dependent linear potential was recently studied in .

The paper is organized as follows. In Section 1 we perform an analytical study of the effects on bright and dark BEC solitons induced by spatial and temporal adiabatic variations of the scattering length. For temporal variations we use both a variational approach and perturbation theory to obtain equations for soliton’s parameters both in presence and absence of a confining parabolic potential. In the former case the problem can be linked to an adiabatic invariant of the Kepler problem with results depending only on adiabaticity (smoothness) but not on the smallness of the perturbation. For spatial variations we use standard soliton’s perturbation theory. In Section 2 we compare the results of this analysis with direct numerical integrations of the full 1D GPE both in absence and in presence of an harmonic trap. In the last section we summarize the main results of the paper.

## 2 Analysis

The dynamics of a dilute trapped BEC is described by the Gross-Pitaevskii equation (GPE)

$$i\mathrm{}\varphi _t(r,t)=\frac{\mathrm{}^2}{2m}\mathrm{\Delta }\varphi (r,t)+V_{tr}(r,t)\varphi (r,t)+\mathrm{\Gamma }|\varphi |^2\varphi (r,t),$$
(1)
with $`\mathrm{\Gamma }=\frac{4\pi \mathrm{}^2a_s}{m}`$, $`m`$ being the atomic mass, $`a_s`$ the scattering length, and $`V_{tr}=\frac{1}{2}m(\omega _1^2x^2+\omega _2^2y^2+\omega _3z^2)`$ represents the harmonic trap. The quasi 1D geometry corresponds to the case $`\omega _1^2<<\omega _{2,3}^2`$. To model adiabatic variation of the atomic scattering in 1D (cigar shaped) BEC we consider the following normalized Gross-Pitaevskii equation

$$i\psi _t+\psi _{xx}+\sigma \gamma (x,t)\psi |\psi |^2\omega ^2x^2\psi =0,$$
(2)
where $`\psi `$ is the ground state wavefunction of the condensate, $`\gamma (x,t)`$ is a slowly varying function of space and time, $`\sigma \pm 1`$ corresponds to the case of negative and positive scattering length $`a_s`$ ($`\omega `$ denotes the longitudinal frequency of the trap). In Eq. (2) the space has been normalized with respect to the healing length $`\xi =1/\sqrt{8\pi \rho |a_s|}`$ ($`\rho `$ is the atomic density), while the time with respect to $`t_0=m\xi ^2/\mathrm{}`$. We remark that the 1D approximation is valid for the number of atoms $`N<l_r/|a_s|`$, where $`l_r=\sqrt{\mathrm{}/(m\omega _{2,3})}`$ is the transverse oscillator length. It is also worth to mention that in the strong interaction limit $`a_sN|\psi |^2>>1`$, ($`a_s>0`$), the effective nonlinearity is of the type $`|\psi |\psi `$  so that deviations from Eq. (2) occur (this will not be considered in this paper). Although the analysis can be performed for a generic smooth functions $`\gamma (x,t)`$, we shall restrict to the limiting cases: $`\gamma \gamma (t)`$, and $`\gamma \gamma (x)`$, the former being experimentally more easy to realize.

### 2.1 Case $`\gamma \gamma (t)`$

To study adiabatic temporal variations of the scattering length, case $`\gamma \gamma (t)`$, we remark that, due to the number of atoms conservation, it is natural to look for solution of Eq. (2) with $`\sigma =+1`$ (negative scattering length) in the form of a bright soliton with time dependent parameters

$$\psi (x,t)=A(t)\text{sech}(\frac{x\zeta (t)}{a(t)})e^{ib(t)(x\zeta (t))^2+iC(x\zeta (t))+i\varphi (t)},$$
(3)
where $`A(t)`$, $`a(t)`$, $`\varphi (t)`$, $`\zeta (t)`$, $`C(t)`$, $`b(t)`$, denote, respectively, amplitude, width, phase, center, velocity and ”chirp” oscillation of the perturbed soliton. Taking these parameters (with their derivatives), as collective coordinates for the soliton, one can derive their time evolution from the Euler-Lagrange equations with respect to the space averaged Lagrangian $`\overline{L}=L(x,t)𝑑x`$ of Eq. (2) . This gives for the center of mass: $`\zeta _{tt}+4\omega ^2\zeta =0`$, and for the width :

$$a_{tt}+4\omega ^2a=\frac{16}{\pi ^2a^3}\frac{4\gamma (t)N}{\pi ^2a^2}.$$
(4)
Notice that Eq. (4) has the form of a perturbed Kepler problem . This analogy can be used to reformulate the problem in terms of the evolution of a unit mass moving in an adiabatically varying Kepler potential. Since action-angle variables for the unperturbed case are known, one can describe the adiabatic change of the soliton width by means of the adiabatic invariant of the Kepler problem. As the result we find for $`\omega ^2<1`$ that

$$a\frac{4}{\gamma (t)N}[1\frac{64\omega ^2}{(\gamma (t)N)^4}].$$
(5)
We see that for $`\omega =0`$ one recovers the law $`a=4/(\gamma (t)N)`$ known from nonlinear fiber optics. In this case, one can also show that the frequency of the width oscillation is given by $`\omega _a\gamma ^2(t)N^2/4\pi `$. The analysis can be easily extended to two dimensions (2D) for the radially symmetric GPE. In this case the variational equation for the width becomes $`a_{tt}+2\omega ^2a=\frac{2(2N_1\gamma _1(t))}{a^3}`$, where $`N_1=N/(2\pi ),\gamma _1=\gamma /(\sqrt{2}\pi l_x)`$. Applying the above arguments, one finds that the width of a 2D soliton changes as $`a(2N_1\gamma _1(t))^{1/4}/\omega ^{1/2}`$. Eq. (5) predicts for small $`\omega `$ a dependence $`1/\gamma `$ of the soliton width . This result can also be obtained by perturbation theory. To this end we assume in the following that the range of variation of the potential is very large in comparison with the size of the condensate, so that we can neglect the $`\omega `$ term in Eq. (2) (this is particularly appropriate for static bright solitons located at the bottom of the trap, as we will show later). In this case the transformation $`\psi =u/\sqrt{\gamma }`$ (for $`\gamma >0`$), allows to rewrite Eq. (2) in the form of a perturbed nonlinear Schrödinger equation

$$iu_t+u_{xx}+|u|^2u=i(\mathrm{ln}\sqrt{\gamma })_tu+(\mathrm{ln}(\gamma ))_xu_x+0(\frac{T_s^2}{T^2},\frac{L_s^2}{L^2})=R(u),$$
(6)
where $`T`$ and $`L`$ are, respectively, the temporal and spatial characteristic scales of the inhomogeneity, while $`T_s,L_s`$ are characteristic time and space soliton scales, assumed to be much smaller than the corresponding inhomogeneity scales. In the case $`\gamma _x=0`$, the right hand side of Eq. (6) represents a small linear time dependent amplification (damping). If we consider the initial condition as a single soliton

$$u(x,t)=\sqrt{2}A(t)\text{sech}(A(t)(x\zeta ))e^{i\mathrm{\Theta }(t)},\mathrm{\Theta }=k(t)+C(t)(x\zeta (t)),$$
(7)
then, using the equation for the energy $`N=|u|^2𝑑x`$, we find $`A=A_0\gamma (t)`$ this giving, for the amplitude and width of the soliton, $`A_\psi =A_0\sqrt{\gamma },a=1/(A_0\gamma )`$, in agreement with the variational approach.

From this we conclude that an adiabatic increase of the scattering length can be used to narrowing the width of a bright soliton matter wave. Soliton compression phenomena induced by linear damping amplification are also known from nonlinear optics . A similar analysis can be performed for the adiabatic evolution of a dark soliton ($`\sigma =1`$ in Eq.(2)). In this case, the adiabatic evolution of the dark soliton $`\psi _d`$ and of the background $`u_B`$ can be found as

$$\psi _d=u_0\mathrm{tanh}(\frac{u_0\sqrt{\gamma }}{\sqrt{2}}(x))e^{i\theta },u_B=u_0\sqrt{\gamma }\mathrm{exp}(i\theta ),\theta =iu_0^2_0^t\gamma (t^{})𝑑t^{},$$
(8)
where $`u_0=\sqrt{\mu }`$, with $`\mu `$ the chemical potential of the condensate (notice that this last expression for $`u_B`$ is exact). From this expression we see that the width of the dark soliton is changing according to $`a_d=a_{d0}/\sqrt{\gamma (t)}`$, i.e. a factor $`1/\sqrt{\gamma (t)}`$ less in comparison with bright solitons.

### 2.2 Case $`\gamma \gamma (x)`$

The perturbative approach can be used to study also the case of spatial variations of the scattering length i.e. $`\gamma \gamma (x)`$. We have from Eq. (6): $`R(u)=F(x)u_x`$, $`F=(\mathrm{ln}\gamma (x))_x`$, with $`R(u)`$ considered to be a small perturbation. Using the perturbation theory for solitons , we find that the equations for the soliton’s parameters (7) are:

$`A_t`$ $`=`$ $`AC{\displaystyle _{\mathrm{}}^{\mathrm{}}}F({\displaystyle \frac{y}{A}}+\zeta )\text{sech}^2(y)𝑑y2ACF(\zeta ),`$
$`C_t`$ $`=`$ $`A^2{\displaystyle _{\mathrm{}}^{\mathrm{}}}F({\displaystyle \frac{y}{A}}+\zeta )(\text{sech}^2(y)\text{sech}^4(y))𝑑y{\displaystyle \frac{2}{3}}A^2F(\zeta ),`$ (9)
$`\zeta _t`$ $`=`$ $`2C+{\displaystyle \frac{C}{A}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}F({\displaystyle \frac{y}{A}}+\zeta )y\text{sech}^2(y)𝑑y2C(1+O({\displaystyle \frac{1}{L^2}})).`$
From these equations it follows that ($`F(\mathrm{})=1`$)

$$A=A_0\gamma (\zeta ),C_{fin}=\sqrt{C_{in}^2+\frac{1}{3}A_0^2(\gamma ^2(\zeta )\gamma _{ini}^2)}.$$
(10)
In the next section we shall compare these predictions with direct numerical integrations of Eq.(2).

## 3 Numerical results

To check these results we have numerically integrated Eq. (2) with the function $`\gamma (t)`$ given by

$$\gamma (t)=2+A_s\left[\frac{1}{2}+\frac{1}{\pi }\mathrm{tan}^1(s\pi )(t\frac{T_f}{2})\right]$$
(11)
This function models, for small $`s`$, an adiabatic change of the scattering length, while for large $`s`$ it reduces to the step function of amplitude $`A_s`$ centered at $`T_f/2`$ (to have the function starting from $`2`$, the quantity $`\delta =\gamma (0)2`$ must be subtracted). In the following we shall use this function to model also spatial variations of the scattering length. In Fig. 1 we depict the width of a bright (panel a) and a dark soliton (panel b) as a function of time for different values of $`A_s`$, in the case $`\omega =0`$. From this figure we see that the oscillations of the soliton width after the transition, are more pronounced for dark solitons than for bright ones. This is a consequence of the fact that in the dark case the background is also excited by the perturbation and contribute for a relevant part to the dynamics. In Fig. 2a we show the average final value of the width of bright and dark solitons as a function of $`A_s`$. It is evident that, for the same absolute variation of the scattering length, bright solitons are more compressed than dark ones, the difference being just a factor proportional to $`1/\sqrt{\gamma }`$, as expected from our analysis. The agreement with the theoretical predictions is indeed rather good especially for higher values of $`A_s`$. Notice that in this figure the influence of a parabolic trap on the compression phenomenon is also shown for the case of a bright soliton initially at rest in the bottom of the trap (open circles). We see that the results are almost identical to those in absence of trap (this is quite expected, since the pulse is localized in the bottom of a low intensity trap). The panel (b) of the same figure shows the amplitude of a bright soliton as a function of the position of the center of mass for the case of a space dependent variation of the scattering length in absence of the parabolic trap. The soliton, is initially at rest, is sucked into the higher scattering length region, and reaches a constant velocity after passing the inhomogeneity. Also in this case we find an excellent agreement with our analysis (see the inset of the figure).

The effect of the parabolic trap on pulse compression is further investigated in Fig. 3. In panels (a) and (b) of this figure, we show the evolution of a bright soliton oscillating in the trap and subject, respectively, to adiabatic and abrupt variation of the scattering length. We see that in contrast with the adiabatic case, which always leads to compression, abrupt variations of the scattering length may induce the splitting of the soliton into sub-solitons (notice in panel (b), that three solitons, indicated by the arrows, are formed at time $`t=54`$). The splitting phenomenon is similar to the one observed in absence of parabolic potential . The number of sub-solitons created can also be determined in analogy with this last case (we omit details for brevity). Panels (c) and (d) of Fig. 3, show the time evolution of a dark soliton in absence and in presence of a parabolic trap. We remark that, while in the bright case the compression of the soliton is obviously accompanied by an increase of its amplitude, in the dark case the amplitude cannot decrease below zero, so that the compression must be combined with a decrease of the background level (the matter from the background move inside the hole so to reduce its size). This is particularly evident in panel c for the case $`\omega =0`$. We remark that this case, although physically non realisable for BEC, it is of interest for nonlinear optics, and shows the role played by the confinement potential in the case of BEC with positive scattering lengths (it should be compared with panel(d) discussed below).

In analogy with bright solitons, one could expect that an abrupt change in time of the scattering length may split a dark soliton ground state into dark and grey sub-components if the amplitude of the perturbation is big enough. This is indeed what we numerically observed. For the case $`\omega =0`$, using the Inverse Scattering Transform, one can predict that out of an initial profile of the form $`\psi (x,t=0)=b\mathrm{tanh}(cx)`$, one black soliton and a number $`N_0`$ of pairs of grey solitons (with $`N_0=n1`$, $`n`$ being the integer part of the ratio $`b/c`$) , should be generated. By denoting with $`\gamma _1`$ and $`\gamma _2`$ the initial and final values of $`\gamma `$, one can write the initial condition for the equation: $`iv_t+v_{xx}|v|^2v=0`$, as $`v(0)=u_0\sqrt{\gamma _2}\mathrm{tanh}(u_0\sqrt{2\gamma _1}x)`$, with $`v=\psi \sqrt{\gamma _2}`$, and the ratio $`b/c`$ equal to $`\sqrt{\gamma _2/2\gamma _1}`$. We checked that this analysis indeed gives correct predictions of splitting in absence of the parabolic trap. The influence of the parabolic trap on a dark soliton initially at the rest in its bottom, is reported in panel (d) of Fig. 3. Although the phenomenon resemble the one observed in absence of the trap (notice the decrease of the background after the transition), the solution is more complicated due to the boundary effects introduced by the trap. In particular we see that the compression of the hole in the center is less effective than for $`\omega =0`$, and the decrease of the background is accompanied by an expansion of the condensate. This is a consequence of the non constancy of the background (it obviously must go to zero at large distances) which permits the matter to escape up in the trap (the pressure on the hole is thus partially lost in the expansion). The same phenomenon exists also for dark solitons initially oscillating in the trap. In analogy with the case $`\omega =0`$, we also found that abrupt variations in presence of parabolic potential induce dark soliton splitting.

In Fig. 4 we show the averaged amplitude of the soliton (i.e. the time average of the modulo squared at the center of the profile) versus the parameter $`s`$ controlling the rapidity of change of the function $`\gamma (t)`$. From this figure we see that a crossover from adiabatic to fast variations behavior occurs at $`s1.2`$. Notice that this crossover value of $`s`$ corresponds to a variation function $`\gamma (t)`$ which is already quite rapid (see top inset of the figure). The adiabatic region is characterized by the fact that the final profiles of the compressed solitons are almost independent on $`s`$ (they depend only on the amplitude of the scattering length variation). This can be seen from the bottom inset of Fig. 4 in which the initial ($`t=0`$) and final ($`t=100`$) soliton profiles corresponding to points a and b in the adiabatic region, are reported. Very little oscillations of the soliton profile are found in these cases. On the contrary, the non-adiabatic (or fast oscillation) region is characterized by the presence of large oscillations of the soliton profile, induced by the rapid variation of the scattering length. This is illustrated in point c of the lower inset of Fig. 4, in which we depicted, besides the initial profile at $`t=0`$, two snapshots of the soliton (max and min of the profile) to evidenciate its oscillations at the final time $`T_f=100`$. These oscillations increase by increasing $`s`$ and eventually lead to the splitting of the soliton as illustrated in panel b of Fig. 3 (the splitting, however, depends also on the amplitude of the scattering length variation, and usually occurs much above the interface between adiabatic and fast oscillation behavior). From this analysis we conclude that the adiabatic compression of matter wave solitons is possible in a wide range of the parameter $`s`$ (i.e. $`0<s1`$) and is compatible with relatively fast variations of the scattering length, as one can see from the top inset of Fig. 4.

## 4 Conclusions

In conclusion, we have shown that adiabatic changes of the scattering length, both in space and time, can be effectively used to control parameters of BEC matter waves of soliton type. In particular we showed the effect of pulse compression both on dark and bright solitons in presence and absence of a confining parabolic potential. The influence of the smoothness of the transition was also considered. We found that deviations from adiabatic predictions occur only for very large values of $`s`$, i.e. for almost step-like variations. This implies that the results of this paper are valid for generic smooth variations of the scattering length and therefore are expected to be found in real BEC matter waves experiments. The possibility of sharpening bright and dark solitons by means of smooth variations of the scattering length, beside providing another confirmation of the soliton nature of the BEC ground state, may be very useful for concrete applications. We hope that experiments in this direction will be soon performed.

We are grateful to V.V. Konotop for interesting discussions. F.Kh.A. acknowledges the hospitality received at the Physics Department of the University of Salerno, where this work was done, and partial financial support from NATO-Linkage grant No. PST.CLG.978177. M. S. thanks the European grant LOCNET no. HPRN-CT-1999-00163 for partial support.

## References


text:
# Effect of Impurities in the Channel of a Spin Field Effect Transistor (SPINFET)

## 1 Introduction

In a seminal paper published in 1990, Datta and Das  proposed a gate controlled electron spin interferometer which is an analog of the standard electro-optic light modulator. Their device (later dubbed ”Spin Field Effect Transistor” or SPINFET) consists of a one-dimensional semiconductor channel with ferromagnetic source and drain contacts (Fig. 1(a)). Electrons are injected into the channel from the ferromagnetic source with a definite spin orientation, which is then controllably precessed in the channel with a gate-controlled Rashba interaction , and finally sensed at the drain. At the drain end, the electron’s transmission probability depends on the relative alignment of its spin with the drain’s (fixed) magnetization. By controlling the angle of spin precession in the channel with a gate voltage, one can modulate the relative spin alignment at the drain end, and hence control the source-to-drain current (or conductance). This is the principle of the SPINFET.

There have been some studies of spin transport in such a device. Mireles and Kirczenow  carried out a study of ballistic spin transport, but they overlooked two crucial features that are always present in a real device structure. First, there is an axial magnetic field along the channel caused by the ferromagnetic source and drain contacts which are magnetized in the same direction. This field can be quite strong ($``$ 1 Tesla)  and dramatically alters the dispersion relations of the subbands in the channel, causes spin mixing, and has a serious effect on spin transport. Second, there will always be a few impurities in the channel (even if they are remote impurities) associated with channel doping or unintentional defects. We show that even a single (non-magnetic) impurity can cause spin relaxation in the presence of the axial magnetic field.

## 2 Theoretical Approach

We first consider the quasi one-dimensional semiconductor channel of a SPINFET in the absence of any impurities. The channel is along the x-axis (Fig. 1(a)) and the gate electric field is applied along the y-direction to induce a Rashba spin-orbit coupling in the channel. This system is described by the single particle effective-mass Hamiltonian

$`={\displaystyle \frac{1}{2m^{}}}\left(\stackrel{}{p}+e\stackrel{}{A}\right)^2+V_I(x)+V_1(y)+V_2(z)`$
$`(g^{}/2)\mu _B\stackrel{}{B}\stackrel{}{\sigma }+{\displaystyle \frac{\alpha _R}{\mathrm{}}}\widehat{y}\left[\stackrel{}{\sigma }\times (\stackrel{}{p}+e\stackrel{}{A})\right]`$ (1)
where $`\widehat{y}`$ is the unit vector along the y-direction in Fig. 1(a) and $`\stackrel{}{A}`$ is the vector potential due to the axial magnetic field $`\stackrel{}{B}`$ along the channel (x-direction) caused by the ferromagnetic contacts. In (1), $`\mu _B`$ is the Bohr magneton ( $`e\mathrm{}/2m_0`$) and $`g^{}`$ is the effective Landé g-factor of the electron in the channel. The quantity $`\alpha _R`$ is the Rashba spin-orbit coupling strength which can be varied with the gate potential. The confining potentials along the y- and z-directions are denoted by $`V_1(y)`$ and $`V_2(z)`$, with the latter being parabolic in space.

In (1), $`V_I(x)`$ represents an interfacial potential barrier between the ferromagnetic contacts and the semiconducting channel. If the contact neighborhood consists of heavily doped semiconductor material in close proximity to a metallic ferromagnet, the Schottky barriers at the interface will be very narrow and electrons from the contacts can tunnel fairly easily into semiconducting channel resulting in a nearly-ohmic contact. We model these ultra-narrow Schottky barriers as delta-barriers given by:

$$V_I(x)=V_L\delta (x)+V_R\delta (xL)$$
(2)
where $`V_L`$ and $`V_R`$ are assumed equal ($`V_L`$ = $`V_R`$ = $`V_0`$). In practice, the strength of the barrier depends on the ferromagnetic materials and also on the doping level in the channel. These barriers have a beneficial effect; they can facilitate coherent spin injection across a metallic ferromagnet and a semiconducting paramagnet interface  which is crucial for the SPINFET.

In (1), we have neglected a few effects for the sake of simplicity. We have neglected the normal Elliott-Yafet interaction  because it is weak in quasi one-dimensional structures (where elastic scattering is strongly suppressed). We have also neglected the Dresselhaus interaction  since it does not relax spin when the initial spin polarization is along the axis of the wire . The Dresselhaus interaction can however be easily included in the Hamiltonian and is left for future work.

The choice of the Landau gauge $`\stackrel{}{A}`$ = (0, -Bz, 0) allows us to decouple the y-component of the Hamiltonian in (1) from the x-z component. Furthermore, in the absence of any impurity scattering potential ($`V_{imp}=0`$), the Hamiltonian in the semiconducting channel is translationally invariant in the x-direction and the wavevector $`k_x`$ is a good quantum number. The eigenstates of the system can then be determined using plane waves traveling in the x-direction . The two-dimensional Hamiltonian in the plane of the channel (x-z plane) is therefore given by

$`H_{xz}`$ $`=`$ $`{\displaystyle \frac{p_z^2}{2m^{}}}+\mathrm{\Delta }E_c+V_I(x)+{\displaystyle \frac{1}{2}}m^{}\left(\omega _0^2+\omega _c^2\right)z^2+{\displaystyle \frac{\mathrm{}^2k_x^2}{2m^{}}}`$ (3)
$`+{\displaystyle \frac{\mathrm{}^2k_Rk_x}{m^{}}}\sigma _z(g^{}/2)\mu _BB\sigma _x{\displaystyle \frac{\mathrm{}k_Rp_z}{m^{}}}\sigma _x`$
where $`\omega _0`$ is the curvature of the confining potential in the z-direction, $`k_R=m^{}\alpha _R/\mathrm{}^2`$, and $`\mathrm{\Delta }E_c`$ is the potential barrier between the ferromagnet and semiconductor. We assume that $`\mathrm{\Delta }E_c`$ includes the effects of the quantum confinement in the y-direction (Fig.1(a)).

To model localized non-magnetic impurities (i.e., which do not by themselves flip the spin) we use a standard model of delta-scatterers

$$V_{imp}=\mathrm{\Gamma }_{imp}\delta (xx_i)$$
(4)
for an impurity located at a distance $`x_i`$ from the left ferromagnet/semiconductor interface with scattering strength $`\mathrm{\Gamma }_{imp}`$ (assumed to be spin independent). In our numerical examples, we consider the case of both attractive ($`\mathrm{\Gamma }_{imp}`$ negative) and repulsive ($`\mathrm{\Gamma }_{imp}`$ positive) impurities. While (1) and (3) represent a ballistic channel with no scattering, addition of the scattering potential in (4) to (1) or (3) will result in a Hamiltonian describing a weakly disordered channel in which impurity scattering takes place. The eigenstates of this (spin-dependent) Hamiltonian can then be found using a transfer matrix technique to extract the electron wavefunction. From this wavefunction, we can calculate the (spin-dependent) transmission probability through the channel and ultimately the (spin-dependent) channel conductance. The details of the calculations have been presented elsewhere . The linear response conductance of the spin interferometer (for injection from either the +x or -x polarized bands in the left contact) is found from the Landauer formula.

We model the ferromagnetic contacts by the Stoner-Wohlfarth model (Fig.1.(b)). The +x-polarized spin (majority carrier) and -x-polarized spin (minority carrier) band bottoms are offset by an exchange splitting energy $`\mathrm{\Delta }`$. The solutions of the Schrödinger equation for injection from minority and majority spins from the left magnetic contact then be written throughout the entire device leading to the corresponding transmission probabilities and the conductance of the interferometer based on the Landauer formula (see ). The strength of the barrier at the ferromagnet/semiconductor interface is characterized by the following parameter

$$Z=\frac{2m_{f}^{}{}_{}{}^{}V_0}{\mathrm{}^2}.$$
(5)
Typical values of $`Z`$ vary in the range of 0 to 2 . Using $`m_{f}^{}{}_{}{}^{}=m_0`$ and $`k_F`$ = 1.05x$`10^8`$ cm<sup>-1</sup>, we get a barrier strength $`V_0`$ = 16 eV$`\AA `$ for Z = 2.

## 3 Numerical Examples

We consider a spin interferometer consisting of a quasi one-dimensional InAs channel between two ferromagnetic contacts. The electrostatic potential in the z-direction is assumed to be harmonic (with $`\mathrm{}\omega _0`$ = 10 meV in (3). A Zeeman splitting energy of 0.34 meV is used in the semiconductor channel assuming a magnetic field B = 1 Tesla along the channel. This corresponds to a $`g^{}`$ factor of 3 and an electron effective mass $`m^{}=0.036m_o`$ which is typical of InAs-based channels . The Fermi level $`E_f`$ and the exchange splitting energy $`\mathrm{\Delta }`$ in the ferromagnetic contacts are set equal to 4.2 and 3.46 eV, respectively .

The Rashba spin-orbit coupling strength $`\alpha _R`$ is typically derived from low-temperature magnetoresistance measurements (Shubnikov-de Haas oscillations) in 2DEG created at the interface of semiconductor heterostructures . To date, the largest reported experimental values of the Rashba spin-orbit coupling strength $`\alpha _R`$ has been found in InAs-based semiconductor heterojunctions. For a normal HEMT $`In_{0.75}Al_{0.25}As/In_{0.75}Ga_{0.25}As`$ heterojunction, Sato et al. have reported variation of $`\alpha _R`$ from 30- to 15 $`\times 10^{12}`$ eV-m when the external gate voltage is swept from 0 to -6 V (the total electron concentration in the 2DEG is found to be reduced from 5- to 4.5$`\times 10^{11}/cm^2`$ over the same range of bias).

Tuning the gate voltage varies both the potential energy barrier $`\mathrm{\Delta }E_c`$ and the Rashba spin-orbit coupling strength $`\alpha _R`$. Both of these variations lead to distinct types of conductance oscillations. The variation of $`\mathrm{\Delta }E_c`$ causes the Fermi-level in the channel to sweep through the resonant energies in the channel, causing the conductance to oscillate. These are known as Ramsauer oscillations and have been examined by us in detail in . The variation of $`\alpha _R`$, on the other hand, causes spin precession in the channel leading to the type of conductance oscillation which is the basis of the spin interferometer, as originally visualized by Datta and Das . In  we found that the Ramsauer oscillations are much stronger and can mask the oscillations due to spin precession, unless the structure is designed with particular care to eliminate (or reduce) the Ramsauer oscillations.

In the numerical examples below, we assumed that the Rashba spin-orbit coupling strength $`\alpha _R`$ is constant (i.e., independent of the gate potential) and equal to the maximum reported value of 30 $`\times 10^{12}`$ eV-m (best case scenario). We also used a value of Z = 0.25 corresponding to a value of $`V_L`$ and $`V_R`$ in (2) equal to 2 eV-$`\AA `$. The scattering strength of the single impurity in the channel was set equal to $`\mathrm{\Gamma }`$ = $`\pm `$ 0.25 eV-$`\AA `$ (the plus and minus sign corresponding to repulsive and attractive impurity, respectively).

First, we consider the case of a single repulsive scatterer (impurity) whose location is varied from the left to the rigth side of the channel in steps of 10 $`\AA `$. Figure 2 shows that variation of the conductance with the impurity location for three different values of $`\mathrm{\Delta }E_c`$ in the channel, i.e, gate potential on the SPINFET. It is seen that the condutance of the interferometer is a strong function of the impurity location. Figure 2 shows that the condutance modulation (difference between maximum and minimum conductance values) of the interferometer is about 0.9, 0.21, and 0.2 $`e^2/h`$, for $`\mathrm{\Delta }E_c`$ equal to 4.191, 4.188, and 4.185 eV, respectively. The same calculation were repeated for the case of an attractive impurity with the same magnitude of the scattering strength and the conductance modulation as a function of the impurity location are shown in Fig.3. In this case, the condutance modulation (difference between maximum and minimum conductance values) of the interferometer is about 0.76, 0.16, and 0.17$`e^2/h`$, for $`\mathrm{\Delta }E_c`$ equal to 4.191, 4.188, and 4.185 eV, respectively.

## 4 Conclusions:

In this paper, we have shown how the conductance of gate controlled spin interferometers proposed in  is strongly dependent on the location of a single impurity in the semiconducting channel (Figs. 2 and 3). The extreme sensitivity of the amplitude and phase of conductance oscillations of a SPINFET to impurity location is reminiscent of the phenomenon of universal conductance fluctuations of mesoscopic samples . This will hinder practical applications of electron spin interferometers since it will lead to such problems as large threshold variability, random device characteristics, and general irreproducibility.

The work of S. B. was supported by the National Science Foundation.


text:
# Comment on “London Theory for Superconducting Phase Transitions in External Magnetic Fields: Application to "UPt"₃”

The theory of the equilibrium flux lattice in $`UPt_3`$ developed in Ref.  is valid for the superconducting A-phase near the phase transition to the normal state, where the symmetry of the superconducting state is mostly dictated by one of the two components of the order parameter $`𝜼=(\eta _1,\eta _2)`$. An attempt to spread the analysis  towards the phase transition between the superconducting A and B states has been undertaken in the recent Letter . In our opinion, the approach and the conclusions of Ref.  are wrong and here we present why it is so.

The Ginzburg-Landau equations for $`Hz`$ after taking just for brevity $`\beta _2=0`$ are

$`\left(\alpha _0\tau _1+2\beta _1\left(\left|\eta _1\right|^2+\left|\eta _2\right|^2\right)+\left(K_{123}D_x^2+K_1D_y^2\right)\right)\eta _1=`$ (1)
$`\left(K_2D_xD_y+K_3D_yD_x\right)\eta _2,`$ (2)
$`\left(\alpha _0\tau _2+2\beta _1\left(\left|\eta _1\right|^2+\left|\eta _2\right|^2\right)+\left(K_{123}D_y^2+K_1D_x^2\right)\right)\eta _2=`$ (3)
$`\left(K_3D_xD_y+K_2D_yD_x\right)\eta _1,`$ (4)
where $`\tau _1=(TT_1)`$, $`\tau _2=(TT_2)`$ with the temperatures $`T_1>T_2`$. Near the phase transition to the normal state, that is when $`HH_{c2}=\mathrm{\Phi }_0\alpha _0\tau _1/2\pi \sqrt{K_1K_{123}}`$ the component $`\eta _2`$ is expressed through the first component $`\eta _1`$ :

$$\eta _2=\frac{\left(K_3D_xD_y+K_2D_yD_x\right)\eta _1}{\alpha _0\left(TT_{c2}\right)},$$
(5)
and $`\eta _1`$ is solution of the nonlinear Eq. (1) at $`\eta _2=0`$ (here $`T_{c2}=T_22\beta _1\left|\eta _1\right|^2/\alpha _0`$).

It is clear that in the region of the neighbourhood of the AB transition ($`TT_{c2}`$), the simple linear but nonlocal relationship (3) between the two order parameter components is divergent and out of applicability. The grow of the amplitude of $`\eta _2`$ near the transition between the A and B states is limited in fact principally by the nonlinear (and local) terms in the free energy. The mixing gradient terms (right hand side of Eq (1)-(2)) transform the AB and AC phase transition lines into regions of crossover between the superconducting states with the same symmetry. The crossovers are certainly not characterized by a divergence of the corresponding coherence length as it has to be in the case of a real second order phase transition.

The authors of the Letter  disregard the contribution of the nonlinear term $`2\beta _1\left|\eta _2\right|^2\eta _2`$ in Eq. (2) by considering for $`HH_{c1}`$ (with $`K_2=K_3K_1`$) the relationship

$$\left(\alpha _0\left(TT_{c2}\right)+K_1𝐃^\mathrm{𝟐}\right)\eta _2=K_2\left(D_xD_y+D_yD_x\right)\eta _1.$$
(6)
The divergence at $`TT_{c2}`$ is avoided here by the presence of the differential operator acting on $`\eta _2`$. Substituting Eq. (4) back to the gradient part of the free energy, Agterberg and Dodgson have found an effective nonlocal (4-th order in gradients) London energy functional for fields and currents originating from $`\eta _1`$. These additional nonlocal terms in the London energy starts to be important when $`TT_{c2}`$. As the results, the authors of Ref.  have derived many conclusions concerning the structures of the Abrikosov lattice in vicinity of the AB transition.

It is possible to estimate under which conditions the nonlinear term is negligible at $`TT_{c2}`$. Using Eq. (4), we have

$$\frac{\beta _1\left|\eta _2\right|^4}{K_1\left|𝐃\eta _2\right|^2}\left(\frac{K_2}{K_1}\right)^2\kappa ^2\left(\frac{\lambda }{d}\right)^4,$$
(7)
where $`\kappa 60`$ is the Ginzburg-Landau parameter for $`UPt_3`$, $`d`$ the distance between vortices and $`\lambda `$ the London penetration depth. As $`HH_{c1}`$, $`d\lambda (\mathrm{ln}\kappa )^{1/2}`$ so that the nonlinear term is found negligible when

$$\frac{K_2}{K_1}<\left(\kappa \mathrm{ln}\kappa \right)^1.$$
(8)
On the other hand, the approach of Ref.  does not consider the usual nonlocal terms originating from the 4-th order gradient terms acting on $`\eta _1`$ in the free energy, which is valid if

$$K_2/K_1>\kappa ^1.$$
(9)
Taking account of (6) and (7), we thus conclude that the nonlinear terms are important in the vicinity of the AB crossover.


text:
# Oblate-prolate transition in odd-mass light mercury isotopes

## 1 Introduction

The sharp increase in the mean-square charge radius is observed in the odd mass 181-185Hg isotopes with respect to their even neighbors . A hint that a very different configuration may constitute the ground state of these nuclei was given by an isomer shift measurement of 185Hg , where they found an excited state of 185Hg which has a charge radius consistent with the charge radius of heavier neighbors. The sharp increase in the charge radius was interpreted as the transition from a nearly spherical to a deformed shape in the ground state of the odd mass 181-185Hg isotopes.

Frauendorf and Pashkevich  calculated the deformation energy of Hg isotopes by means of Strutinsky’s shell correction method. They found that the light Hg isotope has an oblate and prolate minimum of nearly the same depth at the respective deformations $`ϵ=0.12`$ and 0.22. They discussed a mechanism which causes the transition from a small oblate to a large prolate deformation in the odd mass 181-185Hg isotopes. They noticed that the neutron level density is very high in the almost spherical oblate minimum, whereas the neutron level density is low in the well deformed prolate minimum. They concluded that the loss of pairing energy due to the blocking of the Fermi level by the odd neutron is larger in the oblate minimum than in the prolate one, which gives rise to the oblate-prolate transition in the ground state of these isotopes. As the prolate minimum corresponds to deformation twice as large as the oblate one, the sharp increase in the mean-square radius is observed in the odd-mass light Hg isotopes.

There are several mean field calculations of Hg isotopes. However, all of them are limited to even isotopes . In order to see if the mechanism of Frauendorf and Pashkevich works well, we perform Hartree-Fock-Bogoliubov (HFB) calculations of Hg isotopes with effective interactions SIII , SkI3  and SLy4 . We calculate Lipkin-Nogami corrections  to both mean field and pairing energies . We investigate the influence of the time-odd mean-field terms as well as the blocking of a level by the odd neutron on the oblate-prolate shape transition.

In Section 2, we present the Skyrme HFB method with Lipkin-Nogami corrections to describe both even and odd Hg isotopes. In Section 3, we present calculations with effective interactions SIII, SkI3 and SLy4. Conclusions are given in Section 4.

## 2 Method of calculation

### 2.1 Skyrme Hartree-Fock-Bogoliubov method

The HFB equation is derived by the variation of energy with respect to the quasi-particle vacuum ,

$$H=\underset{\alpha \beta }{}t_{\alpha \beta }c_{\alpha }^{}{}_{}{}^{}c_\beta +\frac{1}{4}\underset{\alpha \beta \gamma \delta }{}\overline{v}_{\alpha \beta \gamma \delta }c_{\alpha }^{}{}_{}{}^{}c_{\beta }^{}{}_{}{}^{}c_\delta c_\gamma ,$$
(1)
$$\delta \frac{\mathrm{\Phi }|H|\mathrm{\Phi }}{\mathrm{\Phi }|\mathrm{\Phi }}=0.$$
(2)
Here the quasi-particle vacuum $`|\mathrm{\Phi }`$ is given by

$$\beta _k=\underset{\alpha }{}U_{\alpha k}^{}c_\alpha +V_{\alpha k}^{}c_{\alpha }^{}{}_{}{}^{},$$
(3)
$$\beta _k|\mathrm{\Phi }=0$$
(4)
for all quasi-particle states $`k`$. In the present study, we employ Skyrme forces in the particle-hole channel and the zero-range density-dependent pairing force in the particle-particle channel,

$$V_p=\frac{V_0}{2}\left(1P_\sigma \right)\left(1\frac{\rho (𝒓_\mathrm{𝟏})}{\rho _\mathrm{c}}\right)\delta (𝒓_\mathrm{𝟏}𝒓_\mathrm{𝟐}).$$
(5)
The HFB energy is then given by

$$E^{\mathrm{HFB}}=𝑑𝒓+\frac{1}{2}\mathrm{\Delta }(𝒓)\kappa ^{}(𝒓)𝑑𝒓,$$
(6)
where $``$ is the Skyrme energy density  and $`\mathrm{\Delta }(𝒓)`$ is the gap field,

$$\mathrm{\Delta }(𝒓)=\frac{V_0}{2}\left(1\frac{\rho (𝒓)}{\rho _\mathrm{c}}\right)\kappa (𝒓).$$
(7)
The pairing strength $`V_0`$ is determined from the odd-even mass difference and the density parameter $`\rho _\mathrm{c}`$ is taken as the saturation density of symmetric nuclear matter which is obtained by using the respective Skyrme parameters. The particle density $`\rho (𝒓)`$ and the anomalous density $`\kappa (𝒓)`$ are written as

$$\rho (𝒓\sigma ,𝒓\sigma ^{})=\underset{\alpha \beta }{}\rho _{\alpha \beta }\phi _\alpha (𝒓,\sigma )\phi _\beta ^{}(𝒓,\sigma ^{}),$$
(8)
$$\kappa (𝒓\sigma ,𝒓\sigma ^{})=\underset{\alpha \beta }{}\kappa _{\alpha \beta }\phi _\alpha (𝒓,\sigma )\phi _\beta (𝒓,\sigma ^{})$$
(9)
with the density matrix and pairing tensor,

$$\rho _{\alpha \beta }=\mathrm{\Phi }|c_{\beta }^{}{}_{}{}^{}c_\alpha |\mathrm{\Phi }=\underset{k}{}V_{\alpha k}^{}V_{\beta k},$$
(10)
$$\kappa _{\alpha \beta }=\mathrm{\Phi }|c_\beta c_\alpha |\mathrm{\Phi }=\underset{k}{}V_{\alpha k}^{}U_{\beta k}.$$
(11)
Variation of Eq. (6) with respect to $`\rho _{\alpha \beta }`$ and $`\kappa _{\alpha \beta }`$ yields the HFB equation,

$$\left(\begin{array}{cc}h_{\alpha \beta }\lambda \delta _{\alpha \beta }& \mathrm{\Delta }_{\alpha \beta }\\ \mathrm{\Delta }_{\alpha \beta }^{}& h_{\alpha \beta }^{}+\lambda \delta _{\alpha \beta }\end{array}\right)\left(\begin{array}{c}U_{\beta k}\\ V_{\beta k}\end{array}\right)=E_k\left(\begin{array}{c}U_{\alpha k}\\ V_{\alpha k}\end{array}\right).$$
(12)
Diagonalization is carried out in a deformed oscillator basis with axial symmetry . A Greek index in Eq. (12) denotes a set of quantum numbers {$`n_rn_z\mathrm{\Lambda }\mathrm{\Sigma }`$} of a deformed oscillator wave function.

In the present study we will restrict our basis space as

$$E_\alpha =(2n_r+|\mathrm{\Lambda }|+1)\mathrm{}\omega _{}+\left(n_z+\frac{1}{2}\right)\mathrm{}\omega _z\left(N_0+\frac{3}{2}\right)\mathrm{}\omega _0$$
(13)
with $`N_0=18`$. The oscillator constant is held fixed as $`\mathrm{}\omega _0=\mathrm{}(\omega _{}^2\omega _z)^{1/3}=1.21\times \frac{\mathrm{}^2}{m}A^{1/3}`$. The deformation parameter $`q=\omega _{}/\omega _z`$ is varied to obtain the maximum deformation energy in the oblate and prolate minimum, respectively. The chemical potential $`\lambda `$ is adjusted by the particle number condition, $`\widehat{N}=N`$. Using the occupation amplitude $`N_k=_\alpha V_{\alpha k}^{}V_{\alpha k}`$ and quasi-particle energy $`E_k`$, we calculate an auxiliary single-particle energy and pairing gap in order to find the optimum value of $`\lambda `$ ,

$$\overline{e}_k=E_k(2N_k1)+\lambda ,$$
(14)
$$\overline{\mathrm{\Delta }}_k=2E_k\sqrt{N_k(1N_k)}.$$
(15)
As in Ref. , the cut-off procedure of quasi-particle states is imposed on the auxiliary single-particle energy. The value of $`\overline{e}_{\mathrm{max}}=60`$ MeV is assumed in the present study.

### 2.2 Description of an odd nucleus

The above formalism can be straightforwardly extended to the case of one quasi-particle state in odd mass nuclei. All the equation retain their form. Only the density matrix $`\rho _{\alpha \beta }`$ and pairing tensor $`\kappa _{\alpha \beta }`$ should be modified due to the blocked quasi-particle $`\beta _k^{}|\mathrm{\Phi }`$ in the following way ;

$`\stackrel{~}{\rho }_{\alpha \beta }`$ $`=`$ $`\rho _{\alpha \beta }+U_{\alpha k}U_{\beta k}^{}V_{\alpha k}^{}V_{\beta k},`$ (16)
$`\stackrel{~}{\kappa }_{\alpha \beta }`$ $`=`$ $`\kappa _{\alpha \beta }+U_{\alpha k}V_{\beta k}^{}V_{\alpha k}^{}U_{\beta k}.`$ (17)
We introduce a reference state $`|\mathrm{\Phi }_{\mathrm{HFBE}}`$ and its energy expectation value $`E^{\mathrm{HFBE}}`$. The state is constructed as an even vacuum without quasi-particle creation and without breaking time-reversal invariance but with an odd average particle number. The difference of energy $`E^{\mathrm{HFB}}E^{\mathrm{HFBE}}`$ is thus composed of the loss of pairing energy due to the blocking of the Fermi level and the polarization energy of the mean field by the presence of an odd neutron .

A proper description of an odd nucleus by mean field theories requires to break the time-reversal symmetry . When the time-reversal symmetry is broken, the Skyrme energy density $``$ is written by the sum of kinetic energy and potential energy of isoscalar $`_0`$ and isovector $`_1`$ terms;

$``$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}\tau _0+_0+_1,`$
$`_t`$ $`=`$ $`C_t^\rho \rho _t^2+C_t^s𝒔_t^2+C_t^\mathrm{}\rho \rho _t\mathrm{}\rho _t+C_t^\mathrm{}s𝒔_t\mathrm{}𝒔_t`$ (18)
$`+C_t^\tau (\rho _t\tau _t𝒋_t^2)+C_t^J(𝑱_t^2𝒔_t𝑻_t)+C_t^J[\rho _t𝑱_t+𝒔_t(\times 𝒋_t)].`$
The coupling constants $`C`$ are expressed by the Skyrme force parameters . The isospin index $`t`$ can have values 0 or 1. The isoscalar and isovector parts of the particle densities $`\rho _t`$ are defined as

$$\rho _0=\rho _n+\rho _p,\rho _1=\rho _n\rho _p,$$
(19)
and analogous expressions are used to define other densities. It should be noted that since we consider the Skyrme force as the density-dependent two-body interaction, the time-odd terms $`C_t^s`$ in Eq. (18) are different from the original expressions  which were derived from the three-body interaction.

The characteristic combinations of $`(\rho _t\tau _t𝒋_t^2)`$, $`(𝑱_t^2𝒔_t𝑻_t)`$ and $`[\rho _t𝑱_t+𝒔_t(\times 𝒋_t)]`$ in Eq. (18) are the results of the local gauge invariance of the Skyrme force . The time-odd coupling constants are uniquely determined from the original Skyrme-force parameters and therefore the time-odd mean field has no adjustable parameters once the time-even mean field is given. However, as discussed by Dobaczewski and Dudek , one may consider the energy density to be a more fundamental construction than the Skyrme interaction itself. In such a case, all 20 coupling constants in Eq. (18) can be treated and adjusted independently. However, the gauge-invariant conditions restrict the values of six time-odd coupling constants and leave the freedom to modify the values of $`C_t^s`$ and $`C_t^{\mathrm{\Delta }s}`$. We are going to vary $`C_t^s`$ and $`C_t^{\mathrm{\Delta }s}`$ to see the influence of these terms on the oblate-prolate transition in the odd-mass 181-185Hg isotopes. We will omit the terms $`(𝑱_t^2𝒔_t𝑻_t)`$ because $`𝑱_t^2`$ are omitted in the parameterization of the present Skyrme forces.

### 2.3 Lipkin-Nogami corrections

In light Hg isotopes, the neutron level density is high in the oblate minimum, while it is low in the prolate minimum. Because of the difference of level densities in the respective minima, it happens that neutrons in the oblate minimum are superconductive, whereas they are normal in the prolate minimum within the same isotope. The situation is extremely inconvenient because we want discuss the energy difference between the oblate and prolate minimum. Besides, a sudden phase-transition from the normal to superconductive state or vice versa at certain neutron numbers may be unphysical and it may be due to the defect of the number non-conserving pairing scheme. In the present study we employ the approximate number-projection method of Lipkin and Nogami (LN)  and calculate its corrections to both mean field and pairing energies . The number projected energy is expressed as

$$E_{\mathrm{HFBLN}}=E_{\mathrm{HFB}}+\lambda _1(N\widehat{N})+\lambda _2(N^2\widehat{N}^2)+\mathrm{},$$
(20)
where $`\lambda _1`$ and $`\lambda _2`$ are constant parameters . The second order expansion is considered in the Lipkin-Nogami method. Variation with respect to the density matrix and pairing tensor as in the case of the HFB equation yields the HFBLN equation applicable to both even and odd nuclei (hereafter simply called HFB). The parameter $`\lambda _1`$ is given by the particle number condition, whereas the parameter $`\lambda _2`$ is held fixed during the variation and is determined after the variation from the additional condition,

$$|\xi =\mathrm{exp}(i\xi \widehat{N}_{20})|\mathrm{\Phi },$$
(21)
$$\lambda _2=\frac{_\xi ^2\mathrm{\Phi }|H|\xi |_{\xi =0}}{_\xi ^2\mathrm{\Phi }|\widehat{N}^2|\xi |_{\xi =0}}.$$
(22)
Equation (22) is calculated by the following density matrix and pairing tensor ;

$`\rho _{\alpha \beta }^{(\xi )}`$ $`=`$ $`\mathrm{\Phi }|c_\beta ^{}c_\alpha |\xi =\left(\rho +2i\xi (1\rho )\rho \right)_{\alpha \beta },`$ (23)
$`\kappa _{\alpha \beta }^{(\xi )}`$ $`=`$ $`\mathrm{\Phi }|c_\beta c_\alpha |\xi =\left(\kappa +2i\xi (1\rho )\kappa \right)_{\alpha \beta },`$ (24)
$`\overline{\kappa }_{\alpha \beta }^{(\xi )}`$ $`=`$ $`\mathrm{\Phi }|c_\alpha ^{}c_\beta ^{}|\xi =(\kappa ^{}2i\xi \kappa ^{}\rho )_{\alpha \beta }.`$ (25)
For an odd nucleus, the density matrix and pairing tensor of Eqs. (16) and (17) are used on the right hand side of the above equations .

## 3 Results and discussion

We have analysed the oblate-prolate transition in the odd-mass light Hg isotopes using the Skyrme Hartree-Fock-Bogoliubov method with effective interactions SIII, SkI3 and SLy4. Shape transition occurs only in the odd-mass isotopes as the result of a subtle balance between the time-even mean-field energy, the loss of pairing energy by the odd neutron and the time-odd mean-field energy. We have examined effective interactions SIII, SGII , SkM , SKX , SkI1, SkI2, SkI3, SkI4 and SkI5  and SLy4 whether they can predict the oblate-prolate shape coexistence in light Hg isotopes. We have found that both SIII and SkI3 interactions have desirable isospin dependence and deformation energies.

We have assumed the zero-range density-dependent pairing interaction. The pairing strength $`V_0`$ was adjusted to reproduce the one-particle separation energy . The density parameter $`\rho _c`$ is taken as the saturation density of symmetric nuclear matter which is obtained from Skyrme parameters. Adjusted values are given in the respective figure captions.

We have investigated the effect of time-odd mean-field terms on the HFB energy, first by using the time-odd coupling constants derived from the Skyrme force parameters and second by switching on and off the coupling constants $`C_t^{\mathrm{\Delta }s}`$ and $`C_t^s`$, respectively. These coupling constants in $`C`$-representation  are given in Table 1. They vary very much with different Skyrme forces, e.g., the coupling constant $`C_0^s(\rho =0)`$ is 14, 84 and –208 for the respective SIII, SkI3 and SLy4 forces.

The polarization energy depends strongly on the blocked level by an odd neutron. Even the sign of energy changes with different choice of the blocked level. In the present study we simply choose as the blocked level the lowest quasi-particle state that is obtained from the HFBE calculation.

Mercury isotopes are very interesting because the shape of the ground state varies very much as the neutron number decreases from the closed shell. Figure 1 shows the deformation energy of even Hg isotopes calculated by the constraint SIII-HF method. We see the spherical shape in the neutron-closed 206Hg. The oblate minimum with deformation $`\beta _20.17`$ then develops in 200Hg. A rather complicated prolate minimum appears in 190Hg. Finally the deep prolate minimum with deformation $`\beta _20.29`$ develops in 184Hg. The oblate-prolate shape coexistence in light Hg isotopes has been predicted by many authors, e.g., .

Figure 2 shows the SIII calculation of the oblate-prolate energy difference. The HF calculation (white circle) predicts prolate deformation in the ground state of all Hg isotopes. The HFBE calculation (black circle), the fully paired state without the blocking effect, predicts oblate deformation in the ground state of all Hg isotopes. The shape changes due to the large gain of pairing energy in the oblate minimum because of the high level density as compared to the level density in the prolate minimum. In the SIII-HFB calculation (white square), on the other hand, the shape of odd mass 181-185Hg isotopes returns to the prolate shape. This is because the loss of pairing energy due to the blocking by the odd neutron is also large in the oblate minimum than in the prolate one. The effect of the time-odd mean fields on the oblate-prolate transition is weak and a few tenths of the blocking effect.

Figure 3 shows the SIII-HFB calculation of (a) the mean-square charge radius and (b) one-neutron separation energy of 180-206Hg isotopes. The sharp increase in the mean-square radius of the odd-mass 181-185Hg isotopes is very well explained by the calculation. The one-neutron separation energy is also well reproduced in the calculation. This is so because the depth of the pairing force has been adjusted to reproduce the observed separation energies. It is noted that the white circle in Fig. 3(b) is a derived value from systematic trends  and not an experiment.

As is well known, a proper description of odd nuclei by mean field theory requires to break the time-reversal symmetry and introduces the time-odd components (TOC) in the mean field. The effect of an odd neutron on the HFB energy is therefore twofold; (a) The loss of pairing energy due to the blocking of the Fermi level by the odd neutron, and (b) the time-odd mean-field energy arising from the blocked level. The time-odd mean fields have been studied both in the analysis of super-deformed rotational bands  and the odd-even nuclear mass difference .

In order to see the blocking energy (the loss of pairing energy) and the polarization energy (the time-odd mean-field energy) separately, we define the blocking energy $`E^{\mathrm{block}}`$ as the energy difference between $`E^{\mathrm{HFB}}`$ with no TOC and $`E^{\mathrm{HFBE}}`$. We also define the polarization energy $`E^{\mathrm{pol}}`$ as the energy difference between $`E^{\mathrm{HFB}}`$ and $`E^{\mathrm{HFB}}`$ with no TOC, i.e.,

$`E^{\mathrm{block}}`$ $`=`$ $`E^{\mathrm{HFB}}(\text{no TOC})E^{\mathrm{HFBE}},`$ (26)
$`E^{\mathrm{pol}}`$ $`=`$ $`E^{\mathrm{HFB}}E^{\mathrm{HFB}}(\text{no TOC}).`$ (27)
Figure 4 shows the SIII-HFB calculation of (a) the blocking energy in the oblate and prolate minimum, (b) the polarization energy in the oblate minimum, (c) the polarization energy in the prolate minimum and (d) the oblate-prolate difference of HFB energy. In Fig. 4(a), the blocking energy in the oblate minimum is larger by $``$200 keV than the one in the prolate minimum. This results from higher neutron level-density in the former than in the latter. The blocking energies are reversed in 189Hg because of the small deformation $`\beta =0.08`$ in the prolate minimum of this isotope and the level density in the minimum is high accordingly.

As mentioned in Section 2.2, the $`C_t^s`$ and $`C_t^\mathrm{}s`$ terms are free from the time-even parameterization. To investigate their effects on the HFB energy, we have further divided the polarization energy into contributions with (i) the full time-odd terms, (ii) the time-odd terms without $`C_t^s`$, (iii) the time-odd terms without $`C_t^{\mathrm{\Delta }s}`$, and (iv) the time-odd terms without both $`C_t^s`$ and $`C_t^{\mathrm{\Delta }s}`$. In the following we assume simultaneous modifications of the isoscalar and isovector coupling constants of a given species.

Figures 4 (b) and (c) display the contributions to the polarization energy. The polarization energy with full time-odd terms (white circle) is quite small. It is attractive and less than –50 keV in the oblate minimum, while it is repulsive and less than 50 keV in the prolate minimum. When the $`C_t^s`$ terms are omitted (black circle), the polarization energy becomes attractive and more than –250 keV in the oblate minimum and it is attractive around –150 keV in the prolate minimum. The $`C_t^s`$ terms have repulsive contribution to the HFB energy. When the $`C_t^{\mathrm{\Delta }s}`$ terms are omitted (white square), on the other hand, the polarization energy becomes repulsive and $``$50 keV in both minima. The $`C_t^{\mathrm{\Delta }s}`$ terms have large attractive contributions to the HFB energy. When both $`C_t^s`$ and $`C_t^{\mathrm{\Delta }s}`$ terms are omitted (black square), the polarization energy becomes small. The effect of the $`C_t^s`$ and $`C_t^{\mathrm{\Delta }s}`$ terms is opposite and cancels with each other.

Figure 4(d) shows the difference of energy in the oblate and prolate minimum. The figure shows the calculation with the full time-odd terms (white circle), without the $`C_t^s`$ terms (black circle), without the $`C_t^{\mathrm{\Delta }s}`$ terms (white square), without both $`C_t^s`$ and $`C_t^{\mathrm{\Delta }s}`$ terms (black square) and without all time-odd terms (white triangle). When both $`C_t^s`$ and $`C_t^{\mathrm{\Delta }s}`$ terms are omitted, the oblate shape is predicted in 181Hg, while the prolate shape is calculated in 185Hg with marginal energy difference between the oblate and prolate minima.

Table 2 summarizes the oblate-prolate difference of the HFBE energy $`\mathrm{\Delta }E^{\mathrm{HFBE}}`$, the blocking energy $`\mathrm{\Delta }E^{\mathrm{block}}`$, the polarization energy $`\mathrm{\Delta }E^{\mathrm{pol}}`$ and the sum of these contributions $`\mathrm{\Delta }E^{\mathrm{HFB}}`$. The polarization energy $`\mathrm{\Delta }E^{\mathrm{pol}}`$ is calculated with four different time-odd mean fields. They are the HFB calculation (i) with full time-odd terms, (ii) without the $`C_t^s`$ terms, (iii) without the $`C_t^{\mathrm{\Delta }s}`$ terms and (iv) without both $`C_t^s`$ and $`C_t^{\mathrm{\Delta }s}`$ terms.

The energy difference $`\mathrm{\Delta }E^{\mathrm{HFBE}}`$ shows that the minimum is found at the neutron number $`N=102`$ and the oblate minimum gets deeper than the prolate one as the neutron number is away from $`N=102`$ in both directions. For the blocking energy $`E^{\mathrm{block}}`$, it is repulsive and larger by 400 keV in the oblate minimum than in the prolate one. The blocked level is labeled with the quantum number $`\mathrm{\Omega }`$, the component of neutron angular momentum along the symmetry axis.

For the odd mass 181-185Hg, the blocked level has $`\mathrm{\Omega }=1/2`$ in the oblate minimum, while it has $`\mathrm{\Omega }=5/2`$ or $`7/2`$ in the prolate minimum. In the calculation with full time-odd terms (case (i)), the polarization energy is attractive in the oblate minimum, while it is repulsive in the prolate minimum, leading to the the attractive energy difference of $`\mathrm{\Delta }E^{\mathrm{pol}}`$. The attractive $`\mathrm{\Delta }E^{\mathrm{pol}}`$ weakens the repulsive $`\mathrm{\Delta }E^{\mathrm{block}}`$ and prevents the oblate-prolate shape transition in these isotopes. In the calculation without the $`C_t^s`$ terms (case (ii)), the polarization energy becomes attractive because of the absence of the repulsive $`C_t^s`$ terms. In this case, the polarization energy $`\mathrm{\Delta }E^{\mathrm{pol}}`$ almost cancels the blocking energy $`\mathrm{\Delta }E^{\mathrm{block}}`$. In the calculation without the $`C_t^{\mathrm{\Delta }\mathrm{s}}`$ terms (case (iii)), the polarization energy is weakly repulsive because of the absence of strongly attractive $`C_t^{\mathrm{\Delta }s}`$ terms. Finally in the calculation without both $`C_t^s`$ and $`C_t^{\mathrm{\Delta }s}`$ terms, the polarization energy is weakly attractive in both minima. The effects of the $`C_t^s`$ and $`C_t^{\mathrm{\Delta }\mathrm{s}}`$ terms on the polarization energy will be discussed more generally in the last paragraph of this section.

Figure 5 shows the results of the SkI3-HFB calculation. We had the difficulty of convergence in this calculation when the full time-odd components are considered. The difficulty arises mainly due to the strong coupling-constant $`C_0^{\mathrm{\Delta }s}`$ derived from the SkI3 parameters (see Table 1) and partly due to the large HF basis of the present study ($`N_0=18`$). The difficulty occurs when the energy splitting between the time-reversed partners becomes large and a few of the quasi-particle states have negative energy. Since we have no way to correct this difficulty, we simply omit the $`C_t^{\mathrm{\Delta }s}`$ terms and display the results in Fig. 5. The blocking energy $`E^{\mathrm{block}}`$ in Fig. 5(a) is similar to the one in the SIII-HFB calculation. Figures 5(b) and 5(c) show the polarization energy (iii) without $`C_t^{\mathrm{\Delta }s}`$ terms (white square) and (iv) without both $`C_t^{\mathrm{\Delta }s}`$ and $`C_t^s`$ terms (black square). The $`C_t^s`$ terms have repulsive contributions of $``$150 keV in the oblate minimum while it is $``$200 keV in the prolate minimum. In Fig. 5(d), the oblate-prolate transition is very well explained when the terms $`C_t^{\mathrm{\Delta }s}`$ are omitted. It is also well explained by omitting the whole time-odd components.

Table 3 summarizes the results of the SkI3-HFB calculation. The energy difference $`\mathrm{\Delta }E^{\mathrm{HFBE}}`$ is nearly the same for even and odd 180-185Hg isotopes. This fact helps very much in reproducing the oblate-prolate transition in the odd mass 181-185Hg isotopes. The SkI3 force is very well suited to describe the mean fields of light Hg isotopes.

Figure 6 shows the results of the SLy4-HFB calculation. The blocking energy $`E^{\mathrm{block}}`$ in the odd-mass 181-185Hg is larger by $``$200 keV in the oblate minimum than in the prolate one. However, their magnitudes are reversed at the neutron number $`N=107`$ and $`N=109`$. The reason is the same as the SIII calculation. The equilibrium deformation in the prolate minimum is reduced to $`\beta =0.12`$ in these isotopes. Figures 6(b) and 6(c) show the terms $`C_t^s`$ have repulsive contributions to the polarization energy. In addition, the terms $`C_t^s`$ and $`C_t^{\mathrm{\Delta }s}`$ have opposite signs and nearly cancel with each other.

The oblate-prolate energy difference in Fig. 6(d) shows that the SLy4 force is not very well suited to describe deformation properties of light Hg isotopes. Here we also had the difficulty of convergence when the $`C_t^s`$ terms are not included in the calculation. The difficulty may be due to the $`C_t^{\mathrm{\Delta }s}`$ terms. The repulsive $`C_t^s`$ terms cancel the attractive $`C_t^{\mathrm{\Delta }s}`$ terms. Without the $`C_t^s`$ terms, the attractive $`C_t^{\mathrm{\Delta }s}`$ terms are quite strong and cause the difficulty of convergence.

Table 4 summarizes the results of the SLy4-HFB calculation. The energy difference $`\mathrm{\Delta }E^{\mathrm{HFBE}}`$ shows that the SLy4 force favors the oblate shape in light Hg isotopes much more than the SIII and SkI3 forces. This fact results in the oblate-prolate transition only in the ground state of 181Hg isotope. The oblate shape is predicted for all heavier isotopes.

To see the effects of the time-odd $`C_t^s`$ and $`C_t^{\mathrm{\Delta }s}`$ terms more closely, we have analysed the odd-mass Sn, Ba, Yb, Hg and Pb isotopes by the SkI3-HFB method with the basis space of $`N_0=10`$. A small basis space was chosen in order to avoid the difficulty of convergence. Results of the SkI3-HF calculation are shown in Fig. 7. The polarization energy is plotted against $`2\times \mathrm{\Omega }^\pi `$ of the blocked level. Polarization energies are calculated (a) with the full time-odd terms, (b) without $`C_t^{\mathrm{\Delta }s}`$ terms and (c) without both $`C_t^{\mathrm{\Delta }s}`$ and $`C_t^s`$ terms. We may conclude that the $`C_t^{\mathrm{\Delta }s}`$ terms are attractive and have a large effect on the blocked level with low $`\mathrm{\Omega }`$, in particular $`\mathrm{\Omega }^\pi =1/2^\pm `$. The terms $`C_t^s`$ are repulsive for levels with any value of $`\mathrm{\Omega }`$. The remaining time-odd terms are attractive for levels with high $`\mathrm{\Omega }`$ and are negligible for levels with low $`\mathrm{\Omega }`$. The same features of the time-odd terms are observed also for the SIII and SLy4 interactions. For the odd-mass 181-185Hg in Tables 2, 3 and 4, we have seen that the time-odd terms have attractive contributions to $`\mathrm{\Delta }E^{\mathrm{HFB}}`$ in the cases (i) and (ii), while the terms have repulsive contributions in the case (iv).

## 4 Conclusions

The shape and the mean-square radius of Hg isotopes have been investigated by using the Skyrme Hartree-Fock-Bogoliubov method. We have first examined effective interactions SIII, SGII, SkM, SKX, SkI1, SkI2, SkI3, SkI4, SkI5 and SLy4 to see if they can predict the oblate-prolate shape coexistence in light Hg isotopes. We have then examined the loss of pairing energy due to the odd neutron by employing the zero-range density-dependent pairing force. Lipkin-Nogami corrections to both even and odd nuclei are calculated within the HFB formalism. We have also investigated the effect of the time-odd $`C_t^{\mathrm{\Delta }s}`$ and $`C_t^s`$ terms on the HFB energy in relation to the blocked level by the odd neutron.

We have obtained the following conclusions. (1) The effective interactions SIII and SkI3 have desirable deformation properties in predicting the oblate-prolate shape coexistence in light Hg isotopes. (2) The blocking energy of the pairing is the most important factor to change the shape of the odd-mass 181-185Hg isotopes. (3) The effect of the time-odd terms on the prolate-oblate transition is weak and cancels the effect of the blocking energy by a few tenths.

We have also investigated the effects of the time-odd terms $`C_t^s`$ and $`C_t^{\mathrm{\Delta }s}`$ by blocking several levels near the Fermi surface of Sn, Ba, Yb, Hg and Pb isotopes. We have found that the terms $`C_t^{\mathrm{\Delta }s}`$ are attractive for the blocked levels with low $`\mathrm{\Omega }`$, in particular $`\mathrm{\Omega }^\pi =1/2^\pm `$. On the other hand, the terms $`C_t^s`$ are repulsive for the blocked levels with all values of $`\mathrm{\Omega }`$. The remaining time-odd terms are attractive for the blocked levels with high $`\mathrm{\Omega }`$ and negligible for the blocked levels with low $`\mathrm{\Omega }`$.

We had the difficulty of convergence in the SkI3-HFB calculation when the $`C_t^{\mathrm{\Delta }s}`$ terms are included. The difficulty also arises from the SLy4-HFB calculation when the $`C_t^s`$ terms are omitted. These difficulties are mainly due to the strong $`C_t^{\mathrm{\Delta }s}`$ terms and partly due to the large basis space of the calculation ($`N_0=18`$). For large $`C_t^{\mathrm{\Delta }s}`$, the energy splitting between the time-reversed partners becomes large and a few of the quasi-particle states have negative energy.

## Acknowledgments

We would like to thank K. Arita for many illuminating discussions. Numerical computation in this work was partly carried out at the Yukawa-Institute Computer Facility, Kyoto University.


text:
# 1 Introduction

## 1 Introduction

The usual calculations of statistical errors for entries in histograms and the use of these errors in the fitting procedure can bias the measurements if there are several entries from the same event. Traditionally, in the past this problem was ignored. The effect is small for low multiplicity events. However, for LEP and especially for the future LHC, RHIC, etc… experiments this effect is not small at all. These entries are correlated creating bin-to-bin and inside bin correlations. Neglecting these correlations, as we will show below, leads to a remarkable underestimation of the errors in the measured quantities and less precise estimation of the quantities themselves.

An effective approach has been proposed in . The exact method decribed here includes all possible correlation effects.

The method was applied to the two-particle correlation functions of $`e^+e^{}`$W+W- events.

## 2 The Method

In this section a new method is presented for construction of the covariance matrix between the bins of the histogram for an unbiased measurement of the fitted quantities.

The presence of bin-to-bin correlations in two particle distributions is unavoidable. If there are $`N`$ positive tracks in the events, each of them has $`(N1)`$ entries in the two-particle density $`P`$, contributing to different bins of the histogram. Also, due to the finite bin width, the same track can also enter several times in the same bin, which is a source of inside bin correlations.

The method is based on classical statistics. Let us consider the $`i`$-th event from the set of $`n`$ events and the two-particle density $`P`$ which is presented in the histogram $`h^i`$ with $`N_p`$ bins.

The histogram $`H=_{i=1}^nh^i`$ and values

$`c_{jk}={\displaystyle _{i=1}^n}(h_j^iH_j/n)(h_k^iH_k/n)(1+1/n)`$
were calculated event by event. Here $`j`$ and $`k`$ are the bin numbers for the histograms. We do not know the correlations and errors for one event. But we know that the different events are uncorrelated. Let us consider bin values of the histogram made for one event as an random vector with unknown distribution. We have an uncorrelated ensemble of these vectors and hence we can estimate the covariance matrix statistically. It is important to note, that this algorithm computes both the “technical” and the physical correlations. For all events we have the resulting histogram H for the two-particle density $`P`$ and $`V_{jk}=c_{jk}n/(n1)`$ covariance matrix for this histogram. Note, that the expresion for unbiased estimation of sample covariance was used for the calculation of covariance between two values. The diagonal terms $`V_{jj}`$ of the covariance matrix are assumed to be the estimate for the squares of the error $`\sigma _j`$ of the $`j`$-th of histogram $`H`$.

## 3 The Test of the Method Using the Simulated Events

The precise measurement of the correlations between particles became important for $`e^+e^{}`$W+W- events due to possible large impact of these correlations on the measured W mass . Thus, we applied the method to the correlation functions in fully hadronic and semileptonic WW events.

The correlation function $`R`$ is used to study the enhanced probability for emission of particles. For pairs of particles, it is defined as

$$R(p_1,p_2)=\frac{P(p_1,p_2)}{P_0(p_1,p_2)},$$
(1)
where $`P(p_1,p_2)`$ is the two-particle probability density, $`p_i`$ is the four-momentum of particle $`i`$, and $`P_0(p_1,p_2)`$ is a reference two-particle distribution which, ideally, resembles $`P(p_1,p_2)`$ in all respects, apart from the lack of Bose-Einstein symmetrization. The effect is usually described in terms of the variable $`Q`$, defined by $`Q^2=(p_1p_2)^2=M^2(\pi \pi )4m_\pi ^2`$, where $`M`$ is the invariant mass of the two pions. The correlation function can then be written as

$$R(Q)=\frac{P(Q)}{P_0(Q)},$$
(2)
which is frequently parametrised by the function

$$R(Q)=1+\lambda e^{r^2Q^2}.$$
(3)
In the above equation, in the hypothesis of a spherically symmetric pion source, the parameter $`r`$ gives the RMS radius of the source and $`\lambda `$ is the strength of the correlation between the pions.

The method described in section 2 was tested on the JETSET simulation . We choose the fully hadronic and semileptonic decay of WW pairs with Bose-Einstein correlations included using the LUBOEI code. The value of $`\lambda =0.85`$ and $`r=0.5`$ fm was used. In case of fully hadronic channel the Bose-Einstein correlations were switched on for all pions(full Bose-Einstein correlations).

The correlation matrices $`\rho _{jk}=V_{jk}/(\sigma _j\sigma _k)`$ for like-sign pairs of WW semileptonic channel (refered as ($`2q`$) mode) and for fully hadronic channel (referred as ($`4q`$) mode), computed using 100 000 simulated events for each sample, are shown in Fig. 1 and Fig. 2. The correlations between bins for the WW fully hadronic channel(Fig. 2) are larger than in the mixed hadronic and leptonic channel(Fig. 1). Thus, the bin-to-bin correlations are increased with multiplicity, as expected. Notice that the correlations are nearly independent of Q, which shows why the effective approach of  is a good approximation in this case.

For the future analysis the 500 samples of 3000 events each for ($`4q`$) channel and 500 samples of 1500 events each for ($`2q`$) channel were simulated. For each of these samples a histogram of the correlation function $`R(Q)`$ was built, using the 25 bins of 100 MeV from 0 to 2.5 GeV. The simulated $`R`$ distributions were normalized to unity in the region $`Q>`$ 0.8 GeV/$`c^2`$. We performed a $`\chi ^2`$ fit to the $`R(Q)`$ to the form (3) for each of 500 samles.

The average values of $`\lambda `$ and $`r`$ from these “naive” fits were:

$`\lambda _{2q}`$ $`=`$ $`0.333\pm 0.029,`$ (4)
$`r_{2q}`$ $`=`$ $`0.562\pm 0.033\text{ fm}`$ (5)
for ($`2q`$) events and

$`\lambda _{4q}`$ $`=`$ $`0.416\pm 0.013,`$ (6)
$`r_{4q}`$ $`=`$ $`0.561\pm 0.012\text{ fm}`$ (7)
for ($`4q`$) events. The statistical errors correspond to the average of the errors for 500 samples. The “pull” of the fitted values of $`\lambda `$ and $`r`$ for the simulated ($`4q`$) events are shown in Fig. 3. A Gaussian fits gave that the errors in parameres $`\lambda `$ and $`r`$ are underestimated by a factor 1.20 $`\pm `$ 0.05 and 1.30 $`\pm `$ 0.05 for the ($`2q`$) events and by a factor 1.42 $`\pm `$ 0.06 and 1.53 $`\pm `$ 0.06 for the ($`4q`$) events.

The average values of $`\lambda `$ and $`r`$ from the 500 fits using the inverted $`V_{jk}`$ matrix (calculated for each of 500 samples) were:

$`\lambda _{2q}`$ $`=`$ $`0.332\pm 0.034,`$ (8)
$`r_{2q}`$ $`=`$ $`0.556\pm 0.040\text{ fm}`$ (9)
for ($`2q`$) events and

$`\lambda _{4q}`$ $`=`$ $`0.403\pm 0.017,`$ (10)
$`r_{4q}`$ $`=`$ $`0.565\pm 0.016\text{ fm}`$ (11)
for ($`4q`$) events. The “pull” of the fitted values of $`\lambda `$ and $`r`$ for the simulated ($`4q`$) events are shown in Fig. 4. A Gaussian fits gave $`\sigma _{(\lambda <\lambda >)}/\sigma _\lambda `$=1.06 $`\pm `$ 0.05 and $`\sigma _{(r<r>)}/\sigma _r`$=1.02 $`\pm `$ 0.05 for (2q) channel, and $`\sigma _{(\lambda <\lambda >)}/\sigma _\lambda `$=0.96 $`\pm `$ 0.05 and $`\sigma _{(r<r>)}/\sigma _r`$=1.08 $`\pm `$ 0.05 for (4q) channel. The above values are in a good agreement with unity and thus the errors are correctly estimated.

## 4 Summary

A model independent method for measuring the quantities influenced by bin-to-bin and inside bin correlations is described. A package, as an addition to HBOOK, was written to support this new functionality. The method was tested using the simulated WW events.