text:
# Faint Radio Sources and Star Formation History

## 1 Introduction

In the last few years, a variety of observational methods have been used to study the global star formation history of the universe at a range of redshifts. Figure 1 compiles the results from several of these studies, scaling all to the same cosmology and IMF. In plotting the points, we used the corrections for extinction by dust calculated by Steidel et al. (1999, their figure 9). The diagram shows significant scatter in the star formation density at each redshift. Most studies agree, however, that the star formation density rises rapidly from $`z=0`$ to $`z=1`$. Beyond a redshift of 1 it is unclear whether the star formation density decreases significantly (as suggested at first by Madau et al., 1996) or stays roughly constant (for example, Steidel et al., 1999).

Radio observations have important advantages in determining the global star formation history, and are a useful complement to studies at other wavelengths. Unlike calculations based on ultraviolet and optical observations, there is no need to make uncertain corrections for dust extinction since the radio emission at $`\nu 1`$ GHz passes freely through dust. Compared to far-infrared and submillimeter studies, interferometric radio observations typically have better positional accuracy, allowing for more reliable identifications with objects detected at other wavelengths (Richards, 1999; Downes et al., 1999). Finally, since far-infrared emission is due to re-heated dust, the original source of energy (whether star formation or AGN) can be unclear; radio properties, such as spectral index and morphology, can help distinguish between these. We do note, however, that contamination of the radio flux by emission from AGN is a problem (addressed in §2.5). In addition, relatively few high redshift star-forming radio sources have as yet been detected (though current and planned deep radio surveys will rapidly change that). As a consequence, the statistical sample used in this work is small.

Our strategy in this paper is as follows. We use very sensitive radio surveys to detect star forming galaxies at high redshift. Not all of the sources in these surveys are star-forming (some are probably AGN), but we deal with this problem by defining data samples that give lower and upper limits to the star formation history, as described in §2. Next, in §3, we use these data to determine the evolving luminosity function for star-forming radio sources. The redshift and flux distribution of the sources, as well as the total extragalactic radio background, are used to constrain the evolution. In §4, we use the well-known relationship between radio luminosity and star formation rate (and discuss the assumption that this relationship holds for all redshifts) to find the star formation history directly from the observed radio sources. The evolving luminosity function is used to correct for faint sources below the observational detection limits. We discuss our conclusions in §5.

Unless stated otherwise, we assume a cosmology of $`H_0=50`$, $`\mathrm{\Omega }_m=1`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$, and a non-evolving Salpeter Initial Mass Function with a stellar mass range of $`0.1100M_{}`$. We use a radio spectral index of $`\alpha =0.4`$ (where $`S\nu ^\alpha `$), which is appropriate for faint sources selected at 5 or 8 GHz (Windhorst et al., 1993; Richards, 2000).

## 2 The Data Sets

### 2.1 Surveys

Three fields have been observed to micro-Jansky sensitivities at centimeter wavelengths and also have extensive photometric and spectroscopic data: the Hubble Deep Field (HDF) at 8 GHz (Table 2), the SSA13 field at 8 GHz (Table 3), and the V15 field at 5 GHz (Table 2). Table 1 summarizes information on the three fields, all of which were observed at the Very Large Array (VLA111The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.). Tables 2, 3, and 4 list the individual sources in each field; the columns are as follows:
Column 1: The source name.
Columns 2,3: $`S_8`$, $`S_5`$, and $`S_{1.4}`$ are the radio flux densities at 8 GHz, 5 GHz, and 1.4 GHz respectively. If a measured value of $`S_{1.4}`$ is not available, we use the spectral index shown in Column 4 to calculate $`S_{1.4}`$ and list the value in parentheses in Column 3.
Column 4: The radio spectral index, $`\alpha `$, defined as $`S\nu ^\alpha `$. If 1.4 GHz measurements are not available, we use assume a spectral index of 0.4 unless this violates the survey detection threshold at 1.4 GHz.
Column 5: The primary beam correction factor, $`B_i`$ (see §2.2, eq. 1).
Column 6: The galaxy type (see §2.4): sim refers to spiral, irregular, or merger; el refers to elliptical; fr refers to faint ($`I>25`$) or red ($`IK>4`$); and un refers to unknown type or undetected. The galaxy type in parentheses is the assumed galaxy type in the case of un.
Column 7: The redshift type (see §2.3): sp refers to spectroscopic, ph refers to rough photometric (based on I or HK magnitudes), a refers to random assignment, and afr refers to random assignments for galaxies of type faint/red.
Column 8: The redshift used in calculations.
Columns 9,10: The I and HK magnitudes. K magnitudes are converted to HK magnitudes by $`K=HK^{}0.3`$.
Column 11: The maximum redshift, $`z_{max}`$, at which the galaxy would have been detected, based on its emitted luminosity $`L_{e,1.4}`$.
Column 12: The log of the luminosity emitted at 1.4 GHz.
Column 13: The star formation rate for each individual galaxy, derived from its radio luminosity ($`L_{e,1.4}`$) and equation 15. If the source is elliptical (or assumed elliptical), the emission is probably contaminated by AGN and the calculated star formation rate is only an upper limit. For some of these sources, the AGN contamination causes the calculated star formation rate to be unphysically large (greater than 5000$`M_{}/yr`$), so it does not provide an interesting upper limit and we do not list it; we do, however, include these sources in our “upper” sample to give conservative upper limits on our results.
Column 14: The samples to which the source was assigned (see §2.5): U refers to the upper sample, M to the middle sample, and L to the lower sample.

For the first time we have a sample of microJy radio sources with nearly complete optical identifications and nearly 60% complete spectroscopic redshift measurements. Others are doing similar work on faint star forming radio galaxies (Hopkins et al., 1999; Mobasher et al., 1999; Benn et al., 1993; Gruppioni, Mignoli, & Zamorani, 1999) with larger catalogs of sources. While our survey samples have fewer sources, we have generally more sensitive radio flux limits and more complete optical followup. In principle, this allows us to probe higher redshifts, and to increase the fraction of sources identified with star-forming galaxies (see §2.4).

### 2.2 Primary beam corrections

In each of the three radio surveys, the flux threshold varies significantly across the field due to the shape of the primary beam response of the VLA antennas. The flux limit listed in Table 1 is for the center of the field; the limiting flux $`S_{lim}`$ increases to the edge of the field by about 1.5 for SSA13, and by about 10 for the larger HDF and V15 fields. To find the total surface density of sources $`n`$ (the number of sources per angular area on sky), we determine the contribution of each source by considering the portion of the field in which that source could have been detected. Therefore, each source $`i`$ contributes

$$B_i=\frac{1}{A_i(S_i,S_{lim})}$$
(1)
to the surface density of sources, where $`A_i`$ is the solid angle on the sky in which the flux $`S_i`$ of source $`i`$ would be greater than the sensitivity limit of the radio survey $`S_{lim}`$. A faint source which could only be detected at the center of the field (small $`A_i`$) contributes more to the average surface density $`n`$ than a strong source which could be detected over the entire primary beam area (large $`A_i`$). To determine $`A_i`$ for each source, we used the shape of the VLA primary beam.222The gain of the VLA primary beam is well-matched by $`g(r)=[cos(\frac{0.23226+74.567639r\nu }{57.2957795})]^6`$ where $`r`$ is the distance from the beam center in degrees and $`\nu `$ is the observing frequency in GHz (Oort & Windhorst, 1985). The total surface density of sources for a survey is then

$$n=\mathrm{\Sigma }B_i$$
(2)
Other instrumental effects that affect the point source sensitivity across the field, such as bandwidth smearing, time delay smearing, and geometrical smearing (Richards, 2000), are negligible in the three VLA surveys. Pascarelle et al. (1998) discuss the importance of surface brigtness corrections in determining star formation history. However, the relatively low resolution of these radio surveys (3′′ for the HDF and 3-10′′ for the SSA13 and V15 fields), combined with the resolution correction as a function of flux density (Windhorst et al., 1990, 1993), suggests that few of these sources are resolved and thus no correction for surface brightness biases have been made.

### 2.3 Redshifts

The redshifts for the sources in the sample are either spectroscopic measurements, estimates from I or HK-band magnitudes, or random assignments (see Tables 1 through 4 and Figure 3). About 58% of the sources have spectroscopic redshifts. The highest spectroscopic redshift in the sample is a source at $`z=4.42`$ in the HDF (Waddington et al., 1999); however, there is some evidence that this source may contain an AGN, so that its radio flux is not dominated by star formation (see §2.5 for how this source is treated in the calculations).

For 13% of the sources, approximate redshifts were found from I and HK band magnitudes. Windhorst et al. (1994b) used Bruzual-Charlot models (1993) to find the dependence of I and HK magnitude on redshift for milliJy radio sources; these models are plotted in Figure 2. When comparing these models with the 45 spectroscopic redshifts in our sample of fainter microJy radio sources, we found a significant dependence on radio flux density. We compensated by shifting the magnitude scale of these models to fit the I or HK values of our sources that do have spectroscopic redshifts, with different shifts for different radio flux ranges (as listed in the caption of Figure 2). We then used the revised curves to estimate redshifts for our sources that have I or HK band magnitudes but not spectroscopic redshifts. The resulting redshifts are crude, but are better than random assignments. We converted $`K`$ magnitudes to $`HK^{}`$ using $`K=HK^{}0.3`$ (Barger et al., 1999). The $`I(z)`$ model is double valued for $`z1`$, and the $`HK(z)`$ model increases sharply for $`z3`$, so when the I or HK magnitude of a source indicated a redshift above these limits, we randomly assigned a redshift instead (see below).

For the remaining 29% of the sources, neither spectroscopic redshifts nor rough photometric redshifts were available. Rather than removing these sources from the sample, we assigned redshifts in the following manner. We first separated the sources into two groups: those with very faint (I$`>25`$) or red (I-K$`>4`$) optical identifications, and those with brighter optical identifications. The group with brighter identifications (10 sources) were assigned redshifts randomly selected from the list of spectroscopic redshifts of star-forming galaxies in the sample. The group with very faint or red optical identifications (12 sources) are probably star-forming galaxies at redshifts greater than 1 (Richards et al., 1999; Barger, Cowie, & Richards, 2000). These sources were assigned redshifts randomly in the range $`z=13`$.

The redshift distributions and total source densities of the three surveys are strikingly different (see Figure 3). The HDF and SSA13 surveys were both performed at 8 GHz with similar flux limits, yet the average source density (including all redshifts) is 1.3 sources/arcmin2 in the HDF, but 2.7 sources/arcmin2 in the SSA13 field, over twice as great. The total source density of the V15 field at 5 GHz is 0.7 sources/arcmin2, which is nearly the same as the HDF field when the differences in flux limit and observing frequency are taken into account. The redshift distribution in the three fields peaks at somewhat different redshifts (see Figure 3), possibly due to galaxy superclustering or other high-redshift structure, although all three fields peak at $`z<1`$ and have a long tail which extends to $`z3`$. The differences in the redshift distribution of the fields are probably due to cosmic variance (note that each field is sampling only a small solid angle). Since these fields were generally chosen to be free of bright sources, the number counts may be too low in the HDF and V15 fields rather than too high in the SSA13 field (indeed, Richards, 2000 reports a deficit of radio sources detected at 1.4 GHz in the HDF). To deal with the differences between fields, we average the three data sets together in our calculations.

### 2.4 Optical Identifications

Figure 4 indicates the nature of the available optical identifications of the radio sources. Known QSOs were removed from the sample (two from the SSA13 field, one from V15), as was one star in the V15 field and are not shown in the figure. For the remaining sources, the three fields appear to be significantly different. In the HDF, about 50% of the radio sources have star-forming counterparts (spirals, mergers, and irregulars; Richards et al., 1998), another 30% are in the red/optically faint category discussed above (several of which are identified with bright sub-mm objects and may include star-forming galaxies; Barger et al., 2000), and only 20% are identified with elliptical galaxies which presumably are associated with low luminosity AGN. In the SSA13 field, 50% are identified with star-forming or red/optically faint galaxies, and 50% are of unknown type. In the V15 field, 15% of the sources have unknown galaxy types, and 40% of the sources are likely to be star-forming or in the red/optically faint population (Hammer et al., 1995). A further 35% are claimed to be elliptical/AGN counterparts based on deep I-band images from the Canada-France Hawaii Telescope, and another 15% are classified by Hammer et al. as AGN based on emission line studies. Thus, Hammer et al. (1995) report a larger fraction of low-luminosity AGN in the V15 field (50%) than observed in the HDF (20%). However, HST/WFPC2 images of these identifications from the Groth Strip survey (Hammer, 1996) show some ambiguity in the optical identifications, indicating a higher fraction of late type galaxies than originally reported by Hammer et al. (1995); thus some uncertainty remains. Also, the V15 field lies only 20 arcmin from the cluster associated with 3C295, and there is another supercluster or redshift structure at $`z=0.98`$ within the field (Le Fevre et al., 1994). These structures have probably caused some bias in the identifications, increasing the fraction of early type radio galaxies. Thus we adopt the statistics from the HDF and SSA13 surveys, which imply that about 50% of the sources have disk or late type galaxy counterparts, 30% have red/optically faint identifications, and 20% are associated with ellipticals and low luminosity FR I type AGN.

### 2.5 Strategy for dealing with incomplete identifications and redshifts

The goal of this work is to calculate the global star formation history based on star-forming radio sources. Since not all of the sources are star-forming, and only 58% have spectroscopic redshifts, we must define the target population carefully. To do this, we separate the data into subsets in order to calculate lower and upper limits on the luminosity function and hence the star formation history (see Figure 5):

* a “lower limit” sample (23 sources): only those sources which are both identified with spiral/irr/mergers and have spectroscopic redshifts. These are the sources which definitely belong in the population of interest.
* a “middle value” sample (37 sources): only those sources for which two criteria are met:

+ redshift is spectroscopic or based on I- or HK-band magnitude (no randomly assigned redshifts)
    + galaxy type is spiral/irr/merger, or faint/red. In addition, about 80% of the “unknown” identifications are assumed to be spiral/irr/merger or faint/red and are included here. The highest redshift source ($`z=4.42`$ in the HDF) is not included here because Waddington et al. (1999) argue that it contains an AGN component.

This sample is our best estimate of the true redshift distribution of star-forming radio galaxies.
* an “upper limit” sample (all 77 sources): all sources (including identifications with elliptical, emission line, and a few Seyfert galaxies, but not verified QSOs which were removed from sample). Redshifts were assigned for those sources without spectroscopic or I- or HK\- band estimates. This sample shows the maximum star formation rate that the data would allow, assuming all radio flux from all detected sources is due to star formation and that the redshift estimates are correct.

## 3 Evolution of the luminosity function

Next we determine the evolution of the luminosity function for this population of faint star-forming radio galaxies, using the data described in §2. In §4.3 we will use this evolving luminosity function to build a model of the star formation history. Since our data contains very few low-redshift objects, we can not fit for the shape of the local luminosity function. Instead, we use the local luminosity function found by Condon (1989), and fit for the evolution of that function in luminosity and number density. In §3.1 we calculate the luminosity function directly from the data, in §3.2 we describe the evolution model, and in §3.3 we describe the observational constraints on that model and the resulting best fit.

We convert all observed luminosities to a rest-frame frequency of 1.4 GHz, since most of the work on the local luminosity function has been done at this frequency. Our samples are defined at 5 and 8 GHz, but some sources have also been detected at 1.4 GHz. We use the observed 1.4 GHz flux densities when available, and for the remaining sources we assume a spectral index of $`\alpha =0.4`$ (see §1) to convert to 1.4 GHz, unless this violates an observed limit on the 1.4 GHz flux density. All source luminosities are then converted from the observed 1.4 GHz value to their rest-frame 1.4 GHz value. Thus, the observed luminosity of each galaxy $`L_{o,\nu }`$ at an observing frequency $`\nu `$ and redshift $`z`$, is converted to the emitted luminosity at 1.4 GHz rest-frame frequency using

$$L_{e,1.4}=L_{o,\nu }\left(\frac{\nu }{1.4\mathrm{GHz}}\right)^\alpha (1+z)^\alpha $$
(3)
We define the luminosity function $`\varphi (L_{e,1.4})`$ as the number per comoving Mpc3 per $`d\mathrm{log}_{10}L`$ of star-forming radio sources with emitted luminosity $`L_{e,1.4}`$(W/Hz)at 1.4 GHz.

### 3.1 Luminosity function estimated from the data

We can calculate the luminosity function directly from the detected sources for those luminosity and redshift ranges which are sampled by the data sets described in §2. For each bin in luminosity ($`L_{min}<L_{e,1.4}<L_{max}`$) and redshift ($`z_{min}<z<z_{max}`$), the luminosity function is

$$\varphi (L_{e,1.4},z)d\mathrm{log}_{10}L=\mathrm{\Sigma }_i\frac{B_i}{V_c[z_{min},z_{max}(L_i)]}$$
(4)
where $`B_i`$ is the surface density, corrected for the primary beam (eq. 1), and $`z_{max}(L_i)`$ is the largest $`z`$ in the bin for which the luminosity of the source $`L_i`$ was above detection limit of survey. $`V_c`$ is the comoving volume (in Mpc3) between $`z_{min}`$ and $`z_{max}`$ for solid angle $`\mathrm{\Delta }\mathrm{\Omega }`$,

$`V_c(z_{min},z_{max},\mathrm{\Delta }\mathrm{\Omega })`$ $`=`$ $`{\displaystyle 𝑑\mathrm{\Omega }r^2(z)𝑑r}`$ (5)
$`=`$ $`{\displaystyle \frac{\mathrm{\Delta }\mathrm{\Omega }}{\mathrm{ster}}}\left({\displaystyle \frac{\mathrm{ster}}{1.18\times 10^7\mathrm{arcmin}^2}}\right){\displaystyle \frac{[r^3(z_{max})r^3(z_{min})]}{3}},`$
where the comoving distance is

$$r(z)=\frac{2c}{H_0}\left(1\frac{1}{\sqrt{1+z}}\right)$$
(6)
for our assumed cosmology (see §1).

The binned luminosity function was calculated using eq. 4 for the lower, middle, and upper samples described in §2.5, using the average of the three surveys (bottom panel of Fig. 5). The result is shown in Figure 6, where the data points are from the “middle” sample, and are plotted at the average of the luminosities in each bin. Vertical error bars are either the lower and upper limits (from the samples described in §2.5), or the Poisson errors ($`1/\sqrt{N}`$ weights from the number of galaxies per bin), whichever is larger (generally the Poisson errors dominate for the low redshift data, and the lower/upper sample limits dominate for the high redshift data). Horizontal error bars are the range of source luminosities in each bin. Bins were chosen such that each contains 4 to 6 galaxies (except for the lowest redshift bin which has only 2 galaxies).

### 3.2 Description of evolving luminosity function model

We now build a model of the evolving luminosity function in order to compare it to several observables. In §3.3 we describe the observational constraints, the fitting process, and the resulting fit of this model to the observed data. Here we describe the model and its free parameters.

We use the local 1.4 GHz luminosity function for star-forming/spiral galaxies from Condon (1989, eq. 8 and discussion after eq. 7), adopting different notation,

$`\mathrm{log}_{10}[\varphi (L_{e,1.4})]d\mathrm{log}_{10}L=28.83+Y1.5\mathrm{log}_{10}L_{e,1.4}`$
$`\left[B^2+{\displaystyle \frac{1}{W^2}}(\mathrm{log}_{10}L_{e,1.4}X)^2\right]^{1/2}d\mathrm{log}_{10}L,`$ (7)
with the fitted parameters for star-forming galaxies of $`Y=2.88`$, $`X=22.40`$, $`W=2/3`$, and $`B=1.5`$. The factor of 28.83 includes unit conversions and the conversion from magnitudes ($`d\mathrm{log}_{2.5}L`$) to base 10 ($`d\mathrm{log}_{10}L`$).

To describe the evolution of the luminosity function, we use the functional form suggested by Condon (1984b, eq. 24), a power-law in $`(1+z)`$ with an exponential cut-off at high redshift. The luminosity evolves as

$$f(z)=(1+z)^Q\mathrm{exp}\left[\left(\frac{z}{z_q}\right)^q\right],$$
(8)
and the number density evolves as

$$g(z)=(1+z)^P\mathrm{exp}\left[\left(\frac{z}{z_p}\right)^p\right].$$
(9)
This gives six free parameters $`\{Q,q,z_q,P,p,z_p\}`$ to use in describing the evolution. Thus the general expression for the evolving luminosity function is (Condon, 1984a)

$$\varphi (L_{e,1.4},z)=g(z)\varphi (\frac{L_{e,1.4}}{f(z)},0).$$
(10)
Once we know the evolving luminosity function, it can be used to predict the observed redshift distribution, $`n(z)`$. The number of sources between $`z_{min}`$ and $`z_{max}`$ that could be detected in a survey of angular area $`\mathrm{\Delta }\mathrm{\Omega }`$ and flux limit $`S_{lim}`$ at frequency $`\nu `$ is

$$n(z)=V_c(z_{min},z_{max},\mathrm{\Delta }\mathrm{\Omega })_{L^{}(S_{lim},z)}^{inf}\varphi (L_{e,1.4},z)d\mathrm{log}_{10}L$$
(11)
where the lower limit of the integral is the luminosity corresponding to the flux limit $`S_{lim}`$ at the redshift $`z`$. The comoving volume $`V_c`$ is defined in eq. 5.

The evolving luminosity function can also be used to predict the extragalactic background due to this population. The background intensity at observing frequency $`\nu _0`$ is (Dwek et al., 1998)

$$I(\nu _0)=\frac{1}{4\pi }\rho (\nu ,z)\left|\frac{cdt}{dz}\right|𝑑z$$
(12)
where the luminosity density, $`\rho `$, emitted at redshift $`z`$ and frequency $`\nu `$ is found from the luminosity function,

$$\rho (\nu ,z)=L_{e,1.4}\varphi (L_{e,1.4},z)d\mathrm{log}_{10}L,$$
(13)
and

$$\left|\frac{dt}{dz}\right|=\frac{1}{H_0(1+z)^{5/2}}$$
(14)
for the assumed cosmology (§1).

### 3.3 Fitting the evolving luminosity function model to the data

We can now compare the evolution model to the observed data in order to fit for the evolution parameters. The model is constrained by three observables:

1. The redshift distribution, $`n(z)`$. We use the “middle” sample (defined in §2.5) for the average of three surveys, shown as the hashed area in the bottom panel of Figure 5.
2. The observed luminosity function, shown as the data points in Figure 6. We use the error bars shown in the figure (the larger of lower/upper limits and Poisson errors).
3. The extragalactic radio background, which is an important constraint on the integral of the luminosity function.

The first two constraints are not independent from each other, but both are needed. In order for the observed luminosity function to have 4-6 sources per bin, only coarse redshift resolution is possible; the $`n(z)`$ distribution allows for more detailed redshift information, but does not include the luminosity information.

The extragalactic radio background is about half due to star-formation activity and half due to AGN (see discussion in Haarsma & Partridge, 1998). In order to isolate the part of the radio background due to star formation, we use the far-infrared (FIR) background found by DIRBE of $`1.15\pm 0.20\times 10^{20}\mathrm{Wm}^2\mathrm{ster}^1\mathrm{Hz}^1`$ near 200$`\mu \mathrm{m}`$ (Hauser et al., 1998), assumed to be due primarily to star formation. The FIR-radio correlation (Helou, Soifer, & Rowan-Robinson, 1985; Condon, Anderson, & Helou, 1991) can then be used to predict the portion of the radio background due to star-formation, which at 1.4 GHz is $`\rho =3.2\pm 0.6\times 10^{23}\mathrm{Wm}^2\mathrm{ster}^1\mathrm{Hz}^1`$.

For each trial set of evolution parameters $`\{Q,q,z_q,P,p,z_p\}`$, we calculate the model prediction for $`n(z)`$, the evolving luminosity function, and the radio background due to star formation. The evolution parameters are adjusted to improve the model fit to the three data constraints, using a downhill simplex algorithm (Press et al., 1992) to find the global $`\chi ^2`$ minimum. Since the $`n(z)`$ and luminosity function constraints are not independent from each other, the reduced $`\chi ^2`$ can not be used to calculate the “goodness of fit” or to quantitatively compare the quality of different fits, but its minimum still indicates the parameters of the best available fit.

Our best fit is $`\{Q=3.97,q=1.02,z_q=1.39,P=0.0579,p=23.1,z_p=14.3\}`$. The resulting evolution factors $`f(z)`$ and $`g(z)`$ are plotted in Figure 7. The $`n(z)`$ distribution predicted by the model is shown in Figure 5 (bottom panel). The peak of the model $`n(z)`$ distribution falls at a lower redshift ($`z0.3`$) than the peak in the data ($`z0.5`$), but the tail of the distribution is reasonable and the total number density under the curve is similar (within 5%) for the data and the model (models with a peak at higher redshift tend to have a much shorter tail or a larger total number of sources $`n`$ and thus a larger discrepancy with the observed total). The luminosity function predicted by the model is shown as the curves in Figure 6, and is a good fit to the data points, except for the $`z=0.050.2`$ bin (which includes only two galaxies), and the $`z=0.60.9`$ bin (where one point is too high and the other is too low). The model predicted star-forming radio background is $`3.0\times 10^{23}\mathrm{Wm}^2\mathrm{ster}^1\mathrm{Hz}^1`$, which is a reasonable fit to the observed value. While the model does not perfectly match the three data constraints, it is the best compromise between them. Models that give a better fit to the observed $`n(z)`$ shape result in a poor fit to the other two data constraints. For instance, some evolution models can produce a longer tail on the $`n(z)`$ distribution, but that raises the total background and the total surface density $`n`$ significantly above the observed levels.

In the early stages of this work, it seemed that the full 6 parameters of our evolution model (eqs. 8, 9) were necessary to achieve a good fit. In the end, the best fit model shows virtually no number density evolution, and only a mild turn-over in luminosity evolution. Pure luminosity evolution {i.e. $`f(z)=(1+z)^Q`$ and $`g(z)=1`$} has often been suggested in the literature (Rowan-Robinson et al., 1993; Hopkins et al., 1998). It turns out that pure luminosity evolution with $`Q=2.74`$ gives a similar fit to the luminosity function data points, but predicts a tail on the $`n(z)`$ distribution which extends beyond $`z=5`$, and a star-forming radio background that is slightly too high ($`4.1\times 10^{23}\mathrm{Wm}^2\mathrm{ster}^1\mathrm{Hz}^1`$).

The 6 parameters used here are very interdependent. As we progressed through this work, adding data and model features, we fit our model to the data numerous times. The resulting fits occured in a wide range of this parameter space, with the turnover at high redshift sometimes occuring in $`f(z)`$ (luminosity evolution) and sometimes occuring in $`g(z)`$ (number density evolution). Some fits had a steeper increase in $`f(z)`$ and a decrease in $`g(z)`$, others had a shallower increase in $`f(z)`$ and an increase in $`g(z)`$. This reminds us that the 6 parameters are degenerate, and most likely a different parameterization with fewer free parameters could describe the data as well. Finding a new parameterization, however, would go beyond the scope of this present work, given the limited sample size currently available. Although the shapes of $`f(z)`$ and $`g(z)`$ varied greatly between different fits, all of the fits predicted generally similar shapes for the observables (the $`n(z)`$ distribution and the luminosity function), and for the predicted global star formation history (see §4.3).

## 4 Star Formation History

Now that we have a model of the evolving luminosity function, we can use it to determine the star formation history. First we describe the relationship between star formation rate and radio luminosity (§4.1), then calculate the star formation history directly from the data with minimal model dependence (§4.2), and finally calculate it from the model (§4.3).

### 4.1 Star formation rate from radio luminosity

For an individual star-forming galaxy, the star formation rate is directly proportional to its radio luminosity (Condon, 1992):

$$\mathrm{SFR}=Q\left(\frac{L_\nu /(\mathrm{WHz}^1)}{5.3\times 10^{21}\left(\nu /\mathrm{GHz}\right)^{0.8}+5.5\times 10^{20}\left(\nu /\mathrm{GHz}\right)^{0.1}}\right)\mathrm{M}_{}\mathrm{yr}^1$$
(15)
Condon (1992) derives this relation by calculating the synchrotron radio emission from supernova remnants (the first term in the denominator) and the thermal radio emission from HII regions (the second term). The spectral index of 0.8 is typical for the non-thermal component of a radio source at 1.4 GHz. This relation is derived purely from radio considerations. Cram et al. (1998), however, compares this relation to H$`\alpha `$ studies, and finds that they give similar star formation rates for local individual galaxies, with the except of galaxies with extremely large star formation rates. (Hopkins et al., 2000 has found SFR-dependent dust corrections that shift the optical results to match 15 and radio observations.)

Both the thermal and non-thermal components of the radio expression are proportional to the formation rate of high-mass stars ($`M>5M_{}`$) which produce supernova and large HII regions, so the factor $`Q`$ is included to account for the mass of all stars in the interval $`0.1100M_{}`$,

$$Q=\frac{_{0.1M_{}}^{100M_{}}M\psi (M)𝑑M}{_{5M_{}}^{100M_{}}M\psi (M)𝑑M},$$
(16)
where $`\psi (M)M^x`$ is the initial mass function (IMF). We have assumed throughout a Salpeter IMF ($`x=2.35`$), for which $`Q=5.5`$. If an upper limit of $`125M_{}`$ is used, then $`Q=5.4`$. If we use a range of mass $`0.25100M_{}`$, as suggested by Gould, Bahcall, & Flynn (1996) then $`Q=3.9`$. We will use $`Q=5.5`$ in the following.

Condon’s relationship (eq. 15) uses the emitted source luminosity at a frequency of 1.4 GHz, and thus the corrections given in eq. 3 must be applied. We should also consider whether there are other ways in which the connection between SFR and radio luminosity might evolve with redshift. At 1.4 GHz, the thermal term in eq. 15 is much smaller than the synchrotron term, so evolution in the thermal term will have little effect. In the synchrotron term, the dependence of the emitted flux on the supernova environment is weak (Condon, 1992), so little evolution is expected. However, at high redshifts, relativistic electrons may experience significant inverse-Compton cooling from the intense FIR energy density or the cosmic microwave background. Another effect that might cause significant evolution in eq. 15 is an evolving IMF, entering through the factor $`Q`$. In active starbursts, the IMF may be weighted to high-mass stars (Elmegreen, 1999), which would result in a smaller value of $`Q`$. However, the smallest $`Q`$ is unity (when virtually all mass occurs in high-mass stars), so the strongest decrease in our calculated star formation history from a radical change in the IMF would be roughly a factor of five. Note that evolution of the IMF would affect optical estimates of the star formation rate as well. In the following calculations we assume that eq. 15 does not evolve.

To determine the star formation rate per comoving volume, we simply substitute the radio luminosity density (such as eq. 13 or 18) for the source luminosity $`L_\nu `$ in eq. 15, giving

$$\mathrm{\Psi }(z)=Q\left(\frac{\rho _{e,1.4}(z)}{4.6\times 10^{21}\mathrm{WHz}^1\mathrm{Mpc}^3}\right)\mathrm{M}_{}\mathrm{yr}^1\mathrm{Mpc}^3$$
(17)
where 1.4$`\mathrm{GHz}`$ is used in the denominator of eq. 15, as all data have already been converted to 1.4 GHz in the rest frame.

### 4.2 Star formation history estimated from the data

We now calculate the star formation history directly from the survey data in §2 by using the luminosity density of the detected sources. For each redshift bin ($`z_{min}<z<z_{max}`$), the luminosity density is

$$\rho _{e,1.4}(z)=\mathrm{\Sigma }_i\frac{L_iB_iC(z)}{V_c(z_{min},z_{max}(L_i))}$$
(18)
where $`B_i`$ is the surface density given in eq. 1, $`C(z)`$ is a correction for faint sources described below, $`z_{max}(L_i)`$ is the largest $`z`$ in the bin for which the luminosity of the source $`L_i`$ was above detection limit of survey, and $`V_c`$ is the comoving volume given in eq. 5. This luminosity density can then be used in eq. 17 to find the evolving star formation density. Without the correction factor $`C(z)`$, the luminosity density includes only individual sources brighter than the flux limit of the survey. It does not include the luminosity density of sources too faint to be detected individually, and so it clearly underestimates the star formation rate (but calculations without $`C(z)`$ have the advantage of being independent of our evolution model and provide a lower limit).

To account for these faint sources, we use the evolving luminosity function found in §3.3. Figure 8 illustrates the faint source correction, using a redshift of 1.6 as an example. The integral under the curve in the figure is proportional to the total luminosity density. The flux limit of the survey, however, only allows detection of individual sources above a certain luminosity, i.e. in the cross-hash area. The faint source correction factor $`C(z)`$ in eq. 18 would then be the ratio of the total area to the cross-hash area. We have argued, however, that the slope of the source number counts changes below about 1$`\mu \mathrm{Jy}`$ (Haarsma & Partridge, 1998), so that in fact most of the sources will occupy only the hashed area of fig. 8. Thus, a more realistic correction to eq. 18 is the ratio of the hashed region to the cross-hashed region, i.e. the ratio of the luminosity density due to sources brighter than 1$`\mu \mathrm{Jy}`$ to the luminosity density from sources brighter than the flux limit of the survey,

$$C(z)=\frac{_{L(S_{lim},z)}^{inf}L_{e,1.4}\varphi (L_{e,1.4},z)d\mathrm{log}_{10}L}{_{L(1\mu \mathrm{Jy},z)}^{inf}L_{e,1.4}\varphi (L_{e,1.4},z)d\mathrm{log}_{10}L}$$
(19)
A list of corrections for several redshifts is given in Table 5. Note that if the slope of the number counts of radio sources were assumed to stay the same below 1$`\mu \mathrm{Jy}`$, these corrections would be even larger, and so would the calculated star formation density.

We calculated the star formation density using eqs. 17, 18, 19 for the lower, middle, and upper samples described in §2.5; the results are shown in Figures 9 and 10 and listed in Table 6. Recall that the “lower” sample includes only sources with spectroscopic redshifts and definite identifications with spirals, irregulars, or mergers, and thus is the minimum amount of star formation activity consistent with the data. The “middle” value includes some sources with ambiguous identifications and rough photometric redshifts, but is our best guess at the total radio-selected star-forming population. The “upper” sample includes all the sources and is the maximum possible star formation activity allowed by the radio data. In Figure 9, the points are values calculated from the middle sample, plotted at the average of the source redshifts in the bin. The vertical error bars are either the limits from the lower and upper samples, or the Poisson errors (from the number of galaxies per bin), whichever is larger (in most cases, the lower limits and upper limits are larger than the Poisson errors). To make the lower sample a true lower limit, we did not include the faint source correction $`C(z)`$. Thus the lower limit is for only those sources clearly identified with star forming systems and having spectroscopic redshifts, with no allowance made for evolution of the luminosity function or for sources below the survey flux limits. This is surely a gross underestimate of the true star formation density value, since faint sources, unidentified sources, and sources without spectroscopic redshifts are all missing, but it provides a firm lower limit to the true star-formation rate. The “middle” and “upper” samples do have the correction $`C(z)`$ for faint sources applied.

Mobasher et al. (1999) have done a similar calculation of star formation density from a survey of faint radio sources. They find no evidence for evolution from $`z=0`$ to $`z=1`$, and a decrease in star formation density above $`z=0.3`$. Their results, however, are based on a radio survey sample that is only 50% complete, and thus their results are highly dependent on assumptions made when correcting for incompleteness. Similarly, their optical identifications and spectroscopic redshifts are much less complete than ours. Finally, the radio surveys we use extend to much fainter flux densities where star formation is more likely to dominate the radio emission. The fainter flux limit also allows us to detect more high redshift sources. Thus, we believe our results for the star formation density are more reliable than those of Mobasher et al. (1999).

### 4.3 Star formation history predicted by model

The star formation history can also be determined directly from the evolution model found in §3.3. We simply calculate the luminosity density emitted at 1.4 GHz (eq. 13), and use eq. 17 to find the star formation density. The resulting star formation density prediction is the curve plotted in Figures 9 and 10.

The model curve in Figures 9 and 10 falls somewhat below the averaged data points (heavy crosses). Note that the model curve was not fitted to these averaged data points (which were calculated using equations 17, 18, 19), but rather the model was found from the evolving luminosity function alone (using equations 13 and 17). The evolving luminosity function, in turn, was fitted to the $`n(z)`$ data, the luminosity function data, and the radio background (see §3.3). As discussed in §3.3, our model fits these data well but not perfectly, so small differences between the model prediction and the star formation data points are not unreasonable.

This sort of calculation was previously done by Cram (1998), using the Condon (1989) luminosity function and Condon (1984a) evolution model. Please note the typographical error in eq. 2 of Cram (1998), which differs by a factor of 28.2 from our eq. 15; the correct values were used in their calculations (L. Cram, private communication). We agree with Cram’s calculation of 0.026$`\mathrm{M}_{}\mathrm{yr}^1\mathrm{Mpc}^3`$ for the local star formation density, and calculation of star formation history from Condon’s early model.

Note that these methods and results are much improved over our very preliminary work (Haarsma & Partridge, 1999), which assumed that the majority of detected faint radio sources lie at the redshift of peak star formation activity. In fact, the peak of the observed redshift distribution (Figures 3-5) is at a lower redshift than the peak star formation activity (Figures 9 and 10), due to cosmological factors, such as the dependence of the comoving volume on redshift.

## 5 Discussion and Conclusions

Radio wavelength determinations of the Universal star formation history have the important advantage of being independent of the dust content of galaxies. Additionally, it is possible to cull relatively clean samples of star-forming objects using radio properties such as spectral index, morphology, and variability. Our results are shown in Figure 10, overlaid with the star formation histories found in several other studies. In Figure 10, the thick crosses show the star formation density of our middle sample, with error bars indicating Poisson errors. Although it is possible that our middle sample may include some low-luminosity AGN (Seyferts, etc.), our careful definition of the sample (§2.5) and the large fraction of sources with clear optical identifications reduces this contamination. The thin curve is our model prediction, found in §4.3 by fitting to the luminosity function, redshift distribution, and radio background (not to the thick crosses). The thick lines indicate our lower and upper limits on the star formation density, which are calculated in §4.2 using samples defined in §2.5. The lower limit is very firm, since it includes only those sources with spectroscopic redshifts and identifications with spirals, irregulars and mergers, and does not include the star formation in galaxies fainter than the survey detection limit. The upper limit is also firm, since it includes all detected radio sources, even those not associated with star forming galaxies, but is more uncertain than the lower limit since it includes sources without spectroscopic redshifts.

If we were to assume a different cosmology, our results would change somewhat. The values of $`H_0`$, $`\mathrm{\Omega }_m`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ affect the calculation of distance from redshift, and luminosity from flux density. The star formation density is proportional to luminosity/volume, so it is inversely proportional to distance and directly proportional to $`H_0`$. If we had assumed $`H_0`$=100$`\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ instead of $`H_0`$=50$`\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, our data points and model for the star formation density would be twice as large. The values of $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ also affect the distance measurement. As an extreme example, in a nearly empty but flat universe ($`\mathrm{\Omega }_m`$=0.1, $`\mathrm{\Omega }_\mathrm{\Lambda }`$=0.9), the distance to a $`z=1`$ object is about 1.5 times larger than the distance to it in our assumed cosmology ($`\mathrm{\Omega }_m`$=1, $`\mathrm{\Omega }_\mathrm{\Lambda }`$=0), and thus the star formation density would be about 2/3 of our listed value. Note that all data points in Figures 1 and 10 depend similarly on these cosmological parameters.

At low redshifts, we agree with the findings of many studies that star formation density increases rapidly from the local universe to $`z=1`$. We disagree with Cowie, Songaila, & Barger (1999) who find a gradual (rather than steep) increase from $`z=0`$ to $`z=1`$, and with Mobasher et al. (1999) who find a decrease in star formation density $`z=0`$ to $`z=1`$. Our firm lower limit is significantly higher $`z<1`$ than the extinction corrected optical results of Lilly et al. (1996) at $`z0.7`$, indicating that some star formation has been obscured by dust. Our “middle sample” data points fall above all the optical and ultraviolet studies shown, indicating that these studies have probably missed some star formation by underestimating the dust extinction (see Hopkins et al., 2000 for SFR-dependent dust corrections that bring these data more into agreement with our radio results). Our results are similar to the star formation density from the Infrared Space Observatory (ISO) survey of the HDF (Rowan-Robinson et al., 1997).

At redshifts above $`z=1`$, we cannot draw strong conclusions. There are few sources with spectroscopic redshifts in this range, so our calculations are based in large part on less secure photometric redshift estimates and random redshift assignments for the very red objects (§2.3). Our assumption that the relationship between radio luminosity and star formation rate does not evolve also becomes less sure as we move to higher redshift (see discussion in §4.1). The Initial Mass Function may also be evolving, although this would affect optical and ultraviolet estimates of star formation history as well. Finally, the faint source corrections (eq. 19, Table 5) become larger at high redshift, and thus depend more strongly on the assumed shape of the luminosity function. In fact, at redshifts above 1.5, the current radio survey limits only probe the extreme end of the luminosity function (SFR per galaxy $`>1000M_{}`$). Deeper surveys are needed to detect radio counterparts to typical high-redshift optical objects, e.g. Lyman break galaxies (for instance, the predicted radio fluxes for even the most luminous Lyman break galaxies in the HDF are only a few $`\mu \mathrm{Jy}`$ at 1.4 GHz, e.g. Meurer, Heckman, & Calzetti, 1999). Planned improvements to the VLA will allow future surveys to reach this sensitivity.

Still, our calculations at high redshift show that if even a small number of star-forming radio sources exist beyond $`z1.5`$, they would indicate a large, optically-hidden fraction of star-formation density. In particular, the population of radio sources with faint, red optical counterparts may be dust-enshrouded (Richards et al., 1999; Barger et al., 2000; Waddington et al., 1999, §2.4), and hence missed even in the deepest optical and ultraviolet studies. Deeper high resolution radio observations, accompanied by close to complete spectroscopic identifications, are needed to determine accurately the amount of “hidden” star formation in the early Universe.

## Acknowledgments

D.B.H. and R.B.P. acknowledge the support of NSF grant AST 96-16971 and the Keck Northeast Astronomy Consortium. D.B.H. acknowledges the support of a Cottrell College Science Award from Research Corporation. E.A.R. acknowledges support from Hubble Fellowship Grant HG-0123.01-99A, and from STScI which is operated by AURA under NASA contract NAS 5-26555. R.A.W. acknowledges support from AST-8821016 and AST-9802963 from the National Science Foundation and NASA grants AR-6948.04-96A and GO-7452.01-96A from STScI under NASA contract NAS5-26555. We are grateful to James Lowenthal for providing us with spectroscopic redshifts in advance of publication. We thank Ian Waddington, Ed Fomalont, Ken Kellermann, and Andrew Hopkins for helpful comments and discussions, and our anonymous referee for detailed and helpful comments on the manuscript.


format: text

text:
# Relative Entropy: Free Energy Associated with Equilibrium Fluctuations and Nonequilibrium Deviations

## I Introduction

Entropy is the most important concept in both statistical mechanics and information theory. In the latter , entropy is a quantity associated with any discrete probability distribution $`\{p_k\}`$

$$S[\{p_k\}]=\underset{k}{}p_k\mathrm{ln}p_k$$
(1)
and its generalization to continuous probability distributions is straightforward. In statistical mechanics, the same Eq. 1 gives the entropy, in units of $`k_B`$, for a canonical ensemble of a molecular system at constant temperature. This is know as Gibbs entropy. The probability of a molecular configuration, say $`\{x_k\}`$, is related to its free energy $`F(\{x_k\})`$ according to Boltzmann’s law $`P(\{x_k\})e^{F(\{x_k\})/k_BT}`$ where $`k_B`$ is the Boltzmann’s constant and $`T`$ is the temperature in Kelvin . It is generally accepted, however, that the entropy in information theory and that in statistical mechanics, though share the same name, are not the same.

This paper is to demonstrate a second connection between the two fields, both are based on the theory of probability. There is another important concept in information theory called relative entropy $`H[\{p_k\}|\{p_k^{}\}]`$. It is associated with probability distributions $`\{p_k\}`$ and $`\{p_k^{}\}`$:

$$H[\{p_k\}|\{p_k^{}\}]=\underset{k}{}p_k\mathrm{ln}\left(\frac{p_k}{p_k^{}}\right).$$
(2)
The relative entropy has many important mathematical properties, for example, it is positive and equals zero if and only if $`p_k=p_k^{}`$, and furthermore, it is a convex function of $`p_k`$.

Until now, however, the relative entropy has not found naturally a physical interpretation in statistical mechanics, despite its important role, as a mathematical device, in the stability analysis of master equations  and Fokker-Planck equations . In this note, we demonstrate, using two specific examples, that the relative entropy is in fact the free energy associated with isothermal equilibrium fluctuations and a generalized free energy associated with transient nonequilibrium deviations. The two examples we use are from the recently developed stochastic theory of macromolecular mechanics , however, the generalization of our result is straightforward, and in fact is known in mathematical literature .

The statistical mechanical system we discuss is a single macromolecule in aqueous solution at a constant temperature $`T`$. Interest on such system is motivated by the recent experimental studies on single biological molecules . We point out that such systems provide a unique type of nonequilibrium problems in which the momentum distribution is in rapid equilibrium due to collisions with the solvent molecules. Hence the nonequilibrium problem is only for the conformational, stochastic dynamics of the macromolecule. A separation of the time scales for momentum and conformation (over-damped mechanics) is assumed, which leads to the Smoluchowski equation.

The Smoluchowski approach to nonequilibrium statistical mechanics of single macromolecules in aqueous solution is in parallel and complementary to the approach based on Boltzmann equation for gases and liquids. Both approaches are based on Newtonian mechanics, but both invoke a stochastic element a priori in dealing with collisions . There is an extensive literature on the nonequilibrium statistical mechanics, including studies on relative entropy, based on the Boltzmann’s framework . This work illustrates the approach based on the Smoluchowski’s framework , which is simpler, conceptually straightforward, and applicable to biomolecular applications. We note that the rate of uncompensated heat in  seems to correspond to the entropy production rate in the Smoluchowski’s framework. The precise and concrete relation between these two quantities remains to be established .

## II Equilibrium Fluctuations of a Polymer Chain

Let’s consider a one-dimensional polymer chain with $`N`$ identical subunits. Each subunit is an elastic element with free energy function (potential of mean force) $`\varphi (x)`$. Let $`x_1,x_2,\mathrm{},x_{N1}`$ be the junctions between successive subunits. We assume that the end $`x_0=0`$ is anchored and the end $`x_N`$ is freely fluctuating. This model is motivated by the mechanical studies on giant muscle protein titin . According to Boltzmann’s law, the joint probability for all $`\{x_k\}`$ is

$$P_{eq}(\{x_k\})=Z^1\mathrm{exp}\left[\frac{\varphi (x_1x_0)+\varphi (x_2x_1)+\mathrm{}+\varphi (x_Nx_{N1})}{k_BT}\right]$$
(3)
where

$`Z=\left({\displaystyle _0^{\mathrm{}}}\mathrm{exp}\left[{\displaystyle \frac{\varphi (z)}{k_BT}}\right]𝑑z\right)^N.`$

The partition function of the equilibrium state of such system can be obtained analytically using Laplace transform . The theory of large deviation of level-1 serves as its rigorous mathematical foundation .

In a laboratory, usually only $`x_N`$ is observable . For this case, we obtain the marginal distribution of (3)

$$P_{eq}(x_N)=_0^{\mathrm{}}\mathrm{}_0^{\mathrm{}}P_{eq}(\{x_i\})𝑑x_1\mathrm{}𝑑x_{N1}.$$
(4)
Combining the probability distribution in Eq. 4 with Boltzmann’s law, we obtain a free energy function for the entire polymer under the condition that the end of the chain is at $`x_N`$

$`F(x_N)=k_BT\mathrm{ln}P_{eq}(x_N).`$

In a completely parallel fashion, one can obtain the free energy function $`F(x_i,x_j)`$ from a marginal distribution of (3).

A different type of laboratory measurements is to obtain the density for subunit length. Optical spectroscopy is sensitive to the length of subunits, hence it provides a measurement on the density function

$$\nu (x)=\frac{1}{N}\underset{k=1}{\overset{N}{}}\delta \left(xx_k+x_{k1}\right),$$
(5)
which is known as the empirical measure in the theory of large deviation of level-2 .

Clearly, the function $`\nu (x)`$ fluctuates in an equilibrium state since $`(x_kx_{k1})`$ fluctuates. It has an expectation

$`E\left[\nu (x)\right]=p(x)={\displaystyle \frac{e^{\varphi (x)/k_BT}}{_0^{\mathrm{}}e^{\varphi (x)/k_BT}𝑑x}}.`$

For large $`N`$, the fluctuations of $`\nu (x)`$ around $`p(x)`$ are so small that one seldomly considers their existence. Nevertheless, there is a free energy associated with each $`\nu (x)`$, and we now show that

$$F\left[\nu (x)\right]F\left[p(x)\right]=Nk_BT_0^{\mathrm{}}\nu (x)\mathrm{ln}\left(\frac{\nu (x)}{p(x)}\right)𝑑x.$$
(6)
This is in fact a mathematical result for large deviation of level-2 . We give only a heuristic proof below and leave the rigorous treatment to the mathematical literature.

Let’s denote

$`\nu _n={\displaystyle _{n\delta }^{(n+1)\delta }}\nu (x)𝑑x,\text{and}p_n={\displaystyle _{n\delta }^{(n+1)\delta }}p(x)𝑑x.`$

Because all the subunits are independent, $`\nu _n`$ is a multinomial distribution

$$P(\nu _1,\nu _2,\mathrm{}\nu _m)=\frac{N!}{(N\nu _1)!(N\nu _2)!\mathrm{}(N\nu _m)!}p_1^{N\nu _1}p_1^{N\nu _2}\mathrm{}p_1^{N\nu _m}.$$
(7)
According to Boltzmann’s law, the free energy difference between configuration $`\{\nu _n\}`$ and $`\{p_n\}`$ is

$`F\left[\{\nu _n\}\right]F\left[\{p_n\}\right]`$ $`=`$ $`k_BT\mathrm{ln}\left({\displaystyle \frac{P\left(\{\nu _n\}\right)}{P\left(\{p_n\}\right)}}\right)`$
$``$ $`Nk_BT{\displaystyle \underset{n=1}{\overset{m}{}}}\nu _n\mathrm{ln}\left({\displaystyle \frac{\nu _n}{p_n}}\right)`$
$``$ $`Nk_BT{\displaystyle _0^{\mathrm{}}}\nu (x)\mathrm{ln}\left({\displaystyle \frac{\nu (x)}{p(x)}}\right)𝑑x(\delta 0)`$
Therefore, the relative entropy in Eq. 6 is the free energy difference between the distribution $`\{\nu (x)\}`$ and its (equilibrium) expectation $`\{p(x)\}`$. Relative entropy is the free energy associated with a fluctuating density at equilibrium !

## III Nonequilibrium Relaxation of a Polymer Chain

Our second example extends the concept of free energy beyond an equilibrium state, and reveals its central role in the transient, isothermal, relaxation processes to equilibria. The dynamic model for the polymer chain in an aqueous solution at constant temperature $`T`$ is a Smoluchowski equation :

$$\frac{P(𝒙,t)}{t}=\frac{k_BT}{\eta }^2P+\frac{1}{\eta }\left[U(𝒙)P\right]$$
(8)
where $`\eta `$ is a frictional coefficient, $`𝒙`$ = $`(x_1,x_2,\mathrm{},x_N)`$, $``$ = $`(\frac{}{x_1},\frac{}{x_2},\mathrm{},\frac{}{x_N})`$, and

$$U(𝒙)=\varphi (x_1x_0)+\varphi (x_2x_1)+\mathrm{}+\varphi (x_Nx_{N1}).$$
(9)
It is easy to verify that Boltzmann’s distribution $`P_{eq}(𝒙)`$ in Eq. 3 is the stationary solution to the Eq. 8. In fact, the steady-state solution of (8) defines a stationary, time-reversible, stochastic process with equilibrium fluctuations .

How is an arbitrary distribution $`P(𝒙)P_{eq}(𝒙)`$ changing with time and approaching to $`P_{eq}(𝒙)`$? We now show that a free energy functional can be introduced, and it is in fact the relative entropy. Let’s define a $`\mathrm{\Psi }`$-function

$`\mathrm{\Psi }[P(𝒙)]`$ $`=`$ $`{\displaystyle \left(U(𝒙)P(𝒙)+k_BTP(𝒙)\mathrm{ln}P(𝒙)\right)𝑑𝒙}`$ (10)
$`=`$ $`k_BT\mathrm{ln}Z+k_BT{\displaystyle P(𝒙)\mathrm{ln}\left(\frac{P(𝒙)}{P_{eq}(𝒙)}\right)𝑑𝒙}.`$ (11)
The first term is the Helmhotz’s free energy of the entire polymer chain in its equilibrium state, and the second term should be interpreted as the free energy difference between the arbitrary distribution $`P(𝒙)`$ and the equilibrium distribution $`P_{eq}(𝒙)`$.

The mathematical properties of relative entropy immediately lead to the following statement: The system reaches equilibrium if and only if the free energy functional (the $`\mathrm{\Psi }`$-function) is minimized. Furthermore, if $`P(𝒙,t)`$ changes with $`t`$ in a transient process according to Eq. 8, then

$$\dot{\mathrm{\Psi }}[P(𝒙,t)]=J(𝒙,t)\mathrm{\Phi }(𝒙,t)𝑑𝒙0$$
(12)
where

$`𝑱(𝒙,t)={\displaystyle \frac{k_BT}{\eta }}P(𝒙,t){\displaystyle \frac{1}{\eta }}U(𝒙)P(𝒙,t)`$

and

$`𝚽(𝒙,t)=k_BT\mathrm{ln}P(𝒙,t)U(𝒙)`$

are fluxes and forces, which are both zero at equilibrium . Furthermore, $`𝑱=P𝚽`$ and Eq. 12 is related to the entropy production rate . Eq. 12 immediately leads to a second statement: The dynamic of $`P(𝒙,t)`$ follows a path of decreasing free energy. The free energy functional in (10) is a Lyapunov function  for the stochastic dynamics of the polymer. The relative entropy is the free energy difference between an arbitrary distribution and the equilibrium distribution. It is associated with the nonequilibrium deviation from the equilibrium state. However, it is interesting to note that the dynamics does not follow the steepest descent of the free energy functional. The interpretation and significance of this observation are not clear to us at present time. We also note that Eq. 12 corresponds to the $`H`$-theorem in Boltzmann’s framework.

It is of course not a coincidence that the relative entropy appears as a free energy difference in both the equilibrium and nonequilibrium situations. Onsager  has pointed out that the force driving the nonequilibrium relaxations is in fact the same force causing the equilibrium fluctuations to return to its mean. Our result, therefore, firmly relates the force to a free energy difference in terms of the relative entropy.

In summary, we have shown that the relative entropy in information theory has a natural physical meaning in equilibrium and nonequilibrium statistical mechanics. It is in fact the free energy difference associated with the equilibrium fluctuations of a density function, a result known in the theory of large deviations. In nonequilibrium systems, it is the free energy difference between an arbitrary distribution and the equilibrium distribution. Again the latter result is known in the theory of Markov processes. What is novel of the present note is to identify the relative entropy with the Helmholtz’s free energy for isothermal systems. Conversely, the mathematical theorems mentioned above become an integral part of the statistical physics of macromolecules.


text:
# Paramagnetic Meissner Effect in Multiply-Connected Superconductors

\[

## Abstract

We have measured a paramagnetic Meissner effect in $`\mathrm{Nb}\mathrm{Al}_2\mathrm{O}_3\mathrm{Nb}`$ Josephson junction arrays using a scanning SQUID microscope. The arrays exhibit diamagnetism for some cooling fields and paramagnetism for other cooling fields. The measured mean magnetization is always small, $`\mathrm{\Phi }_{\mathrm{tot}}\mathrm{\Phi }_{\mathrm{ext}}<0.3\mathrm{\Phi }_0`$ (in terms of flux per unit cell of the array, where $`\mathrm{\Phi }_0`$ is the flux quantum) for the range of cooling fields investigated ($`12\mathrm{\Phi }_0`$ to $`12\mathrm{\Phi }_0`$). We demonstrate that a new model of magnetic screening, valid for multiply-connected superconductors, reproduces all of the essential features of paramagnetism that we observe and that no exotic mechanism, such as $`d`$-wave superconductivity, is needed for paramagnetism.

\]

A paramagnetic response was first observed in BSCCO  and subsequently reported in many other high $`T_c`$ materials. This response is in striking contrast to the diamagnetic response which is a hallmark of superconductivity. It has been argued that the paramagnetic Meissner effect (PME) results from a $`d`$-wave symmetry in the BSSCO order parameter. The observation of PME in niobium and aluminum, known to be $`s`$-wave superconductors, questioned the validity of the $`d`$-wave argument.

There are several proposed theories for the PME, although no experiment has been yet able to provide conclusive evidence for one or another; and there is disagreement as to the validity of the $`d`$-wave explanation. In high $`T_c`$ materials one possible explanation is the existence of $`\pi `$-junctions which result from the misalignment of grains in materials with $`d`$-wave symmetry.  In the conventional materials, which are known to be $`s`$-wave, PME has been attributed to different mechanisms causing non-equilibrium flux configurations. A non-equilibrium flux state may result from flux-compression, the similar giant vortex state, pinning randomness, or a surface barrier at the edge of the sample.

To create a controlled experiment we have chosen to look at Josephson-junction arrays using a scanning SQUID microscope (SSM), using a technique similar to that of Kirtley et al. The SSM provides spatially resolved magnetization images of the sample. Arrays eliminate many of the complications in studying the underlying cause of the PME. Since our arrays are lithographically fabricated from $`\mathrm{Nb}`$ and $`\mathrm{Al}_2\mathrm{O}_3`$, we know that there are neither $`\pi `$-junctions nor significant randomness. Additionally, our sample is uniformly attached to the thermal bath, eliminating thermal gradients across the sample during cooling. We find that arrays can be either paramagnetic or diamagnetic depending upon the cooling field. Additionally, we propose a new model, valid for multiply-connected supercondutors, which predicts a paramagnetic response and agrees with our scanning SQUID measurements.

The arrays consist of a square array of niobium crosses in two layers, with a unit cell size of $`a=46\mu \mathrm{m}`$. Junctions are formed at the cross overlap. The calculated self-inductance of each loop of four junctions is $`L=64\mathrm{p}\mathrm{H}`$ and the measured critical current density is $`J_c=600\mathrm{A}/\mathrm{cm}^2`$ at $`4.2\mathrm{K}`$, with a junction area of $`5\times 5\mu \mathrm{m}^2`$. The entire array consists of $`30\times 100`$ junctions. The unit-cell dimensions compare to the typical grain size seen in BSCCO samples  which exhibit PME, to which similar fields are applied. The array is placed in an SSM sample stage, sitting at the end of a sapphire rod of $`10\mathrm{m}\mathrm{m}`$ diameter, around which a solenoid has been wound. Before any data is collected we first remove any flux trapped within the pick up loop of the SQUID to avoid interaction with the sample. Additionally, we measure and correct for any background field present.

We cooled our arrays to 4.2K in fixed external magnetic field, and then used the SSM to measure the magnetic flux threading the SQUID from the sample, $`\mathrm{\Phi }_{\mathrm{tot}}`$, as a function of position. Our spatial resolution was limited by our SQUID-sample separation distance, which we estimate to be between $`40`$ and $`60\mu \mathrm{m}`$. Data were taken every $`5\mu \mathrm{m}`$ in the $`x`$ direction and every $`50\mu \mathrm{m}`$ in the $`y`$ direction, and converted from flux through the SQUID to normalized flux per unit cell. We measure the external flux, $`\mathrm{\Phi }_{\mathrm{ext}}`$, directly by warming the sample above $`T_c`$ and imaging. From the data we calculate both a local normalized magnetization $`(\mathrm{\Phi }_{\mathrm{tot}}\mathrm{\Phi }_{\mathrm{ext}})/\mathrm{\Phi }_0`$ and an average overall magnetization $`\mathrm{\Phi }_{\mathrm{tot}}\mathrm{\Phi }_{\mathrm{ext}}/\mathrm{\Phi }_0`$. We were unable to resolve individual unit cells in the array,  but we were still able to resolve large scale flux distributions.

Our measurement and analysis proceeded in a fashion similar to Ref. . We first make an image of the magnetization of the array (a typical image is shown in Fig. 1(a)). We then make a magnetization histogram and compute the mean magnetization. Fig. 1(b) shows the magnetization histogram generated from the image in Fig. 1(a). The histogram is computed by looking at the lower third of the array and counting the number of points at a given magnetization. We find that the array may be either diamagnetic or paramagnetic, but surprisingly, appears to be preferentially paramagnetic. Additionally, we always observe diamagnetic screening currents around the outside edges of the sample. The array shown in Fig. 1(a) was cooled in an external field of $`4.8\mathrm{\Phi }_0`$ and clearly is paramagnetic (i.e. it is darker than the background). We find a mean paramagnetic response of $`\mathrm{\Phi }_{\mathrm{tot}}\mathrm{\Phi }_{\mathrm{ext}}=(0.063\pm 0.005)\mathrm{\Phi }_0`$ for this cooling field.

The data plotted in Fig. 2 show the mean magnetization versus external cooling field. Each data point represents a single field-cooled average measurement. The error in the mean values plotted, such as discussed above, is less than the size of each data point as drawn. The sign of the magnetization is reproducible, and the magnitude is reproducible to within a factor of two when the sample is warmed above 10K and recooled in the same field, indicating that our external field is consistent to within $`\pm 0.05\mathrm{\Phi }_0`$ The external field has a spatial variation of $`\pm 0.5\mathrm{\Phi }_0`$, which is limited by the gradient of the background field in our probe. Although there is some variation from one cooling to the next in the actual magnetic images for the same external field, the paramagnetism is remarkably reproducible.

We may consider a single loop of the array as a jumping off point for modeling the entire array’s magnetization. Originally used to describe the observed magnetization of a single loop with a single Josephson-junction, Sigrist and Rice  later used this model to argue that $`\pi `$-junctions could cause PME in $`d`$-wave superconductors. The single-loop model gives a relationship between the external flux applied to a four-junction loop and the total measured flux

$$\mathrm{\Phi }_{\mathrm{tot}}=\mathrm{\Phi }_{\mathrm{ext}}+LI_c\mathrm{sin}\left(\frac{\pi }{2}n\frac{\pi }{2}\frac{\mathrm{\Phi }_{\mathrm{tot}}}{\mathrm{\Phi }_0}\right).$$
(1)
Here $`I_c`$ is the critical current of a single junction. For all values of $`\frac{LI_c}{\mathrm{\Phi }_0}>\frac{2}{\pi }`$, $`\mathrm{\Phi }_{\mathrm{tot}}`$ is a multi-valued function of $`\mathrm{\Phi }_{\mathrm{ext}}`$ depending on an integer $`n`$. We note here that Ref.  used a model with a $`\frac{LI_c}{\mathrm{\Phi }_0}`$ slightly greater than one to demonstrate the absence of potential paramagnetism for a single loop with a single junction. Had they used a value of $`\frac{LI_c}{\mathrm{\Phi }_0}1`$ they would have noticed the potential for paramagnetism even in the absence of external flux, although such a situation is not energetically favorable. For a four junction loop, $`\beta _{\mathrm{}}=2\pi \frac{LI_c}{\mathrm{\Phi }_0}=30`$ allows for the possibility of paramagnetism (or diamagnetism) at $`\mathrm{\Phi }_{\mathrm{ext}}=0`$.

Computing the Gibbs free energy for the various possible solutions to (1) we find that for $`m<\frac{\mathrm{\Phi }_{\mathrm{ext}}}{\mathrm{\Phi }_0}<m+\frac{1}{2}`$ the single loop will be diamagnetic and for $`m+\frac{1}{2}<\frac{\mathrm{\Phi }_{\mathrm{ext}}}{\mathrm{\Phi }_0}<m+1`$ will be paramagnetic, where $`m`$ is any integer. If one applies a positive offset flux $`\mathrm{\Phi }_{\mathrm{off}}`$ to the four-junction loop (i.e. the total applied flux is $`\mathrm{\Phi }_{\mathrm{ext}}+\mathrm{\Phi }_{\mathrm{off}}`$) the solutions $`\mathrm{\Phi }_{\mathrm{tot}}\mathrm{\Phi }_{\mathrm{ext}}`$ vs. $`\mathrm{\Phi }_{\mathrm{ext}}`$ become more paramagnetic and conversely if $`\mathrm{\Phi }_{\mathrm{off}}<0`$ the solutions become more diamagnetic. If, e.g., a simple diamagnetic screening current flows around the outside edge of the superconductor (negative offset), the measured magnetization will be more diamagnetic.

A conventional superconductor is a perfect diamagnet: An applied flux causes screening currents to flow around the outside. In an ideal array of Josephson junctions, one may decompose this edge current, $`I`$, into loop currents of magnitude $`I`$, flowing around each plaquette, as shown in Fig. 3(a). In a real array, this is not likely to be the complete picture. The currents will flow only within one penetration depth in the superconductor. For our array, the penetration depth is $`900\mathrm{\AA }`$ while the wire size is $`10\mu \mathrm{m}`$ so that the interior loop currents do not cancel. This means that the array must create an energy $`\frac{1}{2}LI^2`$ for each plaquette in order to generate a screening current $`I`$ for the array. In an overlap geometry, the currents through the Josephson-junctions are uniform, so the currents through the junctions still cancel and we do not have any increase in the total Josephson energy, $`E_J`$. For an array of $`N\times N`$ junctions Fig. 3(a) costs a total energy of

$$E_a=\frac{1}{2}LI^2N^2+4N(1\mathrm{cos}\gamma )E_J$$
(2)
where $`\gamma `$ is the gauge invariant phase difference.

We propose that the array actually screens in a much different fashion. The array may generate a diamagnetic screening current around the outside perimeter if each exterior plaquette maintains a screening current. as shown schematically in Fig. 3(b). In addition to the diamagnetic screening current, we then also have a paramagnetic current, of the same magnitude, flowing just inside the array. The array must generate the $`\frac{1}{2}LI^2`$ energy only for each exterior plaquette, so the energy cost to the array in this case is

$$E_b=\frac{1}{2}LI^2(4N)+\{4N+4(N1)\}(1\mathrm{cos}\gamma )E_J.$$
(3)
It is clear from (2) and (3) that $`E_aN^2`$ while $`E_bN`$ so that for a large enough array, the conventional picture of screening will *not* be energetically favorable. Fig. 3(b) has lower energy than Fig. 3(a) when, assuming $`\gamma `$ is proportional to $`\mathrm{\Phi }_{\mathrm{tot}}`$,

$$\beta _{\mathrm{}}>\frac{4N1}{N^24N},N>2.$$
(4)
For our array the criteria is easily satisfied and we expect screening to take place as in Fig. 3(b).

In Fig. 3(a) we can treat the outer ring of Josephson-junctions as a single loop of $`4N`$ junctions. This loop will screen up to $`\frac{LI_c}{\mathrm{\Phi }_0}=4.8`$ flux quanta from the interior of the array, when $`I_c`$ is flowing around the exterior of the array. In Fig. 3(b) each exterior plaquette can screen up to $`\frac{LI_c}{\mathrm{\Phi }_0}`$ flux quanta *from itself*. This is a much larger absolute field than in the conventional case and the difference can be easily observed experimentally. We have observed, that flux will not penetrate the array until the external flux exceeds at least $`4.0\mathrm{\Phi }_0`$ per unit cell. Additionally, ac susceptibility work demonstrated that little effect was seen until the amplitude of the driving ac flux exceeded $`4.6\mathrm{\Phi }_0`$ per unit cell of the array.

The interior plaquettes see the external field, a diamagnetic flux from the extreme outside loop of the array, and a paramagnetic flux from just inside the array. Because the loop generating the paramagnetic flux is closer, the interior loops will see a net paramagnetic flux in addition to the external field. Considered individually, the interior plaquettes obey a modified Eq. (1),

$$\mathrm{\Phi }_{\mathrm{tot}}=\mathrm{\Phi }_{\mathrm{ext}}+\mathrm{\Phi }_{\mathrm{scr}}+LI_c\mathrm{sin}\left(\frac{\pi }{2}n\frac{\pi }{2}\frac{\mathrm{\Phi }_{\mathrm{tot}}}{\mathrm{\Phi }_0}\right)$$
(5)
where $`\mathrm{\Phi }_{\mathrm{scr}}`$ is the flux seen by the inner plaquette due to the screening currents produced by the exterior plaquettes. When $`\mathrm{\Phi }_{\mathrm{scr}}>0`$ the interior plaquette sees a paramagnetic offset and its magnetization will be more paramagnetic. This causes the array to be paramagnetic more often than diamagnetic in agreement with experiment.

If we consider only the mutual inductance between a test plaquette in the interior of the array and the exterior plaquettes which provide the screening, we can make a crude estimate of the paramagnetic response. The maximum current that the junctions can sustain is the critical current, $`I_c`$; this sets the upper bound for the magnetization offset. We have numerically computed the flux induced from the screening currents in Fig. 3(b) onto the interior plaquettes for our $`30\times 100`$ junction array. When the screening currents are $`I_c`$, the maximum flux is induced, and the *minimum* $`\mathrm{\Phi }_{\mathrm{scr}}`$ we compute is $`0.15\mathrm{\Phi }_0`$ per unit cell, near the center of the array. This means that the center plaquette should have its magnetization shifted upward by $`0.15\mathrm{\Phi }_0`$. Plaquettes nearer to the exterior see $`\mathrm{\Phi }_{\mathrm{scr}}`$ values as large as several $`\mathrm{\Phi }_0`$. Additionally, we compute an average flux induced into the interior of the array to be $`\mathrm{\Phi }_{\mathrm{scr}}=0.277\mathrm{\Phi }_0`$ per unit cell. Considering the crudeness of the model, this is remarkable agreement with our measured magnetization (see Fig. 2).

We expect that the screening current will increase monotonically with the external flux, up to $`I_c`$, and therefore expect that the paramagnetic offset in the interior plaquette magnetization will increase as the external flux increases. This is also consistent with our observations in Fig. 2.

The model described here does not account for many features of the real array. We have not accounted for variations in the Josephson-junction parameters and we have ignored the mutual induction between interior plaquettes. Phillips et al. showed the importance of including the mutual inductance matrix when computing array magnetization. Chandran has, in the overdamped junction limit, used a model including the full inductance matrix and computed a paramagnetic moment in a Josephson-junction array. Unfortunately, our arrays are highly underdamped, which makes a model utilizing the full-inductance and dynamics much more computationally difficult.  The simplicity of our model does not describe all of the features that we observe in our array. The observed flux pattern is very complicated, with both paramagnetic and diamagnetic regions.

We have shown that in a multiply connected superconductor such as our array, that the lowest energy configuration that provides diamagnetic screening will also tend to make the interior of the array paramagnetic. Furthermore, the predictions of this model are consistent with our measurements, that the array is more often paramagnetic than diamagnetic. We believe that the paramagnetism in Josephson-junction arrays is best described by the effects detailed above.

We would like to thank J. R. Kirtley, M. Sigrist and M. Chandran for useful communications. Additionally we would like to thank R. L. Greene for his critical reading of our manuscript. We acknowledge support from AFOSR Grant F49620-98-1-0072, NSF Grant DMR-9732800 and NSF Grant DMR-9801825.


text:
# BEPPOSAX OBSERVATIONS OF MARKARIAN 501 IN JUNE 1999

## 1. Scientific goal of the program

The radio-to-$`\gamma `$-ray spectra of blazars ($`\nu F_\nu `$) are typically “double-humped”, with the first peak commonly attributed to synchrotron radiation within a relativistic jet and the second to inverse Compton scattering of relativistic electrons off synchrotron or ambient soft photons (Ulrich, Maraschi & Urry 1997, and references therein). In X-ray bright BL Lacs (XBL) the synchrotron maximum generally occurs in or close to the X-ray band, and the inverse Compton emission peaks above the GeV spectral region extending in some cases to the TeV band, as observed so far for several sources by ground based Cherenkov telescopes. Among these is the nearby ($`z=0.034`$) BL Lac object Mkn 501, one of the brightest blazars at UV, X- and $`\gamma `$-ray energies, and a typical XBL according to the numerous multiwavelength data taken prior to the BeppoSAX launch. Repeated observations with the satellites ASCA and XTE have detected large amplitude X-ray variability at different time scales, often correlated with strong activity in the TeV band (e.g., Kataoka et al. 1999; Sambruna et al. 2000).

BeppoSAX observations of Mkn 501 in April 1997, during an outburst, revealed a completely new behavior. In fact, the joint LECS, MECS and PDS spectra showed that at that epoch the synchrotron component peaked at 100 keV or higher energies (Fig. 1). Correspondingly the source was extremely bright in the TeV band and exhibited rapid flares (Catanese et al. 1997; Aharonian et al. 1999). The flux at 10 keV was an order of magnitude higher than the historical level, while around 1 keV the flux was only moderately brighter than usual. The X-ray spectrum hardens with increasing intensity, and the peak energy of the synchrotron component varies by more than two decades with respect to the quiescent state, a unique behavior in blazars, if compared with the variations of the peak energy exhibited by similar sources, never exceeding an order of magnitude with respect to quiescence (Fossati et al. 2000; Giommi et al. 1999). A discussion of the SED of Mkn 501 in the frame of the SSC model is reported in Tavecchio & Maraschi (this volume).

Further BeppoSAX observations in April-May 1998 showed that the synchrotron peak energy was located at $``$20 keV (Fig. 1), indicating a decrease of an order of magnitude with respect to the previous year (Pian et al. 1999). The simultaneously measured TeV flux was also much lower than in 1997 (Krawczinsky, priv comm.) However, the fact that, despite the radiation losses, the synchrotron peak was still at such high energies one year after the huge outburst, clearly indicates the presence of very powerful, efficient and continuously active mechanisms of particle energization and acceleration in this source.

## 2. BeppoSAX Observations: June 1999 Campaign

We have carried out observations with BeppoSAX in June 1999, simultaneously with TeV Cherenkov telescopes and coordinated with XTE, to investigate the

* X-ray variability of the continuum on short (hours and sub-hour) time scales;
* the presence of possible time lags between soft and hard X-ray emission;
* the correlation between X-ray and gamma-ray flux and spectral index variations.

A single, 180 ks long BeppoSAX pointing has been performed between 1999 June 10, 23:01:18 UT and June 16, 02:11:08 UT, during one of the longest observing windows of the Cherenkov telescopes Whipple, HEGRA, CAT and CELESTE. The X-ray data have been cleaned and linearized at the BeppoSAX Science Data Center (SDC). Spectra and light curves have been extracted from the images with the standard XSELECT package. For the spectral analysis, we used the background files and response matrices distributed by the SDC.

## 3. Results

The source has been clearly detected by the BeppoSAX LECS and MECS instruments, and by the PDS up to 40 keV. The flux exhibits variability of up to 20-30% in amplitude on time scales of 10-12 hours or more (Fig. 1).

The June 1999 flux level at 1 keV is similar to the one observed for Mkn 501 prior to 1997 (“historical” state), namely more than a factor of 2 fainter than detected by BeppoSAX in April-May 1998 and in 1997 April 7, and an order of magnitude fainter than observed by BeppoSAX during the outburst of 1997 April 16 (Fig. 1).

The 0.1-40 keV spectrum is steeper than found in 1997 and 1998 and is progressively steepening with energy. It is not well fitted ($`\chi _r^25`$) by a single power-law plus Galactic absorption ($`N_H=1.73\times 10^{20}`$ cm<sup>-2</sup>), therefore two power-laws have been used to fit the data. The fitted energy break is around 1 keV and the spectral indices are $`\alpha _1=0.89\pm 0.03`$ and $`\alpha _2=1.44\pm 0.02`$ (errors are at 90% confidence level; $`\chi _r^2=1.4`$). Therefore, the X-ray spectrum is consistent with being produced with a unique emission component, which we identify with synchrotron radiation.

In a $`\nu f_\nu `$ representation, the X-ray spectrum appears to peak at the energy of $``$1 keV, which can be identified with the maximum of the synchrotron component. The comparison with the BeppoSAX spectra of the previous years (see Fig. 1) indicates that in 1999 the synchrotron peak has shifted toward lower energies, following the trend already noted in 1997 and 1998, when the peak was observed at $`>`$100 keV and $``$20 keV, respectively.

Optical and TeV coverage simultaneous with the present BeppoSAX campaign was limited due to bad weather. The optical flux level is similar to that usually observed for Mkn 501. At TeV energies, only marginal detections have been obtained on each night. The analysis of these data is still underway.

Our preliminary conclusions are that

* The energy at which the synchrotron component peaks can vary on long term by a large amount (more than a factor of $`200`$), in correspondence with long time scale large amplitude flux variations (2 orders of magnitude at 10 keV).
* This “shift” in energy takes place on much longer time scales (years) than the synchrotron cooling times at X-ray energies, estimated from multiwavelength energy distribution fitting (e.g. Tavecchio & Maraschi, this volume). This further indicates that electrons are continuously accelerated, while the high state is gradually turning to quiescence.
* The low TeV flux in June 1999 indicates that the TeV flux variations are probably correlated with the X-ray variations.

## ACKNOWLEDGEMENTS

We thank the BeppoSAX SDC and Mission Planning Team for their support of this project.

## REFERENCES

Aharonian F.A., Akhperjanian A.G., Barrio J.A., et al., 1999, A&A 342, 69

Catanese, M., et al. 1997; ApJ, 487, L143

Fossati, G., et al. 2000, ApJ in press (astro-ph/0005067)

Giommi, P., et al., 1998, Nucl. Phys. B Proc. Suppl., 69, 407

Kataoka, J., et al. 1999, ApJ, 514, 138

Pian, E., et al. 1998, ApJ, 492, L17

Pian, E., et al. 1999, in BL Lac Phenomenon, ASP Conf. Ser. 159, ed. L. Takalo (San Francisco: ASP), p. 180

Sambruna, R. M., et al. 2000, ApJ in press (astro-ph/0002215)

Ulrich, M.-H., Maraschi, L., & Urry, C. M. 1997, ARAA, 35, 445


text:
# Anomalous Larmour Frequency Dependence of Proton Spin-Lattice Relaxation Time (T1) in the Ferroelectric Glycine Phosphite

## I Introduction

Glycine phosphite (abbreviated as GPI) is one of the many inorganic addition compounds of amino acid group. Most of these addition compounds show interesting physical properties. Some well known ferroelectric addition compounds formed with the amino acid Glycine (+NH3CH2COO-) are tri glycine sulphate (TGS), selenate (TGSe), tetrafluoroberrylate (TGFBe) and diglycine nitrate (DGN) \[1-4\]. Another interesting group of amino acid addition compounds, which is formed by the $`\alpha `$ amino acid betaine ((CH3)3N+CH2COO-) includes among others, betaine phosphate (BP) which is an antiferroelectric, betaine phosphite (BPI) and betaine arsenate (BA) which are ferroelectrics and betaine calciumchloride dihydrate (BCCD) which exhibits a number of phases including incommensurate ones \[5-8\].

GPI (also the other inorganic addition compounds mentioned above) belongs to the family of hydrogen bonded crystals. At room temperature it is monoclinic with space group P21/a . The structure consists of a layered arrangement of planes containing infinite \[H(HPO3)\]$`{}_{}{}^{n}{}_{n}{}^{}`$ chains, alternating with planes built by the organic cations, both layers developing parallel to the bc plane. The chains of HPO3 tetrahedra which are connected by strong hydrogen bonds are directed along the c-axis.

GPI undergoes a ferroelectric phase transition at about 224 K . The spontaneous polarisation appears along the b direction which is perpendicular to the hydrogen-bonded chains of the phosphite anions. This is not so in the case of BPI or caesium dihydrogen phosphate where the hydrogen-bonded chains are parallel to the ferroelectric axis . However, in ferroelectric potassium dihydrogen phosphate (KDP), the hydrogen bonds are perpendicular to the ferroelectric axis.

Dielectric, infrared, Raman, EPR and DSC measurements have indicated that the ferroelectric phase transition in GPI is second order and of order-disorder type . Tritt et.al.  have performed proton NMR spin-lattice relaxation time measurements at 90 MHz to understand the role played by glycinium cations in the ferroelectric transition in GPI . The main feature in the variation of spin-lattice relaxation time (T1) with temperature is the T1 minimum observed at 294 K.

We have taken up the 1 H NMR spin-lattice relaxation time measurements in GPI as a part of our investigation to study proton dynamics in the “Orientational glass” state exhibited by the mixed systems of glycine phosphite and betaine phosphate (abbreviated as BP). Surprisingly, the spin-lattice relaxation time versus 1000/T ( K-1) plot in GPI has exhibited a frequency dispersion in the high temperature side of the T1 minimum. We present here our experimental results and a model to explain the observed T1 dispersion.

## II Experimental Details

Polycrystalline samples of GPI were grown at 308 K from aqueous solution containing glycine and phosphorous acid in the ratio 1:1.by slow evaporation. The samples were purified by recrystallisation. The d-spacings obtained using the powder X-ray diffraction technique agree well with the published crystal structure data . The crystal structure as mentioned earlier, is monoclinic with space group P21/a and the unit cell dimensions are a = 0.9792nm, b = 0.84787nm, c = 0.7411nm, $`\beta `$ = 100.43 Åand Z = 4.

For our experiments, the polycrystalline sample was powdered and sealed in a glass tube after evacuation. 1 H NMR spin- lattice relaxation time were measured at two Larmor frequencies of 11.40MHz and 23.56MHz using a home made pulsed NMR spectrometer working in the range 3-30MHz. A $`\pi \tau \frac{\pi }{2}`$ pulse sequence was employed for spin-lattice relaxation time measurements. The recovery of magnetization followed a single exponential in the temperature range studied . Above 419K the sample melts. The melting temperature observed matched with that obtained from DSC experiments. The temperature of the sample was controlled using a nitrogen gas flow arrangement and measured using a Pt-100 sensor.

## III Results and Discussion

Fig:1 shows a plot of spin-lattice relaxation time (T1) with 1000/T at the two Larmour frequencies of 11.40 MHz and 23.56 MHz. T1 passes through a minimum as the temperature is decreased. The minimum in T1 is at T = 262K (1000/T = 3.82 K-1) and has a value of 10.8 mS at 11.40 MHz. At 23.56 MHz the T1 minimum occurs at 276 K (1000/T = 3.62 K-1) and has a value of 21.7 mS. An interesting feature seen in Fig:1 is the Larmour frequency dependence of T1 at temperatures above the T1 minimum till the highest temperature ( 419 K). The T1 increases with decreasing temperature below the T1 minimum till the lowest temperature studied. The T1 dispersion on the low temperature side of the T1 minimum is as expected from the BPP theory .

In the earlier 1H NMR spin-lattice relaxation work in the same system by Tritt et.al. at a single Larmour frequency of 90 MHz, the spin-lattice relaxation behaviour has been explained by assuming that the reorientation of the NH$`{}_{}{}^{+}{}_{3}{}^{}`$ relaxes all the protons via spin-diffusion in the entire temperature range from 195K to 420 K. However, the model used by Tritt et.al. could not explain our T1 results, especially the T1 dispersion observed on the high temperature side of the T1 minimum. The expression for the proton spin-lattice relaxation rate used by Tritt et.al.  is given as,

$$\frac{1}{T_{1NH_3}}=\frac{3}{8}D\left(\frac{\tau _c}{1+\omega _0^2\tau _c^2}+\frac{4\tau _c}{1+4\omega _0^2\tau _c^2}\right),$$
(1)
where $`\tau _c=\tau _{c0}exp\left(\frac{E_c}{kT}\right)`$ and $`D=\frac{9}{20}\left(\frac{\mu _0}{4\pi }\right)^2\gamma ^4\mathrm{}^2r^6`$. The factor $`\frac{3}{8}`$ in the Eqn:1 arises from spin -diffusion and in this particular case there are three relaxing protons belonging to the NH$`{}_{}{}^{+}{}_{3}{}^{}`$ group, relaxing the total number of eight protons in the molecule. $`\tau _c`$ is the correlation time of the NH$`{}_{}{}^{+}{}_{3}{}^{}`$ group reorientation which is assumed to follow the Arrhenius behaviour, Ec is the activation energy and D is the relaxation constant for the NH$`{}_{}{}^{+}{}_{3}{}^{}`$ group reorientation. The Tritt model predicts a Larmour frequency independent behaviour for T1 on the high temperature side of T1 minimum . The T1 minimum due to reorientation of the NH$`{}_{}{}^{+}{}_{3}{}^{}`$ group alone without spin-diffusion will be 5.5 mS and 11.4 mS at the two frequencies of our study and these are half of the values that have been observed. This can occur if there are inequivalent reorienting groups in the system contributing to relaxation . It is known that if there are Na number of “a” type groups and Nb number of “b” type groups, then the relaxation rate due to these inequivalent groups is given as,

$$\frac{1}{T_1}=\left(N_a\frac{1}{T_{1a}}+N_b\frac{1}{T_{1b}}\right)\frac{1}{N_a+N_b}$$
(2)
We have fitted our data at 11.40 MHz and 23.56 MHz to Eqn:2 and the best fit T1 values do not agree with the observed values for both the cases of Na = Nb and Na$``$Nb. We have also considered the possibility of inequivalence with spin-diffusion but in vein. These failed attempts to fit the data to inequivalent glycinium cations agree with the observations from the vibrational studies in the paraelectric state of GPI .

We have modified the Tritt’s model to explain our observations ( Fig:1) in the following way. We propose that in addition to the reorientation of the NH$`{}_{}{}^{+}{}_{3}{}^{}`$ group, the flip motion of the glycine group (NH2CH2COOH) about the long axis of the GPI molecule may also contribute to the relaxation of the protons particularly at higher temperatures. The relaxation will be dominated by the NH$`{}_{}{}^{+}{}_{3}{}^{}`$ group reorientation at lower temperatures and at higher temperatures where $`\omega \tau _{NH_3}1`$, the 1200 flip of the NH2CH2COOH group could contribute significantly to the relaxation. The 1200 flip in GPI is supported by the second moment experiments carried out by Tritt et.al. .Similar observations where the flip of the entire molecule contributing to the relaxation have been reported earlier in the literature . For example in the case of thiourea CS(NH2)2 the high temperature T1 minimum has been attributed to the flip motion of NH2 . We give below the expression for the proton NMR spin-lattice relaxation rate due to the combined effect of NH$`{}_{}{}^{+}{}_{3}{}^{}`$ reorientation and 1200 flip of NH3CH2COOH group along with spin-diffusion.

$`{\displaystyle \frac{1}{T_1}}`$ $`=`$ $`{\displaystyle \frac{3}{8}}D\left({\displaystyle \frac{\tau _c}{1+\omega _0^2\tau _c^2}}+{\displaystyle \frac{4\tau _c}{1+4\omega _0^2\tau _c^2}}\right)`$ (4)
$`+{\displaystyle \frac{2}{8}}F\left({\displaystyle \frac{\tau _i}{1+\omega _0^2\tau _i^2}}+{\displaystyle \frac{4\tau _i}{1+4\omega _0^2\tau _i^2}}\right)`$
where $`\tau _c`$ and $`\tau _i`$ are the correlation times for the NH$`{}_{}{}^{+}{}_{3}{}^{}`$ group reorientation and NH3CH2COOH group flip motion respectively. Both the correlation times are assumed to follow the Arrhenius behaviour. D and F represent the relaxation constants for the NH$`{}_{}{}^{+}{}_{3}{}^{}`$ group reorientation and NH3CH2COOH group flip motion respectively. The factors $`\frac{3}{8}`$ and $`\frac{2}{8}`$ come from spin-diffusion. The constant D is given as in Eqn:1 and F is given by ,

$$F=\frac{27}{40}\left(\frac{\mu _0}{4\pi }\right)^2\frac{\gamma ^4\mathrm{}^2}{R^6}\mathrm{sin}^2\alpha $$
(5)
where $`\alpha `$ is the flip angle , R is the inter -proton distance between the protons attached to the CH2 group and the NH3 group.

We have done a least square fitting of the T1 vs 1000/T data to Eqn: 4 and the best fit is shown in Fig:1 and the best fit parameters obtained at 11.40 MHz and 23.56 MHz are given in the table I. The correlation times obtained for the reorientation and flip motion compare well with those reported in the literature . The activation energy of 26.8 kJ/mol for the NH$`{}_{}{}^{+}{}_{3}{}^{}`$ group reorientation (which agrees with that of Tritt et.al.) is relatively high compared to that in the other glycine compounds like TGS, TGSe etc . Tritt et.al.  have attributed this high value to the relatively strong hydrogen bonds in GPI compared to the other glycine systems. The high value of 45.6 kJ/mol for the activation energy of the flip of NH3CH2COOH group is also to be expected and therefore this motion comes into play only at higher temperatures. Similar values of activation energy for the flip motion have been reported in other systems like thiourea, hydraziniumnitrate etc . The value of inter-proton distance in the NH$`{}_{}{}^{+}{}_{3}{}^{}`$ group calculated using the relaxation constant D is 1.69 $`\AA `$ and it matches well with the value of 1.70 $`\AA `$ obtained from neutron diffraction experiments in amino acids . We have not observed any distinct T1 minimum due to the flip motion. But a careful look at the spin-lattice relaxation behaviour (Fig:1) in the higher temperature region at 11.40 MHz , shows a weak indication around 370K. As the flip motion is a high activation energy motion, the motional frequencies may not be very high and it is expected to show up significantly only at lower Larmour frequencies. Fig:2 shows the simulated T1 vs 1000/T curve using the motional parameters from the table at 11.40 MHz for reorientational motion, flip motion and the combined motion of reorientation and flip. The arrow indicates the approximate temperature at which T1 minimum is expected due to the flip motion. As the Larmour frequency increases the T1 minimum from the flip motion shifts to higher temperatures. The data of Tritt et.al. was obtained at a higher Larmour frequency of 90 MHz where the flip motion of the glycine group as a whole is not very effective and hence could be explained on the basis of NH$`{}_{}{}^{+}{}_{3}{}^{}`$ reorientation alone in the range of temperature studied. $`\omega _0\tau _i`$ at 11.40 MHz and 23.56 MHz do not satisfy the condition $`\omega _0\tau _i1`$ at higher temperatures and this explains the T1 dispersion in the data at higher temperatures.

It would be interesting to compare the behaviour of T1 near Tc in GPI (Tc = 224 K) and other glycine compounds which exhibit ferroelectric transition. In triglycine sulphate (TGS, Tc = 320 K),selenate (TGSe, Tc = 295 K), tetrafluoroberrylate (TGFBe, Tc = 343 K), partially deuterated TGS (Tc = 333 K) and potassium dihydrogen phosphate (KDP, Tc = 123 K) the transition from the paraelectric phase to the ferroelectric phase has been found to be accompanied by a rather sudden decrease in the proton spin-lattice relaxation time . KDP is compared here due to its similarity with GPI in the direction of spontaneous polarisation. In the TG systems a sudden change in T1 near the transition temeperature indicates a direct involvement of protons in the ferroelectric transition. In GPI also a large isotopic effect (98 K) on Tc, points towards the essential role of the interphosphite hydrogen bonds in the phase transition mechanism. However the T1 varies smoothly with temperature near the ferroelectric transition temperature of 224 K  indicating that the dynamics of the protons in the hygrogen bonds are fast for the NMR T1 scale. This fast proton jump between the interphosphite hydrogen bond also could be the reason for a lower ferroelectric transition temperature in GPI compared to the other ferroelectrics mentioned above . Vibrational spectroscopic studies show that normal H2PO$`{}_{}{}^{1}{}_{3}{}^{}`$ anions appear in the paraelectric phase and the ferroelectric phase, confirming that the protons in the interphosphite hydrogen bonds are dynamically disorderd in the time scale of vibrational spectroscopy . These studies also suggest that the dynamical disorder of interphosphite hydrogen bonds is coupled to the motions of glycinium cations . The above statement is well supported by the splitting of some glycinium bands in the ferroelectric phase where the interphosphite hydrogen bonds settle to one of the positions in the time scale of vibrational spectroscopy. As stated above the motion of these protons appears to be too fast ($`\omega _0\tau _c1`$) to be observed by the NMR T1 measurements.

## IV Conclusion

We have observed a Larmour frequency dependence in the spin-lattice relaxation time behaviour in GPI on the high temperature side of the T1 minimum . Presence of inequivalent glycinium cations could not explain our observed results and an attempt has been made to explain this by invoking the 1200 flip motion of NH2CH2COOH group about it’s long axis, in addition to the NH$`{}_{}{}^{+}{}_{3}{}^{}`$ reorientation . The motional parameters obtained match well with those obtained in the earlier studies. No significant change in T1 behaviour could be observed at Tc.

We are grateful to Prof. H. L. Bhat and Ms.Deepthy. A. Pillai for providing us the samples of Glycine phosphite used in these investigations.