text:
# Ruling out chaos in compact binary systems

## Abstract

We investigate the orbits of compact binary systems during the final inspiral period before coalescence by integrating numerically the second-order post-Newtonian equations of motion. We include spin-orbit and spin-spin coupling terms, which, according to a recent study by Levin \[J. Levin, Phys. Rev. Lett. 84, 3515 (2000)\], may cause the orbits to become chaotic. To examine this claim, we study the divergence of initially nearby phase-space trajectories and attempt to measure the Lyapunov exponent $`\gamma `$. Even for systems with maximally spinning objects and large spin-orbit misalignment angles, we find no chaotic behavior. For all the systems we consider, we can place a strict lower limit on the divergence time $`t_\mathrm{L}1/\gamma `$ that is many times greater than the typical inspiral time, suggesting that chaos should not adversely affect the detection of inspiral events by upcoming gravitational-wave detectors.

For the current generation of ground-based gravitational-wave (GW) detectors, such as LIGO, one of the most promising sources is the final inspiral of two compact stars (black holes or neutron stars in a relativistic binary orbit) thorne ; cutler . To detect such an event successfully, one must be able to match theoretical GW templates to experimental data containing a great deal of instrumental noise. This “matched filtering” technique has the potential of greatly increasing the effective signal-to-noise ratio for the detector thorne ; owen . Even to specify quantitative upper limits on inspiral events (a major objective of early LIGO runs), it is critical to have strong confidence in the physical accuracy of the templates allen .

Recent work has suggested that the orbits of two rapidly spinning compact objects may be chaotic levinprl ; levingrqc . The presence of chaos on the inspiral timescale (typically $`100`$s as the frequency sweeps up from $`10`$ to $`10^3`$Hz) could significantly reduce the probability of detection, even with GW templates that include spin effects kidder95 . In the matched filtering technique, if the signal and template are off by as little as half a wavelength (over $`10^3`$ cycles), the event could be missed because of destructive interference owen . In order to quantify this threat, we attempt to measure Lyapunov divergence times for a broad sample of initial conditions.

As in the Newtonian two-body problem, a relativistic binary system can be expressed in terms of a reduced mass $`\mu =m_1m_2/(m_1+m_2)`$ orbiting a fixed mass $`m=m_1+m_2`$ with a separation $`\stackrel{}{𝐫}=\stackrel{}{𝐫}_\mathrm{𝟏}\stackrel{}{𝐫}_\mathrm{𝟐}`$. We adopt units with $`G=c=m=1`$ and use the same notation and post-Newtonian (PN) equations of motion as described by Kidder kidder93 , expanded to second order (2PN) terms in $`(v/c)^2`$ and $`Gm/r`$, and including spin-orbit (SO) and spin-spin (SS) effects as well as the 2.5PN radiation reaction (RR) term. The individual spins precess because of frame-dragging and the Lens-Thirring effect according to $`\dot{\stackrel{}{𝐒}}_i=\stackrel{}{𝛀}_i\times \stackrel{}{𝐒}_i`$. Here $`i=1,2`$, $`\stackrel{}{𝐒}_i`$ is the spin of the compact object, and $`\stackrel{}{𝛀}_i`$ is a variable axis of precession as defined in eqn. (2.4) of kidder95 . The magnitude of $`\stackrel{}{𝛀}_i`$ is the instantaneous spin precession frequency. We also define a mass ratio $`\beta m_2/m_1`$, and a spin-orbit misalignment angle $`\theta _i`$ for each object as the angle between the spin vector $`\stackrel{}{𝐒}_i`$ and the Newtonian angular momentum $`\stackrel{}{𝐋}_N=\mu (\stackrel{}{𝐫}\times \dot{\stackrel{}{𝐫}})`$. The relative position vector $`\stackrel{}{𝐫}`$ evolves according to a second-order ordinary differential equation of the type $`\ddot{\stackrel{}{𝐫}}=\stackrel{}{𝐚}_{\mathrm{𝐏𝐍}}+\stackrel{}{𝐚}_{\mathrm{𝐒𝐎}}+\stackrel{}{𝐚}_{\mathrm{𝐒𝐒}}+\stackrel{}{𝐚}_{\mathrm{𝐑𝐑}}`$. The full expressions for these terms can be found in Refs. levingrqc  and kidder95 .

We integrate these equations numerically in double precision using a 5th-order Runge-Kutta scheme with an adaptive time step. The robust nature of the Runge-Kutta algorithm makes it particularly attractive for measuring exponential divergence of nearby trajectories. Indeed, since the Runge-Kutta integration can only introduce errors that grow at a polynomial rate, any exponential divergence, if it occurs, will rapidly dominate the evolution and should be easily distinguishable from numerical effects sussman . To quantify chaotic behavior in a dynamical system, the equations of motion must be conservative in the sense that there should be no dissipative terms that could act as attractors in phase space (which would eliminate the possibility of the system being formally chaotic; see cornish ). For this reason, when integrating the equations of motion to calculate Lyapunov exponents, the radiation reaction terms are not included. This also allows for the option of very long integrations to test whether the system is formally chaotic (but on a time scale much longer than the actual inspiral time).

Historically, an important tool for identifying chaos in celestial mechanics has been the Poincaré surface of section murray . By plotting the position in phase space only at certain values of independent coordinate variables, one can reduce a complicated four-dimensional (4-D) phase space trajectory to a 2-D scatter plot. Conserved quantities such as energy or angular momentum generally are held constant for different initial conditions plotted on a single section. For many systems, some initial conditions behave regularly, producing 1-D curves as if there were additional integrals of the motion, while others fill out 2-D regions of the section. This spreading away from the invariant curves is evidence for chaotic behavior in the system. However, this method is only successful in reducing the phase space by one additional dimension, so for systems with \[# of degrees of freedom\]$``$\[# of integrals of motion\] greater than two, the surface of section technique is not very useful for identifying chaos wisdom . For higher-dimensional systems, the projection of non-chaotic orbits onto a two-dimensional section in general will not generate confined curves, and thus a spreading of points in this pseudo-section is not necessarily an indication of chaos. For the problem of two spinning compact objects, the phase space is 10-D (3 from $`\stackrel{}{𝐫}`$, 3 from $`\dot{\stackrel{}{𝐫}}`$, and 2 each from the spin vectors $`\stackrel{}{𝐒}_\mathrm{𝟏}`$ and $`\stackrel{}{𝐒}_\mathrm{𝟐}`$, which can precess in arbitrary orientations but maintain constant magnitude phase\_dim ) while there are only 4 obvious constraints, corresponding to angular momentum and energy conservation due to the invariance of the PN Lagrangian under rotations and time translations kidder93 . This almost guarantees that the projection of the trajectory onto a 2-D section of phase space will not be constrained to a 1-D curve, so regular behavior could easily be misinterpreted as chaos (cf. levinprl ; levingrqc ).

A more quantitative method for identifying and measuring chaos is to calculate the maximum Lyapunov exponent, defined as the divergence rate between initially nearby trajectories:

$$\gamma (t)\frac{1}{t}\mathrm{ln}\left(\frac{dX(t)}{dX(0)}\right).$$
(1)
Here the difference $`dX`$ between two points in phase space is simply the Cartesian distance between the dimensionless 12-component coordinate vectors $`[\stackrel{}{𝐫},\dot{\stackrel{}{𝐫}},\stackrel{}{𝐒}_\mathrm{𝟏},\stackrel{}{𝐒}_\mathrm{𝟐}]`$ and $`[\stackrel{}{𝐫}^{},\dot{\stackrel{}{𝐫}^{}},\stackrel{}{𝐒}_{\mathrm{𝟏}}^{}{}_{}{}^{},\stackrel{}{𝐒}_{\mathrm{𝟐}}^{}{}_{}{}^{}]`$ of two nearby trajectories. In our numerical integrations, the initial separation is a small displacement in phase space with a random orientation and a magnitude of $`dX(0)=10^{10}`$. The two trajectories are then integrated forward in time, recording the separation $`dX(t)`$, from which $`\gamma `$ is computed. In chaotic systems, the divergence will be exponential in time with a roughly constant (positive) exponent: $`dX(t)dX(0)\mathrm{exp}(\gamma t)`$. To be precise, there are actually many different Lyapunov exponents, one for each dimension of phase space. The maximum Lyapunov exponent is automatically selected because, like a classical eigenvector problem, any vector (such as the random initial displacement vector) when multiplied repeatedly by the same matrix will grow fastest in the direction corresponding to the largest eigenvalue. Similarly, for a chaotic system, any random initial displacement in phase space is expected to grow as fast as the maximum Lyapunov exponent. However, for a regular system, the divergence generally grows linearly or at most as a power law in time. In that case, the Lyapunov exponent $`\gamma `$ approaches zero for large $`t`$. We define the Lyapunov time $`t_\mathrm{L}1/\gamma `$ as the time scale on which nearby trajectories separate by a factor of $`e`$ (the “e-folding time”). While it is difficult, on the basis of numerical integrations, to claim that a system is categorically non-chaotic on all time scales, we can make the more precise claim that on a particular time scale the system shows no chaotic behavior, that is, we can set a lower limit on $`t_\mathrm{L}`$. For the problem of coalescing compact binaries, the relevant time scale is the inspiral time $`t_{\mathrm{insp}}`$, so if $`t_{\mathrm{insp}}t_\mathrm{L}`$, chaos will not affect the dynamics.

To establish the validity of our numerical approach, we have re-examined the problem studied by Suzuki and Maeda suzuki , where both chaotic and regular trajectories have clearly been identified for super-maximally spinning test particles ($`S_2m_2^2`$; $`m\mu `$) orbiting around a Schwarzschild black hole ($`S_1=0`$). While physically unrealistic, these conditions introduce no mathematical singularities into the equations of motion and are therefore fine examples for identifying chaotic or regular behavior in this dynamical system. Here we show numerical results for two different sets of initial conditions: both begin with $`r/m=4.5`$, $`\beta =10^4`$, $`S_1=0`$, $`S_2=1.4\times 10^4m_2^2`$, and spin-orbit misalignment angle $`\theta _2=\pi /4`$, but one trajectory has an initial Newtonian orbital angular momentum of $`L_N/m=1.485\mu `$ while the other begins with $`L_N/m=1.53\mu `$. In Fig. 1 we show the measured value of $`\gamma `$ as a function of time (measured in orbital periods). The orbit with $`L_N/m=1.485\mu `$ (a) is clearly chaotic with $`\gamma 0.1`$ (corresponding to a divergence time scale of $`t_\mathrm{L}10`$ orbital periods). The other trajectory (b) exhibits very different qualitative behavior, with no evidence for chaos even on the time scale of thousands of orbits. Our results for these and many other trajectories that we have computed agree well with Suzuki and Maeda’s results, confirming that both chaotic and regular orbits exist for rapidly spinning test particles around a Schwarzschild black hole.

Another method for distinguishing between chaotic and quasi-periodic orbits is to look for certain qualitative features in the power spectrum of one of the dynamical variables sussman ; sussman92 . For the compact binary problem, the obvious choice is to analyze the spectrum of the GW signal, defined as the squared amplitude of the Fourier transform of the GW strain $`h(t)`$. With RR turned off, the spectrum should exhibit sharp lines for quasi-periodic orbits, but broad-band noise for chaotic orbits. Including only the leading quadrupole radiation terms, the components of the transverse-traceless GW tensor observed from the $`\widehat{z}`$ direction are

$$h_+=\frac{4\mu }{D}\left(v_x^2\frac{m}{r^3}x^2\right),$$
(2)
$$h_\times =\frac{4\mu }{D}\left(v_xv_y\frac{m}{r^3}xy\right),$$
(3)
where $`D`$ is the distance to the source and $`v_i=\dot{x}_i`$ kidder95 ; mtw . Figure 2 shows the power spectra in $`h_+(t)`$ for the same two test-particle orbits as in Fig. 1. Both spectra show peaks around $`f_{\mathrm{GW}}25`$Hz, the fundamental quadrupole frequency (twice the orbital frequency) for these initial conditions. As expected, the chaotic orbit (a) produces primarily broad-band noise with few distinguishable features, whereas the spectrum for the quasi-periodic orbit (b) shows many sharp lines.

Having established the ability of our method to distinguish between regular and chaotic trajectories described by the PN equations of motion, we now apply the same techniques to the astrophysically relevant systems expected to be detectable by ground-based laser interferometers. Here we will show results for an illustrative case of two maximally spinning $`10M_{}`$ black holes ($`S_i=m_i^2`$) in an initially circular orbit with spin-orbit misalignment angles of $`\theta _1=38^o`$ and $`\theta _2=70^o`$. To measure the Lyapunov exponent at different points during the final inspiral, 4 different trajectories were integrated, corresponding to GW frequencies of 10, 40, 100, and 400 Hz, emitted at orbital separations $`r/m`$ 50, 20, 10, and 5, respectively initnote . For each separation, we then integrate the PN equations of motion (with RR turned off) over a period much longer than the inspiral time. If no exponential divergence is observed, we can safely conclude that the system is, for all practical purposes, not chaotic.

The measured values for $`\gamma (t)`$ are plotted in Fig. 3 for each of the selected stages of the inspiral. For any power-law divergence in phase space with $`dX(t)t^\alpha `$, on a log-log plot of $`\gamma `$ versus time the slope is $`d\mathrm{ln}\gamma /d\mathrm{ln}t=1/\mathrm{ln}t1`$, so that for large times the curve should be linear with slope -1 \[as we see in Fig. 1 (b)\]. All plots in Fig. 3 clearly show this behavior, characteristic of regular, non-chaotic orbits. Also shown for comparison in Fig. 3 are the inspiral times corresponding to each separation. In all four cases, the lower limit on $`t_\mathrm{L}`$ far exceeds $`t_{\mathrm{insp}}`$.

Just as in the test particle example presented above, we also calculate the gravitational wave power spectrum for this black hole binary. Here the results are even more striking: Fig. 4 shows the $`h_+(t)`$ power spectrum at the point in the inspiral where the GW frequency is $`100`$Hz. The spectrum has a few very sharp lines with no evidence for broad-band noise. The strongest features correspond to the fundamental GW quadrupole and spin-orbit precessional frequencies ($`100`$Hz and $`7`$Hz, respectively) and their harmonics. The spectra from the other points in the inspiral show the same qualitative features, all indicative of quasi-periodic motion.

In addition to this prototypical binary black hole system, we have investigated the behavior of a large number of other systems at different stages throughout the inspiral. Varying the binary mass ratio, spin magnitudes (always with $`S_im_i^2`$), misalignment angles, eccentricity, and initial separation over wide ranges, we consistently find the same regular, non-chaotic behavior for all trajectories. In particular, we have calculated orbits for the same high-eccentricity systems considered by Levin in Ref. levinprl  and again measure only linear divergence of nearby initial conditions, finding $`t_{\mathrm{insp}}t_\mathrm{L}`$ belczynski . We have also looked at many binary systems in which either one or both objects are solar-mass neutron stars. Because of their smaller mass (and thus slower radiation loss), these systems can spend significantly longer times in the inspiral band, completing tens of thousands of orbits before merger. Yet, even on the longest relevant time scales ($`t_{\mathrm{insp}}10^4`$s), we again find no chaotic behavior.

Recent work in binary star evolution kalogera  suggests that systems with large spin misalignment angles are likely to be astrophysically relevant sources, giving rise to complicated nonlinear SO and SS interactions. It is possible that there are undiscovered narrow regions in phase space where specific resonant conditions could give rise to chaotic behavior, yet even if a binary system passes through them during inspiral, the time spent in these resonant bands will most likely be much shorter than the Lyapunov time. Even though the Lyapunov exponent approaches zero for all the systems we have considered, their orbits can still appear quite irregular, especially when both objects are spinning rapidly levinprl . As mentioned above, this can lead to the spreading of the surface of section away from a 1-D curve, but in dynamical systems with many degrees of freedom, this spreading alone does not imply the presence of chaos. Nevertheless, the detection of GW signals from binaries with rapidly spinning components may be very challenging. Systems with similar initial conditions still may produce waveforms that look quite different even without exponential divergence. Orbital precession will add many complicated features to the basic inspiral waveforms. The current template database will most likely need to be extended (perhaps by orders of magnitude) to include at least a portion of parameter space where the signals are significantly modified by spin effects apost . Furthermore, for the late stages of inspiral, the 2.5PN expansion becomes invalid and one should include higher order PN terms (or even perform calculations in full general relativity buonanno ). While these terms introduce additional features in the dynamics, at this late stage the orbit is certainly decaying too rapidly for the final GW signal to be affected by formal chaotic behavior.

###### Acknowledgements.

We are very grateful to J. Wisdom for many useful comments and discussions. We also thank V. Kalogera and J. Levin for helpful discussions. This work was supported by NSF Grants AST-9618116 and PHY-0070918 and the NSF Graduate Research Fellowship Program.


format: text

text:
# Two-particle pairing and phase separation in a two-dimensional Bose-gas with one or two sorts of bosons

## I Introduction

It is well known that, in contrast with two-particle Cooper pairing in Fermi systems, the essence of a superfluidity in Bose systems is one particle Bose-Einstein condensation (BEC). This asymmetry between Fermi and Bose ( two- particle versus one-particle condensation) was challenged in a pioneering paper by Valatin and Butler . They proposed a BCS-like variational function for the description of an attractive Bose gas. The most difficult problem with the validity of their description is connected with a tendency toward phase separation which arises in attractive Bose systems. Later on Nozieres and Saint-James  conjectured that in a Bose system with a short-ranged, hard core repulsion and a van der Waals attractive tail, in principle, it is possible to create a two-particle bosonic bound state and to escape collapse. Unfortunately their calculations in a three-dimensional system showed that at least for one sort of structureless boson either standard one-particle BEC is more energetically beneficial, or that phase separation takes place earlier than the two-particle condensation. Note that the same result was obtained earlier by Jordanski  for the case of weak van der Waals attraction.

The important development of the ideas of Nozieres and Saint James belongs to Rice and Wang . The authors claimed that in two-dimensional (where already an infinitely small attraction leads to the bound state in a symmetrical potential well) it is possible to realize a two-particle boson pairing. Moreover, this two-particle pairing results, for small momenta $`p\xi _0<1`$ in a linear, soundlike, dispersion law of quasiparticles at $`T=0`$ in an analogy with a standard one-particle Bose-condensation.

To escape a collapse in a 2D attractive Bose gas, the authors of Ref. introduced in their model a Hartree-Fock shift of the chemical potential $`\mu _BUn`$, connected with the short-range repulsion $`U`$. This shift in the case of $`U>2V`$, where $`V`$ is the magnitude of the van der Waals tail, leads to a positive compressibility in the system $`\kappa ^1=d\mu _B/dn=U2V>0`$.

The main goal of a present paper is to construct a phase diagram of a 2D dilute Bose gas with the van der Waals interaction between particles, by taking into account on equal grounds the full contribution of a hard-core repulsion $`U`$ and a van der Waals tail $`V`$. Throughout the paper we will consider the lattice model, and will base our results on the exact solution of the two-particle $`T`$-matrix problem presented in Refs. and . We will study the possibility of different two-boson pairings, as well as the possibility of a total phase separation in the system. We will also consider the two sorts of structureless bosons described by the two-band Hubbard model. In the case of attraction between bosons of two different sorts, we will find a possibility of an $`s`$-wave two-boson pairing $`b_1b_20`$.

## II Theoretical model

The model under consideration is described by the following Hamiltonian on the 2D square lattice:

$$H=t\underset{ij}{}b_i^{}b_j+\frac{U}{2}\underset{i}{}n_i^2\frac{V}{2}\underset{ij}{}n_in_j,$$
(1)
where $`n_i=b_i^{}b_i`$ is a 2D boson density. We will work in the limit of strong hard-core repulsion $`U\{V;t\}`$, and restrict ourselves mostly to a low-density limit $`n_Bd^21`$, $`d`$ being the interatomic distance. Further on we put $`d=1`$.

Note that in the case of $`V=0`$ model (1) is just the Bose-Hubbard model, extensively studied in the literature for the case of 2D 4He submonolayers, as well as for the flux lattices and Josephson arrays in the type-II superconductors (see Refs.).

We would like to emphasize that model (1) is, to some extent, a Bose analog of the fermionic $`tJ`$ model with released constraint, considered by Kagan and Rice . After the Fourier transformation from Eq. (1) we obtain:

$`H`$ $`=`$ $`{\displaystyle \underset{𝐩}{}}\epsilon _𝐩b_𝐩^{}b_𝐩+{\displaystyle \frac{U}{2}}{\displaystyle \underset{𝐤_\mathrm{𝟏}𝐤_\mathrm{𝟐}𝐪}{}}b_{𝐤_\mathrm{𝟏}}^{}b_{𝐤_\mathrm{𝟐}}^{}b_{𝐤_\mathrm{𝟐}𝐪}b_{𝐤_\mathrm{𝟏}+𝐪}`$ (2)
$``$ $`{\displaystyle \underset{𝐤_\mathrm{𝟏}𝐤_\mathrm{𝟐}𝐪}{}}V(𝐪)b_{𝐤_\mathrm{𝟏}}^{}b_{𝐤_\mathrm{𝟏}+𝐪}b_{𝐤_\mathrm{𝟐}}^{}b_{𝐤_\mathrm{𝟐}𝐪},`$ (3)
where $`\epsilon _𝐩=2t(\mathrm{cos}p_x+\mathrm{cos}p_y)`$ is a bosonic spectrum on the square lattice and $`V(𝐪)=V(\mathrm{cos}(q_x)+\mathrm{cos}(q_y))`$ is a Fourier transform of the Van der Waals tail. As a result, a total interaction in the momentum space is given by the formula:

$$V_{eff}(𝐪)=\frac{U}{2}V(𝐪).$$
(4)
## III The $`T`$-matrix problem.

An instability toward a two-particle boson pairing manifests itself in the appearance of a pole at a temperature $`T=T_c`$ in the solution of the Bethe-Salpeter equation for the two-particle vertex $`\mathrm{\Gamma }`$ for zero total momentum of the two bosons $`(𝐩,𝐩)`$. To proceed to the solution of this equation, we must solve at first the $`T`$-matrix problem for the two particles in vacuum. Here we can use the results of Ref., because the solution of the two-particle problem does not depend upon statistics of colliding particles. For the $`T`$-matrix problem it is convenient to expand $`V_{eff}(𝐪)`$ in Eq.(4) in series with the eigenfunctions of the irreducible representation of the lattice symmetry group $`D_4`$. This yields

$$\begin{array}{c}V_{eff}(\text{extended }s\text{-wave})=\frac{U}{2}\frac{V}{2}(\mathrm{cos}p_x+\mathrm{cos}p_y)(\mathrm{cos}p_x^{}+\mathrm{cos}p_y^{}),\hfill \\ V_{eff}(p\text{-wave})=\frac{V}{2}(\mathrm{sin}p_x\mathrm{sin}p_x^{}+\mathrm{sin}p_y\mathrm{sin}p_y^{}),\hfill \\ V_{eff}(d_{x^2y^2})=\frac{V}{2}(\mathrm{cos}p_x\mathrm{cos}p_y)(\mathrm{cos}p_x^{}\mathrm{cos}p_y^{}),\hfill \end{array}$$
(5)
where $`𝐪=𝐩𝐩^{}`$ is a transfered momentum. Note that, for spinless bosons, which we formally consider in Eq. (1), the total spin of the Bose pair is zero. Hence only $`s`$\- and a $`d`$-wave pairings are allowed by the symmetry of the pair $`\mathrm{\Psi }`$ function. A $`p`$-wave pairing is allowed only for an odd total spin of the two bosons. Nevertheless we will conserve the results for the $`p`$-wave pairing in our paper because the generalization of Eq.(1) for the case of bosons with internal degrees of freedom is straightforward. The $`T`$-matrix problems for $`p`$\- and $`d`$-wave channels are very simple. Solutions of these problems for the two particles with a total momentum zero and a total energy $`E`$ yield

$$T_{d;p}(E)=\frac{\frac{1}{2}V}{\left(1+\frac{1}{2}VI_{d;p}\right)},$$
(6)
where

$$I_{d;p}=\underset{0}{\overset{2\pi }{}}\underset{0}{\overset{2\pi }{}}\frac{dp_xdp_y}{(2\pi )^2}\frac{|\varphi _{d;p}|^2}{E+4t(\mathrm{cos}p_x+\mathrm{cos}p_y)},$$
(7)
and the functions $`\varphi _d=\mathrm{cos}p_x\mathrm{cos}p_y`$ and $`\varphi _p=\mathrm{sin}p_x+i\mathrm{sin}p_y`$ correspond, respectively, to the $`d`$\- and $`p`$-wave channels.

Let us find the thresholds for the bound states in the $`d`$\- and the $`p`$-wave channels. The appearance of a bound state means that $`E=W|\stackrel{~}{E}|`$, where $`W=8t`$ is a bandwidth. For the threshold $`\stackrel{~}{E}=0`$. An exact solution of Eqs. (6) and (7), which involves the calculation of elliptic integrals of first and second order, yields

$$\begin{array}{c}\left(\frac{V_c}{4t}\right)_{p\text{wave}}=\frac{1}{1{\displaystyle \frac{2}{\pi }}}2.8,\hfill \\ \left(\frac{V_c}{4t}\right)_{d\text{wave}}=\frac{1}{{\displaystyle \frac{4}{\pi }}1}3.7.\hfill \end{array}$$
(8)
Note that a threshold for a $`p`$-wave pairing is lower. Now let us proceed to an $`s`$-wave channel. Here an ordinary $`s`$-wave pairing is suppressed by large hard-core repulsion $`U`$, however an extended $`s`$-wave pairing with a symmetry of the order parameter $`\mathrm{\Delta }_s=\mathrm{\Delta }_0(\mathrm{cos}p_x+\mathrm{cos}p_y)`$ is allowed. In real space this pairing corresponds to the particles on the neighboring sites. Moreover the pair $`\mathrm{\Psi }`$ function is zero in the region of a hard core ($`r<r_0`$) and is centered (has a maximum) in the region of a van der Waals attraction for $`rr_1`$ (see Fig. 1). On the lattice $`r_0d/2`$ and $`r_1d`$.

One can see that the $`\mathrm{\Psi }`$ function has a region of zero value zero. But it has no nodes because it does not change its sign for all values of $`r`$. The rigorous calculation of the threshold for an extended $`s`$-wave pairing gives

$$\left(\frac{V_c}{4t}\right)_{s\text{wave}}=1.$$
(9)
Moreover for $`V>V_c`$ an energy of the bound state has the form

$$|E_b|=|\stackrel{~}{E}|=8W\mathrm{exp}\left\{\frac{\pi V}{(VV_{crit})}\right\}.$$
(10)
Of course, in a strong coupling case $`VW`$

$$|E_b|V.$$

The $`T`$ matrix in an $`s`$-wave channel for small and intermediate values of $`V`$ is given by

$$T_s(|\stackrel{~}{E}|)=\frac{W\left(1{\displaystyle \frac{V}{4t}}\right)}{{\displaystyle \frac{1}{\pi }}\left(1{\displaystyle \frac{V}{4t}}\right)\mathrm{ln}{\displaystyle \frac{8W}{|\stackrel{~}{E}|}}{\displaystyle \frac{V}{4t}}}.$$
(11)
The most important is that a strong Hubbard repulsion $`U`$ acts only as an excluded volume and effectively drops out from Eq. (11) at low energies. For $`V4t`$ the $`T`$ matrix

$$T_s(|\stackrel{~}{E}|)\frac{\pi W}{\mathrm{ln}\frac{8W}{|\stackrel{~}{E}|}}$$

corresponds to repulsion and coincides with the $`T`$ matrix for the 2D Bose-Hubbard model at low density. Of course, at very high energies the $`T`$ matrix contains an additional pole $`EU`$ which corresponds to an antibound state in a total analogy with the fermionic Hubbard model. For $`V=4t`$: $`T_s(|\stackrel{~}{E}|)=0`$ and there is no interaction at all. Finally, for $`V>4t`$, $`T_s(|\stackrel{~}{E}|)<0`$ corresponds to an attraction and reflects the appearance of the bound state.

## IV Bethe-Salpeter equation for an $`s`$-wave pairing of the two bosons

Let us consider at first the most interesting case of $`V>4t`$ and find the critical temperature for an extended $`s`$-wave pairing of the two bosons. The solution of the Bethe-Salpeter equation for bosonic systems reads

$$\mathrm{\Gamma }_s=\frac{T_s}{1+T_s{\displaystyle \frac{dp_xdp_y}{(2\pi )^2}\frac{\mathrm{coth}\frac{\epsilon _𝐩\mu }{2T}}{2(\epsilon _𝐩\mu )}}}.$$
(12)
For a low density of bosons $`n_B1`$ one has

$$\epsilon _𝐩4t+\frac{p^2}{2m},$$

where $`m=1/2t`$ is a boson mass. Accordingly $`\mu =4t+\stackrel{~}{\mu }`$ and $`\xi _𝐩=\epsilon _𝐩\mu =p^2/2m\stackrel{~}{\mu }.`$

The most substantial difference of Eq. (12) from an analogous fermionic equation is the replacement of $`\mathrm{tanh}\xi /2T`$ by $`\mathrm{coth}\xi /2T`$ in its kernel. Moreover, as shown in Ref. for the 2D Fermi gas $`\stackrel{~}{\mu }=\epsilon _F|E_b|/2`$. So, in a weak-coupling case, when $`\epsilon _F|E_b|`$, the chemical potential $`\stackrel{~}{\mu }\epsilon _F>0`$ is positive. In contrast to this we shall see below that a bosonic chemical potential $`\stackrel{~}{\mu }`$ is always negative even in the weak-coupling case, when a binding energy is much smaller than a degeneracy temperature:

$$|E_b|T_0=\frac{2\pi n}{m}.$$

Another very important point is that the $`T`$ matrix, which enters into the Bethe-Salpeter equation, must be calculated for a total energy $`\stackrel{~}{E}=2\stackrel{~}{\mu }`$ of colliding bosons. The chemical potential $`\stackrel{~}{\mu }`$ can be determined from the requirement of the number of particle conservation. This requirement yields

$$n_B=\frac{d^2p}{(2\pi )^2}\frac{1}{\mathrm{exp}\left\{\frac{p^2/2m\stackrel{~}{\mu }}{T}\right\}1}.$$
(13)
From Eq.(13) for the temperatures $`|E_b|<T<T_0<W`$ we obtain

$$\stackrel{~}{\mu }=T\mathrm{exp}\left(\frac{T_0}{T}\right)<0.$$
(14)
Note that a standard Hartree-Fock shift $`nU`$ drops out from the expression for $`\xi _𝐩=\epsilon _𝐩\mu `$ both in the Bethe-Salpeter equation (12) and in the equation for the number of particle conservation (13). Now we are ready to solve the Bethe-Salpeter equation (12). The critical temperature $`T_c`$ corresponds to the pole in Eq. (12)

$$1+\frac{mT_s(2\stackrel{~}{\mu })}{2\pi }I=0,$$
(15)
where

$$I=\underset{0}{\overset{W/T_c}{}}𝑑y\frac{\mathrm{coth}\left(y+{\displaystyle \frac{|\stackrel{~}{\mu }|}{2T_c}}\right)}{y+{\displaystyle \frac{|\stackrel{~}{\mu }|}{2T_c}}},$$
(16)
and

$$y=\frac{p^2}{4mT_c}.$$

An analysis of Eq.(16) shows that the main contribution to the integral comes from the lower limit of integration.

Hence providing $`|\stackrel{~}{\mu }|/T_c1`$ we have

$$I\underset{0}{\overset{W/T_c}{}}\frac{dy}{\left(y+{\displaystyle \frac{|\stackrel{~}{\mu }|}{2T_c}}\right)^2}\frac{2T_c}{|\stackrel{~}{\mu }|}.$$
(17)
As a result Eq. (15) can be represented in the following form:

$$\frac{T_c}{|\stackrel{~}{\mu }|}=\frac{\pi }{mT_s(2\stackrel{~}{\mu })}.$$
(18)
It is useful now to represent $`T_s(2\stackrel{~}{\mu })`$ in terms of the binding energy $`E_b`$. Utilizing Eqs. (10) and (11) we can write

$$T_s(2\stackrel{~}{\mu })=\frac{\pi W}{\mathrm{ln}\frac{2|\stackrel{~}{\mu }|}{|E_b|}}=\frac{4\pi }{m\mathrm{ln}\frac{2\stackrel{~}{|\mu |}}{|E_b|}}$$
(19)
It is important to mention here that $`\stackrel{~}{\mu }<0`$, and hence the $`T`$ matrix in Eq. (19) does not contain an imaginary part. In the fermionic case $`\stackrel{~}{\mu }=\epsilon _F>0`$, and the $`T`$ matrix contains an imaginary part corresponding to the resonant scattering. As a result, from Eq. (19) we obtain:

$$4\mathrm{ln}\frac{2|\stackrel{~}{\mu }|}{E_b}=\frac{|\stackrel{~}{\mu }|}{T_c}.$$
(20)
Assuming that $`|E_b|T_cT_0`$ , we obtain

$$\stackrel{~}{\mu }(T_c)=T_c\mathrm{exp}\left(\frac{T_0}{T_c}\right)$$

and

$$\frac{|\stackrel{~}{\mu }|}{T_c}=\mathrm{exp}\left\{\frac{T_0}{T_c}\right\}.$$

Later on we will justify this assumption.

As a result from Eq. (20) we will obtain

$$T_c=\frac{T_0}{\mathrm{ln}\left(\frac{1}{4}\mathrm{ln}\frac{2|\stackrel{~}{\mu }|}{|E_b|}\right)}.$$
(21)
Recall that in the case of the fermionic pairing in two dimensionsa critical temperature reads

$$T_c=\sqrt{2\epsilon _F|E_b|}.$$

Let us analyze expression (21). As we already know

$$|E_b|=8W\mathrm{exp}\left\{\frac{1}{\lambda }\right\},$$

where

$$\lambda =\frac{(VV_{cs})}{\pi V}.$$
(22)
Then a condition $`|E_b|T_0`$ means

$$\lambda \frac{1}{\mathrm{ln}\frac{W}{T_0}}1.$$
(23)
Hence $`\mathrm{ln}T_0/|E_b|=1/\lambda \mathrm{ln}W/T_01/\lambda `$, and

$$T_c\frac{T_0}{\mathrm{ln}\frac{1}{4}\mathrm{ln}\frac{T_0}{|E_b|}}\frac{T_0}{\mathrm{ln}\frac{1}{4\lambda }},$$
(24)
which is in an agreement with Ref.. Note that $`T_c`$ from Eq. (24) satisfies the conditions $`|E_b|T_cT_0`$, so an assumption used for the derivation of $`T_c`$ is justified.

For $`T<T_c`$ the spectrum of the quasiparticles acquires a gap:

$$E_𝐩=\sqrt{(\epsilon _𝐩+|\stackrel{~}{\mu }|)^2\mathrm{\Delta }^2}.$$
(25)
Note that at low densities of bosons a gap $`\mathrm{\Delta }`$ becomes isotropic in the principal approximation.

The gap $`\mathrm{\Delta }`$ together with the chemical potential $`\stackrel{~}{\mu }`$ must be defined self-consistently from the two coupled equations

$$1=\frac{\lambda }{4}\underset{\frac{|\stackrel{~}{\mu }|}{2T}}{\overset{\frac{W}{2T}}{}}𝑑z\frac{\mathrm{coth}\sqrt{z^2\frac{\mathrm{\Delta }^2}{4T^2}}}{\sqrt{z^2\frac{\mathrm{\Delta }^2}{4T^2}}},$$
(26)
$$n_B=\frac{\lambda }{4}\underset{|\stackrel{~}{\mu }|}{\overset{W}{}}𝑑\xi \frac{1}{\mathrm{exp}\left\{\frac{\sqrt{\xi ^2\mathrm{\Delta }^2}}{T}\right\}1},$$
(27)
where $`\xi =\epsilon +|\stackrel{~}{\mu }|`$ and $`z=\xi /2T`$.

Of course, the solution of the system of equations Eqs. (26) and (27) exists only if $`|\stackrel{~}{\mu }|\mathrm{\Delta }`$, or, in other words, only if $`E_g^2=|\stackrel{~}{\mu }|^2\mathrm{\Delta }^2>0`$. The exact solution of these equations yields for zero temperature in an agreement with

$$|\stackrel{~}{\mu }(T=0)|=\mathrm{\Delta }=\frac{|E_b|}{2}.$$
(28)
This result is very important. It justifies our scenario, leading to a linear, soundlike spectrum of the quasiparticles for a small momenta $`p`$. Indeed

$$E_𝐩=\sqrt{\epsilon _𝐩^2+2\epsilon _𝐩|\stackrel{~}{\mu }|}=\sqrt{\epsilon _𝐩^2+\epsilon _𝐩|E_b|}.$$
(29)
From Eq. (29) for the case $`p\xi _01`$, where $`\xi _0=1/\sqrt{2m|E_b|}`$ is the coherence length of the boson pair, we immediately obtain a linear dispersion law:

$$E_p=cp.$$
(30)
In (30) $`c^2=|E_b|/2m`$ is a sound velocity squared. This means that an inverse compressibility of the system $`\kappa ^1=c^2`$ is positive. This fact proves the stability of a superfluid paired state and excludes the possibility of the collapse of the pairs in the system. Note also that close to $`T_c`$ one has

$$\mathrm{\Delta }(T)\mathrm{\Delta }(0)\sqrt{\frac{T_cT}{T_c}},$$
(31)
which is similar to the BCS theory. We would like to mention that bosonic pairs in the limit $`|E_b|T_0`$ are extended in full analogy with the BCS theory. That is the coherence length in this limit,

$$\xi _0\frac{1}{\sqrt{n}}1,$$
(32)
is larger than the mean distance between the bosons. The Bose pairs are strongly overlapping in this limit. The pairing takes place in the momentum space in an analogy with the Cooper pairing in the BCS picture of superconductivity.

In the opposite limit $`|E_b|T_0`$, the situation closely resembles the bipolaronic limit for the fermionic systems . That is, the creation of the bosonic bound pairs is associated with the crossover temperature

$$T^{}=\frac{|E_b|}{\mathrm{ln}(1/n)}.$$
(33)
The Bose condensation of the pairs occurs at lower temperature :

$$T_c=\frac{T_0}{\mathrm{ln}\mathrm{ln}(1/n)}.$$
(34)
Note that this temperature is obtained from the condition $`\mu _b(T_c)=T_c\mathrm{exp}\left(T_0/T_c\right)+f_0T_c=0`$, where $`f_0=1/\mathrm{ln}(1/n)`$ is a repulsive interaction between the local pairs. So the superfluid transition takes place only for a residual repulsion between the pairs. Also note that in a dilute Bose gas the Berezinski-Kosterlitz-Thouless contribution of vortices  is important only very close to $`T_c`$, so the mean field expression (34) gives a very good estimate for the exact BKT critical temperature:

$$(T_cT_{BKT})/T_c1/\mathrm{ln}\mathrm{ln}(1/n)1.$$

In the case of the local pairs the coherence length is small:

$$\xi _0\frac{1}{\sqrt{n}}.$$
(35)
The pairs are compact, and the pairing takes place in the real space.

## V Possibility of $`p`$-wave and $`d`$-wave pairing of the two bosons

Now let us analyze the solution of the Bethe-Salpeter equation for $`p`$\- and $`d`$-wave boson pairings. Here the critical temperatures should be found from the conditions

$$1+T_{p;d}(2\stackrel{~}{\mu })\stackrel{~}{I}_{p;d}=0,$$
(36)
where

$$\stackrel{~}{I}_{p;d}=\underset{0}{\overset{2\pi }{}}\underset{0}{\overset{2\pi }{}}\frac{dp_xdp_y}{(2\pi )^2}\frac{\mathrm{coth}\frac{\epsilon \mu }{2T_c}}{2(\epsilon \mu )}|\varphi _{p;d}|^2.$$
(37)
In a low-density limit the $`\varphi `$ functions can be approximated by the following expressions:

$$\varphi _pp_x+ip_y=pe^{i\varphi },$$

$$\varphi _d\frac{1}{2}(p_x^2p_y^2)=\frac{1}{2}p^2\mathrm{cos}2\varphi .$$

Hence after an angular integration we obtain:

$$\begin{array}{c}\stackrel{~}{I}_p=\frac{m}{2\pi }p𝑑p\frac{\mathrm{coth}\frac{\xi }{2T_c}}{2\xi }p^2,\hfill \\ \stackrel{~}{I}_d=\frac{m}{16\pi }p𝑑p\frac{\mathrm{coth}\frac{\xi }{2T_c}}{2\xi }p^4,\hfill \end{array}$$
(38)
where again $`\xi =p^2/2m+|\stackrel{~}{\mu }|`$.

Additional factors $`p^2`$ and $`p^4`$ in the integral expressions for $`\stackrel{~}{I}_p`$ and $`\stackrel{~}{I}_d`$ reflect a well-known fact, that for slow 2D particles in vacuum an $`s`$-wave harmonics of the scattering amplitude behaves as $`f_0\mathrm{ln}1/p^2`$, whereas for a magnetic number $`m0`$, the scattering amplitude vanishes for $`p`$ goes to zero as $`f_mp^{2m}`$. The additional factor $`p^4`$ leads to the absence of an infra-red singularity for $`\epsilon 0`$ in $`\stackrel{~}{I}_d`$:

$$\stackrel{~}{I}_d\frac{d\epsilon \epsilon ^2}{\epsilon ^2}\epsilon 0.$$
(39)
For the $`p`$-wave channel the infra-red singularity becomes logarithmically weak:

$$\stackrel{~}{I}_p\frac{d\epsilon \epsilon }{\epsilon ^2}\mathrm{ln}\epsilon .$$
(40)
This means that the Bethe-Salpeter equation has no solutions in $`p`$\- and $`d`$-wave channels for $`|V|/t<1`$.

Hence the boson pairing with a large coherence length $`\xi _0>1/\sqrt{n}`$ is absent in a $`p`$-wave channel as well as in a $`d`$-wave channel. Here only the limit of the local pairs is possible. For $`p`$\- and $`d`$\- wave channels local pairs are created at the crossover temperatures given by Eq. (33). In this case the binding energies for $`p`$\- and $`d`$-wave channels read, correspondingly,

$$|E_b^p|W\frac{(VV_{cp})}{V}\frac{1}{\mathrm{ln}\frac{V}{(VV_{cp})}},$$
(41)
$$|E_b^d|W\frac{(VV_{cd})}{V}.$$
(42)
Note that $`V_{cp}`$ and $`V_{cd}`$ in Eqs. (41) and (42) are the thresholds for the $`p`$\- and $`d`$-wave pairings given by Eq. (8). We would also like to mention that for a fixed $`V`$: $`|E_b^p|>|E_b^d|`$. Providing that the interaction between the local pairs is repulsive, the temperature of the Bose condensation of the local pairs in these two cases is again given by Eq. (34). Summarizing this case we see that, while the crossover temperatures $`T^{}`$ are different for the $`s`$-, $`p`$\- and $`d`$-wave local pairs, their critical temperatures coincide with each other in a principal order.

## VI Total phase separation

As we discussed in the Sec. 1, the real collapse is prohibited in our system by large Hubbard repulsion $`U`$. However the phase separation on the two large clusters is allowed. The first cluster corresponds to the Mott-Hubbard Bose solid. In this cluster $`n_b1`$, that is, each site is practically occupied by one boson. Such a cluster is localized due to Mott-Hubbard considerations. It has no kinetic energy. However it has a potential energy of the order of $`2V`$ for one particle. A second cluster has a very small boson density $`n0`$. In this cluster for $`V<4t`$ the energy per particle is $`W/2+f_0n`$, where $`f_0=m/4\pi \mathrm{ln}1/n`$ is a two-particle $`T`$ matrix in the absence of a bound state. Rigorously speaking, at a given bosonic density $`n`$ the phase separation results in the formation of the two clusters with the densities $`n_1>n`$ and $`n_2<n`$, where $`n_1`$ is close to or identically equal to 1. The phase separation for $`V<4t`$ takes place if the energy per particle in the cluster with the density $`n_1`$ becomes smaller than the energy per particle in the cluster with the density $`n_2`$

$$2\eta VW/2,$$
(43)
where $`\eta `$ is an unknown numerical coefficient of the order of $`1`$. Note that in the fermionic $`tJ`$ model, considered in Ref., the Mott-Hubbard cluster with $`n=1`$ has an antiferromagnetic order. Hence instead of $`\eta V`$ in (43) one should write $`1.18J`$ – the energy per spin of the 2D antiferromagnetic. As a result, in a fermionic case $`J_{ps}=3.8t`$. In our system $`V_{ps}2t`$, due to the absence of kinetic energy and zero point energy in the case of structureless bosons. In the same time for $`n1`$ the phase separation between the Bose solid and the one-particle BEC takes place. According to Ref. for $`n_B1`$, the phase separation takes place already for small values of $`|V|/t`$.

In principle, another scenario of the phase separation connected with the creation of quartets  is also possible in our system. It requires an evaluation of the four-particle vertex which is impossible to do analytically. However we think that our scenario of total phase separation takes place for smaller values of $`V/t`$ than the quartet formation. This is in agreement with numerical calculations  for the 2D fermionic $`tJ`$ model.

## VII Phase diagram of the system

In this section we will complete the phase diagram of the system. At first, note that for $`V<V_{ps}`$ (when the $`T`$ matrix for an $`s`$-wave channel is repulsive) we have at low density a standard Bose gas with a hard-core repulsion. It will be unstable toward a standard one-particle BEC at a critical temperature given again by Eq. (34). For $`V>2t`$ a total phase separation on two large clusters takes place in our system. One of these clusters contains a Mott-Hubbard Bose solid, another one contains a Bose gas with one particle condensation.

For large densities $`n=n_c1`$ ($`n_c1`$ in Ref. for structureless bosons) the system will undergo a transition to the Mott-Hubbard Bose solid. As a result, on a qualitative level the phase diagram for our system has the form, presented in Fig.2. Note that our phase diagram should be important for the understanding of the physics of the gas of kinks and steps on a solid interface of 4He. It will be interesting in this context to obtain phase separation in a system of kinks on a quantum atomically rough surface of 4He. Another possible application of our results can be connected with the study of biexcitonic pairing in semiconductors . In this context we must mention an interesting recent paper of Lozovik et al . .

## VIII Two-band Hubbard model for the two sorts of bosons

Let us consider the two-band Hubbard model for the two sorts of structureless bosons. The Hamiltonian of the system has the form

$`H=`$ $`t_a{\displaystyle \underset{ij}{}}a_i^{}a_jt_b{\displaystyle \underset{ij}{}}b_i^{}b_j`$ (45)
$`+{\displaystyle \frac{U_{aa}}{2}}{\displaystyle \underset{i}{}}n_{ia}^2+{\displaystyle \frac{U_{bb}}{2}}{\displaystyle \underset{i}{}}n_{ib}^2{\displaystyle \frac{U_{ab}}{2}}{\displaystyle \underset{i}{}}n_{ia}n_{ib},`$
where $`t_a`$ and $`t_b`$ and $`n_a`$ and $`n_b`$ are, respectively, the hopping matrix elements and densities for bosons of sorts $`a`$ and $`b`$. For simplicity we will consider the case $`t_a=t_b`$, corresponding to the equal masses $`m_a=1/(2t_a)=m_b`$. We also assume that the bottoms of the bands coincide. In the Hamiltonian $`U_{aa}`$ and $`U_{bb}`$ are Hubbard onsite repulsions for bosons of sorts $`a`$ and $`b`$. Finally $`U_{ab}`$ is an onsite attraction between bosons of two different sorts.

Let us consider the low-density limit, when both $`n_a1`$ and $`n_b1`$. In this limit we must replace the Hubbard interaction $`U_{ab}`$ by the corresponding $`T`$ matrix. The relevant expression for the $`T`$ matrix $`T_{ab}`$ is given by

$$T_{ab}(\stackrel{~}{E})=\frac{U_{ab}}{1U_{ab}\frac{d^2p}{4\pi ^2}\frac{1}{p^2/m+|\stackrel{~}{E}|}},$$
(46)
where $`\stackrel{~}{E}`$ is given again by $`\stackrel{~}{E}=E+W`$. The $`T`$ matrix has the pole for the energy

$$|\stackrel{~}{E}|=|E_b|=W\mathrm{exp}\left\{\frac{4\pi }{mU_{ab}}\right\}.$$
(47)
In the extremely strong coupling case $`U_{ab}>W`$ the pole corresponds to the energy $`|E|=|E_b|=U_{ab}`$. The pole in the $`T`$ matrix reflects the appearance of the bound state of the two bosons of different sorts.

Now we can solve the two-particle problem in the presence of the bosonic background. A simple analysis shows that only local bosonic pairs (bipolarons) are possible in our case. They are formed at a high crossover temperature $`T^{}=|E_b|/\mathrm{ln}(W/T_0)`$, where $`T_0=\mathrm{min}\{T_0^a,T_0^b\}`$, and $`T_0^a`$ and $`T_0^b`$ are degeneracy temperatures for bosons of the two sorts. Correspondingly the pairs are bose condensed at lower temperature $`T_c^{ab}=T_0/\mathrm{ln}(1/\lambda )`$, where $`\lambda =mU_{ab}/(4\pi )`$. Our results are valid is the case $`|E_b|>\{T_0^a,T_0^b\}`$. In the opposite case of higher densities, when at least one of the temperatures $`T_0^a`$ or $`T_0^b`$ is larger then $`|E_b|`$, we have at first a standard one particle condensation for bosons with higher density. As a result the two-particle pairing between bosons of different sort can take place only as a second superfluid transition inside the superfluid phase with one particle BEC. Let us now analyze the stability of our system with respect to collapse. For simplicity we consider an extremely strong coupling case $`\{U_{aa},U_{bb},U_{ab}\}>W`$. In this case the local pairs of two bosons of different sorts have onsite character $`a_ib_i0`$. To escape collapse in this case we must satisfy the following stability criteria: namely, we must prevent the collapse of $`ab`$ pairs on one site (a creation of the quartets $`aabb`$). This requires the inequality $`U_{aa}+U_{bb}4U_{ab}>0`$. These results can be important for slave-boson and Schwinger-boson theories.

## IX Slave-boson formulation of the $`tJ`$ model.
Application to HTSC.

We conclude our paper with a brief discussion of the problems which arise in a slave-boson formulation of the $`tJ`$ model. In this formulation close to half-filling ($`n1`$) an electron is represented as a product of holon and spinon :

$$c_{i\sigma }^{}=f_{i\sigma }^{}b_i.$$

A superconductive $`d`$-wave gap

$$\mathrm{\Delta }_d=c_ic_jc_ic_j$$

is a direct product

$$\mathrm{\Delta }_d=\mathrm{\Delta }_{sp}\mathrm{\Delta }_h^{}$$

of a spinon $`d`$-wave gap

$$\mathrm{\Delta }_{sp}=f_if_jf_if_j$$

and a holon $`s`$-wave gap

$$\mathrm{\Delta }_h=b_ib_j.$$

Then a natural question arises whether $`b_i0`$ and, accordingly, $`\mathrm{\Delta }_h=b_ib_j`$, or $`b_i=0`$ but $`b_ib_j0`$. In other words, whether a one particle or two-particle condensations of the holons takes place in our system .

This problem is a very difficult one and, surely deserves a very careful analysis. Our preliminary considerations show, however, that the more beneficial conditions for the two-particle condensation may arise in the SU(2) formulation of the $`tJ`$ model . In the standard U(1) formulation of the model  an effective potential of the two-holon interaction on neighboring sites appearing after the Hubbard-Stratonovich transformation has a form

$$\left(\frac{8t^2}{J}\frac{J}{4}\right)\underset{ij}{}b_i^{}b_j^{}b_ib_j,$$
(48)
and thus corresponds to the repulsion for $`t>J`$. This observation excludes the possibility of the two-holon pairing in the U(1) formulation of the $`tJ`$ model.

In the SU(2) case it will be desirable to derive conditions when $`b_1=b_2=0`$ but $`b_1b_20`$. For such a nondiagonal pairing, as already discussed above it is easier to satisfy the stability criteria . Note also that the same situation with two sorts of bosons and a possible attraction between them can be realized for 2D magnetic systems. We obtain the corresponding bosonic Hamiltonian here after a Schwinger transformation of spins  in extended Heisenberg models. Work along these lines is now in progress.

## X Conclusions

In conclusion we analyzed the possibility of the formation of boson pairs with $`s`$-wave symmetry and an appearance of total phase separation in a 2D Bose gas. In addition we considered the case of boson pairs with symmetry of $`p`$\- and $`d`$-wave type. We also constructed the qualitative phase diagram for the 2D Bose gas with the van der Waals interaction between the particles, which, besides a standard one particle BEC, contains regions of the Mott-Hubbard Bose solid and a total phase separation. We also consider the situation for two sorts of bosons described by the two-band Hubbard model, and found the conditions for the two-particle pairing between bosons of different sorts. We discuss the applicability of our results for the different physical systems ranging from submonolayers of 4He and excitons in the semiconductors till Schwinger bosons in magnetic systems and holons in the slave-boson theories of HTSC.

The authors acknowledge helpful discussions with P. Fulde, O. Fischer, Yu. Kagan, N.M. Plakida, G. Khaliullin, S.H.S. Salk, M.A. Baranov, D. Ivanov, G. Jackeli, G.V. Shlyapnikov and V.I. Marchenko. This work was also supported by Russian Foundation for Basic Research (grant 02-02-17520), Russian President Program for Science Support (grant 96-15-9694) and by Grant of Russian Academy of Sciences for young scientists.


text:
# Low temperature resistance minimum in non-superconducting 3R-Nb1+xS2 and 3R-GaxNbS2

## 1 Introduction

The layered transition metal dichalcogenides (LTMDs) MX2 of the Group V metals (M = V, Nb, Ta; X = S, Se) and their intercalation compounds have been the subject of numerous studies on the inter-relationship between superconductivity and charge density waves (CDW), both of which arise from the strong electron-phonon ($`e`$-$`ph`$) coupling within the layers . Parameters such as stoichiometry, polymorphism, disorder and intercalation have been extensively used to study the physical properties of these low dimensional compounds. There is however, no satisfactory explanation for their effect on the above transitions.

Amongst the binary compounds, all di-Selenides (V, Nb, Ta) and all polymorphs of TaS2 show CDW transitions, while the 2H and 4H Nb and Ta compounds are also superconducting. 2H-NbS2 and 1T-VS2 are unusual – In the former, any CDW is suppressed below the superconducting $`T_\mathrm{c}`$ 6.2 K due to $`e`$-$`e`$ interactions in the narrow unhybridised Nb $`d`$ bands. 1T-VS2 is structurally metastable, supposedly due to reduced covalency which destabilizes the layered structure in its stoichiometric composition. The effects of $`e`$-$`e`$ interactions are also significant in 1T-VSe2 in which the normal to incommensurate CDW transition temperature $`T_{}`$ increases from 110 K under pressure induced broadening of the narrow V $`d`$ bands. It also shows a small Curie-like contribution in $`\chi `$ from a small excess of V between the layers. 2H-NbSe2 ($`T_{}`$ 33 K), 2H-TaS(Se)2 ($`T_{}`$ 75(122) K) and 4Hb-TaS(Se)2 ($`T_{}`$ 22(75) K) also exhibit $`T_\mathrm{c}`$ 7.2 K, $``$ .6 K and $``$ 1 K respectively. High pressure lowers $`T_{}`$ and raises $`T_\mathrm{c}`$ towards its ‘undistorted’ values ( and references therein). Lattice-dynamics calculations suggest that in these materials, $`e`$-$`ph`$ renormalisation effects on phonon frequencies are crucial in causing lattice instability, as well as raising $`T_\mathrm{c}`$ . The above transitions are absent in di-Tellurides in which strong intermetallic bonding and metal atom clustering leads to more stable but distorted structures.

Stoichiometry, intercalation and disorder significantly affect the CDW and superconductivity in these compounds. In 1T-TaS2 a mere 50 ppm of isoelectronic Nb impurities were observed to destroy the long range phase coherence of the commensurate CDW and completely suppress the associated metal-insulator (MI) transition . On intercalation, the $`T_\mathrm{c}`$ in Ta compounds increases from $`<`$ 1 K upto 5 K while it is lowered for Nb . In self-intercalated 2H-Nb1.05Se2 the $`T_\mathrm{c}`$ was suppressed below 2.2 K , while 2H-NbSe2(EDA)1/4 showed no superconductivity but a resistance minimum at $``$ 25 K . Post-transition metal intercalated 2H-NbS2 also did not show superconductivity . ‘Kondolike’ resistance minima at 20 K along with a CDW were observed for 2H-Fe.05Nb(Ta)Se2 . LixNbS2 ($`0x.5`$) has a complex dependence of $`T_\mathrm{c}`$ on $`x`$ due to interplay of polymorphic changes (2H–3R) and electron transfer effects upon intercalation . In fact, comparable disorder induced resistance minima at low temperatures are also seen in structurally similar Graphite intercalation compounds, and explained using weak localisation and $`e`$-$`e`$ interactions .

A less studied polymorph of NbS2 is the 3R phase (space group $`R3m`$). While the inherently stoichiometric 2H-NbS2 forms at high temperatures ($`950`$C) under high S vapour pressure, any excess Nb (at low S vapour pressure) results in the metal-rich 3R phase . The stoichiometry limit for single phase 3R-Nb1+xS2 depends on the preparation temperature, and a minimum of $`x=.03`$ has been reported at 650. Samples prepared at high temperature and reported closer to stoichiometry are essentially a 2H–3R mixture with properties correspondingly in between. The 3R phase can thus be considered a self-intercalated phase with the excess Nb in octahedral interlayer vacancy sites . While no superconductivity has been observed in 3R-Nb1+xS2 down to 1.7 K, there is also no evidence of a CDW which could depress the same. Conduction in LTMDs is essentially intralayer. Hence the role of even small amounts of interlayer metal in suppressing superconductivity between different polymorphs needs to be investigated to understand better the effects of polymorphism, stoichiometry, disorder and dimensionality on the properties of these compounds.

We have prepared 2H-NbS2, 3R-Nb1+xS2 ($`x.07`$) and 3R-GaxNbS2 ($`x=.1,.25,.33`$). 2H-NbS2 has $`T_\mathrm{c}`$ $``$ 6.2 K, as expected. The metal-rich 3R-Nb1+xS2 shows no superconductivity; rather, a resistance minimum $``$ 20 K. Non-magnetic Ga-intercalated 3R phases exhibit similar low temperature resistance minima (20 K $`T_{\mathrm{min}}<`$ 60 K). Thermopower also shows anomalies around $`T_{\mathrm{min}}`$. We observe a clear correspondence of $`T_{\mathrm{min}}`$ and extent of localisation to the cation-to-anion ratio and the preparation temperature, i.e., to the extent of structural distortion from the ideal layered host.

## 2 Experimental Details

### 2.1 Preparation and characterisation

The compounds were prepared by heating stoichiometric mixtures of the pure elements (Ga 99.999%, Nb 99.7%, S 99.999% pure, CERAC) in evacuated quartz ampoules. 2H-NbS2 was obtained using 10% excess Sulfur. Chemical analysis by ICP-AES and EDX techniques gave Fe $`<`$ 240 ppm, Mn $`<`$ 15 ppm and Cr $`<`$ 10 ppm. Pressed polycrystalline pellets were sintered at different temperatures – 950C for the binary compound, and 850C (LT phase) & 1100C (HT phase) for Ga-intercalated compounds – and quenched to room temperature. Single crystal flakes were obtained by vapour transport (950 $``$ 900C) for 2H-NbS2 and 3R-Ga.1NbS2. Some Sulfur rejection was observed at high temperature in all cases, indicating a metal rich composition.

The compounds were characterised by room temperature powder XRD. We compare the XRD patterns of 2H-NbS2, 3R-Nb1+xS2 and 3R-Ga.1NbS2(LT) in Fig. 1. The 2H phase ($`a`$ = 3.32 Å, $`c`$ = 11.97 Å) had broad ($`10l`$) lines indicating well known disorder in the stacking of the weakly coupled antiparallel layers . The metal rich 3R phases on the other hand had sharp lines, the HT phases more ordered than the LT ones. The excess metal (Nb / Ga) occupying the octahedral interlayer vacancies  pins the layers together and prevents stacking faults. 3R-Nb1+xS2 ($`a`$ = 3.32 Å, $`c`$ = 17.88 Å) and 3R-GaxNbS2 ($`x`$ = .1, .25, .33(LT phase)) ($`a`$ = 3.335 Å, $`c`$ = 17.905 Å) have similar lattice parameters. However, the Ga-intercalated phases exhibit superlattice reflections proportional in number and intensity to Ga content and the preparation temperature, indicating progressive distortion. The .33 Ga (HT) phase ($`a`$ = 7.19 Å, $`c`$ = 17.30 Å) is structurally very different, with about 10% reduction in volume and a doubling of the $`a`$-axis. The superlattice lines order into a doubled $`a`$-axis and show much lower $`c/a1.603`$. Thus with increasing cation concentration and preparation temperature, the trigonal prismatic close-packed layers transform towards a distorted octahedral coordination of Nb by staggered Sulfur layers. Progressive metal clustering and vacancy formation within the Nb layers finally leads to a well ordered phase with a different structure. This is similar to cation-rich Nb2Se3, Nb3S4, Cu.33NbS2, etc. which have strong Nb-Nb bonds giving rise to zig-zag chains / clusters. These changes are also reflected in the electronic properties. The exact nature of the Nb atom clustering in our compounds would require a more detailed structural study.

### 2.2 Electronic properties

The electronic properties were studied by 4-probe d.c.-resistance (4.2–300 K), thermopower ($`S(T)`$) (14–300 K) and magnetic susceptibility ($`\chi (T)`$)(80–300 K) measurements. The d.c.-resistance of pellets as well as single crystal flakes was measured in van der Pauw geometry  using Ag-paste contacts. The $`S(T)`$ of pellets pressed between Cu stubs was measured in the differential mode. Absolute $`S`$ was determined by calibrating with Pb , and correcting for the Cu leads. The $`\chi (T)`$ was measured on compacted powders under a 9.7 kOe field in a vibrating sample magnetometer (VSM) . The results are summarised in Table 1.

The samples showing resistance minimum obviously do not follow Mattheissen’s rule for scattering from dilute impurities under the quasi-free electron approximation. Therefore, we consider the $`\rho _{\mathrm{min}}`$ as characteristic of the residual resistance, and in subsequent discussions the residual resistance ratio (RRR) is defined by $`\rho _{300}/\rho _{\mathrm{min}}`$; while, for the superconducting samples the RRR = $`\rho _{300}/\rho _{T_\mathrm{c}}`$, where $`\rho _{T_\mathrm{c}}`$ is the value of the resistivity at $`T_\mathrm{c}`$.

The resistivity results on single crystal flakes and polycrystalline pellets of both 2H and 3R polymorphs are shown in Fig. 2 and 3. Grain boundary scattering in the pellets and the anisotropy of conduction in flakes would prevent their inter-comparison. Figure 3 therefore shows the normalised resistance behaviour after subtracting low temperature residual resistances. The 2H phase $`T_\mathrm{c}`$ $`6.2`$ K and its large RRR of $`69`$ for in-plane conduction in the flakes compare well with earlier studies . Thus inspite of a considerable degree of stacking disorder, the essentially in-plane conduction and also the superconductivity is not affected in 2H polymorphs. The 3R phases on the other hand have no stacking disorder as inferred from their sharp XRD pattern, but gave a resistance minimum around and above 20 K. The temperature and extent of the minimum increases with the amount of intercalate atoms.

The overall conduction behaviour of the 3R phases showing a resistance minimum is also significantly different. Firstly, we notice a large increase in residual resistance upon intercalation (Table 1). For example, the .1 Ga crystal flake measured along the plane has a RRR of only 1.65, although its resistivity is actually much lower than superconducting 2H-NbS2 (RRR $``$ 69). The scattering of electrons at higher temperatures is also significantly different. We observe in the 3R phases a comparatively faster drop in resistance of flakes as well as pellets and also a larger contribution of T2 \- term above their resistance minimum (Fig. 3). The reduction in saturation effects in the high temperature conduction of the polycrystalline 3R-Nb1+xS2 phase probably indicates an increase in its isotropy of conduction. A proper explanation of this difference would require a detailed investigation of the effects of non-stoichiometry and polymorphic changes on the band structure and $`e`$-$`ph`$ interaction effects. The observed differences, however, cannot be simply related to doping of carriers since extra Nb or Ga atoms have similar effects. We presently conclude that in the intrinsically non-stoichiometric 3R-Nb1+xS2 phases the intralayer scattering potentials are significantly increased, leading to large residual resistance, increased $`e`$-$`e`$ scattering effects, and consequently low temperature localisation of carriers.

The $`S(T)`$ and $`\chi (T)`$ results on polycrystalline 2H and 3R phases are shown in Fig. 4. The $`S(T)`$ of 3R phases of self intercalated and .1 Ga intercalated phase are very similar but differ markedly from that of 2H-NbS2. For the former a large, fairly constant, negative $`S(T)`$ at high temperature changing rapidly towards positive below 100 K can be seen. For $`\chi (T)`$ we have subtracted a small saturated contribution which was observed in the $`M(H)`$ behaviour at low fields ($`<1`$ kOe). Our data differ from previous studies, as we observe a continuous reduction in $`\chi (T)`$ instead of a slight increase on cooling . These are the typical dependences observed in compounds showing CDW instabilities.

## 3 Discussion

We have earlier mentioned that group V TMD compounds of other than the present studies, i.e., layered 1T-VSe2, 2H-NbSe2 and various polymorphs of TaS(Se)2, show varying degrees of CDW formation. The associated anomalies in their transport and magnetic properties are much weaker than in 1D CDW structures such as NbS3, TaS3 etc. For example, except for 1T-TaS2 which has a CDW lead MI transition, the usual increase in resistance at $`T_{}`$, is barely observable and only a steep increase in its slope is observed below $`T_{}`$. Recent optical studies on 2H-TaSe2 confirm the absence of any abrupt formation of a charge-excitation gap at $`T_{}`$ . The steep increase in the resistance slope below $`T_{}`$ is found to be a consequence of freezing-out of scattering channels since the Drude scattering peak in $`\sigma (\omega )`$ ($`\omega `$ $``$ 0) becomes narrower below $`T_{}`$. In the light of the above observations, we cannot rule out CDW-correlations in our 3R phases, more so since superstructural distortions are clearly observed in intercalated phases. The thermopower and magnetic susceptibility variation also suggest some non-magnetic electronic correlations developing on cooling. A careful structural study of the 3R polymorphs of NbS2 at low temperature is therefore required to ascertain the relation of the superlattice with the observed transport and magnetic properties.

In Fig. 5 we show a scaled plot of resistance behaviour with temperature for the compounds of our study between $`.2<T/T_{\mathrm{min}}<2`$. There is seen to be a close similarity for the crystal flakes and pellets of compounds having different intercalate (Ga/Nb) concentration and showing vastly different resistivity values. It clearly indicates the role of defects in these compounds since $`\rho _{\mathrm{min}}`$ and $`T_{\mathrm{min}}`$ show a systematic increase with the intercalate concentration.

The resistance minimum may be sought to be explained by Kondo scattering from dilute magnetic impurities such as Fe or Mn . However, as stated earlier, the maximum concentration of such impurities in our compounds ($`<250`$ ppm) is too low for significant spin-flip scattering contribution. Neither do we observe the expected logarithmic dependence of resistivity upon temperature for $`T<T_{min}`$. The temperature of the Kondo minimum plotted as $`\rho /\rho _{\mathrm{min}}`$ vs $`T`$ is in fact only weakly dependent upon the extrinsic impurity concentration. Instead, we find a strong dependence of $`T_{\mathrm{min}}`$ on intercalate concentration which scales as shown in Fig. 5. Moreover, our $`S(T)`$ is smoothly varying in the temperature range of resistance minimum instead of the broad, large maxima expected in Kondo systems. The Pauli like $`\chi (T)`$ of our compounds at high temperature is also in contrast to the Curie-Weiss like behaviour of Kondo systems for $`TT_{\mathrm{min}}`$. Thus Kondo scattering cannot be invoked to explain the observed behaviour of our compounds.

On the other hand, disorder and Coulomb interaction effects on the scattering of electrons in a narrow anisotropic band of these compounds are expected to be quite important. For $`k_\mathrm{F}l1`$, ($`k_\mathrm{F}`$ – Fermi momentum, $`l`$ – elastic mean free path), the quasi-classical treatment of elastic scattering leading to Mattheissen’s rule in dilute alloys breaks down. The rise in resistance on cooling can be understood in terms of either increasing interference corrections for the elastically scattered waves from the static disorder or subtle changes in the excitation spectrum at the Fermi level caused by $`e`$-$`e`$ interaction effects in the presence of ionic disorder, as first suggested by Altshuler and Aronov . The relative importance of these quantum correlations and their temperature dependence depends on dimensionality and the details of band structure of the conduction electrons.

The plots of $`\mathrm{\Delta }\sigma /\sigma _0`$ at low temperature are shown in Fig. 6, where $`\sigma _0`$ is the value of conductivity at 0 K obtained by extrapolating below 4.2 K. Our preliminary results show a nearly $`T^{1/2}`$ dependence of the conductivity far below its maximum, the coefficient of the $`T^{1/2}`$ term depending upon the stoichiometry and structural details. The observed behaviour suggests quantum corrections to the DOS due to long range Coulomb interactions between conduction electrons . The metal rich .33 Ga (HT) phase also exhibits similar behaviour, though with a much larger slope. Here the localisation effects are dramatically increased due to the clustered nature of the metallic lattice.

We mention here an interesting possibility which may increase the long range Coulomb interaction effects between charge carriers. The theory of Altshuler and Aronov gives a depression in DOS($`E_\mathrm{F}`$) . The present compounds are prone to CDW correlations which also reduce DOS($`E_\mathrm{F}`$). However, for CDW the $`e`$-$`e`$ interactions are very special since they require $`e`$-$`ph`$ coupling across a nesting wave vector. Therefore, the role of incipient CDW fluctuations in the presence of disorder to give the observed resistance minimum behaviour in these compounds should be seriously explored.

To summarise, we have observed low temperature resistance minima in non-superconducting metal-rich NbS2 derivatives – non-stoichiometric 3R-Nb1+xS2 ($`T_{\mathrm{min}}`$$``$ 20 K) as well as 3R-GaxNbS2 (20 K $`T_{\mathrm{min}}<`$ 60 K). A common physical origin of these minima is evident in the scaling of resistance as $`\rho /\rho _{\mathrm{min}}(T/T_{\mathrm{min}})`$ between $`.2<T/T_{\mathrm{min}}<2`$ for different phases whose resistivity differs by an order of magnitude. The low temperature behaviour of these non-magnetic compounds cannot be explained by Kondo scattering effects. Instead, we find that the conductivity varies as $`T^{1/2}`$ below its maximum (i.e., below $`T_{\mathrm{min}}`$($`\rho _{\mathrm{min}}`$)). We, therefore, propose that a possible cause for the observed behaviour in these narrow band anisotropic systems is a correction to the DOS($`E_\mathrm{F}`$) due to $`e`$-$`e`$ interaction effects in the presence of ionic disorder.

Acknowledgements: We are grateful to Prof. Deepak Kumar for illuminating discussions. AN thanks CSIR, New Delhi for financial support during his doctoral work.


text:
# 1 Introduction

## 1 Introduction

There are extensions of the Standard Model in which Higgs bosons have suppressed couplings into $`b`$-quarks. This can occur for specific parameters of the Two Higgs Doublet Model , or of the Minimal SuperSymmetric Model , as well as for some composite models . Standard Model Higgs searches  would have a reduced sensitivity in such cases, because of their strong reliance on the identification of the $`b`$-quarks from the Higgs boson decay to maximize the separating power. It is important to also cover such scenarios experimentally with dedicated searches in which the information from the flavour of the quarks in the Higgs boson decay is not exploited, so that the model-dependance of our final Higgs results can be reduced.

All four LEP collaborations have pursued such flavour independent searches in recent years , analysing the four-jet ($`q\overline{q}q\overline{q}`$), missing energy ($`q\overline{q}\nu \overline{\nu }`$) and leptonic ($`q\overline{q}l^+l^{}`$) topologies. None has found evidence for any signal. In this note, the first combination of the results obtained is presented, for the $`e^+e^{}Z^{}hZ`$ production mechanism, in terms of upper limits on the corresponding cross-section as a function of the Higgs boson mass, and of a lower mass limit in the assumption of production cross-sections equal to those in the Standard Model and a Higgs boson decaying into hadrons with a 100 % branching fraction. Higgs mass assumptions from 60 to 115 GeV were tested.

Some systematic effects have been studied by individual experiments, but were not considered nor included in the evaluation presented here. All results are preliminary.

## 2 Experimental Analyses

The analyses used by the four collaborations to search for $`hZ`$ production in the flavour independent hypothesis were to a large extent adaptations of existing LEP-2 searches or measurements. In all channels except four-jets, events were selected exactly, or almost exactly, as in the corresponding Standard Model Higgs search, removing the $`b`$-tagging from the final selection. In the four-jet channel, on the other hand, all collaborations have used dedicated test-mass dependent selections, to exploit maximally the kinematic features available and the mass reconstruction. This was necessary to enable reducing as much as possible the dominant backgrounds from $`WW`$ and $`ZZ`$ production and from QCD processes $`qqgg`$, in the absence of $`b`$-tagging.

Because the searches performed by the four collaborations have been developed relatively, or in some cases very, recently, they are less optimised, slightly less comprehensive, and have used more simplifying assumptions to extract the results than the corresponding searches for the Standard Model Higgs boson.

In the ALEPH search , the data collected in 1998, 1999 and 2000 were used, clustering the data from the last year in seven energy bins. All the topologies corresponding to possible decay products of the $`Z`$ boson were investigated. In the four-jet channel, a dedicated search based on a neural network method was developed. The compatibility of the data with the signal hypothesis was tested in the Higgs boson mass range from 60 up to 115 GeV/c2, with 5 GeV/c2 steps, interpolating the final discriminant variable to test intermediate values.

In the DELPHI search , only the data collected in 1999 and 2000 were used, clustering the data from the last year in two energy bins. All the topologies corresponding to possible decay products of the $`Z`$ boson were investigated, except the $`q\overline{q}\tau ^+\tau ^{}`$ final state. The compatibility of the data with the signal hypothesis was tested in the Higgs boson mass range from 50 up to 110 GeV/c2, with 5 GeV/c2 steps. A small degradation of the performance was introduced in the final evaluation, between each test mass, to account for mass resolution effects.

In the L3 search , only the data collected in 1999 and 2000 were used, clustering the data from the last year in five energy bins. All the topologies corresponding to possible decay products of the $`Z`$ boson were investigated. The compatibility of the data with the signal hypothesis was tested in the Higgs boson mass range from 60 up to 115 GeV/c2, with 1 GeV/c2 steps.

In the OPAL search , the data collected in 1998, 1999 and 2000 were used, clustering the data from the last year in a single energy bin. Dedicated analyses were conducted for all the topologies corresponding to possible decay products of the $`Z`$ boson. For the four-jet channel the test-masses were chosen in the Higgs boson mass range from 60 up to 115 GeV/c2, with 1 GeV/c2 steps.

Finally, in spite of not taking advantage of the $`b`$-tagging of the jets from the Higgs boson decay, every experiment has found small but still significant differences in performance between the different possible decay products of the Higgs boson, arising from slight differences in mass resolution and jet structure. In particular, Higgs boson decays into gluon pairs have larger multiplicities, but at the same time coarser dijet mass resolution, than decays into light quarks. In order to enable quoting genuine flavour independent results, the samples used for the final evaluation were conservatively chosen as those which gave the weakest expected performance in each channel, and for each value of the Higgs boson mass.

## 3 Results

The common evaluation of the results used the standard statistical procedures based on the likelihood ratio technique, as applied in the other combinations performed by the LEP working group for Higgs boson searches .

Results were first obtained for each of the four collaborations, with two independent implementations of the combination software , and compared with the results obtained within each collaboration. Although some slight differences were found, the general shapes and features of the observed and expected confidence levels for the signal, $`CL_s`$ and background-only, $`CL_b`$ hypotheses were very consistent.

The calculations of the observed and expected 95% CL lower limits on the mass of the Higgs boson, assuming production cross-sections equal to those in the Standard Model and BR(h$``$ hadrons) = 1.0, all agreed within a few hundred MeV, or better. The values111Because the evaluations performed did not include systematic uncertainties, the values obtained differ in some cases from those quoted by the individual colloborations which included them in their evaluations. This is most notably the case for the OPAL results . from one of the implementations  are shown in Table 1, together with the combined results using the data from all four collaborations. The combined observed and median expected limits were 112.9 and 113.0 GeV/c2, respectively.

The confidence levels $`CL_s`$ and $`CL_b`$ obtained from the full combination in the signal and background-only hypotheses are shown as a function of the mass in Figure 1.a and 1.b, respectively. Good overall agreement can be seen between the observation and expectation in the absence of a signal. The slightly depressed values for $`1CL_b`$ for masses below 80 GeV/c2, witnessing excesses of data in this region, may be the results of some statistical or systematic effect, and should be investigated further. A 5 sigma discovery corresponds to a value of 5.7 $`10^7`$, as indicated by the horizontal line. The maximal sensitivity for such a discovery is reached for an assumed Higgs boson mass of 107 GeV/c2, when this line is intersected by the expected median confidence in the background-only hypothesis.

Upper limits on the production cross-section as a function of mass were also determined for each individual collaboration, and combined using the data from all four collaborations. Typically, cross-sections larger than about 10-60 % of the expected Standard Model value were excluded at 95% CL in the mass range 60-100 GeV/c2, by each collaboration alone. The combined exclusion from the four collaborations is shown as a function of the mass in Figure 2. In the same mass range, it is possible to exclude cross-sections larger than about a few to 30 % of the expected Standard Model value with the full LEP-2 data set. Good overall agreement between the observation and expectation can be seen also here, except in the lower mass range, where the weaker observed limit is resulting from some excesses of data in this region, as was already noted.

## Acknowledgements

We are greatly indebted to our technical collaborators, to the members of the CERN-SL Division for the excellent performance of the LEP collider, and to the funding agencies for their support in building and operating the ALEPH, DELPHI, L3 and OPAL detectors. We also acknowledge the fruitful collaboration between all four collaborations, and the open spirit which prevailed and enabled to implement efficiently the established combination procedures to the new channels described in this note.


text:
# Defect Turbulence in Inclined Layer Convection

## Abstract

We report experimental results on the defect turbulent state of undulation chaos in inclined layer convection of a fluid with Prandtl number $`1`$. By measuring defect density and undulation wavenumber, we find that the onset of undulation chaos coincides with the theoretically predicted onset for stable, stationary undulations. At stronger driving, we observe a competition between ordered undulations and undulation chaos, suggesting bistability between a fixed-point attractor and spatiotemporal chaos. In the defect turbulent regime, we measured the defect creation, annihilation, entering, leaving, and rates. We show that entering and leaving rates through boundaries must be considered in order to describe the observed statistics. We derive a universal probability distribution function which agrees with the experimental findings.

Weakly driven, dissipative pattern-forming systems often exhibit the spatiotemporally chaotic state of defect turbulence, where the dynamics of a pattern is dominated by the perpetual nucleation, motion, and annihilation of point defects (or dislocations) . Examples can be found within striped patterns in wind driven sand, electroconvection in liquid crystals , nonlinear optics , fluid convection , autocatalytic chemical reactions , and Langmuir circulation in the oceans . The hope is that the dynamics of these very different systems can be characterized by a universal description which is based only on the defect dynamics.

A first description of defect turbulence was given by Gil et al.  for a spatiotemporally chaotic state of the complex Ginzburg-Landau equation (CGLE). They postulated that the nucleation rate for defect pairs is independent of the number of pairs $`M`$, and based on the topological nature of defects the annihilation rate is proportional to $`M^2`$. Through detailed balance, they showed that these assumptions lead to a squared Poisson distribution for the probability distribution function (PDF) of $`M`$. They also found agreement with this PDF in simulations of the CGLE with periodic boundary conditions. Rehberg et al.  measured the PDF of $`M`$ for defect turbulence in electroconvection of nematic liquid crystals and found it to be consistent with the predicted squared Poisson distribution. Later, Ramazza et al. investigated a defect turbulent state in optical patterns and found that their data was not conclusive. To date, studies in both simulation and experiment have relied purely on comparisons of the PDFs. The gain and loss rates, fundamental to the universal description of defect turbulence, have not been measured. In addition, effects due to boundaries were not considered.

In this Letter, we report experimental results on the defect turbulent state of undulation chaos in inclined layer convection of a fluid of Prandtl number $`1`$. By tracking all defects in a finite area of the convection cell we measured, for the first time, defect creation, annihilation, leaving, and entering rates for a defect turbulent state. The observed pair creation and annihilation rates agree with the predictions . In an experimental system where periodic boundary conditions do not apply, pairs of defects are no longer an appropriate description. To describe the statistics correctly, single defect leaving and entering rates through the boundaries must be considered. Our data show that the entering rates are approximately independent of the number of positive/negative defects $`N_\pm `$ and the leaving rates are proportional to $`N_\pm `$. We derive a universal PDF which reduces to the earlier predicted squared Poisson distribution  when boundary effects are negligible. Our measurements agree with this new distribution. In addition, by measuring both $`N_\pm `$ and the undulation wavenumber we determined that the onset of undulation chaos coincides with the theoretically predicted onset for stable, stationary undulations. At higher driving, we observe competition between regions of ordered undulations and undulation chaos, suggesting bistability between a fixed point attractor and spatiotemporal chaos.

When an inclined fluid layer is heated from below and cooled from above for small temperature differences $`\mathrm{\Delta }T`$ the fluid layer experiences not only a linear temperature gradient (as in Rayleigh-Bénard convection) but also a symmetry-breaking large scale shear flow, with fluid rising at the warmer and falling at the cooler plate . Above a critical temperature difference $`\mathrm{\Delta }T_c`$, a pattern of longitudinal rolls aligned with the uphill/downhill direction is observed for intermediate inclination angles. For $`ϵ0.02`$, where $`ϵ\mathrm{\Delta }T/\mathrm{\Delta }T_c1`$, these longitudinal rolls are unstable to undulations . While it was predicted theoretically that a perfect pattern of undulations is stable , a defect turbulent state of undulation chaos was found experimentally . Recently, undulation chaos was also observed in numerical simulations of the Navier-Stokes equations when random initial conditions were used , while with controlled initial conditions undulations were found to be stable. As shown in the snapshot in Fig. 1a, undulation chaos is characterized by a pattern of disordered undulating rolls, punctuated by point defects which carry a topological charge of $`\pm 2\pi `$ . Defects are nucleated pairwise in regions of low convective amplitude (green), enter and leave the system at the boundaries, preferentially move at right angles to the rolls, and annihilate as pairs. For $`ϵ0.10`$, the defect turbulent state appears intermittently with ordered regions of undulations, suggesting a competition between the fixed point attractor (undulations) and spatiotemporal chaos (undulation chaos). This is shown by the two snapshots in Fig. 1b, c. The dynamics of the state is best viewed in the movies available at Ref. .

The experiments were conducted in the same apparatus as Ref.  for an inclination angle of $`\theta =30^{}`$. The fluid was compressed CO2 at a pressure of $`(56.5\pm 0.01)`$ bar regulated to $`\pm 0.005`$ bar with a mean temperature of $`(28\pm 0.05)^{}`$C regulated to $`\pm 0.0003^{}`$C. For these parameters, convection appeared at $`\mathrm{\Delta }T_\mathrm{c}=(1.763\pm 0.005)^{}`$C. The planform of the convection pattern was observed by the usual shadowgraph technique . A cell of height $`d=(388\pm 2)\mu `$m and dimensions $`101d\times 50d`$ was used, for which the vertical diffusion time was $`\tau _v1.3`$ sec. All data was measured in a homogeneous central subregion of size $`32d\times 25d`$ . To track defects, runs of 100 shadowgraph images were taken at a rate of 3 frames/sec ($`\tau _v/4`$) with a 10 bit 1K $`\times `$ 1K digital CCD camera. For 17 values of $`ϵ`$ between 0.04 and 0.22 this was repeated for 600 (500, 400) runs at $`ϵ0.07`$ ($`0.08ϵ0.10`$, $`ϵ0.12`$). Runs at the same $`ϵ`$ were separated by at least $`100\tau _v`$ for statistical independence. A total of 1.5 terabytes of data were analyzed. Each value of $`ϵ`$ was reached via a quasistatic temperature increase in steps of $`0.001^{}`$C, waiting $`700\tau _v`$ between steps for transients to die out. We also conducted a sequence of measurements with quasistatic temperature decreases to check for possible hysteresis, which was not observed.

To lowest order, the planform of undulation chaos can be captured by the real part of the two dimensional field $`\mathrm{\Psi }(x,y)`$ with

$$\mathrm{\Psi }=\psi _0e^{i(qx+\varphi )}+i\psi _1\left(e^{i(qx+py)}+e^{i(qxpy)}\right)+\mathrm{}$$
(1)
where the wavenumbers $`q`$ and $`p`$ describe the stripe spacing and waviness, respectively. Both the wavenumbers and the amplitudes $`\psi _0`$ and $`\psi _1`$ vary slowly in space and time. For a perfect undulation pattern $`\varphi `$ is constant everywhere, while for undulation chaos $`\varphi (x,y,t)`$ approximates the deformation of the pattern due to the defects. The ansatz $`\varphi (x,y,t)=c_i\mathrm{arctan}(\frac{yy_i}{xx_i})`$ describes the observed pattern well for constant $`q,p,`$ and $`\psi _1/\psi _0`$, where $`(x_i,y_i)`$ are the locations of defects and $`c_i=\pm 1`$ are their corresponding topological charges.

In Fig. 2 the peak values for the distributions of $`p`$, $`q`$, and $`\psi _1/\psi _0`$ are shown as a function of $`ϵ`$ . The values for $`q`$ change little with $`ϵ`$, while the undulation wavenumbers can be fit by $`p\sqrt{ϵϵ_u}`$, shown as a solid black line. This fit provides a value for the onset of undulation chaos at $`ϵ_u=0.017\pm 0.001`$, which is in good agreement with $`ϵ_u=0.016`$ predicted theoretically for the onset of stable stationary undulations .

In order to better characterize the generic properties of defect turbulence, we determined defect locations and the corresponding topological charges in each shadowgraph image. A time trace of the number of defects over $`2\times 10^5\tau _v`$ at $`ϵ=0.13`$ ($`10^4`$ defects) revealed that the undulation chaos state is statistically stationary and neutral with respect to topological charge. Nonetheless, any single realization of the pattern does not necessarily satisfy $`N_+=N_{}`$ due to the boundaries. From the list of defect locations in each frame, we determined the mean $`N_\pm `$ and the variance $`\sigma _\pm ^2=N_\pm ^2N_\pm ^2`$ in the chosen subregion, plotted in Fig. 3 as a function of $`ϵ`$. Here and below, we report results only for $`N_+`$; equivalent results were obtained for $`N_{}`$.

For $`ϵ0.10`$, $`N_\pm (ϵ)`$ is linearly increasing. Extrapolation to $`N_+=0`$ gives $`ϵ_d=0.025\pm 0.002`$, which is close to the onset of undulations ($`ϵ_u=0.017\pm 0.001`$) determined above from wavenumber measurements. The difference may be attributable to the finite size of the system. For $`ϵ>0.10`$, $`N_\pm (ϵ)`$ decreases in the bistability regime, while no such decrease is seen in $`p`$.

Two predicted generic features of defect turbulence  are the constant creation rate and the annihilation rate quadratic in the number of defect pairs $`M`$. By connecting defect locations in adjacent frames into defect tracks, and then associating tracks of opposite sign to locate creation and annihilation events, we were able to perform the first direct test of the theoretical postulates, as well as the entering and leaving rates to/from the observation area. For a system where single defects can enter and leave through the boundaries a pair-based description is no longer justified. Thus, we separately consider the PDFs for $`N_\pm `$.

The creation rate $`C`$ was observed to be approximately independent of $`N_\pm `$ and the annihilation rate agreed with $`AN_+N_{}`$, by extension from $`AM^2`$ for the topologically neutral case. Fig. 4 shows the four gain/loss rates, tabulated using positive defects and events only. To lowest order, the entering rate is independent of the number of defects already in the subregion of the cell, suggesting that both creation and entering are random processes. The leaving rate was found to increase in proportion to $`N_\pm `$ as expected, since unlike annihilations such events do not depend on the presence of two oppositely-charged defects.

The observed rates are approximately given by

$`E(N_\pm )`$ $`=`$ $`E_0`$ (2)
$`C(N_\pm )`$ $`=`$ $`C_0`$ (3)
$`L(N_\pm )`$ $`=`$ $`L_0N_\pm `$ (4)
$`A(N_\pm )`$ $`=`$ $`A_0N_\pm ^2`$ (5)
which consider defects of a single sign only. The fits corresponding to Eqn. 4 are shown in Fig. 4, for three values of $`ϵ`$. It is important to note that this is a simplification. $`C(N_\pm )`$ and $`E(N_\pm )`$ show a weak, negative slope: defects are slightly less likely to enter or be created if there are already defects present. We also observe a weak dependence of these rates on the number of oppositely-charged defects: $`E(N_+)`$ is slightly diminished and $`L(N_+)`$ is slightly elevated for $`N_+>N_{}`$. Thus, to fully describe the system, we would need to consider rates and distributions for both $`N_+`$ and $`N_{}`$ simultaneously. In order to make analytical progress we will assume below that Eqn. 4 are a complete description. As will be shown below, these assumptions approximate the data well.

Using Eqn. 4, we can construct the PDF for $`NN_\pm `$. Assuming a stationary distribution, detailed balance requires $`\mathrm{l}oss(N)𝒫(N)=\mathrm{g}ain(N1)𝒫(N1)`$ to describe the probabilities at adjacent $`N`$ . For the relevant rates,

$`𝒫(N)`$ $`=`$ $`{\displaystyle \frac{E(N)+C(N)}{L(N)+A(N)}}𝒫(N1)`$ (6)
$`=`$ $`{\displaystyle \frac{\alpha }{\beta N+N^2}}𝒫(N1).`$ (7)
where $`\alpha \frac{E_0+C_0}{A_0}`$ and $`\beta \frac{L_0}{A_0}`$. By performing the recursion and properly normalizing the distribution, we find a modified Poisson distribution

$$𝒫_{\alpha ,\beta }(N)=\frac{\alpha ^{\frac{\beta }{2}+N}}{I_\beta (2\sqrt{\alpha })\mathrm{\Gamma }(1+\beta +N)N!}$$
(8)
where $`I_\beta `$ is the modified Bessel function. For $`\beta =0`$ ($`\alpha =N^2`$) this PDF reduces to the squared Poisson distribution of Ref. . In our experiments, $`\beta 3`$ .

Fig. 5b shows that the squared Poisson and Poisson distributions match the data poorly. Both the modified Poisson distribution and a distribution obtained by performing the recursion of Eqn. 6 on the raw experimental data from Fig. 4 agree well with the experimental data. The modified Poisson distribution also correctly captures the mean and width of the PDF. Using the already normalized modified Poisson distribution (Eqn. 8), we calculate $`N=N𝒫_{\alpha ,\beta }(N)`$ and $`N^2=N^2𝒫_{\alpha ,\beta }(N)`$ and obtain the analytic expressions

$`N`$ $`=`$ $`\sqrt{\alpha }{\displaystyle \frac{I_{\beta +1}(2\sqrt{\alpha })}{I_\beta (2\sqrt{\alpha })}}`$ (9)
$`N^2`$ $`=`$ $`\sqrt{\alpha }{\displaystyle \frac{I_{\beta +1}(2\sqrt{\alpha })+\sqrt{\alpha }I_{\beta +2}(2\sqrt{\alpha })}{I_\beta (2\sqrt{\alpha })}}.`$ (10)
Again, the above equations reduce to the corresponding equations for the squared Poisson distribution in the limit $`\beta =0`$. In Fig. 5a the dependence of the mean on the variance is plotted for the raw data and the modified and squared Poisson distributions. The modified Poisson distribution (Eqn. 5) shows excellent agreement, while the squared Poisson distribution clearly fails to describe the experimental data.

In summary, we find that the state of undulation chaos, and not the fixed point attractor to ordered undulations, is selected above the predicted onset for undulations. At higher driving, we find an apparent competition between the two attractors.

Our result for the $`N`$-dependence of the defect creation, annihilation, leaving, and entering rates and the modified Poisson (Eqn. 8) should be generic to any defect turbulent system. We have found agreement with the theoretical predictions for creation and annihilation rates  and have extended the analysis to the experimentally relevant case of a finite system without periodic boundary conditions where it is important to consider both positive and negative defects (rather than pairs). Although we observed weak fluctuations away from topological charge neutrality, a description based purely on a single topological charge was sufficient to describe the observed behavior .

We thank W. Pesch and O. Brausch for sharing the results of simulations and stability analysis which have enhanced our understanding of the bistability. We would like to thank the NSF for support under DMR-0072077.