format: text

text:
# THE COSMIC RADIO AND INFRARED BACKGROUNDS CONNECTION

## 1 INTRODUCTION

Star forming galaxies exhibit a remarkable correlation between their radio and infrared (IR) fluxes covering over four decades of IR flux intensities (Helou, Soifer, & Rowan–Robinson 1985, and references therein). The correlation results from the fact that the radio and infrared fluxes are different manifestations of the physical processes associated with the lifecycle of the same stellar objects (Helou & Bicay 1993, Lisenfeld, Völk, & Xu 1996). Massive stars heat the dust in the interstellar medium, form H II regions, and after their explosive deaths accelerate particles to cosmic ray energies, processes that, respectively, give rise to the observed galactic thermal infrared, radio thermal and synchrotron emission. The correlation may not hold for active galactic nuclei (AGN), in which the radio emission is not associated with the life cycle of massive stars.

If the radio–IR correlation for individual star forming galaxies in the local universe holds for all redshifts, then their cumulative contributions to the cosmic radio and infrared backgrounds should be related. Haarsma & Partridge (1998, hereafter HP98) used this correlation to estimate the intensity of the cosmic radio background (CRB) that can be attributed to star–forming galaxies. At the time of their analysis only the CIB at wavelengths larger than $``$ 120 $`\mu `$m was definitively detected by the Diffuse Infrared Background Experiment (DIRBE) and Far Infrared Spectrophotometer (FIRAS) instruments on board the Cosmic Background Explorer (COBE) satellite (Hauser et al. 1998, Fixsen et al. 1998). For simplicity, HP98 assumed that all the sources giving rise to the CIB release their energy instantaneously at redshift $`z`$ = 1. With this simple assumption about the formation epoch of the CIB, most of its detected intensity is determined by the $``$ 80 $`\mu `$m spectral luminosity density of galaxies at the assumed redshift of its creation. Adopting a linear relation, $`S_\nu (20cm)S_\nu (80\mu `$m), between the radio and IR fluxes from star–forming galaxies, HP98 derived their contribution to the 40 cm CRB intensity. To compare this CRB estimate to the observed Bridle (1967) datum point at 178 MHz, HP98 adopted a $`\nu ^{0.7}`$ power law spectrum for the radio sources. They derived a brightness temperature of $`T_{\mathrm{crb}}`$(178 MHz) $``$ 15 K, about half the value inferred from the observations. HP98 attributed the remaining CRB intensity to the contribution of AGN which do not contribute significantly to the CIB (Barger et al. 2001, Fadda et al. 2001).

In this paper we extend the analysis of HP98 to the more general case in which the cosmic star formation rate (CSFR) of the CIB sources is a general function of redshift. First, in §2, we briefly review the various presentations of the radio–IR correlation and rederive this correlation in terms of the galaxies’ 8 – 1000 $`\mu `$m IR fluxes. Previous presentations expressed this correlation as a function of only the far–IR (FIR) fluxes, derived from the galaxies 60 and 100 $`\mu `$m detections by the Infrared Astronomy Satellite (IRAS). Such a correlation may systematically underestimate the radio flux from more luminous star forming galaxies in which an increasing fraction of the IR flux may be emitted in the shorter (12 and 25 $`\mu `$m) IRAS bands. We then explore the functional form of the radio–IR correlation. HP98 adopted a linear relation between the radio and IR fluxes, which may not be an accurate fit to the data. Significant deviations from linearity will require detailed knowledge of the evolution of the galaxies’ IR luminosity function with redshift. However, we argue that the data are consistent with a linear relation over most of the range of the IR luminosities that contribute to the local IR luminosity density. In §3 we first derive the expressions for the CIB and CRB intensities produced by star–forming galaxies for a general power–law correlation between their radio and IR fluxes. We then derive a simple analytic expression relating these background intensities when their radio–IR correlation is linear. In §4 we apply these results to different cosmic star formation histories with their associated CIB intensity limits. A brief summary of our paper is presented in §5.

## 2 THE RADIO–INFRARED FLUX CORRELATION

The radio–IR correlation can be written in a general form as:

$$P_\nu (\nu ,\mathrm{L}_{\mathrm{IR}})(\mathrm{W}\mathrm{Hz}^1)=\kappa (\nu )\left(\frac{\mathrm{L}_{\mathrm{IR}}}{\mathrm{L}_{}}\right)^{\beta (\nu )}$$
(1)
where LIR is the IR luminosity of the galaxy, $`P_\nu (\nu )`$ is the specific radio luminosity at the radio frequency $`\nu `$, and $`\kappa `$ and $`\beta `$ are coefficients that depend on the frequency at which the correlation is expressed. The coefficient $`\kappa `$ is dimensional, and its value depends on the units used in the correlation. The presentation above is often referred to as the {Pν–LIR} presentation of the radio–IR correlation. An alternative is the {Sν–FIR} presentation, in which the radio and IR outputs are expressed in units of W m-2 Hz-1 and W m-2, respectively.

The two presentations are, of course, equivalent for individual galaxies. However, for determining $`\beta `$, the choice of presentation is important if $`\beta `$ is a priori known to be different from unity (Cox et al. 1988). If all galaxies lie on a line of slope $`\beta `$ = 1 in a log–log Sν vs. FIR plot, then the conversion to a {Pν–LIR} presentation will move them up or down along a diagonal line, leaving the overall slope of the correlation unchanged. However, if the slope of the correlation is different from unity in a {Sν–FIR} presentation, then fainter galaxies, which tend to be more distant than brighter ones, may move systematically away or toward the $`\beta `$ = 1 slope line in a {Pν–LIR} presentation. The overal effect will be to flatten the slope if $`\beta <`$ 1, and to steepen it if $`\beta >`$ 1 in a transition from a {Sν–FIR} to a {Pν–LIR} presentation.

Various investigators examined the radio–IR correlation using a variety of sample selection criteria. All use FIR fluxes measured by the IRAS satellite and construct galaxy FIR luminosities (fluxes) from the 60 and 100 $`\mu `$m detections using the following relation (Sanders & Mirabel 1996):

$`F_{\mathrm{FIR}}(\mathrm{Wm}^2)`$ $`=`$ $`1.26\times 10^{14}[2.58S_\nu (60\mu m)+S_\nu (100\mu m)]`$ (2)
$`\mathrm{L}_{\mathrm{FIR}}(\mathrm{L}_{})`$ $`=`$ $`1.26\times 10^{12}[2.58\mathrm{L}_\nu (60\mu m)+\mathrm{L}_\nu (100\mu m)]`$
where $`S_\nu `$ and Lν are flux and luminosity densities expressed in units of Jy ( = 10-26 W m-2 Hz-1) and L/Hz, respectively.

As we will show in §3, the value of $`\beta `$ plays an important role in quantifying the CIB–CRB connection. A value of $`\beta `$ that is significantly different from unity will require detailed knowledge of the evolution of the IR luminosity function with redshift, whereas a value of unity will considerably simplify the relation between the two background emissions. We therefore briefly review previous determinations of $`\beta `$.

Initial investigations, using a relatively small sample of 91 radio–selected galaxies observed at $`\nu `$ = 4.8 GHz (de Jong et al. 1985) and 38 optically– and radio– selected galaxies observed at 1.4 GHz (Helou, Soifer, & Rowan–Robinson 1985), found the radio–IR correlation to be linear. Wunderlich & Klein (1988) extended the analysis of de Jong et al. (1985) to a wider range of galaxy types, covering about 4 orders of magnitudes in FIR luminosities. They found that the correlation at $`\nu `$ = 4.8 GHz is approximately linear ($`\beta =0.99\pm 0.38`$) up to IR luminosities of $`9.1\times 10^9`$ L, with $`\beta `$ increasing to a value of $`1.26\pm 0.33`$ at higher luminosities. However, we found that their sample of galaxies could also be fitted with a single power law with a slope of $`\beta `$ = 1.05$`\pm `$0.03.

Devereux and Eales (1989) correlated the 1.49 GHz power with the FIR output from an optically selected sample of galaxies with luminosities between $`10^810^{11}`$ L. They found the slope of the correlation to be significantly larger than unity with a value of $`\beta `$ = 1.28. A similar conclusion was reached by Cox et al. (1988) who studied the correlation for a flux limited sample of 74 radio galaxies at 151 MHz. They found $`\beta `$ = 1.32$`\pm `$0.06 in a {Sν–FIR} presentation, and $`\beta `$ = 1.15$`\pm `$0.04 in a {Pν–LIR} presentation of the data, and adopted an average slope of $`\beta `$ = 1.23.

Chi & Wolfendale (1990) provided theoretical arguments why the slope of the radio–IR correlation should not be unity. They predicted a break in the slope of the correlation at an IR luminosity above which most of the electrons producing the radio synchrotron emission remained trapped in the galaxy. For these galaxies the slope of the correlation should be unity, but for the smaller and less IR luminous galaxies they argued that $`\beta `$ should be larger than unity. They claimed to have found evidence for this effect in the data and reported values of $`\beta `$ = 1.37 for low luminosity galaxies, a trend exactly opposite to that found by Wunderlich & Klein (1988). Our analysis of the Chi & Wolfendale sample did not confirm either trend. We found that their data could be represented by a single slope with a value of $`\beta =1.27\pm 0.06`$.

Condon, Anderson, & Helou (1991; hereafter CAH91) reexamined the radio–IR correlation using the IRAS revised bright galaxies sample (BSG) from which spectroscopically identified AGN were subtracted, and which were detected with the VLA at 1.49 GHz. They found the slope of the correlation to be significantly greater than unity with a smaller value of $`\beta =1.11\pm 0.02`$ over the $`10^9`$ to 1013 L luminosity range.

For the purpose of the present analysis we reexamined the radio–IR correlation starting from the same galaxy sample used by CAH91, the revised IRAS Bright Galaxy Sample, elimating galaxies that were not detected in all 4 IRAS bands, and removing AGN using a more recent catalog of spectroscopically confirmed AGN galaxies (Veron-Setty & Veron 2001, catalog of Quasars and AGN). The total number of galaxies left in the sample is 222, compared to 258 galaxies used in the analysis of CAH91. We used the additional 12 and 25 $`\mu `$m IRAS bands to derive the total 8–1000 $`\mu `$m IR luminosity from the relation (Sanders & Mirabel 1996):

$$\mathrm{L}_{\mathrm{IR}}(\mathrm{L}_{})=1.8\times 10^{12}[13.48\mathrm{L}_\nu (12\mu m)+5.16\mathrm{L}_\nu (25\mu m)+2.58\mathrm{L}_\nu (60\mu m)+\mathrm{L}_\nu (100\mu m)]$$
(3)
Figure 1 depicts the radio–IR correlation for a {Pν–LIR} presentation of the data, and the best power-law fit to the data. We derive an almost identical slope $`\beta `$ for the same sample of galaxies as CAH91, with a smaller value for $`\kappa `$, reflecting the fact that the 8 – 40 $`\mu `$m emission was not included in their galactic IR energy budget. The values of the fit for radio–IR correlation are: {$`\kappa `$$`\beta `$} = {2.96$`\times 10^{10}`$$`1.083\pm 0.02`$}. An almost identical slope was obtained for the {Sν–FIR} presentation, suggesting that there were no significant flux related systematic errors in determining the distances of the galaxies in the sample. Also shown in Figure 1 is a forced linear fit to the data with {$`\kappa `$, $`\beta `$} = {2.47$`\times 10^{11}`$, 1.0}. The linear fit is almost indistinguishable from the best fit for galaxies with luminosities above $``$ 1010 L, which produce most of the IR luminosity density in the local and high–$`z`$ universe. The same linear fit will provide a similarly good representation of the radio–IR correlation even if the AGN were included in the data. The AGN population however shows more dispersion around the fit.

## 3 THE CIB and CRB INTENSITIES PRODUCED BY STAR–FORMING GALAXIES

The specific intensity $`I_\nu ^{\mathrm{crb}}(\nu _R)`$ of the CRB at the observed radio frequency $`\nu _R`$ is given by the integral over sources (e.g. Peebles 1993):

$$I_\nu ^{\mathrm{crb}}(\nu _R)=\left(\frac{c}{4\pi }\right)_0^{\mathrm{}}_\nu (\nu ,z)\left|\frac{dt}{dz}\right|𝑑z$$
(4)
where $`_\nu (\nu ,z`$) is the specific radio luminosity per comoving volume element at redshift $`z`$, $`\nu =\nu _R(1+z)`$ is the frequency in the rest frame of the luminous objects and $`\left|dt/dz\right|`$ is given by (Hogg 1999)

$$\left|dt/dz\right|^1=H_0(1+z)\left[(1+z)^2(1+\mathrm{\Omega }_mz)z(2+z)\mathrm{\Omega }_\mathrm{\Lambda }\right]^{1/2}$$
(5)
where $`H_0`$ is the Hubble constant, $`h`$ its value in units of 100 km s-1 Mpc-1, $`\mathrm{\Omega }_m\rho _m/\rho _c`$ is the present mass density of the universe normalized to the critical density $`\rho _c=1.88\times 10^{29}h^2`$ g cm-3, and $`\mathrm{\Omega }_\mathrm{\Lambda }\mathrm{\Lambda }/3H_0^2`$ is the dimensionless cosmological constant.

We will assume that the radio–IR correlation, expressed at a given rest frame frequency $`\nu _0`$, holds for all redshifts and that the radio spectrum of the individual galaxies follows a $`\nu ^\alpha `$ power law with spectral index $`\alpha `$, so that $`P_\nu (\nu ,\mathrm{L}_{\mathrm{IR}},z)=(\nu /\nu _0)^\alpha P_\nu (\nu _0,\mathrm{L}_{\mathrm{IR}},z)`$. The spectral luminosity density $`_\nu `$ can then be expressed as:

$`_\nu (\nu ,z)`$ $`=`$ $`{\displaystyle _{\mathrm{L}_{\mathrm{min}}}^{\mathrm{L}_{\mathrm{max}}}}P_\nu (\nu ,\mathrm{L}_{\mathrm{IR}},z)\mathrm{\Phi }(\mathrm{L}_{\mathrm{IR}},z)\mathrm{dL}_{\mathrm{IR}}`$
$`=`$ $`\kappa _0\left({\displaystyle \frac{\nu }{\nu _0}}\right)^\alpha {\displaystyle _{\mathrm{L}_{\mathrm{min}}}^{\mathrm{L}_{\mathrm{max}}}}\mathrm{L}_{\mathrm{IR}}^{}{}_{}{}^{\beta _0}\mathrm{\Phi }(\mathrm{L}_{\mathrm{IR}},z)\mathrm{dL}_{\mathrm{IR}}`$
where $`\mathrm{\Phi }(\mathrm{L}_{\mathrm{IR}},z)[\mathrm{d}\phi (\mathrm{L}_{\mathrm{IR}},z)/\mathrm{dL}_{\mathrm{IR}}]`$ is the differential IR luminosity function, d$`\phi (\mathrm{L}_{\mathrm{IR}},z)`$ is the comoving number density of galaxies with luminosities between LIR and LIR \+ dLIR in the {$`\mathrm{L}_{\mathrm{min}},\mathrm{L}_{\mathrm{max}}`$} luminosity interval, and $`\kappa _0`$, $`\beta _0`$ are the parameters of the radio–IR correlation at frequency $`\nu _0`$.

The specific intensity of the CRB can now be written as:

$$I_\nu ^{\mathrm{crb}}(\nu __R)=\left(\frac{c}{4\pi }\right)\left(\frac{\nu __R}{\nu _0}\right)^\alpha _0^{\mathrm{}}\frac{f(z)_{\mathrm{IR}}(z)}{(1+z)^\alpha }\left|\frac{dt}{dz}\right|𝑑z$$
(7)
where $`_{\mathrm{IR}}`$(z) is the comoving IR luminosity density at redshift $`z`$:

$$_{\mathrm{IR}}(z)_{\mathrm{L}_{\mathrm{min}}}^{\mathrm{L}_{\mathrm{max}}}\mathrm{L}_{\mathrm{IR}}\varphi (\mathrm{L}_{\mathrm{IR}},z)\mathrm{dL}_{\mathrm{IR}}$$
(8)
and $`f(z)`$ is defined as:

$$\frac{f(z)\kappa _0_{\mathrm{L}_{\mathrm{min}}}^{\mathrm{L}_{\mathrm{max}}}\mathrm{L}_{\mathrm{IR}}^{}{}_{}{}^{\beta _0}\varphi (\mathrm{L}_{\mathrm{IR}},z)\mathrm{dL}_{\mathrm{IR}}}{_{\mathrm{IR}}(z).}$$
(9)
The same star forming galaxies give rise to the CIB with an intensity $`I_{\mathrm{CIB}}`$ that is given by:

$$I_{\mathrm{CIB}}=\left(\frac{c}{4\pi }\right)_0^{\mathrm{}}\frac{_{\mathrm{IR}}(z)}{1+z}\left|\frac{dt}{dz}\right|𝑑z$$
(10)
For a linear correlation, $`\beta _0`$ = 1, and $`f(z)=\kappa _0`$(Hz-1). The spectral intensity of the CRB then becomes:

$$I_\nu ^{\mathrm{crb}}(\nu __R)=\left(\frac{c}{4\pi }\right)\left(\frac{\nu __R}{\nu _0}\right)^\alpha \kappa _0_0^{\mathrm{}}\frac{_{\mathrm{IR}}(z)}{(1+z)^\alpha }\left|\frac{dt}{dz}\right|𝑑z$$
(11)
In this case the relation between $`I_{\nu _R}^{\mathrm{crb}}`$ and $`I_{\mathrm{CIB}}`$ is considerably simplified and given by:

$`I_{\mathrm{CIB}}`$ $`=`$ $`\left({\displaystyle \frac{\nu _R}{\nu _0}}\right)^\alpha {\displaystyle \frac{g(\alpha )}{\kappa _0}}I_\nu ^{\mathrm{crb}}(\nu __R)`$ (12)
$`=`$ $`1.18\times 10^8\left({\displaystyle \frac{\nu _R}{\nu _0}}\right)^\alpha {\displaystyle \frac{g(\alpha )}{\kappa _0}}\nu _R^2(\mathrm{MHz})\mathrm{T}_{\mathrm{crb}}(\nu _\mathrm{R})\mathrm{nW}\mathrm{m}^2\mathrm{sr}^1`$
where

$$g(\alpha )\left\{_0^{\mathrm{}}\frac{_{\mathrm{IR}}(z)}{1+z}\right|\frac{dt}{dz}|dz\}/\left\{_0^{\mathrm{}}\frac{_{\mathrm{IR}}(z)}{(1+z)^\alpha }\right|\frac{dt}{dz}|dz\},$$
(13)
and $`T_{\mathrm{crb}}(\nu _R)`$ is the CRB brightness temperature at frequency $`\nu _R`$.

## 4 THE CIB–CRB CONNECTION FOR VARIOUS COSMIC STAR FORMATION HISTORIES

Several models have been recently developed for calculating the evolution of $`_{\mathrm{IR}}(z)`$ with redshift. Here we will focus on some select representative cases: Malkan & Stecker (2001, baseline model), Chary & Elbaz (2001, pure luminosity evolution model), and Xu et al. (2001, peak model). These models were developed to explain the limits and detections of the CIB spectrum and galaxy number counts obtained with the IRAS, the Infrared Space Observatory (ISO), and the Submillimeter Common User Bolometric Array (SCUBA) instrument on the James Clerk Maxwell Telescope (JCMT) at various IR and submillimeter wavelengths. Figure 2 presents the evolution of $`_{\mathrm{IR}}(z)`$ for these select models. Also plotted are observational estimates of the comoving cosmic star formation rate (CSFR) at different redshifts. Star formation rates were converted to IR luminosity densities using the relation (Kennicutt 1998): $`_{\mathrm{IR}}`$(L Mpc-3) = $`6\times 10^9\times \rho _{}`$(M yr-1 Mpc-3). References to the CSFR data can be found in the papers listed above. Also shown in the figure are approximate upper and lower limits to the IR luminosity density. The top–hat function centered on redshift $`z`$ = 1 in the figure represents the instantaneous energy injection model used by Haarsma & Partridge (1998) to calculate the CRB. The amplitude of the function was chosen so that it reproduced the nominal $``$ 5–1000 $`\mu `$m CIB intensity of 50 nW m-2 sr-1 (Hauser & Dwek 2001).

Figure 3 shows the value of $`g(\alpha )`$ as a function of $`\alpha `$ for the different star formation histories depicted in Figure 2. The function $`g(\alpha )`$ is insensitive to the cosmic star formation history, a direct consequence of the fact that the integrands in eq. (13) differ only by a factor of $`(1+z)^{\alpha 1}`$. For $`\alpha `$ 0.6–1.2, the most probable range of values for the radio spectral index, $`g(\alpha )`$ is well approximated by:

$$g(\alpha )=\alpha ^{0.688}\mathrm{for}\alpha 0.61.2$$
(14)
to an accuracy better than 2%.

Using eqs. (12), (14) and the linear radio–IR correlation at $`\nu _0`$ = 1.49 GHz: $`P_\nu (\mathrm{W}\mathrm{Hz}^1)=2.47\times 10^{11}\mathrm{L}_{\mathrm{IR}}`$(L), we can express the CIB intensity in terms of the radio brightness temperature at frequency $`\nu _R`$ = 178 MHz as:

$`I_{\mathrm{CIB}}(\mathrm{nWm}^2\mathrm{sr}^1)`$ $`=`$ $`𝒜(\alpha )T_{\mathrm{crb}}`$
$`\mathrm{where}𝒜(\alpha )`$ $`=`$ $`\{3.0,2.7,2.4,2.1\}\mathrm{for}\alpha =\{0.6,0.7,0.8,0.9\}`$ (15)
Table 1 lists the values of the CIB intensity and the 178 MHz brightness temperature obtained from the respective use of eqs. (10) and (11) for the different star formation histories and radio spectral index $`\alpha `$. The entries in the table satisfy the analytical approximation for the relation between $`I_{\mathrm{CIB}}`$, $`T_{\mathrm{crb}}`$, and $`\alpha `$ to an accuracy of a few percent. For comparison we also listed the observed limits and detections of the CIB (Hauser & Dwek 2001) and the 178 MHz CRB brightness temperature inferred from the Bridle data (Bridle 1967) for the different values of $`\alpha `$. Also shown in the table are the observed $`T_{\mathrm{crb}}/I_{\mathrm{CIB}}`$ ratios and those predicted by the CIB-CRB correlation \[eq. (15)\].

It is interesting to compare our model predictions with the result obtained by Haarsma & Partridge (1998). HP98 found that the star–forming galaxies that produce the detected 120–260 $`\mu `$m CIB intensity of $``$ 22 nW m-2 sr-1 contribute about 15 K to the 178 MHz brightness temperature (they adopted a value of $`\alpha `$ = 0.7). This result seems at first glance to be in disagreement with the entry for the HP98 model in the table. The reason for this apparent ”discrepancy” is that HP98 expressed the radio–IR correlation in terms of the FIR luminosity of galaxies. Had they expressed this correlation in terms of their IR luminosity, which is about twice the FIR value, they would have derived a radio brightness temperature of $``$ 17 K for a 3.5–1000 $`\mu `$m background intensity of 50 nW m-2 sr-1, almost identical to the value listed in the table. The analytical expression presented in this paper reproduces the results derived by HP98 for their specific CIB production scenario, and generalizes their treatment to any cosmic star formation history.

## 5 DISCUSSION AND SUMMARY

In this paper we derived an analytical expression for the correlation between the CIB and CRB intensities as a function of $`\alpha `$, the spectral index of the radio sources \[see eq. (15)\]. This correlation is summarized for select galaxy number counts and CIB models in Table 1. Several conclusions can be drawn from a simple examination of the table: (1) the minimal CSFR,defined by the lower envelope of the data, is definitely ruled out as a viable representation of the CSFR since it falls short of providing the observed intensity of the CIB. The maximal CSFR,defined by the upper envelope of the data, is consistent with the CIB and CRB limits only for values of $`\alpha `$ 0.8; (2) the baseline model of Malkan & Stecker is barely consistent with the lower limit on the CIB intensity and consequently, for $`\alpha `$ = 0.7, it requires AGN to contribute more than 70% of the 178 MHz background; (3) for $`\alpha `$ = 0.9, all models, with the exception of those of Malkan & Stecker and Chary & Elbaz, predict CRB temperatures that are 2$`\sigma `$ above the observed value. The latter two models leave no room for a significant contribution of AGN to the CRB; (4) for values of $`\alpha `$ 0.6–0.8, the calculated $`T_{\mathrm{crb}}/I_{\mathrm{CIB}}`$ ratio is lower than the nominal observed ratios, suggesting that $``$60 to 20% of the CRB at 178 MHz must arise from AGN, regardles of the radiative history of star–forming galaxies.

All conclusions listed above assume that the radio–IR correlation observed in the local universe can extended to galaxies at all redshifts. Future observations will provide a larger sample of galaxies with which to study the radio–IR correlation and place tighter limits on the CIB, advances that will lead to a better understanding of the relative contribution of star–forming galaxies and AGN to the CRB.

Acknowledgement: We thank Rick Arendt for useful discussions and his insightful comments on an earlier version of the manuscript. ED acknowledges NASA’s Astrophysics Theory Program NRA 99-OSS-01 for support of this work. MB’s summer student internship at NASA/GSFC was supported by the ”Research Opportunities for Undergraduates in the Laboratory for Astronomy and Solar Physics” program.


text:
# 1 The potentials 𝑉_𝑝(𝑥) as a function of 𝑥 for 𝑝={0,1,2}. Note the asymptotic linear behavior discussed in the text.

The study of one-dimensional fermionic systems has recently attracted much attention due to progress in manufacturing of mesoscopic devices. Although much progress has been made in the theoretical understanding of the dynamics of these strongly correlated systems, many interesting questions remain still unanswered. One theoretical tool that is successfully used is the bosonization of the fermionic degrees of freedom : for a given fermionic system in the thermodynamic limit, one considers the low-lying particle-hole excitations around the Fermi points. Linearization of the dispersion relation around these leads to a relativistic conformal field theory  of these gapless effective degrees of freedom. More recently, an extended bosonization procedure incorporating the subsequent non-linear corrections of the actual dispersion relation has been proposed . The algebraic bosonization technique is an all-order bosonization approach based on the $`W_{1+\mathrm{}}`$ algebra, a characteristic fermionic symmetry  extending the conformal symmetry of the standard bosonization. The main advantage of this method is that it yields a systematic expansion of the full hamiltonian (including interactions) and all physical observables in powers of the momentum fluctuations around the Fermi points, involving algebraic procedures only.

In this letter, we present a one-parameter family of hermitian one-matrix models which originate in the $`W_{1+\mathrm{}}`$ algebra and are intimately connected to the algebraic bosonization of one-dimensional fermionic systems. Consider the Fock space of a $`(1+1)`$-dimensional relativistic Weyl fermionic field with creation and annihilation operators $`a_n^{}`$ and $`a_n`$, respectively ($`n`$ is an integer that labels the momentum), satisfying $`\{a_n^{},a_m\}=\delta _{n,m}`$. We define the $`W_{1+\mathrm{}}`$ operators $`V_n^k`$ ($`k=0,1,2,\mathrm{}`$) in this specific basis :

$$V_n^k=\underset{r𝐙}{}P_n^k(r):a_{rn}^{}a_r:,$$
(1)
which satisfy the $`W_{1+\mathrm{}}`$ algebra ,

$$[V_n^i,V_m^j]=(jnim)V_{n+m}^{i+j1}+q(i,j,n,m)V_{n+m}^{i+j3}+\mathrm{}+\delta ^{ij}\delta _{n+m,0}cd(i,n),$$
(2)
where the structure constants $`q(i,j,n,m)`$ and $`d(i,n)`$ are polynomial in their arguments, $`c`$ is the central charge ($`c=1`$ for the single Weyl fermion), and the dots stand for a finite number of terms involving the operators $`V_{n+m}^{i+j12k}`$, with $`k=0,1,\mathrm{},[(i+j1)/2]`$, where the brackets denote the integer part function (the complete expression of (2) is given in ). The meaning of the operators $`V_n^k`$ is as follows: they parametrize particle-hole excitations near the Fermi points, such that the index $`n`$ accounts for the fixed momentum transfer, whereas the index $`k`$ denotes the geometric moments (i.e.,“multipole moments”) of the deformation of the Fermi-Dirac sea . The explicit expression of the polynomials $`P_k^n(x)`$ can be inferred directly from the $`W_{1+\mathrm{}}`$ algebra :

$$P_n^k(x)=\frac{(1)^k}{\left(\begin{array}{c}2k\\ k\end{array}\right)}\underset{j=0}{\overset{n}{}}(1)^j\left(\begin{array}{c}n\\ j\end{array}\right)^2\left(n+1x\right)_{kj}\left(x\right)_j,$$
(3)
where $`(a)_n=_{\mathrm{}=0}^{n1}(a+\mathrm{})`$ denotes the Pochhammer symbol. Some properties can be deduced from (3): i) the polynomials have a definite parity with respect to the reflection about $`x=1/2`$: even $`n`$ (resp. odd) correspond to even (odd) polynomials. ii) for $`n`$ large enough, the zeroes of $`P_n^k(x)`$ are located on the line $`x=1/2+i\zeta `$, with $`\zeta `$ real, and for pure imaginary values of $`n`$, all zeroes fall along this line. For $`k`$ even, all zeroes come in complex conjugate pairs, and for $`k`$ odd the same is true except for a single zero at $`x=1/2`$. However, the most important consequence of (3) is that the $`P_n^k(x)`$ satisfy a three-term recurrence relation:

$$xP_n^k(x)=P_n^{k+1}(x)+\frac{1}{2}P_n^k(x)+\frac{k^2(n^2k^2)}{4(4k^21)}P_n^{k1}(x).$$
(4)
This property has the important consequence of leading to a new class of random matrix models associated to the $`W_{1+\mathrm{}}`$ algebra, as we shall show in the following.

Matrix models constitute another important tool when studying systems displaying universal behavior in the thermodynamic limit . There are many know examples in which the symmetries that characterize the relevant effective degrees of freedom of a given system are simply encoded in a specific random matrix model, which posses the advantage of being a simpler $`(0+1)`$-dimensional theory. Consider a generic random one-matrix model, defined by its partition function:

$$Z(N)=d^{N^2}\mathrm{\Phi }\mathrm{exp}\left[\mathrm{Tr}V(\mathrm{\Phi })\right],$$
(5)
where $`\mathrm{\Phi }`$ is an $`N\times N`$ Hermitian matrix, $`\mathrm{Tr}`$ is the trace, ane the real function $`V(\mathrm{\Phi })`$ is the potential. Eventually, we shall be interested in the limit $`N\mathrm{}`$, which is the relevant thermodynamic limit. The matrix $`\mathrm{\Phi }`$ has real eigenvalues $`x_i`$, $`i=1,\mathrm{},N`$, in terms of which the partition function can be written as :

$$Z(N)=\mathrm{\Omega }(N)\mathrm{}\underset{i=1}{\overset{N}{}}dx_i\mathrm{\Delta }_{N}^{}{}_{}{}^{2}(x)\mathrm{exp}\left(\underset{i=1}{\overset{N}{}}V(x_i)\right),$$
(6)
where $`\mathrm{\Delta }_N(x)=_{i<j}(x_jx_i)`$ and $`\mathrm{\Omega }(N)=(2\pi )^{N(N1)/2}/_{k=1}^Nk!`$ is the volume of the “angular” (non-diagonal) entries of $`\mathrm{\Phi }`$ . A general strategy to compute $`Z(N)`$ is to factorize the $`N`$ integrals on the r.h.s. of (6), which can be done provided there exists a set of orthogonal polynomials $`𝒫_k(x)`$ satisfying:

$$_{\mathrm{}}^{\mathrm{}}𝑑x\mathrm{e}^{V(x)}𝒫_k(x)𝒫_l(x)=h_k\delta _{k,l}.$$
(7)
In turn, this condition can be met provided the polynomials $`𝒫_k(x)`$ satisfy a three-term recurrence relation given by:

$$x𝒫_k(x)=𝒫_{k+1}(x)+S_k𝒫_k(x)+R_k𝒫_{k1}(x),$$
(8)
where $`S_k`$ and $`R_k`$ are specific to the matrix model . A crucial equality that relates the potential with the recursion relation is $`R_k=h_k/h_{k1}`$, and the partition function is given by:

$$Z(N)=\mathrm{\Omega }(N)N!\underset{k=0}{\overset{N1}{}}h_k=\mathrm{\Omega }(N)N!h_0^N\underset{k=0}{\overset{N1}{}}R_k^{Nk}.$$
(9)
In the following, we shall focus in the set of polynomials (3), and show that they are orthogonal with respect to an specific measure which is then used to define the associated matrix model. One anticipates, therefore, an infinite family of models parametrized by the momentum $`n`$, and defined by the family of potentials $`V_n(x)`$. Given the location of the zeroes of the polynomials (3), it is natural to perform a ”Wick rotation” in $`x`$ and define $`𝒫_n^k(x)=i^kP_n^k\left(ix+(n+1)/2\right)`$ as the set of $`W_{1+\mathrm{}}`$ orthogonal polynomials. They posses a definite parity under reflection in $`x`$; the first few are: $`𝒫_n^0(x)=1`$, $`𝒫_n^1(x)=x`$, $`𝒫_n^2(x)=x^2+(n^21)/12`$,$`𝒫_n^3(x)=x^3+x(3n^27)/20`$. The recursion relation (4) now takes the form:

$$x𝒫_n^k(x)=𝒫_n^{k+1}(x)+\frac{k^2(k^2n^2)}{4(4k^21)}𝒫_n^{k1}(x),$$
(10)
which implies the existence of the set of potentials $`V_n(x)`$ making them orthogonal. These functions are readily found once one recognizes that the polynomials $`𝒫_n^k(x)`$ are Hahn polynomials, for which the orthogonality relations are known. There are, in fact, two classes of Hahn polynomials: discrete and continuous. The first case arises for $`n`$ integer, such that there is a finite set of polynomials, since the recursion relation (10) stops at $`k=n`$. However, given (1), we are interested in the case for which the family is infinite: this is the continuous case, that arises for $`n`$ not an integer. For our purposes, it is convenient to consider the values $`n=ip`$, $`p=0,1,2,\mathrm{}`$ (i.e., imaginary momentum) since this yields a positive and non-degenerate inner product. This condition can be viewed as the Wick rotation in momentum space that corresponds to the one previously performed in $`x`$. Under these conditions, the orthogonality relations are given by :

$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x`$ $`g_p(x)`$ $`𝒫_{ip}^k(x)𝒫_{ip}^l(x)={\displaystyle \frac{|\mathrm{\Gamma }(k+1+ip)|^2}{(2k+1)\left(\begin{array}{c}2k\\ k\end{array}\right)^2}}\delta _{k,l},`$ (13)
$`g_p(x)`$ $`={\displaystyle \frac{\pi }{\mathrm{cosh}(2\pi x)+\mathrm{cosh}(p\pi )}}.`$ (14)
From this expression, one extracts the family of potentials $`V_p(x)`$ defining the matrix models as in (5):

$$V_p(x)=\mathrm{ln}g_p(x)=\mathrm{ln}\left[\frac{\mathrm{cosh}(2\pi x)+\mathrm{cosh}(p\pi )}{\pi }\right],p=0,1,2,\mathrm{},$$
(15)
which asymptotically behave as $`V_p(x)2\pi |x|2\mathrm{ln}(2\pi )`$ for $`|x|>\gamma _p/\pi `$, $`\gamma _p=\mathrm{ln}[2\pi \mathrm{cosh}(p\pi )]`$ (see figure 1).

From (14), one also knows that the norm of the $`kth`$ polynomial in the $`pth`$ matrix model is:

$$h_p^k=\frac{|\mathrm{\Gamma }(k+1+ip)|^2}{(2k+1)\left(\begin{array}{c}2k\\ k\end{array}\right)^2},$$
(16)
which is in agreement with the recursion relation (10) through the equality $`R_p^k=h_p^k/h_p^{k1}=k^2(k^2+p^2)/4(4k^21)`$. From (16), one can therefore define the partition function $`Z_p(N)`$ and free energy density $`f_p(N)`$ for the $`pth`$ model as:

$$Z_p(N)=\mathrm{\Omega }(N)N!\underset{k=0}{\overset{N1}{}}h_p^k=\mathrm{exp}(N^2f_p(N)).$$
(17)
The asymptotic form of $`f_p(N)`$ in the thermodynamic limit $`N\mathrm{}`$ is:

$$f_p(N)=\frac{1}{2}\mathrm{ln}N\left(\frac{3}{4}+\frac{1}{2}\mathrm{ln}\frac{8}{\pi }\right)+O\left(\mathrm{ln}N/N\right).,$$
(18)
independent of $`p`$ to the given order.

The thermodynamic limit leads also to a simple physical picture in terms of a fictitious one-dimensional plasma: consider the partition function (5) for the potential (15) (we consider for simplicity the case $`p=0`$). The asymptotic form of $`V_0(x)`$ is linear, and suggests the definition of a new scaling variable $`y=x/N`$, appropriate in the limit $`N\mathrm{}`$. This implies that $`V_0(x)2\pi N|y|2\mathrm{ln}(2\pi )`$ for $`|y|>1/\mathrm{\Lambda }`$, with $`\mathrm{\Lambda }=\pi N/\gamma _0`$. In the scaling limit, $`y`$ is defined in the entire real axis, except for a small region of the size of the cutoff $`1/\mathrm{\Lambda }`$ around the origin. We can write the partition function approximately as:

$$Z_0(N)\frac{\mathrm{\Omega }(N)}{N^N}\underset{i=1}{\overset{N}{}}dy_i\mathrm{exp}[N(2\pi \underset{i=1}{\overset{N}{}}|y_i|\frac{1}{N}\underset{i<j}{}\mathrm{ln}(y_iy_j)^2+f(N))],$$
(19)
where $`f(N)=(N1)\mathrm{ln}N+2\mathrm{ln}(2\pi )`$. In the large $`N`$ limit, $`Z_0(N)`$ is dominated by the saddle-point configurations $`\{y_1,\mathrm{},y_N\}`$ satisfying the equations:

$$2\pi \mathrm{sg}(y_i)=\frac{2}{N}\underset{j,ji}{}\frac{1}{y_iy_j},i=1,\mathrm{},N,$$
(20)
where $`\mathrm{sg}(y)`$ is the sign function. The scaled variable $`y`$ is appropriate for taking the continuum limit: $`y_iy`$, $`_i(\mathrm{})/N𝑑y\rho (y)(\mathrm{})`$, where $`\rho (y)`$ is the density function. The saddle-point equations are:

$`\pi \mathrm{sg}(y)`$ $`=`$ $`𝒫{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{\rho (y^{})}{yy^{}}}𝑑y^{},`$ (21)
$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\rho (y^{})𝑑y^{}`$ $`=`$ $`1,`$ (22)
where $`𝒫`$ denotes the Cauchy principal part of the integral, and the second equality is the normalization condition. These equations give a semiclassical physical interpretation to our system: they describe a system of fictitious equally charged “particles”, located on a one-dimensional rod, and subject to a neutralizing background. The particles repel each other with a two-dimensional Coulomb force, while they are attracted with constant force to an infinite two-dimensional uniformly charged plane, that is perpendicular to the rod. One can find a solution to equations (22) by assuming that the system is symmetric under $`yy`$, and therefore the total number of “particles” is even, with $`N/2`$ of them on each side of the infinite charged plane. Therefore, one seeks for a solution of the form $`\rho (y)=\rho _+(y)E(y)+\rho _{}(y)E(y)`$, with $`E(y)`$ the Heaviside function. In the following we shall focus on one side only, e.g. $`(+)`$. As it is customary, one introduces the function

$$F_+(y)=_{\mathrm{}}^{\mathrm{}}\frac{\rho (y^{})}{yy^{}}𝑑y^{},$$
(23)
defined on the complex half-plane $`\mathrm{Re}(y)>0`$. Assuming that this function has one cut running along the real axis from $`0`$ to $`\alpha `$, one can solve (22) by a standard procedure , yielding:

$$F_+(y)=\pi \pi \sqrt{\frac{y\alpha }{y}},$$
(24)
From (24) one determines the density:

$$\rho _+(y)=\sqrt{\frac{\alpha y}{y}},$$
(25)
where $`\alpha =1/\pi `$ is determined from the normalization condition (22). One can verify that along the cut and close to its extreme $`y\alpha `$, the behavior of $`\rho _+(y)`$ is in agreement with the expected one for a linear background potential .

Let us now focus on the application of the previous results to one-dimensional condensed matter models. Consider a free non-relativistic fermionic field on a one-dimensional lattice of unit spacing. We shall outline the effective field theory formulation of this problem in terms of the $`W_{1+\mathrm{}}`$ symmetry generators . While this description does not appear as necessary in the free theory, it becomes rather advantageous when including the effects of the interactions, as, for example, in the $`XXZ`$ Heisenberg model. We shall postpone the full discussion of the effects of the interactions for a forthcoming publication, but remark here that the scope of our approach is to develop a bosonization scheme valid to all orders in the expansion parameter, which is the inverse size of the system, given by $`N`$. In the thermodynamic limit, the non-relativistic field can be expanded in the basis of the relativistic Weyl field, and the Fermi sea interpreted as a Dirac sea . Moreover, the free hamiltonian $``$ can be written as the sum of two independent terms, $`_+`$ and $`_{}`$, describing the fluctuations around the right $`(+)`$ and left $`()`$ Fermi points at momenta $`\pm n_F`$, respectively, with $`n_F=(N2)/4`$; i.e., we consider the ground state of the system to be at half-filling, for simplicity. Both terms are decoupled, up to global (zero-mode) constraints. Therefore, we shall consider $`_+`$ only in the following:

$$_+=\underset{r=\mathrm{}}{\overset{+\mathrm{}}{}}f(n_F+r):a_r^{}a_r:,f(n)=\mathrm{cos}\left(\frac{2\pi }{N}n\right),$$
(26)
where $`a_r`$ are the fermionic Fock operators of the Weyl field, obtained from the original ones by shifting the momentum index, and $`r`$ is integer. The dispersion relation around the Fermi point given by $`f(n_F+r)`$ is non-linear, and can be expanded in terms of the zero-momentum $`W_{1+\mathrm{}}`$ polynomials (3) as :

$$_+=\underset{k=0}{\overset{\mathrm{}}{}}c_kV_0^k,$$
(27)
where the coefficients $`c_k`$ are such that:

$$\mathrm{sin}\left[\frac{2\pi }{N}\left(r\frac{1}{2}\right)\right]=\underset{k=0}{\overset{\mathrm{}}{}}c_kP_0^k(r).$$
(28)
Making the change of variables $`r(ir+1/2)`$ and passing to the basis $`𝒫_0^k(r)`$ allows us to determine $`c_k`$ using the orthogonality relation (7):

$$c_k=\frac{i^k}{h_0^k}_{\mathrm{}}^+\mathrm{}𝑑xg_0(x)\mathrm{sin}\left(\frac{2\pi }{N}r\right)𝒫_0^k(x)k=0,1,2,\mathrm{},$$
(29)
where for convenience we have made the redefinition $`2\pi /Ni2\pi /N`$, so as to deal with regular expressions; this analytic continuation should be undoed at the end of the computations. From (29) one immediately concludes that $`c_{2n}=0`$ for $`n`$ integer, and defines $`C_k=c_{2k+1}`$ ($`k=0,1,2,\mathrm{}`$) for the other cases, which can be worked out by considering the set of integrals:

$$I_n(y)=_{\mathrm{}}^+\mathrm{}𝑑xg_0(x)\mathrm{e}^{ixy}𝒫_0^n(x),n=0,1,2,\mathrm{}.$$
(30)
As a consequence of (10), they satisfy the recursion relations:

$$I_{n+1}(y)=i\frac{dI_n(y)}{dy}R_0^nI_{n1}(y),I_0(y)=\frac{y}{2\mathrm{sinh}(y/2)},$$
(31)
which are formally solved by $`I_n(y)=𝒫_0^n\left(i\frac{d}{dy}\right)I_0(y)`$. The coefficients (29) are then readily determined yielding:

$$C_k=\frac{(1)^k}{h_0^{2k+1}}I_{2k+1}\left(\frac{2\pi }{N}\right),k=0,1,2,\mathrm{}.$$
(32)
As and example, after analytically continuing back $`\pi /Ni\pi /N`$ as explained above,

$$C_1=\frac{6}{\mathrm{sin}(\pi /N)}\left[1\frac{\pi }{N\mathrm{tan}(\pi /N)}\right]=\frac{2\pi }{N}+\frac{7}{120}\left(\frac{2\pi }{N}\right)^3+O\left\{\left(\frac{2\pi }{N}\right)^5\right\},$$
(33)
which reproduces the results of . We have also verified the leading term in $`C_2=(2\pi /N)^3/6+\mathrm{}`$ yields the correct result.

Finally, we remark that the matrix models considered in this letter belong to the class of models with non-polynomial potentials, which include also the Penner  and Frustrated Spherical  models. According to the ‘t Hooft correspondence, one can associate to each matrix model a tiling of a two-dimensional random surface, such that the the number of sides and the ways of gluing the elementary poligons are determined by the algebraic powers of the monomials that define the potential . The class of matrix models with non-polynomial potentials possess the characteristic property of having an unbounded set of polygons and gluing vertices among them. We conclude by remarking that it would be interesting to explore whether other infinite-dimensional Lie algebras ecould give rise to further new classes of random matrix models.

GZ acknowledges the hospitality of the Abdus Salam Center for Theoretical Physics (Italy) and the Physics Department and I.N.F.N. of Florence (Italy).


text:
# On the existence of a bound tetraneutron

(Received: )

## Abstract

Following recent work in which events which may correspond to a bound tetraneutron (4n) were observed, it is pointed out that from the theoretical perspective the two-body nucleon-nucleon force cannot by itself bind four neutrons, even if it can bind a dineutron. A very strong phenomenological four-nucleon (4N) force is needed in order to bind the tetraneutron. Such a 4N force, if it existed, would bind 4He by about 100 MeV. Alternative experiments such as (8He,4n) are proposed to search for the tetraneutron.

In a recently reported experiment marques ; Orr  events were observed that exhibit the characteristics of a multineutron cluster liberated in the breakup of 14Be, most probably in the channel 10Be+4n. The lifetime of order 100 ns or longer suggested by this measurement, would indicate that the tetraneutron is particle stable. If confirmed, the existence of a bound tetraneutron is of great importance as it would challenge our understanding of nuclear few-body systems and nucleon-nucleon (NN) interactions. Here I would like to make several observations in this context.

I

As discussed in Refs marques ; Orr  the breakup 14Be $``$ 10Be+4n represents an attractive channel to search for a tetraneutron. The nucleus 10Be is strongly bound and the four neutron separation energy for 14Be is only about 5 MeV with respect to the 14Be $`^{10}`$Be + 4n breakup. It is known that in an external potential well neutrons can form drops ndrops , ndrops1 . In the attractive field of 10Be the last four neutrons could form a tetraneutron which might be liberated in the breakup.

In earlier experiments, the tetraneutron was searched using heavy-ion transfer reactions such as 7Li($`{}_{}{}^{11}B,^{14}O`$)4n Bel , 7Li($`{}_{}{}^{7}Li,^{10}`$C)4n Al  and double exchange reaction 4He($`\pi ^{},\pi ^+`$)4n U , G . These reactions should be strongly suppressed by spectroscopic reasons and a negative outcome from them could be easily anticipated.

II

No proper ab-initio four-body calculations of the tetraneutron with realistic two-body and three-body NN forces are known to the author. However, several theoretical calculations of the tetraneutron are available.
$`(i)`$ In Ref. bev  the tetraneutron, studied in the translationally invariant, symmetrized oscillator basis using the $`4\mathrm{}\omega `$ model space, was found to be unbound by 18.5 MeV. However, the oscillator basis is not appropriate for the description of unbound systems.
$`(ii)`$ No bound tetraneutron was found in Ref. gorb  within the angular potential functions method with semirealistic NN interactions. No search for the four-body resonance state was carried out there.
$`(iii)`$ No bound tetraneutron was found in Ref. varga  within the stochastic variational method on a correlated Gaussian basis for a range of simple effective NN potentials. No search for the four-body resonance state was carried out there.
$`(iv)`$ A search of the four-body resonance in the lowest order of the hyperspherical functions method (HSFM) gave a null result BBKE , SRV .
$`(v)`$ The energy behaviour of the eigen phases, studied in Ref. GNO  within the HSFM using the $`K_{max}`$ = 6 model space, led the authors to the conclusion that the tetraneutron may exist as a resonance in the four-body continuum at an energy of about $`13`$ MeV. However, such a conclusion is not convincing because no convergency of the eigen phases with an increase of the hyperangular momentum has been achieved in these calculations. Besides, a clear indication of the resonance has been seen with only one of the NN potentials used in the calculations, namely with the Volkov effective NN force V1 volkov . Volkov effective NN interactions reproduce the experimental binding energy of another four-body system, 4He. However, their singlet even and triplet even components are equal, and therefore with these potentials a singlet dineutron has exactly the same binding energy as the deuteron. In the particular case of V1, the dineutron is bound by 0.547 MeV. Therefore, V1 cannot be used in calculations of the multineutron systems. Another NN potential, used in Ref. GNO , namely that of Reichstein and Tang (RT) RT , reproduces the n-p triplet and p-p singlet scattering lengths and does not bind a dineutron. With this potential the energy derivatives of the eigen phases monotonically decrease with increasing energy. This means that resonances in the four neutron system are absent, at least within the model space considered.

To understand whether the RT potential can produce any resonance or bound state if the model space of the HSFM is increased, I have calculated the hyperradial potentials for the tetraneutron up to $`K_{max}`$ = 16 using the technique developed in Ref. T . The RT potential has only central even components $`V_{10}`$ and $`V_{01}`$ and the odd potentials are usually obtained from them using an arbitrary parameter $`u`$ so that $`V_{00}=(u1)V_{01}`$ and $`V_{11}=(u1)V_{10}`$. The present calculations have been performed with $`u=1`$ which corresponds to no-interaction in the odd partial waves. The diagonalised hyperradial potentials $`V_{diag}(\rho )`$, calculated for different model spaces $`K_{max}`$, are plotted in Fig.1 as a function of hyperradius $`\rho `$. Such potentials can be used to obtain a quick approximate solution for the binding energy in the extreme adiabatic approximation of the HSFM fabre . One can see that the $`V_{diag}(\rho )`$ have almost converged. They are purely repulsive, monotonically decreasing with the hyperradius, and do not show any sign of local attractive pockets. Therefore, the RT potential can neither bind the tetraneutron nor produce any resonances. To get a bound tetraneutron, the value of $`u`$ should be increased up to $`u=2.3`$. However, with this value, the isotopes 4H and 5He must become bound with respect to the t + n and 4He + n thresholds by about 6.5 MeV and 32 MeV respectively. In reality, these nuclei are unbound by 3 and 0.9 MeV.

For comparison, the calculations of the hyperradial potentials have been performed with Volkov NN potential V1 as well. The diagonalised hyperradial potentials $`V_{diag}(\rho )`$, shown in Fig. 2, reveal local attractive pockets for $`K_{max}>6`$. These pockets become negative for $`K_{max}>12`$, but they are too shallow to form a bound state with respect to the four-body decay. Although with $`K_{max}`$ = 16 convergence has not yet been reached, the general trend seen in Fig. 2, suggests that it is unlikely, that a tetraneutron, bound with respect to the 2n + 2n decay, may exist. To get a bound tetraneutron, the Majorana parameter $`m`$ should be changed from its standard value of 0.6 to $`m=0.2`$. Such a change provides E(4n) = $``$1.2 MeV and does not influence the binding energy of 4He, however, it binds 4H and 5He with respect to the t + n and 4He + n thresholds by 7 MeV 25 MeV respectively.

And finally, calculations with the realistic two-body isospin-conserving central NN interaction Argonne $`v_4`$ wiringa  have been performed as well and the corresponding hyperradial potentials $`V_{diag}(\rho )`$ are shown in Fig.1. These purely repulsive potentials have almost converged and their behaviour excludes any possibility of either a bound or a resonance state.

III

The results of the theoretical calculations suggest that the tetraneutron must be unbound and most likely should not exist as a resonance, at least when only two-body central forces are considered. It is remarkable that even when an effective NN potential binds the dineutron (the case of V1), it still cannot bind two dineutrons, although a resonance in such a system should exist. In order for the two-body central force to be able to bind the tetraneutron, a huge unphysical attraction should be introduced in the triplet odd potential. Such an attraction would strongly overbind 4H, 5He and other known $`A>4`$ nuclei. The reason for necessity of huge attraction lies in the relative number of nucleon pairs in even and odd states allowed by the Pauli principle. For the tetraneutron with $`L=0`$ and $`S=0`$ the probabilities $`P`$ to find a pair of nucleon in the singlet even (01) and triplet odd (11) states are equal, $`P_{01}=P_{11}`$ = $`\frac{1}{2}`$. For the tetraneutron with $`L=1`$ and $`S=1`$, these probabilities are $`P_{01}=\frac{1}{3}`$ and $`P_{11}=\frac{2}{3}`$. The NN force is not sufficiently attractive in the singlet even state, and is repulsive in the triplet odd state, which leads to the absence of the tetraneutron. For comparison, for the $`L=0`$, $`S=0`$ state of the closed shell nucleus 4He $`P_{10}=P_{01}=\frac{1}{2}`$ and for the $`L=1`$, $`S=1`$ state of 4H $`P_{10}=P_{01}=P_{11}=\frac{1}{3}`$. The binding energies of 4He and 4H are $``$28.3 MeV and $``$5.2 MeV respectively. It is clear that the binding of a nucleus is strongly correlated with the probabilities to find a pair of nucleons in even states, especially in the triplet even state (10) where the $`np`$ bound state exists.

IV

The fact that, despite the theoretical predictions made with central two-body forces, six events 14Be $`^{10}`$Be + 4n have possibly been observed, means that either non-central interactions and/or three-nucleon (3N) forces, and/or a strong four-nucleon (4N) force may bind the tetraneutron. If not, then a different interpretation of the six events from Refs marques ; Orr  should be looked for.

The first possibility should be rejected, because the calculations of Ref. gorb  indeed included the spin-orbit and tensor forces. The second possibility does not look convincing either because the contribution of the 3N potential to the binding energies of the neutron drops 7n and 8n, calculated with Green’s function Monte Carlo Methods in Ref. pieper , is about 1 to 5 MeV which is small compared to 4He and to other $`A=7`$ and $`A=8`$ nuclei. A similar contribution from the 3N force would be too small to bind the tetraneutron. Therefore, if what has been seen in the 14Be experiment was indeed the tetraneutron, it could indicate at presence of a 4N force.

At present, the results of the ab-initio calculations of the $`3A8`$ nuclei suggest that a 4N force is not needed to fit the observed binding energies of light nuclei at the 1$`\%`$ level pieper . Therefore, either the 4N contributions are smaller than 1$`\%`$ for these nuclei or parts of the 3N forces are mocking up their effects pieper . On the other hand, the latest models of the 3N force used to calculated the binding energies of these nuclei have not yet been tested in the description of polarization observables of the low-energy $`Nd`$ scattering where the contribution of the 3N force is very important kiev , cadman . Besides, the latest numerically accurate calculations employing a 3N force do not reproduce the proton analyzing power for the p-3He scattering viviani . The simultaneous fit of the $`Nd`$ and $`p^3He`$ polarization observables and of binding energies of the lightest nuclei may leave room for a 4N force.

In order to get an idea of what strength is needed for the 4N force, the HSFM calculations of this paper have been repeated with the RT potential using a 4N force simulated by the potential $`V_{4N}(\rho )=W_0e^{\alpha \rho }`$. The values $`\alpha =0.7,1.2`$ and 1.5 $`fm^1`$, used in calculations, were the same as the range of the phenomenological term of the 3N spin-orbit force introduced in Ref. kiev . The tetraneutron becomes bound if the corresponding values $`W_0`$ are equal to $`410`$, $`1460`$ and $`2530`$ MeV respectively. These values are two orders of magnitude larger as compared to the values $`1`$, $`10`$, and $`20`$ MeV obtained for the 3N spin-orbit force kiev . For the Volkov potential V1, these strengths are similar: $`W_0`$ = $`320`$, $`1565`$ and $`2900`$ MeV. If the same 4N force existed in 4He, the binding energy of 4He, calculated with V1, would be $`88`$, $`82`$ and $`134`$ MeV respectively. Therefore, the 4N force, if it exists, should be strongly $`T`$-dependent. In principle, such a strong $`T=2`$ 4N force could be noticable in $`AZ4`$ nuclei, especially near the neutron edge of stability, where the number of the $`T=2`$ four-nucleon states increases. However, HSFM calculations of the neutron-rich helium isotopes 6,8He, performed with V1 T , do not indicate a need for a strong 4N force.

To get a real understanding of the role of the 4N force, ab-initio calculations of the tetraneutron with realistic 2N and 3N forces must be performed. The HSFM is the best method for such calculations because it is well suited to systems without bound subsystems, which is the case for the tetraneutron, and because it can also treat a four-body continuum.

V

Since the discovery of the tetraneutron would most certainly require a revision of modern realistic models of the three-nucleon interaction, it is extremely important to confirm or reject the results reported in Ref. marques . Therefore, independent searches for the tetraneutron in alternative reactions should continue. Among such reactions, $`\alpha `$ transfer reactions $`A(^8He,^4n)A+\alpha `$ deserve special attention. Using a light target and detecting a product nucleus $`A+\alpha `$, one can study not only a hypothetical tetraneutron bound state but also the four-neutron continuum as well.

The spectroscopic factor $`S`$ for the $`^8He|^4He^4n`$ overlap can be estimated using the translation invariant shell model TS . If the probability amplitudes of the $`|L=0,S=0`$ and $`|L=1,S=1`$ configurations are $`\alpha _0`$ and $`\alpha _1`$ in 8He and $`\beta _0`$ and $`\beta _1`$ in 4n respectively, then $`S=\left(\sqrt{\frac{2}{3}}\alpha _0\beta _0+\sqrt{\frac{5}{6}}\alpha _1\beta _1\right)^2`$.

The choice of the target $`A`$ must be determined by kinematic and spectroscopic considerations. It should be light enough for residual nucleus $`A+\alpha `$ to come out of the target and its wave function should overlap strongly with the $`A+\alpha `$ wave function. The obvious targets of this type are the deuteron, triton, 3He and 4He. The list of other light targets with large spectroscopic factors $`S_\alpha `$ is given in Ref. kurath .

As an example, in Figs. 3 and 4 the DWBA $`\alpha `$-transfer cross sections of the d($`{}_{}{}^{8}He,^4n`$)6Li and 12C($`{}_{}{}^{8}He,^4n`$)16O reactions are shown as a function of the four-body binding energy of the hypothetical tetraneutron. The energies of these reactions have been chosen based on availability of optical potentials in the entrance and exit channels, which have been taken from perey  and bec . In addition, the energy choice for the 12C target was also motivated by experimental measurements of the similar reaction 12C($`{}_{}{}^{8}Li,\alpha `$)16N which revealed very large cross sections at small angles, about 25 mb/sr bec . The latest theoretical value $`S_\alpha `$ = 0.87 NWS  has been used for 6Li and $`S_\alpha `$ for 8He was taken to be $`\frac{2}{3}`$. According to Figs. 3 and 4, the d($`{}_{}{}^{8}He,^4n`$)6Li reaction is better suited for the tetraneutron search because its cross sections are not strongly influenced by unknown tetraneutron binding energy. The cross sections and 8He beam intensities make this experiment feasible in the very near future. In contrast, the cross sections of the 12C($`{}_{}{}^{8}He,^4n`$)16O reaction decrease dramatically with increasing the binding of the tetraneutron and become small if the tetraneutron is bound.

Summarizing, due to the large probability for a pair of neutrons to be in the triplet odd state, the two-body nucleon-nucleon force cannot by itself bind four neutrons, even if it can bind a dineutron. A very strong phenomenological $`T`$-dependent four-nucleon force is needed in order to bind the tetraneutron. Alternative experiments such as (8He,4n) are proposed to search for the tetraneutron.

I am grateful to W. Catford for providing me with the text of Ref. marques  prior to the publication. I am also grateful to N.Orr and W. Catford for useful comments concerning my paper. I also thank K. Varga for providing me with the gaussian expansion of the Argonne NN potential.