text:
# 1 𝑝_𝑇 distribution of double spin asymmetries at √𝑠=50 GeV for Λ_𝑐⁺ production (left panel) and for Λ̄_𝑐⁺ production (right panel).

Charmed Hadron Production
in Neutrino Reactions
and Polarized $`s`$–quark Distribution111Talk presented by K. Sudoh at the NuFACT’01 workshop, Tsukuba, Japan, May 24–30 (2001).

Kazutaka SUDOH

Graduate School of Science and Technology,
Kobe University, Nada, Kobe 657–8501, Japan

E-mail:sudou@radix.h.kobe-u.ac.jp

Toshiyuki MORII

Faculty of Human Development,
Kobe University, Nada, Kobe 657–8501, Japan

E-mail:morii@kobe-u.ac.jp

Abstract

In order to obtain information about the polarized strange quark distribution, we studied the semi–inclusive $`\mathrm{\Lambda }_c^+`$/$`\overline{\mathrm{\Lambda }}_c^+`$ production in neutrino and polarized proton collisions. We found that these reactions are effective to clearly extract the polarization of strange quark by measuring the spin correlation between the target proton and the produced $`\mathrm{\Lambda }_c^+`$/$`\overline{\mathrm{\Lambda }}_c^+`$ baryon.

PACS numbers: 13.15.g, 13.88.e, 14.20.Lq

Proton spin puzzle is currently one of the most challenging topics in high energy spin physics. As is well known, proton spin is composed of the spin and orbital angular momentum of quarks and gluons which constitute proton. The polarized parton distribution plays an important role on the spin structure of proton. However, knowledge of the polarized sea quark distribution remain still poor. In order to understand the spin structure of proton, we need more information about the polarized sea quark distribution functions. Recently, HERMES group has reported  that direct measurement of the strange sea is required to explain the violation of the Ellis–Jaffe sum rule . So far, there are several parametrization models of the strange quark distribution. Though the most simplest case is to assume the flavor $`\mathrm{SU}(3)_f`$, a new parametrization taken account of the violation of $`\mathrm{SU}(3)_f`$ is also recently proposed .

Here we focus on the polarized $`s`$/$`\overline{s}`$ quark distribution. To examine the polarized $`s`$/$`\overline{s}`$ quark distribution, we have studied semi–inclusive $`\mathrm{\Lambda }_c^+`$/$`\overline{\mathrm{\Lambda }}_c^+`$ neutrino production; $`\nu +\stackrel{}{p}l^{}+\stackrel{}{\mathrm{\Lambda }}_c^++X`$, $`\overline{\nu }+\stackrel{}{p}l^++\stackrel{}{\overline{\mathrm{\Lambda }}}_c^__++X`$, which might be observed at the planned factory, where arrows attached to particles mean that these particles are polarized. Since in the naive quark model, the $`\mathrm{\Lambda }_c^+`$ baryon is composed of a heavy $`c`$–quark and anti–symmetrically combined light $`u`$ and $`d`$–quarks, we can assume the polarization of $`\mathrm{\Lambda }_c^+`$ baryon to be the one of $`c`$–quark. In addition, $`\mathrm{\Lambda }_c^+`$ is dominantly produced by the fragmentation of $`c`$–quark which is originated from $`s`$–quark through the $`t`$–channel $`W`$ exchange at the leading order. Therefore, there can be a correlation between the $`s`$–quark polarization and the produced $`c`$–quark polarization. Hence we can expect that the measurement of the spin correlation between the incident proton and the produced $`\mathrm{\Lambda }_c^+`$/$`\overline{\mathrm{\Lambda }}_c^+`$ gives us information about the polarized $`s`$/$`\overline{s}`$–quark distribution in proton.

For above processes, we calculated the double spin asymmetry $`A_{LL}^H`$ for final state hadron specified by $`H`$ ($`\mathrm{\Lambda }_c^+`$ or $`\overline{\mathrm{\Lambda }}_c^+`$), which is given by

$$A_{LL}^H=\frac{d\mathrm{\Delta }\sigma /dp_T}{d\sigma /dp_T},$$
(1)
where $`p_T`$ is a transverse momentum of final hadron $`H`$. $`d\mathrm{\Delta }\sigma `$ represents the spin–dependent differential cross section and is defined in terms of $`d\sigma _{hh^{}}`$ as

$$d\mathrm{\Delta }\sigma \frac{1}{4}[d\sigma _{++}d\sigma _++d\sigma _{}d\sigma _+],$$
(2)
with definite helicities $`h`$ and $`h^{}`$ for incoming proton and outgoing $`\mathrm{\Lambda }_c^+`$/$`\overline{\mathrm{\Lambda }}_c^+`$, respectively.

For $`\mathrm{\Lambda }_c^+`$ production, the spin–dependent differential cross section is given

$$d\mathrm{\Delta }\sigma (\nu pl^{}\mathrm{\Lambda }_c^+X)=\left\{\mathrm{U}_{cs}^2\mathrm{\Delta }s(x)+\mathrm{U}_{cd}^2\mathrm{\Delta }d(x)\right\}dx\left(\frac{d\mathrm{\Delta }\widehat{\sigma }}{d\widehat{t}}\right)d\widehat{t}\mathrm{\Delta }D_c^{\mathrm{\Lambda }_c^+}(z)dz,$$
(3)
where $`\mathrm{\Delta }s(x)`$ and $`\mathrm{\Delta }d(x)`$ are the polarized $`s`$–quark and $`d`$–quark distribution functions, respectively. $`\mathrm{U}_{cs}`$ and $`\mathrm{U}_{cd}`$ are CKM parameters. $`\mathrm{\Delta }D_c^{\mathrm{\Lambda }_c^+}(z)`$ is the polarized fragmentation function of outgoing charm quark decaying into $`\mathrm{\Lambda }_c^+`$. We used the model of Peterson et al.  as the unpolarized fragmentation function, which is given by

$$D_c^{\mathrm{\Lambda }_c^+}(z)=\frac{1}{z[1\frac{1}{z}\frac{ϵ_p}{1z}]^2}ϵ_p0.25\text{for}\mathrm{\Lambda }_c^+.$$
(4)
Unfortunately, the polarized fragmentation function is not yet established because of lack of experimental data. By analogy with the study on $`\mathrm{\Lambda }`$ polarization , we took the following ansatz:

$$\mathrm{\Delta }D_c^{\mathrm{\Lambda }_c^+}(z)=C_c^{\mathrm{\Lambda }_c^+}(z)D_c^{\mathrm{\Lambda }_c^+}(z),$$
(5)
where $`C_c^{\mathrm{\Lambda }_c^+}(z)`$ is the scale–independent spin transfer coefficient. Here we apply the analysis on $`\mathrm{\Lambda }`$ production to $`\mathrm{\Lambda }_c^+`$ production and choose the following two models; $`C_c^{\mathrm{\Lambda }_c^+}(z)=1`$ (the naive nonrelativistic quark model) and $`C_c^{\mathrm{\Lambda }_c^+}(z)=z`$ (the jet fragmentation models ). The double spin asymmetry $`A_{LL}^{\mathrm{\Lambda }_c^+}`$ is described in terms of the spin transfer coefficient $`C_c^{\mathrm{\Lambda }_c^+}(z)`$ and the ratio of the parton distribution functions as

$$A_{LL}^{\mathrm{\Lambda }_c^+}C_c^{\mathrm{\Lambda }_c^+}(z)\frac{\left\{\mathrm{U}_{cs}^2\mathrm{\Delta }s(x)+\mathrm{U}_{cd}^2\mathrm{\Delta }d(x)\right\}}{\left\{\mathrm{U}_{cs}^2s(x)+\mathrm{U}_{cd}^2d(x)\right\}}.$$
(6)
Thus $`A_{LL}^{\mathrm{\Lambda }_c^+}`$ is proportional to linear combination of the polarized $`s`$–quark and $`d`$–quark distribution function. Therefore, $`\mathrm{\Lambda }_c^+`$ production is not so good process for clearly extracting the $`s`$–quark distribution, since the contribution from the valence $`d`$–quark is large.

On the contrary, $`A_{LL}^{\overline{\mathrm{\Lambda }}_c^+}`$ in $`\overline{\mathrm{\Lambda }}_c^+`$ production can be represented as

$$A_{LL}^{\overline{\mathrm{\Lambda }}_c^+}C_c^{\mathrm{\Lambda }_c^+}(z)\frac{\left\{\mathrm{U}_{cs}^2\mathrm{\Delta }\overline{s}(x)+\mathrm{U}_{cd}^2\mathrm{\Delta }\overline{d}(x)\right\}}{\left\{\mathrm{U}_{cs}^2\overline{s}(x)+\mathrm{U}_{cd}^2\overline{d}(x)\right\}}=C_c^{\mathrm{\Lambda }_c^+}(z)\frac{\mathrm{\Delta }\overline{s}(x)}{\overline{s}(x)}=C_c^{\mathrm{\Lambda }_c^+}(z)\frac{\mathrm{\Delta }s(x)}{s(x)},$$
(7)
in which the $`\overline{d}`$–quark contribution can be eliminated and thus $`A_{LL}^{\overline{\mathrm{\Lambda }}_c^+}`$ is directly proportional to the strange quark distribution function. Therefore, we can clearly extract the polarized strange quark distribution $`\mathrm{\Delta }s(x)`$. Note that above equation is derived in both the flavor $`\mathrm{SU}(3)_f`$ case ($`\mathrm{\Delta }\overline{u}(x)=\mathrm{\Delta }\overline{d}(x)=\mathrm{\Delta }\overline{s}(x)`$) and the non–$`\mathrm{SU}(3)_f`$ case ($`\mathrm{\Delta }\overline{u}(x)=\mathrm{\Delta }\overline{d}(x)=\lambda \mathrm{\Delta }\overline{s}(x)`$***$`\lambda `$ represent a degree of $`\mathrm{SU}(3)_f`$ violation and is a parameter which should be determined from experiments .).

Setting a charm quark mass $`m_c=1.5`$ GeV and the relevant collider energy $`\sqrt{s}=50`$ GeV, we numerically calculated the spin–independent and dependent differential cross sections and the double spin asymmetry. As for the parton distribution functions, we used the GRV98  parametrization as the unpolarized parton distribution function, and the AAC  and “standard scenario” of GRSV00  parametrizations for the polarized one.

We show the double spin asymmetry in Fig. 1 as a function of $`p_T`$ at $`\sqrt{s}=50`$ GeV for $`\mathrm{\Lambda }_c^+`$ production (left panel) and for $`\overline{\mathrm{\Lambda }}_c^+`$ production (right panel). In order to suppress the contributions from the diffractive process and higher twist corrections, we have imposed the kinematical cut on $`p_T`$ as $`p_T>2`$ GeV in numerical calculations. In figures, the bold and normal lines show the case of AAC parametrization and “standard scenario” of GRSV00 parametrization, respectively. The solid lines represent the spin transfer coefficient $`C_c^{\mathrm{\Lambda }_c^+}(z)=1`$ case, while the dashed lines represent $`C_c^{\mathrm{\Lambda }_c^+}(z)=z`$ case.

As shown in figures, $`A_{LL}^{\mathrm{\Lambda }_c^+}`$ in the smaller $`p_T`$ regions does not depend on the model of polarized parton distribution functions, and is strongly affected by the shape of the spin transfer coefficient $`C_c^{\mathrm{\Lambda }_c^+}(z)`$. Therefore, The $`\mathrm{\Lambda }_c^+`$ production is effective for extracting information about the polarized fragmentation function $`\mathrm{\Delta }D_c^{\mathrm{\Lambda }_c^+}(z)`$. On the other hand, measuring $`A_{LL}^{\overline{\mathrm{\Lambda }}_c^+}`$ in larger $`p_T`$ regions is quite effective for testing the model of polarized parton distribution functions, since the ambiguity of the polarized fragmentation function is small and we see a big difference between two parametrization models.

In summary, to extract information about the polarized strange quark distribution in proton, the semi–inclusive $`\mathrm{\Lambda }_c^+`$/$`\overline{\mathrm{\Lambda }}_c^+`$ production in neutrino–polarized proton collisions was studied. $`\overline{\mathrm{\Lambda }}_c^+`$ production is most promising not only for testing but also for directly extracting the polarized strange quark distribution $`\mathrm{\Delta }s(x)`$ by measuring $`A_{LL}^{\overline{\mathrm{\Lambda }}_c^+}`$.


format: text

text:
# 1 Introduction

## 1 Introduction

About 20 years ago Koch and Pietarinen performed an energy-independent partial-wave analysis on pion-nucleon elastic and charge-exchange differential cross sections and elastic polarisations for laboratory momenta below 500 MeV/c incorporating the constraints from fixed-$`t`$ dispersion relations as well as crossing and unitarity (the KH78 and KH80 analyses) . Since then, however, new low-energy data have emerged in all charge channels: examples of recent high precision results for the differential cross sections are given in refs. , polarisation parameter in refs.  and the spin rotation parameter in ref. . There are also some new measurements of integrated cross sections . Especially the high precision measurements of the hadronic level shift and width on pionic hydrogen and deuterium , giving information of the pion-nucleon interaction just at the threshold, have opened a completely new chapter in the study of the low-energy pion-nucleon interaction. Another direction where significant advances have taken place is the theoretical framework where we study the low-energy pion-nucleon interaction. The tool, chiral perturbation theory ($`\chi `$PT), has been developed in the 80’s and 90’s and the work continues. The development of $`\chi `$PT motivates a new partial wave analysis from the theoretical point of view — on one hand strong interaction physics is becoming a precision science also in the low-energy region, on the other hand, there is need for $`\pi N`$ phenomenology to fix some of the low-energy constants appearing in the meson-baryon lagrangian. See, for example, the talk of Meißner in these proceedings .

It is our goal to make use of the new data in an analysis which, in addition to the requirements of analyticity, crossing and unitarity, includes the constraints from chiral symmetry.

## 2 The old KH78/80 analysis

The aim of the KH analysis was to determine the amplitudes satisfying several conditions:

1. The amplitudes had to reproduce all the data which were available at that time, which means $`\frac{d\sigma }{d\mathrm{\Omega }}`$, $`\sigma _{\mathrm{tot}}`$ and $`P`$ within their associated errors.
2. The solution had to fulfil the isospin invariance.
3. All partial waves had to satisfy the unitarity condition.
4. The crossing symmetry was implicitly assumed, because Mandelstam variables were used.
5. The invariant amplitudes were going to have the correct analyticity properties in $`s`$ at fixed-$`t`$.
6. The amplitudes at fixed-$`s`$ were to be analytic in $`\mathrm{cos}\theta `$ in the small Lehman ellipse.

The experimental data is not enough to fix a unique partial wave solution, but further theoretical constraints are needed. The constraints from fixed-$`t`$ analyticity and from the isospin invariance are strong enough to resolve the ambiguities .

## 3 Three stages of the analysis

### 3.1 Fixed-$`t`$ analysis

The old KH analyses consist of three phases: fixed-$`t`$ analysis, fixed-$`\theta `$ analysis and fixed-$`s`$ analysis, which were iterated until the results agreed up to about 3 %. The fixed-$`t`$ analysis was carried out at 40 $`t`$-values in the range from zero to $`1.0`$ GeV2. The analysis would be too complicated, if one were working with dispersion integrals, so the expansion

$$C(\nu ,t)=C_N(\nu ,t)+H(Z,t)\underset{i=0}{\overset{n}{}}c_iZ^i$$
(1)
was used for $`t`$-values smaller than $`4m\mu `$, ($`\nu =(su)/4m_N`$). In the expansion the nucleon pole term $`C_N(\nu ,t)`$ is treated separately, and the sum is multiplied by a factor $`H(Z,t)`$, which describes the expected asymptotic behaviour. The essence of the expansion is that the sum is written in terms of functions $`Z`$, which have the correct analytic behaviour, i.e. it is not a polynomial approximation, but a series presentation of an analytic function, which is just truncated at some reasonable point (ca. $`n=50`$ or $`n=100`$), because infinite accuracy is impossible. The condition of smoothness and the compatibility with the data constrain the terms with large index $`i`$ to be negligible .

The expansion coefficients $`c_i`$ are determined by minimising

$$\chi ^2=\chi _{\mathrm{data}}^2+\chi _{\mathrm{pw}}^2+\chi _{\mathrm{penalty}}^2.$$
(2)
Here $`\chi _{\mathrm{data}}^2`$ comes from the experimental errors, $`\chi _{\mathrm{pw}}^2`$ belongs to the deviation from the fixed-$`s`$ partial wave solution and the last term is used to suppress large values of the higher coefficients of the expansion. In practice, the analyticity constraints cannot be used without smoothing the data. The aim is, of course, to smooth out the statistical fluctuations without distorting the physically relevant structures.

### 3.2 Fixed-$`\theta `$ analysis

The fixed-$`t`$ constraint is often used only for $`t`$ values from zero to about $`0.5`$ GeV2, because the partial wave expansions for the imaginary parts of the invariant amplitudes do not converge for large $`|t|`$. In the range $`t(0.5,1.5]`$ GeV2 the truncated partial wave expansions can still be reasonable approximations, but for $`t`$ values below ca. $`1`$ or $`1.5`$ GeV2 the fixed-$`t`$ analyticity cannot be applied anymore. So another analyticity constraint is used to cover the rest of the angular range at intermediate and at high energies. The calculation was made at 18 angles between $`\mathrm{cos}\theta =0.9\mathrm{}0.8`$. Analysing methods are the same as in the fixed-$`t`$ analysis: the expansion method is used and the coefficients are fixed by minimising $`\chi ^2`$, i.e. by fitting to the data and to the fixed-$`t`$ solution.

### 3.3 Fixed-$`s`$ analysis

The third stage, the fixed-$`s`$ analysis, is a standard phase shift analysis in the sense that the partial waves are fitted to the data. On the other hand, it is not the usual one, because the partial waves are fitted also to the fixed-$`t`$ and to the fixed-$`\theta `$ amplitudes. Now 92 momentum values were selected from the energy range $`0\mathrm{}200`$ GeV/c, 6 of the momenta were above 6 GeV/c. Again, the coefficients were fixed by minimising $`\chi ^2`$ which now included also a term suitable to enforce unitarity.

## 4 Treatment of the data

The electromagnetic corrections proposed by Tromborg et al. at momentum values below $`0.65`$ GeV/c , were applied to the data. At higher momentum, only the one-photon exchange correction was applied, taking into account the Coulomb phase. In all three different analyses, the data were shifted to the selected energy bins (i.e. the selected values of $`s`$, $`t`$ or $`\theta `$) using the previous solution of the iteration to calculate the correction. Some data points requiring too large a momentum shift were omitted. The normalisation of some data sets had to be corrected to guarantee a smooth extrapolation to the forward direction, and to the input for the forward amplitude.

## 5 Life after KH80

The latest KH phase shift analysis was finished in 1980. After that there has been very accurate measurements of pionic hydrogen level shift and width, $`\frac{d\sigma }{d\mathrm{\Omega }}`$ has been measured with good accuracy, many spin rotation parameter measurements has been done as well as polarisation parameter measurements. Also, some integrated cross section measurements have been performed. The newer data has never been analysed by the methods of Koch and Pietarinen, and for example the results of Pavan et al.  are not compatible with the results of the old analyses. So, an updated version of the analysis is certainly needed.

## 6 The code of Pietarinen

The original code was made for the NDP Fortran compiler, which runs under MS-DOS. The code needs to be ported to UNIX. We have tested the code in an old MS-DOS machine, and most of the main tasks seem to be working correctly. What is still needed, is a modification of the $`s`$-plane conformal mapping. This has effects on all routines, which are related to the fixed-$`t`$ analysis. Also some modifications are needed to be able to study the isospin analysis.

The code base is divided into several parts:

* There is a program for comparing the partial wave solution to existing data. It simply plots the data and the solution in the same picture, and allows the comparison of different data sets to the solution.
* The second part is for shifting the experimental data points into the fixed-$`t`$ bins. The earlier solution is used for interpolation, and those data points which are to be shifted too much are rejected.
* The next part is for making a starting value for the fixed-$`t`$ expansion.
* One program is for the iteration to adjust the fixed-$`t`$ amplitude to the experimental data and to the current solution.
* The main part of the program makes the actual partial wave analysis and adjusts the solution simultaneously to the data and to the fixed-$`t`$ amplitudes.

## 7 The current status

### 7.1 Porting the code

During the porting process the most extensive work is needed for writing the graphical user interface, the data base engine and the plotting routines. The graphics routines of the original code were impossible to get working under UNIX, so we decided to use the Python/Tkinter combination for the GUI, and the Python/Gnuplot combination for plotting routines. The old code was reused as much as possible, but many parts still needed almost complete rewriting. We decided to write all the new code in Fortran 95, so, at the moment, most of the calculation engine is written in Fortran 77 and some parts in Fortran 95. All routines of the old code were modified to take almost all input from the stdin and to write output to stdout, so it should now be possible to change the whole GUI with a reasonable amount of work, whenever it becomes necessary.

### 7.2 The current status

At present, the plotting program, the interpolation program and the program calculating the starting value of the fixed-$`t`$ analysis are ported to UNIX, and all the functionality of the original versions is implemented (fig. 1). The part, which iterates to adjust the fixed-$`t`$ amplitudes, the partial wave solution, and the experimental data, compiles OK but the GUI is still under construction. The heart of the whole program, the part making the actual analysis, is still in a completely untested stage.

### 7.3 Comparison of the numerics

We have compiled the code with different compilers running on different operating systems in order to check the stability of the mathematical subroutines111We have tried the old NDP Fortran compiler for MS-DOS, the Compaq Fortran in Digital UNIX, the Lahey/Fujitsu Fortran for Linux, and different versions of the GNU/Linux Fortran.. Comparing the results of the routines compiled by different compilers showed that the routines are not currently stable enough for production use.

For illustration, the interpolator part of the program calculates 173670 interpolated data points. When comparing the results of routines compiled by GNU Fortran to those calculated by NDP Fortran, one notices that 96% of the new data points agree up to 0.01%, but in some cases there are significant discrepancies: namely 73 data points of the 173670 differ by more than 1%, and in the worst case the difference is 28%. The cause of these discrepancies is unknown when writing this.

## 8 The next phase

After chasing the bugs and finding the reasons for the numerical unstabilities, we are hoping to find a better way to handle the propagation of the experimental errors than that used in the old KH analysis.

Because of the size of the data base, and because of the discrepancies in the different data sets, the visualisation of the data and the partial wave solutions is essential. Also, during the analysis, the program is used a lot, so the graphical user interface has to be easy and efficient to use. For making the choice of the data sets as easy as possible, we have plans to convert our text file data bases to MySQL format.

We intend to get the first preliminary results by the end of the year.

Acknowledgments: We wish to thank A.M. Green for useful comments on the manuscript. We gratefully acknowledge financial support by the Academy of Finland grant 47678, the TMR EC-contract CT980169, and the Magnus Ehrnrooth foundation.


text:
# Pair annihilation and radio emission from galactic jet sources: The case of Nova Muscae

## 1 Introduction

It has been known for some time now that a handful of Galactic X-ray transients has exhibited radio jets. The first X-ray binary for which a kinematic model involving radio jets was proposed to explain its bizarre observed spectrum was SS 433 in 1979 (Abell & Margon 1979). Subsequently, in 1994, apparent superluminal motion of ejected material was observed for the first time in our Galaxy from GRS 1915$`+`$105 (Mirabel & Rodríguez 1994). The velocity of the ejecta was estimated to be 0.92$`c`$ for an assumed distance of 12.5 kpc. Since that first discovery, four other transients have shown apparent superluminal motion: GRO J1655$``$40 (Hjellming & Rupen 1994; Tingay et al. 1995); XTE J1748$``$288 (Rupen, Hjellming & Mioduszewski 1998; V4641 Sgr (Hjellming et al. 2000); and XTE J1550$``$564 (Hannikainen et al. 2001).

In addition to the transient jet sources, there has also been evidence of steady compact radio jets during the low hard state of Galactic black hole candidate X-ray binaries from, for example, Cyg X-1 (Stirling et al. 2001) and 1E 1740.7$``$2942, or “The Great Annihilator” (e.g. Mirabel et al. 1992; see Fender 2001 for a full review on radio jets from X-ray binaries in the low hard state). The presence of jets in X-ray binaries gives rise to radio emission and in this paper we will take the observation of radio emission to indicate the presence of jets. Note however, that in the absence of resolved radio observations, this emission may also stem from other regions within the source.

The composition of the jets, either electron-positron or electron-proton plasma, is still not fully established. For example, Gliozzi, Bodo & Ghisellini (1999) investigated the role of electron-positron pairs (amongst other possibilities) as energy carriers from the inner jet regions, and have excluded a pair plasma as a viable possibility. They argue that either hot or cold pairs cannot survive the annihilation.

In this paper we propose that the annihilation line features observed in the hard X-ray spectra of some X-ray transients arise from pairs in a bipolar outflow. At the time of annihilation, this outflow is already accelerated to relativistic bulk speeds causing a significant Doppler shift of the frequency of the annihilation lines. We also show that the subsequent emission of radio synchrotron radiation from the outflow may be caused by only a small fraction of pairs escaping from the annihilation region. We apply this idea to radio and hard X-ray observations of Nova Muscae 1991, thus constraining the properties of the possible outflow from this system. The energy requirements for this source rule out a large contribution of protons to the outflow. We would like to point out that 1E 1740.7$``$2942 has exhibited annihilation features near 511 keV (e.g. Mandrou et al. 1990), and that there is some evidence of the same in Cyg X-1 (Ling & Wheaton 1989). We are currently working on expanding the arguments presented in this paper to encompass the other two sources (Hannikainen & Kaiser, in preparation).

In Section 2 we discuss under which plasma conditions we can expect Doppler-shifted annihilation lines. Section 3 briefly reviews radio synchrotron emission. The case of Nova Muscae 1991 is investigated in Section 4 and in Section 5 we summarise our results.

## 2 Doppler-shifted annihilation lines

### 2.1 Conditions for a strong, narrow annihilation line

The direct annihilation of an electron-positron pair results in the production of two $`\gamma `$-ray photons, each with an energy of 511 keV in the rest-frame of the annihilating particles. In the case of a hot pair plasma the resulting emission line is broadened by the random thermal motion of the pairs (Ramaty & Mészáros 1981). Also, the annihilation of a thermal pair plasma is accompanied by bremsstrahlung from the pairs. The bremsstrahlung emissivity around 511 keV exceeds that of the annihilation process at temperatures above $`\mathrm{\Theta }=kT/(m_\mathrm{e}c^2)=3`$ (Svensson 1982; Maciołek-Niedźwiecki, Zdziarski & Coppi 1995, hereafter MZC95). Below a temperature of about $`10^6`$ K the annihilation proceeds via the formation of positronium leading to a distinctive low-energy wing of the annihilation line (e.g. Longair 1994). Thus in the case of a thermal pair plasma the plasma temperature must lie roughly in the range $`10^6`$ K to $`4\times 10^9`$ K for a narrow annihilation line to be observable.

In a plasma, pairs may be produced by various mechanisms. Unless there is a strong external radiation field containing many photons with energies in excess of 511 keV, all of these require the presence of highly energetic electrons and soft seed photons. The electrons upscatter the soft photons beyond the pair creation threshold and thus the plasma starts producing pairs. As a by-product of the pair production, a strong Comptonization spectrum is emitted by the plasma. MZC95 point out that this Comptonization spectrum completely swamps the annihilation line in thermal plasmas. An observable narrow line can only be produced if a large fraction of the pairs escapes from the plasma and has time to cool before annihilating. However, this implies the existence of a rather sharp boundary between a fairly hot pair plasma bathed in an intense radiation field of soft photons and a region virtually devoid of any soft photons. Also, although the pairs escape the plasma, the volume containing the pairs must be confined at least for a time longer than the cooling time. Any expansion would lead to enormously reduced annihilation rates. Although it is certainly possible to construct such geometries, the existence of non-thermal leptons within the plasma allows for a simpler way of producing a narrow annihilation line.

In non-thermal plasmas relativistic electrons or pairs may be injected into the plasma. The injected pairs and those produced in the plasma may cool to sub-relativistic energies and thermalize before annihilating, thus leading to a narrow annihilation line (Lightman & Zdziarski 1987, hereafter LZ87). If the line is strong, it can rise above the Comptonization spectrum and becomes detectable. This requires a high pair yield, $`Y`$, defined as the ratio of the energy converted to pairs and the energy supplied to the plasma. The highest pair yields can be achieved when the plasma is ‘photon-starved’, i.e. when the number of injected relativistic photons strongly exceeds that of the injected soft photons (Zdziarski, Coppi & Lamb 1990). In this case, $`Y0.25`$ and a strong, narrow annihilation line above the Comptonization continuum becomes observable.

We conclude that the observation of a narrow annihilation line most likely indicates a plasma with strong injection of non-thermal electrons or pairs. The injection of leptons into a spherical volume of radius $`R`$ is characterised by the compactness (e.g. LZ87)

$$l_\mathrm{e}=\frac{L_\mathrm{e}\sigma _\mathrm{T}}{Rm_\mathrm{e}c^3},$$
(1)
where $`\sigma _\mathrm{T}`$ is the Thomson cross-section and $`L_\mathrm{e}`$ is the power of the electron injection,

$$L_\mathrm{e}=\frac{4\pi R^3}{3}m_\mathrm{e}c^2Q_\mathrm{l}(\gamma )\left(\gamma 1\right)𝑑\gamma .$$
(2)
Here, $`Q_\mathrm{l}(\gamma )`$ is the rate of injection of leptons with Lorentz factor $`\gamma `$ per unit volume per unit time per unit $`\gamma `$. If the compactness $`l_\mathrm{e}`$ can be inferred from the observations of an annihilation line, then Equations (1) and (2) can be used to constrain $`Q_\mathrm{l}(\gamma )`$.

### 2.2 Relativistic Doppler-shifts

Any emission of relativistically moving material is Doppler-shifted in its frequency. For material moving with bulk velocity $`v_\mathrm{b}=\beta c`$ at an angle $`\theta `$ to the line of sight to the observer, the observed frequency, $`\nu `$, of radiation emitted at frequency $`\nu ^{}`$ in the rest-frame of the material is given as

$$\nu =\frac{\nu ^{}}{\gamma _\mathrm{b}\left(1\pm \beta \mathrm{cos}\theta \right)}=\nu ^{}\delta _\pm ,$$
(3)
where $`\gamma _\mathrm{b}`$ is the Lorentz factor corresponding to the velocity $`\beta `$ (e.g. Rybicki & Lightman 1979). The upper signs correspond to material receding along the line of sight to the observer while the lower signs indicate approaching material. From Equation (3) it is clear that radiation of material receding from an observer is always redshifted. However, for approaching material the emission may be blueshifted or redshifted, depending on the combination of $`\beta `$ and $`\theta `$. Figure 1 shows this effect for the example of positron annihilation radiation emitted at a rest-frame energy of 511 keV. Any combination of $`\beta `$ and $`\theta `$ below the dashed line in this figure corresponds to redshifted emission from the approaching material. The range of velocities which result in such a Doppler redshift is largest for angles to the line of sight close to $`90^{}`$.

For a source with a bipolar outflow with components travelling in opposite directions there will always be one approaching and one receding component. If both components emit a strong annihilation line, we expect to observe two such lines at frequencies $`\nu _{}`$ and $`\nu _+`$. Assuming the same bulk velocity for both components, Equation (3) can then be solved for the bulk velocity of the components,

$$\gamma _\mathrm{b}=\frac{\nu ^{}}{2}\left(\frac{1}{\nu _{}}+\frac{1}{\nu _+}\right),$$
(4)
and for the angle of the bulk velocity to the line of sight

$$\mathrm{cos}\theta =\frac{1}{\beta }\left(\pm \frac{\nu ^{}}{\gamma _\mathrm{b}\nu _\pm }1\right),$$
(5)
where again the upper signs correspond to the receding component and the lower signs to the approaching one.

In practice it may be difficult to observe the annihilation line from the receding component as this will be strongly deboosted. In fact, the ratio of the line fluxes provides for an independent check on the quantity $`\beta \mathrm{cos}\theta `$ as (Rybicki & Lightman 1979)

$$\frac{F_{}}{F_+}=\left(\frac{1+\beta \mathrm{cos}\theta }{1\beta \mathrm{cos}\theta }\right)^k,$$
(6)
with $`k=3`$ for discrete components and $`k=2`$ for continuous jets. It is usually difficult to determine the line fluxes accurately from observations. This, combined with the ambiguity of the value of $`k`$, implies that the line flux ratio is of limited use.

## 3 Radio synchrotron emission and absorption

Synchrotron emission arises from relativistic electrons and/or positrons spiralling in a magnetic field. It is usually described by a power law spectrum with $`L_\nu \nu ^\alpha `$. Both the relativistic particles and the magnetic field store energy, but how the total energy is divided between the two cannot be determined from the observed radiation alone. However, it is well known that there is a minimum total energy required to produce a given monochromatic synchrotron luminosity, $`L_\nu `$. This is usually expressed in terms of the magnetic field corresponding to this condition

$$B_{\mathrm{min}}=\left[\frac{3\mu _0G(\alpha )L_\nu }{2V}\right]^{2/7},$$
(7)
where $`V`$ is the volume of the emitting region and $`G(\alpha )`$ is a numerical constant depending on the slope of the synchrotron spectrum, $`\alpha `$, and on $`\nu _{\mathrm{min}}`$ and $`\nu _{\mathrm{max}}`$, the limits of the spectrum. It is usually assumed that $`\nu _{\mathrm{min}}10`$ MHz and $`\nu _{\mathrm{max}}100`$ GHz. Note here that the exact choice for the spectral cut-offs, and particularly the choice for $`\nu _{\mathrm{max}}`$, has very little influence on the value of $`G(\alpha )`$ for $`\alpha 0.50.9`$, typical for optically thin synchrotron radiation (e.g. Longair 1994). The energy in the relativistic particles is then given by

$$W_{\mathrm{part}}=\frac{2VB_{\mathrm{min}}^2}{3\mu _0}=Vm_\mathrm{e}c^2P(\gamma )(\gamma 1)𝑑\gamma ,$$
(8)
where $`P(\gamma )`$ is the number of relativistic particles with Lorentz factor $`\gamma `$ per unit volume per unit $`\gamma `$. For a given synchrotron luminosity and a known volume of the emission region we can therefore constrain $`P(\gamma )`$. Obviously, this is just an estimate as the distribution of the energy between the relativistic particles and the magnetic field may depart strongly from the minimum energy values. However, such a departure would lead to much higher overall energy requirements.

The estimates presented above implicitly assume that the emission region is optically thin to radio waves. This is not necessarily the case, particularly in the inner, dense regions of the source. For $`P(\gamma )=P_0\gamma ^p`$ the optical depth is

$$\tau =F(p)n_0B_{\mathrm{min}}^{(p+2)/2}\nu ^{(p+4)/2}R,$$
(9)
where $`F(p)`$ is a function of $`p`$ only (e.g. Longair 1994) and $`R`$ is the typical path length photons have to travel through the emission region. If the magnetic field is tangled on scales significantly smaller than the physical size of the emission region, its energy density obeys $`u_{\mathrm{mag}}V^{4/3}`$. So for a spherical volume of radius $`R`$ and $`VR^3`$ we find

$$\tau \left(R^2\nu ^{}\right)^{(p+4)/2}.$$
(10)
## 4 The case of Nova Muscae

Nova Muscae (GRS 1124$``$684) was discovered on 1991 January 8 by Granat (Lund & Brandt 1991; Sunyaev 1991) and Ginga (Makino 1991). An extensive monitoring campaign with Granat covering most of January and February 1991 followed this initial discovery (Gil’fanov et al. 1991, hereafter GSC91; Sunyaev et al. 1992; Gil’fanov et al. 1993). During observations on January 20–21, a strong, narrow line near 500 keV was detected. This line had not been noted before and the line flux was observed to increase during the observation within the space of a few hours. Unfortunately, the observations stopped before the line flux decreased again, implying a lifetime of the line emission of at least 10 hours. During subsequent observations on February 1–2 it was not detected. Simultaneously to the strong line near 500 keV, there was also increased emission detected near 200 keV. The spectrum obtained with Granat during 1991 January 20–21 is shown in Figure 2. GSC91 fitted the spectrum with a power-law continuum and two Gaussian line profiles. The lines peak at $`474\pm 30`$ keV and $`194\pm 13`$ keV. The line fluxes are given as $`(6.6\pm 3.6)\times 10^3`$ photons s-1 cm-2 and $`(1.5\pm 1.3)\times 10^3`$ photons s-1 cm-2 respectively. The FWHM of the line centred at 474 keV is $`70\pm 70`$ keV.

The original discovery also triggered a radio monitoring programme using the Molongolo Observatory Synthesis Telescope (MOST) and the Australia Telescope Compact Array (ATCA) at frequencies ranging from 843 MHz to 8.6 GHz (Kesteven & Turtle 1991; Ball et al. 1995). The campaign started on January 17 and continued into early February. The radio lightcurve at 843 MHz obtained with the MOST is shown in Figure 3. Ball et al. (1995) note four distinct features in the lightcurve. There is a general decay of the radio flux from the very beginning of the programme continuing until about January 22. This is interrupted by a brief flare around January 18 observed with ATCA (not shown in Figure 3). There is another flare observed from January 31 until around February 5. Finally, a short flare was observed by MOST at 843 MHz on January 24 with a measured flux density of 24 mJy. The event lasted only one day. There were no observations at other frequencies.

In the following we will concentrate on the detection of the $`\gamma `$-ray lines on January 21 and the brief radio flare on January 24. We speculate that the same ejection event is responsible for the emission at opposite ends of the electromagnetic spectrum observed on the two days.

### 4.1 Constraints from the annihilation line

We assume here that the two $`\gamma `$-ray lines observed in the spectrum of Nova Muscae arise from pair annihilation in a relativistic bipolar outflow from the very centre of the system. In this model, the line at 474 keV is associated with the approaching component while that at 194 keV arises from the receding component.

Using Equations (4) and (5) we then find $`\theta =60^{}\pm 7^{}`$ for the angle to the line of sight of the component motion and $`\beta =0.84\pm 0.02`$ for the bulk velocity of the components. This corresponds to a Lorentz factor of $`\gamma _\mathrm{b}=1.86\pm 0.12`$. The errors are determined from the uncertainties in the frequency of the line peaks as given by GSC91. With these values we expect the ratio of the line fluxes of the approaching and the receding components to be $`F_{}/F_+=15\pm 10`$ for discrete components and $`F_{}/F_+=6\pm 3`$ for continuous jets (see Equation 6). The observed line fluxes give $`F_{}/F_+=4\pm 6`$. The velocity and angle to the line of sight calculated from the frequencies of the line peaks are consistent with this within the errors. However, the uncertainties are very large.

The companion star in Nova Muscae is tidally distorted by the black hole. The resulting variation in the lightcurve has been modelled by Orosz et al. (1997) who constrain the inclination of the system to $`54^{}i65^{}`$, and also by Shahbaz, Naylor & Charles (1997) who find $`i=54_{15}^{+20}`$. It is reasonable to expect that $`i`$ is also the inclination angle of the accretion disc around the black hole and that the bipolar outflow moves in a direction perpendicular to this disc. However, note that the jets in GRO J1655-40 appear to be misaligned with respect to the direction perpendicular to the disc (e.g. Orosz & Baylin 1997). With these assumptions we should find $`\theta =i`$, i.e. the angle of the outflow direction to the line of sight is equal to the inclination of the system. Our value for $`\theta `$ is fully consistent with this interpretation.

We now take a closer look at the profile of the line centred at 474 keV. The line at 194 keV is too weak to allow any detailed investigation of its profile, and hence we only fit the line at 474 keV for this present study. The FWHM of the line at 474 keV is only 70 keV which implies a pair temperature of only about $`4\times 10^7`$ K at the time of annihilation (Ramaty & Mészáros 1981; GSC91).

The rate of annihilation implied by the line flux is $`\dot{N}_+2\times 10^{43}`$ pairs s-1 for a distance of 5.5 kpc (Orosz et al. 1996). This enormous rate is sustained for at least 10 hours (Sunyaev et al. 1992). It is therefore very unlikely that this feature at 474 keV arises from a large number of pairs formed practically instantaneously and then slowly annihilating away. This would imply a short-lived flash of annihilation photons with a fast, exponential decay contrary to the lifetime of the annihilation features of at least 10 hours. A more promising approach is to assume that the pair producing plasma is in equilibrium, i.e. the annihilation losses are balanced by pair creation. MZC95 failed to fit the line with such a model assuming a thermal plasma with strong pair escape. They also tried a model with injection of non-thermal electrons based on the work of LZ87. Again, the fit was poor. However, they did not take into account a possible Doppler-redshift of the line.

After correcting the spectrum for the redshift derived above, we attempted a fit essentially identical to that of MZC95. We used the model NTEEA within the XSPEC package which is an implementation of the model developed in LZ87. We disabled any contribution to the spectrum by reflection off the disc. The temperature of the black body providing the soft seed photons for the pair creation process was set to 1 keV which gives a reasonable fit to the low energy spectrum of Nova Muscae (see MZC95). Further, we assumed that there was no thermal heating of the plasma, no escape of pairs from the source and that the optical depth due to protons was zero. The energy spectrum of the injected relativistic electrons was described by $`Q_\mathrm{l}=Q_0\gamma ^{2.2}`$, extending from $`\gamma _{\mathrm{min}}=1`$ to $`\gamma _{\mathrm{max}}=1000`$. We have chosen the slope of the power law to be equal to that of the radio emission (see Section 4.2). Variation of this parameter as well as changing $`\gamma _{\mathrm{min}}`$ and $`\gamma _{\mathrm{max}}`$ was found to cause only minor differences in the result. Finally, the model assumes that the plasma is contained within a spherical volume. The radius of this we set to $`R=2\times 10^5`$ m corresponding to the Schwarzschild radius of a black hole of mass $`6`$ M. This is the minimum conceivable size of the annihilation region. In any case, the exact value of $`R`$ is only important for the determination of energy losses due to Coulomb and bremsstrahlung processes. It does not significantly influence the model fit of the annihilation line. The only free parameters in this fit were the compactness of the non-thermal electron injection, $`l_\mathrm{e}`$, and the compactness of injected soft photons, $`l_\mathrm{s}`$. The latter is defined to be analogous to $`l_\mathrm{e}`$ with $`L_\mathrm{e}`$ in Equation (1) replaced by the luminosity of the soft photon injection. As we were only interested in fitting the annihilation line at 474 keV, we only used data points in the range 195–699 keV. We also excluded all negative flux measurements within this range which resulted in 23 data points to be fitted.

Figure 2 shows our best fit with a reduced $`\chi ^2`$-value of 1.35 (for 84 d.o.f). The free parameters, $`l_\mathrm{e}3000`$ and $`l_\mathrm{s}50`$, imply a strongly photon-starved plasma. The parameters are not well constrained as under these conditions the pair yield does not depend strongly on their exact values (e.g. LZ87). We find a lower limit for $`l_\mathrm{e}`$ of about 100, below which the pair yield is too small to give a good fit. Despite these uncertainties, the correction of the data for a Doppler-redshift clearly makes a reasonable fit of the annihilation line possible.

Re-arranging Equation (1) yields

$$L_\mathrm{e}=\frac{Rm_\mathrm{e}c^3l_\mathrm{e}}{\sigma _\mathrm{T}}.$$
(11)
Substituting in our lower limits for $`R`$ and $`l_\mathrm{e}`$, we find $`L_\mathrm{e}7\times 10^{30}`$ W. From Equation (2) it follows that $`Q_0=9\times 10^{26}`$ m-3 s-1 for this lower limit. The rate at which relativistic electrons are injected is then $`\dot{N}_{\mathrm{inj}}=3\times 10^{43}`$ particles s-1. This is comparable to the observed annihilation rate of $`2\times 10^{43}`$ particles s-1.

### 4.2 Constraints from the radio emission

The radio flare occured on 1991 January 24, about three days after the observation of the annihilation lines. If this emission arose from the same material that was responsible for the annihilation features at the very centre of the source, then these ejection components were by then located roughly $`7\times 10^{13}`$ m further out. This assumes that the components did not decelerate. When resolved in the radio, the ejections of X-ray transients are quite collimated (e.g. Mirabel & Rodríguez 1994; Hjellming & Rupen 1995). In order to estimate the size of the radio-emitting region, we assume a half-opening angle of the ejection of $`1^{}`$. This implies a radius of $`10^{12}`$ m for a spherical component. The measured flux density was $`23.9`$ mJy at 843 MHz which we associate entirely with the approaching component. For the values found in Section 4.1 for the bulk velocity, $`\beta =0.84`$, and the angle to the line of sight, $`\theta =60^{}`$, this converts to 29.9 mJy emitted at 909 MHz in the rest frame of the approaching component. This assumes a discrete component rather than a continuous jet. We have no information on the slope of the radio spectrum for 24 January. Therefore, we adopted a spectral index of $`\alpha =0.6`$ as derived from the multi-frequency observations of the earlier flares. This implies that $`P(\gamma )=P_0\gamma ^{2.2}`$. For a distance to the source of 5.5 kpc (Orosz et al. 1996) and using Equation (7), we find the magnetic field strength correponding to minimum energy conditions to be $`B_{\mathrm{min}}=0.015`$ mT. Using Equation (8) this implies $`W_{\mathrm{part}}=5\times 10^{32}`$ J and $`P_0=8\times 10^8`$ m-3. The number of relativistic particles in the emitting volume is $`N_{\mathrm{rel}}=2\times 10^{44}`$.

Adopting the radius of the emission region at the time of the observation as a scale radius $`R_0=10^{12}`$ m, Equation (9) for the optical depth at $`\nu ^{}=909`$ MHz yields

$$\tau 0.1\left(\frac{R}{R_0}\right)^{6.2}.$$
(12)
### 4.3 Implications of a bipolar flow

The annihilation feature observed persists for at least 10 hours (Sunyaev et al. 1992). The production of an outflow must have lasted at least for that length of time. In the case of a shorter ejection event the annihilation rate would decrease dramatically as the ejected material travels outwards and expands. Of course, the annihilation could also proceed in a region close to the centre of the source without moving outwards. In this case, the displacement of the line at 474 keV could be explained as arising from gravitational redshift. This requires the annihilation region to be inside about 7 Schwarzschild radii of the 6 M black hole. However, in this case it is not clear what causes the second, simultaneous line at 194 keV.

An entirely different explanation for the line at 474 keV is given by GSC91 and also by Martín et al. (1994). This involves the emission of $`\gamma `$-ray photons of energy 468 keV during the formation of 7Li close to the black hole. This is supported by the observation of optical lithium lines. However, GSC91 point out that the proximity of the lithium production region to the black hole should result in a strong rotational broadening of the line. Furthermore, the secondary line at 194 keV is again not explained in this scenario.

Misra & Melia (1993) presented a model for electron-positron pair jets to explain the annihilation line observed in the spectrum of the Galactic centre source 1E 1740.7-2942. Their model is based on a radiative acceleration mechanism for the jets. The great majority of pairs annihilate within the acceleration zone with a large variation in the thermal pair temperature. This leads to the formation of quite a broad line. Also, as mentioned before, a model with purely thermal pairs is unable to account for the annihilation line in the case of Nova Muscae (MZC95). In the model presented here we argue that the bulk of the outflow containing the pairs is already accelerated when they annihilate, thus explaining the redshifts of the two observed lines. The velocity of the outflow is then 0.84 c, while the radiative acceleration allows only for a terminal velocity of about 0.7 c. Because of this, the acceleration of the outflow must, at least partly, be due to some other mechanism(s).

The energy required to drive the outflows is enormous. As the relativistic electrons necessary for pair production are highly relativistic, their mass in the rest frame of the outflow material is given by $`\gamma m_\mathrm{e}`$. Therefore, the total kinetic power of the relativistic electrons, as measured in the source rest frame, injected into the outflow can be approximated as

$$\dot{E}_{\mathrm{kin}}=VQ_0m_\mathrm{e}c^2\left(\gamma _\mathrm{b}1\right)\gamma ^{1p}𝑑\gamma 6\times 10^{31}\mathrm{W}.$$
(13)
In this estimate we have used the limiting values for the electron compactness, $`l_\mathrm{e}=100`$, and the size of the emission region, $`R=2\times 10^5`$ m. Any increase in these parameters will also cause the energy estimate to rise. Therefore, the lower limit on $`\dot{E}_{\mathrm{kin}}`$ corresponds already to about 80% of the Eddington luminosity of a 6 M black hole. Balance of electrical charge requires the presence of positively charged particles in the outflow. If these are protons with negligible thermal energy but travelling at the necessary bulk velocity, then another $`4\times 10^{33}`$ W in kinetic energy are required. This is more than a factor 50 more than the Eddington luminosity. This energy injection into the outflow has to be sustained for more than 10 hours and thus makes a large proton content in the jet very unlikely.

The only alternative is then that from the very beginning of the ejection the outflow material in the annihilation region consists of virtually a pure pair plasma. This removes the requirement of relatively inefficient pair production from relativistic electrons. The pairs must be injected into the relativistic bulk flow with a relativistic velocity distribution to explain the strength of the annihilation line (MCZ95), the cooling within the outflow and then annihilation at a rate of $`2\times 10^{43}`$ s-1. In this case, the required kinetic power is about $`6\times 10^{30}`$ W, or 8% of the Eddington luminosity.

Assuming that the outflows consist entirely of pairs, at least $`7\times 10^{47}`$ pairs are injected into the outflow approaching us. Most of these annihilate but a small fraction may escape the annihilation region. If the radio flare discussed above is associated with the annihilation line feature and the radio-emitting region only contains the relativistic pairs required by minimum energy arguments, then only one in 7000 pairs need escape. In the model of Misra & Melia (1993), roughly 10% of all pairs escape, which would imply a large thermal population underlying the radio-emitting relativistic pairs. Nevertheless, the kinetic energy of this thermal population is, at $`3\times 10^{29}`$ J, small compared to the kinetic energy of $`10^{33}`$ J of the relativistic and therefore heavy radio-emitting particles. In any case, the energy required to explain the radio observations is at most only a fraction of 5% of the originally injected energy. Thus the radio emission observed in X-ray transients may represent only an almost negligible fraction of the energy initially available to the outflow.

The time between the detection of the annihilation line features from Nova Muscae and the first observations of a radio flare is roughly three days. During this time, a number of observations at various frequencies showed no sign of enhanced radio emission (Ball et al. 1995). This would argue against a connection between the two events. However, Equation (12) shows that the outflow became optically thin to emission around 1 GHz only when $`R0.7R_0`$. For a constant opening angle of the outflows as assumed throughout, this implies that radio observations would have detected this flare only two days after its start at the earliest.

The existence of highly relativistic electrons far out from the centre of the source implies a secondary acceleration mechanism. The annihilation line at 474 keV is consistent with a pair temperature of only 5 keV. Therefore the escaping pairs are sub-relativistic and must be re-accelerated to explain the radio emission. One possibility is the existence of internal shocks within the outflow (e.g. Kaiser, Sunyaev & Spruit 2000). The duration of the outflow production is easily long enough to allow for variations in the outflow speed which may lead to the formation of such shocks. Alternatively, the particles may be re-accelerated at a working surface at the end of the flow, as is the case in many extragalactic jet sources.

## 5 Conclusions

From the above discussion it is likely that Nova Muscae underwent a number of bursts resulting in relativistic outflows. The production of these outflows would then always be accompanied by strong annihilation line emission followed by radio synchrotron emission a few days later. Here we are assuming that the radio emission arises in jets. Unfortunately, only one of the rather minor outbursts of Nova Muscae in 1991 was covered by both $`\gamma `$-ray and radio observations. The strongest outburst leading to the first detection of the system on 1991 January 8 was observed with insufficient energy resolution to establish the presence of an annihilation line. Subsequent high-resolution observations on January 9 did not detect any lines. There were two more minor radio bursts. One of them occurred on January 18 and a longer burst started around January 31 (Ball et al. 1995). If the three days’ delay between the appearance of the annihilation lines and the start of the radio emission is typical, then the two radio flares should have been preceeded by annihilation lines on January 15 and around January 28 respectively. There were no $`\gamma `$-ray observations during the period from January 11 until January 16 and from January 22 until February 1 (GSC91). The non-detections of January 9 and January 16 may be taken as evidence that the outflow production, when it occurs, is rather short-lived.

If it is true that the annihilation features observed originate from accelerated plasma, then protons as a contributing factor to the matter in the plasma are ruled out. The scenario we have discussed favors a pair plasma — even if only a tiny fraction of the pairs survives, this is sufficient to produce the observed radio emission. Furthermore, the detection of Doppler-shifted annihilation lines is a direct method, independent of uncertainties arising from other observational parameters such as fluxes, by which to directly measure jet speeds.

However, to confirm the outflow scenario linking $`\gamma `$-ray line observations with radio emission, further detailed, time-resolved observations are needed. The spectral and temporal resolution of integral will enable the detection and follow the evolution of possible annihilation features in other jet sources. This will prove extremely valuable in the case of, for example, GRS 1915$`+`$105 which is heavily obscured at visible wavelengths. In addition, this study was applied to Nova Muscae which is a transient, and hence it will be interesting to observe the behaviour of permanent jet sources, such as 1E 1740.7$``$2942, from which an annihilation line has already been detected.

## Acknowledgements

We would like to thank Marat Gil’fanov for providing us with the original GRANAT data. It is also a pleasure to thank Andrzej Zdziarski for help with NTEEA. We also thank Duncan Campbell-Wilson and Lewis Ball for so generously handing us the MOST and ATCA data for this project. The MOST is operated by the University of Sdyney and funded by grants from the Australian Research Council. We thank the anonymous referee for helpful comments. DCH acknowledges the support of a PPARC postdoctoral research grant to the University of Southampton.


text:
# Modified BCS mechanism of Cooper pair formation in narrow energy bands of special symmetry III. Physical interpretation

## I Introduction

In Part I of this paper josn  a modified BCS mechanism of Cooper pair formation was proposed that operates in narrow “superconducting” energy bands ($`\sigma `$ bands) and can be established within a nonadiabatic extension of the Heisenberg model of magnetism hei . This “nonadiabatic Heisenberg model” (NHM) is defined by three postulates given, e.g., in enhm .

The present paper gives a physical interpretation of this mechanism. In Section II the special properties of the spin-dependent Wannier functions (spin-dependent Wfs) in narrow $`\sigma `$ bands are interpreted in terms of spin-flipping processes. In Sections III and IV the spin-phonon coupling in narrow $`\sigma `$ bands and the resulting mechanism of Cooper pair formation, respectively, are illustrated in terms of these spin processes.

In Section V the important difference between the new mechanism in narrow $`\sigma `$ bands and the conventional BCS mechanism bcs  is interpreted in terms of “spring-mounted” Cooper pairs.

## II Spin-flipping processes

Assume a metal to be given that possesses a narrow, roughly half-filled $`\sigma `$-band complex in its calculated band structure. As in Part I we assume that this metal has only one atom in the unit cell. That means that the $`\sigma `$-band complex consists only of a single $`\sigma `$ band.

The Bloch functions of the $`\sigma `$ band can be transformed into spin-dependent Wfs $`w_m(\stackrel{}{r}\stackrel{}{T},t)`$ \[as given in Eq. (A14) of Part I\] with the following properties:

(i) they are centered on the atomic positions $`\stackrel{}{T}`$;

(ii) they are obtained by a unitary transformation from the Bloch functions of the $`\sigma `$ band;

(iii) they are adapted to the symmetry of the considered metal; and

(iv) they are as well localized as possible.

($`m=\pm \frac{1}{2}`$ is the quantum number of the crystal spin; $`\stackrel{}{r}`$ and $`t`$ stand for the local and spin coordinate, respectively, of the electron.) However, the spin-dependent Wfs cannot be written as a product of a local function $`w(\stackrel{}{r}\stackrel{}{T})`$ with the spin functions $`u_s(t)`$ ($`s`$ denotes the spin quantum number) because the matrix $`[f_{sm}(\stackrel{}{k})]`$ in Eq. (A9) of Part I cannot be chosen independent of $`\stackrel{}{k}`$ in a $`\sigma `$ band. Thus, these functions represent localized states $`|\stackrel{}{T},m`$ with spin directions depending on the positions $`\stackrel{}{r}`$ of the electrons. This motion may be described in terms of “spin-flipping processes” as, e.g., illustrated in Fig. 1. Here two electrons with parallel spins align their spins antiparallel when they approach within the region of overlap.

It is essential that within the NHM the states $`|\stackrel{}{T},m`$ are really occupied by electrons. This is a consequence of the second postulate of the NHM, i.e., of the special electron-electron interaction in the nonadiabatic system which is forced by this postulate. Hence, within the NHM the spin-flipping processes in a narrow $`\sigma `$ band are produced by a real electron-electron interaction that operates in the nonadiabatic system.

The flippings of the spins are a consequence of the special mathematical properties of the spin-dependent Wfs $`w_m(\stackrel{}{r}\stackrel{}{T},t)`$ which are determined by the four above conditions (i) - (iv) that must be satisfied within the NHM. Hence, we understand the physical character of this interaction if we understand why the $`w_m(\stackrel{}{r}\stackrel{}{T},t)`$ must comply with these conditions.

The first and the fourth condition (i) and (iv) are required by the second postulate. Hence, these conditions allow the electrons to take the state of lowest Coulomb repulsion energy. The second condition (ii) yields a nonadiabatic Hamiltonian with the correct symmetry and the third condition (iii) guarantees that the spin-dependent Wfs are orthogonal. Hence, the third condition is required that electrons occupying the states $`|\stackrel{}{T},m`$ comply with the Pauli principle.

Consequently, in the ground state of a narrow $`\sigma `$ band the electrons change permanently their spin directions in order to comply with the Pauli principle.

## III Spin-phonon interaction

In a narrow $`\sigma `$ band, the operator $`H_{\mathrm{Cb}}^{n\sigma }`$ of Coulomb interaction has the rather complicated form given in Eq. (3.13) of Part I. $`H_{\mathrm{Cb}}^{n\sigma }`$ depends on boson as well as on fermion operators and represents a special “spin-phonon interaction” which may be understood as follows.

On a spin-flipping process as indicated in Fig. 1 the spin angular momentum of the two electrons is not conserved. During this process the electrons behave like gyroscopes: inertial forces act on the electrons that deform their orbitals within their localized states. As indicated in Fig. 2, the localized charge distributions become asymmetric with respect to the lattice and accelerate the positive ions in such a way that phonons are emitted or absorbed. Eq. (3.14) of Part I shows that $`H_{\mathrm{Cb}}^{n\sigma }`$ conserves the total crystal spin angular momentum.

## IV Cooper pair formation

As a consequence of the spin-phonon coupling in a narrow $`\sigma `$ band, crystal-spin-1 phonons are excited in the ground state of the nonadiabatic system. We may assume that at zero temperature these phonons are virtually excited, i.e., each phonon pair is reabsorbed immediately after its generation. Therefore, in Part I we approximated the nonadiabatic Hamiltonian by a purely electronic Hamiltonian $`H^0`$ in order to describe the electron-phonon system at zero temperature; see Eq. (3.17) of Part I.

$`H^0`$ satisfies the equation of constraint (3.36) in Part I. This equation distinguishes the mechanism of Cooper pair formation within the NHM from the conventional BCS mechanism. As a consequence of this equation of constraint, in all the eigenstates of $`H^0`$ the electrons form Cooper pairs. This may be understood as follows.

The system represented by $`H^0`$ is purely electronic though the spin-flipping processes still occur. The electron spins still couple to the motion of the atomic cores and, hence, to the phonons. These phonons, however, cannot store crystal spin angular momenta but must be reabsorbed immediately after their generation. In this way, the phonons mediate a special spin-spin interaction between the electrons.

Thus, we have to look for a state with the total spin angular momentum being conserved within the electron system alone. That means, within the system represented by $`H^0`$ every spin-flipping process of two $`\sigma `$-band electrons must be accompanied by the opposite flipping process of two other $`\sigma `$-band electrons. This opposite process may occur somewhere in the crystal since the electron spins interact via virtual crystal-spin-1 phonons.

In this context, “opposite” means “time inverted” because the time-inversion symmetry is the only symmetry which makes such an opposite motion possible. Hence, the total spin angular momentum is conserved within the electron system if all the electrons form pairs that are invariant under time inversion. That means, whenever a certain one-electron state is occupied, then its time-inverted state must be occupied, too. In other words, all the electrons must form Cooper pairs.

## V Constraining forces

We obtain further information on the spin-flipping processes in a narrow $`\sigma `$ band by transforming the operator $`H_{\mathrm{Cb}}^{n\sigma }`$ of Coulomb interaction into the $`\stackrel{}{k}`$ representation,

$$H_{\mathrm{Cb}}^{n\sigma }=\underset{\stackrel{}{k},m}{}\stackrel{}{k}_1^{\prime \prime },l_1;\stackrel{}{k}_2^{\prime \prime },l_2;\stackrel{}{k}_1,m_1,n;\stackrel{}{k}_2,m_2,n|H_{\mathrm{Cb}}|\stackrel{}{k}_1^{},m_1^{},n;\stackrel{}{k}_2^{},m_2^{},nb_{\stackrel{}{k}_1^{\prime \prime }l_1}^{}b_{\stackrel{}{k}_2^{\prime \prime }l_2}^{}c_{\stackrel{}{k}_1m_1}^nc_{\stackrel{}{k}_2m_2}^nc_{\stackrel{}{k}_2^{}m_2^{}}^nc_{\stackrel{}{k}_1^{}m_1^{}}^n+\text{H.c.}$$
(1)
with

$$\stackrel{}{k}_1^{\prime \prime }+\stackrel{}{k}_2^{\prime \prime }+\stackrel{}{k}_1+\stackrel{}{k}_2=\stackrel{}{k}_1^{}+\stackrel{}{k}_2^{}+\stackrel{}{K}$$
(2)
where $`\stackrel{}{K}`$ stands for a vector of the reciprocal lattice. As already mentioned above, in addition $`H_{\mathrm{Cb}}^{n\sigma }`$ conserves the crystal spin angular momentum. In bcc materials this conservation law may be expressed by the analogous equation es3

$$l_1+l_2+m_1+m_2=m_1^{}+m_2^{}+4n$$
(3)
where $`n`$ is an integer. The boson operators $`b_{\stackrel{}{k}l}^{}`$ and $`b_{\stackrel{}{k}l}`$ create and annihilate phonons with wave vector $`\stackrel{}{k}`$ and crystal spin $`l=1,0,+1`$; the fermion operators $`c_{\stackrel{}{k}m}^n`$ and $`c_{\stackrel{}{k}m}^n`$ create and annihilate nonadiabatic Bloch states with wave vector $`\stackrel{}{k}`$ and crystal spin $`m=\pm \frac{1}{2}`$.

Remember that $`H_{\mathrm{Cb}}^{n\sigma }`$ depends on boson operators because during a spin-flipping process of two electrons a phonon pair is emitted (or absorbed) in a narrow $`\sigma `$ band. Eq. (1) shows that these flipping processes are closely connected with scattering processes of the Bloch electrons: a flipping of two spins effects a scattering of these electrons within the $`\stackrel{}{k}`$ space or, vice versa, a scattering of two Bloch electrons effects a flipping of the spins.

The mechanism of Cooper pair formation via crystal-spin-1 phonons in a narrow $`\sigma `$ band is characterized by the equation of constraint (3.36) of Part I. The physical origin of the related “constraining forces” may be understood as follows.

Assume all the $`N`$ electrons (with $`N`$ being even) in a narrow $`\sigma `$ band to form (symmetrized) Cooper pairs of the form

$$\beta _{\stackrel{}{k}m}^n=c_{\stackrel{}{k}m}^nc_{\stackrel{}{k}m}^nc_{\stackrel{}{k}m}^nc_{\stackrel{}{k}m}^n;$$
(4)
cf. Eq. (3.25) of Part I. Two of these Cooper pairs are depicted in Fig. 3(a). The total crystal spin of each Cooper pair is zero.

Let be $`𝒫^0`$ the subspace of the Hilbert space spanned by the $`N`$-electron states in the $`\sigma `$ band in which all the electrons form Cooper pairs. In $`𝒫^0`$ the electrons possess only $`2N`$ degrees of freedom because always two electrons are characterized by four quantum numbers $`k_x`$, $`k_y`$, $`k_z`$ and $`m`$ (cf. Section 4 of Part I). However, the electrons cannot stay in $`𝒫^0`$ because in reality they have $`4N`$ degrees of freedom. Thus, they necessarily become scattered out of this subspace creating then a phonon pair \[satisfying Eqs. (2) and (3)\] because any scattering process is connected with a flipping of the spins.

In Fig. 3(b) the two electrons with the wave vectors $`\stackrel{}{k}_1^{}`$ and $`\stackrel{}{k}_2^{}`$ become scattered emitting then a phonon pair with the wave vectors $`\stackrel{}{k}_1^{\prime \prime }`$ and $`\stackrel{}{k}_2^{\prime \prime }`$. The total crystal spin of the four electrons and the two phonons still is exactly zero (since $`H_{\mathrm{Cb}}^{n\sigma }`$ conserves the crystal spin angular momentum), but the phonons and electrons are strongly coupled in the sense that neither the two phonons nor the four electrons alone possess a defined total crystal spin. The strength of this coupling as well as the wave vectors of the scattered electrons and the emitted phonons strongly depend on the special features of the preceding spin-flipping process. As stated above, the spins flip in a narrow $`\sigma `$ band in order that the localized electrons comply with the Pauli principle. Hence, the features of the spin-flipping processes are determined by the wave vectors of the flipping electrons as well as by the wave vectors of all the other electrons that produce the potential in which the scattered electrons move.

The two phonons cannot be reabsorbed spontaneously because their complete absorption can only be the result of an appropriate motion of two other electrons in the $`\sigma `$ band. Such a motion is given when these two other electrons are forced by the Pauli principle to change their spin directions in such a way that the phonon pair ($`\stackrel{}{k}_1^{\prime \prime }`$, $`\stackrel{}{k}_2^{\prime \prime }`$) is completely reabsorbed. Such a motion is performed by the time-inverted electron pair ($`\stackrel{}{k}_1^{}`$, $`\stackrel{}{k}_2^{}`$): these two electrons carry out exactly the time inverted flipping process leading to a complete absorption of the phonon pair. Hence, after the absorption of the phonon pair the electrons again form Cooper pairs with zero total crystal spin, as illustrated in Fig. 3(c).

This mechanism may be interpreted in terms of “spring-mounted” Cooper pairs: There exist “springs” in the nonadiabatic system that push the electrons back into the subspace $`𝒫^0`$. These springs become tense when phonons get excited, and they become slacken again when these phonons are reabsorbed.

The operator $`H^0`$ represents a system with these springs being rigid. The effect of these rigid springs may be interpreted in terms of constraining forces that produce frozen Cooper pairs possessing half the degrees of freedom of unpaired electrons.

In analogy to classical physics, I suppose that only these quantum mechanical constraining forces produced by the described “springs” in narrow $`\sigma `$ bands can produce stable Cooper pairs (with half the degrees of freedom of unpaired electrons). This supposition is corroborated by the calculated band structures of those metals that experimentally prove to be superconductors; see Part II of this paper josm .

###### Acknowledgements.

I thank Ernst Helmut Brandt for critical comments on the manuscript.


text:
# Finite-Temperature Phase Diagram of Hard-Core Bosons in Two Dimensions

\[

## Abstract

We determine the finite-temperature phase diagram of the square lattice hard-core boson Hubbard model with nearest neighbor repulsion using quantum Monte Carlo simulations. This model is equivalent to an anisotropic spin-1/2 $`XXZ`$ model in a magnetic field. We present the rich phase diagram with a first order transition between solid and superfluid phase, instead of a previously conjectured supersolid and a tricritical end point to phase separation. Unusual reentrant behavior with ordering upon increasing the temperature is found, similar to the Pomeranchuk effect in 3He.

\] A nearly universal feature of strongly correlated systems is a phase transition between a correlation-induced insulating phase with localized charge carriers, to an itinerant phase. High temperature superconductors , manganites  and the controversial two dimensional (2D) “metal-insulator-transition”  are just a few examples of this phenomenon. The 2D hard-core boson Hubbard model provides the simplest example of such a transition from a correlation-induced charge density wave insulator near half filling to a superfluid (SF). It is a prototypical model for preformed Cooper pairs , of spin flops in anisotropic quantum magnets , of SF Helium films and of supersolids .

In simulations on this model, which does not suffer from the negative sign problem of fermionic simulations, we can investigate some of the pertinent questions about such phase transitions: what is the order of the quantum phase transitions in the ground state and the finite temperature phase transitions? Are there special points with dynamically enhanced symmetry ? Can there be coexistence of two types of order (such as a supersolid – coexisting solid and superfluid order)? Answers to these questions provide insight also for the other problems alluded to above.

The Hamiltonian of the hard-core boson Hubbard model we study is

$`H=`$ $``$ $`t{\displaystyle \underset{𝐢,𝐣}{}}(a_𝐢^{}a_𝐣+a_𝐣^{}a_𝐢)+V{\displaystyle \underset{𝐢,𝐣}{}}n_𝐢n_𝐣\mu {\displaystyle \underset{𝐢}{}}n_𝐢,`$ (1)
where $`a_𝐢^{}`$ $`(a_𝐢)`$ is the creation (annihilation) operator for hard-core bosons, $`n_𝐢=a_𝐢^{}a_𝐢`$ the number operator, $`V`$ the nearest neighbor Coulomb repulsion and $`\mu `$ the chemical potential. This model is equivalent to an anisotropic spin-1/2 $`XXZ`$ model with $`J_z=V`$ and $`|J_{xy}|=2t`$ in a magnetic field $`h=2V\mu `$. The zero field (and zero magnetization $`m_z=0`$) case of the spin model corresponds to the half filled bosonic model (density $`\rho =m_z=1/2`$) at $`\mu =2V`$. Throughout this Letter we will use the bosonic language, and refer to the corresponding quantities in the spin model where appropriate. Due to the absence of efficient Monte Carlo algorithms for classical magnets in a magnetic field there are still many open questions even in the classical version of this model, which was only studied by a local update method .

In Fig. 1 we show the ground-state phase diagram . For dominating chemical potential $`\mu `$ the system is in a band insulating state ($`\rho =0`$ and $`\rho =1`$ respectively), while it shows staggered checkerboard charge order ($`\rho =1/2`$) for dominating repulsion $`V`$. These solid phases are separated from each other by a SF. Earlier indications for a region of supersolid phase between checkerboard solid and SF phase turned out to be due to phase separation at this transition which is of first order at $`T=0`$ .

All of these phases extend to finite temperatures. On the strong repulsion side the hard-core boson Hubbard model is equivalent to an antiferromagnetic Ising model at $`t=0`$ and the insulating behavior extends up to a finite-temperature phase transition of the Ising universality class ($`T_c0.567V`$ at half filling). In case of zero repulsion the system is equivalent to an XY model in a magnetic field and the SF extends up to a transition of Kosterlitz Thouless (KT) type at finite temperature. The point $`V=2t`$ and $`\mu =2V`$ is also of special interest, it corresponds to the Heisenberg point in the spin model and has an enlarged $`O(3)`$ symmetry, with long range order only at zero temperature in 2D .

We focus our attention on the finite-temperature phase diagram in the vicinity of the first order quantum phase transition separating the SF and solid ground states away from the Heisenberg point. It is instructive to first consider the three dimensional (3D) version of this model, with a similar ground-state phase diagram. There, Fisher and Nelson  determined that the first order line terminates at a finite-temperature bicritical point at which the transition temperatures of the solid and SF phases meet. This bicritical point has a dynamically enhanced $`O(3)`$ symmetry – the symmetry of the Hamiltonian is only $`O(2)`$.

This feature alone dictates that despite the similar ground-state phase diagrams, the finite-temperature phase diagrams have to be very different in 2D, as an $`O(3)`$ critical point in 2D can exist only at zero temperature . Even with the absence of a supersolid in the ground state several alternative scenarios are still plausible: (i) the bicritical point could be pulled down to zero temperature, and the transition temperatures $`T_{KT}`$ and $`T_c`$ approach zero temperature without a direct phase transition between SF and solid at finite temperatures, (ii) a supersolid phase could exist at finite temperatures, although it is absent in the ground state, or (iii) the phase transition lines meet at finite temperature, but without dynamically enhanced symmetries.

In order to decide between these possibilities and determine the finite-temperature phase diagram we use quantum Monte Carlo simulations, employing a variant  of the worm algorithm  in the stochastic series expansion (SSE) representation . The simulations were performed in the grand canonical ensemble using lattices up to size $`96\times 96`$ and do not suffer from any systematic error apart from finite size effects. In contrast to the loop algorithm  the worm algorithm, while not as efficient at half filling (a spin model in zero magnetic field), remains effective away from half filling (in a finite field), where the loop algorithm slows down exponentially.

From now on we will restrict our discussion to $`V=3t`$, which corresponds to the dashed line in Fig. 1. We present the finite-temperature phase diagram for this representative value in Fig. 2. Due to particle-hole symmetry it is sufficient to show densities up to $`1/2`$ and chemical potentials up to $`\mu =2V=6t`$. We find that scenario (iii) is realized: the first order phase transition extends to a tricritical point (D) at around $`T_30.4t`$, where it meets with the second order melting transition of the solid. The KT transition temperature $`T_{KT}`$ of the SF meets the first order line at a lower temperature (B) – in contrast to the 3D case, where all the three lines meet at the same bicritical point (B$`=`$D) with enhanced $`O(3)`$ symmetry. We want to emphasize one unusual feature before discussing in detail how the phase diagram was determined: there is reentrant behavior in the vicinity of the tricritical point (D) and the system melts upon decreasing the temperature at constant chemical potential.

The first order phase transition is characterized by a finite jump $`\mathrm{\Delta }\rho `$ in the density. As discussed in Ref. , where we present the ground-state phase diagram, two peaks can be observed in the density distribution at the critical chemical potential. Since a two-peak structure can also be seen at a second order transition in a finite system, a careful finite size scaling analysis is needed. Figure 3a shows the finite size scaling of the densities of the two coexisting phases. For $`T0.39t`$ the jump $`\mathrm{\Delta }\rho `$ remains finite under extrapolation, thus confirming the first order nature up to this temperature (line from C to D in Fig. 2), while it scales to zero for $`T0.45t`$, an indication for a second order transition (line from D to E). Additionally we have checked the finite size scaling of the fourth order cumulant ratio $`C_4=1m_s^4/3m_s^2^2`$, of the staggered charge order parameter $`m_s`$ in the solid phase. Consistency with the expected size independent universal value of the Ising universality class $`C_4^c=0.6106900(1)`$  at the transition confirms the second order nature. We have thus determined that the first order transition extends up to a critical point (D) at a temperature larger than $`0.39t`$ and smaller than $`0.45t`$. In the next paragraph we show that the critical point where the KT line meets the first order line (point B) is below $`0.39t`$. Therefore the end point D of the first order line is a tricritical point.

The KT transition temperatures $`T_{KT}`$ (line from A to B in Fig. 2) can be determined using the SF number density $`\rho _s=(mT/2)W_x^2+W_y^2`$ . Here $`W_x`$ and $`W_y`$ are the winding numbers of the world lines in $`x`$ and $`y`$ direction respectively, and $`m=1/(2t)`$ is the mass of a boson. For system size $`L=\mathrm{}`$, $`\rho _s`$ has a universal value at $`T_{KT}`$, $`\rho _s(T_{KT})=(2/\pi )mT_{KT}`$, while it is always larger than this universal value for $`T<T_{KT}`$. At $`T=T_{KT}`$, $`\rho _s`$ jumps from the universal value to zero and we have $`\rho _s(T)=0`$ for all $`T>T_{KT}`$. However, for finite $`L`$, $`\rho _s`$ is nonzero at all $`T`$ and the finite size corrections of $`\rho _s`$ at $`T=T_{KT}`$ are given by  $`\rho _s\pi =2mT_{KT}(1+\{2\mathrm{log}[L/L_0(T_{KT})]\}^1)`$. Although these finite size corrections are known to be notoriously difficult, we can obtain strict upper and lower bounds for $`T_{KT}`$ by plotting $`1/(\pi \rho _s/mT2)\mathrm{log}(L)`$ as function of $`\mathrm{log}L`$. As $`L`$ is increased, this quantity converges to the constant $`(\mathrm{log}L_0)`$ at $`T_{KT}`$ and diverges to $`\pm \mathrm{}`$ for $`T>T_{KT}`$ and $`T<T_{KT}`$, respectively . The error bars we obtain using this method are reliable and small enough for our purposes.

To investigate how the KT transition line meets the phase separation line we perform canonical measurements of the SF density $`\rho _s`$ at the density $`\rho _f`$ of the fluid on this coexistence line. In Fig. 3b we show the finite size scaling of $`\rho _s(\rho _f)`$. At low temperatures $`T=0.3t`$ $`\rho _s`$ of the SF at the coexistence line is well above the universal value for a KT transition. Therefore the transition between SF and solid at temperatures $`T0.3t`$ does not include a KT transition. The system remains SF right up to a direct first order phase transition from the SF to the solid. At $`T=0.39t`$ on the other hand, $`\rho _s`$ as a function of the system size $`L`$ drops below the universal value and should scale to zero according to . There is thus a first order transition from solid to normal fluid at $`T=0.39t`$ and the KT line has to meet the first order line at a finite temperature between $`0.3t`$ and $`0.39t`$.

As can be seen in Fig. 2, the checkerboard solid is melting in a first order transition to a SF upon decreasing the temperature around $`\mu 2.35t`$. Even more special is the melting to a normal fluid in the vicinity of the tricritical point (D) both upon increasing or decreasing the temperature. To proof this reentrant behavior, we first determine the critical chemical potential $`\mu _c(T=0.39t)`$ of the first order transition at $`T=0.39t`$. We then show that $`\mu _c(T=0.5t)<\mu _c(T=0.39t)`$ On a finite system we define $`\mu _c`$ as the chemical potential where there is an equal probability to be in the solid and fluid phase (equal peak areas in the density probabilities). Performing a finite size extrapolation (Fig. 4a) we obtain $`\mu _c(0.39t)=2.3015(10)t`$. Next we calculate the upper bound for the critical chemical potential $`\mu _c(T)`$ at $`T=0.5t`$. Observing that the fourth order cumulant ratio $`C_4`$ scales above the universal critical value $`C_4^c`$ at $`\mu =2.295t`$ (Fig. 4b), we thus have $`\mu _c(T=0.5t)<2.295t`$. As $`\mu _c(T=0.5t)<\mu _c(T=0.39t)`$ we have pinpointed reentrant behavior where the solid melts to a normal fluid both when increasing and decreasing the temperature.

To summarize, we have investigated the finite-temperature phase diagram of the hard-core boson Hubbard model with nearest neighbor repulsion using quantum Monte Carlo simulations. The first order phase transition in the ground state between SF and solid phases extends to finite temperatures and turns into a second order melting line of the solid phase at a finite-temperature tricritical point. As the KT line meets the first order line at a lower temperature, we have a qualitatively different phase diagram than in 3D , where there is a bicritical point with dynamically enhanced $`O(3)`$ symmetry.

We expect this phase diagram, obtained for a representative value of $`V=3t`$, to be valid also for other couplings $`V>2t`$. On approaching the Heisenberg point $`V=2t`$ the temperature scale of all the transitions drops to zero. In the other limit $`t/V0`$, we approach the antiferromagnetic Ising model. The KT transition temperatures $`T_{KT}`$ and the tricritical point $`T_3`$ should scale with $`t`$, since in the Ising model the phase transition is first order at zero temperature, but second order at any finite temperature , consistent with $`T_3/Vt/V0`$.

The melting line of the solid shows several unusual features. As a function of chemical potential there are three different phase transitions: at low temperatures the melting transition is first order into a SF. At intermediate temperatures the melting transition is still first order, but to a normal fluid, while it is a second order transition to a normal fluid at even higher temperatures.

There is reentrant behavior in the vicinity of the tricritical point: as a function of increasing temperature there are three phase transitions: first the SF turns into a normal fluid. Upon increasing the temperature further this normal fluid solidifies, and finally melts again. The reason for this reentrant behavior is conjectured to be a large entropy of the weakly interacting vacancies (or interstitials at $`\rho >1/2`$) in the solid, larger than the entropy of the normal fluid. Such an ordering transition upon increasing temperature is rare but has been observed before, e.g. in the bcc Ising antiferromagnet .

Finally, we want to mention the striking qualitative similarity of the presented phase diagram to those of fermionic 3He and bosonic 4He in two dimensions. In 3He there is also a reentrant melting line (the Pomeranchuk effect ) and a SF phase at low temperatures . In 3He however the large entropy of the solid is due to the nuclear spin degrees of freedom and the effect even more pronounced. In 4He films, the solid phase melts to a SF in a small range for the pressure upon decreasing the temperature , like in our phase diagram where there is a direct transition form SF to solid for constant values of $`\mu `$ around $`\mu 2.35t`$. Note that in our phase diagram we have pinned down the even more subtle effect of a reentrant melting line between solid and normal fluid.

We wish to thank G.G. Batrouni and D.M. Ceperley for stimulating discussions. The simulations were performed on the Asgard Beowulf cluster at ETH Zürich, employing the Alea library for Monte Carlo simulations . G.S., S.T. and M.T. acknowledge support of the Swiss National Science Foundation.