text:
# Quantum discreteness and fundamental 1/f noise in tunnel junctions, nano-composites and other many-electron systems

## I Introduction

1. The low-frequency flicker noise (1/f-noise) known in a large variety of systems is actual problem of applied and theoretical physics \[1-5\]. This noise manifests itself as fluctuations in rates of relaxation and transport processes. In electronics thats are conductance fluctuations in practice often representing most bad type of noise. As a rule, these fluctuations are much higher sensitive to a structure of materials and external influences than conductivity itself and hence may bring a delicate information about mechanisms of charge transport. But the problem is that there is no adequate comprehension of the problem.

For many years theory tries to reduce flicker noise to large temporal scales (the concepts of hierarchycal kinetics  and self-organized criticality  involve also large spatial scales). Usually, flicker noise is thought be result of slow thermally activated fluctuations in structural disorder or in occupancy of electron traps, etc. \[1-3,5-15,29,37\]. Spectrum 1/f is composed by Lorentzians corresponding to ”fluctuators” with different activation energies, under assumption that sufficiently broad and uniform distribution of the energies takes place \[1,5-15,37\]. But in fact either qualitative origin of fluctuators remains unclear or quantitative treatment of experiments remains made at a stretch. Many important facts in no way could be kept within this theory, for instance, 1/f mobility fluctuations of carriers in clear intrinsic semiconductors , or 1/f-noise in liquid metals .

At present, the noise in systems with bad (narrow-band, hopping, tunnel, percolative) conductivity is under common attention, namely, in strongly doped or defected or amorphous semiconductors \[5,7,11-14\] and metals , in oxides , in manganites with colossal magneto-resistance \[16-20\], etc. Since in such systems long-range Coulombian forces can be significant, the charge redistributions were suggested to the role of slow fluctuators \[5,11-15\]. But again it is hard to explain the absence of saturation of real 1/f spectra at low frequencies . Perhaply, it was the true prophecy  that the game ”1/f from Lorentzians” may turn to a disaster for 1/f-noise theory.

2. With no doubts of the importance of thermally activated ”slow” fluctuations, it would be reasonable to separate them from the fundamental 1/f-noise whose direct source is nothing but the same ”fast” kinetic events (interactions and collisions of particles and quanta) which produce the electric resistivity itself \[4,23-27\].

Let us recall that generally the phenomenon of relaxation and irreversibility (in particular resistivity) in one or another dynamical system is equivalent to its ability to constantly forget its own history. Most easy the oblivion spreads to the number of kinetic events which took place in the past (for instance, more than a few inelastic free-path times ago). Hence, the system loses both stimulus and means for keeping some definite ”number (or probability) of events per unit time” (more accurately, only the ratio of numbers of time-reversed events is under thermodynamical control, but not their sum or difference) . Therefore, the fluctuations in ”thickness” (the term from ) of kinetic events and correspondingly in the rates of relaxation and transport do not cause some back compensating reaction, hence, they themselves do not relax and have no upper time scale.

Such statistical behaviour results just in scale-invariant non-saturating spectrum 1/f (for details see ). Thus, the extremes converge: long-living correlations associated with 1/f spectrum reflect not a long memory but, in opposite, indifference with respect to the past (it is usefull, as recommended by Krylov , to get rid of the prejudice that any statistical correlation testifies some actual causal correlation).

The kinetic theory throws out this child at the same time as turbid water if postulating quite certain ”probabilities per unit time” (i.e. collision integrals etc.). But if the path from statistical mechanics of a gas to gas kinetics is overcomed without such an ansatz  then it leads to 1/f fluctuations of diffusivities and mobilities of molecules. The origin of these fluctuations is not long relaxation times (surely absent here) but merely indifference of the system with respect to impact parameters of previous collisions (and thus to their effective cross-section). The flicker fluctuations in dissipation and in light scattering in quartz also can be explained as statistical property of phonon kinetics . The latter does not reduce to a certain marginal collision integrals, because phononic kinetic events are strongly entangled in space and time and parametrically interplay one with another.

3. Hence, the general principle of ”1/f-noise from lossing memory” realizes by different concrete ways. In this work (representing development of the preliminary study ) we argue that in many-electron systems it may realize like in phononic systems by means of mutual entangling kinetic events, i.e. electron transfers and jumps. We will see that flicker conductance fluctuations naturally come to light if a theory takes into account factual finite duration of kinetic events and besides real discreteness of electron energy states.

Because of its finite duration any one-electron transfer becomes the part of a complex many-particle process. If cutting off this part across bosonic lines one gets that quantum amplitude of given transfer evolves under influence by effective time-varying fields representing other components of the process. For instance, currently realizing electron jumps through a tunnel junction, as well as thermal charge motion in its sides, induce fast voltage fluctuations at the junction. In its turn, this voltage noise causes random phase changes in the increments of quantum amplitude of the maturing transfer. Such a picture with mutual influence of one-electron events mediated by many soft photons was considered in the theory of Coulomb blockade and low-temperature anomalies of current-voltage cvharacteristics (CVC) of small tunnel junctions \[31-33\]. The mathematically close interpretation of many-phonon processes was used in the theory of mobility of strongly coupled polarons .

Obviously, this picture must foresee not only renormalization of transport characteristics but in addition their specific fluctuations. If quantum amplitude of any electron transfer becomes a functional of enviromental noise (electric, elastic, exchange fields, etc.) then corresponding quantum probabilities become random values, and eventually the conductance as well. As far as we know, previously this phenomenon was not under consideration. It must be well expressed if correlation time of the noise is small as compared with typical duration (expectation time) of electron jumps. We will demonstrate that such a situation is easily realizable, and the result of fast enviromental noise, regardless of its concrete mechanisms and statistics, is just flicker scaleless fluctuations of transfer rates (probabilities) and conductance (therefore, flicker noise what accompanies hopping conductivity does not need in slowness of charge redistributions, instead such thats are sufficient which are as fast as diffusion of carriers).

As the experiments prompt, in relative units conductance fluctuations are inversely proportional to a number of movable carriers in a noisy region, that is to its volume. From here the principal statement does follow: flicker noise can not be cathed by a theory which underestimates role of the discreteness and too hurries to turn sums over energy states into integrals.

4. For visuality, we will concentrate at ”ideal” tunnel junction. Previously various aspects of many-particle processes and the discreteness in tunnel junctions were touched (see e.g.  and refs. therein), but with no respect to the effect of our interest. This effect, i.e. conductance fluctuations caused by fast thermal noise, principally differs of the low-temperature conductance fluctuations in tunnel micro-junctions due to quasi-elastic electron scattering by static structural disorder. . The 1/f-noise observed in ”material” junctions  usually is attributed to structural fluctuators, so called two-level systems. The possibilities of this theory (developing general approach ) were analysed in . This theory surely explains Lorentzian contributions to the whole noise but avoids an interpretation of the ”residual” 1/f component observed in . Unlikely one might categorically say that this component (and then almost all high-temperature noise possessing smooth 1/f spectrum) can not have another origin. Such a suspicion is supported by the experiment  as discussed below.

## II Time scales of tunnelling

1. If a voltage $`U<T/e`$ ($`T`$ is temperature) is applied to a tunnel junction then the mean charge transported during time $`\mathrm{\Delta }t`$ and corresponding conductance can be expressed as

$$\mathrm{\Delta }Q=e\frac{Ue}{\delta E}\frac{\mathrm{\Delta }t}{\tau _t},G=\frac{\mathrm{\Delta }Q}{U\mathrm{\Delta }t}=e^2\nu \gamma ,\nu =\frac{1}{\delta E},\gamma =\frac{1}{\tau _t}$$
(1)
Here $`\delta E`$ is mean spacing of electron energy levels, thus $`\nu `$ is density of states; $`Ue/\delta E`$ is the number of ”active” levels effectively contributing to electron transport between outward leads; $`\tau _t`$ is mean transmission time required for one electron jump from a given level at one side to wherever at opposite side, in other words, this is the time of accumulation of quantum jump probability to a value $`1`$ (though quantum jumps are momentary, the moment is unpredictable and may be expected up to $`\tau _t`$ or longer); $`\gamma `$ plays the role of mean jump probability per unit time (thus $`e/\tau _t`$ is the mean current per active level).

Any real junction has a finite capacity $`C`$ and thus characteristic time $`\tau _cRCC/G`$ . Its physical sense is the relaxation time of charge captured by $`C`$ and the correlation time of thermal voltage fluctuations at the junction. Compare the defined time scales and let be convinced of that tunneling is very long event, that is

$$\frac{\tau _t}{\tau _c}=\frac{e^2\nu }{C}=\frac{E_c}{\delta E}>>1$$
(2)
even if Coulombian effects are weak in the trivial sense $`E_c<<T`$ . For certainty, consider flat metallic electrodes with thickness $`w`$ , correspondingly flat barrier with thickness $`d`$ and typical dielectricity $`ϵ20`$ . Then taking from textbooks standard estimates for capacity and Fermi energy and velocity one obtains

$`{\displaystyle \frac{e^2\nu }{C}}{\displaystyle \frac{e^2}{\mathrm{}v_F}}{\displaystyle \frac{4\pi dw}{ϵa^2}}{\displaystyle \frac{dw}{a^2}},`$where $`a`$ is atomic size, $`a310^8cm`$ .

This expression clearly shows that inequality (2) is always satisfied. Consequently, while tunneling an electron has time to virtually feel many variations of the inter-side voltage noise, $`u(t)`$ , produced by many other jumps in both directions. At one-electron language, this means that quantum-mechanical transfer probabilities behaves as random processes. From the point of view of rigorous field theory of many particles, the description of corresponding excess fluctuations of transport current would need in at least four-particle Green functions . Since the necessary formal technique is still absent, we will try to understand the essence of the matter by means of archaic tunnel Hamiltonian method and quantum-mechanical non-stationary perturbation theory. We confine ourselves by the simplest variant: the sides are identical, all tunnel matrix elements $`g_{kq}g`$ are approximately equal and transparency is small.

2. Let us recall the underlying reasons of standard tunnel Hamiltonian scheme (for a time forgetting voltage noise). Let $`p_{kq}(\mathrm{\Delta }t,U)`$ be the probability of electron transfer from a state ”k” at left side to a state ”q” at right side during an observation time $`\mathrm{\Delta }t`$ (or probability of the time-reversed event). The probability of electron jump from left level ”k” to wherever at right-hand side is given by sum of these elementary probabilities, and jump duration is the inverse jump probability per unit time:

$`\gamma =\tau _t^1=p_k(\mathrm{\Delta }t,U)/\mathrm{\Delta }t,p_k(\mathrm{\Delta }t,U){\displaystyle \underset{q}{}}p_{kq}(\mathrm{\Delta }t,U)`$(by above assumptions, these values are weakly dependent on ”k”).

Further, the standard scheme recommends (as usually when constructing kinetics) to resort to the Fermi ”golden rule”,

$`p_{kq}(\mathrm{\Delta }t,U)2\pi g_{kq}^2\delta (E_{kq})\mathrm{\Delta }t/\mathrm{},E_{kq}E_q^+E_k+eU`$(plus relates to right-hand side). This recipe ensures the time-linear increase of the jump probability and thus the existence of certain rate of jumps,

$$\gamma =\tau _t^12\pi g^2\nu /\mathrm{}$$
(3)
The using golden rule (in other words, continuous spectrum approximation) implies that the duration of watching on the evolution of quantum amplitudes, $`\mathrm{\Delta }t`$ , what is sufficient for adequate evaluation of jump probabilities is in the frames $`\mathrm{}/T<<\mathrm{\Delta }t<<\tau _g`$ . Here the characteristic time is introduced, $`\tau _g=2\pi \mathrm{}/\delta E`$ , involved by discreteness of energy spectrum.

Clearly, if one wants to account for the influence of voltage noise $`u(t)`$ then the adequate watching time should be much greater than its correlation time, and in view of (2) it is desirable to continue the watching time up to a value of order of factual duration of jumps. Hence, the standard scheme needs in the condition $`\tau _g>\tau _t`$ . But the Eqs.1 and 3 lead to the relation

$$\frac{\tau _t}{\tau _g}=\frac{R}{R_0}\left(\frac{\delta E}{2\pi g}\right)^2,R_0\frac{2\pi \mathrm{}}{e^2},$$
(4)
which shows that in a bad transparent junction just the opposite situation takes place,

$$\tau _t/\tau _g>>1$$
(5)
Thus it comes out that the golden rule is inapplicable, and one runs into the problem: ”probabilities per unit time” remain uncertain.

We may take in mind the case of well expressed discreteness (5), since this is principally most interesting case and besides at $`R<R_0`$ compatibility of the tunnel Hamiltonian method and perturbation theory could be under question  (nevetheless, junctions with $`R<R_0`$ also will be described if thats can be represented as parallel sum of autonomous high-resistance junctions). Under condition (5), the perturbation theory is undoubtly applicable even in presence of voltage noise.

## III Fluctuations of probabilities

1. Consider quantum transfers influenced by time-varying potential difference, $`u(t)`$ , between initial and final states. According to modern theory of quantum chaos, stochastic behaviour is typical for quantum systems regardless of energy discreteness . Therefore, we will treat $`u(t)`$ like classical random process (though remembering that in rigorous theory $`u(t)`$ is an operator entangled with particle operators). At $`\mathrm{\Delta }t\tau _t`$ , the perturbation theory is sufficient to approximately solve the Shroedinger equation for wave function of tunneling particle. The result reads as

$$p_{kq}\left|A_{kq}\right|^2,A_{kq}\frac{g_{kq}}{\mathrm{}}_0^{\mathrm{\Delta }t}\mathrm{exp}(iE_{kq}t/\mathrm{})Z(t)𝑑t,Z(t)=\mathrm{exp}[i\phi (t)],$$
(6)
where the diffusively accumulating random phase is introduced,

$`\phi (t)={\displaystyle \frac{e}{\mathrm{}}}{\displaystyle _0^t}u(t^{})𝑑t^{}`$At $`u(t)=0`$ these formulas reduce to usual ones, otherwise they describe chaotic parametrical excitation or damping of probabilities by the random phase.

Let us introduce phase correlation function, corresponding coherence time of quantum amplitudes and ”coherence region” by formulas

$$K(t_1t_2)=Z(t_1)Z^{}(t_2),\tau _{coh}=_0^{\mathrm{}}\left|K(\tau )\right|𝑑\tau ,\mathrm{\Delta }E=2\pi \mathrm{}/\tau _{coh},$$
(7)
where angle brackets denote the averaging with respect to $`u(t)`$ . Generally, because of multiplicative structure of $`Z(t)`$ calculations of even simplest statistical characteristics of the solution require a vast statistical information about $`u(t)`$ . But if the coherence time is much smaller than observation time then factor $`Z(t)`$ under integral in (6) acts as complex fast (”white”) noise. Consequently, the transfer amplitudes $`A_{kq}`$ behave as (complex) Brownian walks, and the minimal information presented by characteristics (7) is already sufficient.

2. Discuss briefly the coherence time, assuming for simplicity that shunting influence by external circuit on $`\tau _c`$ is negligible. Note that $`K(t)`$ represents characteristic function of the phase, in the sense of probability theory. This function can be easy written at $`E_c<<T`$ if one treats voltage noise as Gaussian random process. Corresponding calculations give

$`\tau _{coh}(\mathrm{}/e)(C/T)^{1/2}<<\tau _g`$(perhaply, that is minimally possible value of $`\tau _{coh}`$ ). At $`E_cT`$ , charge quantization is essential and it is better to treat $`u(t)`$ like triple-valued random process, $`u(t)=0,\pm e/C`$ . Corresponding analysis finishes with simple estimate

$`\tau _{coh}\tau _c`$In this case coherence time may be comparable with $`\tau _g`$ .

Hence, there are all grounds to consider $`Z(t)`$ as white noise and consequently the amplitudes as Brownian walks. This allows to write

$$p_{kq}^2=\left|A_{kq}\right|^42\left|A_{kq}\right|^2^2=2p_{kq}^2,p_{kq},p_{kq}p_{kq}^2$$
(8)
We use Malakhov cumulant brackets,

$`x,yxyxy`$At first, we see that after a time longer than coherence time the quantum amplitudes and probabilities acquire 100-percent uncertainty. Secondly, the mean probabilities grow linearly with time:

$$p_{kq}\mathrm{\Delta }t(g_{kq}/\mathrm{})^2K(\tau )\mathrm{exp}(iE_{kq}\tau /\mathrm{})𝑑\tau \mathrm{\Delta }t,$$
(9)
Now, the fixed scale $`\mathrm{\Delta }E`$ instead of $`2\pi \mathrm{}/\mathrm{\Delta }t`$ determines energy region available for transfers from a given level.

3. We are most interested in the summary jump probability. It can be represented as

$$p_k\underset{q}{}p_{kq}=_0^{\mathrm{\Delta }t}\mathrm{\Gamma }_k(t_1t_2)Z(t_1)Z^{}(t_2)𝑑t_1𝑑t_2,$$
(10)
with the integral kernel

$`\mathrm{\Gamma }_k(\tau )={\displaystyle \underset{q}{}}({\displaystyle \frac{g_{kq}}{\mathrm{}}})^2\mathrm{exp}(i\tau E_{kq}/\mathrm{})`$Analytical properties of such kernels dictated by quantum discreteness play a key role in the theory. In the continuous spectrum approximation this kernel turns into a function which quickly (integrably) and irreversibly decays to zero, for instance, delta-function. But in reality this kernel is quite non-local and never decays (instead, sometimes it returns to a value of order of its value at $`\tau =0`$). In principle, this property is nothing but reflection of the unitarity of quantum dynamics. If take, for visuality, the equidistant spectrum at side accepting jumps, $`E_{kq}=n\delta E+\epsilon _k`$ , with integer $`n`$ , then

$$\mathrm{\Gamma }_k(\tau )=\frac{1}{\tau _t}\mathrm{exp}(i\epsilon _k\tau /\mathrm{})\underset{n}{}\delta (\tau n\tau _g)$$
(11)
Now the third point is obvious: if $`\mathrm{\Delta }E>\delta E`$ then the mean jump probability $`p_k`$ $`\mathrm{\Delta }t/\tau _t`$ practically coinsides with what is used in kinetics. Hence, due to the noise transfers into discrete spectrum realize as successfully as to continuous spectrum.

But, of course, the noise can not make the ”probabilities per unit time” be better certain. The most inportant point is that due to the discreteness the summary jump probabilities also have non-zero fluctuations. With accounting for the white-noise character of $`Z(t)`$ , variancies of jump probabilities can be expressed as

$$p_k,p_k\mathrm{}_0^{\mathrm{\Delta }t}\mathrm{\Gamma }_k(t_1t_2)\mathrm{\Gamma }_k(t_3t_4)K(t_1t_4)K(t_3t_2)𝑑t_1\mathrm{}𝑑t_4$$
(12)
Here under condition $`\mathrm{\Delta }E>\delta E`$ only regions $`t_1t_4,`$ $`t_3t_2`$ are significant but many delta functions from (11) contribute. The estimate of this integral yields

$$p_k,p_k\frac{\mathrm{\Delta }t^2}{\tau _t^2\tau _g}\left|K(\tau )\right|^2𝑑\tau \frac{\tau _{coh}}{\tau _g}p_k^2=\frac{\delta E}{\mathrm{\Delta }E}p_k^2$$
(13)
(we took into account that ”width” of delta functions determined by inverse energy band is wittingly smaller than $`\tau _{coh}`$ ). Notice that conributions to $`p_k(\mathrm{\Delta }t,U)`$ from different parts of the observation time are completely statistically correlated although produced by independent pieces of realization $`u(t)`$ . Such unusual transport statistics were considered in . Here, it is due to that transport is governed by subsequent time summation of the amplitudes, not probabilities.

Below mutual correlation between fluctuations of probabilities of jumps from different levels will play an important role. It reads as

$$p_{k_1},p_{k_2}p_{k_1}p_{k_2}\frac{\delta E}{\mathrm{\Delta }E}S(E_{k_1}E_{k_2}),$$
(14)
where function

$`S(E)={\displaystyle \mathrm{exp}(iE\tau /\mathrm{})\left|K(\tau )\right|^2𝑑\tau \left[\left|K(\tau )\right|^2𝑑\tau \right]^1}`$describes its dependence on energy distance between two levels, $`E`$ . According to these formulas, the ”uncertainty principle” takes place: with increasing $`\mathrm{\Delta }E`$ all fluctuations of jump probabilities decrease but instead their correlations spread to more and more wide energy distancies. The rates of charge injection from levels belonging to the same ”coherence region” $`\mathrm{\Delta }E`$ are closely correlated.

At $`\mathrm{\Delta }E<\delta E`$ fluctuations of jump probabilities may increase up to 100-percent value or even higher. In this extremal case quantitative estimates are sensible to details (commensurability etc.) of energy spectra, therefore one needs in some statistics of energy levels. Besides, in this case the noise induced renormalization of mean jump probabilities and consequently average transport current and CVC becomes significant. All the above results can be extended to when inequality (5) is not satisfied.

## IV Conductance fluctuations

1. Consider fluctuations of the charge transported through the junction between outward leads. We will confine ourselves by the case of not too low temperature, $`T>>\delta E`$ , and not too high external voltage, $`U<T/e`$ .

Below let $`\mathrm{\Delta }Q`$ denote the random value. It consists of two parts,

$`\mathrm{\Delta }Q=\mathrm{\Delta }Q_{th}+\mathrm{\Delta }Q_{ex}`$Here first term is contribution of fast thermal (shot) current noise caused by uncertainty of jump instants. This part contributes to variance of $`\mathrm{\Delta }Q`$ even in equilibrium (at $`U=0`$ ) and can be estimated with the help of fluctuation-dissipation theorem (Nyquist formula),

$`\mathrm{\Delta }Q_{th}^22TG\mathrm{\Delta }t`$We are more interested in the second term which besides mean value of transport current includes also excess transport fluctuations caused by the above discussed uncertainty of jump probabilities.

Within statistical language, $`\mathrm{\Delta }Q_{ex}`$ may be defined as conditional average value of $`\mathrm{\Delta }Q`$ under fixed $`p_{kq}`$ . It is clear, from simple thermodynamical reasonings, that this conditional average has the same sign as $`U`$ and disappears at $`U=0`$ . Therefore, it can be represented as the sum of definitely directed jumps (for instance, from left to right at $`U>0`$ ). This circumstance allows us to concentrate at the excess contribution without special analysing correlations between oppositely directed electron transfers. These correlations create no obstacles to arbitrary low-frequency transport fluctuations (distributed at times longer $`\tau _c`$ ) since such fluctuations do not touch inner statistical balance in the junction . Thus we can assume that all the correlations are already included into statistics of voltage noise $`u(t)`$ . All the more this is reasonable in view of the appointed fact that details of the noise statistics are unimportant.

2. Therefore, let us attribute formulas of previous Section to the probabilities of one-directed excess jumps responsible for transport. For a first, one may neglect correlations between random energy positions of ”active” levels which inject the excess jumps, and preliminarily perfom the averaging with respect to the postions. Then the start expression for the excess transport looks as

$$\mathrm{\Delta }Q_{ex}=e\underset{k}{}[f(E_k)f(E_k+eU)]p_k(\mathrm{\Delta }t,U),$$
(15)
with $`f(E)`$ being the Fermi distribution function. Below with no noticable losses we may put on $`p_k(\mathrm{\Delta }t,U)`$ $`p_k(\mathrm{\Delta }t,0)`$ . If one additionally neglected the randomness of probabbilities, then Eq.15 would reduce to the well known expression for mean transport current. It must be underlined that Eq.15 completely accounts for the Fermi statistics in both electrodes. The result of averaging (15) coinsides with Eq.1.

Consider fluctuations of the conductance $`G=\mathrm{\Delta }Q_{ex}/U\mathrm{\Delta }t`$ implied by Eq.15, omitting details of calculations. At small coherence time, what corresponds to ”large” junction, dense levels, and $`\mathrm{\Delta }E>\delta E`$ , with the help of formulas of previous section the variance of (15) can be transformed into

$$\mathrm{\Delta }Q_{ex},\mathrm{\Delta }Q_{ex}(GU\mathrm{\Delta }t)^2\frac{\delta E}{\mathrm{\Delta }E}W(E^{})W(E^{\prime \prime })S(E^{}E^{\prime \prime })𝑑E^{}𝑑E^{\prime \prime }$$
(16)
Here function

$`W(E)=[f(E)f(E+eU)]/eUf(E)/E`$has the sense of the ”one-particle” probability density distribution of energy of active levels. The Eq.16 yields rather universal estimate

$$\delta G^2\frac{G,G}{G^2}\frac{\delta E}{T},(\tau _{coh}<\tau _g)$$
(17)
At large coherence time, what corresponds to ”small” junction, rarefied levels, $`\mathrm{\Delta }E<\delta E`$ , and weakly correlated jumps, i.e. to the above mentioned extremal situation, so definite estimate is impossible, because conductance fluctuations are very sensible to statistics of energy levels, first of all to their commensurability. In this case ”everything is possible”,

$$\frac{\delta E}{T}<\delta G^2<1,(\tau _{coh}>\tau _g),$$
(18)
up to relative fluctuations of order of unit.

3. As it is seen, in general the discreteness directly serves as the measure of conductance fluctuations. At the same time, factor $`\mathrm{\Delta }E`$ which characterizes intensity of ”enviromental noise” does not contribute to (17). This is likely if the mean number of active levels, $`NeU/\delta E`$ , is not large and these levels make take arbitrary relative positions, in particular, all occur in one and the same coherence region $`\mathrm{\Delta }E`$ . But the latter case is impossible when $`N`$ exceeds total number of levels in such a region, $`\mathrm{\Delta }E/\delta E`$ . Consequently, at $`eU>\mathrm{\Delta }E`$ permissible energy distributions of active levels are more uniform, therefore jumps from them are less correlated. Under increase of applied voltage, this may result in supressing conductance fluctuations even at $`eU<<T`$ when average conductance is still more or less constant.

This effect is not scoped by the approximation (15). Formally, as the Eq.16 shows, the matter is that it corresponds to the factored pair (”two-particle”) energy distribution $`W(E^{})W(E^{\prime \prime })`$ . But it is not hard to improve our analysis. Let us enumerate active levels in the increasing energy order. The typical distance between two levels with numbers $`j>i`$ can not be smaller than $`(ji)\delta E`$ . Therefore, any pair of numbers implies its own pair distribution,

$$W_{ij}(E^{},E^{\prime \prime })W(E^{})W(E^{\prime \prime })\vartheta (\left|E^{\prime \prime }E^{}\right|\left|ji\right|\delta E)$$
(19)
wher $`\vartheta ()`$ is step-like function. Further, instead of (15), we should write

$$\mathrm{\Delta }Q_{ex}=e\underset{j=1}{\overset{n}{}}p_{k_j}(\mathrm{\Delta }t,U),$$
(20)
thus taking into consideration all possible random positions of active levels at energy scale (i.e. substituting average occupancies by factual unit). Here $`n`$ is random total number of active levels whose mean value equals to $`N`$ . Then calculate variance of this expression. Since at $`N1`$ the answer must coinside with previous result, we may put on $`N>>1`$ , $`n=N`$ . Performing the averaging with respect to both random energy positions and jump probabilities yields

$$\mathrm{\Delta }Q_{ex},\mathrm{\Delta }Q_{ex}e^2(\gamma \mathrm{\Delta }t)^2\frac{\delta E}{\mathrm{\Delta }E}\underset{i,j=1}{\overset{N}{}}W_{ij}(E^{},E^{\prime \prime })S(E^{}E^{\prime \prime })𝑑E^{}𝑑E^{\prime \prime }$$
(21)
After elementary manipulations we obtain instead of (17) the corrected estimate of conductance fluctuations,

$$\delta G^2\frac{\delta E}{T}D(eU),$$
(22)
where function

$$D(X)=\frac{1}{X}_0^X𝑑E_E^{\mathrm{}}S(E^{})𝑑E^{}\left[_0^{\mathrm{}}S(E)𝑑E\right]^1$$
(23)
describes their dependence on external voltage. For example, in the case of exponential phase correlation we have

$`K(\tau )=\mathrm{exp}(\left|\tau \right|/\tau _{coh}),D(X)={\displaystyle \frac{1}{X}}{\displaystyle _0^X}\left(1{\displaystyle \frac{2}{\pi }}\mathrm{arctan}{\displaystyle \frac{E}{\mathrm{\Delta }E}}\right)𝑑E`$The esimate (18) modifies in analogous fashion.

According to these formulas, at $`eU>\mathrm{\Delta }E`$ the fluctuations decrease approximately as $`U^1`$ . One can say, that the increase of applied voltage results in increase of effective number of statistically independent (energetical) tunneling channels, and relative conductance fluctuations decrease inversely proportionally to the number of channels. Thus the degree of discreteness, $`\delta E`$ , as well as the degree of voltage noise, $`\mathrm{\Delta }E`$ , both have direct reflection in ”current-noise characteristics” if not in current-voltage characteristics. The clear separation of duties of these two parameters gives hope that their significance will be maintained in a more rigorous theory.

4. What is for the transparency of junction, in relative units it slips out from the estimates. In this sense it can not be treated as ”small expansion parameter”, as well as in the sense of relation between excess transport noise and thermal (shot) noise. Though the shot noise contribution to variance of the transported charge is proportional to the transparency, while excess contribution is proportional to it squared, the latter grows as the observation time squared, while the former grows linearly with time. Therefore, excess transport noise inevitably dominates at large times and low frequencies. It is easy to verify that at $`eUT`$ this takes place already after a time of order of one jump duration, $`\tau _t`$ . In fact, in experiments and practice excess noise noticably dominates at low frequencies at any realistic operational currents. Besides, naively formal reasonings on ”powers of transparency” are wrong in view of that the random phase $`\phi (t)`$ involves in some complicated non-perturbative way all orders of the transparency.

## V Discussion and comparison with experiments

1. Experimental confirmation of the theory can be found in the work  devoted to 1/f noise in films of the cermet (nano-composite) formed by $`Ni`$ nano-particles in matrix $`Al_2O_3`$ . In this system the parameters of a typical elementary tunnel junction between neighbouring particles were $`\delta E0.2`$meV , $`d2`$nm, $`C510^6`$cm, $`E_CT`$ (at room temperature), $`R30`$MOhm, what means $`\tau _g310^{11}`$s, $`\tau _c1.510^{10}`$s and $`\tau _t310^8`$s. In this system very intensive conductance fluctuastions were observed with relative spectral power $`S_{\delta G}(f)\alpha /N_gf`$, where $`\alpha 610^3`$ , and $`N_g`$ is number of metal granules in the sample. Since $`E_cT`$ , then $`N_g`$ approximately coinsides with the number of movable (simultaneously transported) electrons in the sample , therefore it was almost ”standard” 1/f-noise with ”classical” Hooge constant \[1-4\].

These data correspond to the noise $`S_{\delta G}(f)\alpha /f`$ in a separate elementary junction. Since both the inequalities (2) and (5) are well satisfied, we may assume that this noise is due to the mechanism under consideration. For stationary 1/f noise, the above calculated conductance variance related to a definite smoothing time can be connected with 1/f spectrum as

$`\delta G^2fS_{\delta G}(f)\mathrm{ln}(\mathrm{\Delta }t/\tau _c)\alpha \mathrm{ln}(\mathrm{\Delta }t/\tau _c)`$At $`\mathrm{\Delta }t\tau _t`$ this gives the value $`0.03`$ . Formula (17) at $`\delta E=0.2`$meV and room temperature yields $`0.008`$ . This is good agreement, in view of that strictly speaking this case must be subjected to formula (18). Indeed, due to above estimates, here $`\tau _{coh}\tau _c>\tau _g`$ , hence $`\mathrm{\Delta }E`$ is of order of $`\delta E`$ or even smaller what corresponds to the extremal situation.

The remarkable observation of  is 1/f-noise sensibility to the discreteness of electron energy spectra in metal granules. When applied voltage per elementary junction exceeds $`\delta E/e`$ , i.e. when more than one level in a granule is active, relative 1/f noise intensity decreases as inverse voltage, although average CVC is still Ohmic up to $`T/\delta E100`$ times larger voltages. In view of $`\mathrm{\Delta }E\delta E`$ , it is clear that in respect to this effect our theory is in full agreement with the experiment.

2. In case of cermet the discreteness of levels is directly determined by a volume of junction sides, i.e. metal nano-particles. Obviously, in a junction with massive sides (electrodes) the quantity $`\delta E`$ also must be determined by volume of the space region which is physically available for electron jumps, that is by geometry of junction and processes of electron interactions and scatterings in electrodes. If temperature is not too low then one has no many variants. The available region must be restricted by junction area, $`A`$ , and non-elastic free path in sides, $`\lambda `$ . In other words, this is the region in which level spacing has the order of broadening of levels due to non-elastic relaxation. Of course, now it would be more precise to say not about literally levels but about statistics of energy repulsion of spatially close electron states , which has a sense and acts also at non-zero temperatures . Thus the term discreteness serves as synonym of uncertainty: nuances of electron states depend on the ”enviromental noise” and can not be controlled with precision better than $`\delta E`$ .

3. Consequently, in case of massive metallic electrodes one can estimate $`\delta EE_Fa^3/A\lambda `$ , with $`E_F`$ being Fermi energy and $`a`$ atomic size. Then, if roughly relate $`\lambda `$ to conductivity in electrodes, $`\sigma \lambda /a^2R_0`$ , the estimate (17) transforms into

$$\delta G^2\frac{E_F}{T}\frac{a^2}{A}\frac{\sigma _{min}}{\sigma }$$
(24)
where $`\sigma _{min}(aR_0)^1`$ is minimal metallic conductivity. In particular, let the metal be so clear that the electron-phonon mechanism of relaxation is dominating. Then, according to ,

$`\sigma (\mathrm{}e^2n_e/m^{}T)(T_D/T)^4`$(all notations are standard), and one may expect that at temperatures lower than Debay temperature conductance fluctuations are proportional to $`T^4`$ .

If say about possible orders of value only, then at standard $`E_F`$ the Eq.24 yields

$$\delta G^2(a^2/A)(T/T_D)^4$$
(25)
From the other hand, left side of (25) may be connected, like it was made above, with spectral power of 1/f-noise. Then in the case of micro-junction with area $`10^9cm^2`$ investigated in , under assumption $`TT_D`$ we obtain the estimate $`fS_{\delta G}(f)10^7`$ . This value is in agrtement with meausurements in  at 260 K . Besides, the estimate (25) naturally explain sharp (approximately by two orders of magnitude) increase of 1/f-noise under increasing temperature from 100 K to 300 K as was observed in . To explain this fact the authors supposed sharp step in the activation energy distribution of fluctuators (although this assumption contradicts to uniformity of the distribution wished for composing spectrum 1/f). From our look, likely in  two types of low-frequency noise were observed: one has structural origin and dominates below 100 K, while another is just what was considered and dominates at higher temperatures.

4. In this connection, it is useful to notice that generally summation of structural fluctuations and fundamental flicker noise under our attention may not obey trivial adding rules . Clear 1/f spectrum corresponds to fast forgetting the past , while entangling between the fundamental noise and actually slow relaxation processes (for instance, heat transport) may lead to somehow deformed spectrum (for instanse, $`1/f^\gamma `$ with $`\gamma >1`$ ). In particular, structural defects can affect 1/f noise although not being its origin. For example, elastic electron scattering by impurities supresses mobility 1/f noise . Experiments with metallic films show that 1/f noise there is strongly influenced by vacancies . But this effect can not be explained in terms of thermally activated diffusion of vacancies. Indeed, unlikely separate diffusive jumps of vacancies have very large variance of activation energies. What is for the diffusion as the whole it is described by the same equations as temperature diffusion, but the latter far ago was removed as hypothetical source of 1/f spectrum \[1-4\]. Perhaply, vacancies promote non-elastic scattering of electrons and in this way become mediator for the ”1/f-noise from lossing memory”.

## VI Conclusion

As it was demonstrated, the careful treatment of quantum discreteness displays that amplitudes of elementary quantum transfers being sensitive to surrounding noise of the whole system behave like random Brownian walks what leads eventually to low-frequency flicker fluctuations of transport rates. These fluctuations have no characteristic time scale and are indifferent to duration of observation and time averaging.

Formally, in our approximation restricted by time intervals $`\mathrm{\Delta }t\tau _t`$ we obtained fluctuations with perfectly undecaying correlations corresponding to spectrum $`\delta (f)`$ . This resembles static universal conductance fluctuations , although by essence thats are quite different things: one starts from inelastic relaxation and decoherence while another finishes at them. We expect that rigorous many-particle theory will extend our results (and first of all the scaleless character of conductance fluctuations) to arbitrary long time intervals. In this theory quantities kindred to above $`p_k`$ will be measure of number of electron jumps but will remain proportional to both time and transparency (highest orders of transparency will be absorbed by the random phase, - see notion at end of Sec.3), while spectrum $`\delta (f)`$ must transform to a spectrum $`1/f`$ with the same frequency dimensionality. Formal grounds of such a confidence will be considered elsewhere. Its principal motivation is that the only source of our results is the trivial quantum mechanical rule that even at presence of noise and decoherence any kinetics is governed by summation or other play of amplitudes, not of intermediate probabilities.

REFERENCES

1. P.Dutta, P.Horn, Rev.Mod.Phys., 53(3), 497 (1981).

2. F.N.Hooge, T.G.M.Kleinpenning, L.K.J.Vandamme, Rep.Prog.Phys., 44, 481 (1981).

3. M.B.Weissman, Rev.Mod.Phys., 60, 537 (1988).

4. G.N.Bochkov, Yu.E.Kuzovlev, Sov.Phys.-Usp., 26, 829 (1983).

5. M.B.Weissman, Rev.Mod.Phys., 65, 829 (1993).

6. B.Raquet, J.M.D.Coey, S.Wirth and S. Von Molnar, Phys.Rev., B59, 12435 (1999).

7. J.G.Massey and M.Lee, Phys.Rev.Lett., 79, 3986 (1997).

8. M.J.C. van den Homberg, A.H.Verbruggen et al., Phys.Rev., B 57, 53 (1998).

9. A.Ghosh, A.K.Raychaudhuri et al., Phys.Rev., B 58, R14666 (1998).

10 X.Y.Chen, P.M.Koenrad, F.N.Hooge et al., Phys.Rev., B 55, 5290 (1997).

11. G.M.Khera and J.Kakalios, Phys.Rev., B 56, 1918 (1997).

12. M.Gunes, R.E.Johanson and S.O.Kasap, Phys.Rev., B 60, 1477 (1999).

13. K.M.Abkemeier, Phys.Rev., B 55, 7005 (1997).

14. G.Snyder, M.B.Weisman and H.T.Hardner, Phys.Rev., B 56, 9205 (1997).

15. V.I.Kozub, Solid State Commun., 97, 843 (1996).

16. A.Lisauskas, S.I.Khartsev and A.M.Grishin, J.Low.Temp.Phys., MOS-99 Proceedings.

17. A.Lisauskas, S.I.Khartsev, A.M.Grishin and V.Palenskis, Mat.Res.Soc.Proc., Spring-99 Meeting.

18. V.Podzorov, M.Uehara, M.E.Gershenson and S.-W.Cheong, lanl arXiv cond-mat/9912064.

19. M.Viret, L.Ranno and J.M.D.Coey, Phys.Rev., B 55, 8067 (1997).

20. J.M.D.Coey, M.Viret and S. von Molnar, Adv. Phys., 48, 167 (1999).

21. P.Bak. Self-organized criticality: why nature is complex. Springer, N-Y, 1996.

22. J.L.Tandon and H.P.Bilger, J.Appl.Phys., 47, 1697 (1976).

23. Yu.E.Kuzovlev, G.N.Bochkov, Izv. VUZov. Radiofizika, 26, 310 (1983); 27, 1151 (1984) (transl. in Radiophysics and Quantum Electronics, N.3, 1983, and N.9, 1984).

24. Yu.E.Kuzovlev, Sov.Phys.-JETP, 67, N.12, 2469 (1988).

25. Yu.E.Kuzovlev, JETP, 84, N.6, 1138 (1997).

26. Yu.E.Kuzovlev, Phys.Lett., A 194, 285 (1994).

27. Yu.E.Kuzovlev, lanl arXiv cond-mat/9903350.

28. N.S.Krylov. Works on the foundations of statistical physics. Princeton, 1979.

29. J.P.Hessling, Yu.Galperin, M.Jonson, R.I.Shekhter and A.M.Zagoskin, Aplied Physics Report 94-38.

30. Yu.M.Galperin, V.G.Karpov, V.I.Kozub, Zh. Eksp. Teor. Fiz., 95, 1123 (1989).

31. Yu.V.Nazarov, Zh. Eksp. Teor. Fiz., 95, 975 (1989).

32. M.H.Devoret, D.Esteve, H.Grabert et al., Phys.Rev.Lett., 64, 1824 (1990).

33. S.M.Girvin, L.I.Glazman, M.Jonson et al., Phys.Rev.Lett., 64, 3183 (1990).

34. I.G.Lang, Yu.A.Firsov, Zh. Eksp. Teor. Fiz., 43, 1843 (1962).

35. Yu.A.Genenko and Yu.M.Ivanchenko, Phys.Lett., 126, 201 (1987).

36. Yu.M.Ivanchenko, Yu.V.Medevedev, Fiz. Nizk. Temp., 2, 142 (1976).

37. C.T.Rogers and R.A.Buhrman, Phys.Rev.Lett., 53, 1272 (1984).

38. Yu.A.Genenko, Yu.M.Ivanchenko, Teor. Mat. Fiz., 69, 142 (1986).

39. G.Casati and B.Chirikov. Fluctuations in quantum chaos. Preprint, Budker Inst. of Nuclear Physics SB RAS, 1993.

40. C.W.J.Beenakker, Rev.Mod.Phys., 69, N 3, 731 (1997).

41. J.V.Mantese, W.I.Goldburg, D.H.Darling et al., Solid State Commun., 37, 353 (1981).

42. G.P.Zhigalski, Usp. Fiz. Nauk, 167, 623 (1997).

43. E.M.Lifshitz, L.P.Pitaevski. Physical kinetics. Moscow, 1974.

44. A. van Oudenaarden, M.H.Devoret et al., Phys.Rev.Lett., 78, 3539 (1997).

45. Yu.V.Nazarov, Sov.Phys.-JETP, 71, 171 (1990).

46. Yu.E.Kuzovlev, lanl arXiv cond-mat/0004398.


format: text

text:
# Solitons in a one-dimensional interacting Bose-Einstein system

## I Introduction

There has been a number of theoretical studies on Bose-Einstein condensation (BEC) in a trap that is effectively lower dimensional . A three-dimensional (axially symmetric) anisotropic harmonic trap may be regarded as (quasi) two- or one-dimensional when $`\omega _{}\omega _x`$, or $`\omega _{}\omega _x`$ respectively, at temperatures much smaller than the larger of the oscillator energy gaps. Recently, Burger et al.  have experimentally created dark solitons in such a quasi one-dimensional system of dilute $`Rb^{87}`$ atoms. In this experiment, $`\omega _{}`$ was taken to be about $`30`$ times larger than $`\omega _x`$. After the creation of the soliton by phase-imprinting on the cigar-shaped condensate, the trap is switched off, and the soliton is observed to propagate along the $`x`$-direction. The $`s`$-wave inter-atomic interaction between $`Rb^{87}`$ atoms is repulsive, and stable dark solitons are expected to exist in one dimension, as seen in this experiment. We note, however, that for $`Rb^{85}`$ and $`Li^7`$, the interatomic interaction is attractive, and bright solitons may form. Such solitons have been studied for the quasi one-dimensional case by several groups , using the Gross-Pitaevskii (GP) mean-field theory . The standard GP theory replaces the bosonic field operator by a classical order parameter field $`\mathrm{\Phi }(𝐫,t)`$ in three dimensions, and the potential energy is given by a fourth-degree $`|\mathrm{\Phi }|^4`$ term in the energy-density functional. Perez-Garcia et al.  have shown how this may be reduced, for extreme anisotropy ($`\omega _{}\omega _x`$), to the quasi one-dimensional geometry with a fourth-degree potential energy term. It has been also shown that for the one-dimensional case with this attractive fourth-degree interaction potential, the solitons are dynamically stable ( for the repulsive case, this has been known for a long time ). One has also considered the limiting case for $`\omega _x=0`$, when there is no confinement in the $`x`$-direction, and studied the free propagation of the soliton.

A different approach has been taken recently by Kolomeisky et al. . Instead of starting with the usual three-dimensional GP energy functional and reducing it to the quasi one-dimensional case as described above, they have derived the GP functional for truly lower dimensions, $`d2`$, using the renormalization group approach . In particular, they proposed, for a repulsive interatomic interaction, a sixth-degree potential energy term replacing the usual fourth-degree term in the order parameter. With this modification they obtained a (dark) soliton solution in the repulsive case.

In this paper, we first discuss the solutions of the modified GP formalism in one dimension in the stationary frame. We show that the soliton solution of Kolomeisky et al.  may be reproduced in the stationary frame with a repulsive interaction. Since the number of particles (see the definition in the next section) is generally not finite unless regulated, the phase cannot be obtained by a Galilean boost even though the system has Galilean invariance to begin with.

Next we focus on an attractive sixth-degree interaction term in the energy functional for bosons. Unlike the repulsive sixth-degree interaction, the attractive interaction may not be obtainable as a modified GP functional in one dimension. However, the motivations for taking such a phenomenological attractive potential energy functional are two-fold. In three dimensions, using the standard fourth-degree attractive interaction energy term in the GP functional, it was predicted that the soliton-like solution would be stable only up to a limited number of bosons in the condensate , and this was verified experimentally . Of course, in one dimension with an attractive fourth-degree term, the soliton is variationally stable as we will show in Sec. III D below. We want to examine the analytical solutions for the attractive sixth-degree interaction in the one-dimensional case. Analytical solutions for a more general case which includes the sixth-degree interaction have been derived earlier, in a different context, by Pelinovsky et al. . They have also discussed numerically the stability of these solutions. Here we use the density functional approach to derive these analytical solutions in the stationary frame. We show that there exists a limitation on the particle number for stability. Since the particle number is finite in the attractive case, we may use the Galilean invariance to obtain the general solution in an arbitrary frame.

It should also be mentioned that the fourth- and sixth-degree interaction (in the field variable) in the energy functional has been used extensively in various applications of physics , with a recent one  examining the stability of the $`Li^7`$ condensate with an attractive two-body and a repulsive three-body contact interaction in three dimensions. More recently, the stability of an attractive BEC in the presence of a sixth-degree interaction with confinement has been analyzed by Zhang. The purely sixth-degree energy functional in one dimension, without the confining potential considered by Kolomeisky et al. , is special because it is scale-invariant, and it has a peculiar condition on the particle number for the soliton in the attractive case. We further show that there are soliton as well as stationary wave solutions if one uses the phenomenological model for an attractive interaction.

In Sec. II, we discuss the equations of motion in the stationary frame and their symmetries. In Sec. III, we obtain the solutions for repulsive and attractive interactions, and discuss the stability of the soliton solution in the attractive case. Sec. IV provides a general discussion of the stationary solutions which includes all the above solutions as special cases.

## II Equations of Motion

Our starting point is Eqs. (3) and (4) of Kolomeisky et al.  without a trap potential. The modified GP functional for a one-dimensional Bose-Einstein system is given by

$$E=\frac{\mathrm{}^2}{2m}𝑑x[\frac{d\mathrm{\Phi }}{dx}\frac{d\mathrm{\Phi }^{}}{dx}+\frac{g\pi ^2}{3}(\mathrm{\Phi }\mathrm{\Phi }^{})^3],$$
(1)
The first term is the kinetic energy density and the second term arises from interactions; here the dimensionless coupling $`g`$ may be positive or negative. Note that the interaction term in (1) is cubic in the density $`\rho =\mathrm{\Phi }\mathrm{\Phi }^{}`$, and not quadratic as in the standard theory. A rigorous derivation of this term for the repulsive case may be found in Ref. . It suffices to state here that in one dimension, where particles are not allowed to cross, such a term naturally arises when the interaction is short-range. This has been demonstrated for the hard-core interaction by Girardeau  in the thermodynamic limit, and by Lieb and Liniger  for a repulsive delta function potential. Even with a two-body interaction which varies inversely as the distance between particles, the leading term in the density functional has a cubic dependence to the density .

A few remarks may be made here in connection with the comment made by Bhaduri and Sen  recently. The Lieb and Liniger analysis  yields the result that in the limit of very low density of very strong repulsive interaction, the system becomes fermionic with $`g=1`$ in the functional given in Eq. (1). In the ground state, the density is constant and the kinetic term vanishes. Hence one obtains the free fermion result for the energy per unit length,

$$E=\frac{\mathrm{}^2}{2m}\frac{\pi ^2}{3}\rho ^3.$$
(2)
Obviously, in the free fermion limit no soliton solutions occur. In what follows, we keep both the kinetic term and the interaction term between the bosons with a general coupling $`g`$, and therefore are not in the free fermion limit.

The equation of motion may be obtained from the energy functional in Eq. (1), and is given by

$$\frac{\mathrm{}^2}{2m}[\frac{^2}{x^2}+g\pi ^2|\mathrm{\Phi }|^4]\mathrm{\Phi }(x,t)=i\mathrm{}\frac{\mathrm{\Phi }(x,t)}{t}.$$
(3)
Henceforth we shall denote the partial derivatives with respect to $`x`$ by a prime and with respect to $`t`$ by a dot. It is useful to multiply Eq. (3) by $`\mathrm{\Phi }^{}(x,t)`$,

$$\frac{\mathrm{}^2}{2m}\mathrm{\Phi }^{}\mathrm{\Phi }^{\prime \prime }+\mathrm{\Phi }^{}W(x)\mathrm{\Phi }=i\mathrm{}\mathrm{\Phi }^{}\dot{\mathrm{\Phi }},$$
(4)
and take the complex conjugate,

$$\frac{\mathrm{}^2}{2m}\mathrm{\Phi }\mathrm{\Phi }^{\prime \prime }+\mathrm{\Phi }W(x)\mathrm{\Phi }^{}=i\mathrm{}\mathrm{\Phi }\dot{\mathrm{\Phi }}^{}.$$
(5)
where $`W(x)=(\mathrm{}^2/2m)g\pi ^2|\mathrm{\Phi }|^4`$.

We now set, quite generally,

$$\mathrm{\Phi }(x,t)=\sqrt{\rho (x,t)}\mathrm{e}^{i\theta (x,t)},$$
(6)
where $`\rho `$ and $`\theta `$ are real. It is useful to define the total particle number

$$N=𝑑x\rho ,$$
(7)
and the momentum functional

$$P=\frac{i\mathrm{}}{2}𝑑x[\mathrm{\Phi }^{}\mathrm{\Phi }^{}\mathrm{\Phi }^{}\mathrm{\Phi }]=\mathrm{}𝑑x\rho \theta ^{}.$$
(8)
We then obtain the equations of motion as follows. The difference between Eqs. (4) and (5) reduces to the continuity equation for the number density $`\rho `$,

$$\dot{\rho }+\frac{\mathrm{}}{m}(\rho \theta ^{})^{}=\mathrm{\hspace{0.33em}0}.$$
(9)
Thus it is legitimate to regard the velocity of the “fluid” to be

$$v(x,t)=\frac{\mathrm{}}{m}\theta ^{}(x,t).$$
(10)
The energy density equation is obtained by adding Eqs. (4) and (5) and dividing by 2,

$$\tau (x,t)+W(x)\rho =\mathrm{}\rho \dot{\theta },$$
(11)
where the kinetic energy density is

$`\tau (x,t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{4m}}(\mathrm{\Phi }^{}\mathrm{\Phi }^{\prime \prime }+\mathrm{\Phi }^{\prime \prime }\mathrm{\Phi })`$ (12)
$`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}\left({\displaystyle \frac{1}{4}}{\displaystyle \frac{(\rho ^{})^2}{\rho }}{\displaystyle \frac{1}{2}}\rho ^{\prime \prime }+\rho (\theta ^{})^2\right).`$ (13)
Note that the Eqs. (9) and (11) are invariant under scaling and Galilean boosts . Under scaling by a factor $`\alpha `$,

$`\rho (x,t)`$ $``$ $`\alpha \rho (\alpha x,\alpha ^2t),`$ (14)
$`\theta (x,t)`$ $``$ $`\theta (\alpha x,\alpha ^2t).`$ (15)
Then the particle number $`N`$ remains invariant, while the energy and momentum transform as

$$E\alpha ^2E\mathrm{and}P\alpha P.$$
(16)
Under a Galilean boost by velocity $`v`$,

$`\rho (x,t)`$ $``$ $`\rho (xvt,t),`$ (17)
$`\theta (x,t)`$ $``$ $`\theta (xvt,t)+{\displaystyle \frac{mv}{\mathrm{}}}(x{\displaystyle \frac{1}{2}}vt).`$ (18)
Then the particle number $`N`$ remains invariant while the energy and momentum transform as

$`E`$ $``$ $`E+Pv+{\displaystyle \frac{1}{2}}mNv^2,`$ (19)
$`P`$ $``$ $`P+mNv.`$ (20)
We should mention a subtlety here. If we are interested in solutions involving a single localized object (such as the solitons discussed in the next section), then it is reasonable to demand that its energy and momentum should only change by a finite amount under a boost. From Eq. (20), we see that this is true only if the particle number $`N`$ is finite. Therefore we will perform boosts only when $`N`$ is finite.

## III Stationary Solutions: Some Analytical Results

We will now find some stationary solutions of Eqs. (9) and (11). For such solutions we make the following ansatz,

$$\mathrm{\Phi }(x,t)=\sqrt{\rho (x)}\mathrm{exp}(i\mu t/\mathrm{}),$$
(21)
namely, we set $`\theta (x,t)=\mu t/\mathrm{}`$. Note that this ansatz is consistent with the continuity equation since both $`\dot{\rho }`$ and $`\theta ^{}`$ are zero. The equation for $`\rho (x)`$ is then given by

$$\mu \rho =\frac{\mathrm{}^2}{2m}\left(\frac{1}{4}\frac{(\rho ^{})^2}{\rho }\frac{1}{2}\rho ^{\prime \prime }+g\pi ^2\rho ^3\right).$$
(22)
Then the following solutions exist.

### A Thomas-Fermi limit

Assuming the bulk density to be a constant, we set $`\rho (x)=\rho _0`$. This is an allowed solution with

$$\rho _0=\sqrt{2m\mu /g\pi ^2\mathrm{}^2},$$
(23)
provided $`g`$ is positive. This is the Thomas-Fermi result . Computing the energy, we reproduce the result given in Eq. (2) when $`g=1`$.

### B Repulsive interaction

We assume $`g`$ is positive and look for solutions when the density is not a constant. The Kolomeisky et al. solution given in their Eq. (12) may be written in the stationary frame with their velocity parameter $`\beta =0`$ . The solution is given by

$$\rho (x)=\frac{k}{\pi \sqrt{g}}\frac{\mathrm{cosh}(2kx)1}{\mathrm{cosh}(2kx)+2}$$
(24)
with

$$\mu =\frac{\mathrm{}^2k^2}{2m}.$$
(25)
Note that the density tends to $`\rho _0=1`$ at $`x\pm \mathrm{}`$; hence it is better to define the particle number of the soliton to be

$$N_{soliton}=𝑑x[\rho _0\rho (x)]=\frac{3}{2\pi \sqrt{g}}_{\mathrm{}}^{\mathrm{}}\frac{dt}{\mathrm{cosh}t+2}=\frac{\sqrt{3}}{\pi \sqrt{g}}\mathrm{ln}\left[\frac{\sqrt{3}+1}{\sqrt{3}1}\right].$$
(26)
In non-linear optics this number corresponds to the total intensity of light with a constant background subtracted out. Since $`N_{soliton}`$ is related to $`g`$, there exists exactly one solution for the system for a given value of $`g`$ (up to a scaling of $`k`$).

The above solution is called a dark soliton in optics  because the density goes to zero at $`x=0`$, and the soliton physically corresponds to a dip in the density distribution. Note that the above solution is valid when $`g`$ is positive, i.e., for repulsive interactions.

One can get the energy, $`E`$ of the soliton by integrating the energy density  given in Eq. (1). Using the soliton profile given in Eq. (24) and using Eq. (22) it is easy to see that the energy of the soliton is infinite. However if the energy of the soliton is defined as the difference

$$E_{soliton}=EE_0,$$
(27)
where $`E_0`$ is the energy of the solution with $`\rho (x)=\rho _0=k/\pi \sqrt{g}`$, then it turns out that the energy of the soliton is zero. It is interesting that even in the presence of a repulsive interaction, it does not cost any energy to create a soliton.

Until now we have checked the results given in Ref.  in the stationary frame. Next we consider some new solutions in the presence of an attractive interaction.

### C Attractive interaction

The equation of motion for $`\mathrm{\Phi }(x,t)`$ in the stationary frame for $`\rho `$ may also be solved by the following ansatz

$$\rho =\frac{2kN}{\pi }\mathrm{sech}(2kx)\mathrm{and}\theta =\frac{\mu t}{\mathrm{}}.$$
(28)
Note that the density is normalized so that

$$_{\mathrm{}}^{\mathrm{}}𝑑x\rho (x)=N.$$
(29)
The ansatz in (28) satisfies the equation of motion provided that

$$N^2=\frac{3}{4g}\mathrm{and}\mu =\frac{\mathrm{}^2k^2}{2m}.$$
(30)
Since $`N`$ is real we are forced to choose $`g`$ to be negative. We may therefore write

$$N^2=\frac{3}{4|g|}.$$
(31)
Following Ref.  we call this the bright soliton. We note in passing that a spatial soliton of the secant hyperbolic type is well-known in nonlinear optics, and has been experimentally realized long back , although the differential equation obeyed by the solitonic electrical field in optics is very different from ours.

Since $`N`$ and $`g`$ are related to each other by Eq. (31), there exists exactly one solution (up to scaling) for the system given a value for either $`g`$ or $`N`$. The soliton solution in this case is a lump, and it has its maximum at $`x=0`$. Note the peculiar fact that the larger the magnitude of $`g`$, the smaller is the particle number $`N`$ of the soliton. The reason for this is explained in the next subsection where we will show that a configuration whose particle number is too large is unstable to collapse.

Interestingly, we find that the energy of the soliton in this case is zero. Since $`N`$ is finite, we can use Eq. (18), to boost the solution given in Eq. (28). Thus

$`\rho `$ $`=`$ $`{\displaystyle \frac{2kN}{\pi }}\mathrm{sech}(2k(xvt)),`$ (32)
$`\theta `$ $`=`$ $`{\displaystyle \frac{\mu t}{\mathrm{}}}+{\displaystyle \frac{mv}{\mathrm{}}}(x{\displaystyle \frac{1}{2}}vt).`$ (33)
We then find that

$$E=\frac{1}{2}mNv^2\mathrm{and}P=mNv.$$
(34)
Hence it is natural to define the mass of this soliton to be $`M=mN`$.

We further note that there also exist stationary wave solutions apart from the above. Starting from the dark soliton solution given in Eq. (24), we can transform $`kik`$ and $`\rho \rho `$ to obtain the following stationary wave solution for negative values of $`g`$:

$$\rho (x)=\frac{k}{\pi \sqrt{g}}\frac{1\mathrm{cos}(2kx)}{2+\mathrm{cos}(2kx)},$$
(35)
with

$$\mu =\frac{\mathrm{}^2k^2}{2m}.$$
(36)
Finally, note that we can apply the scaling transformation given in Eq. (15) for all the solutions given above; that would change the normalization of the density in all cases. But we can apply the boost transformation only for the bright soliton since that is the only case in which the total particle number is finite.

### D Stability of the soliton solutions

We will now study the stability of the solitons, particularly the bright soliton given in Eq. (28), using both variational and small fluctuation analysis. Before doing that, let us make some preliminary comments.

The soliton solutions obtained above for both attractive and repulsive cases are scale invariant. Hence the type of instability which is associated with the soliton solutions for certain types of nonlinear Schrödinger equations (see Barashenkov et al. ) are not expected here. Another form of instability which can occur is the breaking of the soliton into some other fundamental constituents . To check the possibility of that, let us note that our energy functional has a global $`U(1)`$ symmetry, and hence a conserved charge associated with such a soliton solution. For the attractive case, the charge is given by the expression

$$Q=i\mathrm{}𝑑x(\mathrm{\Phi }^{}\dot{\mathrm{\Phi }}\dot{\mathrm{\Phi }}^{}\mathrm{\Phi })=2\mu N=\frac{\mathrm{}^2k^2}{2m}\sqrt{\frac{3}{|g|}}.$$
(37)
This is of course proportional to the integral of the density $`\rho (x)`$ in our case which is the same in any frame. In that case the soliton solution is classically stable against the aforesaid complete dissociation into more fundamental constituents (mfc), because its rest energy which is $`0`$ here must be lower than $`m_0|Q|`$ if we assume that the mass $`m_0`$ of the hypothetical mfc must be positive.

We will now show that in the attractive case $`g<0`$, the bright soliton can collapse to a point if the particle number of a configuration is bigger than a critical number $`N_c`$ which is a function of the parameter $`g`$. We can show this as follows. Consider a trial wave function of the form

$$\mathrm{\Phi }(x)=\left(\frac{N\alpha }{\pi \mathrm{cosh}(\alpha x)}\right)^{1/2},$$
(38)
where $`\alpha `$ is a free parameter (whose inverse determines the width of the configuration), and we set $`\theta =0`$. The particle number for this configuration is given by $`N`$. We then find that the energy functional in Eq. (1) is given by

$$E=\frac{\alpha ^2\mathrm{}^2N}{16m}[1\frac{4}{3}|g|N^2].$$
(39)
for this configuration. We observe that this expression does not have a minimum at any finite and nonzero value of $`\alpha `$. If

$$N^2>N_c^2\frac{3}{4|g|},$$
(40)
then the energy of the configuration given in Eq. (38) can be made arbitrarily large and negative by letting $`\alpha \mathrm{}`$; this corresponds to the density of the configuration collapsing to a $`\delta `$-function at the origin. We thus see that there is a critical value of $`N_c`$ which is proportional to $`1/\sqrt{|g|}`$ beyond which there is no lower bound to the energy. We observe that the particle number of the bright soliton in Eq. (31) is exactly equal to this critical value. Since this argument is only based on a trial wave function, it does not prove anything about the stability of the bright soliton against a collapse.

It is instructive to contrast the situation here with a more general form of the GP functional where the interaction term in Eq. (1) is taken to be of the form $`g\rho ^n`$ instead of $`g\rho ^3`$, where we assume $`n>1`$. If the exponent of $`\rho `$ is taken to be $`n<3`$ (instead of $`n=3`$ as we have considered so far), and $`g`$ is negative, we find from the above argument that there is no instability to a collapse, no matter what the values of $`g`$ and $`N`$ are. Namely, for any value of $`g<0`$ and $`N`$, the energy of the trial configuration given in (38) has a minimum at some finite value of the inverse width $`\alpha `$. On the other hand, if we take the exponent $`n>3`$ and $`g<0`$, the above argument indicates an instability to collapse for any value of $`N`$, no matter how small. Thus, an attractive interaction of the form $`\rho ^3`$ is rather special; this is the only value of the exponent $`n`$ for which there is a finite critical number $`N_c`$, such that the trial configuration is stable against a collapse if $`N<N_c`$ and is unstable if $`N>N_c`$. This special property of $`n=3`$ can be traced to the fact that the coefficient of the interaction $`g`$ is dimensionless only for this value of $`n`$.

Let us now return to the bright soliton given in Eq. (28), and directly study its stability by studying small fluctuations around it. Our main conclusion will be that there is at least one unstable mode. To do the analysis, let us set $`k=1/2`$ for simplicity in (28). Then the soliton takes the form

$`\rho _c`$ $`=`$ $`{\displaystyle \frac{N}{\pi \mathrm{cosh}x}},`$ (41)
$`\theta _c`$ $`=`$ $`{\displaystyle \frac{\mathrm{}t}{8m}},`$ (42)
with $`N^2g=3/4`$. We now assume small fluctuations of the form

$`\rho `$ $`=`$ $`{\displaystyle \frac{N}{\pi \mathrm{cosh}x}}+\delta \rho (x,t),`$ (43)
$`\theta `$ $`=`$ $`{\displaystyle \frac{\mathrm{}t}{8m}}+\delta \theta (x,t).`$ (44)
The equations of motion then imply that

$`{\displaystyle \frac{\delta \rho }{t}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}N}{\pi m}}({\displaystyle \frac{1}{\mathrm{cosh}x}}{\displaystyle \frac{^2}{x^2}}{\displaystyle \frac{\mathrm{sinh}x}{\mathrm{cosh}^2x}}{\displaystyle \frac{}{x}})\delta \theta ,`$ (45)
$`{\displaystyle \frac{\delta \theta }{t}}`$ $`=`$ $`{\displaystyle \frac{\pi \mathrm{}}{4Nm}}\mathrm{cosh}x({\displaystyle \frac{^2}{x^2}}+{\displaystyle \frac{\mathrm{sinh}x}{\mathrm{cosh}x}}{\displaystyle \frac{}{x}}+{\displaystyle \frac{\mathrm{sinh}^2x}{2\mathrm{cosh}^2x}}+{\displaystyle \frac{9}{2\mathrm{cosh}^2x}}{\displaystyle \frac{1}{2}})\delta \rho .`$ (46)
To solve these equations, we eliminate $`\delta \theta `$ and define a function $`\psi (x)`$ as

$$\delta \rho (x,t)=(\mathrm{cosh}x)^{1/2}\psi (x)e^{i\omega t}.$$
(47)
We then obtain the eigenvalue equation

$$A\psi (x)=\omega ^2\psi (x),$$

where

$$A=\frac{\mathrm{}^2}{4m^2}(\frac{d^2}{dx^2}\frac{1}{4}+\frac{3}{4\mathrm{cosh}^2x})(\frac{d^2}{dx^2}\frac{1}{4}+\frac{15}{4\mathrm{cosh}^2x})\psi .$$
(48)
The operator $`A`$ is hermitian with the weight $`w(x)=1`$; namely, for any two normalizable functions $`\psi _1(x)`$ and $`\psi _2(x)`$,

$$_{\mathrm{}}^{\mathrm{}}𝑑xw(x)\psi _1^{}(x)A\psi _2(x)=_{\mathrm{}}^{\mathrm{}}𝑑xw(x)(A\psi _1(x))^{}\psi _2(x),$$
(49)
with $`w(x)=1`$. (It is for this reason that we introduced the factor of $`(\mathrm{cosh}x)^{1/2}`$ in the definition of $`\psi `$ in Eq. (47)). Hence the eigenvalues of $`A`$, $`\omega ^2`$, must be real. However, $`\omega ^2`$ could be negative; in that case, there would be a solution in which $`\omega `$ is purely imaginary and $`i\omega >0`$. That would be an unstable mode with the function $`\delta \rho (x,t)`$ blowing up exponentially as $`t`$ goes to $`\mathrm{}`$. We therefore need to know if $`A`$ has any negative eigenvalues.

We already know two eigenfunctions of $`A`$ with $`\omega =0`$ corresponding to the translation and scaling symmetries. From infinitesimal versions of these two symmetries, we find that the eigenfunctions of the operator $`A`$ with zero eigenvalues are given by

$$\psi =\frac{\mathrm{sinh}x}{(\mathrm{cosh}x)^{3/2}},$$
(50)
and

$$\psi =\frac{1}{(\mathrm{cosh}x)^{1/2}}\frac{x\mathrm{sinh}x}{(\mathrm{cosh}x)^{3/2}}.$$
(51)
Note that these two functions are normalizable (i.e., they would be called bound states in the language of quantum mechanics), and that they have one and two nodes respectively. Thus the operator $`A`$ has quite different properties from the usual Hamiltonian operator $`H=d^2/dx^2+V(x)`$ which appears in a one-dimensional Schrödinger equation; a second derivative operator like $`H`$ always has non-degenerate eigenvalues for bound states, and the corresponding eigenvalues strictly increase with the number of nodes.

We have not been able to find an analytical expression for the eigenfunction of $`A`$ corresponding to an eigenvalue $`\omega ^2<0`$. However it is possible to show its existence by a variational argument. We use the result from quantum mechanics that for a hermitian operator $`A`$, if $`\psi _v(x)`$ is a normalizable variational function, and

$$E_v=\frac{<\psi _v|A|\psi _v>}{<\psi _v|\psi _v>}$$
(52)
is the corresponding variational energy, then there must be an eigenvalue of $`A`$ which is less than or equal to $`E_v`$. We now try a function of the form

$$\psi _\alpha (x)=\frac{1}{(\mathrm{cosh}x)^{1/2}}+\alpha \frac{x\mathrm{sinh}x}{(\mathrm{cosh}x)^{3/2}},$$
(53)
where $`\alpha `$ is a variational parameter which we assume to be real for convenience. Then the variational energy given in (52) takes the form

$$E_\alpha =\frac{3\mathrm{}^2}{10m^2}\frac{\alpha (1+\alpha )}{1+2\alpha +\alpha ^2(1+\pi ^2/8)}.$$
(54)
We find numerically that this has a minimum value of $`E=0.1350\mathrm{}^2/m^2`$ at $`\alpha =0.4738`$. Since this minimum value is negative, we know that there must be an eigenfunction of the operator $`A`$ for which $`\omega ^2<0`$. This shows that the bright soliton has at least one small fluctuation mode which blows up as $`t\mathrm{}`$. Since we do not know the explicit form of this eigenfunction, the physical interpretation of this mode remains unclear.

Although we have found an instability classically, it is possible that the instability may disappear if the theory is quantized, i.e., if we impose the equal-time commutation relation $`[\rho (x),\theta (y)]=i\delta (xy)`$ . There are examples known in the literature (for instance, skyrmions in three dimensions) where the quantization of some collective coordinates eliminates a classical instability . However, we will not examine the question of the quantum stability of the bright soliton in this paper where our entire analysis has been classical.

## IV General treatment of stationary solutions

We will now see how the general form of the stationary solution of Eq. (3) may be found. We do this by mapping the one-dimensional problem considered in previous sections to that of a pseudo-particle problem in two dimensions. To do this let us write

$$\mathrm{\Phi }(x,t)=e^{i\mu t/\mathrm{}}[\xi _1(x)+i\xi _2(x)].$$
(55)
Thus the density given by $`\rho =\xi _1^2+\xi _2^2`$ is independent of time. Eq. (3) then takes the form

$`{\displaystyle \frac{d^2\xi _1}{dx^2}}`$ $`=`$ $`{\displaystyle \frac{2m\mu }{\mathrm{}^2}}\xi _1+g\pi ^2\rho ^2\xi _1,`$ (56)
$`{\displaystyle \frac{d^2\xi _2}{dx^2}}`$ $`=`$ $`{\displaystyle \frac{2m\mu }{\mathrm{}^2}}\xi _2+g\pi ^2\rho ^2\xi _2.`$ (57)
These are precisely the equations of motion of a classical particle of unit mass moving in two dimensions with the spatial coordinates $`(\xi _1,\xi _2)`$ and a time coordinate $`x`$. Note that in the one-dimensional problem $`x`$ refers to position coordinate. The total energy and pseudo-angular momentum of the particle are given by

$``$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{d\xi _1}{dx}}\right)^2+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{d\xi _2}{dx}}\right)^2+{\displaystyle \frac{m\mu }{\mathrm{}^2}}\rho {\displaystyle \frac{g\pi ^2}{6}}\rho ^3,`$ (58)
$`L`$ $`=`$ $`\xi _1{\displaystyle \frac{\xi _2}{x}}\xi _2{\displaystyle \frac{\xi _1}{x}}.`$ (59)
The pseudo-angular momentum $`L`$ is in fact the same, up to a factor, as the momentum density given in Eq. (8). These are conserved since the particle is moving in a central potential

$$V(\xi _1,\xi _2)=\frac{m\mu }{\mathrm{}^2}\rho \frac{g\pi ^2}{6}\rho ^3.$$
(60)
Since there are two conserved quantities and the particle has two degrees of freedom, the motion can be found exactly if the initial conditions are given. For instance, $`\rho `$ is given in terms of $`x`$ by the equation

$$\left(\frac{d\rho }{dx}\right)^2=4[2\rho L^2\frac{2m\mu }{\mathrm{}^2}\rho ^2+\frac{g\pi ^2}{3}\rho ^4].$$
(61)
We are interested in solutions of Eqs. (57) and (61) in which the particle density $`\rho `$ remains finite as $`x\pm \mathrm{}`$. In that case, the solution for $`\rho `$ will generally be a periodic function of $`x`$; this corresponds to a stationary wave. From the form of the potential energy in (59), we see that this will happen in the following situations.

(i) If $`g>0`$, the potential has a negative sixth-degree term. Hence we must have $`\mu >0`$ and the particle energy $``$ must be equal to or less than the potential barrier (this can be derived from Eq. (60)) in order to have $`\rho `$ finite for all times $`x`$. Then $`\rho `$ will always remain bounded by the location of the potential barrier given by $`\rho _0=\sqrt{2m\mu /g\pi ^2\mathrm{}^2}`$ as in Eq. (23).

(ii) If $`g<0`$, the potential has a positive sixth-degree term, so $`\rho `$ will be finite regardless of whether $`\mu `$ is positive or negative. However the possible motions are different for $`\mu >0`$ and $`\mu <0`$ as we will see.

We will now discuss how the various solutions discussed in the previous section can be obtained as special cases from this general discussion. First of all, let us set the pseudo-angular momentum $`L=0`$ and therefore $`\xi _2=0`$. Let us now consider the cases with $`g`$ positive and negative separately.

(i) In the repulsive case $`g>0`$, we have already seen that we need $`\mu >0`$ and the energy must lie below the potential barrier so that the particle does not escape to infinitely large values of $`\rho `$. If the energy is exactly equal to the potential barrier, then we have two possible solutions. Either the particle stays at the top of the potential barrier for all times $`x`$ (this corresponds to the Thomas-Fermi result in (23)), or we have a solution in which the particle begins at $`\xi _1=\sqrt{\rho _0}`$ at time $`x=\mathrm{}`$ and ends at $`\xi _1=\sqrt{\rho _0}`$ at time $`x=\mathrm{}`$. The latter solution corresponds to the dark soliton of Kolomeisky et al. On the other hand, if the energy is less than the potential barrier, then the particle will perform a periodic motion which repeatedly passes through $`\xi _1=0`$; this corresponds to a stationary wave with the minimum density being zero at an infinite number of values of $`x`$. This is a new stationary wave solution.

(ii) In the attractive case $`g<0`$, there is only one kind of motion possible if $`\mu >0`$, i.e. the particle will perform a periodic motion passing repeatedly through $`\xi _1=0`$, thus giving a new stationary wave. If $`\mu <0`$, we have a double well potential (see Eq. (60), and there are various kinds of motion possible depending on the energy $``$. If $`>0`$, then the particle again performs a periodic motion passing through $`\xi _1=0`$, thereby giving the stationary wave described in Sec. II C. If $`=0`$, then we have a bounce solution in which the particle begins at $`\xi _1=0`$ at $`x=\mathrm{}`$, goes to some maximum (or minimum) value of $`\xi _1`$ at some finite time and then returns to $`\xi _1=0`$ at $`x=\mathrm{}`$. This corresponds to the bright soliton described in Sec. II C; note that such a solution necessarily has $`N^2g=3/4`$ no matter what the value of $`\mu `$ is. (Thus there is no bright soliton solution possible if $`N^2g3/4`$). If $``$ is negative but is greater than the bottom of the potential, the particle performs a periodic motion in which $`\xi _1`$ always remains positive (or negative) and never reaches zero. This is a new solution corresponding to a stationary wave whose minimum density is nonzero. Finally, if $``$ is equal to the bottom of the potential located at $`\rho _0=\sqrt{2m\mu /g\pi ^2\mathrm{}^2}`$, then we have a solution with uniform density $`\rho (x)=\rho _0`$.

## V Comments

We have shown that the modified one-dimensional GP mean-field theory has many solutions apart from the one outlined by Kolomeisky et al. . The soliton solution in the attractive case is particularly simple since the full time dependent solutions in an arbitrary moving frame may be obtained simply by boosting the solutions in the stationary frame. The question of the quantum stability of this soliton remains an open problem.

S. G. thanks R. Rajaraman for discussions and comments. R.K.B. acknowledges financial support from NSERC, Canada, for this research, and the hospitality of I.M.Sc., Chennai, India.


text:
# On the Mass-Donating Star in GRS 1915+105

## 1 Introduction

As the first Galactic source of superluminal jets, GRS 1915+105 has undergone intensive study during recent years. Despite this level of investigation, one of the system’s key parameters – the approximate mass of the mass-donating star in the binary – remains a subject of controversy. No binary modulation in the system, such as eclipses, have been observed at any wavelength, preventing insights based on orbital properties and evolution. Mirabel et al., (1997) observed near-infrared emission lines of He I and Br$`\gamma `$, which they interpreted as evidence for a massive Oe- or Be-type companion star. However, Castro-Tirado et al. (1996) and Eikenberry et al. (1998) reported the observation of a weak, time-variable He II emission line (later confirmed by Marti et al., 2000). Since normal Oe and Be stars do not typically doubly-ionize He, they suggested that the emission lines arise in the accretion disk surrounding the compact object, leaving the companion star type unconstrained. Furthermore, Eikenberry et al. (1998) observed the emission lines changing in strength by factors of several on timescales of $`510`$ minutes during jet-producing activity. They interpreted the constancy of the line profiles during these rapid changes as further evidence of the accretion disk origin of the lines.

More recently, Marti et al. (2000) have presented VLT spectra of GRS 1915+105, in which they confirm the weak He II emission line and identify another emission line due to NaI ($`2.21\mu `$m). They interpret this feature as evidence that the mass-donating star in GRS 1915+105 is a super-massive star, such as a luminous blue variable (LBV) star. However, they note that this interpretation does not seem able to explain the line variability observed by Eikenberry et al. (1998).

In this paper, we will also address the nature of the mass-donating star. Based on the known inclination angle of the jet axis, its stability over timescales of years, and the lack of X-ray eclipses in GRS 1915+105, we place constraints on the ratio of the mass-donating star’s radius $`R_s`$ to the semi-major axis of the binary system $`a`$. We then show, based on estimates of the black hole mass in GRS 1915+105, its X-ray luminosity, and the observed X-ray absorption column towards GRS 1915+105, that the accretion in the system cannot be provided by the wind of a massive star. We then show that for a $`30M_{\mathrm{}}`$ black hole undergoing accretion via Roche lobe overflow, the mass-donating star is constrained to have $`M_s<19.4M_{\mathrm{}}`$. Next, we use limits on X-ray reprocessing on the stellar surface to place further constraints on the mass-donating star. Finally, we discuss the implications of these results for GRS 1915+105 and present our conclusions.

## 2 Geometry of GRS 1915+105

In the paper announcing the discovery of superluminal motion from GRS 1915+105, Mirabel & Rodriguez (1994) were able to constrain the inclination angle of the jet axis to be $`\theta =70\mathrm{°}\pm 2\mathrm{°}`$. Furthermore, in the observed jet outbursts since then, this angle has remained stable to within a few degrees on timescales of years (e.g. Fender et al., 1999). Thus, the orbital motion of the compact object does not significantly alter the direction of the jet axis. This constancy can be provided by one of two possibilites: either the mass ratio of the system is small ($`q=M_s/M_{BH}1`$), or else the orbital axis is aligned with the jet axis. The first case in and of itself strongly constrains the nature of the companion star. In the second case, further analysis is required to see if the nature of the companion star can be constrained.

The first component in this analysis is the lack of X-ray eclipses in GRS 1915+105. This object has been monitored by the Rossi X-ray Timing Explorer (RXTE) All-Sky Monitor since early 1996, and with a huge number of pointed observations using the RXTE Proportional Counter Array (e.g. Muno et al., 2000). It has also been extensively studied using ASCA (e.g. Kotani et al., 2000), BeppoSAX (e.g. Feroci et al., 2000), Granat-Sigma (Castro-Tirado, et al., 1994), and the Compton Gamma-Ray Observatory (e.g. Harmon, et al., 1997). Despite this wealth of coverage, no X-ray eclipse has ever been observed from GRS 1915+105. For an orbit aligned with the jet axis at $`\theta =70`$ degrees, this constrains the ratio of the star radius to the orbital separation to be

$$\frac{R_s}{a}\mathrm{cos}(70\mathrm{°})=0.342$$

## 3 Accretion via Stellar Wind

### 3.1 Limits on $`\dot{M}`$ and the Stellar Wind

If the mass-donating star in GRS 1915+105 is in fact an early-type star, one possibility for providing the accreted mass is through a stellar wind. The accretion rate must be able to account for the high X-ray luminosity observed from GRS 1915+105, with $`L_x3\times 10^{39}\mathrm{ergs}\mathrm{s}^1`$ (e.g. Belloni et al., 1997). If we assume a canonical efficiency of 15% for converting accreted mass into radiated energy (e.g. Frank et al., 1992), we arrive at an accretion rate of $`\dot{M}3\times 10^7M_{\mathrm{}}\mathrm{yr}^1`$.

Wind accretion is known to be rather inefficent compared to Roche lobe overflow – only the wind material directed towards the compact object is accreted. More quantitatively, the accretion rate from the wind $`\dot{M}_{wind}`$ is related to the stellar mass loss rate $`\dot{M}_s`$ according to

$$\dot{M}_{wind}\frac{1}{4}\dot{M}_sq^2(R_s/a)^2$$

The mass of the black hole in GRS 1915+105 is not well-determined. However, Zhang, et al. (1997) and Nowak, et al. (1997) suggest $`M_{BH}30M_{\mathrm{}}`$ based on the 67 Hz quasi-periodic oscillation observed from the system (Morgan et al., 1997). Furthermore, by requiring the peak X-ray luminosity to be less than the Eddington luminosity, we find $`M_{BH}>23M_{\mathrm{}}`$ – consistent with the result of Zhang, et al. (1997). Therefore, we adopt $`M_{BH}=30M_{\mathrm{}}`$ here. For a super-massive star as suggested by Marti et al., we have $`M_s50M_{\mathrm{}}`$ (Humphreys & Davidson, 1994), giving $`q1.6`$. For an accretion rate $`\dot{M}3\times 10^7M_{\mathrm{}}\mathrm{yr}^1`$, we then arrive at $`\dot{M}_s2.7\times 10^5M_{\mathrm{}}\mathrm{yr}^1`$. This mass-loss rate is very high – even the LBV stars have an upper end to their mass loss between outbursts of $`10^5M_{\mathrm{}}\mathrm{yr}^1`$ (e.g. Davidson & Humphreys, 1997). Since GRS 1915+105 has shown repeated patterns of X-ray emission over timescales of several years, implying steady accretion over this interval as compared to the variable aperiodic outflows typical of LBV stars in outburst (e.g. Humphreys & Davidson, 1994), it seems likely that any LBV-type mass-donating star is not in outburst. Thus, the wind mass-loss rate required to produce the observed X-ray emission is marginally consistent with the mass-loss rates known to occur in such stars.

The required wind mass-loss rate decreases with decreasing $`q`$ until $`q\frac{R_s}{a}=0.34`$, at which point essentially all of the stellar wind not initially flowing away from the compact object is captured. For $`q=0.34`$, we have $`\dot{M}_s1.2\times 10^6M_{\mathrm{}}\mathrm{yr}^1`$ for a $`10M_{\mathrm{}}`$ star. This mass-loss rate is again consistent with the rates inferred for the O/B companions in wind-fed high-mass X-ray binaries (e.g. Frank et al., 1992), though this stellar mass range is at the lower limit for such stars.

### 3.2 Limits from X-ray Absorption

The tremendous stellar wind required to provide the X-ray luminosity of GRS 1915+105 should also have an impact on the observed absorbing column density towards the X-ray source. Given a minimum stellar mass-loss rate from above and the wind velocity $`v_{wind}`$, we can place a lower limit on the wind density as a function of distance from the star’s center $`r`$ according to

$$\rho (r)=\frac{\dot{M}_s}{4\pi r^2v_{wind}}$$

We take $`v_{wind}=10^3\mathrm{km}\mathrm{s}^1`$ here – somewhat larger than typical for these systems, making the above density a firm lower limit.

For this minimum $`\dot{M}_s`$, the ratio $`R_s/a`$ is fixed at 0.342. Furthermore, for supermassive LBV stars, an upper limit on their size is $`R_s1`$ AU, giving an upper limit on the orbital separation of $`a_{max}=2.9`$ AU. For a less-massive O/B-type star, we have $`R_s0.1`$ AU, giving $`a_{max}=0.29`$ AU. Given these system parameters and assuming a uniform spherical wind, we can calculate the line-of-sight column density through the wind for an observer at a $`20\mathrm{°}`$ angle above the orbital plane. It is important to note that the assumption of uniform spherical outflow is not generally critical for the results of this calculation. Rather, we only require that the line of sight sample the typical density regimes of any more-complicated outflow – a requirement that will be satisfied for all but the most contrived geometries in this case. Based on this, we estimate the lower limit to the column density to be $`n_H1.5\times 10^{23}\mathrm{cm}^2`$ for an LBV-type star and $`n_H7\times 10^{22}\mathrm{cm}^2`$ for an O/B-type star when the mass-donating star lies between the observer and the X-ray source.

For comparison, the measured total column density towards GRS 1915+105 is $`n_H6\times 10^{22}\mathrm{cm}^2`$ (e.g. Belloni et al., 1997)111This measurement assumes solar abundances. For enriched winds (commonly seen in these stars), the apparent absorption would be even higher, increasing the contrast between observations and the wind-accretion model., and the vast majority of this absorption is attributable to the intervening ISM (Chaty, et al., 1996). Furthermore, as the orbital orientation changes, any wind-induced absorption column would vary systematically, dropping by a factor of $`>5`$ at opposition – see Figure 1. However, spectra collected over long time intervals show no such variability in absorption for GRS 1915+105(e.g. Muno et al., 1999). So, even the lowest estimates for stellar winds capable of providing the accreting matter produce X-ray absorption greatly in excess of that observed. Therefore, we conclude that a stellar wind cannot provide the accretion observed in GRS 1915+105.

## 4 Accretion via Roche-lobe Overflow

### 4.1 Maximum stellar mass

Roche-lobe overflow is generally a much more efficient mode of mass-transfer than stellar winds, with essentially all of the mass-loss transferring onto the compact object. Thus, we have $`\dot{M}_s3\times 10^7M_{\mathrm{}}\mathrm{yr}^1`$. Moreover, as Eggleton (1983) has shown, Roche lobe overflow requires that the mass ratio $`q`$ and the ratio of the stellar radius to the binary separation follow

$$\frac{R_s}{a}=\frac{0.49q^{2/3}}{0.6q^{2/3}+\mathrm{ln}(1+q^{1/3})}$$

For our upper limit of $`R_s/a0.342`$ above, we arrive at an upper limit of $`q0.647`$, corresponding to a limit on the stellar mass of $`M_s19.4M_{\mathrm{}}`$ (assuming $`M_{BH}=30M_{\mathrm{}}`$). This upper limit is more than a factor of 2 below the lower limit for supermassive stars of the type suggested by Marti et al. (2000) as the companion for GRS 1915+105.

### 4.2 Limits from X-ray reprocessing

We can also place limits on Roche-lobe-filling companion stars to GRS 1915+105 due to the apparent lack of significant infrared (IR) flux from reprocessed X-rays in the system (Eikenberry et al., 2000). They place limits on the IR flux density due to reprocessing of $`\mathrm{\Delta }F_\nu <3\times 10^{26}\mathrm{ergs}\mathrm{s}^1\mathrm{cm}^2\mathrm{Hz}^1`$ at $`\nu =1.4\times 10^{14}`$ Hz ($`2.2\mu `$m), in the presence of $`\mathrm{\Delta }L_x10^{39}\mathrm{ergs}\mathrm{s}^1`$. For thermal reprocessing on a blackbody of temperature $`T`$, we have a reprocessing efficiency $`ϵ`$ (ratio of the bolometric reprocessed luminosity to the emitted X-ray luminosity) given by

$$ϵ=\frac{c^2\sigma T^3}{2\pi \nu ^2k}\frac{\mathrm{\Delta }L_\nu }{\mathrm{\Delta }L_x}$$

assuming the change in temperature is small compared to the temperature ($`\mathrm{\Delta }TT`$), changes in the surface area of the blackbody are negligible, and $`\mathrm{\Delta }L_\nu =4\pi d^2\mathrm{\Delta }F_\nu `$ where we take $`d=12.5`$ kpc (Mirabel & Rodriguez, 1994). From simple geometry, we can see that $`ϵ`$ should also have the form

$$ϵ=(1\eta )\frac{\pi R_s^2}{4\pi a^2}$$

where $`\eta `$ is the X-ray albedo of the stellar surface (assumed to be $`\eta 0.5`$ here – e.g. Anderson, 1981).

Given this relationship between $`ϵ`$ and $`\frac{R_s}{a}`$ at a given temperature $`T`$ and the relationship between $`q`$ and $`\frac{R_s}{a}`$ for Roche-lobe overflow, we can then calculate an upper limit on $`q`$ as a function of stellar surface temperature. Figure 2 shows a plot of the parameter space (mass, temperature) allowed/excluded for Roche-lobe-filling stars222Note that for stars with $`M_s3M_{\mathrm{}}`$ ($`q0.1`$), the constraint on the alignment of the orbital axis with the jet axis may be relaxed, so that the stellar temperature is not constrained by this method. (assuming $`M_{BH}=30M_{\mathrm{}}`$).

One important consideration here is that we have ignored the possible effects of the accretion disk casting an X-ray shadow on the star, thus reducing the amount of reprocessing on the stellar surface. However, this intercepted radiation does not simply disappear – rather it is reprocessed on the disk itself. While the efficiency of this reprocessing depends on the temperature and albedo of the disk region intercepting the radiation, the star subtends only an azimuthal range of $`\frac{2R_s}{a}`$ radians as seen from the center of the accretion disk, while the disk extends a full $`2\pi `$ radians. Thus, if the disk intercepts any fraction of the X-ray light that would otherwise reach the star, it will in fact intercept at least $`\frac{\pi a}{R_s}`$ times that amount of light over its total azimuthal range. Given our limit of $`\frac{R_s}{a}<0.342`$, this means that as long as the disk reprocessing efficiency is at least 10% of the stellar reprocessing efficiency per unit area, any disk shadowing would actually increase the total amount of reprocessed light over the above estimate, and thus tighten the constraint given the upper limits on reprocessed radiation. Therefore, we consider the constraints in Figure 2 to be fairly robust against the effects of disk shadowing.

## 5 Discussion

### 5.1 The companion star mass

As shown above, any reasonable parameters for accretion via stellar wind produce large X-ray absorption which is inconsistent with observations of GRS 1915+105. Furthermore, for accretion via Roche-lobe overflow, we can rule out any star with $`M_s>25M_{\mathrm{}}`$ as well as all but the hottest stars with $`M_s3M_{\mathrm{}}`$. Thus, a super-massive companion, as suggested by Marti et al. (2000), does not seem possible for GRS 1915+105.

The likeliest “normal” candidates for the mass-donating star are either a main-sequence or giant (luminosity class III) early-B star, or a low-mass ($`M_s3M_{\mathrm{}}`$) star. For the expected stellar radii and masses of massive stars and assuming accretion via Roche lobe overflow, we would then expect binary periods for massive donor stars in GRS 1915+105 to range over $`1015`$ days. For low-mass companions, the range of expected binary periods is very poorly constrained.

Given the $`\dot{M}3\times 10^7M_{\mathrm{}}\mathrm{yr}^1`$ accretion rate required to produce the X-ray luminosity, one might expect a low-mass star to imply a short lifetime for the “active” phase of GRS 1915+105 – $`\tau M/\dot{M}10^610^7`$ years. However, while GRS 1915+105 has been active for most of the last 8 years, we know it was quiescent for several decades prior to this. Therefore, the duty cycle for active mass transfer in this system may be $`1`$, so that the lifetime can be considerably larger than the above estimate.

### 5.2 Whence the IR emission lines?

The primary motivation for the initial suggestion of a super-massive companion by Marti et al. (2000) was the IR emission line spectrum. The stars allowed under the above analysis could not produce sufficient ionizing radiation to explain the observed IR spectrum from GRS 1915+105. However, it is important to note that the IR lines observed from LBV and other such stars are not produced directly in the stellar photosphere/chromosphere. Rather, they are the result of the dense stellar wind being irradiated by UV emission from the luminous hot photosphere.

In the case of GRS 1915+105, we know that the system contains a hot, highly-luminous blackbody – the inner accretion disk, which we observe directly in X-rays, should also produce tremendous ionizing UV luminosity. Furthermore, the outer portions of such a disk are widely considered to be the sources of strong winds, and we are most interested in this particular system due to its strong outflows. Recent X-ray spectral analysis by Kotani et al. (2000)also indicate the presence of such a disk wind. Therefore, we hypothesize that the emission spectrum observed by Marti et al. (2000) (as well as others previous to them) derives not from a super-massive companion star, but from the accretion disk itself.

This latter explanation also fits well with the observations of Eikenberry et al. (1998). The observed rapid variability of the IR emission lines is correlated with flaring behavior in GRS 1915+105 in a manner suggesting radiative pumping. Since the line profiles are constant, the line-emitting region must be uniformly irradiated by the flares which appear to come from the inner accretion disk. This is not possible for a line emission from the envelope of a massive companion – the star itself will shield a considerable fraction of its envelope from the flare. However, if the lines arise from the disk, the line-emitting region is axisymmetric about the flaring region, guaranteeing uniform irradiation and maintaining the constancy of the emission line profiles.

### 5.3 The black hole mass

A key assumption for much of the above analysis is the mass of the black hole in GRS 1915+105. As noted above, both arguments based on the Eddington luminosity and (more speculative) interpretations of quasi-periodic oscillations imply $`M_{BH}30M_{\mathrm{}}`$. This estimate would make the black hole in GRS 1915+105 the most massive of the “stellar-mass” black hole candidates in the Galaxy – roughly 3 times more massive than the black hole in Cyg X-1. Based on the above arguments, $`M_{BH}=30M_{\mathrm{}}`$ eliminates the possibility of a super-massive companion star ($`M_s50M_{\mathrm{}}`$) as suggested by Marti et al. (2000). Furthermore, in order to allow a companion consistent with even the lowest portion of the LBV mass range would require $`M_{BH}75M_{\mathrm{}}`$.

## 6 Conclusions

We have presented an analysis of the possible companion star mass ranges for GRS 1915+105, and come to the following conclusions:

* Based on the $`70\mathrm{°}`$ inclination angle and stability of the jet axis, we find that either the mass ratio of the system is $`q1`$, or the orbital axis is aligned with the jet axis.
* Based on the lack of X-ray eclipses, for massive companions ($`q0.1`$), the radius of the star must be $`R_s<0.342a`$, where $`a`$ is the semi-major axis of the orbit.
* Given these constraints and the X-ray luminosity of GRS 1915+105, all reasonable parameters assuming accretion via a stellar wind produce an X-ray aborption in the wind which is incompatible with observations. This indicates that accretion must occur via Roche-lobe overflow in GRS 1915+105.
* For Roche-lobe overflow, the constraint on $`R_s/a`$ requires the companion star to have $`M_s<19.4M_{\mathrm{}}`$ (assuming $`M_{BH}=30M_{\mathrm{}}`$).
* Based on limits for IR flux from reprocessed X-ray in GRS 1915+105 (Eikenberry et al., 2000), any companion star with $`3M_{\mathrm{}}<M_s<25M_{\mathrm{}}`$ must have a surface temperature $`T12\times 10^4`$ K. The “normal” main sequence and giant stars which fall into this range are early-B stars.
* The IR emission lines do not arise in the envelope of a massive companion, but instead are likely to be produced by the outer accretion disk and/or a disk wind being irradiated by the hot inner accretion disk.

SE is supported in part by an NSF Faculty Early Career Development (CAREER) award (NSF-9983830) and NASA grant NAG5-9178. This work was performed while RMB held a National Research Council Research Associateship Award at NRL.