text:
# FS Aur - a new class of Cataclysmic Variables or the missing link between Intermediate Polars and SW Sex objects?

## 1 Introduction

FS Aur was discovered back in 1949 Hoffmeister (1949) as a dwarf nova. Howell & Szkody (1988) made an attempt to determine its orbital period, among several other systems, by means of CCD photometry. They cite a spectroscopic period of 85 min, but they failed to find any definite period from photometry, although they have detected 0.15 mag variations which could be consistent with 87 to 105 min periods. Later, Thorstensen et al. (1996) refined the orbital period from spectroscopy, and suggested a SU UMa classification based on the short orbital period of 0.0595 days (85.7 min).

Nevertheless, long-term monitoring of the system by several groups failed to detect any superoutburst/superhumps in its light curve. The publicly available VSOLJ and AAVSO light curves show instead a steady-cyclic outburst pattern that is more similar to a SS Cyg type dwarf nova lightcurve. Andronov (1991) found 12 outbursts on 180 archival photographic plates of the Sternberg State Astronomical Institute. But he also failed to find any super-outbursts in the archival data and notes the regular character of the outbursts, similar to dwarf novae.

Neustroev (2002), hereafter N02, among others, for a long time hunted for superoutbursts and also conducted a new study of the system in quiescence. In his paper he noted several interesting spectral characteristics of the object, roughly determined the system parameters and also discovered variability with a period longer than the orbital period.

Here we present the results of a new long-term photometric campaign and additional spectroscopic observations of FS Aur (§2). They are combined with data previously obtained by N02, separated by a 3 year time gap. We compare the properties of FS Aur with SW Sex stars (§3.2), and describe our discovery of the existence of a 3.42 hour photometric period in the light curve of FS Aur (§3.3). Finally, we discuss possible causes of this unusual behavior of the system (§4) and possible implications on the classification of Cataclysmic Variables(§5) .

## 2 Observations

### 2.1 Optical Observations

New observations of FS Aur were obtained at the Observatorio Astronomico Nacional (OAN SPM) in Mexico. The 2.1 m telescope with the B&Ch spectrograph, equipped with a SITe $`1024\times 1024`$ pixel CCD with a pixel size of 24 $`\mu \mathrm{m}`$, was used for the spectroscopic observations. Spectra were obtained in the second order of a 400 l/mm grating, covering the 3800–5200Å region and reaching 2.5 Å FWHM spectral resolution with the used slit width of 1.5′′. Spectral observations were done in two consecutive nights (31 Jan/1 Feb 2001) with a total coverage of 12 hours, consisting of 10 min duration individual exposures. He-Ar arc lamp exposures were taken before and after the object observations for wavelength calibration, and the standard spectrophotometric star G191B2B (Oke, 1990) was observed each night for flux calibration. The reduction of the images and the extraction of the spectra was done using standard IRAF111IRAF is distributed by the National Optical Astronomical Observatories, which is operated by the Associated Universities for Research in Astronomy, Inc., under contract to the National Science Foundation. routines customly used for long-slit spectra.

Photometric observations were performed with the 1.5 m telescope of OAN SPM, and cover a longer time span. Differential photometry was obtained with the three broad band Johnson-Cousins filters B V R in long time-series with individual exposures ranging from 20 to 180 sec. Selected comparison and check star form the field were used for differential photometry performed using the DAOPHOT package within IRAF.

Additional photometry was acquired by members of the Center of Backyard Astronomy (CBA)222http://cba.phys.columbia.edu/. For details of their activity and procedures used to obtain and reduce the data we refer to Patterson et al. (2000).

The log of all optical observations is presented in Table 1.

### 2.2 ROSAT X-Ray Observations

We also retrieved all FS Aur X-ray observations from the ROSAT archive at MPE Garching. FS Aur was observed in the course of the all-sky survey in August/September 1990 for a total exposure time of 441 sec, and is listed in the ROSAT All-Sky Survey Bright Source Catalog as 1RXS J054748.5+283512. It is detected at a count rate of 0.13$`\pm `$0.02 cts/s.

FS Aur was serendipitously observed again in pointed mode on 23 September 1993 with the position-sensitive proportional counter (PSPC, Pfeffermann et al. 1986) with a total exposure time of 9509 sec. Due to its off-axis angle, the effective exposure of FS Aur is 6130 sec. The mean count rate is 0.128$`\pm `$0.005 cts/s, identical to that of the all-sky survey observation in 1990.

The spectrum is identical in shape and brightness during both observations. The hardness ratio of the object inferred from these observations is HR1 = 0.911$`\pm `$0.027 and HR2 = 0.371$`\pm `$0.034, These hardness ratios are defined as HR1 = $`(BA)/(B+A)`$ and HR2 = $`(DC)/(D+C)`$, where $`A(0.10.4`$ keV), $`B(0.52.0`$ keV), $`C(0.50.9`$ keV), and $`D(0.92.0`$ keV) are the counts in the given energy range. Given the galactic coordinates of FS Aur of LII = 180$`.^{}`$55, BII = +0$`.^{}`$23, these hardness ratios can be interpreted as a rather hard source spectrum, with strong interstellar absorption which practically leaves no counts below 0.5 keV.

Due to this strong absorption and the restricted spectral coverage, several spectral models fit the data. The spectral fit using a thermal bremsstrahlung model is presented on Figure 1, and gives the following parameters: absorbing column $`N_\mathrm{H}`$ = (3.0$`\pm `$0.3)$`\times `$1021 cm-2, temperature $`kT`$ = 0.59$`\pm `$0.30 keV, and a normalization of (1.3$`\pm `$0.1)$`\times `$10-3 ph/cm2/s/keV. The absorbing column corresponds to about half the total galactic column in this direction (Dickey & Lockman 1990). The unabsorbed bolometric flux for this model is 1.0$`\times `$10-11 erg/cm2/s, and the luminosity is 1.2$`\times `$1031 (D/100 pc)2 erg/s.

## 3 Results

### 3.1 Spectroscopy

Our spectroscopic data basically repeat all characteristics described by N02. For radial velocity (RV) measurements we simply did one gaussian fitting to the mostly single peaked Hβ emission line profiles. Using RV measurements from our observations and combining them with the measurements from N02, we derived the spectroscopic (orbital) period of FS Aur. We used the DFT method to obtain the power spectrum, the spectral window corresponding to the data and the time distribution of our observations. The period analysis resulted in a strong peak confirming previous determination of the orbital period. Using the CLEAN procedure (Roberts et al., 1987) that allows to suppress alias frequency peaks caused by uneven data sampling, we obtained the strongest peak at the frequency corresponding to 85.65 (0$`\stackrel{\mathrm{d}}{\mathrm{.}}`$059479) minutes. This estimate is within the range of periods and margins of error cited by previous authors (Thorstensen et al., 1996). There is another peak at around 10 day-1. We will discuss this later in the text after considering our photometric data (§4). The superior conjunction of the emission line source (i.e. red-blue crossing) corresponds to HJD T0=2450101.2656. However, since the source of emission is not defined, it is not clear whether it is a superior conjunction of the WD or not.

Other emission lines show similar results with no significant phase shifts between the lines.

N02 estimates the inclination angle of the system to be within 51-65. While it is based on dynamical interpretation of emission lines which might be highly unreliable indicator of the motions of the underlying stars, it is in a good agreement with the velocities obtained from our new observations. It is surprising though not to see broad double peaked emission lines from this short-period system (Dwarf Nova) at this inclination. It is noteworthy, that He ii $`\lambda 4686`$ Å shows up strongly. It is unusual to have He ii this strong in an ordinary dwarf nova. It appears that the emission lines of higher excitation show central absorption cores at some phases while Hβ remains mostly single-peaked. Probably because there is no systematics behind this behaviour, we could not find any firm phase dependent pattern. We averaged the spectrum of FS Aur by co-adding all spectra after correcting for the orbital radial velocity (Figure 2). Individual lines are blown up in separate panels to display the abovementioned features. These and some other characteristics noted in FS Aur bring us to the idea that it may belong to the class of SW Sex objects.

### 3.2 Is FS Aur an SW Sex object?

Careful examination of spectra of FS Aur shows that:

* Inspite of the estimate of a relatively high inclination angle of the binary system of $`51<i<65`$, emission lines of FS Aur are single peaked. However, distorted wings can be seen in form of blue shifted shoulders around phase 0.8 (N02).
* Central small absorption can be seen in the higher Balmer lines starting from Hγ, and in the He i 4471 and He ii lines at some phases. The absorption core of He lines could be seen even in the integrated spectrum (see Figure 2).
* Although there are no eclipses, high excitation lines do notably diminish at some phases (N02).

These features are proper to SW Sex objects (Thorstensen et al., 1991; Hellier, 2000). Although most of the SW Sex stars are eclipsing systems (having been one of the original class identifiers), recently also low-inclination systems were discovered pertaining to that class. In particular, FS Aur shows remarkable similarities in the line profiles to LS Peg (Taylor et al., 1999) and V442 Oph (Hoard et al., 2000). The shoulder in the blue wing of Hβ around phase 0.8 is similar to those observed in LS Peg and V442 Oph, where it is identified with the high velocity component in these systems. We were not able to separate the high velocity component, probably because of the low S/N ratio of our data. Also worth mentioning is the relatively small intensity of the He ii line in both these low-inclination systems. Finally, there is a striking resemblance of the Doppler tomograms of FS Aur with those of SW Sex stars. The ring in the maps, corresponding to the accretion disk, is rarely seen in SW Sex systems. Instead, bright spots are often detected in the lower left quadrant or below the location of the WD, and they may vary from line to line. The velocity maps, the trailed and reconstructed spectra of Hβ and Hγ of FS Aur are presented in Figure 3. Similar tomograms were obtained by Neustroev (N02). Since the phasing in the absence of eclipses or any other clear indication of zero phase is arbitrary, the location of the spots on velocity map is arbitrary too.

Nevertheless, it should be noted that other characteristics of FS Aur do not fit well into the SW Sex definition. The latter are mostly novalike systems exhibiting high and low states of luminosity, while FS Aur undergoes $`2`$ mag dwarf nova like outbursts. Another significant difference between FS Aur and the SW Sex stars is its short orbital period. However, the spectral (orbital) period, unambiguously detected in FS Aur, appears to be in large contrast with the photometric period discovered by us.

### 3.3 Photometric Period

The extensive photometric coverage reveals a most unexpected result. The long-term variability noted in N02 in reality appears to be strictly periodic. All data that we were able to collect since then show the same period of 3$`.^\mathrm{h}`$4 and easily can be phase-folded over the time exceeding 3 years of observations. In Figure 4 we present several thousands of measurements collected from a variety of CBA small telescopes, in addition to SPM observations reported earlier in Tovmassian et al. (2002) and folded with the photometric period. The data were normalized to each night’s mean magnitude, since they were obtained with various detectors and comparison stars. In spite of the large 0$`.^\mathrm{m}`$2 scatter the sinusoidal shape of the light curve is mostly undisturbed. Neither the phase nor the amplitude of variations seem to be altered by the cyclic outbursts that the system underwent in between the different sets of observations. This demonstrates the very coherent nature of the periodic signal in the light curves.

However, when we look closely at each night’s data, we can see that the above mentioned periodic signal in the light curve varies in its strength and appearance. Since the observed phenomenon is of unusual nature we decided to present them all for clarity. In Figures 5, 6, 7, 8 each night’s data is presented in its corresponding panel. The overlayed solid curve represents the photometric period. At nights when the photometric period is not seen clearly or there are faster variations present, we have marked also the spectroscopic period with a dashed line. Its phase was estimated by a least-square fitting of the curve to the moment of minima and maxima measured at these nights. We are not sure that these rapid variations are the result of orbital modulations and these curves are presented solely for demonstration purposes.

In the broad-band filters obtained at the 1.5 m telescope of OAN SPM, the system shows smooth sinusoidal light curves. The data collected from a variety of CBA small telescopes is not as smooth, and additionally spontaneous wiggles are detected in the light curves of individual nights. However, the 3.4 hours periodic signal can be clearly distinguished in the light curves and unambiguously stands-up in the power spectrum. The most important characteristics of the light curves are:

* The period of variation in the light curve is 3$`.^\mathrm{h}`$4247$`\pm `$0$`.^\mathrm{h}`$0005 as estimated from the combined data (differential magnitudes of all dates and colors were mean-subtracted and combined in one large dataset). The power spectra of individual nights were also calculated and checked for consistency. Depending on the duration and data quality they all show a similar pattern within errors between individual nights. The CLEANed power spectrum of the combined data and a blowout of the most prominent peak are presented on Figure 9.
* The B V R light curves from OAN SPM are smooth and almost perfectly sinusoidal (see Figure 10).
* The amplitude of 0.24 mag is approximately equal in all three bands (B V R).
* The individual light curves from the unfiltered CBA observations are more “noisy”. By this we do not mean photon statistics or a data scatter due to the flickering which is usual for CVs. Instead, the light curves just show more erratic wiggling imposed on a longer term variation due to the 3$`.^\mathrm{h}`$4 period. However, the period is coherent over about 3 years. The phase-folded light curve of all combined data, presented in Figure 4, clearly demonstrates this coherence.
* The examination of other peaks in the power spectrum of the radial velocities shows an unidentified peak at 9.88 day-1, which is very close to the side-band frequency of

$$\frac{1}{\mathrm{P}_{sb}}=\frac{1}{\mathrm{P}_{orb}}\frac{1}{\mathrm{P}_{phot}}$$
(1)

The photometric period exceeds the spectroscopic by 2.4 times! This is an unprecedented case for a low mass binary system unambiguously identified as a Cataclysmic Variable. In order to show that the two periods do not match we plotted the combined photometric data folded with the spectroscopic (orbital) period (Figure 11), along with the radial velocity curve obtained from measurements of the Hβ line from our observations. The zero point for spectroscopic as well as photometric phase, the $`+/`$ crossing of the radial velocity, was taken from N02.

## 4 Possible causes of the large photometric period

1. The possibility of rotational effects in the reference frame of the secondary star as well as a triple/multiple component system can be effectively ruled out based on the canonical sinusoidal shape of the light curves, the equality of the amplitude of variation in the three bands and the absence of any spectral evidence of such kind.
2. As a dwarf nova with an orbital period below the period gap, FS Aur was expected to be a SU UMa type object and show superhumps due to an eccentric disc (see the extensive discussion in the recent paper by Patterson (2001) and references therein). The number of superhumpers is statistically significant, theory and models of superhumps are well developed. Based on our knowledge of this phenomenon one would expect to detect superhumps with periods exceeding the orbital only by a few percent. Noting that permanent superhumps were detected in novalike systems with short periods, a dwarf nova of SU UMa type produces superhumps only during superoutbursts. FS Aur in turn never showed superhumps during decades of optical monitoring. We think that it is extremely unlikely that the photometric period in FS Aur is due to an eccentric disc precession because the period difference is much too large.
3. The long period could be the precession period of a warped disc. One possible model to explain the difference between spectral (orbital) period and the photometric variations in FS Aur is the presence of a tilted accretion disc which is rigidly precessing with the observed photometric period. Similar models have been proposed to explain long-term variations in X-ray binaries. The analytical estimates of the induced precession period are given by Papaloizou & Terquem (1995) and Larwood (1998). The ratio between the binary orbital period and the forced precession period is

$$\frac{P}{P_p}=\frac{3}{7}q\left(\frac{1}{1+q}\right)^{1/2}\left(\frac{R_o}{a}\right)^{3/2}\mathrm{cos}\delta $$
(2)
where $`q`$ = Ms/Mp is the mass ratio, $`R_\mathrm{o}`$ is the outer radius of the disc, the orbital separation is $`a`$, and $`\delta `$ is the orbital inclination with respect to the disc. For mass ratios $`q=0.220`$, Larwood found that period ratios for precessing tilted discs is in the range $`P{}_{\mathrm{p}}{}^{}/P1840`$. If the mean Papaloizou & Pringle (1977) disc truncation radius applies, then the lower value is $`10`$.

From the absence of superoutbursts in FS Aur we are not bound to the 3:1 resonance criterium, so the mass ratio can be different to the $`0.05<q<0.25`$ requirement for the eccentric disc resonance. However, it is hard to find reasonable values of masses in order to reach a rigid precession of a gaseous Keplerian disc in a system similar to FS Aur. The secondary in the system should be late M dwarf according to the relation between orbital period and spectral class of the secondary. The absence of any contribution in the near infra-red spectrum from the secondary confirms that it is at least later than M 0 with M$`{}_{2}{}^{}<`$ 0.5M. On the other hand, rough estimates of M1 yield values somewhere between 0.34$`<M_1<0.46`$ (N02). Thus, the mass ratio $`q`$ can reach at maximum 1.5, if we assume an M0 secondary with 0.5 M and the lower limit on M1. With these parameters $`P_\mathrm{p}/P_{\mathrm{o}rb}>14`$ according to Larwood (1998) and the precession period in FS Aur should have been of the order of a day ($`20`$ hours)

Radiation driven warping of accretion discs (Oglive & Dubus, 2001) imposes additional restrictions on the parameters of binary systems in order to maintain stable precessing discs.

Thus, from a purely theoretical point of view it is very difficult to justify a disc precession model. From the observational point of view the contribution of the disc in emission is inferior to the spot, whatever the origin of the spot is. Thus, the spectroscopic data do not provide strong evidence of a luminous disc, whose precession could lead to the observed photometric variations.
4. Could it be due to the rotation of the magnetic poles of the white dwarf? The colors and shapes of the optical light curve strongly support the magnetic spin version of the photometric period; there are numerous similar light curves in the literature on intermediate polars (IPs). For example we refer to Norton (1999) and references therein. The intermediate polar scenario is so far the most plausible way to immerse unrelated to the orbital period in a close binary system.

With this interpretation, we run basically into one important issue: why is the period of rotation so slow? None of the dozens of known intermediate polars show white dwarf spin periods less than the orbital one. Yet, there are no theoretical restrictions on the velocity of rotation of the white dwarf in CVs. Generally, one would expect the rotation periods of white dwarfs in binary systems not to be very different from single WDs. According to Spruit (2001), the distribution of rotational velocities of WDs permits periods as the one observed in FS Aur.

However, the presence of a magnetic field anchored on a WD makes almost impossible the existance of a stable spin periods less than 0.68 Porb according to King & Wynn (1999) unless the magnetic field is high enough, so that the magnetosphere of the WD overfills its respective Roche lobe and locks the system into synchronous rotation. Somerscales et al. (2002) considered the spin periods of magnetic CVs and show that there are very strong restrictions on the spin periods in magnetic CVs. But we should note here that there are a few examples of asyncroneous polars, that are supposed to be jolted out of synchronism by a nova explosion. The period difference is again not so big as in the case of FS Aur, and any system would rapidly tend to establish synchronization, but little is known what happens during such cataclysms and how far it may take the spin period. Therefore, it is difficult to directly apply the intermediate polar scenario to the case of FS Aur in this pretext.

An assumption alternative to the longer than orbital spin period might be the possibility that the real white dwarf spin period is of the order of 50–100 sec and FS Aur is a fast rotator as AE Aqr and DQ Her. If we assume also that the white dwarf in FS Aur freely precesses then such fast rotation will result in the wobbling of the X-ray projector together with the WD. In a certain orientation of the magnetic pole angle to the orbital plane (colatitude $`\beta `$), the X-ray beam may or may not sweep across the disk or its inner edge. What, or how much that fan beam perpendicular to the magnetic field lines will illuminate may explain the nature of this so far mysteriuos photometric period. The idea of free precession of WD was explored by Leins et al. (1992) They mention intermadiate polars as a possible laboratory to test it. According to Leins et al. (1992), Euler’s frequency of free precession $`\omega _e`$ of a rigid and axially symmetric body will be:

$$\omega _e=\alpha (\rho _c)\mathrm{\Omega }^3$$
(3)
In case of an elastic body the Chandler frequency $`\omega _c`$ is used instead:

$$\omega _c=\beta (\rho _c)\mathrm{\Omega }^3$$
(4)
where $`\alpha _c(\rho _c)`$ and $`\beta _c(\rho )`$ are coefficients that are determined by the central density $`\rho _c`$ or a mass M, radius R, and moment of inertia A of the white dwarf. The results of numerical calculations for the realistic models of a white dwarf with different central densities are presented in the Table 2 of Leins et al. (1992). For a mass of a white dwarf between 0.55-0.63 M, $`\alpha _c(\rho _c)`$ ranges between 2.109 and 1.299, while $`\beta _c(\rho )`$ varies from 0.219 to 0.124, correspondingly. Thus, a white dwarf rotating with $`50`$ P$`{}_{spin}{}^{}100`$ sec would have a free precession period of P$`{}_{prec}{}^{}=205`$ min depending on which model (Euler or Chandler) we adopt.

A good test for this intermediate polar scenario is the X-ray light curve. The spectral hardness ratio is in agreement with X-ray spectra from other known IPs. We have folded the X-ray data of the ROSAT pointed observation with the photometric and spectroscopic period, respectively. The corresponding curves are shown in Figure 12. None of the two figures contradicts the presence of a periodic signal with the corresponding period in the light curve. FS Aur was observed 5 times with 10-20 min each, at time intervals of 3-4 hrs. The periodic pattern in both panels most probably is the observing window. This is what causes the left panel in the paper to show the ”coherent” pattern - it is just the observing window pattern. The statistics are very poor and the limit on the pulse fraction is only $`<`$80%. Therefore, depending on the binning of a handful of photons one can reach misleading results. New X-ray observations may are needed and may help to sort this out.

## 5 Conclusions

There were arguments in the literature as to whether SW Sex objects owe their unusual characteristics to the phenomenon known as intermediate polars (Warner, 1995). V795 Her has been proposed to be intermediate polar and SW Sex object at the same time (see Casares et al., (1996)). The newly discovered KUV 03580+0614 (Szkody et al., 2001) is a good candidate to qualify for both categories. There are other, less obvious examples. These objects lack the detection of a stable, clearly detectable spin period. In FS Aur, we were able to observe two distinct, unrelated periods similar to the intermediate polars. There are very firm restrictions on upper limits of the spin periods of IPs, and 205 min is certainly excluded for FS Aur as the period of rotation of the WD. Therefore, we propose that the WD in FS Aur is a fast rotator with Pspin around 50 to 100 sec and freely precesses with P$`{}_{prec}{}^{}=205`$ min. The cause of optical modulations with spin and side-band periods in IPs is still a matter of debate, but regardless of its nature the amplitude of the modulation is simply a geometrical effect of the orientation of the X-ray beam relative to the line of sight and/or orbital plane. Therefore, it may also explain why the observed phenomenon is unique so far. The number of IPs is still statistically small and the detection of this effect might be merely a matter of a distinct magnetic pole colatitude. We suggest that this phenomenological model can be observationally tested with fast optical photometry and/or more sensitive X-ray observations (Chandra/XMM).

While the combination of the SW Sex and intermediate polar models is appealing, there are big difficulties in squeezing FS Aur into each of these. Thus, FS Aur may be a rare, new kind of cataclysmic variable which can be defined as a dwarf nova, which in addition to its orbital period shows a 2-3 times longer photometric period in its light curve. Recently Woudt & Warner (2002) discoverd another similar system. GW Lib is another short period CV, that apart from its 1.28 hour spectroscopic period shows well defined photometric modulation of light with 2.09 hour periodicity. As in the case with FS Aur, the reasons for such behaviour are not yet clear. However we can talk of a new phenomenon within Cataclysmic Variables. FS Aur and GW Lib can be classified as new type of CVs, which still needs explanation of its nature.

This work was supported in part by CONACYT projects 25454-E, 36585-E.

| HJD start | Duration | Time of exposure | Band | Telescope |
| --- | --- | --- | --- | --- |
| 2450000+ | min | sec. |  |  |
| 1621.65 | 164 | 100 | R | 1.5m/SPM |
| 1622.64 | 184 | 180 | V | 1.5m/SPM |
| 1864.84 | 322 | 150 | V | 1.5m/SPM |
| 1865.80 | 298 | 150 | V | 1.5m/SPM |
| 1928.63 | 409 | 180 | B | 1.5m/SPM |
| 1941.63 | 368 | 600 | 3800-5200Å | 2.1m/SPM |
| 1942.62 | 299 | 600 | 3800-5200Å | 2.1m/SPM |
| 1946.65 | 315 | 30,20 | B,R | 1.5m/SPM |
| 2206.84 | 289 | 30 | B | 1.5m/SPM |
| 2259.78 | 157 | 30 | B | 1.5m/SPM |
| 2260.76 | 178 | 30 | B | 1.5m/SPM |

| SPM - San Piedro Martir Observatory, Mexico |
| --- |
| MDM - MDM Observatory |
| CBA - Center of Backyard Astronomy |


format: text

text:
# Stability Issues for 𝑤<-1 Dark Energy

## Abstract

Precision cosmological data hint that a dark energy with equation of state $`w=P/\rho <1`$ and hence dubious stability is viable. Here we discuss for any $`w`$ nucleation from $`\mathrm{\Lambda }>0`$ to $`\mathrm{\Lambda }=0`$ in a first-order phase transition. The critical radius is argued to be at least of galactic size and the corresponding nucleation rate glacial, thus underwriting the dark energy’s stability and rendering remote any microscopic effect.

preprint: April 2003
Introduction.

The equation of state for the dark energy component in cosmology has been the subject of much recent discussion  Present data are consistent with a constant $`w(Z)=1`$ corresponding to a cosmological constant. But the data allow a present value for $`w(Z=0)`$ in the range $`1.38<w(Z=0)<0.82`$ . If one assumes, more generally, that $`w(Z)`$ depends on $`Z`$ then the allowed range for $`w(Z=0)`$ is approximately the same. In the present article we shall forgo this greater generality as not relevant. Instead, in the present article we address the question of stability for a dark energy with constant $`w(Z)<1`$.

Interpretation as a limiting velocity

Consider making a Lorentz boost along the 1-direction with velocity $`V`$ (put c = 1). Then the stress-energy tensor which in the dark energy rest frame has the form:

$$T_{\mu \nu }=\mathrm{\Lambda }\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& w& 0& 0\\ 0& 0& w& 0\\ 0& 0& 0& w\end{array}\right)$$
(1)
is boosted to $`T_{\mu \nu }^{^{}}`$ given by

$$T_{\mu \nu }^{^{}}=\mathrm{\Lambda }\left(\begin{array}{cccc}1& V& 0& 0\\ V& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& w& 0& 0\\ 0& 0& w& 0\\ 0& 0& 0& w\end{array}\right)\left(\begin{array}{cccc}1& V& 0& 0\\ V& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right)=\mathrm{\Lambda }\left(\begin{array}{cccc}1+V^2w& V(1+w)& 0& 0\\ V(1+w)& V^2+w& 0& 0\\ 0& 0& w& 0\\ 0& 0& 0& w\end{array}\right)$$
(2)
We learn several things by studying Eq.(2). First, consider the energy component $`T_{00}^{^{}}=1+V^2w`$. Since $`V<1`$ we see that for $`w>1`$ this is positive $`T_{00}^{^{}}>0`$ and the Weak Energy Condition (WEC) is respected. For $`w=1`$, $`T_{00}^{^{}}0`$ as $`V1`$ and is still never negative. For $`w<1`$, however, we see that $`T_{00}^{^{}}<0`$ if $`V^2>(1/w)`$ and this violates the WEC and is the first sign that the case $`w<1`$ must be studied with great care. Looking at the pressure component $`T_{11}^{^{}}`$ we see the special role of the case $`w=1`$ because $`w=T_{11}^{^{}}/T_{00}^{^{}}`$ remains Lorentz invariant as expected for a cosmological constant. Similarly the off-diagonal components $`T_{01}^{^{}}`$ remain vanishing only in this case. The main concern is the negativity of $`T_{00}^{^{}}<0`$ which appears for $`V^2>(1/w)`$. One possibility is that it is impossible for $`V^2>(1/w)`$. The highest velocities known are those for the highest-energy cosmic rays which are protons with energy $`10^{20}eV`$. These have $`\gamma =(1V^2)^{1/2}10^{11}`$ corresponding to $`V110^{22}`$. This would imply that:

$$w>110^{22}$$
(3)
which is one possible conclusion.

First-Order Phase Transition and Nucleation Rate.

But let us suppose, as hinted at by  that more precise cosmological data reveals a dark matter which violates Eq.(3). Then, by boosting to an inertial frame with $`V^2>(1/w)`$, one arrives at $`T_{00}^{^{}}<0`$ and this would be a signal for vacuum instability. If the cosmological background is a Friedmann-Robertson-Walker (FRW) metric the physics is Lorentz invariant and so one should be able to see evidence for the instability already in the preferred frame where $`T_{\mu \nu }`$ is given by Eq.(1). This goes back to work in the 1960’s and 1970’s where one compares the unstable vacuum to a superheated liquid. As an example, at one atmospheric pressure water can be heated carefully to above $`100^0`$ C without boiling. The superheated water is metastable and attempts to nucleate bubbles containing steam. However, there is an energy balance for a three-dimensional bubble between the positive surface energy $`R^2`$ and the negative latent heat energy of the interior $`R^3`$ which leads to a critical radius below which the bubble shrinks away and above which the bubble expands and precipitates boiling. For the vacuum the first idea in  was to treat the spacetime vacuum as a four-dimensional material medium just like superheated water. The second idea in the same paper was to notice that a hyperspherical bubble expanding at the speed of light is the same to all inertial observers. This Lorentz invariance provided the mathematical relationship between the lifetime for unstable vacuum decay and the critical radius of the four-dimensional bubble or instanton.

In the rest frame, the energy density is

$$T_{00}=\mathrm{\Lambda }(10^3eV)^4(1\mathrm{m}\mathrm{m})^4$$
(4)
since $`10^3eV(1\mathrm{m}\mathrm{m})^1`$.

In order to make an estimate of the dark energy decay lifetime in the absence of a known potential, we can proceed by assuming (without motivation from observation) that there is a first-order phase transition possible between the $`\mathrm{\Lambda }=(10^3eV)^4`$ “phase” and a $`\mathrm{\Lambda }=0`$ “phase”. This hypothesized decay is the Lorentz invariant process of a hyperspherical bubble expanding at the speed of light, the same for all inertial observers. Let the radius of this hypersphere be R, its energy density be $`ϵ`$ and its surface tension be $`S_1`$. Then according to  the relevant instanton action is

$$A=\frac{1}{2}\pi ^2R^4ϵ+2\pi ^2R^3S_1$$
(5)
where $`ϵ`$ and $`S_1`$ are the volume and surface energy densities, respectively. The stationary value of this action is

$$A_m=\frac{27}{2}\pi ^2S_1^4/ϵ^3$$
(6)
corresponding to the critical radius

$$R_m=3S_1/ϵ$$
(7)
We shall assume that the wall thickness is negligible compared to the bubble radius. The number of vacuum nucleations in the past lightcone is estimated as

$$N=(V_u\mathrm{\Delta }^4)exp(A_m)$$
(8)
where $`V_u`$ is the 4-volume of the past and $`\mathrm{\Delta }`$ is the mass scale relevant to the problem. This vacuum decay picture led to the proposals of inflation, for solving the horizon, flatness and monopole problems (only the horizon problem was generally known at the time of ). None of that work addressed why the true vacuum has zero energy. Now that the observed vacuum has non-zero energy density $`+ϵ(10^3eV)^4`$ we may interpret it as a false vacuum lying above the true vacuum with $`ϵ=0`$. In order to use the full power of Eq.(8) taken from  and the requirement $`N1`$ we need to estimate the three mass-dimension parameters $`ϵ^{1/4},S_1^{1/3}`$ and $`\mathrm{\Delta }`$ therein and so we discuss these three scales in turn.

The easiest of the three to select is $`ϵ`$. If we imagine a tunneling through a barrier between a false vacuum with energy density $`ϵ`$ to a true vacuum at energy density zero then the energy density inside the bubble will be $`ϵ=\mathrm{\Lambda }=(10^3eV)^4`$. No other choice is reasonable. As for the typical mass scale $`\mathrm{\Delta }`$ in the prefactor of Eq. (8), the value of $`\mathrm{\Delta }`$ does not matter very much because it appears in a power rather than an exponential so let us put (the reader can check that the conclusions do not depend on this choice) $`\mathrm{\Delta }=ϵ^{1/4}=(1mm)^1`$ whereupon the prefactor in Eq.(8) is $`10^{116}`$. The third and final scale to discuss is the surface tension, $`S_1`$. Here we appeal to comparison of spontaneous decay to stimulated decay. The former dictates $`N1`$ in Eq.(8): the latter requires further discussion.

Spontaneous dark energy decay brings us to the question of whether such decay can be initiated in an environment existing within our Universe. The question is analogous to one of electroweak phase transition in high energy particle collision. This was first raised in  and revisited for cosmic-ray collisions in . That was in the context of the standard-model Higgs vacuum and the conclusion is that high-energy colliders are safe at all present and planned foreseeable energies because much more severe conditions have already occurred (without disaster) in cosmic-ray collisions within our galaxy. More recently, this issue has been addressed in connection with fears that the Relativistic Heavy Ion Collider (RHIC) might initiate a diastrous transition but according to careful analysis there was no such danger.

The energy density involved for dark energy is some 58 orders of magnitude smaller \[$`(10^3eV)^4`$ compared to $`(300GeV)^4`$\] than for the electroweak case and so the nucleation scales are completely different. One is here led away from microscopic towards astronomical size scales.

The energy density of Eq.(4) is so readily exceeded that the critical radius cannot be microsopic. Think first of a macroscopic scale e.g.  1 meter and consider a magnetic field practically-attainable in bulk on Earth such as 10 Tesla. Its energy density is given by

$$\rho _{mag}=\frac{1}{\mu _0}^2$$
(9)
Using the value $`\mu _0=4\pi \times 10^7NA^2`$ and $`1T=6.2\times 10^{12}(MeV.s.m^2)`$ leads to an energy density $`\rho _{mag}=2.5\times 10^{17}eV/(mm)^3`$, over 20 orders above the value of Eq.(4) for the interior of nucleation. Magnetic fields in bulk exist in galaxies with strength $`1\mu G`$ and the rescaling by $`^2`$ then would give $`\rho _{mag}(2.5\times 10^5eV)(mm)^3`$, slightly below the value in Eq.(4). Assuming the dark energy can exchange energy with magnetic energy density the observed absence of stimulated decay would then imply a critical radius of at least galactic size, say, $`10kpc.`$ Using Eq(7) then gives for the surface tension $`S_1>10^{23}(mm)^3`$ and number of nucleations in Eq.(8) $`N<exp(10^{92})`$. The spontaneous decay is thus glacial. Note that the dark energy has appeared only recently in cosmological time and has never interacted with background radiation of comparable energy density. Also, this nucleation argument does not require $`w<1`$.

Discussion

As a first remark, since the critical radius $`R_m`$ for nucleation is astronomical, it appears that the instability cannot be triggered by any microscopic process. While it may be comforting to know that the dark energy is not such a doomsday phenomenon, it also implies at the same time the dreadful conclusion that dark energy may have no microscopic effect. If any such microscopic effect in a terrestrial experiment could be found, it would be crucial in investigating the dark energy phenomenon. We note that the present arguments are less model-dependent than those in .

In closing one may speculate how such stability arguments may evolve. One may expect most conservatively that the value $`w=1`$ will eventually be established empirically in which case both quintessence and the “phantom menace” will be irrelevant. In that case, indeed for any $`w`$, we may still hope that dark energy will provide the first connection between string theory and the real world as in e.g. . Even if precise data do establish $`w<1`$, as in the “phantom menace” scenario, the dark energy stability issue is still under control.

Acknowledgements

We thank Stuart Raby and Tomo Takahashi for useful discussions. This work is supported in part by the US Department of Energy under Grant No. DE-FG02-97ER-41036.


text:
# Boundary conditions for hyperbolic formulations of the Einstein equations

## 1 Introduction

Achieving long time simulations of binary black hole spacetimes remains one of the most pressing issues in numerical relativity. Numerical simulations are plagued with instabilities, the origin of which has not been isolated among a number of possible factors of relevance to stability (see for instance  and references therein). A partial list of suspected sources of instabilities includes ill-posed evolution, ill-posed constraint propagation and poor choices of boundary conditions, as well as poor choices of binary black hole data. The problem of well posedness of the evolution equations and constraint propagation in the analytic (as opposed to numerical) sense has been widely studied (see the Living Review article by Reula ), with the result that there is a great deal of choice of strongly hyperbolic formulations of the Einstein equations available for analytic studies (aside from numerical implementation). On the other hand, appropriate boundary values for evolution remain to be identified. Since a great deal more about the boundary value problem of a strongly hyperbolic set of equations is known than for any other kind, currently analytic investigations of appropriate boundary conditions tend to use some kind of strongly hyperbolic formulation of the Einstein equations. Relevant considerations usually run along the following lines.

The boundary values of a strongly hyperbolic system of partial differential equations split always into two sets : those that are determined by the initial conditions of the problem and cannot be chosen freely, and those that are entirely arbitrary because they are determined by values of the fields outside of the region of interest. In practice, the free boundary data must be specified in order for a unique solution of the equations to be determined in the region of interest.

Suppose we have a boundary delimiting a region where we seek a solution to the Einstein equations in some strongly hyperbolic formulation. In this case one calculates data on the initial slice satisfying the constraints, and then integrates up the hyperbolic equations step by step. But, by the previous paragraph, the initial data alone are not sufficient to find a unique solution in the future of a bounded region in space. Boundary values must be prescribed as well. Suppose that we know exactly which variables’ boundary values are to be considered independent of the initial data. How are we to prescribe them? Since in some sense the solution that we are seeking must satisfy all the Einstein equations, and by construction the evolution equations are being imposed, then it seems that the appropriate question to ask is: What are the boundary data that preserve the constraints?

This question is addressed by Stewart in  within a particular strongly hyperbolic formulation of the Einstein equations . The constraints themselves, considered as real functions of spacetime, propagate according to a strongly hyperbolic system of their own. This implies that vanishing values of the constraints at the initial time will propagate towards the future along characteristics. Some of the constraints propagate towards the boundary and cross out of the region of interest, whereas others propagate into the region of interest by crossing the boundary from outside. Clearly the values of the incoming constraints at the boundary are arbitrary, and one wants to have them vanish. But the vanishing of the constraints cannot be imposed along the boundaries in practice. The constraints involve derivatives of the fields across the boundaries, not just the values of the fields themselves. Stewart argues that the vanishing of the incoming constraints can in fact be used as a boundary condition on the fields when the equations are linearized around flat space. If the fields are expressed in integral form in terms of Fourier-Laplace transforms, the linearity of the differential equations implies algebraic equations for the transforms of the fields. Additionally, the constraints transform into algebraic expressions in terms of the transforms of the fields, thus acquiring an algebraic look. Regardless of whether these seemingly algebraic conditions are practical boundary conditions for a numerical simulation, the argument does not hold up in the non-linear case.

The idea of imposing the vanishing of the incoming constraints as boundary conditions is pursued further in , where space derivatives of the fields are eliminated in favor of time derivatives in the expression of the incoming constraints in terms of the fundamental variables. In this case, a different formulation of the Einstein equations is used , restricted to spherical symmetry.

Clearly not all formulations have the problem of imposing the constraints at the boundary. Suppose there is a formulation of the Einstein equations that has no incoming constraints, that is: a formulation in which the constraints propagate upwards along the boundary. In such a case, the constraints will be satisfied at the boundary by virtue of the initial values alone, and there is no need to impose additional conditions on the boundary values in order to enforce them. The boundary values of the incoming fields in this case must be arbitrary. Such formulations exist at least in symmetry-reduced cases, and some boundary problems for those formulations have been studied .

Here we stray away from the general trend of seeking a way to impose the constraints along the boundary, in order to propose a seemingly unrelated method to write down equations that must hold among the boundary values of many quite generic first-order formulations of the Einstein equations. These consist of the vanishing of the four components of the projection of the Einstein tensor along the normal to the boundary.

In Section 2 we briefly describe the formulation of the Einstein equations that we choose to write down the proposed boundary conditions explicitly in a model case. The discussion in Section 2 hinges heavily on the existence of a complete set of characteristic fields, which is guaranteed by the strong hyperbolicity of the formulation chosen (as opposed to weak hyperbolicity). In Section 3 we show that such boundary conditions are consistent with constraint propagation, realize the goal of guaranteeing the vanishing of the incoming constraints at the boundary and, furthermore, they coincide with the boundary conditions proposed in  by the trading of space and time derivatives alluded to above. The discussion on Section 3 hinges heavily on the existence of a complete set of characteristic constraint fields, namely: on the fact that the constraint propagation is strongly hyperbolic (as opposed to weakly hyperbolic). In Section 4 we state the generality of the arguments with regards to three-dimensional strongly hyperbolic formulations, namely: the aspects of the argument that apply to any strongly hyperbolic formulation of the Einstein equations irrespective of particularities such as the specific characteristic fields and speeds. Formulation-dependent features such as any explicit form of the boundary conditions arising from the projection of the Einstein equations perpendicularly to a boundary are intentionally excluded from Section 4 with the purpose of exposing the true reach and relevance of the argument across the lines delimiting alternative strongly hyperbolic formulations. Concluding remarks appear in Section 5.

## 2 Boundary conditions for the Einstein-Christoffel formulation with spherical symmetry

As usual, in a foliation of spacetime by level surfaces of a time function $`t`$, and using $`t`$ as a coordinate, the spacetime metric can be viewed as a time dependent Riemannian metric evolving in time according to a freely specifiable rate and labeling of the 3-space. From now on, we use the term metric with no qualifier to refer to the evolving Riemannian metric of the slices. In this section we use the Einstein-Christoffel (EC) formulation , restricted to spherical symmetry. The full 3-D equations and their structure of characteristics can be found in , as well as the restriction to spherical symmetry. We borrow the notation directly from . With the restriction to spherical symmetry, the spacetime metric has the form

$$ds^2=N^2dt^2+g_{rr}(dr+\beta ^rdt)^2+g_Tr^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2)$$
(1)
where $`g_{rr},g_T,\beta ^r`$ and $`N`$ are functions only of $`r`$ and $`t`$. The shift $`\beta ^r`$ and the densitized lapse

$$\stackrel{~}{\alpha }\frac{N}{g_T\sqrt{g_{rr}}}$$
(2)
are considered arbitrarily given. The EC formulation of the Einstein equations restricted to spherical symmetry consists of six fundamental variables $`(g_{rr},g_T,K_{rr},K_T,f_{rrr},f_{rT})`$ which evolve according to six evolution equations and whose initial data satisfies four constraints . The evolution equations are

$`_tg_{rr}\beta ^r_rg_{rr}`$ $`=`$ $`2NK_{rr}+2g_{rr}_r\beta ^r,`$ (3a)
$`_tg_T\beta ^r_rg_T`$ $`=`$ $`2NK_T+2{\displaystyle \frac{\beta ^r}{r}}g_T,`$ (3b)
$`_tK_{rr}\beta ^r_rK_{rr}+{\displaystyle \frac{N}{g_{rr}}}_rf_{rrr}`$ $`=`$ $`N[2f^r{}_{rr}{}^{}(f^r{}_{rr}{}^{}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{4f_{rT}}{g_T}})+K_{rr}(2{\displaystyle \frac{K_T}{g_T}}K_r^r)`$ (3c)
$`{\displaystyle \frac{6}{r^2}}6\left({\displaystyle \frac{f_{rT}}{g_T}}\right)^2_r^2\mathrm{ln}\stackrel{~}{\alpha }(_r\mathrm{ln}\stackrel{~}{\alpha })^2`$
$`+({\displaystyle \frac{4}{r}}f^r{}_{rr}{}^{})_r\mathrm{ln}\stackrel{~}{\alpha }]+2K_{rr}_r\beta _r`$
$`+4\pi N(Tg_{rr}2S_{rr}),`$
$`_tK_T\beta ^r_rK_T+{\displaystyle \frac{N}{g_{rr}}}_rf_{rT}`$ $`=`$ $`N\left(K_TK_r^r+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{2f_{rT}^2}{g_{rr}g_T}}{\displaystyle \frac{f_{rT}}{g_{rr}}}_r\mathrm{ln}\stackrel{~}{\alpha }\right)`$ (3d)
$`+{\displaystyle \frac{2\beta ^r}{r}}K_T,`$
$`_tf_{rrr}\beta ^r_rf_{rrr}+{\displaystyle \frac{N}{g_{rr}}}_rK_{rr}`$ $`=`$ $`N[4g_{rr}{\displaystyle \frac{K_T}{g_T}}(3{\displaystyle \frac{f_{rT}}{g_T}}f^r{}_{rr}{}^{}+{\displaystyle \frac{2}{r}}_r\mathrm{ln}\stackrel{~}{\alpha })`$ (3e)
$`K_{rr}(10{\displaystyle \frac{f_{rT}}{g_T}}+f^r{}_{rr}{}^{}{\displaystyle \frac{2}{r}}+_r\mathrm{ln}\stackrel{~}{\alpha })]`$
$`+3f_{rrr}_r\beta ^r+g_{rr}_r^2\beta ^r+16\pi NJ_rg_{rr},`$
$`_tf_{rT}\beta ^r_rf_{rT}+{\displaystyle \frac{N}{g_{rr}}}_rK_T`$ $`=`$ $`NK_T(2{\displaystyle \frac{f_{rT}}{g_T}}f^r{}_{rr}{}^{}_r\mathrm{ln}\stackrel{~}{\alpha })`$ (3f)
$`+\left(_r\beta ^r+{\displaystyle \frac{2\beta ^r}{r}}\right)f_{rT}.`$
The constraints are

$`𝒞`$ $``$ $`{\displaystyle \frac{_rf_{rT}}{g_{rr}g_T}}{\displaystyle \frac{1}{2r^2g_t}}+{\displaystyle \frac{f_{rT}}{g_{rr}g_T}}({\displaystyle \frac{2}{r}}+{\displaystyle \frac{7f_{rT}}{2g_T}}f^r{}_{rr}{}^{}){\displaystyle \frac{K_T}{g_T}}(K_r^r+{\displaystyle \frac{K_T}{2g_T}})+4\pi \rho =0,`$ (3da)
$`𝒞_r`$ $``$ $`{\displaystyle \frac{_rK_T}{g_T}}+{\displaystyle \frac{2K_T}{rg_T}}+{\displaystyle \frac{f_{rT}}{g_T}}\left(K_r^r+{\displaystyle \frac{K_T}{g_T}}\right)+4\pi J_r=0,`$ (3db)
$`𝒞_{rrr}`$ $``$ $`_rg_{rr}+{\displaystyle \frac{8g_{rr}f_{rT}}{g_T}}2f^r{}_{rr}{}^{}=0,`$ (3dc)
$`𝒞_{rT}`$ $``$ $`_rg_T+{\displaystyle \frac{2g_T}{r}}2f_{rT}=0.`$ (3dd)
Here $`𝒞`$ and $`𝒞_r`$ are the scalar and vector constraints, respectively, whereas the vanishing of $`𝒞_{rrr}`$ and $`𝒞_{rT}`$ defines the variables $`f_{rrr}`$ and $`f_{rT}`$. We are here interested only in the vacuum equations, so we set all the source terms to zero ($`T=S_{rr}=J_r=\rho =0`$), and we do not need to consider matter equations.

In essence, the evolution system is strongly hyperbolic because one can find a complete set of linearly independent combinations of the fundamental variables such that the principal part of the equations decouples. Such combinations are referred to as the characteristic variables. In a linear homogeneous system, the characteristic variables propagate their initial values exactly along the characteristic lines, and may be interpreted as waves propagating with the characteristic speed. In our case , the characteristic fields and their characteristic speeds are

$`U_r^0`$ $``$ $`g_{rr}(v_c=\beta ^r)`$ (3dea)
$`U_T^0`$ $``$ $`g_T(v_c=\beta ^r)`$ (3deb)
$`U_r^\pm `$ $``$ $`K_{rr}\pm {\displaystyle \frac{f_{rrr}}{\sqrt{g_{rr}}}}(v_c=\beta ^r\stackrel{~}{\alpha }g_T)`$ (3dec)
$`U_T^\pm `$ $``$ $`K_T\pm {\displaystyle \frac{f_{rT}}{\sqrt{g_{rr}}}}(v_c=\beta ^r\stackrel{~}{\alpha }g_T)`$ (3ded)
Here we have changed the signs of the characteristic speeds relative to  in order to agree with the sign conventions in . The characteristic fields give insights into how much data one can or must provide in addition to the initial data, depending on the region of spacetime where the solution is sought. If the region of interest has a boundary, in general one may need to prescribe additional data on this boundary in order to obtain a unique solution within that region. The additional data are the values of the characteristic fields that enter the region from outside. If such values are not prescribed, then there are many solutions for the same initial data in the region of interest. For each boundary value prescribed arbitrarily – as long as it is consistent with the initial data at the intersection of the boundary with the initial slice –, there is a unique solution. No boundary value may be freely prescribed for the characteristic fields that cross the boundary from the inside, or that run along the boundary, since such values are determined by the initial data already, in principle. In practice, since the equations are not linear nor homogeneous, the actual values of the outgoing characteristic fields at the boundary are not known in advance. Most of the time the values of some or even all the fields at the boundary, regardless of formulation type, are prescribed by means of some type of wave condition, as in  and, more elaborately, in .

In our problem, we restrict attention now to a region inside a fixed value of $`r`$. There is a set of incoming characteristic fields, which are the ones that have positive characteristic speeds in our convention, which we will be able to single out once we choose the values of the densitized lapse and shift. Suppose initial data are prescribed with vanishing values of the constraints. One can now prescribe the values of the incoming fields independently of the initial values. The question is: are such values free? It is worth investigating the possibility that the Einstein equations themselves may restrict the boundary values.

Consider the following argument. The boundary is a surface of constant radius in spacetime. The unit normal vector to the boundary is well defined. We denote it by $`e^a=(e^t,e^r,e^\theta ,e^\varphi )`$ and it can easily be computed from the gradient vector $`r,_a=(0,1,0,0)`$ by $`e^a=^4g^{ab}r,_b/\sqrt{{}_{}{}^{4}g_{}^{ab}r,_br,_a}=^4g^{ar}/\sqrt{{}_{}{}^{4}g_{}^{rr}}`$, where $`{}_{}{}^{4}g_{}^{ab}`$ is the contravariant spacetime metric. In our case, we have

$`e^t`$ $`=`$ $`{\displaystyle \frac{\beta ^r}{\stackrel{~}{\alpha }g_T\sqrt{g_{rr}(\stackrel{~}{\alpha }^2g_T^2(\beta ^r)^2)}}},`$ (3defa)
$`e^r`$ $`=`$ $`{\displaystyle \frac{\sqrt{\stackrel{~}{\alpha }^2g_T^2(\beta ^r)^2}}{\stackrel{~}{\alpha }g_T\sqrt{g_{rr}}}},`$ (3defb)
$`e^\theta `$ $`=`$ $`e^\varphi =0.`$ (3defc)
The projector $`p_{ab}`$ on the boundary surface can be found from the spacetime metric and the unit normal to the surface via

$$p_{ab}{}_{}{}^{4}g_{ab}^{}e_ae_b$$
(3defg)
with latin indices $`a,b,c\mathrm{}`$ raised and lowered with the spacetime metric. In particular, we are interested in the projector with mixed indices

$$p^a{}_{b}{}^{}=\delta ^a{}_{b}{}^{}e^ae_b$$
(3defh)
We saturate the indices of the Einstein equations $`G_{ab}=0`$ in two different ways. One equation is obtained by contracting twice with the normal to the boundary surface, $`e^a`$. We have

$$G_{ab}e^ae^b=G_{tt}(e^t)^2+2G_{tr}e^re^t+G_{rr}(e^r)^2=0$$
(3defi)
Another set of equations is obtained by contracting one idex with $`e^a`$ and the other index with $`p^b_c`$, namely: $`G_{ab}e^ap^b{}_{c}{}^{}=0`$. But because of the spherical symmetry, this contraction is identically zero for $`c=\theta ,\varphi `$. The remaining equations are

$`G_{ab}e^ap^b_t`$ $`=`$ $`G_{tt}e^t+G_{tr}e^r=0,`$ (3defj)
$`G_{ab}e^ap^b_r`$ $`=`$ $`(G_{tt}e^t+G_{tr}e^r)p^t{}_{r}{}^{}+(G_{tr}e^t+G_{rr}e^r)p^r{}_{r}{}^{}=0.`$ (3defk)
Clearly, if (3defj) and (3defi) are satisfied, then (3defk) is an identity. Thus there are only two independent equations, and we pick (3defj) and (3defi). These equations are extremely valuable for the boundary value problem because they contain no second-order derivatives with respect to $`r`$ irrespectively of the choice of $`\beta ^r`$ or $`\stackrel{~}{\alpha }`$. In the case of a first-order formulation, like ours, such equations are said to be interior to the boundary surface – not to confuse with the interior of our region! In the following, we restrict ourselves to $`\beta ^r=0`$ and $`\stackrel{~}{\alpha }=1`$ just for the sake of argument. Non trivial choices of densitized lapse and shift only make the argument more complicated without introducing any real obstacle.

Explicitly, for the case of $`\beta ^r=0`$ and $`\stackrel{~}{\alpha }=1`$, we have

$`G_{ab}e^ae_b=0`$ $`=`$ $`{\displaystyle \frac{1}{r^2g_T{}_{}{}^{2}g_{rr}^{}^2}}(2r^2g_{rr}{}_{}{}^{3/2}\dot{K}_{T}^{}3r^2g_{rr}f_{rT}{}_{}{}^{2}4rg_Tg_{rr}f_{rT}`$ (3defla)
$`+2r^2g_Tf_{rT}f_{rrr}g_{rr}{}_{}{}^{2}g_{T}^{}+r^2g_{rr}K_T{}_{}{}^{2}),`$
$`G_{ab}e^ap^b{}_{t}{}^{}=0`$ $`=`$ $`{\displaystyle \frac{2}{rg_Tg_{rr}^{3/2}}}(rg_{rr}\dot{f}_{rT}+rg_T\sqrt{g_{rr}}f_{rT}K_{rr}\sqrt{g_{rr}}K_T(rg_{rr}f_{rT}`$ (3deflb)
$`+2g_Tg_{rr}rg_Tf_{rrr})).`$
We can put these equations entirely in terms of the characteristic fields, using

$`K_T`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(U_T^++U_T^{}\right),`$ (3deflma)
$`f_{rT}`$ $`=`$ $`{\displaystyle \frac{\sqrt{g_{rr}}}{2}}\left(U_T^+U_T^{}\right),`$ (3deflmb)
$`K_{rr}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(U_r^++U_r^{}\right),`$ (3deflmc)
$`f_{rrr}`$ $`=`$ $`{\displaystyle \frac{\sqrt{g_{rr}}}{2}}\left(U_r^+U_r^{}\right).`$ (3deflmd)
and the fact that the metric components $`g_{rr}`$ and $`g_T`$ are characteristic fields themselves. For the derivatives $`\dot{K}_T`$ and $`\dot{f}_{rT}`$ we take derivatives of (3deflma) and (3deflmb) and use (3a) – written in terms of the characteristic fields– to substitute $`\dot{g}_{rr}`$ which will appear in $`\dot{f}_{rT}`$. Explicitly, for $`\beta ^r=0`$ and $`\stackrel{~}{\alpha }=1`$ we have

$`\dot{K}_T`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\dot{U}_T^++\dot{U}_T^{}\right)`$ (3deflmna)
$`\dot{f}_{rT}`$ $`=`$ $`{\displaystyle \frac{\sqrt{g_{rr}}}{2}}\left(\dot{U}_T^+\dot{U}_T^{}\right){\displaystyle \frac{g_T}{4}}\left(U_r^++U_r^{}\right)\left(U_T^+U_T^{}\right)`$ (3deflmnb)
On account of (3deflmna-3deflmnb), one can see by inspection that when both equations (3defla-3deflb) are written out entirely in terms of characteristic fields, they involve the time derivatives only of $`U_T^\pm `$, but not of $`U_r^\pm `$. Furthermore, $`\dot{U}_T^\pm `$ appear linearly in both equations. We can algebraically solve the two equations for both $`\dot{U}_T^\pm `$ as independent variables in terms of the rest. We readily find:

$`\dot{U}_T^{}`$ $`=`$ $`{\displaystyle \frac{1}{2r^2\sqrt{g_{rr}}}}(g_{rr}g_T+r^2g_{rr}\left(U_T^{}\right)^22r^2g_{rr}U_T^{}U_T^+`$ (3deflmno)
$`4r\sqrt{g_{rr}}g_TU_T^{}+r^2g_TU_T^{}U_r^+r^2g_TU_T^{}U_r^{})`$
and

$`\dot{U}_T^+`$ $`=`$ $`{\displaystyle \frac{1}{2r^2\sqrt{g_{rr}}}}(g_{rr}g_T+r^2g_{rr}\left(U_T^+\right)^22r^2g_{rr}U_T^+U_T^{}`$ (3deflmnp)
$`+4r\sqrt{g_{rr}}g_TU_T^+r^2g_TU_T^+U_r^++r^2g_TU_T^+U_r^{})`$
Both equations (3deflmno) and (3deflmnp) are necessary because they are linearly independent. The question is: how are we to interpret them within the framework of the initial-boundary value problem?

To start with, the fields $`g_{rr}`$ and $`g_T`$ may be considered as known sources, because they travel upwards along the boundary. If the problem was linear and homogeneous, then the values of $`g_{rr}`$ and $`g_T`$ at any point on the boundary would be exactly their initial values. Because the problem is non-linear, the values at any point on the boundary must be integrated, but can still be regarded as sources. A similar thinking may be used for the outgoing characteristic field $`U_r^+`$, with the complication that $`U_r^+`$ travels towards the boundary but not directly upwards, and thus the computation of its values at the boundary will be more involved, but conceptually no different.

Additionally, we may take the point of view that the incoming characteristic field $`U_r^{}`$ is arbitrary, since there are only two equations to satisfy at the boundary and $`\dot{U}_r^{}`$ does not appear in any of them. From this perspective, the value of this field at the boundary is truly free, and must be considered as one true degree of freedom of the boundary value problem for the Einstein equations, in addition to the degrees of freedom contained in the initial data.

Continuing within such an interpretation, Eq. (3deflmno) provides exact boundary values for the incoming characteristic field $`U_T^{}`$ in terms of the free incoming field $`U_r^{}`$ and the fields that either run along the boundary or cross from the inside, all of which are determined by the initial data in principle, as argued before.

There remains Eq. (3deflmnp). From our point of view, this equation predicts boundary values for $`U_T^+`$. There is no doubt that it must hold, in order for the solution that we are seeking to satisfy all the Einstein equations at all points on the boundary surface. But $`U_T^+`$ is, in principle, determined by the initial data by propagation along characteristics, as explained before. There seem to be two possibilities here: either the outgoing field $`U_T^+`$ propagated from the initial data satisfies Eq. (3deflmnp), or it does not. If the boundary values of $`U_T^+`$ do not satisfy Eq. (3deflmnp), then the initial-boundary value problem is inconsistent and can not yield a solution to the Einstein equations in the region of interest.

Alternatively, if the boundary values of $`U_T^+`$ propagated along characteristics do satisfy Eq. (3deflmnp), then the initial-boundary value problem is consistent. In this case, Eq. (3deflmnp) can be used to prescribe the values of $`U_T^+`$ with more accuracy and reliability than the propagation along characteristics discussed above because it is an ordinary differential equation, and should be preferred.

In the following Section we show that the initial-boundary value problem is consistent and we should be using Eq. (3deflmnp) to prescribe the boundary values of $`U_T^+`$ instead of propagating them by characteristics from the initial data, because the propagation by characteristics is inaccurate in itself –even more so in the case that the characteristic speeds depend on the fields themselves and must be calculated at every step. One can anticipate that this is the case based on some intuition. Why should the boundary values of the outgoing fields be “constrained” by an equation such as Eq. (3deflmnp)? Because, in principle, the boundary values of the outgoing fields must reflect in some way the constraints on the intial data that determine them. The initial data are related; their inter-relationships must be propagated along characteristics. Eq. (3deflmnp) may well be a shortcut out of constraint propagation.

## 3 Consistency with constraint propagation

If thought of as functions which could take any real values, the initial constraints $`𝒞,𝒞_r,𝒞_{rrr},𝒞_{rT}`$ evolve in time according to another strongly hyperbolic evolution system, with the following characteristic fields and speeds:

$`C_1`$ $`=`$ $`𝒞+{\displaystyle \frac{𝒞_r}{\sqrt{g_{rr}}}}(v_1^c=\beta ^r\stackrel{~}{\alpha }g_T)`$ (3deflmnqa)
$`C_2`$ $`=`$ $`𝒞{\displaystyle \frac{𝒞_r}{\sqrt{g_{rr}}}}(v_2^c=\beta ^r+\stackrel{~}{\alpha }g_T)`$ (3deflmnqb)
$`C_3`$ $`=`$ $`𝒞_{rrr}(v_3^c=\beta ^r)`$ (3deflmnqc)
$`C_4`$ $`=`$ $`𝒞_{rT}(v_4^c=\beta ^r)`$ (3deflmnqd)
In the case of $`\beta ^r=0`$ and $`\stackrel{~}{\alpha }=1`$, $`C_2`$ is the only incoming characteristic field. If one wishes for all four functions to take the value zero in the region interior to some fixed radius, then three of them will be zero by choice of initial data, but $`C_2`$ must be set to zero at the boundary in order to propagate inwardly from there. But one simply may set the value of $`C_2`$ to zero because one does not evolve the constraint propagation system. It is one thing to impose the value $`0`$ on the function $`C_2`$, but a different thing to expect $`C_2`$ to be vanishing at the boundary as a function of the fundamental variables of evolution. In particular, by inspection one can readily see that $`C_2`$ contains the derivative of $`U_T^{}`$ with respect to $`r`$. Therefore, although it is natural and desirable to have $`C_2=0`$ on the boundary, it is an impossible task to have $`C_2=0`$ as a boundary condition for the fundamental fields of evolution. $`C_2`$ will have to vanish for the evolution to produce a solution of the Einstein equations, and thus the vanishing of $`C_2`$ and all the other constraints must be checked for consistency at the boundary, after the solution has been found and the $`r`$derivatives can be evaluated; but it is not practical as a boundary condition.

A strategy that has been proposed  to circumvent this obstacle is to use the evolution equations to turn the $`r`$derivatives of the characteristic fields that appear in the expression of $`C_2`$ into $`t`$derivatives. The procedure must work in this symmetry-reduced case because, by construction, all characteristic fields have time derivatives proportional to their $`r`$derivatives, up to terms of zeroth order. Therefore, $`C_2=0`$ can be turned into an evolution equation for $`U_T^{}`$ that will be restricted to the boundary. Following , we do that next, maintaining the restrictions $`\beta ^r=0`$ and $`\stackrel{~}{\alpha }=1`$. We have explicitly

$`C_2`$ $`=`$ $`{\displaystyle \frac{1}{2r^2g_T{}_{}{}^{2}g_{rr}^{}^{5/2}}}(2r^2g_Tg_{rr}_rK_T2r^2g_Tg_{rr}{}_{}{}^{3/2}_{r}^{}f_{rT}+g_Tg_{rr}^{5/2}`$ (3deflmnqr)
$`4rg_Tg_{rr}{}_{}{}^{3/2}f_{rT}^{}7r^2g_{rr}{}_{}{}^{3/2}f_{rT}^{}+2r^2g_T\sqrt{g_{rr}}f_{rT}f_{rrr}+2r^2g_Tg_{rr}{}_{}{}^{3/2}K_{T}^{}K_{rr}`$
$`+r^2g_{rr}{}_{}{}^{5/2}K_{T}^{}{}_{}{}^{2}+4rg_tg_{rr}^2K_T2r^2g_Tg_{rr}f_{rT}K_{rr}2r^2g_{rr}^2f_{rT}K_T)`$
Solving for $`_rf_{rT}`$ from the evolution equation (3d), and for $`_rK_T`$ from the evolution equation (3f), substituting into (3deflmnqr), and replacing all appearances of $`f_{rT},f_{rrr},K_T`$ and $`K_{rr}`$ in terms of the characteristic fields using (3deflma-3deflmd), Eq. (3deflmnqr) turns into a new expression, which is not an initial constraint anymore, so we use quotes to make this point explicit:

$`\mathrm{`}\mathrm{`}C_2\mathrm{"}`$ $`=`$ $`{\displaystyle \frac{1}{2r^2g_T{}_{}{}^{2}g_{rr}^{}^{5/2}}}(2r^2g_{rr}{}_{}{}^{2}\dot{U}_{T}^{}r^2g_{rr}{}_{}{}^{5/2}(U_T^{})_{}^{2}r^2g_Tg_{rr}{}_{}{}^{3/2}U_{T}^{}U_T^+`$ (3deflmnqs)
$`+r^2g_Tg_{rr}{}_{}{}^{3/2}U_{T}^{}U_r^{}g_Tg_{rr}{}_{}{}^{5/2}+2r^2g_{rr}{}_{}{}^{5/2}U_{T}^{+}U_T^{}+4rg_Tg_{rr}{}_{}{}^{2}U_{T}^{})`$
Setting $`\mathrm{`}\mathrm{`}C_2\mathrm{"}=0`$ yields a boundary condition for $`U_T^{}`$. The point is that $`C_2=0`$ on the boundary is not the same as $`\mathrm{`}\mathrm{`}C_2\mathrm{"}=0`$, but they are the same along the boundary if evaluated on fundamental fields that satisfy the evolution equations. In other words: they differ by a linear combination of the evolution equations. But more interestingly, $`\mathrm{`}\mathrm{`}C_2\mathrm{"}=0`$ is exactly the same equation as (3deflmno), as can be verified by inspection. Therefore the method of “trading” space derivatives for time derivatives advocated in  is equivalent to solving one of the Einstein equations that are interior to the boundary, there being a two-dimensional set of such equations. This could be expected because the projections of the Einstein equations along the direction normal to a surface of fixed radius are the linear combinations that have no second-order $`r`$derivatives, and clearly the “trading” is equivalent to linearly combining the Einstein equations. Now, because Eq. (3deflmnp) is equivalent to $`C_2=0`$ on the boundary, it guarantees that $`C_2`$ will vanish in the interior by propagation along characteristics. So Eq. (3deflmnp) is not only consistent with the boundary value problem, but is necessary as well.

We can also use the “trading” of space derivatives for time derivatives to write an equation $`\mathrm{`}\mathrm{`}C_1\mathrm{"}=0`$ from $`C_1=0`$ along the boundary. This is exactly Eq. (3deflmnp), the meaning of which is explained in the previous section, and we conclude that Eq. (3deflmnp) and $`C_1`$ are equivalent along the boundary if the evolution equations are satisfied. But how much of this argument proves that the constraint has been “propagated” and where does the propagation start? Clearly $`C_1`$ propagates out from the initial data towards the boundary, therefore $`C_1`$ will vanish along the boundary as a consequence of its vanishing initially. Since the evolution will be satisfied by construction, then $`\mathrm{`}\mathrm{`}C_1\mathrm{"}`$ will vanish as a consequence of $`C_1`$ vanishing initially as well. Seemingly, thus, the boundary values of the fields satisfying Eq. (3deflmnp) must be consistent with the initial values satisfying $`C_1=0`$. For this reason, and because Eq. (3deflmnp) is an ordinary differential equation, it is not only consistent but advisable to use Eq. (3deflmnp) to obtain the values of $`U_T^+`$ instead of propagating them from the initial values along characteristics. At any rate, Eq. (3deflmnp) should be checked for consistency regardless of method.

Two objections may come to mind on a surface glance to Eq. (3deflmnp), which we want to anticipate. It might conceivably be argued that if the values of $`U_T^+`$ at the boundary are obtained from the initial values by propagation along characteristics, then Eq. (3deflmnp) could also be thought of as an equation for the incoming field $`U_r^{}`$ with $`U_T^+`$ given, instead of prescribing $`U_T^+`$ with a free $`U_r^{}`$. But this is a flawed line of thought, for the boundary values of $`U_T^+`$ propagated along characteristics are consistent with Eq. (3deflmnp), and substituting $`U_T^+`$ into Eq. (3deflmnp) must yield an identity, leaving no equation to solve for $`U_r^{}`$, which would thus remain arbitrary.

Secondly, how is it that Eq. (3deflmnp) involving a free boundary field $`U_r^{}`$ can be consistent with propagation of $`U_T^+`$ along characteristics, given that $`U_r^{}`$ is a field which the initial data know nothing about? The initial data do not know about $`U_r^{}`$, but the characteristics do, because the problem is non linear. Thus the propagation of $`U_T^+`$ along characteristics will necessarily yield boundary values of $`U_T^+`$ that know about earlier boundary values of incoming fields. It is fine for Eq. (3deflmnp) to involve $`U_r^{}`$, as long as it does not involve $`\dot{U}_r^{}`$.

For additional support of our scheme to use Eq. (3deflmno) and Eq. (3deflmnp) to prescribe $`U_T^{}`$ and $`U_T^+`$, respectively, assuming that $`U_r^{}`$ is given freely, we can linearize the equations around flat space and verify that the scheme is consistent, since the boundary value problem of linear hyperbolic equations is clear and enlightening. We have

$`g_{rr}`$ $`=`$ $`1+\widehat{g}_{rr}`$ (3deflmnqta)
$`g_T`$ $`=`$ $`1+\widehat{g}_T`$ (3deflmnqtb)
$`K_{rr}`$ $`=`$ $`\widehat{K}_{rr}`$ (3deflmnqtc)
$`K_T`$ $`=`$ $`\widehat{K}_T`$ (3deflmnqtd)
$`f_{rrr}`$ $`=`$ $`{\displaystyle \frac{4}{r}}+\widehat{f}_{rrr}`$ (3deflmnqte)
$`f_{rT}`$ $`=`$ $`{\displaystyle \frac{1}{r}}+\widehat{f}_{rT}`$ (3deflmnqtf)
where all hatted quantities are small. We will keep only linear terms in such quantities. The evolution equations become

$`\dot{\widehat{g}}_{rr}`$ $`=`$ $`2\widehat{K}_{rr}`$ (3deflmnqtua)
$`\dot{\widehat{g}}_T`$ $`=`$ $`2\widehat{K}_T`$ (3deflmnqtub)
$`\dot{\widehat{K}}_{rr}+_r\widehat{f}_{rrr}`$ $`=`$ $`{\displaystyle \frac{42}{r^2}}\widehat{g}_{rr}+{\displaystyle \frac{48}{r^2}}\widehat{g}_T+{\displaystyle \frac{10}{r}}\widehat{f}_{rrr}{\displaystyle \frac{44}{r}}\widehat{f}_{rT},`$ (3deflmnqtuc)
$`\dot{\widehat{K}}_T+_r\widehat{f}_{rT}`$ $`=`$ $`{\displaystyle \frac{1}{r^2}}\widehat{g}_{rr}+{\displaystyle \frac{2}{r^2}}\widehat{g}_T{\displaystyle \frac{4}{r}}\widehat{f}_{rT}`$ (3deflmnqtud)
$`\dot{\widehat{f}}_{rrr}+_r\widehat{K}_{rr}`$ $`=`$ $`{\displaystyle \frac{12}{r}}\widehat{K}_{rr}+{\displaystyle \frac{4}{r}}\widehat{K}_T`$ (3deflmnqtue)
$`\dot{\widehat{f}}_{rT}+_r\widehat{K}_T`$ $`=`$ $`{\displaystyle \frac{2}{r}}\widehat{K}_T`$ (3deflmnqtuf)
Thus the characteristic fields are $`\widehat{g}_{rr},\widehat{g}_T`$ and

$`\widehat{U}_T^\pm `$ $``$ $`\widehat{K}_T\pm \widehat{f}_{rT},`$ (3deflmnqtuva)
$`\widehat{U}_r^\pm `$ $``$ $`\widehat{K}_{rr}\pm \widehat{f}_{rrr}.`$ (3deflmnqtuvb)
Our boundary conditions, Eqs. (3deflmno) and (3deflmnp), linearize to

$$\dot{\widehat{U}}{}_{T}{}^{}=\frac{9}{2r^2}\widehat{g}_{rr}\frac{3}{2r^2}\widehat{g}_T+\frac{1}{r}\widehat{U}_r^{}+\frac{1}{r}\widehat{U}_T^+,$$
(3deflmnqtuvw)
and

$$\dot{\widehat{U}}{}_{T}{}^{+}=\frac{9}{2r^2}\widehat{g}_{rr}\frac{3}{2r^2}\widehat{g}_T\frac{1}{r}\widehat{U}_r^+\frac{1}{r}\widehat{U}_T^{},$$
(3deflmnqtuvx)
respectively. Thus (3deflmnqtuvw) can be used to prescribe values for $`\widehat{U}_T^{}`$ if $`\widehat{U}_r^{}`$ is given arbitrarily, whereas (3deflmnqtuvx) can be used to calculate values for $`\widehat{U}_T^+`$ irrespective of $`\widehat{U}_r^{}`$, as anticipated.

## 4 Boundary conditions for generic three-dimensional strongly hyperbolic formulations

How much of the argument in spherical symmetry actually depends on the symmetry assumption? Not a great deal. Suppose we have a strongly hyperbolic formulation of the Einstein equations in terms of 6+6+18 variables representing the three-metric and all its first derivatives. The evolution equations then look like

$$\dot{u}=A^i_iu+b$$
(3deflmnqtuvy)
where $`u`$ is the 30-dimensional vector of all the fundamental variables, and there are 4+18 constraints on the initial data:

$`𝒞`$ $`=`$ $`0`$ (3deflmnqtuvza)
$`𝒞_i`$ $`=`$ $`0`$ (3deflmnqtuvzb)
$`𝒞_{ijk}`$ $`=`$ $`0`$ (3deflmnqtuvzc)
where $`𝒞`$ and $`𝒞_i`$ are the scalar and vector constraint, respectively, and $`𝒞_{ijk}`$ are the constraints necessary to reduce the equations from second to first order in space (they define the 18 first order variables). Wherever the set of evolution equations and the constraints are satisfied, the ten Einstein equations $`G_{ab}=0`$ for the ten components of the spacetime metric $`g_{ab}`$ are satisfied, equivalently.

Suppose the evolution equations (3deflmnqtuvy) are strongly hyperbolic, and that they imply, in the usual manner , a second (or subsidiary) system of equations for the constraint functions $`𝒞,𝒞_i`$ and $`𝒞_{ijk}`$. Suppose that this second system of equations, for the constraints, is also strongly hyperbolic. Assume we have identified the characteristic variables of both strongly hyperbolic systems.

Suppose now that we seek a solution in the region interior to some fixed value of the coordinate $`x^1`$ (any spacelike coordinate). The unit normal to the boundary is then

$$e^a=\frac{g^{ab}\delta _b^1}{\sqrt{g^{ab}\delta _a^1\delta _b^1}}=\frac{g^{1a}}{\sqrt{g^{11}}}$$
(3deflmnqtuvzaa)
and the projector on the boundary surface is $`p^a{}_{b}{}^{}=\delta _b^a+e^ae_b`$ with $`e_b=\delta _b^1/\sqrt{g^{11}}`$. We can write down $`G_{ab}e^ae^b=0`$ and $`G_abe^ap^b{}_{c}{}^{}=0`$. In general these will be four equations with no second derivatives with respect to $`x^1`$, as can be proven by direct calculation. Contracting the Einstein tensor with $`e^a`$ yields a vector equation with two pieces:

$$G_{ab}e^a=R_{ab}e^a\frac{\delta _b^1}{2\sqrt{g^{11}}}R.$$
(3deflmnqtuvzab)
We next show that the components $`b1`$ of $`R_{ab}e^a`$ have no second derivatives of the metric with respect to $`x^1`$, and that the Ricci scalar term cancels out the second $`x^1`$derivatives that appear in the first term for $`b=1`$. We use the following well-known expression of the Ricci tensor  in which the second-order derivatives of the metric are explicit:

$$R_{ab}=\frac{1}{2}g^{cd}\left(g_{cb,ad}+g_{ad,cb}g_{cd,ab}g_{ab,cd}\right)+g^{cd}\left(\mathrm{\Gamma }_{ad}^e\mathrm{\Gamma }_{edb}\mathrm{\Gamma }_{ab}^e\mathrm{\Gamma }_{ecd}\right).$$
(3deflmnqtuvzac)
Contracting with $`e^a`$ yields

$$R_{ab}e^a=\frac{1}{2\sqrt{g^{11}}}\left(g^{1a}g^{1c}g_{ac,1b}g^{11}g^{cd}g_{cd,1b}\right)+\mathrm{}$$
(3deflmnqtuvzad)
where $`\mathrm{}`$ represents terms that have no second derivatives with respect to $`x^1`$. Thus $`R_{ab}e^a`$ has no second $`x^1`$derivatives except for $`b=1`$. The second part of the Einstein tensor is nonvanishing only for $`b=1`$. Therefore the only component of $`G_{ab}e^a`$ that might contain a second $`x^1`$derivative is $`b=1`$. Now $`R`$ is exactly

$$R=g^{1a}g^{1c}g_{ac,11}g^{11}g^{cd}g_{cd,11}+\mathrm{}$$
(3deflmnqtuvzae)
Therefore, the terms with second $`x^1`$derivatives of the Ricci scalar part of the Einstein tensor cancel exactly with those in the Ricci tensor part, with the consequence that $`G_{ab}e^a`$ has no second $`x^1`$derivatives of the metric for any value of $`b`$. Contracting with $`e^b`$ or with $`p^b_c`$ will not add new second $`x^1`$derivatives, but provides perhaps a convenient splitting of the four equations. The four equations

$$_bG_{ab}e^a=0$$
(3deflmnqtuvzaf)
or equivalently

$$G_{ab}e^ae^b=0\text{ and }G_{ab}e^ap^b{}_{c}{}^{}=0$$
(3deflmnqtuvzag)
must be satisfied by any solution to the Einstein equations and are interior to the boundary surface, which makes them ideal boundary conditions. There are four equations and thirty variables, but of the thirty variables many are outgoing or run along the boundary. In fact, one can always write the evolution system so that the metric components propagate with vanishing characteristic speeds, regardless of the choice of lapse and shift. So we have at most $`(306)/2`$ relevant incoming characteristic fields. Clearly, some of the incoming fields will be left arbitrary at the boundary – perhaps more than 8, since in our example in spherical symmetry we find that one of the boundary conditions applies to an outgoing field instead of an incoming one. The particulars of which characteristic fields will be affected by the boundary equations will depend entirely on the details of the first-order strongly hyperbolic formulation of choice. As an illustration, the direct generalization of our argument to the case of the Einstein-Christoffel formulation in three-dimensions without the restriction of spherical symmetry will be reported elsewhere.

Conceptually, the scheme runs as follows. Choose initial data that satisfy (3deflmnqtuvz). Choose boundary data that satisfy (3deflmnqtuvzaf). The vanishing of the initial constraints guarantees that the outgoing and static constraints will remain vanishing in the region of interest. The boundary equations (3deflmnqtuvzaf) are equivalent to the vanishing of the incoming constraints at the boundary on account of their being related to the incoming constraints by linear combinations with the Einstein equations that are satisfied at the boundary – the evolution is satisfied by construction and the outgoing constraints are satisfied by propagation along characteristics. Therefore the boundary equations (3deflmnqtuvzaf) guarantee that the incoming constraints are vanishing at the boundary, which in turn guarantees that they will remain vanishing in the region of interest.

Although the preceding calculation assumes that the normal to the surface of fixed value of $`x^1`$ is spacelike, clearly Eqs. (3deflmnqtuvzaf) will be interior to the surface of fixed value of $`x^1`$ even if the normal is timelike, in which case the proof that the equations contain no second $`x^1`$derivatives runs the same if we substitute $`\sqrt{g^{11}}`$ with $`\sqrt{|g^{11}|}`$. This remark applies to interior boundaries that lie within the event horizon of a black hole spacetime, and may be relevant to numerical simulations of binary black holes.

On the other hand, in this Section we have restricted the discussion to three-dimensional boundaries taken one at a time. As is the case with any three-dimensional boundary-value problem, the intersection of two boundaries requires special attention. In this respect, consistency issues between the boundary equations that may arise at the intersection of two boundary surfaces remain to be studied, and their resolution may depend on the particulars of the formulation of the initial value problem at hand. The interested reader is referred to  where a particular representation of the projection of the Einstein equations on the boundary is implemented and corners and edges are treated, for the case of an ADM-like formulation of the initial value problem.

## 5 Concluding remarks and outlook

Even though we have, for the most part, developed our arguments explicitly in the case of spherical symmetry for the EC formulation, the fundamental relevance of the argument to the initial-boundary value problem of the Einstein equations does not depend on the symmetry restrictions nor on the particulars of the hyperbolic formulation, as we show in Section 4.

We argue that given any boundary at a fixed value of a coordinate of an initial-boundary value problem for the Einstein equations, the vanishing of the components of the projection of the Einstein tensor $`G_{ab}`$ along the normal $`e^a`$ to the boundary, namely $`G_{ab}e^b_a=0`$, constitute necessary and consistent boundary conditions, for essential reasons, as follows.

First, the components of the projection of the Einstein tensor along the normal to such a boundary, $`_a`$, contain no second-order derivatives sticking out of the boundary. If the initial value problem is stated in first order form, such as any of the hyperbolic formulations available, their vanishing becomes differential equations for the boundary values of the fundamental fields excluding their derivatives across the boundary. Moreover, if the initial value problem is stated in second-order form in the space coordinates, such as ADM form  or conformal form , they constitute a type of mixed Neumann-Dirichlet conditions, the appropriateness of which is worth investigating in full depth. In this regard, some representations of $`_a=0`$ as boundary conditions have been used along with an ADM type formulation of the linearized Einstein equations in  to investigate the consistency of pure Neumann or Dirichlet conditions with the evolution.

Second, in the case of strongly hyperbolic formulations that propagate the constraints in a strongly hyperbolic fashion, the vanishing of all four boundary equations $`_a`$ guarantees the vanishing of the incoming constraints in the region of interest and is consistent with the propagation of the outgoing constraints along characteristics. In practice, this means that one can use $`_a=0`$ to prescribe as many boundary values as possible, without worrying about inconsistencies with the initial values. This necessary analytic consistency should not be mistaken for numerical consistency, however. In fact, it is not known at this stage whether this analytic consistency is stable under small perturbations of either the initial data or the boundary data, a point that is critical to the numerical implementation. Additionally, how to proceed in order to investigate the stability of an initial-boundary-value problem with constraints remains, to our knowledge, an open problem not addressed in the standard reference literature of strongly hyperbolic systems such as .

In the case of formulations that do not guarantee stable constraint propagation, $`_a=0`$ on the boundary still appeals to us as the next best thing. In fact, in such cases there are hardly any handles on boundary conditions. Notice that imposing $`_a=0`$ on the boundary does not involve any manipulation of the choice of evolution equations; it does not affect the evolution equations themselves nor the system of propagation of the constraints. The converse is not true: any manipulation of the evolution equations (by adding linear combinations of the constraints in the usual manner) necessarily affects both the propagation of the constraints and the role of $`_a=0`$ on the boundary. The interested reader is referred to  for a series of numerical studies of the effect of the manipulation of the evolution equations on constraint propagation exclusively, without imposing $`_a=0`$.

Several issues remain open at this time. First and foremost, there remains the issue of whether the boundary equations $`_a=0`$ are consistent with a well-posed initial-boundary-value problem in the analytic sense, that is: in the sense that the solutions at a later time are continuous functions of the initial data and the boundary free data. The answer to this question will vary among the different strongly hyperbolic formulations available. The answer may be relevant to numerical relativity because, even though well-posedness of the initial-boundary value problem does not guarantee a stable numerical implementation, it is considered as a necessary  or at least desirable feature of the continuous equations being implemented. Second, any issues that may arise exclusively in connection with the numerical implementation of $`_a=0`$ –but that are otherwise irrelevant to the analytic initial-boundary value problem of the Einstein equations– are also open at this time and are worth pursuing, their details being strongly dependent on the particulars of the formulation at hand. From the numerical point of view, for instance, the presence of undifferentiated terms, –which is not relevant to our current argument and does not affect the well-posedness of first-order problems– is critical to numerical stability . In this respect, the effect of the undifferentiated terms of the boundary equations $`_a=0`$ on any issues of stability remains to be determined, being, almost certainly, critically dependent on the particulars of the formulation at hand. Additional open issues at this time include checking for the consistency of $`_a=0`$ with other boundary conditions already in use in numerical simulations, and, eventually, how the use of $`_a=0`$ on the boundary may affect the run time of numerical simulations.

This work was supported by the NSF under grants No. PHY-0070624 to Duquesne University and No. PHY-0135390 to Carnegie Mellon University.

## References