format: text

text:
# Exploring the Role of Jets in the Radio/X-ray Correlations of GX 339-4

## 1 Introduction

The Galactic black hole candidate (BHC) X-ray binary (XRB) GX 339$``$4 displays a wide variety of spectral states, and has recently resurfaced at a high flux level after being in an extended “off” state since 1999 (see, e.g., Kong et al. 2000). It has been extensively, and often simultaneously, observed across broad energy bands: from the radio (Hannikainen et al. 1998; Fender et al. 1999; Corbel et al. 2000—hereafter C00; Corbel et al. 2002, to be submitted—hereafter C02), infrared (IR) and optical (Soria et al. 1999; Shahbaz et al. 2001; Corbel & Fender 2002—hereafter CF02), to the X-rays (Wilms et al. 1999—hereafter W99; Nowak et al. 2002—hereafter NWD). These data sets are relatively untainted by contamination from the companion star. This, and the fact that GX 339$``$4 often transitions between the canonical XRB states, makes it an ideal source for comparing the emission processes involved in these various states.

GX 339$``$4 spends the great majority of its time in the low/hard state (LHS) (for a review of BHC states, see Nowak 1995), at flux levels varying over roughly three orders of magnitude in the X-ray waveband (NWD; C02). Its LHS spectra are fairly “typical” in that they show a flat-to-inverted radio component, a weak thermal contribution, and a hard power-law in the soft to hard X-rays (Zdziarski et al. 1998; Fender 2001; Revnivtsev et al. 2001; NWD). Although no jets have yet been resolved, the radio emission shows 2% linear polarization, and the brightness temperature (defined as $`S_\nu =\mathrm{\Omega }B_\nu (T_\mathrm{B})`$ where $`\mathrm{\Omega }`$ is the solid angle of the emitting region) requires a spatial extension greater than the binary separation (W99; C00). At the time of writing, GX 339-4 briefly entered the very high state (VHS) and made a transient-like radio outburst which was five times brighter than the brightest radio emission ever observed in the LHS. The optically thin spectrum had a spectral index $`\alpha 0.6`$ ($`F_\nu \nu ^\alpha `$) (R. Fender, S. Corbel, private comm.), which is characteristic of a transient bright jet. This makes GX 339$``$4 the first source to show its radio properties clearly in the LHS, the high/soft state (HSS) and transient/VHS. Combined with the fact that jets have already been imaged in the other three persistent Galactic BHC which share similar LHS spectra (e.g., 1E1740.7-2942, Mirabel et al. 1992; GRS 1758-258, Rodríguez et al. 1992; and Cygnus X-1, Stirling et al. 2001), the case for a jet in this source is now very strong.

A special feature of GX 339$``$4 is that its radio emission has been observed to correlate tightly with the X-ray emission over the entire range of LHS luminosities, and also down to a very weak level of emission in the “off” state (Hannikainen et al. 1998; C00; C02). Similarly, during changeover to the HSS, the seemingly linked radio and hard (20-100 keV) X-ray emission show an anti-correlation with the assumed thermal disk emission in the softer (2-12 keV) X-ray band (Fender et al. 1999; C00). This suggests not only that the plasma leading to the synchrotron radiation plays some role at higher energy, but also that the jet is either quenched (Fender et al., 1999) or otherwise unable to form when the thermal disk component dominates (Meier, 2001; Merloni & Fabian, 2002).

On the other hand, the LHS X-ray spectra from GX 339$``$4 show additional features which have been interpreted as signs of reflection off colder disk material, and can be fit by one of several disk-corona models as discussed in detail in NWD. In this paper, the authors find correlations between the reflection covering fraction and both the soft X-ray flux (see also Zdziarski et al. 1999; Revnivtsev et al. 2001), and the time lags between the soft and hard X-ray variability. Similar to what Pottschmidt et al. (2000) report for Cyg X-1—which also shows radio/X-ray correlations (Brocksopp et al., 1999)—NWD find that these time lags anti-correlate with the coronal compactness, and peak shortly after the transition from the HSS to the LHS, when the assumed jet emission begins to rise. They suggest an interpretation that the time lags may be associated with the extended synchrotron-emitting plasma.

While the model fits presented in NWD can well explain the X-ray features, and suggest a link between the corona and jet, they do not explore the nature of this link or how this emission is related to that of the radio. One model which can explain the tight radio/X-ray correlations proposes that the jets directly contribute to the LHS X-ray spectrum via synchrotron emission (Markoff et al. 2001—hereafter MFF), as occurs in several AGN (e.g., 3C 66B, Hardcastle et al. 2001 and M87, Marshall et al. 2002; Wilson & Yang 2002). However, this model has only been applied until now to one source, XTE J1118$`+`$480, which may be unusual in that it has extremely weak or absent reflection features (Miller et al., 2001).

In this paper, we discuss the application of a jet model to the collected simultaneous or quasi-simultaneous data sets for GX 339$``$4 in order to investigate whether it succeeds in explaining the observed radio/X-ray correlations. We discuss how the model fares against the broadband data and compare the model parameters for the various data sets to each other, as well as to XTE J1118$`+`$480.

## 2 Data

The 1981 data set was compiled by CF02, where it is explained in further detail. The X-rays come from the Ariel 6 instrument (Ricketts, 1983), with simultaneous optical/IR from Pedersen (1981) and Motch et al. (1981).

All radio observations were performed with the Australia Telescope Compact Array (ATCA) and, with the exception of the 2000 September 15 observation (to appear in C02), have been presented previously by C00. All observations from C00 in which the radio flux density was weak (i.e., $`1`$ mJy) were reanalyzed for this paper.

The X-ray data presented here have been discussed previously by NWD and C02. All data were obtained by the Rossi X-ray Timing Explorer (RXTE), with the exception of the 2000 September 15 data, which were obtained by BeppoSAX.

In this work we are attempting to describe the *broadband features* of the jet model, and to simultaneously explain radio, IR, optical, and X-ray data (taken with a wide variety of instruments). Furthermore, we are not yet attempting to explain detailed features within the X-ray, e.g., Fe lines or reflection features. Therefore, rather than forward fold the jet models through the detector response matrices, we instead compare the models to “unfolded” count flux rates. This is not strictly fitting, as we cannot obtain chi-squared estimates from unfolded data in this way, but for simplicity we will use the word “fit” to describe the comparison to the data throughout the paper.

To unfold the X-ray data, the RXTE data extraction, data binning, response matrix generation, and background subtraction were carried out in an identical manner as presented by NWD. We then fit models comprised of: a multi-temperature disk blackbody (e.g., Mitsuda et al. 1984; Makishima et al. 1986) with peak temperature fixed at 0.25 keV, an exponentially cutoff broken power-law with break energy at $`10`$ keV, and a (potentially) broad gaussian line with energy fixed to 6.4 keV. (This model is essentially identical to that discussed by W99). These model fits were then unfolded through the detector response matrix and multiplied by the ratio of the count rate data to the fitted model folded through the detector response. For the RXTE observations, data was considered only for the 3-200 keV range or less for the faintest data sets. Furthermore, HEXTE data was normalized to the PCA data (e.g., see W99).

As discussed by C02, simultaneous RXTE/ASCA observations indicate a probable faint background source that contaminates the RXTE spectrum— but not the narrow field of view ASCA spectrum—at *extremely* low flux levels. The RXTE observation of 1999 July 29 is assumed to be heavily dominated by this contaminating source, and this spectrum, multiplied by 0.78, is subtracted, before model fitting and data unfolding, from all RXTE observations occurring later than the observation of 1999 May 14. This normalization was determined by comparison of spectra obtained simultaneously by RXTE and ASCA on two separate occasions (the required normalization constants for these two observations were 0.73 and 0.83; see C02). Note that the 2000 September 15 BeppoSAX observation is most likely *not* contaminated by this background source, which could easily be half a degree away and still affect the RXTE observations, as BeppoSAX has a narrower field of view than RXTE. A more detailed discussion of this background source appears in C02.

The dates of the observations are listed in Table 1. Data sets labeled 20181 were first discussed by W99. Data sets 40108-02 were first discussed by NWD. Data sets 40108-03, first appearing in C02, were part of the same series of observations originally discussed by NWD; however, these latter observations were too faint for simultaneous spectral and X-ray variability analyses. Observations 40105-02-02, 40104-01-08 and 21136001 were part of a series of observations to study GX 339$``$4 in quiescence, and are discussed by C02.

## 3 Model

The details of the jet model can be found in MFF, and references therein, and we give a brief summary here. We consider that an accretion disk may contribute to the near infrared (NIR), optical and possibly soft X-ray emission, but that the radio through at least far infrared (FIR) are due to synchrotron emission from the jets, which could possibly extend through the NIR to the X-rays. The inverse Compton (IC) upscattering by the hot jet electrons of both the thermal disk photons, as well as the synchrotron photons, is also included.

Our disk model is relatively simplistic, since we are currently exploring its limits in terms of contribution and there are already many detailed disk models in existence (see Poutanen 1998; NWD for reviews). We here consider one commonly invoked model in which a standard thin, optically thick disk (SD) (Shakura & Sunyaev, 1973) exists down to some transition radius $`r_{\mathrm{tr}}`$, at which point the flow becomes hotter, optically thin and relatively non-radiative (e.g., Esin et al. 1997). This type of non-radiative inner disk has also been invoked for low-luminosity AGN such as NGC 4258 (Yuan et al., 2002). The thermal bump—seen in several of the GX 339$``$4 data sets in the optical band—is assumed to originate in the disk, which also provides seed photons for jet IC as mentioned above. We assume for now that the jet is oriented 90 from the plane of the disk, but in the future we will consider other orientations (see, e.g., Maccarone 2002).

In their study of the BHC XTE J1118+480, MFF interpreted the extreme-ultraviolet (EUV) upper limits (Hynes et al., 2000) as the hottest, and thus innermost, emission from this outer disk. Therefore, they could place limits on its maximum temperature and luminosity. They derived a location for the transition radius of $`r_{\mathrm{tr}}10^210^3r_g`$ which is consistent with other models (e.g., Liu et al. 1999; Esin et al. 2001). However, no EUV data exist for GX 339$``$4 , and so we constrain $`r_{\mathrm{tr}}`$ with the low-frequency X-ray data. We have also made the outer disk slightly more realistic than the simple one-temperature component we used before by considering a multi-temperature blackbody spectrum as its representation (e.g., Shakura & Sunyaev 1973; Mitsuda et al. 1984; Makishima et al. 1986). We then explore its possible geometry using the data sets where its signature seems the most prevalent.

A significant problem for modeling GX 339$``$4 is that most of its system parameters are not well constrained. The only report of an orbital period based on optical photometry is from Callanan et al. (1992), who found 14.8 hours. This value remains unconfirmed (though Cowley et al. 2002 reconsidered the problem and found a similarly tentative period of 16.8 hours). However, it has been adopted over the years to limit the mass function, which is also problematic since the companion star is still unclassified. Finally, because the system has not yet been resolved, the inclination angle is also unknown but limited to $`60^{}`$ because of the lack of eclipsing. Based on reasonable assumptions, a distance of 4 kpc has been adopted (see Zdziarski et al. 1998; W99; but also Shahbaz et al. 2001, who argue for $`5.6`$ kpc) and a mass for the black hole ranging from $`36M_{}`$, where we here use $`5M_{}`$.

We fix the inner temperature of the SD, $`T_{\mathrm{in}}`$, based on the 1981 data set. We then estimate the outer radius, $`r_{\mathrm{out}}`$, of the SD based on the physical parameters above, the best fit inclination angle, $`\theta _i`$, and assuming a companion of $`0.4M_{}`$, consistent with the 16.8 hour period. The exact luminosity contributed by the SD, $`L_\mathrm{d}`$, is not uniquely determined by the data, and so we can only roughly estimate the accretion rate, $`\dot{M}_\mathrm{d}`$. Because of the simplicity of our disk model, this is not a main focus of this work, but we want to be sure we are consistent in terms of the energy budget. We express $`L_\mathrm{d}`$ in units of $`L_{\mathrm{Edd}}=1.25\times 10^{38}\left(\frac{M_{\mathrm{BH}}}{M_{}}\right)`$ erg s-1.

For accreting black holes it has been argued that the jet power is of order $`Q_\mathrm{j}q_\mathrm{j}\dot{M}_\mathrm{d}c^2`$ with an efficiency inferred to be of order $`q_\mathrm{j}=10^310^1`$ (Falcke & Biermann, 1999; Meier, 2001). We consider $`Q_\mathrm{j}`$ a free parameter in our model, also expressed in units of $`L_{\mathrm{Edd}}`$ erg s-1, and which we check for consistency with our limits on $`\dot{M}_\mathrm{d}`$. While the jet formation itself is still a looming question in the field, the physics of calculating most of the jet emission is relatively straightforward because the flat-to-inverted spectrum stems from the part of the jet where it is basically undergoing free expansion. We briefly summarize the jet component as follows (and see also Falcke & Markoff 2000; MFF).

At the inner edge of the optically thin accretion flow, hot plasma is ejected out from symmetric nozzles, where it becomes supersonic. We assume free jets which accelerate along their axes only due to their longitudinal pressure gradients. This is the simplest scenario and provides a lower limit to the final velocity distribution. The velocity field along the jet is thus uniquely determined from the Euler equation (see, e.g., Falcke 1996), and at large distances $`z`$ from the base has a dependence $`\gamma _\mathrm{j}(z)\beta _\mathrm{j}(z)\sqrt{\mathrm{ln}(z)}`$. Thus, the maximum value of the bulk Lorentz factor $`\gamma _\mathrm{j}`$ is roughly determined by the distance at which the lowest frequency radio emission is emitted and is dependent on other parameters in this model but is itself not a fitted parameter. For both GX 339$``$4 and XTE J1118+480, we find $`\gamma _{\mathrm{j},\mathrm{max}}3`$, which is consistent with the low beaming factors suggested by the recent work of Gallo et al. (2002).

We further assume the jets expand sideways with their initial proper sound speed, $`\gamma _\mathrm{s}\beta _\mathrm{s}c0.4c`$ for a hot electron/cold proton plasma. Beyond the nozzle region, the radius $`r(z)`$ is related to the distance $`z`$ by the Mach number along the jet, $`\gamma _\mathrm{j}\beta _\mathrm{j}/\gamma _\mathrm{s}\beta _\mathrm{s}`$.

We choose the simplifying assumption that the internal energy, dominated by the magnetic field, is equal to the bulk kinetic energy of the particles, consistent with a magnetic launching mechanism. The plasma is assumed to originate in the hot accretion flow and therefore contains equal numbers of protons and electrons, with the exact temperature of the electrons at the base of the jet, $`T_e`$, remaining a fitted parameter. A process such as pair-loading via interactions with the disk photon field will not be efficient enough in the LHS to be significant (see, e.g., Sikora & Madejski 2000). Under these assumptions, the magnetic energy density and particle density in each jet can thus be defined as $`B^2(z)/8\pi 0.25Q_\mathrm{j}/(c\gamma _\mathrm{j}\beta _\mathrm{j}\pi r(z)^2)`$ and $`n^2(z)0.25Q_\mathrm{j}/(c\gamma _\mathrm{j}\beta \mathrm{j}m_pc^2\pi r(z)^2)`$. Taking into account the non-constant velocity field, the dependencies of the magnetic field, $`B`$, and density, $`n`$, on distance are then similar to, but slightly stronger than, the canonical $`z^1`$ and $`z^2`$ dependencies for conical jets, respectively (Blandford & Königl, 1979; Hjellming & Johnston, 1988; Falcke & Biermann, 1995). In this way, the basic physical properties governing the emission at each point in the jet are fixed after specifying the jet power $`Q_\mathrm{j}`$ and nozzle scale $`r_0`$, which determine $`B_0`$ and $`n_0`$.

The thermal particles travel along the jet until they encounter an acceleration region which begins at $`z_{\mathrm{sh}}`$. An acceleration process is inferred because of the observed optically thin synchrotron power-law seen in the recent transient radio events, indicating the presence of non-thermal particles. This is a standard process invoked to explain many features of AGN jet emission (Marscher & Gear, 1985). We consider the process to be diffusive shock acceleration, but it could really be one of a few possible scenarios (e.g., stochastic acceleration) which lead to a non-thermal power-law with spectral index $`p`$ ($`\frac{dN}{dE}E^p`$). For standard shock acceleration theory, $`p2.02.5`$ (e.g., Jones & Ellison 1991). In order to maintain the non-thermal particles out along the jet to account for the entire radio spectrum, we assume distributed acceleration after $`z_{\mathrm{sh}}`$. This is because the cooling timescales after acceleration are very fast compared to the dynamical timescales, a problem which has also required the assumption of distributed acceleration in the two best studied AGN jets, M87 and 3C 273 (Meisenheimer et al. 1996; Jester et al. 2001, but see also Perlman et al. 2001).

We assume that the accelerated distribution extends from the peak of the thermal distribution, with Lorentz factor $`\gamma _\mathrm{e}4T_{\mathrm{e},10}`$. As discussed in MFF, the maximum $`\gamma _\mathrm{e}`$ occurs when the acceleration rate (assuming Bohm diffusion) is matched by the cooling rate giving

$$\gamma _{\mathrm{e},\mathrm{max}}10^8\left(\xi B\right)^{0.5}\left(\frac{u_{\mathrm{sh}}}{c}\right),$$
(1)
where $`u_{\mathrm{sh}}`$ is the shock speed in the plasma frame and $`B`$ is the magnetic field in Gauss. The parameter $`\xi `$ is the ratio between the diffusive scattering mean free path and the gyroradius of the particle. It measures how many times the particle gyrates per shock crossing, where the most efficient acceleration possible is when it crosses the shock once per gyroradius, thus giving a lower limit of $`\xi =1`$. For quasi-perpendicular shocks $`\xi c\gamma _e\beta _e/u_{\mathrm{sh}}`$, although in the case of a jet the field is likely quite tangled in which case an upper bound on $`\xi `$ is not known. While we do not know the exact value of $`u_{\mathrm{sh}}`$, it must be larger than the speed of sound at the shock, $`c_\mathrm{s}`$, and no higher than the local bulk velocity. For very relativistic particles $`\xi `$ can be quite large, but conservative values lie in the $`10100`$ range (e.g., Jokipii 1987). If we define as a reference value $`\xi =\xi _2100`$, the maximum synchrotron frequency is

$$\nu _{\mathrm{max}}=0.29\nu _\mathrm{c}1.210^{20}\xi _2^1\left(\frac{u_{\mathrm{sh}}}{c}\right)^2\mathrm{Hz}$$
(2)
where $`\nu _\mathrm{c}\frac{3}{4\pi }\gamma _{\mathrm{e},\mathrm{max}}^2(eB)/(m_\mathrm{e}c)`$ is the critical synchrotron frequency. This maximum corresponds approximately to the rollover of the power-law cutoff and in order to explain the “canonical” 100 keV cutoff for $`\xi 100`$, we need $`u_{\mathrm{sh}}c_\mathrm{s}`$, consistent with our physical expectations. This cutoff is not dependent on the magnetic field, the jet power, or the shock location as long as we are in the synchrotron cooling dominated regime. Because we would expect XRBs to have similar shock structures, we should get roughly similar cutoffs for different sources and accretion rates. In this sense $`\xi `$ and $`u_{\mathrm{sh}}`$ provide a reasonable physical explanation for the location of the cutoff in much the same way the coronal energy fraction and temperature determine it in coronal IC models (e.g., Sunyaev & Trümper 1979).

The location of the initial acceleration region is determined by the frequency where the flat, highly self-absorbed synchrotron spectrum turns over into the optically thin power-law produced at the shock. In MFF, we showed that from back-extrapolating the X-ray power-law, this maximum self-absorption frequency should occur somewhere in the IR/optical regime at $`10^{14}10^{15}`$ Hz for XTE J1118+480. It is interesting to note that this now holds true for several other BHC sources in the LHS including GX 339$``$4 and V404 Cyg (CF02; Brocksopp et al. 2002, in prep., Gallo et al. 2002, Gallo et al., MNRAS, in prep.).

There are thus seven free parameters relating to the jet component, $`Q_\mathrm{j}`$, $`T_e`$, $`z_{\mathrm{sh}}`$, $`r_0`$, $`p`$ and $`\theta _i`$. These are not, however, entirely independent of each other. Furthermore, in our philosophy of comparing broadband spectral features to the jet model, we can choose reasonable constraints for a number of these parameters. For instance, if the X-rays are due to synchrotron emission, the spectral index of the accelerated particles is uniquely determined ($`\alpha _X=\frac{(p1)}{2}`$). Starting with reasonable assumptions of the underlying disk/jet geometry in order to constrain $`T_e`$ and $`r_0`$, $`Q_\mathrm{j}`$ and $`\theta _i`$ can then be determined from the radio spectrum because they fix its normalization and spectral index, respectively. The beginning of the acceleration region, $`z_{\mathrm{sh}}`$, is then uniquely determined by where the turnover must occur to link the X-ray to the radio spectrum. For the purposes of the spectral models that we describe here, the critical fitted parameters are $`Q_\mathrm{j}`$ and $`\theta _i`$, as they most directly determine the radio spectrum. For $`\theta _i`$, the only observational constraint is that it is $`60^{}`$.

As described above, there are also additional parameters describing the cutoff in the X-rays as well as the disk component, and we will also argue below for an additional irradiation component. However, these are not an integral part of our model, in that they do not constrain the conditions of the radiating plasma in the jet. We include them for completeness, and to open the way for more complex, combined disk/jet models.

Some fraction of the emitted X-rays from the jet base will either directly interact with, or be scattered by hot electrons into, the cold disk, resulting in reflection. This calculation is not attempted here, however, since again this is very dependent on the disk geometry and we are focusing on the role of the jet. This, as well as other issues, will be the focus of an upcoming paper.

## 4 Results

### 4.1 1981 data set: constraining the GX 339$``$4 model parameters

Of all the data sets for GX 339$``$4 considered here, the most interesting is that of May 1981 (see Fig. 1), which is at the highest flux so far for which simultaneous broadband data have been published in the LHS. The IR/optical wavebands seem to indicate the clear presence of both jet and disk components, which has not been seen explicitly in any other XRB. The IR band shows a first component with a negative slope suggesting that we are seeing the expected turnover from the optically thick to optically thin regimes. The shape of this spectrum at higher frequencies is hidden under a component which is likely due to thermal emission from the SD, as indicated by the sharp rise in the optical points. However, if the simultaneously measured X-rays are traced back to the IR, they line up with the turnover remarkably well, supporting their interpretation as synchrotron emission.

While no simultaneous radio observation exists, an estimate of the 8.6 GHz radio flux can be found by extrapolating the radio/X-ray correlation curve to higher fluxes (C00; C02 and see Fig. 6), and is shown as a point in Figs. 1 & 2). As mentioned above, several other LHS sources show this interesting “turnover coincidence”. In the compact cores of most flat-spectrum AGN, the same kind of turnover is observed in the mm range (e.g., Bloom et al. 1994), but the more compact scales in XRBs push this to higher frequencies (e.g., Falcke et al. 2001).

The clear presence of the optically thick-to-thin turnover makes the 1981 data set a strong test of any theoretical model. Furthermore the radio/X-ray correlation, illustrated by the full set of simultaneous data shown here, suggest that any model for the X-rays from this source must address the lower-frequency data as well. This means that if—as is usually assumed—the X-rays are due entirely to IC upscattering off hot electrons in the base of the jet/coronal plasma, one has to show that these same electrons can reproduce the flat-to-inverted radio spectrum via their synchrotron emission, or that some other very strong coupling exists between “corona” and “jet”. As Blandford & Königl (1979) explain, to fit the flat-to-inverted radio spectrum particularly in the cm-band, one needs a conical geometry with outflowing plasma, otherwise the emitted synchrotron spectrum will be too inverted. We therefore think that a corona which has no relation to the jet cannot explain the data for GX 339$``$4 . However, the corona may provide the launching point for the jet (Merloni & Fabian, 2002) or even comprise the base of the jet itself (Fender et al., 1999). Interestingly, Beloborodov (1999) found that magnetized plasma moving away from the disk with a velocity of $`0.3c`$ can explain the weak reflection features in Cyg X-1, which would be consistent with this picture.

In order to understand the IR/optical data, we must first consider how the disk may be contributing in this region. Unfortunately, addressing the thermal disk component is not a uniquely defined problem for GX 339$``$4 . If we assume the SD component extends up to the first X-ray data points and then drops off, we find $`T_{\mathrm{in}}=10^6`$ K, which is the same value which W99 used for their fits to this source. We keep $`T_{\mathrm{in}}`$ fixed to this value for all models discussed in this paper. However, a naive fit to the rising optical data with a multi-temperature blackbody would then require almost the entire $`L_{\mathrm{Edd}}`$ being channeled into the SD emission. For the LHS, which is defined by its assumedly low accretion rate, this seems unrealistically high. One way around this is to assume a lower SD luminosity, and that the optical points are due to the irradiation of the outer disk material by intercepted X-rays. For the mass and inclination angle we use here, the outer disk radius would be on the order of $`r_{\mathrm{out}}10^5`$ cm. We use this as the nominal scale for the X-ray irradiated shell, noting that there could easily be a factor of a few leeway due to the large uncertainties in the system’s binary parameters. To fit the optical data, the temperature of this single blackbody needs to be $`T_{\mathrm{BB}}5\times 10^4`$ K. This high temperature is not in conflict with possible disk instability models (e.g., Cannizzo et al. 1988) because we are observing already in outburst, after the effects of the instability would be felt.

For the models presented in Fig. 1, we fix the SD luminosity to a characteristic value of $`L_\mathrm{d}0.1L_{\mathrm{Edd}}`$, and then make up any deficit in the optical region via the irradiation component. The ensuing fraction of the jet power contained in the irradiated component then can be determined due to the fact that an IR/optical break is observed. For the other data sets included in this paper, however, this break is not observationally resolved. For those models we therefore assume that the luminosity of the irradiated emission can be linked to the power in the jet, and we set $`L_{\mathrm{BB}}=0.1Q_\mathrm{j}`$ as a canonical value.

It turns out that the 1981 data set is in fact so luminous that it creates a problem for our assumed black hole mass of $`5M_{}`$. Even assuming $`L_\mathrm{d}0.1L_{\mathrm{Edd}}`$ in fact requires $`\dot{m}0.5`$ (in Eddington units) for the fixed $`T_{\mathrm{in}}=10^6`$ K, which is much larger than one would expect for the LHS according to some models (e.g., Esin et al. 1997). We could solve this problem by choosing a lower $`L_\mathrm{d}`$ and $`T_{\mathrm{in}}`$, thus having the SD contribution truncate before the X-ray data. However, then this only increases the need for more intercepted jet power in the irradiated component, which we already require to be at the upper bound of physical limits. We discuss the parameters in detail for the models below, but feel this problem is too unconstrained at this time to make any conclusions about the exact nature of the disk contribution.

### 4.2 1981 data set, no acceleration

Assuming that jets are responsible for the radio emission, the first issue to explore is how the solution would look if the jet plasma never encounters an acceleration region. This would be the simplest case model, where a portion of the thermal disk or coronal plasma is simply advected into the outflow, retaining its thermal or quasi-thermal distribution. In this case, the synchrotron spectrum from each segment of the jet will have a “hump” shape peaked around a particular self-absorption frequency determined from the physical parameters for the segment. The combination of all these humps, each located at a different frequency, gives a flat-to-inverted radio to at least IR spectrum which can be seen in Fig. 1. The slope of the spectrum depends on the amount of expansion and the velocity profile along the jet, as well as on $`\theta _i`$.

As explained in the previous section, the electron temperature at the base of the jet is a free parameter, and so we explore two extremes for this scenario. The first, shown in Fig. 1a, is at the lowest end of the relativistic scale with electrons entering the jet base with initial $`T_\mathrm{e}2\times 10^9`$ K, and in Fig. 1b we show a more relativistic solution, with $`T_\mathrm{e}=2\times 10^{10}`$ K, the maximum temperature we could accommodate in our modeling of XTE J1118+480 (MFF). In order to bring the highest frequency “hump” from the part of the jet closest to the accretion disk to low enough frequencies to come near the turnover in the IR data, we must vary the nozzle radius $`r_0`$ and the input power $`Q_\mathrm{j}`$, thus controlling the magnetic field $`B_0`$ and density $`n_0`$ at the base of the jet while the temperature remains fixed. Similarly we must adjust $`\theta _i`$ to fit both the radio and the IR together with synchrotron, which changes the power requirements because of the beaming.

With the requirement that the non-accelerated models adequately describe the IR data as well as the maximum self-absorption turnover, the jet parameters for the fit shown in Fig. 1a are then $`r_0=5\times 10^3r_g`$, $`Q_\mathrm{j}=3.7L_{\mathrm{Edd}}`$ and $`\theta _i=8^{}`$, the former two being obviously unlikely values. The extremely high value of the the nozzle radius is necessary to decrease the self-absorption frequency to a value low enough to fit the IR data, which because of the resulting larger self-absorption also leads to a very peaked contribution from the denser nozzle. This does not seem feasible under the assumption that the jet originates at the inner edges of the accretion flow, or in the corona above it, and the match to the IR data itself is not particularly good. The “nozzle bump” is due to our assumption of this element of the geometry, as required for jet models of Sgr A\* and LLAGN (e.g., Falcke & Markoff 2000; Yuan et al. 2002). Without a nozzle region, one could find a solution with a somewhat larger $`\theta _i`$, but then this would require an even higher $`Q_\mathrm{j}`$ because the emission would be beamed away from the observer. The value of $`Q_\mathrm{j}`$ here is already several times higher than $`L_{\mathrm{Edd}}`$ for a 5 $`M_{}`$ central object. The irradiation component can provide the required rise in the optical data if the outer disk intercepts a small fraction of the jet power in X-rays, $`L_{\mathrm{BB}}0.03Q_\mathrm{j}`$. Overall, however, this model is not so satisfactory.

The observations are slightly better described, at least at low-frequency, for the more relativistic case shown in Fig. 1b. The higher temperature electrons lead to less self absorption, giving a broader nozzle bump and requiring a somewhat smaller nozzle radius of $`r_0=1.2\times 10^3r_g`$ with a more reasonable $`Q_\mathrm{j}=1.4\times 10^2L_{\mathrm{Edd}}`$ for $`\theta _i=15^{}`$. This nozzle width is still rather untenable, however, in the context of jet/corona models. For this scenario we also run into more trouble with the disk modeling, because we would require the outer disk to intercept roughly $`6Q_\mathrm{j}`$ in X-ray luminosity to account for the optical data.

For both cases shown in Fig. 1, the IC emission in the X-rays is dominated by the upscattered thermal disk emission but, because of the low densities, lies under the data and also has the wrong shape. In this and all other figures in this paper, we do not include IC of photons from far out in the jet, since they are orders of magnitude under the data and will not affect the fit. Changing the disk geometry so that the thin disk underlies the base of the jet/corona may provide more IC emission, but would still not address the fundamental problem of explaining the radio thru IR data with reasonable parameters.

In summary, the energy budget in the non-accelerated case is not large enough to create enough synchrotron to fit the radio for a non-relativistic plasma, since the required $`Q_\mathrm{j}`$ for the already weakly relativistic temperature of Fig. 1a is super-Eddington. Considering that the radio through IR bands are likely due to jet synchrotron in the LHS, we conclude, based on all the above difficulties, that a model which does not include particle acceleration in the jet is unlikely. The extent to which the energy budget falls short is in part due to assumptions of the equipartition of energy, as described earlier. However, corona or other accretion flow models which attempt to explain the low-frequency data without a standard jet solution do not seem very feasible.

### 4.3 1981 data set, with acceleration

Therefore we want to focus on solutions where the particles in the jet are accelerated into a power-law distribution, as is typical in AGN. This solves the energy budget problem for the jet power, although we still need a very luminous thermal contribution.

We show four representative fits in Fig. 2, all with a nozzle radius of $`3r_g`$, which is more appropriate for the assumption that the jet originates at the edge of an accretion flow. As discussed above, we fix $`L_\mathrm{d}=0.1L_{\mathrm{Edd}}`$ with $`T_{\mathrm{in}}=10^6`$ and $`r_{\mathrm{out}}=10^5r_g`$ and then fit the optical data with the irradiated component, leaving the intercepted fraction of $`Q_\mathrm{j}`$, $`L_{\mathrm{BB}}`$, as a free parameter in the fit.

In panel (a) of Fig. 2, we show the resulting spectrum from a jet with mildly relativistic electrons, similar to the model shown in Fig 1a. The short-dashed line shows their contribution before encountering the acceleration region. Compared to Fig. 1 the nozzle “hump” lies at a higher frequency due to the more compact nozzle scale, and the slope is more inverted primarily to the greatly increased angle to the line-of-sight necessitated by the radio/IR fit, $`\theta _i=55^{}`$. The slope is flatter in the post-accelerated component (shown with the dot-dashed line), because of the effects of reacceleration and particle losses. Synchrotron from the accelerated particles dominates the entire spectrum with the exception of the optical, with the IR and X-rays stemming from the region near the start of the acceleration $`z_{\mathrm{sh}}`$, while the radio comes from much further out along the jet. For a fixed $`T_\mathrm{e}`$, $`Q_\mathrm{j}`$ and $`\theta _i`$ are found by fitting the radio-IR spectrum, which then uniquely determines $`z_{\mathrm{sh}}`$ by fitting the IR turnover and extension to the X-rays. This also determines the slope of the accelerated particles, $`p`$. The jet parameter values for this mildly relativistic case are: $`T_\mathrm{e}=2\times 10^9K`$, $`Q_\mathrm{j}=0.62L_{\mathrm{Edd}}`$, $`\theta _i=55^{}`$, $`z_{\mathrm{sh}}=2.15\times 10^3r_g`$ and $`p=2.15`$. The IC emission is dominated by upscattered thermal disk photons and is well under the synchrotron contribution.

The problem with this model is that—as in the model from Fig. 1a—the required power in the jet $`Q_\mathrm{j}`$ is unrealistically large, in order to compensate for the low electron temperature. This imbalance also leads to the magnetic field energy density in the jet being a few orders of magnitude super-equipartition with the hot, radiating electrons. While many jet formation models assume the magnetic energy density is in equipartition with the kinetic energy near the base of the jet (e.g., Blandford & Payne 1982; Li et al. 1992), which for a maximal jet would lead to magnetic domination, this may be rather extreme.

In Fig. 2b, we consider the highest temperature electrons allowed by this model, above which the X-ray predictions (from either synchrotron or IC from the jet) are too high, $`T_\mathrm{e}=7\times 10^9`$ K. The jet free parameters are now: $`Q_\mathrm{j}=0.27L_{\mathrm{Edd}}`$, $`\theta _i=55^{}`$, $`z_{\mathrm{sh}}=1.75\times 10^3r_g`$ and $`p=2.15`$. Raising the temperature thus serves mainly to reduce the power requirements and changes the location of the shock region slightly. With the same SD parameters as in (a), the irradiation component now comprises $`0.4Q_\mathrm{j}`$, still rather large but again given the uncertainties in the mass and outer radius, more viable. The hotter electrons now also contribute to a greater degree via IC upscattering, contributing some of the soft X-ray emission.

While there is no measurement of the inclination angle, some papers favor a lower value than what we find for our best fit above (e.g., Cowley et al. 2002). In Fig. 2c, we show what our model predicts for $`\theta _i=15^{}`$ and the same temperature as in panel (b). The decreased angle greatly reduces the power requirements in the jet, but at the same time flattens the predicted radio-IR slope significantly, making a good radio/IR fit impossible. The parameters are $`Q_\mathrm{j}=3.2\times 10^2L_{\mathrm{Edd}}`$, $`z_{\mathrm{sh}}=2.15\times 10^3r_g`$ and $`p=2.15`$. Although this radio point is only from an extrapolation of the radio/X-ray correlation relation (C00; C02), the slope is remarkably well-defined at these higher luminosities. On the other hand, $`\theta _i=55^{}`$ does give good fits to most of the other data sets discussed below. Because the beaming reduces the power requirements for the jet, a low inclination angle also means that more power must irradiate the outer disk than is available, making this model problematic.

Finally in Fig. 2d we consider the case of a power-law of particles entering the base of the jet from the accretion flow instead of thermal particles, as may occur if acceleration also takes place at an accretion shock. Instead of defining a temperature, we take the minimum electron Lorentz factor $`\gamma _{\mathrm{e},\mathrm{min}}=1`$, leaving $`\gamma _{\mathrm{e},\mathrm{max}}`$ a free parameter. In order to prevent the nozzle emission from violating the constraints of the soft X-ray data, we find $`\gamma _{\mathrm{e},\mathrm{max}}10`$, which would require extreme radiative cooling to truncate the accelerated power-law for this scenario to be realistic. This type of scenario may be necessary for certain low-luminosity AGN (see, e.g., Yuan et al. 2002). The other jet parameters are $`Q_\mathrm{j}=0.55L_{\mathrm{Edd}}`$, $`z_{\mathrm{sh}}=2.15\times 10^3r_g`$, $`\theta _i=55^{}`$ and $`p=2.15`$. The high required jet power means that the irradiated component needs only $`0.02Q_\mathrm{j}`$ to account for the optical data.

In this last model, the IC emission still lies under the synchrotron component, but not by much, and has the correct slope. In a different geometry where the SD underlies the corona/jet base, this may provide the best alternative for a solution where the IC dominates the synchrotron emission in the X-rays.

In all four model runs shown here, we needed a rather large fraction of the jet power to be intercepted by the outer disk and reradiated. If we allowed the SD to radiate $`0.2L_{\mathrm{Edd}}`$ away, we could lower the necessary $`L_{\mathrm{BB}}`$, but then again we would be basically at the Eddington accretion rate, as discussed for the non-acceleration case. This problem is clearly not very constrained for this source, and a central mass of $`15`$ instead of 5 $`M_{}`$ (as in, e.g., GRS 1915+105; Greiner et al. 2001), would easily resolve these issues.

In Fig. 3 we show the same model as in Fig. 2b with the flux multiplied by the frequency to show the distribution of power across the wavebands. One can see that most power seems to fall in the disk component, with the synchrotron then creating the second peak around the cutoff. The IC component is in fact dominated by the upscattered disk photons, but this is hidden beneath the synchrotron emission. The self-Comptonized synchrotron component peaks at $`10^{23}`$ Hz, or $`400`$ MeV, corresponding to the synchrotron peak upscattered by the hot electrons near the shock. At high frequencies, the shape is similar to the “camel back” spectra seen in blazars (e.g., Fossati et al. 1998), and in fact at this highest observed flux level is above the sensitivity of GLAST for a one year all-sky survey at 100 MeV. The less powerful epochs of this source presented in the next section may not be so observable, however, as they are orders of magnitude lower.

In conclusion, we have shown that in order to explain the low-frequency data from the 1981 observation of GX 339$``$4 , a jet solution must include some form of particle acceleration. And further, once that condition is fulfilled, X-rays from jet synchrotron (in addition to IC) are unavoidable as long as IC emission is not the main cooling channel. But under the assumption of a disk geometry where a SD transforms to an optically thin, hot accretion flow at some inner radius, synchrotron cooling will dominate in a maximal jet. It is possible, however, that in the alternative geometry of a corona/jet base lying directly on top of the SD could lead to so much IC cooling that the synchrotron cutoff would occur below $`10^{17}`$ Hz, thus accommodating a dominant IC contribution in the X-rays. This would not, however, explain the interesting coincidence of IR turnover/X-ray spectrum seen here and now also in several sources in addition to XTE J1118+480. It may also be difficult to explain the radio/X-ray correlation with an IC model (see Sect. 5). We will, however, soon explore combined jet/corona solutions using lower temperature cyclo-synchrotron and IC scattering in different geometries.

We have explored the relativistic electron parameter range for the 1981 data set and find that the most reasonable parameters are obtained for the case of a fully relativistic initial thermal electron distribution. Besides requiring too much power, a mildly relativistic initial electron distribution would also result in the magnetic energy density being much higher (and further out of equipartition with) the accelerated particles at the shock, more so than would likely be expected for magnetic confinement. For these reasons, in this paper we consider only the higher temperature solutions for all fits, corresponding to the type of fit shown in Fig. 2b.

### 4.4 Other data sets

For the 12 data sets compiled during the other LHS episodes of GX 339$``$4 , we can either optimize each fit individually leaving all parameters free, or fix as many parameters as possible to see if essentially the same model can account for all observed spectra. We here choose to take this latter approach, in order to explore the fundamental physics of the model. We will fix some of the parameters to those best constrained by the 1981 set and then study the affect of changing only two jet parameters. Our emphasis is on understanding the radio/X-ray correlations, and thus the jet component. However, we address the SD contribution as well, and attempt to test the jet/disk connection by finding a scaling solution. We set $`L_{\mathrm{BB}}=0.1Q_\mathrm{j}`$ for all fits, and then fit $`L_\mathrm{d}`$ to the optical data, when present.

While it seems clear that there are small variations in the soft X-ray spectra of these data sets (this is discussed in NWD, who fit a broken power law—we are fitting the lower energy segment), we choose to fix the accelerated spectral index to $`p=2.15`$ as in 1981. Similarly, while there are obvious changes in the SD contribution, they seem consistent with the likely underlying accretion changes that also affect the jet since both components decrease in luminosity with increasing time. We keep $`T_{\mathrm{in}}`$ fixed at $`10^6`$ K and $`T_\mathrm{e}`$ fixed at $`7\times 10^9`$ K based on the 1981 set, as a starting point. We will assume that the orientation of the jet remains fixed at $`55^{}`$ (statistically $`57^{}`$ is the most likely angle) and that the base of the jet is fixed at $`r_0=3r_\mathrm{g}`$. Thus the only jet parameters allowed to vary are $`Q_\mathrm{j}`$ and $`z_{\mathrm{sh}}`$, and the disk parameter $`L_\mathrm{d}`$ is then adjusted to fit the thermal data, when present.

Figs. 4a & 4b show the broadband spectral fits to the data, with the corresponding fit parameters given in Table 2. For the data sets without optical data constraining the disk, $`L_\mathrm{d}`$ is rather arbitrary, and we scaled the luminosity to the jet power.

Although the fits are surprisingly good given the limited fitting parameters, there are obvious features in the X-ray waveband such as a flattening at higher frequencies and peaks which our model cannot account for. These features are likely due to the impacting X-rays on the cooler material of the SD, but we have not yet incorporated this into our model. The softest X-rays should not be affected by these additional components, however, so we focus our fitting at the lower X-ray frequencies. For the lowest flux level data sets of 1999, which we label $`99_4`$$`99_8`$, as well as 2000, the error bars are large enough to make the fitting more difficult.

The most important feature to note is that the radio/X-ray correlations detected in this source (C00; C02) can be explained by the variation of only one main parameter in the jet model: the input power. If we plot $`z_{\mathrm{sh}}`$ against the jet power $`Q_\mathrm{j}`$ (Fig. 5), we see that there is a slight dependence ($`z_{\mathrm{sh}}Q_\mathrm{j}^{0.135}`$), but this likely results from our rather artificially freezing of the other parameters. This small dependence is rather amazing when one thinks that $`Q_\mathrm{j}`$ is changing by about three orders of magnitude over the same range, if 1981 is included.

## 5 Understanding the radio/X-ray correlations

It seems likely that the real source of the variations is the power fed into the jet. In Fig. 6 we plot the prediction of our model against the data for the correlation between the 8.6 GHz radio flux and the integrated 3-9 keV X-rays (C02). We fix $`z_{\mathrm{sh}}1.8\times 10^3r_g`$, to match the highest flux data point of 1981, which is not shown on this figure because the radio was not simultaneously measured. Then we adjust only $`Q_\mathrm{j}`$ while holding every other parameter constant, in order to fit the radio flux, and plot the resulting integrated 3-9 keV X-ray flux. In removing the $`z_{\mathrm{sh}}Q_\mathrm{j}^{0.135}`$ dependence mentioned above, we are likely compromising the quality of the fit to some degree, but this could be compensated for if we allowed particularly the spectral index to vary.

By just changing the jet power we obtain a very good fit to the correlation data, with the exception of the two data sets $`99_7`$ and $`99_8`$. In fact, the radio data in $`99_8`$ was taken a few days after the X-rays, and coincided with a small reflare in the hard X-rays (see Fig. 15 in C00). This is the likely explanation of the excess radio flux. Similarly, the reflare may also have an affect on $`99_7`$, if there is, e.g., a lag between the changing fluxes.

The results shown in Fig. 6 provide strong support for a jet synchrotron model, since it correctly predicts the radio/X-ray correlation over orders of magnitude in flux levels, almost exactly matching the data. The fit is better at the higher flux end, where the error bars are smaller. The exact dependence (excluding $`99_7`$ and $`99_8`$) goes as $`F_\mathrm{X}F_\mathrm{R}^m`$, where $`m=1.41`$, $`F_\mathrm{R}`$ is the 8.6 GHz radio flux and $`F_\mathrm{X}`$ is the integrated 3-9 keV X-ray flux.

The value for $`m`$ in fact follows directly from analytic predictions of the jet model. Eq. (52) in Falcke & Biermann (1995) shows111In Falcke & Biermann (1995), the quantities which we label $`\nu _{\mathrm{SSA}}`$, $`Q_\mathrm{j}`$ and $`F_{\mathrm{SSA}}`$ correspond to $`\nu _{\mathrm{s},\mathrm{obs}}`$, $`q_{\mathrm{j}/1}L_{46}`$ and $`L_{\nu _s,\mathrm{obs}}`$, respectively. that if all parameters except the power are fixed, as we have done to produce Fig. 6, the frequency where the optically thick flat-to-inverted spectrum turns over to the optically thin regime depends on the power as

$$\nu _{\mathrm{SSA}}Q_\mathrm{j}^{2/3}.$$
(3)
Similarly Eq. (56) (ibid.) shows that under the same freezing of other parameters, the flux at the turnover depends on the power as

$$F_{\mathrm{SSA}}Q_\mathrm{j}^{17/12}.$$
(4)
These relations can then be normalized using the 1981 data set, where the flux and turnover frequency are visible.

Knowing these relations, and the spectral indices of the optically thick and thin components from the data, one can roughly calculate the expected correlation slope. Because the spectral indices for the radio-IR, $`\alpha _{\mathrm{RIR}}`$, and X-ray, $`\alpha _\mathrm{X}`$, do vary slightly over GX 339$``$4 ’s observation history, taking constant values will make this less exact. We now want to calculate the slope between two points along the correlation curve, analogous to the model plotted in Fig. 6. Each point represents a different value of $`Q_{\mathrm{jet}}`$, corresponding to different values of $`\nu _{\mathrm{SSA}}`$ and $`F_{\mathrm{SSA}}`$. The 8.6 GHz radio flux is then

$$F_{8.6}=F_{\mathrm{SSA}}\left(\frac{\nu _{\mathrm{SSA}}}{8.6\mathrm{GHz}}\right)^{\alpha _{\mathrm{RIR}}},$$
(5)
and then using Eqs. 3 & 4:

$$\mathrm{log}_{10}F_{8.6}=\mathrm{log}_{10}C_1+\frac{17}{12}\mathrm{log}_{10}Q_\mathrm{j}\frac{2}{3}\alpha _{\mathrm{RIR}}\mathrm{log}_{10}Q_\mathrm{j},$$
(6)
where $`C_1`$ absorbs the exact dependences of the flux and frequency at the turnover on $`Q_\mathrm{j}`$, and other constants.

Similarly, the X-ray flux at any frequency (as long as it falls on the power-law) is

$$F_\mathrm{X}=F_{\mathrm{SSA}}\left(\frac{\nu _{\mathrm{SSA}}}{\nu _\mathrm{X}}\right)^{\alpha _\mathrm{X}},$$
(7)
which after using the forms from Eqs. (3) & (4) and integrating from 3-9 keV becomes

$$\begin{array}{c}\mathrm{log}_{10}F_\mathrm{X}𝑑\nu =\mathrm{log}_{10}C_2+\frac{17}{12}\mathrm{log}_{10}Q_\mathrm{j}\hfill \\ \hfill \frac{2}{3}\alpha _\mathrm{X}\mathrm{log}_{10}Q_\mathrm{j},\end{array}$$
(8)
where again $`C_2`$ includes the constants of integration and the dependencies on $`Q_\mathrm{j}`$.

For any two values of the jet power, $`Q_A`$ and $`Q_B`$, the correlation slope is then

$`m`$ $`=`$ $`{\displaystyle \frac{\mathrm{log}_{10}F_{\mathrm{X},\mathrm{B}}\mathrm{log}_{10}F_{\mathrm{X},\mathrm{A}}}{\mathrm{log}_{10}F_{8.6,\mathrm{B}}\mathrm{log}_{10}F_{8.6,\mathrm{A}}}}`$ (9)
$`=`$ $`{\displaystyle \frac{\frac{17}{12}\frac{2}{3}\alpha _\mathrm{X}}{\frac{17}{12}\frac{2}{3}\alpha _{\mathrm{RIR}}}},`$ (10)
after algebraic cancellations. Taking the values from the power-law fits to the 1981 data set from CF02, $`\alpha _{\mathrm{RIR}}0.15`$ and $`\alpha _\mathrm{X}0.6`$, we find $`\alpha _{\mathrm{cor}}=1.38`$, which considering the simplifications in the model and our ignoring of the spectral variations, is surprisingly close to the value of 1.4 found for GX 339$``$4 . In addition, this type of radio/X-ray correlation has also been seen in V404 Cygni, with exactly the same slope (Gallo et al. 2002; Gallo et al., MNRAS, in prep.). Our preliminary models of this source also show a turnover in the IR band and have similar jet parameters (Markoff et al. 2002, in prep.).

The fact that the well-determined slope $`m`$ extrapolates to the lower flux data also gives us greater confidence in the subtraction of the putative “background source” discussed in Sect. 2. Note also that the lowest flux point is from BeppoSAX, and therefore does *not* rely on the subtraction of a background source.

It is possible that an IC corona model could also explain the correlation slope, if the radio flux scales with the power (assumed to be linearly related to $`\dot{M}`$) as expected from the jet and the X-ray flux scales with $`\dot{M}^2`$. This may be the case in some optically thin accretion solutions. However, it has yet to be shown that a self consistent solution can be found which also well fits all the broadband data available for this source. Similarly this would not explain why the X-rays trace back to the turnover in the IR as 1981 data, or why most if not all LHS sources show this turnover coincidence.

## 6 Conclusions

We show that a model comprising a dominantly synchrotron jet component, in combination with an optically thin accretion flow transitioning to a standard thin disk, is able to explain the broadband spectral data from 13 observations of GX 339$``$4 . At the same time, the model can easily explain the $`m1.4`$ slope of the radio/X-ray correlations by changing only the power input into the jet. This input power, $`Q_\mathrm{j}`$, is assumed to be proportional to the accretion rate, and varies from $`L_{\mathrm{Edd}}0.0030.3`$ for a 5$`M_{}`$ black hole. Above the highest luminosity, the source would likely transition to the HSS.

This model assumes an underlying disk geometry in which the standard thin disk exists beyond $`1001000r_g`$, and has a relatively low accretion rate. In a situation where the SD extends all the way down to the last stable orbit, the photon field at the base of the jet may be high enough to affect our conclusions, which are dependent on synchrotron losses dominating the cooling. If the available photon field for inverse Compton upscattering becomes large enough, it may truncate the accelerated particle distribution and lessen the extent of the synchrotron emission.

With the exception of this model, the non-thermal X-ray spectral component in all XRB states has been modeled in terms of inverse Compton processes in a hot corona, located either above or within the standard thin disk. These classes of models have significant success explaining both the spectral features and timing characteristics in the X-ray waveband, but fail to address the radio/IR spectra which in at least GX 339$``$4 , V404 Cyg and Cyg X-1 seem intimately linked to the higher frequencies.

However, the existence of a hot, magnetized plasma at the base of the jet leads to obvious considerations of a unification scheme between it and this concept of a corona, which until now has not been explicitly observed. Both scenarios (jet and corona) address important features from the observations, but either cannot or have not yet attempted to account for everything. We feel that a concatenation of these two components would be fruitful as the next stage of our investigation. It is, however, important to note that under reasonable physical assumptions, the jet can easily produce X-ray emission via synchrotron radiation, and has significant success explaining the data when it does. This possibility must therefore be considered in spectral modeling; whether it really dominates the entire spectrum is a question of the environment and local acceleration conditions. It may be that it indeed dominates only in those sources which show only weak reflection features (e.g., XTE J1118+480; Miller et al. 2001), but it likely has a wider impact in all sources than has so far been appreciated.

Interestingly, the parameters we derive for the power and location of the shock are roughly consistent with earlier results modeling XTE J1118+480, which because of the higher temperature ($`T_e=2\times 10^{10}`$ K) resulted in a lower necessary power input into the jet. In this source we found $`z_{\mathrm{sh}}10^2r_g`$ as compared to the $`10^3r_g`$ found here, however considering that the jet can extend beyond $`10^{10}r_g`$, this range is quite small. In our recent modeling of other LHS BHC sources, we are in fact finding that all sources with simultaneous radio/X-ray data seem to require acceleration to begin in this same range, and we will discuss the physical implications of this elsewhere. For these sources we also find, as in GX 339$``$4 and XTE J1118+480, that if the X-rays are traced back to lower frequencies, the optically thick-to-thin turnover always occurs in the IR range. If this turnover coincidence is real, it offers us the chance to explore the physical conditions at the acceleration zone of XRB jets.

This work illustrates how critical simultaneous multiwavelength observations are to increasing our understanding of the physics in these sources. With the advent of higher energy missions such as GLAST, INTEGRAL and ASTRO E-2, soon we can hopefully begin to probe the contribution of these smaller jets to the hard X-ray and $`\gamma `$-ray bands.

###### Acknowledgements.

This work is partially supported by AUGER Theory Grant O5CU1ERA/3 from the BMBF (S.M.). S.M. would like to thank Tom Maccarone and David Meier for very helpful discussions.


text:
# New Developments in Interval Arithmetic and Their Implications for Floating-Point Standardization Technical report DCS-273-IR, Department of Computer Science, University of Victoria, Victoria, BC, Canada.

## 1 Introduction

Continuing advances in process technology have caused a tremendous increase in the number of transistors available to the designer of a processor chip. As a result, multiple parallel floating-point units become feasible. The time will soon come when interval arithmetic can be done as fast as conventional arithmetic.

However, to properly utilize the newly available number of transistors, chip designers need to spend ever more time iterating through cycles of synthesis, place-and-route, and physical verification than current design methodology allows. This design bottleneck makes it desirable to simplify floating-point units. Recent developments in the theory of interval arithmetic suggest possibilities for simplification. As far as interval arithmetic is concerned, certain parts the 1985 standard are essential, whereas other parts are superfluous or even a liability.

In this paper we present a consolidated, self-contained account of the new developments in interval arithmetic not available elsewhere. In the conclusions we give a sketch of what would be an ideal standard from the point of view of interval arithmetic when arranged in this way.

## 2 Why interval arithmetic?

Conventional, non-interval, numerical analysis is marvelously cheap and it works most of the time. This was exactly what was needed in the 1950s when computers needed a demonstration of feasibility. A lot has changed since that time. Numerical computation no longer needs to be cheap. It has become more important that it always works. As a result interval arithmetic is becoming an increasingly compelling alternative. For example, civil, mechanical, and chemical engineers are liable for damage due to unsound design. So far, they have been able to get away with the use conventional numerical analysis, appealing to what appears to be best practice. As interval methods mature, it is becoming harder to ignore them when defining “best practice”.

Recent developments, which we call “modern interval arithmetic” provide a practical and mathematically compelling basis. In conjunction with this it has become clear that some aspects of IEEE standard 754 are not needed or are detrimental, whereas other aspects are marvelously suited to interval arithmetic. If these latter are preserved in future development of the standard, then interval arithmetic can help bridge the design gap and lead to the situation where all arithmetic can be faster due to interval methods.

## 3 A theory of approximation

Conventional numerical analysis approximates each real by a single floating-point number. It approximates each of the elementary operations on reals by the floating-point operation that has the same name, but not the same effect. Let us call this approach “point approximation”. It has been amply documented that this approach, though often satisfactory, can lead to catastrophical errors. This can happen because it is not known to what degree a floating-point variable approximates its real-valued counterpart.

Interval arithmetic is based on a theory of approximation, set approximation, that ensures that for every real-valued variable $`x`$ in the mathematical model, there is a machine-representable set $`X`$ of reals that contains $`x`$. Such arithmetic is exact in the sense that $`xX`$ is and remains a true statement in the sense of mathematics. It is of course not exact in the sense that $`X`$ typically contains many reals.

Conversely, in case of numerical difficulties, it will turn out that continued iteration does not reduce the size of $`X`$, in which case we have a notification that there is a problem with the algorithm or with the model. Because of this property, we call this manifest approximation: there is always a known lower bound to the quality of the approximation.

In this way the operations can be interpreted as inference rules of logic; for example, $`xX`$ and $`yY`$ imply that $`x+yZ`$, where $`Z`$ is computed from $`X`$ and $`Y`$.

## 4 Approximation structures

We address the situation where we need to solve a mathematical model in which a variable takes on values that are not representable in a computer, but where it is possible to so represent sets of values. We then approximate a variable by a representable set that contains all the values that are possible according to the model. Models with real-valued variables are but one example of such a situation.

As the theory of set approximation applies to sets in general, we first present it this way.

###### Definition 1

Let $`𝒯`$ be the type of a variable $`x`$, that is, the set of the possible values for $`x`$. A finite set $`𝒜`$ of subsets of $`𝒯`$ that contains $`𝒯`$ and that is closed under intersection is called an approximation structure in $`𝒯`$.

In a typical practical application of this theory, $`𝒜`$ is a set of computer-representable subsets of $`𝒯`$.

###### Theorem 4.1

If $`𝒜`$ is an approximation structure for $`𝒯`$, then for every $`S𝒯`$ there exists a unique least (in the sense of the set-inclusion partial order) element $`S^{}`$ of $`𝒜`$ such that $`SS^{}`$.

###### Definition 2

For every $`S𝒯`$, $`\varphi (S)`$ is the unique least element of $`𝒜`$ that exists according to theorem 4.1.

We regard $`\varphi (S)`$ as the approximation of $`S`$. As $`S`$ can be a singleton set, this theory provides approximations both of elements and subsets of $`𝒯`$.

## 5 An approximation structure for the reals

We seek a set of subsets of the reals that can serve as approximation structure. A first step, not yet computer-representable, is that of the closed, connected sets of reals.

###### Theorem 5.1

The closed, connected subsets of $`R`$ are an approximation structure for $`R`$.

Proof According to a well-known result in topology, all closed connected subsets of $`R`$ have one of the following four forms: $`\{xaxb\}`$, $`\{xxb\}`$, $`\{xax\}`$, and $``$ where $`a`$ and $`b`$ are reals. Here we do not exclude $`a>b`$, because the empty set is also included among the closed connected sets. Clearly the conditions for an approximation structure are satisfied: $``$ is included and the intersection of any two sets of this form is again closed and connected. This completes the proof.

The significance of the closed connected sets of reals as approximation structure is that they can be represented as a pair of extended reals.

###### Definition 3

The extended reals are the set obtained by adding to the reals the two infinities. As with the reals, the extended reals are totally ordered. When two extended reals are finite, then they are ordered within the extended reals as they are in the reals. Furthermore, $`\mathrm{}`$ is less than any real and $`+\mathrm{}`$ is greater than any real.

Theorem 5.1 together with definition 3 suggest the following notation for the four forms of the closed connected sets of reals. In this notation we do not include the empty interval. The reason is that if it is found in an interval computation that an interval is empty, any further operations involving the corresponding variable yield the same result, so that the computation can be halted.

###### Definition 4

Let $`a`$ and $`b`$ be reals such that $`ab`$.

$`a,b`$ $`\stackrel{\mathrm{def}}{=}`$ $`\{xaxb\}`$
$`\mathrm{},b`$ $`\stackrel{\mathrm{def}}{=}`$ $`\{xxb\}`$
$`a,+\mathrm{}`$ $`\stackrel{\mathrm{def}}{=}`$ $`\{xax\}`$
$`\mathrm{},+\mathrm{}`$ $`\stackrel{\mathrm{def}}{=}`$ $``$
Note that each of these pairs denote sets of reals, even though in their notation the infinities are used. These are not reals.

## 6 Floating-point intervals: a finite approximation structure for the reals

Let $`F`$ be a finite set of reals.

###### Theorem 6.1

The sets of the form $`\mathrm{}`$, $`\mathrm{},b`$, $`a,b`$, $`a,+\mathrm{}`$, and $`\mathrm{},+\mathrm{}`$ are an approximation structure when $`a`$ and $`b`$ are restricted to elements of $`F`$ such that $`ab`$.

###### Definition 5

The real intervals are sets of the form described in theorem 5.1. The floating-point intervals are sets of the form described in theorem 6.1, where $`a`$ and $`b`$ are finite IEEE754 floating point numbers (according to a choice of format: single-length, double-length, extended) such that $`a0`$, $`b+0`$, and $`ab`$.

From definitions 4 and 5 and theorem 6.1 we conclude:

* The restriction on the sign of zero bounds in definition 5 is there to make the notation unambiguous. We will see that disambiguating the notation in this way has an advantage for interval division.
* $`\{0\}`$ is written as $`+0,0`$.
* When $`a,b`$ is a floating-point interval, then $`a+\mathrm{}`$ and $`b\mathrm{}`$.

Let us take care to distinguish “real intervals” from “floating-point intervals”. Both are sets of reals. The latter are a subset of the former.

From now on we assume the floating-point intervals as approximation structure when we rely on the fact that for any set $`S`$ of reals there is a unique least floating-point interval $`\varphi (S)`$ containing it.

###### Definition 6

For any real $`x`$, $`x^{}`$ ($`x^+`$) is the left (right) bound of $`\varphi (\{x\})`$.

This operation is implemented by performing a floating-point operation that yields $`x`$ in rounding mode toward $`\mathrm{}`$ ($`+\mathrm{}`$).

## 7 Interval Arithmetic

Much of the standard is concerned with defining, signaling, and trapping exceptions caused by overflow, underflow, and undefined operations. What distinguishes modern interval arithmetic from the old is that no exceptions occur. As we will see, no operation can result in Nan. Every operation is defined on all operands. Moreover, it is defined in such a way that the floating-point endpoints bound the set of the real numbers that are the possible values of the associated variable in the mathematical model.

This property is based on the use of set extensions of the arithmetical operations. It is helped by the use of relational definition rather than functional ones of these operations. We discuss these in turn.

#### Set extensions of functions

Whenever a function $`f`$ is defined on a set $`S`$ and has values in a set $`T`$, there exists the canonical set extension $`\widehat{f}`$, which is a function defined on the subsets of $`S`$ and has as values subsets of $`T`$ according to $`\widehat{f}(X)=\{f(x)xX\}`$ for any $`XS`$. This definition is of interest because it also carries over to partial functions and to multivalued functions.

Though $`X`$ may be an approximation of $`x`$, $`\widehat{f}(X)`$ may not be an element of an approximation structure of $`T`$, so is not necessarily an approximation of $`f(x)`$. But $`\varphi (\widehat{f}(X))`$ does approximate $`f(x)`$. Thus $`\varphi `$ induces a transformation among functions. It changes $`f`$ to the function that maps $`x`$ to $`\varphi (\widehat{f}(\{x\}))`$.

The inverse canonical set extension of $`f`$ is defined as $`f^1(Y)=\{xf(x)Y\}`$. This definition is of interest because such an inverse is defined even when $`f`$ itself has no inverse.

By using the canonical set extensions of a function, one ensures that undefined cases never arise. By considering instead of the arithmetical operations on the reals their canonical set extensions to suitably selected sets of reals (namely, floating-point intervals), undefined cases are eliminated.

An example of a set extension for arithmetical operations is $`X+Y=\{x+yxXyY\}`$. Though $`X`$ and $`Y`$ may be floating-point intervals, that is typically not the case for $`\{x+yxXyY\}`$. So to ensure that addition is closed in the set of floating-point intervals, we need to apply $`\varphi `$, as shown below in the formulas for interval operations that go back to R.E. Moore .

$`X+Y`$ $`=`$ $`\varphi (\{x+yxXyY\})`$
$`XY`$ $`=`$ $`\varphi (\{xyxXyY\})`$
$`XY`$ $`=`$ $`\varphi (\{xyxXyY\})`$
$`X/Y`$ $`=`$ $`\varphi (\{x/yxXyY\})`$
Regarded as a set extension, the above definition of $`X/Y`$ is correct and unambiguous: set extensions are defined just as well for partial functions, functions that are not everywhere defined. Yet many authors have subjected it to the condition $`0Y`$, making it useless in practice. Others have taken a less restrictive stance by changing the definition to:

$$X/Y=\varphi (\{x/yxXyYy0\}).$$

#### Relational definitions

Ratz  has avoided such difficulties by using a relational form of the above definitions. Although not necessary, this relational form also makes it possible to define both addition and subtraction with the same ternary relation $`x+y=z`$. This leads to an attractive uniformity in the definition of the interval arithmetic operations.

###### Definition 7

Let $`X`$ and $`Y`$ be non-empty floating-point intervals. Then interval addition, subtraction, multiplication, and division are defined as follows.

$`X+Y`$ $`\stackrel{\mathrm{def}}{=}`$ $`\varphi (\{zxXyY.x+y=z\})`$
$`XY`$ $`\stackrel{\mathrm{def}}{=}`$ $`\varphi (\{zxXyY.z+y=x\})`$
$`XY`$ $`\stackrel{\mathrm{def}}{=}`$ $`\varphi (\{zxXyY.xy=z\})`$
$`XY`$ $`\stackrel{\mathrm{def}}{=}`$ $`\varphi (\{zxXyY.zy=x\})`$
We use the symbol $``$ in $`XY`$ here for interval division rather than the $`X/Y`$ defined earlier. There is only a difference between the two definitions when $`0,0`$ occurs as an operand. For details, see . The difference is immaterial, as intuition fails in these cases, anyway.

The operations thus defined form an interval arithmetic that is sound in the sense that the resulting sets contain all the real values they should contain according the set extension definition. They are closed in the sense that they are defined for all interval arguments and yield only interval results. Such an interval arithmetic never yields an exception.

It remains to show that these definitions can be efficiently computed by IEEE standard floating-point arithmetic while avoiding the undefined floating-point operations $`\mathrm{}\mathrm{}`$, $`\pm \mathrm{}/\pm \mathrm{}`$, $`0\pm \mathrm{}`$, and $`0/0`$. This we do in the next sections.

### 7.1 The algorithm for interval addition and subtraction

###### Theorem 7.1

If $`X=a,b`$ and $`Y=c,d`$ are non-empty floating-point intervals, then $`X+Y`$ and $`XY`$ according to definition 7 are equal to $`(a+c)^{},(b+d)^+`$ and $`(ad)^{},(bc)^+`$, respectively.

See . The interesting part of the proof takes into account that adding $`a`$ and $`c`$ is undefined if they are infinities with opposite signs. As, according to definition 4, $`a`$ and $`c`$ are not $`+\mathrm{}`$, this cannot happen. Similar reasoning shows that $`b+d`$ is always defined and that the formula for subtraction cannot give an undefined result. Thus, in interval addition and subtraction we achieve the ideal: Never a Nan, and this without the need to test.

### 7.2 The algorithm for interval multiplication

If $`a,b`$ and $`c,d`$ are bounded, real intervals, then

$$a,bc,d=\mathrm{min}(S),\mathrm{max}(S),$$

where $`S=\{ac,ad,bc,bd\}`$.

This formula holds for real rather than floating-point intervals. It is several steps away from interval arithmetic. When we allow the bounds to be any floating-point number, we introduce the possibility that they are infinite. In that case we need to be assured that all four products in $`S`$ are defined. Moreover, we want, as much as possible, to perform only two multiplications, one for each bound. The above formula always requires four.

To attain these goals, we classify the intervals $`a,b`$ and $`c,d`$ according to the signs of their elements, as shown in the table in Figure 1. This classification creates many cases in which intervals can be multiplied with only one multiplication for each bound.

The classification yields four cases (for multiplication the subdivision of $`P`$ and $`N`$ do not matter) for each of the operands, giving at first sight 16 cases. However, when at least one of the operands classifies as $`Z`$, several cases collapse. As a result, we are left with 11 cases.

###### Theorem 7.2

If $`a,b`$ and $`c,d`$ are real intervals, then $`a,bc,d`$ is a real interval whose endpoints are given by the expressions, to be evaluated as extended reals, in Figure 2.

In  the cases indicated as such in the table in Figure 2 are proved directly. The other cases can be proved by symmetry from the case proved already. The symmetries applied are based on the identities $`xy=(xy)`$ or similar ones shown in the last column in the table.

The proofs first show the correctness of the scalar products for bounded real intervals. To allow for floating-point intervals, which can be unbounded, we have to consider whether the products are defined. Let us consider as example the top line according to which $`a,bc,d=ac,bd`$. The undefined cases occur when one operand is 0 and the other $`\mathrm{}`$. It is possible for $`a`$ or $`c`$ to equal 0, but neither can be infinite: because of the classification $`P`$, they cannot be $`\mathrm{}`$; because of their being lower bounds, they cannot be $`+\mathrm{}`$.

Let us now consider $`bd`$. It is possible for $`b`$ or $`d`$ to equal $`+\mathrm{}`$, but neither can be 0 because of the classification $`P`$. One may verify that in every case of the table in Figure 2 undefined values are avoided by a combination of definitions 4 and 5 and the classification of the case concerned.

We need tests to identify the right case in the table anyway to minimize the number of multiplications. We obtain as a bonus the saving of tests to avoid undefined values. Thus, in interval multiplication we achieve the ideal: Never a Nan, and this without the need to test.

### 7.3 Division

For interval multiplication the classification of the interval operands in the classes $`P`$, $`M`$, $`N`$, and $`Z`$ is sufficient. For interval division it turns out that the further subdivision of $`P`$ into $`P_0`$ and $`P_1`$ and of $`N`$ into $`N_0`$ and $`N_1`$ (see the table in figure 1) is relevant for the dividend.

###### Theorem 7.3

If $`a,b`$ and $`c,d`$ are real intervals, then $`a,bc,d`$ is the least floating-point interval containing the real interval whose endpoints are given as the “general formula” column in Figure 3 unless the specified condition in the next column holds, in which case the result is given by the exception case in column 5.

In  the cases indicated as such in the table in Figure 3 are proved directly. The other cases can be proved by symmetry from the case proved already. The symmetries used are based on the identities $`x/y=(x/y)`$ (indicated as $`S_1`$) and $`x/y=(x/y)`$ (indicated by $`S_2`$).

The proofs first show the correctness of the scalar products for bounded real intervals. To allow for floating-point intervals, which can be unbounded, we have to consider whether the products are defined. In the column labelled “unless” we find the values for which an undefined value occurs. In the “exception case” column we find the correct value for the exception case. In every case, evaluating the formula in the third column in IEEE standard floating-point arithmetic in the exception case is defined and gives the infinity of the right sign, as shown in column 5. This property depends on a zero lower bound being $`+0`$ and a zero upper bound being $`0`$, as required by definition 5.

Let us now consider potentially undefined cases. In case of division these are $`\mathrm{}/\mathrm{}`$ and $`0/0`$. Consider for example the top line according to which $`a,bc,d=a/d,b/c\{0\}`$. Because of the classification $`P_1`$, $`a`$ can be neither infinite nor zero. This ensures that $`a/d`$ is defined. Because of the $`P_1`$ classification, $`b`$ cannot be zero. It is possible for $`b`$ to be infinite, but not for $`c`$ because of the $`P`$ classification. This ensures that $`b/c`$ is defined.

One may verify that in every case of the table in Figure 3 undefined values are avoided by a combination of definition 5 and the classification of the case concerned. Thus, in relational interval division we achieve the ideal: Never a Nan, and this without the need to test.

## 8 Related work

For most of the time since the beginning of interval arithmetic, two systems have coexisted. One was the official one, where intervals were bounded, and division by an interval containing zero was undefined. Recognizing the unpracticality of this approach, there was also a definition of “extended” interval arithmetic  where these limitations were lifted. Representative of this state of affairs are the monographs by Hansen  and Kearfott . However, here the specification of interval division is quite far from an efficient implementation that takes advantage of the IEEE floating-point standard. The specification is indirect via multiplication by the interval inverse. There is no consideration of the possibility of undefined operations: presumably one is to perform a test before each operation.

Steps beyond this were taken by Older  in connection with the development of BNR Prolog. A different approach has been taken by Walster , who pioneered the idea that intervals are sets of values rather than abstract elements of an interval algebra. In Walster shares our objective to obtain a closed system of arithmetic without exceptions. He attains this objective in a different way: by including the infinities among the possible values of the variables. In our approach, the variables can only take reals as values; the infinities are only used for the representation of unbounded sets of reals. In this way, the conventional framework of calculus, where variables are restricted to the reals, needs no modification.

## 9 Conclusions

We have presented the result of some recent developments in interval arithmetic that lead to a system with the following properties.

* Correctness The interval operations are such that their result includes all real numbers that are possible as values of the variables according to the mathematical model.
* Freedom of exceptions No floating-point operation needs raise an exception. All divisions by zero are defined and give the correct result: an infinity of the correct sign. This is achieved by a zero lower (upper) bound being $`+0`$ ($`0`$). Mathematically speaking, the system is a closed interval algebra. We do not emphasize the algebra aspect, because it is not important whether it has any interesting properties. Other approaches have limited the applicability of interval arithmetic in their pursuit of a presentable algebra.
* Efficiency The system is efficient in that tests are only needed to determine the right case in the tables in Figures 2 and 3. Tests are not necessary to avoid exceptions.

These properties lead to several observations about the floating-point standard from the point of view of interval arithmetic:

* Exceptions Freedom from exceptions has interesting implications for the standard. A considerable part of the definition effort, and presumably also of the implementation effort, is concerned with defining, signaling, and trapping exceptions caused by overflow, underflow and undefined operations. A processor where the floating-point arithmetic is interval arithmetic can omit this as unnecessary ballast.

Let us review the five exceptions. Invalid Operation is prevented by the design of the algorithms. Division by Zero does occur in our interval arithmetic and is designed to yield the correct result. So it should not be an exception. Overflow occurs in the sense that a real $`x`$ can result in real arithmetic such that $`\varphi (x)`$ is the interval between the greatest finite floating-point number and $`+\mathrm{}`$. This result is mathematically correct and therefore the desired one. There is no reason to terminate computation: it should not be an exception. Underflow means that a lower bound zero is substituted for a nonzero bound with very small absolute value. This is correct and no reason to terminate computation. Inexact result: this might be of some use, but is certainly not essential for interval arithmetic.
* Signed zeros Often signed zeros are regarded as an unavoidable, but regrettable artifact of the sign-magnitude format of floating-point numbers. It is fortunate that the drafters of the standard have nonetheless taken them seriously and defined sensible conventions for operations involving zeros. Especially having the right sign of a zero bound turns out to be useful in interval division.
* Denormalized numbers of view of interval arithmetic, denormalized numbers seem to be neither useful nor harmful. It is different from the point of view of interval constraints. This a method  for using interval arithmetic to solve systems of constraints with real-valued variables. Interval arithmetic is used for the basic operations in constraint propagation. This is an iteration that can be slowed down by denormalized numbers when the limit is zero, even when operations on denormalized numbers are performed at normal speed. Thus the presence of denormalized numbers only plays a role as a performance bug that occurs gratuitously, and fortunately rarely, in this special case.

An argument that is advanced in favour of denormalized numbers is that it justifies compiler optimizations that rely on certain mathematical equivalences that hold only in the presence of denormalized numbers. This is of no interest from the point of view of interval constraints. Any mathematically correct transformation can be performed on the set of constraints without changing the set of solutions obtained by a correctly implemented interval constraint system. This correctness is not dependent on the presence of denormalized numbers. In fact, it only depends on the finite floating-point numbers being some subset $`F`$ of the reals, as described in this paper. Because of this independence, elaborate symbolic processing far beyond currently contemplated compiler optimizations is taken for granted in interval constraints.

## 10 Acknowledgments

Many thanks to Belaid Moa for pointing out errors. We acknowledge generous support from the Natural Science and Engineering Research Council NSERC.