text:
# VLBI Polarimetry of 177 Sources from the Caltech-Jodrell Bank Flat-spectrum Survey

## 1 Introduction

Studies of the polarization properties of jets on the milliarcsecond scale can yield vital insights about their formation, collimation, and propagation \[e.g. Gómez et al., 2000; Meier, 2001\]. Based on results from 25 BL Lacs, Gabuzda, Pushkarev & Cawthorne  claim a tendency for the electric vector polarization angle to lie along the jet axis for these objects, with about 30 percent showing a perpendicular orientation. This has been interpreted in terms of oblique and transverse shocks \[Gabuzda, Pushkarev & Cawthorne, 2000\] or as the result of a helical magnetic field confining the jet \[Gabuzda, 2002\].

Although many individual sources and even large surveys have been undertaken with VLBI in total intensity, the number of sources studied with VLBI polarimetry is far less. In the early years of VLBI this was due in part to hardware restrictions (e.g., feeds with only a single available polarization), but also due to software restrictions that made calibration difficult. The VLBI group at Brandeis pioneered efforts to develop techniques and software to make VLBI polarimetry more tractable \[Roberts et al., 1984; Roberts, Brown, & Wardle, 1991\], and this culminated in a study of the mas polarization structure of 24 sources from the Pearson-Readhead sample \[Cawthorne et al., 1993\]. But it wasn’t until the advent of the VLBA, with homogeneous, low-leakage antennas and feeds, that VLBI polarimetry became widely accessible. Full polarization calibration techniques were implemented in AIPS in the 1990s and in 1999 NRAO began providing regular monitoring111http://www.aoc.nrao.edu/$``$smyers/calibration/ of a suite of calibrators that can be used for the absolute polarization angle calibration \[Taylor & Myers, 2000\].

The Caltech-Jodrell Bank flat-spectrum (CJF) survey is a large complete, flux-limited sample of sources imaged with VLBI, consisting of 293 sources \[Taylor et al., 1996\]. This survey has been used to place constraints on gravitational lensing by 105 M black holes \[Wilkinson et al., 2001\], to study the unbeamed synchrotron luminosities of relativistic jets \[Lister & Marscher, 1997\], and recently to study the relation between the linear size of the VLBI jet and the black hole mass \[Cao & Jiang, 2002\].

As part of a proper-motion study of compact objects, all the sources in the CJF sample were observed at least 3 times at intervals of approximately 2 years \[Britzen et al., 1998\]. The early epochs used Global VLBI observations, but starting in 1995 the VLBA was used exclusively since it could provide better ($`u,v`$) coverage and consequently improved image fidelity. The first observing sessions with the VLBA were limited in bandwidth (8 MHz) and used only a single polarization, but as the VLBA matured and its capabilities increased, we began in 1998 to observe with 16 MHz bandwidth in full polarization mode. The last three epochs of CJF were observed in this mode and provide full VLBI polarimetry of 177 sources.

Here we present a statistical analysis of the properties of those 177 sources in CJF observed with full polarimetry. The sample is described in detail in §2, and the calibration procedures are described in §3. We present results in §4, and discuss comparisons of source properties, and evolution of those properties with redshift in §5.

We assume H$`{}_{0}{}^{}=75`$ km s-1 Mpc-1 and $`q_0`$=0.5 throughout.

## 2 Sample Selection

We present 5 GHz Very Long Baseline Interferometric (VLBI) polarimetry observations of 177 sources in the Caltech-Jodrell Bank flat-spectrum (CJF) survey, a complete flux-density-limited sample of 293 flat-spectrum radio sources. The CJF sample requires that the flux density at 4850 MHz be at least 350 mJy, ($`S_{4850}350`$ mJy), and that the spectral index at 4850 MHz and 1400 MHz be at least $`0.5`$ ($`\alpha _{1400}^{4850}0.5`$). Additionally, the sample is limited by B1950 declinations and galactic latitudes such that $`\delta 35^{}`$ and $`|b|`$ $`10^{}`$. A list of the 293 extragalactic radio sources and a more complete description of the selection criteria for the CJF sample can be found in Taylor et al. . The CJF sample has been completely observed with VLBI as part of the PR \[Pearson & Readhead, 1988\], CJ1 \[Polatidis et al., 1995\] and CJ2 \[Taylor et al., 1994\] surveys.

The 177 sources presented here were observed as part of a proper motion study \[Britzen et al., 1998\] so that polarimetric data in multiple epochs are available for 66 of the 177 sources. Where data from multiple epochs was available, we preferred the 2000.958 observations over the 1998.122 observations, and used the 1999.890 observations, with somewhat poorer u,v coverage, only when other epochs were unavailable. All plots shown here include just one epoch of each of the 177 sources. In the future it may be interesting to study variability in polarized intensity and polarization angle using those 66 sources with multiple epochs available, however no effort has been made toward this goal thus far.

Using optical classifications from the literature \[Henstock et al., 1995; Vermeulen, Taylor, Readhead, & Browne, 1996\], from the NASA Extragalactic Database, and from unpublished observations by the CJ collaboration, we find that the final sample of 177 sources consists of 30 galaxies, 106 quasars, 20 BL Lacertae objects and 21 others, including empty fields, red and blue objects and other optically unresolved objects. We note here a few words of caution regarding the above classifications in that low-luminosity BL Lacs can be missidentified as quasars or radio galaxies \[Marchã & Browne, 1995\], and variability in the emission lines of BL Lacs and quasars can cause some movement between these classes \[Vermeulen et al., 1995\]. Redshifts have been measured for 147 of the sources presented here.

## 3 Observations

### 3.1 Calibration

The observations were carried out on 2000 December 16 and 17 (2000.958), 1999 November 21, 23 and 26 (1999.890), and 1998 February 8, 12, 13, 20 and 21 (1998.122). These three epochs provided about 34, 72 and 120 hours of data, respectively, for a total of 226 hours of observations. The 2000.958 and 1998.122 observations used the 10 element Very Long Baseline Array222The National Radio Astronomy Observatory is operated by Associated Universities, Inc., under cooperative agreement with the National Science Foundation. (VLBA) while the 1999.890 observations used only 8 of the VLBA antennas. The Saint Croix, Virgin Islands antenna was lost due to hurricane Lenny, and the North Liberty, Iowa antenna was lost due to a power hardware problem. Because we did not use phase referencing, we failed to image 6 sources with total intensity peaks of less than 55 mJy. Of these, four sources were identified with galaxies, and two sources had no optical identification. Another source, 0954+556, was omitted due to a position error of 30 arcseconds. Finally, the gravitationally lensed object 0218+357 was also omitted. These 8 sources that we did not image have not been included in any of the distributions or analysis, and are not included in the 177 source count.

We observed at 4995 MHz with a total bandwidth of 16 MHz. Right- and left-circular polarizations were recorded using 1 bit sampling at 64 Mbits/sec. Two intermediate frequencies (IFs) of 8 MHz each were used for right- and left- polarization. Amplitude calibration for each antenna was derived from measurements of the antenna gain and system temperatures during each run. Global fringe fitting was performed using the AIPS task FRING, an implementation of the Schwab & Cotton  algorithm. The fringe fitting was performed on each IF and polarization independently using a solution interval of 2 minutes, and a point source model was assumed. Next, a short segment of the cross hand data from the strongly polarized calibrator 3C 279 was fringe fitted in order to determine the right-left delay difference, and the correction obtained was applied to the rest of the data. Once delay and rate solutions were applied the data were averaged in frequency over 8 MHz. All data were then averaged together over a 30 second time interval. Imaging, editing and self-calibration were performed in Difmap \[Shepherd, Pearson & Taylor, 1995\] and AIPS.

Determination of the leakage terms was performed with the AIPS task LPCAL using a strong calibrator with simple source structure such as 0923+392 or 0716+714. Given the nature of the survey, strong sources with good parallactic angle coverage were abundant in every observing session. The absolute electric vector polarization angle (EVPA) calibration was performed using various well-known target sources such as 0923+392, 1803+784, and 2200+420 (BL Lac). These observations were compared to contemporaneous short VLA observations and VLBI polarimetry reported in the literature, or available from the VLA/VLBA Polarization Calibration Page \[Taylor & Myers, 2000\]. We also used component C4 of 3C 279, which was fairly stable during this time period \[Taylor, 2000; Zavala & Taylor, 2001\]. The absolute uncertainty in the EVPA calibration is $``$4, with the dominant source of error being the variability of the EVPAs of the calibrator sources \[Taylor & Myers, 2000\].

### 3.2 Data Analysis

We automated the imaging process by writing a script that ran the AIPS task IMAGR on the self-calibrated data for every source, creating Stokes I, Q and U maps. Using another script, we then ran the AIPS task COMB on each of the sources to produce polarized intensity maps, polarized intensity noise maps, polarization angle maps and polarization angle noise maps. Some sample images covering a range of morphologies are shown in Fig. 1. These FITS files were brought into IDL where we used programs to objectively and uniformly determine the properties of each source.

### 3.3 Source Property Definitions

The RMS noise levels were calculated for the total intensity and linear polarizations by taking the standard deviations of four corner sections of the image and averaging the two middle values. We calculated the integrated flux in Stokes I, Q and U by summing the flux in a particular region that completely enclosed each source. The algorithm used to systematically define the region boundaries creates a rectangular box around those pixels with values at least 15 times the RMS noise level. Due to image artifacts this high cut-off proved necessary to correctly enclose many sources. A border of width three times the FWHM of the convolving beam was added to the initial rectangular area to define the final region to be used in the summing. As a check we compared the integrated flux calculated in IDL to the sum of the clean components found in the AIPS task IMAGR. The total integrated flux measurements of all 177 sources were found to agree with that calculated by the AIPS clean algorithm to within 8%.

The core of each source was assumed to be located at the brightest peak in the total intensity map when no other data were available for study. However, in the case of compact symmetric objects (CSOs) or sources with complex jet structures we took advantage of 15 GHz images \[Taylor, Readhead & Pearson, 1996\] (and unpublished 15 GHz observations) and published motion studies \[Taylor & Vermeulen, 1997\] to more accurately determine core positions. Jet component positions were found by subtracting from the initial total intensity image a Gaussian model centered at the peak flux position and scaled so that the maximum value of the Gaussian equaled the peak value in the map. As the appropriate Gaussian is subtracted from the previous image, the locations of each new peak define the locations of consecutive jet components. We add the restriction that the total intensity at the jet component position found must be greater than 9 times the RMS noise level, where the noise is calculated as discussed above. At most 8 jet component positions were required to describe a source. The position angle of the elliptical Gaussian model used in the subtraction was matched to the position angle of the convolving beam, and both axes of the Gaussian were given a width of 1.5$`\sigma `$, where $`\sigma =\frac{FWHM}{2\sqrt{2\mathrm{ln}2}}`$ and describes the width of the convolving beam. Due to slight differences in the editing process that changed the u,v coverage for different Stokes parameters, the beam sizes and position angles fit by the AIPS task IMAGR differed slightly in I, Q and U. For consistency we forced all beams to have the dimensions and orientation of the Stokes U beam by setting the beam parameters BMAJ, BMIN and BPA in AIPS for the Stokes Q and I images.

The jet properties and source morphologies based on total intensities were easily characterized after each jet component position was defined. We found 34 sources with no jet components which we call naked cores and 134 sources with at least one jet component which we call core-jet sources. In addition, we found 9 sources which we classified as Compact Symmetric Objects (CSOs) or CSO candidates based on their symmetric structures. CSOs are a recently identified class of radio sources smaller than 1 kpc in size with emission on both sides of the central engine, and are thought to be very young objects ($``$1000 y, Readhead et al. ; Owsianik & Conway ). Eight of these sources were previously identified as CSOs or CSO candidates \[Peck & Taylor, 2000\], and the remaining source (0402+379) is currently under study.

We recorded the jet axis angle to be the angle defined by the core and the closest jet component, where North is $`0^{}`$. The only exception to this is the CSO candidate 0402+379 which has a component close to the core but seemingly unrelated to the jet. The jet axis angle for this source was measured by eye. The jet length is defined in most cases as the distance from the core to the farthest jet component from the core, irrespective of jet bend. It is quoted in milliarcseconds and is uncorrected for angular distances. In §5 we model our angular size calculation after the description given in Kellermann , creating a contour line at 2% of the peak in the total intensity map. The maximum distance from the contour line to the core position defines the angular size. To avoid measuring distances to spurious peaks above the 2% level, we confined the contour to be within the same rectangular area that was previously determined to completely enclose the source.

We define integrated core flux densities as the pixel value at the core component, where the pixel values have units of Jy beam-1 and we assume the core itself to be unresolved. Due to complex jet structure, rather than sum the pixel values for the jet components we define the integrated jet flux as the difference between the integrated total source flux and the core flux.

Due to leakage between orthogonal polarization feeds in the antennas, we define detected, polarized cores or jet components as only those components in which the polarized intensity is greater than 0.2% of the peak in the total intensity image (p $`>`$ 0.002 Ipeak). We require detected polarizations to be above the 3$`\sigma `$ level, where $`\sigma `$ is described by the RMS noise in the Stokes Q and U maps by the relation $`\sigma =(0.21\sigma _U^2+0.21\sigma _Q^2)^{1/2}`$ according to Rayleigh statistics. In the case that the core has undetected polarization, we quote mcore as an upper limit equal to the larger of the two detection restrictions (m$`{}_{\mathrm{core}}{}^{}`$ max(0.002 Ipeak, 3$`\sigma `$)). When studying correlations between jet axis angle and polarization angle, we restrict our analysis to those components with $`\sigma _\chi <20^{}`$, where $`\sigma _\chi `$ is determined by the noise image created by the AIPS task COMB.

## 4 Results

Measured properties for each source are presented in Table 1. Additional information, including measured jet properties, is available on the world wide web333www.aoc.nrao.edu/$``$gtaylor/cjftab.text. Images are available upon request to the authors. Below we present a statistical analysis of source properties.

### 4.1 Core Fraction

From Fig. 2 it is clear that the core fractions of galaxies come from a different parent distribution than that of quasars or BL Lacertae objects. While quasars and BL Lacs generally exhibit a core fraction (defined $`R_\mathrm{c}=S_\mathrm{c}/S_{\mathrm{total}}`$) $`R_\mathrm{c}`$80%, the 30 galaxies in our sample show no tendency toward either high or low core fraction. The tendency towards higher core fraction in quasar and BL Lacs is expected from unified schemes since the core fraction is a strong function of orientation. The Kolmogorov-Smirnov (K-S) test disproves the null hypothesis that the distribution of quasars ($`d_\mathrm{Q}`$) and the distribution of BL Lacs ($`d_{\mathrm{BL}}`$) are the same as the distribution of galaxy core fractions ($`d_\mathrm{G}`$). When comparing $`d_\mathrm{Q}`$ to $`d_\mathrm{G}`$ the K-S test outputs the maximum value of the absolute difference between the two cumulative distribution functions (c.d.f.’s), $`D`$=0.45 with a probability, $`p`$, that quasars and galaxies are drawn from the same population of $`7.7\times 10^5`$. Similarly, when comparing $`d_{\mathrm{BL}}`$ to $`d_\mathrm{G}`$, we find a probability of $`2.9\times 10^3`$. Comparing $`d_{\mathrm{BL}}`$ to $`d_\mathrm{Q}`$ we find $`p`$=0.41, indicating that these samples could come from the same parent distribution. Since the cores and jets generally have different spectral indices, the core fraction will depend weekly on redshift. As an example, a source with $`R_\mathrm{c}`$ = 0.75 at a redshift of zero would shift to $`R_c`$ = 0.83 at z=1 and $`R_\mathrm{c}`$ = 0.87 at z=2, assuming a flat spectrum core and steep ($`\alpha =0.7`$) jets. In Fig. 2 we have not attempted any redshift correction so this could account for the minor difference between the shape of the distribution for the quasars and BL Lac objects.

The distribution of $`R_\mathrm{c}`$ in Fig. 2 for the galaxies appears somewhat bimodal. Of the 15 sources with $`R_\mathrm{c}<0.5`$ more than half (9/15) are CSOs. No obvious properties are identified with the high $`R_\mathrm{c}`$ population, although our knowledge of these sources is not complete.

One might be concerned that this comparison uses a sample of quasars four to five times as large as the samples of BL Lacs and radio galaxies, respectively. This should not be a factor if we consider the effective number of data points $`N_e`$ as defined by Press et al. . $`N_e`$ is defined as

$$N_e=\frac{N_1N_2}{N_1+N_2}$$
(1)
where $`N_i`$ are the number of data points in the respective sets. As long as $`N_e`$ is greater than or equal to four the probability estimate should be accurate. For the quasar$``$galaxy comparison $`N_e=23.4`$, quasar$``$BL Lac $`N_e=16.8`$, and for BL Lac$``$galaxy $`N_e=12.0`$.

### 4.2 Core and Jet Fractional Polarizations

Fig. 3 shows the distribution of average fractional polarization in the cores of quasars, galaxies and BL Lac objects. From this distribution we see that quasars have typical core polarizations of a few percent, and most (87%) are detected. Galaxies are clearly less polarized than quasars on the parsec scale with only 3 of 30 (10%) of the sources found to have detectable polarization in their cores. No detectable polarized flux was found from any of the 9 CSOs observed. This is consistent with results from a VLBI follow-up survey of CSOs by Peck & Taylor  who found less than 1.2 mJy of polarized flux from each of 21 sources observed at 8.4 GHz with the VLBA. Generally we find less than 5 mJy of polarized flux density and fractional polarizations of $`<`$0.5%. For a few CSOs with weak cores, our upper limits on the fractional polarization are quite high. In contrast the BL Lacs appear to have slightly more strongly polarized cores than quasars (K-S test gives $`p`$=0.20 that they are drawn from the same population).

The K-S test only compares detections, and ignores the upper limits. Survival analysis \[Feigelson & Nelson, 1985\] takes into account the upper limits in the comparison of two samples. Using the Penn State ASURV 1.2 package \[Lavalley, Isobe, & Feigelson, 1992\] we tested the hypothesis that the BL Lac and quasar core fractional polarization results were drawn from the same parent distribution. Including the 15 upper limit quasar results decreased the likelihood that the two samples were drawn from the same parent distribution. The five two sample tests in ASURV 1.2 give an average probability of 1.5% that the null hypothesis is correct. This is significantly different from the K-S test result, and we find that this result is not affected by eliminating the one rather high upper limit quasar core fractional polarization of 2.7%.

Fig. 4 shows the same set of distributions, but for the average fractional polarizations of the jet components. In this case the fractional polarization of the BL Lacs is significantly greater from that of the quasars (p=0.0037 that they are drawn from the same population). When the jet polarization is detected, sources have higher fractional polarization in the jets than in the cores. Even the galaxies have an average fractional polarization in the jets of 10%, compared to typical fractional core polarizations of $`<`$1%. In Fig. 5 we show the images for all six radio galaxies with some detected polarization. The measurements range from marginal in 0600+442 and 0847+379 to extremely detailed in the bright, one-sided jet of 1807+698.

### 4.3 Comparing Alignments Between Jet Axis and EVPA

Many claims have been made in the literature regarding comparisons of the orientation of the jet axis ($`\theta `$) with the EVPA of the jet ($`\chi _{\mathrm{jet}}`$) or core ($`\chi _{\mathrm{core}}`$). Below we discuss various tests to look for preferred alignments. Ideally these tests should be made after having removed any Faraday rotation due to the ISM of our Galaxy, or the environment local to the sources. Unfortunately, for these single-frequency observations it was not possible to estimate the Rotation Measure (RM). If significant RMs exist towards any of the source components (e.g., as found to be common in quasar cores by Taylor 1998, 2000) then that will tend to smear out any intrinsic correlation.

The rotation measure of the cores of BL Lacs are typically 2$``$300 rad m-2, and those of quasars range from 200 rad m-2 to several 1000 rad m-2 \[Zavala & Taylor, 2003\]. At 5 GHz a rotation measure of 277 rad m-2 is enough to cause a turn of one radian. This would be sufficient to smear out any correlation in EVPA. The RMs in the jets of 16 BL Lacs, radio galaxies and quasars, with one exception, are all approximately 200 rad m-2, and thus would also smear out any correlation in $`\theta \chi _{\mathrm{jet}}`$ \[Zavala & Taylor, 2003\].

In Fig. 6 we present the distribution of $`\theta \chi _{\mathrm{core}}`$ for quasars, galaxies and BL Lac objects. For the quasars we see that the distribution is not quite flat, and hints at an excess of sources near 80 degrees. These nearly misaligned EVPAs imply a magnetic field nearly perpendicular to the jet axis assuming optically thick synchrotron emission. The K-S test indicates only a 0.3% chance that the distribution is the same as a flat distribution, suggesting the possibility of an intrinsic correlation between EVPA and $`\theta `$. There is a somewhat greater probability (3.7%) that the quasars and BL Lacs are drawn from the same population. However, the BL Lac distribution itself is not significantly different from a flat distribution (17% chance that it is the same as flat).

To further investigate the nature of the misaligned core EVPAs we divided the quasars into those with short (length $`<`$ 6 mas) jets and long (length $`>`$ 6 mas) jets. From Fig. 7 we see that a jet length of 6 mas divides the quasars into two groups of roughly equal sizes since there are jet length peaks at 4 and 9 mas. We used ASURV’s correlation tests \[Isobe, Feigelson, & Nelson, 1986\] to examine the correlation that appears in Fig. 8 for the quasars with short jets. Although it appears that short jets are correlated with misaligned core EVPAs, a Kendall’s tau test gives a less than 1 sigma significance of a correlation. Additionally, plots of $`|\theta \chi _{core}|`$ versus jet length show no obvious correlation, regardless of jet length. Although there is no evidence for a correlation of quasar core EVPA alignment with jet length the K$``$S test of quasars versus a uniform distribution indicates that quasars have a preference for misaligned core EVPAs. As one can see from Fig. 8 the cores associated with short jets appear more misaligned (probability of a flat distribution is 0.009%, probability that they have the same distribution as sources with long jets is 1.7%), with a clear excess near 90 degrees. We also divided the quasars into those with a high ($`R_c>85`$%) or low ($`R_c<85`$%) core fraction. Core fraction seems much less important as a discriminator of the misaligned core EVPA population.

One can also look for any dependence on alignment with redshift. In Fig. 9 we plot the alignments of the jet and core EVPAs with jet axis against redshift. No trends with redshift are apparent for either the jets or the cores. Correlation tests performed using the ASURV package verified that no correlation of either core or jet misalignment with redshift is present. This was true for fractional polarization in the cores and jets as well.

In Fig. 10 we present the distribution of $`\theta \chi _{\mathrm{jet}}`$ for quasars, galaxies and BL Lac objects for the jet components both near ($`<`$ 6 mas) and far ($`>`$ 6 mas) from the core. Contrary to claims in the literature \[Cawthorne et al., 1993; Homan et al., 2002\] that quasars have preferred alignments with the magnetic field oriented parallel to the jet axis, we see no such correlation. (The probability of being drawn from a flat distribution is 34%.) We have also tried examining separately the near ($`<`$ 6 mas) and far ($`>`$ 6 mas) jet components and reach the same conclusions. Cawthorne et al. did apply a rotation measure correction, which might reveal a correlation our non-RM corrected results hide. However, their RMs were derived from the VLA, which does not correctly sample the parsec-scale RM of these AGN. Additionally, the RM corrections were primarily (12/17) less than 50 rad m-2, which amounts to a correction of 10 or less at 5 GHz. As shown in Zavala & Taylor  these corrections were too small to account for the known parsec-scale RMs. Homan et al. did not apply a RM correction, but their observations were at 15 and 22 GHz. Their results are thus less susceptible to any rotation measure smearing.

The distribution for BL Lac objects (Fig. 10) appears to have a lack of sources with EVPA perpendicular to the jet. From the K-S test we find an 8.0% chance of this population being drawn from a uniform $`\theta \chi _{\mathrm{jet}}`$ distribution. A bimodal distribution was found by Gabuzda, Pushkarev & Cawthorne  with peaks near 0 and 90. Gabuzda et al. performed a similar analysis on a 1 Jy flux limited sample of BL Lac objects observed with VLBI polarimetry at 5 GHz. This is shown in Fig. 12 of Gabuzda, Pushkarev & Cawthorne . A K-S test of the data in their Fig. 12 yields a 1.1% probability that the data are drawn from a uniform distribution. With assistance from Gabuzda (2002, private communication) we combined our BL Lac $`\theta \chi _{\mathrm{jet}}`$ data with the 35 measurements from Gabuzda, Pushkarev & Cawthorne  and tested this data against a uniform distribution in a K-S test. The combination of these two samples results in a 0.7% probability that the parent distribution is uniform. However, with just the 17 data points from our sample it is difficult to draw a conclusion on the nature of the underlying distribution, and we cannot hope to isolate what type of source gives rise to the slight excess of sources with EVPA at small angles to the jet axis.

In a 43 GHz polarimetry survey of flat-spectrum sources from the Pearson-Readhead survey \[Pearson & Readhead, 1988\], Lister  found that the most strongly polarized quasar cores display EVPAs that are aligned with the jet axis. In Fig. 11 we examine the relationship between alignment and fractional polarization. No preferred alignment for the more polarized sources, as what Lister  found, is evident in our sample. Also the preferred core EVPA orientation in our sample is for EVPAs perpendicular to the jet (Fig. 6). These two claims can be reconciled with the same projected magnetic field orientation if the cores are optically thin at 43 GHz and optically thick at 5 GHz.

Using the ASURV bivariate tests \[Isobe, Feigelson, & Nelson, 1986\] we examined the hypothesis that there is a correlation between the core fractional polarization and the core EVPA alignment. The Kendall’s tau test gives a probability of 82.7% and the Cox regression a probability of 96.9% that a correlation does not exist. Although there appears to be a slight excess of sources with misaligned core EVPAs at fractional polarizations greater than 2.5%, this is not born out by the bivariate tests. The Cox regression gives a probability of 94.7% and the Kendall’s tau test a 43.6% probability that no correlation is present. The disparity probably results from small number statistics as only 20 sources have m$`{}_{\mathrm{core}}{}^{}>`$ 2.5%.

### 4.4 Fractional Polarization and Jet Length

In Fig. 12 the core and average jet fractional polarization is plotted against the projected jet length. The fractional jet polarization seems fairly smoothly distributed, but the fractional core polarization is anti-correlated with jet length. Sources with short jets seem to have higher core polarizations and there is a distinct lack of strongly polarized cores in sources with long jets. The censoring present in the lower panel of Fig. 12 led us to use survival analysis to test for a correlation for core EVPA alignment with jet length. Both the Kendall’s tau and Spearman’s rho tests give a weak (1.6 $`\sigma `$) support to the anti-correlation which seems to be present. A possible explanation for this trend is that sources with long jets are viewed at angles farther from the line-of-sight. Assuming the standard obscuring torus model, the angle to the core then traverses a higher-density region with magnetic fields causing Faraday depolarization. This interpretation is also consistent with the finding that the cores of BL Lacs (thought to be viewed at very small angles to the line-of-sight \[Antonucci, 1993\]) are somewhat more strongly polarized than quasar cores.

## 5 Cosmology Using the Size$``$Redshift Relation

The $`\theta z`$ relation for compact radio sources has been studied in the past to place limits on $`q_{}`$. Kellermann  used a sample of 82 core-jet sources selected from the literature to find that the deceleration parameter is roughly 0.5, without appealing to source evolution. Gurvits, Kellermann, & Frey  did a similar analysis using 330 5 GHz compact sources taken from the literature, finding that $`q_{}=0.21\pm 0.30`$ if no linear size$``$luminosity, $``$redshift or $``$spectral index dependences are assumed. Gurvits, Kellermann, & Frey  applied a regression model to 145 of the 330 sources with additional restrictions on luminosity and spectral index. With this subset of sources a deceleration parameter of $`q_{}=0.33\pm 0.11`$ was found. All of these results are consistent with Friedmann cosmologies with $`0q_{}0.5`$ and $`\mathrm{\Lambda }`$=0.

Figure 13 shows the $`\theta `$z relationship for 103 core-jet sources from our complete sample and mirrors Fig. 1 in Kellermann  and Fig. 5 in Gurvits, Kellermann, & Frey . We measured angular sizes similarly to how Kellermann and Gurvits did, accounting for the decrease in brightness with redshift as described in § 3.3. CSOs have been excluded due to the proposed relation between their age and angular size, and BL Lacertae objects have been excluded because of their proposed preferred inclination angles and to be consistent with the methods of Kellermann . The redshifts range from 0.0172 $`z`$ 3.469. With the exception of the slope between the first two redshift bins, our plot shows no clear dependence between angular size and redshift. Binning the data in redshift space rather than restricting each bin to equal numbers of sources does not greatly alter the plot. We found that plotting median values rather than mean values also does not change the plot’s shape considerably. Thus contrary to previous claims our complete, flux-density limited sample, although not inconsistent with Friedmann cosmologies, cannot definitively rule out a steady state universe or place any limits on $`q_{}`$. The few sources with higher angular sizes and reasonable error bars clustered around z=1.1 are, however, suggestive of a general increase of angular size with redshift after z=1, which supports inflationary cosmologies.

Although these compact radio sources with smaller redshifts and younger ages are not as subject to evolutionary effects, and although the physical parameters of their central engines are thought to outweigh any effects due to possible variations of intergalactic medium with redshift, we find that these objects nevertheless cannot be used to place restrictions on the deceleration parameter and the geometry of our universe. The large error bars in Fig. 13 prove these sources to be poor standard rods, possibly because of their rapid evolution compared to the Hubble time. It is possible that strong Doppler favoritism is skewing our results so that more objects with small inclinations to the line-of-sight and therefore greater, Doppler-boosted luminosities and smaller angular sizes are getting into the sample without being excluded by the flux density requirement. Our neglect of Doppler boosting may also affect our measurement of angular size since we have defined object size by the contour line at 2% of the peak in the total intensity map. Angular size might also be better measured taking into account jet bend, rather than assuming a one-dimensional jet in all cases. Lastly, it is possible that with a multi-frequency study using an angular size$``$spectral index relation, as is shown in Fig. 7 of Gurvits, Kellermann, & Frey , the $`\theta z`$ correlation may become tighter.

### 5.1 Bent-Jet Sources and Other Interesting Properties

Although we did not do an exhaustive statistical analysis of the correlation between jet bending and other source properties, we do note two basic categories into which most bent sources fall. One class of bent jets typically shows little or no polarized intensity in the jet, while there may be polarized intensity at the source core. Some examples of this are 0600+442, 0700+470, 0707+476, 0831+557, and 0843+575. The second and larger class of bent jets shows a tendency for the electric field to follow the curve of the bending axis so that the electric vector position angles are usually, but not always, perpendicular to the bending axis of the jet. Figure 14 shows 0133+479, 0627+532, 1151+408, 1459+480, 1619+491, and 2351+456 which are all typical sources in this second category. A nice example of an $`𝐄`$-field bending with and parallel to the jet axis is the quasar 0627+532.

## 6 Conclusions and Future Work

A new, fairly strong correlation between core EVPA and jet axis angle has been discovered in quasars with projected jet lengths shorter than 6 mas. These objects exhibit $`𝐄`$-fields nearly perpendicular to the axis while those quasars with jet lengths longer than 6 mas show no obvious $`\theta \chi _{\mathrm{core}}`$ correlation. That this correlation exists in spite of the tendancy of quasar cores to have substantial Faraday RMs (Taylor 1998, 2000), implies that the actual correlation is even stronger than we observe. These findings suggest that quasars with axes anti-aligned with the $`𝐄`$-field are not intrinsically different from those without this correlation; rather, the long-jet sources may be at larger angles to the line-of-sight so that the angle to the core traverses a higher-density region with magnetic fields that produce sufficient Faraday rotation to smear out any intrinsic correlation. This reasoning also can explain the sparsely populated region of high projected jet lengths and high fractional core polarizations in Fig. 12. These results can be made consistent with the finding of Lister , of electric vectors aligned with the jet axis at 43 GHz if the cores are optically thick at 5 GHz, but optically thin at 43 GHz. A magnetic field perpendicular to the jet direction could be produced by a strong transverse shock \[e.g. Laing, 1980\].

Unlike the preferred alignment found in quasar cores and contrary to previous claims in the literature \[Cawthorne et al., 1993; Homan et al., 2002\], quasar jets are not found to exhibit a $`\theta \chi `$ relationship. This seems to be the case regardless of jet length, core dominance, or fractional polarization.

In contrast to quasar cores, BL Lac cores in our sample show no strongly preferred alignment. Since BL Lacs are observed to have lower Faraday Rotation Measures than quasars \[Zavala & Taylor, 2003\], it is unlikely that the absence of correlation is due to Faraday rotation. BL Lacs may be intrinsically different from quasars (i.e., have a lower power), or their intrinsic high variability may make it impossible to find a correlation over just one epoch, with only a modest number of sources (17). Future multi-epoch analysis, and/or the analysis of a larger sample, is necessary to conclusively determine the existence of any preferred alignment. An interesting test to make in the future would be to use spectral information to look for preferential alignments in optically thick or thin cores. Contrary to claims made by Gabuzda, Pushkarev & Cawthorne  of a bimodal $`\theta \chi _{\mathrm{jet}}`$ distribution, the EVPA in BL Lac jets similarly show no correlation with the jet axis angle.

We found no evolution of source properties (size, fractional polarization, or polarization angle) with redshift. This lack of evolution is reasonable if the physical parameters of the central engine outweigh any effects due to possible variations of intergalactic medium with redshift.

We thank Rene Vermeulen and Silke Britzen for help in scheduling and for sharing the data. LKP thanks NSF for summer support through the Research Experience for Underdergraduates (REU) program. We thank Travis Rector and Denise Gabuzda for helpful discussions. We also thank an anonymous referee for numerous helpful suggestions. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.


format: text

text:
# Thermodynamics of black holes in finite boxes

## 1 Introduction

In a seminal work of the 70ś, Hawking , showed that black holes are capable of emitting radiation, and evaluated the quantum process by which these black holes lose mass due to vacuum polarization (induced by the gravitational energy of the black hole). This important work and follow-up contributions, together with Bekenstein’s arguments  about the entropy of black holes, prompted a new discipline of black hole research: thermodynamics of black holes. Today, the study of black hole thermodynamics is very active and several arguments indicate that it must embrace the partial (or perhaps complete) marriage of quantum mechanics, thermodynamics and gravitation.

The original Hawking analysis showed that the particles created display a thermal spectrum $`f(E,T)=\frac{1}{e^{\beta E}\pm 1}`$ with $`\beta =1/k_BT`$, and a temperature determined by the black hole mass $`M`$ given by

$$T_{bh}=\frac{\mathrm{}c^3}{8\pi k_BGM}10^7\left(\frac{M_{}}{M}\right)K$$
(1)
see  for details. Since this spectrum is thermal, the luminosity of the black hole can be calculated using the Stefan-Boltzmann law, using $`LT^4r_{g}^{}{}_{}{}^{2}`$, where the gravitational radius $`r_gM`$ squared plays the role of the emitting area. Using these definitions, $`L(M)M^2`$ is readily obtained.

In addition to ”purely academic” cases (e.g. black holes inside ideal boxes), these emission (absorption) relations are in principle useful for investigating the properties of astrophysical and/or primordial black holes (hearafter PBHs) and their cosmological consequences. For example, some short-timescale gamma-ray burst events have been possibly explained by evaporating PBHs in a terminal (explosive) phase. In fact, the derived relations between period, mass and spectra are consistent with the PBH model, see .

Even in the simplest cases (the black holes interacting with ambient radiation) it is interesting to understand what determines the equilibrium features of the system and how is it achieved. We shall perform an analysis restricted to finite volumes with the aim of illustrating some novel features of the thermodynamics of black holes using gedanken experiments. These excercises serve also as a prelude to the more complicated situation, in which evaporation of PBHs must be studied in an expanding universe. However, and quite independently of these further advanced considerations, ”boxed” black holes display interesting features which need careful clarification.

## 2  Classical absorption in an initially cold box

Let us consider the following gedanken experiment: we start with a black hole and put it within a finite box (with linear size $`L`$). In the beginning, the box has no radiation, i.e. its temperature is zero. Then, since the black hole is hotter than the environment, it will emit thermal radiation, and thus this closed box will be increasingly filled with this radiation shortly afterwards.

The evolution of the mass of the black hole is described by the well-known expression which balances the energy loses to the Hawking luminosity $`M^2`$

$$\frac{dM}{dt}=\frac{C}{M^2}$$
(2)
with an initial condition of the box $`T_{rad}(t=t_i)=0`$. The constant $`C`$ depends on the degrees of freedom allowed to be radiated by the black hole  and is set to $`C10^{26}g^3s^1`$ throughout this work, and the initial mass is $`M_i`$. We assume that the box linear size is much larger than the black hole gravitational radius $`Lr_g`$ to avoid a complicated non-linear feedback between the geometry and the radiation field, which would lead us far beyond the purely thermodynamical approach.

Since the box is closed, the total energy is conserved. Therefore it is easy determine a relation for the black hole mass and the thermal radiation contained within the box from

$$\varrho _{rad}(T_{rad})V_{box}+M(t)=M_i$$
(3)
where $`t>t_i`$ (we use natural units to set the speed of light $`c`$ equal to 1). The radiation density is given by $`\varrho _{rad}(T)=a_{}T^4`$ with $`a_{}8\times 10^{36}gcm^3K^4`$. After some time, the box temperature is

$$T_{rad}(t)=\left[\frac{M_iM(t)}{a_{}V_{box}}\right]^{1/4}$$
(4)
At an even later time $`t>t_i`$, eq.(2) should not simply contain the Hawking term, but an absorption term is also needed. The simplest form of term is generally constructed as the product of the incoming flux times the gravitational cross-section of the radiation falling into the hole, $`B\varrho _{rad}(T)M^2`$, with $`B=\frac{27\pi G^2}{c^3}(or\frac{27\pi }{M_{pl}^{}{}_{}{}^{4}}`$ in natural units, see Ref.  for details and discussion). Therefore the complete differential equation for the mass is given by

$$\frac{dM}{dt}=\frac{C}{M^2}+B\varrho _{rad}(T)M^2$$
(5)
From this equation, it is clear that thermodynamical equilibrium between the black hole plus and the ambient radiation will be given by $`\dot{M}=0`$, and thus the mass of the black hole in equilibrium is

$$M_{bh}=M_c(t_{eq})=\frac{D}{T_{rad}(t_{eq})}$$
(6)
where $`D=(C/a_{}B)^{1/4}2\times 10^{26}gK`$. Solving the set of equations above we track the black hole mass evolution as it approaches the equilibrium (and determine the conditions for this equilibrium).

If we substitute the eq.(4) into eq.(5) we obtain

$$\frac{dM}{dt}=\frac{C}{M^2}+B\left(\frac{M_iM}{V_{box}}\right)M^2$$
(7)
The evolution of the mass of eq.(7) will be given by solving the integral

$$_{M_i}^{M(t)}\frac{dMM^2}{\kappa _1+\kappa _2M^4+\kappa _3M^5}=t$$
(8)
and the constants $`\kappa _i`$ are $`\kappa _1=C`$, $`\kappa _2=BM_i/V_{box}`$ and $`\kappa _3=B/V_{box}`$.

In the limit $`V_{box}>>r_{g}^{}{}_{}{}^{3}`$ the solution of eq.(7) approaches $`M(t)M_i[1(t/t_{evap})]^{1/3}`$ with $`t_{evap}=\frac{M_i^3}{3C}`$, as expected.

The numerical solution of the integral eq.(7) is shown in Fig.1 for different box sizes but the same initial black hole mass $`M_i`$. In order to obtain the conditions for the equilibrium we solve the algebraic equation $`T_{rad}(t_{eq})=T_{Haw}(M_{eq})`$, rewritten as

$$\left[\frac{M_iM(t_{eq})}{a_{}V_{box}}\right]^{1/4}=\frac{C}{M(t_{eq})}$$
(9)
Inserting the respective constants above, eq.(9) becomes

$$[1\mu (t_{eq})]\mu ^4(t_{eq})=D(V_{box}/cm^3)\times (M_{}/M_i)^5$$
(10)
where $`\mu (t_{eq})=\frac{M(t_{eq})}{M_i}`$ and $`D4.2\times 10^{97}`$.

The interpretation of the solution is quite evident, when the volume of the box is such that the right side of eq.(9) is larger than one, the black hole evaporates completely in the absence of additional quantum corrections at the Planck scale. But if $`D(V_{box}/cm^3)\times (M_{}/M_i)^5\mathrm{\hspace{0.17em}1}`$, the size of the box is small enough to stop the evaporation, and the system black hole$`+`$radiation achieves thermodynamical equilibrium before its mass vanishes completely.

The linear size of the box below which the black hole achieves its equilibrium is given by

$$L_c(M_i)1.3\times 10^{32}(M_i/M_{})^{5/3}cm$$
(11)
The time for reaching equilibrium $`t_{eq}`$ is determined by the initial mass and the size of the box only, i.e. $`t_{eq}=t_{eq}(M_i,V_{box})`$. For the same initial mass, smaller boxes will contain black holes with larger final masses at equilibrium (see Fig. 1). Our analysis agrees with the previous work by Page . From now on, we define the critical volume by $`V_{crit}(M)=L_c(M)^3`$ using eq.(10), with its clear physical interpretation given above. In the next sections we revisit the derivation of eq.(10), extend it to more complicated cases and consider causality features.

## 3  Classical absorption in a closed box with initial radiation

Let us now evaluate the behavior of the black hole mass when introduced in a box with a non-zero initial radiation content. We expect the black hole to achieve thermodynamical equilibrium earlier than in the case without initial radiation, and with a smaller equilibrium mass (considering boxes with the same size). In this case, is easy show that

$$_{M_i}^{M(t)}\frac{dMM^2}{\kappa _1+\kappa _4M^4(E_iM)}=t$$
(12)
where $`E_i=M_i+a_{}T_{i}^{}{}_{}{}^{4}V_{box}`$ and $`\kappa _4=B/V_{box}`$. The plot of $`M(t)`$ as a function of time is actually similar to the former case.

As before, we impose the equilibrium condition in to obtain the relation between the black hole mass and the box volume. Following the same steps as before, we obtain

$$\mu ^{}(t_{eq})^4[F(T_i,M_i)\mu ^{}(t_{eq})]=D(V_{box}/cm^3)\times (M_{}/M_i)^5$$
(13)
where $`F(T_i,M_i)=\frac{E_i}{M_i}`$, and the $`t_{eq}`$ is different from the previous case. Actually the black hole achieves thermodynamical equilibrium sooner, as expected.

## 4  $`N`$ black holes plus radiation in a closed box

The generalization to $`N`$ black holes seems quite straightforward, although it will become clear that the initial states and other details must be carefully defined to achieve consistent results. First, we shall consider just two black holes immersed in a closed box (with constant volume) in the following initial situation: one black hole has initial mass $`M_1`$ and the other black hole is more massive than the first, $`M_2>M_1`$. Initially the box has no radiation, and therefore these objects begin to evaporate immediately. Moreover, we shall consider that $`V_{box}V_{crit}(M_1)`$. Then at the initial time $`t_i`$ we have

$$\left(\frac{dM_2}{dt}\right)_{t_i}=\frac{C}{M_2^2}$$
(14)
$$\left(\frac{dM_1}{dt}\right)_{t_i}=\frac{C}{M_1^2}$$
(15)
Since $`M_1<M_2`$ one may consider $`|(\dot{M_1})_{evap}||(\dot{M_2})_{evap}|`$ initially. Afterwards, when the box starts to be filled with the emitted radiation, some energy can be absorbed by the black holes. Therefore, their evolution will be given by

$$\frac{dM_2}{dt}=\frac{C}{M_2^2}+B\varrho _{rad}(T)M_{2}^{}{}_{}{}^{2}$$
(16)
$$\frac{dM_1}{dt}=\frac{C}{M_1^2}+B\varrho _{rad}(T)M_{1}^{}{}_{}{}^{2}$$
(17)
$$M_1(t)+M_2(t)+V_{box}\varrho _{rad}(T)=M_1(t_i)+M_2(t_i)$$
(18)
where the last condition displays the conserved energy of the box. Using this constraint, we can describe this system by the following set of equations

$$\frac{dM_2}{dt}=\frac{C}{M_2(t)^2}+\frac{B}{V_{box}}M_2(t)^2[E_iM_1(t)M_2(t)$$
(19)
$$\frac{dM_1}{dt}=\frac{C}{M_1(t)^2}+\frac{B}{V_{box}}M_1(t)^2[E_iM_1(t)M_2(t)]$$
(20)
combining eqs.(17) and (18) and using the conservation of energy yields

$$C\left[\frac{1}{M_1^2}+\frac{1}{M_2^2}\right]+B\varrho _{rad}[M_1^2+M_2^2]+V_{box}\dot{\varrho }_{rad}=0$$
(21)
where $`E_i=M_1+M_2`$ is the initial energy.

Note that the set of coupled equations does not include the effect of black hole motion due to their mutual gravity. Motion is likely to affect the emission/absorption properties of the black holes in the relativistic regime , . Thus, the analysis is strictly valid whenever $`v_{bh}c`$.

Let us take a look at the global behavior of the solutions. First, it is easy to show that the equilibrium between two black holes is possible. If we impose that the thermodynamical equilibrium will be achieved, then asymptotically

$$\frac{dM_2}{dt}=\frac{dM_1}{dt}=0$$
(22)
From these equations $`\frac{C}{B}V_{box}=M_{1}^{}{}_{}{}^{4}[E_iM_1M_2]`$ immediately follows, and the same is true for the other black hole. Therefore, it follows that

$$M_{1,eq}=M_{2,eq}$$
(23)
The same results of these equilibrium masses could have been obtained from the expression of the entropy. In fact, more complicated cases can be worked out either using the above approach or using that the entropy must be an extreme for a system in equilibrium. For instance, let us consider $`N`$ black holes enclosed in the box. The total entropy for this system is given by

$$S_{total}(t)=\frac{1}{4}\underset{i}{}A_i(t)+g_{}T_{rad}^3(t)V_{box}$$
(24)
Since the horizon area for the i-th black hole is $`A_i=4\pi r_{g,i}^{}{}_{}{}^{2}`$ we have

$$S_{total}(t)=\frac{4\pi }{M_{pl}^4}\underset{i}{}M_i(t)^2+g_{}T_{rad}^3(t)V_{box}$$
(25)
In equilibrium $`T_{rad}=T_{bh1}=T_{bh2}=\mathrm{}=T_{bhi}=\frac{M_{pl}^2}{8\pi M_{eq}}`$, and therefore

$$S_{max}(N,V_{box})=\beta _{}N\mu _{eq}^{}{}_{}{}^{2}+l_{}/\mu _{eq}^{}{}_{}{}^{3}(V_{box}/cm^3)$$
(26)
where $`\beta _{}=2.5\times 10^{40}`$, $`l_{}=8\times 10^{33}\alpha _{}`$ and $`\mu _{eq}=(M_{eq}/10^{15}g)`$ the equilibrium mass scaled to a convenient reference value.

Extremizing the entropy $`\left(\frac{dS_{max}}{dM_{eq}}\right)=0`$ and solving for $`M_{eq}`$ yields

$$M_{eq}(N,V_{box})5.4\times 10^{13}\alpha _{}^{}{}_{}{}^{1/5}N^{1/5}(V_{box}/cm^3)^{1/5}g$$
(27)
with $`\alpha _{}=26.7g_{}`$. After substituting back into $`S_{max}`$ we find

$$S_{max}(N,V_{box})1.2\times 10^{38}\alpha _{}^{}{}_{}{}^{2/5}N^{3/5}(V_{box}/cm^3)^{2/5}$$
(28)
Since $`N=nV_{box}`$ by definition, we have $`S_{max}1.2\times 10^{38}n^{3/5}(V_{box}/cm^3)`$. Therefore, the thermodynamic approach leads us to conclude that the equilibrium black hole mass is bigger for bigger boxes. This result is quite counter-intuitive, but it has a reasonable explanation given below. To confirm the dependence first, we resort to another method.

Let us rewrite the equation for the mass of the i-th black hole as

$$\frac{dM_i}{dt}=\frac{C}{M_i(t)^2}+\frac{B}{V_{box}}M_i(t)^2[E_iNM_i(t)]$$
(29)
where we prepare the system with $`N`$ black holes with the same initial mass (this particular condition is not over-restrictive and simplifies the algebra). From $`\dot{M_i}=0`$ we obtain

$$\mu _{eq}^{}{}_{}{}^{4}\frac{1}{\mu _i}\mu _{eq}^{}{}_{}{}^{5}=\frac{K}{N\mu _i}(V_{box}/cm^3)$$
(30)
where $`K3.5\times 10^6`$. This algebraic equation has non-trivial solutions which can be found numerically. For the sake of definiteness we adopt $`N=1`$, $`\mu _i=10`$, and a small $`V_{box}=10^3cm^3`$, which yields $`\mu _{eq}0.137`$. If we enlarge the box, say $`V_{box}=10^4cm^3`$ instead, then $`\mu _{eq}0.244`$, and so on. If we fix the size of the box instead, say $`V_{box}=2.8\times 10^5cm^3`$, and change the number of black holes, the equilibrium mass will reflect the effect of the presence of other black holes accordingly. For example, taking $`\mu _i=10`$ and $`N=1`$ we have $`\mu _{eq}1.027`$ but for $`N=10^3`$ we obtain $`\mu _{eq}0.178`$ (Fig. 2). These results agree with the previous entropy analysis, and confirm that the coupled equations correctly describe the evolution of the system. It is quite interesting to propose the construction of a plot showing the behavior of the solutions and a discussion about the reasons for that behavior.

As a further related work it is proposed to show that if we consider the existence of an initial thermal radiation with energy $`E_i`$, the effective volume becomes

$$V_{eff}(N,M_i)=\frac{V_{box}}{N+E_i/M_i}$$
(31)
which holds for $`N`$ black holes with the same initial masses. This effective volume must be used instead of $`(V_{box}/N)`$ in eq.(26) and $`M_{eq}(N,V_{box})`$ becomes $`[\frac{V_{box}}{N+E_i/M_i}]^{1/5}`$.

The puzzling growth of the equilibrium mass in these systems must be further explained in more detail to improve the understanding. To proceed, let us consider two boxes with volumes $`V_{box,2}`$ and $`V_{box,1}`$ where the second box is larger than the first. Let us immerse two black holes with the same initial masses in both boxes. Therefore, the energy content is forced to be the same. Then, in equilibrium we will have (dropping some inessential numerical factors)

$$E_{box}=2M_{eq,1}+T_{rad,1}^4V_{box,1}=2M_{eq,2}+T_{rad,2}^4V_{box,2}$$
(32)
Since the equilibrium mass is given by $`M_{eq}[V_{box}]^{1/5}`$ for both boxes, then we have

$$V_{box,1}^{}{}_{}{}^{1/5}V_{box,2}^{}{}_{}{}^{1/5}=T_{rad,2}^4V_{box,2}T_{rad,1}^4V_{box,1}$$
(33)
both terms are negative due to the assumption $`V_{box,2}>V_{box,1}`$ then we obtain from the r.h.s. of this equation

$$\left(\frac{T_{rad,2}}{T_{rad,1}}\right)^4<\left(\frac{V_{box,1}}{V_{box,2}}\right)<1$$
(34)
then $`T_{rad,2}<T_{rad,1}`$. Therefore, the bigger box will have smaller radiation temperature. Since in equilibrium $`T_{rad,2}=T_{bh,2}\frac{1}{M_{eq,2}}`$ then the mass of the second black hole must be bigger than the first, in agreement with our previous analysis. Previous work by Page  (specifically his eq.(10)) has explored this problem, which has been solved here using the conservation of energy only, with the absorption terms explicitly present.

## 5 Causality considerations

The above considerations implicitly assumed that the black hole absorption starts instantaneously as the former is created inside the box. One might wonder whether a black hole immersed within an initially empty box is actually able to absorb ambient radiation immediately; since the Hawking radiation emitted by it will be available to absorption only after $`L/c`$ as required by causality.

In order to gain some insight, we use the fact that if the box is large enough, when the radiation comes back into the black hole, the object has evaporated completely. Then, a critical condition is obtained imposing this time interval to be bigger than the timescale for evaporation. Thus

$$\frac{L}{c}>\frac{M^3}{3C}$$
(35)
Therefore, it follows that $`L_{crit}(M)10^{29}cm\left(\frac{M}{10^{15}g}\right)^3`$.

Then, any box (without initial radiation) containing a black hole with mass $`M`$ whose linear size $`L`$ is bigger than the critical size $`L_{crit}`$ above would let the black hole to evaporate completely, thus precluding equilibrium. Note that the functional dependence and numerical value are very different from the naive approach used previously where we have considered that the absorption term turns on instantaneously (that is, as soon we immersed the black hole in the box). If we boldly compare $`L_{crit}`$ with the size of the Hubble horizon today $`10^{27}cm`$, we deduce that a PBHs with masses smaller than $`2\times 10^{14}g`$ can evaporate completely. However, it is clear that in a realistic case other sources of energy must be considered and our comment is intended just for pedagogical purposes.

The result of eq.(35) can be interpreted as follows. If we put one hot black hole (at $`t_i`$) with initial mass $`M_i`$ inside a closed box with volume $`V_{box}`$, it will evolve according to $`\dot{M}=C/M^2`$ for a time $`t_i<t<\frac{L_{box}}{c}`$. But from $`t>L_{box}/c`$ on, the black hole has to absorb radiation from the box (originated by itself!). Then the absorption term $`T^4V_{box}`$ contributes to avoid its complete evaporation. If the box is slightly larger (but smaller than the critical volume defined by eq.(10)) the onset of absorption is delayed, then this black hole will be hotter than the previous case. This black hole is then filling the box with higher temperature radiation, and after some time it will be absorbing radiation at higher temperature, hence the absorption term $`T^4M^2`$ will be larger. This term drives the black hole into the larger equilibrium mass, as given by eq.(27). It should be remembered that a maximal number for black holes enclosed in a box with volume $`V_{box}`$ must exist, a naive estimate of $`N_{max}`$ can be obtained from the close-packing condition $`N_{max}r_g^3(M_{eq})V_{box}`$.

## 6 Conclusions

In this work we have analyzed the behavior of black holes+radiation systems, considering some gedanken experiments within finite boxes (initially with and without radiation) (see  for a clear pioneering discussion of these issues). These gedanken experiments are important to understand the evolution of black holes in more complex situations, since it is hoped that we can always approximate the thermodynamical behavior of the universe with that of a finite box. In these simple cases the energy available for absorption by the black holes is finite at all times and there is no need to consider other sources. Strictly speaking, and with an eye on the actual cosmological formation and evolution of black holes, it may be argued that the problem is oversimplified, since actual primordial black hole masses must be formed with a substantial fraction of the horizon mass , and therefore finite size effects should be considered from the scratch. The actual situation is much more complicated but not impossible to tackle. The real difficulty, which can be viewed as an advanced exercise dealing with the (generalized) second law of thermodynamics not related to laboratory systems, is that (unlike laboratory gases) the black hole ”gas” will have varying masses and therefore a complicated kinetic equation must be used to describe their evolution. This provides a concrete example, possibly realized in nature, of a variable-mass gas needing a deep understanding of statistical mechanics principles for its very formulation. The simplest approach taken here was enough to show that one black hole plus radiation can achieve thermodynamical equilibrium if the box volume is smaller than a critical volume (in agreement with a previous treatment made by Page ). Some surprises arise in more general cases where we have two (or more) black holes (plus radiation) within these boxes.

Even in their simplest versions the cases of black holes in boxes are very instructive to analyze and may serve as an good issues for a graduate course of thermodynamics/statistical mechanics (see  for a discussion of related formal aspects of black hole thermodynamics) as an alternative to traditional problems dealing with the approach to equilibrium. They pose several questions of deep physical meaning, which are also strongly entangled and sometimes puzzling to interpret.

## 7 Appendix

A complementary set of simple problems intended to reinforce the conceptual aspects of black hole evaporation. Except for the problem 4, their mathematical complexity is quite straightforward.

1) Using eq.(2), prove that the timescale for evaporation for one black hole is proportional to $`M^3`$. Verify that for one black hole with the solar mass ($`M2\times 10^{33}g`$); $`t_{evap}(M)10^{65}years`$. Evaluate this for a initial mass $`10^{15}g`$ and compare with the age of the universe.

2) Verify the expression eq.(6) for the critical mass. Evaluate this number today, considering that $`T_{rad}(t_0)3K`$.

3) Derive eq.(11). Use for the linear size of the box the present size of the universe. Evaluate the corresponding mass at equilibrium.

4) Deduce the set of eqs.(19)-(21) and interpret the solutions graphically.

5) Deduce the $`M_{eq}(N,V_{box})`$ for the cases without initial radiation and with a radiation filling the box. Interpret these results and compare with the cosmological situation plugging representative numbers for the latter.

6) Discuss how all these results may be modified if the box expands following a known temporal dependence.


text:
# Note on ingoing coordinates for binary black holes

## I Introduction

In a previous paper alvi , an approximate solution to Einstein’s equations representing two widely separated nonrotating black holes in a circular orbit was constructed by matching a post-Newtonian metric to two perturbed Schwarzschild metrics. The spacetime metric was presented in a single coordinate system valid up to the holes’ apparent horizons. In this paper, I write the binary black hole four-metric from alvi  in coordinates that are corotating and post-Newtonian in the radiation and near zones and smoothly become ingoing near the black holes. This coordinate system is valid through the holes’ horizons and covers the holes’ interiors as well as the near and radiation zones. The metric components in this coordinate system are explicitly nonsingular on the black hole horizons. The metric presented here is promising as a source of initial data for numerical simulations of binary black holes. Since these simulations require the computational grid to extend inside the holes’ horizons, the coordinate system used near the black holes in alvi  is not suitable for numerical relativity.

Let us begin with the metric near the first black hole, that is, in region I in the terminology of alvi  (see Fig. 1 in alvi ). This metric is the Schwarzschild metric plus electric-type and magnetic-type tidal perturbations due to the second black hole, and is given in isotropic coordinates in Eq. (3.22) of alvi . The second black hole’s tidal field rotates with angular velocity $`\mathrm{\Omega }`$ as seen by inertial observers in the first black hole’s local asymptotic rest frame. However, the tidal perturbation’s angular dependence $`\varphi \mathrm{\Omega }\stackrel{~}{t}`$ as given in Eq. (3.22) of alvi  (I have replaced $`T`$ in that equation with $`\stackrel{~}{t}`$) is singular at the first black hole’s horizon. The reason is that the Schwarzschild time coordinate $`\stackrel{~}{t}`$, which is suitable for applying the technique of matched asymptotic expansions in the buffer zone around the black hole (see alvi  for details), is badly behaved at the horizon. Since our goal in this paper is to come up with coordinates valid through the horizon and inside the black hole, we must use a time coordinate $`T`$ with the property that the hypersurfaces of constant time coincide with Schwarzschild time slices in the buffer zone but smoothly transition into ingoing Eddington-Finkelstein time slices which penetrate the horizon. The singular angular dependence $`\varphi \mathrm{\Omega }\stackrel{~}{t}`$ can be simply replaced by the nonsingular $`\varphi \mathrm{\Omega }T`$, with $`T`$ as described above; this is discussed in further detail below.

It turns out that isotropic coordinates are not a good starting point for an ingoing transformation. The analog of the ingoing Eddington-Finkelstein transformation, which is based on ingoing null geodesics of the Schwarzschild spacetime, is unsuccessful when applied to isotropic coordinates: the coordinate system remains singular at the horizon. Indeed, the isotropic radial coordinate is only defined outside the black hole. However, isotropic coordinates were used in alvi  to match a tidally perturbed Schwarzschild metric to the post-Newtonian near zone metric. It is therefore necessary to define a new radial coordinate that is equal to the (tidally distorted) isotropic radial coordinate in the buffer zone but transitions smoothly into the (tidally distorted) Schwarzschild radial coordinate near the black hole.

## II Ingoing coordinates

Following the notation in alvi , I denote the black holes’ masses by $`m_1`$ and $`m_2`$, and their coordinate separation in post-Newtonian harmonic coordinates by $`b`$. Let $`m=m_1+m_2`$, $`ϵ=(m/b)^{1/2}`$, and $`\mathrm{\Omega }=(1m_1m_2/mb)(m/b^3)^{1/2}`$. By assumption, $`ϵ1`$.

Let us begin with the region I metric given in isotropic coordinates in Eq. (3.22) of alvi . Note that, in this paper, $`T`$ and $`R`$ denote the nonsingular time and radial coordinates described in Sec. I, while in alvi , they denoted the isotropic time and radial coordinates. Set $`\mathrm{\Omega }=0`$ in Eq. (3.22) of alvi  and transform to Schwarzschild coordinates $`(\stackrel{~}{t},\stackrel{~}{r},\theta ,\varphi )`$. This yields the metric $`\stackrel{~}{𝐠}=𝐠_S+\stackrel{~}{𝐡}`$; the Schwarzschild metric $`𝐠_S`$ and the stationary tidal perturbation $`\stackrel{~}{𝐡}`$ are given in Schwarzschild coordinates by

$`𝐠_S`$ $`=`$ $`\left(1{\displaystyle \frac{2m_1}{\stackrel{~}{r}}}\right)d\stackrel{~}{t}^2+\left(1{\displaystyle \frac{2m_1}{\stackrel{~}{r}}}\right)^1d\stackrel{~}{r}^2+\stackrel{~}{r}^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2),`$ (1)
$`\stackrel{~}{𝐡}`$ $`=`$ $`{\displaystyle \frac{4m_2}{b^3}}\sqrt{{\displaystyle \frac{m}{b}}}\left(1{\displaystyle \frac{2m_1}{\stackrel{~}{r}}}\right)\stackrel{~}{r}^3dt[\mathrm{cos}\theta \mathrm{sin}\varphi d\theta +\mathrm{sin}\theta \mathrm{cos}(2\theta )\mathrm{cos}\varphi d\varphi ]`$ (2)
$`+{\displaystyle \frac{m_2\stackrel{~}{r}^2}{b^3}}[3\mathrm{sin}^2\theta \mathrm{cos}^2\varphi 1][(1{\displaystyle \frac{2m_1}{\stackrel{~}{r}}})^2dt^2+d\stackrel{~}{r}^2`$
$`+(\stackrel{~}{r}^22m_1^2)(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2)].`$
In this notation, $`d\stackrel{~}{t}`$, $`d\stackrel{~}{r}`$, $`d\theta `$, and $`d\varphi `$ are coordinate one-forms and $`d\stackrel{~}{t}^2`$ denotes the tensor product $`d\stackrel{~}{t}d\stackrel{~}{t}`$.

Let $`\zeta _1`$ and $`\zeta _2`$ be two numbers satisfying $`2<\zeta _1<\zeta _2<(b/m_1)^{1/2}`$. Define the new ingoing coordinates $`(T,R,\theta ,\varphi )`$ by

$`\stackrel{~}{t}`$ $`=`$ $`T2m_1\mathrm{ln}\left({\displaystyle \frac{R}{2m_1}}1\right)\psi (R),`$ (3)
$`\stackrel{~}{r}`$ $`=`$ $`R+m_1\left(1+{\displaystyle \frac{m_1}{4R}}\right)\eta (R).`$ (4)
The functions $`\psi (R)`$ and $`\eta (R)`$ must satisfy the following properties: (i) for $`R\zeta _2m_1`$, $`\psi (R)=0`$ and $`\eta (R)=1`$ so the coordinates are isotropic; (ii) for $`R\zeta _1m_1`$, $`\psi (R)=1`$ and $`\eta (R)=0`$ so the time coordinate is Eddington-Finkelstein and the radial coordinate is Schwarzschild; and (iii) for $`\zeta _1m_1<R<\zeta _2m_1`$, $`\psi (R)`$ and $`\eta (R)`$ smoothly and monotonically vary between their constant values outside this interval. The transition points $`R=\zeta _1m_1`$ and $`R=\zeta _2m_1`$ and the functions $`\psi (R)`$ and $`\eta (R)`$ can be chosen freely as long as the above properties are satisfied. Since $`d\stackrel{~}{r}/dR=1+(m_1+m_1^2/4R)\eta ^{}(R)m_1^2\eta (R)/4R^21m_1^2/4R^2`$, we have $`d\stackrel{~}{r}/dR>0`$ for $`R>2m_1`$, and so Eqs. (3) and (4) define a valid coordinate transformation.

Transforming (1) and (2) using (3) and (4), we obtain

$`𝐠_S`$ $`=`$ $`f(dT^22\mathrm{\Psi }dTdR+\mathrm{\Psi }^2dR^2)+{\displaystyle \frac{H^2}{f}}dR^2+\stackrel{~}{r}^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2),`$ (5)
$`\stackrel{~}{𝐡}`$ $`=`$ $`{\displaystyle \frac{4ϵm_2}{b^3}}f\stackrel{~}{r}^3(dT\mathrm{\Psi }dR)[\mathrm{cos}\theta \mathrm{sin}\varphi d\theta +\mathrm{sin}\theta \mathrm{cos}(2\theta )\mathrm{cos}\varphi d\varphi ]`$ (6)
$`+{\displaystyle \frac{m_2}{b^3}}\stackrel{~}{r}^2(3\mathrm{sin}^2\theta \mathrm{cos}^2\varphi 1)[f^2(dT^22\mathrm{\Psi }dTdR+\mathrm{\Psi }^2dR^2)+H^2dR^2`$
$`+(\stackrel{~}{r}^22m_1^2)(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2)],`$
where

$`\mathrm{\Psi }(R)`$ $`=`$ $`2m_1\left[R^1\left(1{\displaystyle \frac{2m_1}{R}}\right)^1\psi (R)+\psi ^{}(R)\mathrm{ln}\left({\displaystyle \frac{R}{2m_1}}1\right)\right],`$ (7)
$`H(R)`$ $`=`$ $`1+m_1\left(1+{\displaystyle \frac{m_1}{4R}}\right)\eta ^{}(R){\displaystyle \frac{m_1^2}{4R^2}}\eta (R),`$ (8)
$`f(R)`$ $`=`$ $`1{\displaystyle \frac{2m_1}{\stackrel{~}{r}(R)}},`$ (9)
and $`\stackrel{~}{r}`$ is given in terms of $`R`$ in (4). Since the linearized Einstein equation is a tensor equation (see, e.g., Eq. (7.5.15) in wald ) and we have only performed a coordinate transformation, $`\stackrel{~}{𝐡}`$ remains a solution to this equation on the background $`𝐠_S`$. However, the $`\varphi `$-dependence of the perturbation $`\stackrel{~}{𝐡}`$ does not correspond to the second black hole’s rotating tidal field—we set $`\mathrm{\Omega }=0`$ above. To remedy this, we simply change $`\varphi `$ to $`\varphi \mathrm{\Omega }T`$ in the components of $`\stackrel{~}{𝐡}`$, which yields a new perturbation $`𝐡`$. Note that this replacement is not a coordinate transformation; a new tensor $`𝐡`$ is defined. Also note that $`T`$ becomes an ingoing coordinate near the horizon, so the time dependence $`\varphi \mathrm{\Omega }T`$ causes no problems on the horizon. This simple remedy works for the following reason: if we solve the linearized Einstein Eq. order by order in $`ϵ`$, then time derivatives of the components of $`𝐡`$ produce factors of $`m_1\mathrm{\Omega }ϵ^3`$ and can thus be neglected. The perturbation $`𝐡`$ is given by

$`𝐡`$ $`=`$ $`{\displaystyle \frac{4ϵm_2}{b^3}}f\stackrel{~}{r}^3(dT\mathrm{\Psi }dR)[\mathrm{cos}\theta \mathrm{sin}(\varphi \mathrm{\Omega }T)d\theta +\mathrm{sin}\theta \mathrm{cos}(2\theta )\mathrm{cos}(\varphi \mathrm{\Omega }T)d\varphi ]`$ (10)
$`+{\displaystyle \frac{m_2}{b^3}}\stackrel{~}{r}^2[3\mathrm{sin}^2\theta \mathrm{cos}^2(\varphi \mathrm{\Omega }T)1][f^2(dT^22\mathrm{\Psi }dTdR+\mathrm{\Psi }^2dR^2)+H^2dR^2`$
$`+(\stackrel{~}{r}^22m_1^2)(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2)].`$
We have now arrived at a metric $`𝐠=𝐠_S+𝐡`$, where $`𝐠_S`$ is given in (5) and $`𝐡`$ in (10), which is valid from the black hole’s interior up into and through the buffer zone around the hole, and is written in coordinates that are well behaved throughout this region.

## III Transformation to corotating coordinates

The next step is to transform the metric $`𝐠`$ to corotating post-Newtonian coordinates $`(t,x,y,z)`$ using the transformation given in Eqs. (4.22) and (4.23) of alvi . This transformation contains a rotation that can be performed by first defining $`\phi =\varphi \mathrm{\Omega }T`$ and then setting $`\mathrm{\Gamma }=R\mathrm{sin}\theta \mathrm{cos}\phi `$, $`\mathrm{\Lambda }=R\mathrm{sin}\theta \mathrm{sin}\phi `$, and $`Z=R\mathrm{cos}\theta `$. To complete the transformation, define the functions $`P_{\alpha \beta }(x,y,z)`$ for $`\alpha ,\beta =0,..,3`$ to be components of the metric $`𝐠`$ in coordinates $`(T,\mathrm{\Gamma },\mathrm{\Lambda },Z)`$; write these components as functions of $`(x,y,z)`$ using Eqs. (4.22) and (4.23) in alvi . The functions $`P_{\alpha \beta }`$ are given by

$`P_{00}`$ $`=`$ $`f+{\displaystyle \frac{\stackrel{~}{r}^2}{R^2}}\mathrm{\Omega }^2(\mathrm{\Gamma }^2+\mathrm{\Lambda }^2){\displaystyle \frac{4ϵm_2\stackrel{~}{r}^3}{b^3R^3}}(2Z^2R^2)f\mathrm{\Omega }\mathrm{\Gamma }`$
$`+{\displaystyle \frac{m_2\stackrel{~}{r}^2}{b^3R^2}}(3\mathrm{\Gamma }^2R^2)\left[f^2+{\displaystyle \frac{\mathrm{\Omega }^2}{R^2}}(\stackrel{~}{r}^22m_1^2)(\mathrm{\Gamma }^2+\mathrm{\Lambda }^2)\right],`$
$`P_{01}`$ $`=`$ $`P_{10}={\displaystyle \frac{\mathrm{\Gamma }}{R}}f\mathrm{\Psi }{\displaystyle \frac{\stackrel{~}{r}^2}{R^2}}\mathrm{\Omega }\mathrm{\Lambda }{\displaystyle \frac{2ϵm_2\stackrel{~}{r}^3}{b^3R^3}}f\left[\mathrm{\Gamma }\mathrm{\Lambda }{\displaystyle \frac{\mathrm{\Omega }}{R}}\mathrm{\Psi }\mathrm{\Gamma }^2(2Z^2R^2)\right]`$
$`+{\displaystyle \frac{m_2\stackrel{~}{r}^2}{b^3R^2}}(3\mathrm{\Gamma }^2R^2)\left[{\displaystyle \frac{\mathrm{\Gamma }}{R}}f^2\mathrm{\Psi }{\displaystyle \frac{\mathrm{\Omega }\mathrm{\Lambda }}{R^2}}(\stackrel{~}{r}^22m_1^2)\right],`$
$`P_{02}`$ $`=`$ $`P_{20}={\displaystyle \frac{\mathrm{\Lambda }}{R}}f\mathrm{\Psi }+{\displaystyle \frac{\stackrel{~}{r}^2}{R^2}}\mathrm{\Omega }\mathrm{\Gamma }{\displaystyle \frac{2ϵm_2\stackrel{~}{r}^3}{b^3R^3}}f\left[Z^2\mathrm{\Gamma }^2{\displaystyle \frac{\mathrm{\Omega }}{R}}\mathrm{\Psi }\mathrm{\Gamma }\mathrm{\Lambda }(2Z^2R^2)\right]`$
$`+{\displaystyle \frac{m_2\stackrel{~}{r}^2}{b^3R^2}}(3\mathrm{\Gamma }^2R^2)\left[{\displaystyle \frac{\mathrm{\Lambda }}{R}}f^2\mathrm{\Psi }+{\displaystyle \frac{\mathrm{\Omega }\mathrm{\Gamma }}{R^2}}(\stackrel{~}{r}^22m_1^2)\right],`$
$`P_{03}`$ $`=`$ $`P_{30}={\displaystyle \frac{Z}{R}}f\mathrm{\Psi }{\displaystyle \frac{2ϵm_2\stackrel{~}{r}^3}{b^3R^3}}f\left[\mathrm{\Lambda }Z{\displaystyle \frac{\mathrm{\Omega }}{R}}\mathrm{\Psi }\mathrm{\Gamma }Z(2Z^2R^2)\right]`$
$`+{\displaystyle \frac{m_2\stackrel{~}{r}^2}{b^3R^2}}(3\mathrm{\Gamma }^2R^2)\left({\displaystyle \frac{Z}{R}}f^2\mathrm{\Psi }\right),`$
$`P_{11}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }^2}{R^2}}\left(f\mathrm{\Psi }^2+{\displaystyle \frac{H^2}{f}}{\displaystyle \frac{\stackrel{~}{r}^2}{R^2}}\right)+{\displaystyle \frac{\stackrel{~}{r}^2}{R^2}}+{\displaystyle \frac{4ϵm_2\stackrel{~}{r}^3}{b^3R^4}}f\mathrm{\Psi }\mathrm{\Gamma }^2\mathrm{\Lambda }`$
$`+{\displaystyle \frac{m_2\stackrel{~}{r}^2}{b^3R^2}}(3\mathrm{\Gamma }^2R^2)\left[{\displaystyle \frac{\mathrm{\Gamma }^2}{R^2}}(f^2\mathrm{\Psi }^2+H^2)+R^2(\stackrel{~}{r}^22m_1^2)\left(1{\displaystyle \frac{\mathrm{\Gamma }^2}{R^2}}\right)\right],`$
$`P_{22}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Lambda }^2}{R^2}}\left(f\mathrm{\Psi }^2+{\displaystyle \frac{H^2}{f}}{\displaystyle \frac{\stackrel{~}{r}^2}{R^2}}\right)+{\displaystyle \frac{\stackrel{~}{r}^2}{R^2}}+{\displaystyle \frac{4ϵm_2\stackrel{~}{r}^3}{b^3R^4}}f\mathrm{\Psi }\mathrm{\Lambda }(Z^2\mathrm{\Gamma }^2)`$
$`+{\displaystyle \frac{m_2\stackrel{~}{r}^2}{b^3R^2}}(3\mathrm{\Gamma }^2R^2)\left[{\displaystyle \frac{\mathrm{\Lambda }^2}{R^2}}(f^2\mathrm{\Psi }^2+H^2)+R^2(\stackrel{~}{r}^22m_1^2)\left(1{\displaystyle \frac{\mathrm{\Lambda }^2}{R^2}}\right)\right],`$
$`P_{33}`$ $`=`$ $`{\displaystyle \frac{Z^2}{R^2}}\left(f\mathrm{\Psi }^2+{\displaystyle \frac{H^2}{f}}{\displaystyle \frac{\stackrel{~}{r}^2}{R^2}}\right)+{\displaystyle \frac{\stackrel{~}{r}^2}{R^2}}{\displaystyle \frac{4ϵm_2\stackrel{~}{r}^3}{b^3R^4}}f\mathrm{\Psi }Z^2\mathrm{\Lambda }`$
$`+{\displaystyle \frac{m_2\stackrel{~}{r}^2}{b^3R^2}}(3\mathrm{\Gamma }^2R^2)\left[{\displaystyle \frac{Z^2}{R^2}}(f^2\mathrm{\Psi }^2+H^2)+R^2(\stackrel{~}{r}^22m_1^2)\left(1{\displaystyle \frac{Z^2}{R^2}}\right)\right],`$
$`P_{12}`$ $`=`$ $`P_{21}={\displaystyle \frac{\mathrm{\Gamma }\mathrm{\Lambda }}{R^2}}\left(f\mathrm{\Psi }^2+{\displaystyle \frac{H^2}{f}}{\displaystyle \frac{\stackrel{~}{r}^2}{R^2}}\right)+{\displaystyle \frac{2ϵm_2\stackrel{~}{r}^3}{b^3R^4}}f\mathrm{\Psi }\mathrm{\Gamma }(\mathrm{\Gamma }^2+\mathrm{\Lambda }^2+Z^2)`$
$`+{\displaystyle \frac{m_2\stackrel{~}{r}^2}{b^3R^4}}(3\mathrm{\Gamma }^2R^2)\left[f^2\mathrm{\Psi }^2+H^2R^2(\stackrel{~}{r}^22m_1^2)\right]\mathrm{\Gamma }\mathrm{\Lambda },`$
$`P_{13}`$ $`=`$ $`P_{31}={\displaystyle \frac{\mathrm{\Gamma }Z}{R^2}}\left(f\mathrm{\Psi }^2+{\displaystyle \frac{H^2}{f}}{\displaystyle \frac{\stackrel{~}{r}^2}{R^2}}\right)`$
$`+{\displaystyle \frac{m_2\stackrel{~}{r}^2}{b^3R^4}}(3\mathrm{\Gamma }^2R^2)\left[f^2\mathrm{\Psi }^2+H^2R^2(\stackrel{~}{r}^22m_1^2)\right]\mathrm{\Gamma }Z,`$
$`P_{23}`$ $`=`$ $`P_{32}={\displaystyle \frac{\mathrm{\Lambda }Z}{R^2}}\left(f\mathrm{\Psi }^2+{\displaystyle \frac{H^2}{f}}{\displaystyle \frac{\stackrel{~}{r}^2}{R^2}}\right)+{\displaystyle \frac{2ϵm_2\stackrel{~}{r}^3}{b^3R^4}}f\mathrm{\Psi }Z(\mathrm{\Gamma }^2\mathrm{\Lambda }^2+Z^2)`$
$`+{\displaystyle \frac{m_2\stackrel{~}{r}^2}{b^3R^4}}(3\mathrm{\Gamma }^2R^2)\left[f^2\mathrm{\Psi }^2+H^2R^2(\stackrel{~}{r}^22m_1^2)\right]\mathrm{\Lambda }Z,`$
where $`ϵ=(m/b)^{1/2}`$, $`\mathrm{\Omega }=(1m_1m_2/mb)(m/b^3)^{1/2}`$, $`R=(\mathrm{\Gamma }^2+\mathrm{\Lambda }^2+Z^2)^{1/2}`$, and $`\mathrm{\Gamma }`$, $`\mathrm{\Lambda }`$, and $`Z`$ are to be expressed in terms of $`(x,y,z)`$ via Eqs. (4.22) and (4.23) in alvi . The functions $`\stackrel{~}{r}(R)`$, $`\mathrm{\Psi }(R)`$, $`H(R)`$, and $`f(R)`$ are given in Eqs. (4), (7), (8), and (9).

The remainder of the coordinate transformation from black hole coordinates $`(T,R,\theta ,\varphi )`$ to corotating post-Newtonian coordinates $`(t,x,y,z)`$ can be done exactly as in alvi . The final metric in region I (see Fig. 1 in alvi ) is given by Eq. (4.27) of alvi , but with $`P_{\alpha \beta }`$ taken from Eq. (LABEL:ingP) above. Define $`\overline{P}_{\alpha \beta }`$ to be $`P_{\alpha \beta }`$ (as given in (LABEL:ingP)) with $`m_1`$ and $`m_2`$ exchanged. Then the metric in region II (that is, near the second black hole; see Fig. 1 in alvi ) is given by Eq. (4.28) in alvi , but with $`\overline{P}_{\alpha \beta }`$ taken from here. Note that the final metric components are everywhere explicitly independent of time $`t`$.

To summarize, the expressions for the metric components given in Sec. V of alvi  are valid in ingoing coordinates if $`P_{\alpha \beta }`$ (and $`\overline{P}_{\alpha \beta }`$) are taken from Eq. (LABEL:ingP) and the functions $`\stackrel{~}{r}(R)`$, $`\mathrm{\Psi }(R)`$, $`H(R)`$, and $`f(R)`$ are taken from Eqs. (4), (7), (8), and (9) with $`\psi (R)`$ and $`\eta (R)`$ having properties (i)–(iii) given below Eq. (4). In order to show explicitly that these components are nonsingular at the black holes’ horizons, I write out $`P_{\alpha \beta }`$ for $`R\zeta _1m_1`$. In this region, $`\psi =1`$, $`\mathrm{\Psi }(R)=2m_1R^1(12m_1/R)^1`$, $`\eta =0`$, $`H=1`$, $`\stackrel{~}{r}=R`$, and $`f(R)=12m_1/R`$. Therefore, for $`R\zeta _1m_1`$, $`P_{\alpha \beta }`$ are given by

$`P_{00}`$ $`=`$ $`1+{\displaystyle \frac{2m_1}{R}}+\mathrm{\Omega }^2(\mathrm{\Gamma }^2+\mathrm{\Lambda }^2){\displaystyle \frac{4ϵm_2}{b^3}}\left(1{\displaystyle \frac{2m_1}{R}}\right)(2Z^2R^2)\mathrm{\Omega }\mathrm{\Gamma }`$
$`+{\displaystyle \frac{m_2}{b^3}}(3\mathrm{\Gamma }^2R^2)\left[\left(1{\displaystyle \frac{2m_1}{R}}\right)^2+\mathrm{\Omega }^2\left(1{\displaystyle \frac{2m_1^2}{R^2}}\right)(\mathrm{\Gamma }^2+\mathrm{\Lambda }^2)\right],`$
$`P_{01}`$ $`=`$ $`P_{10}={\displaystyle \frac{2m_1}{R^2}}\mathrm{\Gamma }\mathrm{\Omega }\mathrm{\Lambda }{\displaystyle \frac{2ϵm_2}{b^3}}\left[\left(1{\displaystyle \frac{2m_1}{R}}\right)\mathrm{\Gamma }\mathrm{\Lambda }{\displaystyle \frac{2m_1}{R^2}}\mathrm{\Omega }\mathrm{\Gamma }^2(2Z^2R^2)\right]`$
$`{\displaystyle \frac{m_2}{b^3}}(3\mathrm{\Gamma }^2R^2)\left[{\displaystyle \frac{2m_1}{R^2}}\left(1{\displaystyle \frac{2m_1}{R}}\right)\mathrm{\Gamma }+\mathrm{\Omega }\mathrm{\Lambda }\left(1{\displaystyle \frac{2m_1^2}{R^2}}\right)\right],`$
$`P_{02}`$ $`=`$ $`P_{20}={\displaystyle \frac{2m_1}{R^2}}\mathrm{\Lambda }+\mathrm{\Omega }\mathrm{\Gamma }{\displaystyle \frac{2ϵm_2}{b^3}}\left[\left(1{\displaystyle \frac{2m_1}{R}}\right)(Z^2\mathrm{\Gamma }^2){\displaystyle \frac{2m_1}{R^2}}\mathrm{\Omega }\mathrm{\Gamma }\mathrm{\Lambda }(2Z^2R^2)\right]`$
$`{\displaystyle \frac{m_2}{b^3}}(3\mathrm{\Gamma }^2R^2)\left[{\displaystyle \frac{2m_1}{R^2}}\left(1{\displaystyle \frac{2m_1}{R}}\right)\mathrm{\Lambda }\mathrm{\Omega }\mathrm{\Gamma }\left(1{\displaystyle \frac{2m_1^2}{R^2}}\right)\right],`$
$`P_{03}`$ $`=`$ $`P_{30}={\displaystyle \frac{2m_1}{R^2}}Z+{\displaystyle \frac{2ϵm_2}{b^3}}\left[\left(1{\displaystyle \frac{2m_1}{R}}\right)\mathrm{\Lambda }Z+{\displaystyle \frac{2m_1}{R^2}}\mathrm{\Omega }\mathrm{\Gamma }Z(2Z^2R^2)\right]`$
$`{\displaystyle \frac{2m_1m_2}{b^3R^2}}(3\mathrm{\Gamma }^2R^2)\left(1{\displaystyle \frac{2m_1}{R}}\right)Z,`$
$`P_{11}`$ $`=`$ $`1+{\displaystyle \frac{2m_1\mathrm{\Gamma }^2}{R^3}}+{\displaystyle \frac{8ϵm_1m_2}{b^3R^2}}\mathrm{\Gamma }^2\mathrm{\Lambda }+{\displaystyle \frac{m_2}{b^3}}(3\mathrm{\Gamma }^2R^2)\left[1{\displaystyle \frac{2m_1^2}{R^2}}\left(1{\displaystyle \frac{3\mathrm{\Gamma }^2}{R^2}}\right)\right],`$
$`P_{22}`$ $`=`$ $`1+{\displaystyle \frac{2m_1\mathrm{\Lambda }^2}{R^3}}+{\displaystyle \frac{8ϵm_1m_2}{b^3R^2}}\mathrm{\Lambda }(Z^2\mathrm{\Gamma }^2)+{\displaystyle \frac{m_2}{b^3}}(3\mathrm{\Gamma }^2R^2)\left[1{\displaystyle \frac{2m_1^2}{R^2}}\left(1{\displaystyle \frac{3\mathrm{\Lambda }^2}{R^2}}\right)\right],`$
$`P_{33}`$ $`=`$ $`1+{\displaystyle \frac{2m_1Z^2}{R^3}}{\displaystyle \frac{8ϵm_1m_2}{b^3R^2}}Z^2\mathrm{\Lambda }+{\displaystyle \frac{m_2}{b^3}}(3\mathrm{\Gamma }^2R^2)\left[1{\displaystyle \frac{2m_1^2}{R^2}}\left(1{\displaystyle \frac{3Z^2}{R^2}}\right)\right],`$
$`P_{12}`$ $`=`$ $`P_{21}={\displaystyle \frac{2m_1\mathrm{\Gamma }\mathrm{\Lambda }}{R^3}}+{\displaystyle \frac{4ϵm_1m_2}{b^3R^2}}\mathrm{\Gamma }(R^22\mathrm{\Gamma }^2)+{\displaystyle \frac{6m_1^2m_2}{b^3R^4}}(3\mathrm{\Gamma }^2R^2)\mathrm{\Gamma }\mathrm{\Lambda },`$
$`P_{13}`$ $`=`$ $`P_{31}={\displaystyle \frac{2m_1\mathrm{\Gamma }Z}{R^3}}+{\displaystyle \frac{6m_1^2m_2}{b^3R^4}}(3\mathrm{\Gamma }^2R^2)\mathrm{\Gamma }Z,`$
$`P_{23}`$ $`=`$ $`P_{32}={\displaystyle \frac{2m_1\mathrm{\Lambda }Z}{R^3}}+{\displaystyle \frac{4ϵm_1m_2}{b^3R^2}}Z(2Z^2R^2)+{\displaystyle \frac{6m_1^2m_2}{b^3R^4}}(3\mathrm{\Gamma }^2R^2)\mathrm{\Lambda }Z,`$ (12)
where $`ϵ=(m/b)^{1/2}`$, $`\mathrm{\Omega }=(1m_1m_2/mb)(m/b^3)^{1/2}`$, $`R=(\mathrm{\Gamma }^2+\mathrm{\Lambda }^2+Z^2)^{1/2}`$, and $`\mathrm{\Gamma }`$, $`\mathrm{\Lambda }`$, and $`Z`$ are to be expressed in terms of $`(x,y,z)`$ via Eqs. (4.22) and (4.23) in alvi . Note that the quantities $`P_{\alpha \beta }`$ in Eq. (12) are all finite at the horizon $`R=2m_1`$.

## Acknowledgments

I thank Lee Lindblom, Mark Scheel, Kip Thorne, and Michele Vallisneri for useful discussions. This research was supported in part by NSF grant PHY-9900776.


text:
# General relativistic gravitational field of a rigidly rotating disk of dust: Solution in terms of ultraelliptic functions

## a Introduction.

An infinitesimally thin disk of dust rotating uniformly around its symmetry axis is the simplest model of a rotating self–gravitating body. Within Newton’s theory of gravitation this model is described by the Maclaurin solution of the Laplace equation. In the present letter we show that the model of a uniformly rotating disk of dust allows for an explicit analytic solution in Einstein’s theory of gravitation as well. This seems to be the first rigorous global solution of the rotating body problem (for perfect fluids with $`p=0`$) in general relativity.

An approximate solution of the problem was given by Bardeen and Wagoner , . By applying the Inverse (Scattering) Method to the related boundary value problem of the Ernst equation we found the solution in terms of two linear integral equations, a ‘small’ and a ‘big’ one . The analytic solution to the small integral equation in terms of the Weierstrass function led to explicit expressions for the Ernst potential on the symmetry axis, the disk metric, and the surface mass density . In the following we will show that the complete solution of the problem may be represented, up to quadratures, in terms of ultraelliptic functions.

## b The metric.

The line element reads in Weyl–Lewis–Papapetrou–coordinates:

$$ds^2=e^{2U}[e^{2k}(d\rho ^2+d\zeta ^2)+\rho ^2d\phi ^2]e^{2U}(dt+ad\phi )^2.$$
(1)
(We use units where Newton’s gravitational constant $`G`$ as well as the velocity of light $`c`$ are equal to 1.) The metric functions $`e^{2U}`$, $`e^{2k}`$, and $`a/\rho _0`$ depend uniquely on the normalized coordinates $`\rho /\rho _0`$, $`\zeta /\rho _0`$, and the parameter

$$\mu =2\mathrm{\Omega }^2\rho _0^2e^{2V_0},$$
(2)
where $`\mathrm{\Omega }`$, $`\rho _0`$ and $`V_0`$ are the angular velocity, the coordinate radius, and the ‘surface potential’ $`V_o(\mu )U(\rho =0,\zeta =0,\mu `$), respectively. The disk is characterized by $`\zeta =0`$, $`0\rho \rho _0`$.

## c The Ernst potential.

The (complex) Ernst potential is defined by

$$f=e^{2U}+\mathrm{i}b,$$
(3)
where the imaginary part $`b`$ is related to the metric function $`a`$ according to

$$a_{,\rho }=\rho e^{4U}b_{,\zeta };a_{,\zeta }=\rho e^{4U}b_{,\rho }$$
(4)
i. e., $`a`$ may be calculated from a given Ernst potential (and its derivatives) by a quadrature. The same holds for the metric function $`e^{2k}`$. Hence, an axisymmetric stationary solution to the vacuum Einstein equations is sufficiently characterized by its Ernst potential. Note that our global solution is completely described by the vacuum solution outside the disk including the boundary data on the disk.

The Ernst potential for the dust disk solution is given by the expression

$$f=\mathrm{exp}\left\{\mu \left[\underset{X_1}{\overset{X_a}{}}\frac{X^2dX}{W}+\underset{X_2}{\overset{X_b}{}}\frac{X^2dX}{W}\underset{\mathrm{i}}{\overset{\mathrm{i}}{}}\frac{hX^2dX}{W_1}\right]\right\},$$
(5)
with

$$W=W_1W_2,W_1=\sqrt{(X\zeta /\rho _0)^2+(\rho /\rho _0)^2},W_2=\sqrt{1+\mu ^2(1+X^2)^2},$$
(6)
$$h=\frac{\mathrm{ln}\left(\sqrt{1+\mu ^2(1+X^2)^2}+\mu (1+X^2)\right)}{\pi \mathrm{i}\sqrt{1+\mu ^2(1+X^2)^2}},$$
(7)
$$X_1^2=\frac{\mathrm{i}\mu }{\mu },X_2^2=\frac{\mathrm{i}+\mu }{\mu }(\mathrm{}X_1<0,\mathrm{}X_2>0).$$
(8)
The upper integration limits $`X_a`$ and $`X_b`$ in the first two integrals are hyperelliptic functions of two variables $`u`$ and $`v`$ (i. e. ‘ultraelliptic’ functions) defined by

$$\underset{X_1}{\overset{X_a}{}}\frac{dX}{W}+\underset{X_2}{\overset{X_b}{}}\frac{dX}{W}=u,\underset{X_1}{\overset{X_a}{}}\frac{XdX}{W}+\underset{X_2}{\overset{X_b}{}}\frac{XdX}{W}=v,$$
(9)
where the variables $`u`$ and $`v`$ are given by

$$u=\underset{\mathrm{i}}{\overset{\mathrm{i}}{}}\frac{hdX}{W_1},v=\underset{\mathrm{i}}{\overset{\mathrm{i}}{}}\frac{hXdX}{W_1}.$$
(10)
The third integral in (5) as well as the integrals in (10) have to be taken along the imaginary axis in the complex X–plane with $`W_1`$ and $`h`$ fixed according to $`\mathrm{}W_1<0`$ (for $`\rho `$, $`\zeta `$ outside the disk) and $`\mathrm{}h=0`$.

Eqs. (9) have exactly the form of Jacobi’s famous inversion problem , which was solved by Göpel  and Rosenhain  who were able to express $`X_a(u,v)`$ and $`X_b(u,v)`$ in terms of (ultraelliptic) theta functions. In our case the moduli of the theta functions depend on $`\rho /\rho _0`$, $`\zeta /\rho _0`$, and $`\mu `$.

It should be noted that the integrations from $`X_1`$ to $`X_a`$ resp. $`X_2`$ to $`X_b`$ in (5) and (9) have to be performed along the same curve in the two–sheeted Riemann surface associated with $`W(X)`$. This leads to a unique Ernst potential $`f(\rho /\rho _0,\zeta /\rho _0,\mu )`$. More technical details and final expressions for the metric in terms of theta functions will be presented in a subsequent paper.

## d The Newtonian limit $`\mu 1`$.

At $`\mu =0`$, the solution may be expanded into a power series in $`\mu ^{1/2}`$:

$$f=1+\underset{n=1}{\overset{\mathrm{}}{}}f_n\mu ^{(n+1)/2}.$$
(11)
The coefficients $`f_n(\rho /\rho _0,\zeta /\rho _0)`$ turn out to be elementary functions. A particularly simple result is obtained for $`f_1`$ to $`f_4`$. From Eqs. (5) to (10) we get

$$X_aX_1=𝒪(\mu ^{3/2}),X_bX_2=𝒪(\mu ^{3/2}),$$
(12)
and

$$f=\mathrm{exp}\left\{\mu \underset{\mathrm{i}}{\overset{\mathrm{i}}{}}\frac{h(XX_1)(XX_2)dX}{W_1}+𝒪(\mu ^3)\right\}.$$
(13)
Expanding $`h`$, $`X_1`$ and $`X_2`$ one can easily determine $`f_1`$ to $`f_4`$. (Note that $`X=𝒪(1)`$ in this integral.) The result for $`f_1`$ is

$$f_1=\frac{1}{\pi \mathrm{i}}\underset{\mathrm{i}}{\overset{\mathrm{i}}{}}\frac{(1+X^2)dX}{W_1}=\frac{1}{\pi }\left\{\frac{4}{3}\mathrm{cot}^1\xi +\left[\xi (\xi ^2+\frac{1}{3})\mathrm{cot}^1\xi \right](13\eta ^2)\right\}$$
(14)
with elliptic coordinates $`\xi `$, $`\eta `$ defined by

$$\rho =\rho _0\sqrt{1+\xi ^2}\sqrt{1\eta ^2},\zeta =\rho _0\xi \eta ,0\xi <\mathrm{},1\eta 1.$$
(15)
Now, the Newtonian limit $`U_{Ml}`$ can be read off from the first order relation

$$g_{tt}=\mathrm{}f=(1+\mu f_1)=(1+\frac{2U_{Ml}}{c^2}).$$
(16)
Obviously, $`U_{Ml}=c^2\mu f_1/2`$ is the Maclaurin solution, i. e. the Newtonian solution of the problem of a rigidly rotating disk of dust. Note that in this limit

$$\mu =\frac{2U_{Ml}(\rho =0,\zeta =0)}{c^2}=\frac{2\mathrm{\Omega }^2\rho _0^2}{c^2}.$$
(17)
(For the sake of clarity we have reintroduced the velocity of light $`c`$ in the last two equations.)

## e Discussion.

In the parameter range $`0<\mu <\mu _0=4.62966184\mathrm{}`$  the solution is asymptotically flat and regular everywhere outside the disk. Its baryonic surface mass density is positive and vanishes at the rim of the disk. For $`\mu \mu _0`$ ($`\mathrm{\Omega }`$ finite) and $`\rho ^2+\zeta ^20`$ it approaches exactly the extreme Kerr solution (mass $`M=1/\mathrm{\hspace{0.17em}2}\mathrm{\Omega }`$). The meaning of the solution for $`\mu >\mu _0`$ deserves further discussion.

The analytic expressions above allow for a detailed investigation of characteristic properties, as for instance the formation of toroidal ergospheres, dragging effects, motion of test particles, and stability.