text:
# Extragalactic Gamma-Ray Absorption and the Intrinsic Spectrum of Mkn 501 during the 1997 Flare.

## 1 Introduction

The spectra of high energy $`\gamma `$-rays from extragalactic sources are predicted to be modified by strongly redshift dependent absorption effects caused by interactions of these $`\gamma `$-rays with photons of the intergalactic IR-UV background radiation (Stecker, De Jager & Salamon 1992). Several attempts have been made to infer the IR SED (spectral energy distribution) either from model calculations or observations. (See Hauser & Dwek (2001) for the latest review.) Such information can be used to calculate the optical depth for TeV range photons as a function of energy and redshift so that the intrinsic spectrum of an observed source can be derived.

Of particular interest is the spectrum of Mkn 501 which was observed while strongly flaring in 1997. The spectrum observed at that time by the HEGRA air Cherenkov telescope system (Aharonian, et al. 1999, 2001a) extended to energies greater than 20 TeV, the highest energies yet observed from an extragalactic source.

Application of the Stecker & de Jager (1998) calculations for the optical depth to the observed spectrum as first reported reported by Aharonian, et al. (1999) resulted in an intrinsic spectrum which was consistent with a differential power-law form for the photon spectrum $`E^2`$ between 0.5 and 20 TeV, corresponding to a flat spectral energy distribution (SED), i.e., $`E^2[dN/dE]`$ constant (Konopelko, Kirk, Stecker & Mastichiadas 1999). In this paper, we will rederive the intrinsic spectrum of Mkn 501 using the method employed by Konopelko, et al. (1999) by using the recalibrated observed spectrum as given by Aharonian et al. (2001a) and using a consistent hybrid model for the extragalactic background radiation. We correct for absorption by first recalculating the opacity of intergalactic space to $`\gamma `$-rays, $`\tau (E,z)`$, using the new calculations of the IR background spectrum of Malkan & Stecker (2001) extended into the optical domain as described later in this paper. The consistency between the Whipple telescope and HEGRA spectra of Mkn 501 allows us to derive the intrinsic spectrum of Mkn 501 over two decades of energy.

The intrinsic TeV $`\gamma `$-ray SED which we derive here, $`\nu F(\nu )=E^2[dN/dE]`$, is again quite flat. However, it is slightly convex, with a broad peak in the multi-TeV range. The spectrum can be interpreted as a Compton peak corresponding to a synchrotron peak which has been observed in the X-ray range. This allows us to determine jet parameters which may explain the variability and spectrum of Mkn 501. We discuss the implications of this spectrum in the context of the synchrotoron self-Compton source hypothesis (Stecker, de Jager & Salamon 1996) and for constraints on the breaking of Lorentz invariance.

## 2 A Hybrid SED for the IR-optical Intergalactic Radiation Field

In order to recalculate the $`\gamma `$-ray opacity of intergalactic space from photon-photon pair production interactions, we adopt a hybrid model which uses the near-IR to far-IR background SED calculated by Malkan & Stecker (2001) combined with a reasonable, smooth extention into the optical-UV wavelength range derived from recent HDF (Hubble Deep Field) galaxy counts by Madau & Pozzetti (2001). We will refer to the entire FIR-UV background as the EBL (extragalactic background light).

Malkan & Stecker (2001) have used their empirically based model (Malkan & Stecker 1998) to predict infrared galaxy luminosity functions and deep infrared galaxy counts at various wavelengths. They have also examined their predictions for the IR SED for comparison with the subsequent determinations from the COBE-DIRBE data analysis. Using the assumption of pure luminosity evolution proportional to $`(1+z)^Q`$ out to a redshift of $`z_{flat}`$ and constant (no evolution) for $`z_{flat}<z<z_{max}=4`$, they find that a comparison of their predictions with current ISO galaxy counts at 15 and 175 $`\mu `$m favor their “baseline model” with $`Q=3.1`$ and $`z_{flat}=2`$ (the lower curve in Figure 1). The mid-IR $`\gamma `$-ray upper limits of $``$ 4–5 nW m-2sr-1 (Stecker & De Jager 1997; Stanev & Franceschini 1998; Renault, et al. 2001) also favor $`Q3`$.

On the other hand, the COBE-DIRBE far infrared background flux determinations of Hauser et al. (1998) seem to favor a stronger evolution with $`Q>4`$ up to $`z_{flat}=1`$.111There was a tentative derivation of the background at 60 $`\mu `$m and 100 $`\mu `$m by Finkbeiner, Davis & Schlegel (2001). However, this result now appears to have suffered from contamination by local solar system dust emission (Finkbeiner 2001). The upper curve in Figure 1, which shows the “fast evolution” model of Malkan & Stecker (2001), assumes $`Q=4.1`$ and $`z_{flat}=1.3`$. The lower curve in Figure 1 shows the baseline model.

If one rederives the 140 $`\mu `$m and 240 $`\mu `$m background fluxes from COBE-DIRBE using the COBE-FIRAS calibration, which suffered from smaller systematic errors than the COBE-DIRBE calibration (Fixsen et al. 1997), the fluxes are lowered to the point where they are consistent with both the baseline and fast evolution SEDs of Malkan & Stecker (2001) (see Figure 1). The newly derived flux at 170 $`\mu `$m obtained from ISOPHOT maps by Kiss et al. (2001), also shown in Figure 1, is also consistent with the both Malkan-Stecker SEDs. In this regard, one should also note that the results reported by Hauser et al. (1998) were at the 4$`\sigma `$ level. (See Malkan & Stecker (2001) for further discussion.)

Our hybrid approach extends the IR SEDs of Malkan & Stecker (2001) to shorter wavelengths by exploiting the fact that the Hubble Deep Field (HDF) was able to resolve the lowest intensity galaxies. In other words, a graph of the integral number of galaxies above a given intensity as a function of that intensity shows convergence to a constant value, which implies that the total amount of extragalactic background light (EBL) has been captured in the wavebands of interest. Thus, unless there are significant non-stellar, intergalactic contributions to the optical background or very low surface brightness galaxies missed by the HST (Hubble Space Telescope), the flux derived from the HDF should represent the true EBL.

Galaxy counts provide a strict lower limit on the intensity of the EBL, since faint or unresolved galaxies not included in the counts are additional contributors to the EBL. The relative importance of the unresolved sources compared to that of the resolved sources can be estimated from fluctuation analyses of the EBL. A lack of fluctuations will suggest a certain degree of completeness in the source counts. Based on such fluctuation analyses, the intensities reported by Madau & Pozzetti (2001) can be considered as actual detections of the EBL at the wavelengths studied. Their results are shown as such in Figure 1. These authors also give a lower limit of 15 nW m-2sr-1 for the integrated light between 0.2 and 2.2 $`\mu `$m. 222There is evidence for an upturn in the SED at the shortest observed UV wavelengts near $`2\times 10^{15}`$ Hz as derived from the HST imaging spectrograph (Gardner, Brown and Ferguson 2000), but this feature does not significantly affect our absorption calculations; the UV background only affects $`50`$ to 100 GeV $`\gamma `$-rays from sources at high redshifts (Salamon & Stecker 1998).

The HDF K–band (2.2 $`\mu `$m) flux is further complemented by ground based observations. Gardner (1996) gives a lower limit of 7.4 nW m-2 sr-1 on the EBL intensity, but deeper K-counts (Bershady, Lowenthal & Koo 1998) suggest an increase in the K $`>`$ 23 counts over previous surveys (see e.g. Gardner 1996). The slope in the K-band counts, $`dN/dm`$ is about 0.36. Adopting a slope of 0.4 down to K = 30 magnitudes, E. Dwek (1999, personal communication) derived a conservative upper limit of 13 nW m-2 sr-1 for the EBL intensity in the K-band. These two limits are shown as a thick vertical bar in Figure 1 and provide an important constraint when attempting to connect the far-IR to UV parts of the EBL in any hybrid model. From Figure 1 we also see that the ground based limit on the K-band flux is consistent with the HST galaxy counts in the HDF. Dwek and Arendt (1998) have used COBE-DIRBE data to derive a tentative estimate of the EBL flux at a wavelength of 3.5 $`\mu `$m of 9.9 $`\pm `$ 5.8 (2$`\sigma `$) nW m-2sr-1, shown by the solid diamond point in Figure 1, which is also seen to be in good agreement with our models.

It can be seen from Figure 1 that the Malkan & Stecker (2001) models connect smoothly with the near-IR to UV lower limits for the EBL based on HDF galaxy counts (Madau & Pozzetti 2001). This fact allows us to construct a reasonable hybrid model, smoothly extending the infrared model SEDs of Malkan & Stecker (2001) to UV wavelengths.

The thick curves in Figure 1 show the SEDs predicted by the baseline (lower curve) and fast evolution models (upper curve) of Malkan & Stecker (2001). The baseline SED of Malkan & Stecker (2001) connects smoothly with the mean stellar component (open circles) derived from galaxy counts (HST and ground–based), which is shown as the thin dashed line (lower curve in Figure 1). The integral of this spectrum (extending through the mean UV to K–band fluxes) between 0.2 $`\mu `$m and 2.2 $`\mu `$m gives a value of 16 nW m-2sr-1, which is just above the galaxy count lower limit of 15 nW m-2sr-1 derived by Madau & Pozzetti (2001). The fast evolution model of Malkan & Stecker (2001) also runs smoothly through the $`2\sigma `$ galaxy count limits and gives us a reasonable upper limit for the EBL in the UV to far-IR.

Figure 1 summarizes the observationally derived values for the optical-UV, near-IR and far-IR fluxes which now exist. Unfortunately, foreground emission prevents the direct detection of the EBL in the mid-IR wavelength range (See discussion in Hauser & Dwek 2001.) However, we note that other theoretical models such as those of Tan, Silk & Balland (1999), Rowan-Robinson (2000) and Xu (2000) predict fairly flat SEDs in the mid-infrared range with average flux levels in the 3 to 4 nW m-2sr -1 as do the Malkan & Stecker (2001) models used here. These flux levels are also consistent with the indirect mid-IR constraints indicated by the box in Figure 1 (These constraints are summarized by Stecker (2001). Fine detailed differences between these theoretical SEDs are not important in determining the optical depth of intergalactic space to high energy $`\gamma `$-rays, $`\tau (E,z)`$, (see next section) because its determination involves integrating over a significant range of EBL wavelengths (Stecker, et al. 1992). As an example of this, we note that Vassiliev (2000) has shown that for a flat SED, 90% of the contribution to the $`\gamma `$-ray optical depth is contributed by EBL photons of wavelengths between .28 and 2.7 times the optimum wavelength corresponding to the peak in the pair production cross section. This is roughly an order of magnitude in integrated EBL wavelengths.

We therefore assert that the two Malkan & Stecker curves, extended into the optical-UV range by our hybrid model, give a reasonable representation of the EBL in the UV to far-IR. Other EBL models in the literature whose flux levels roughly fit the present data and have the same spectral characteristics, i.e., a stellar optical peak, a far-IR dust emission peak, and a mid-IR valley which allows for some warm dust emission (see review of Hauser & Dwek 2001), should give similar results on the optical depth of the near Universe $`\tau (E,z)`$ to high energy $`\gamma `$-rays (see next section).

## 3 The Optical Depth for Very High Energy $`\gamma `$-rays from sources at $`z<0.3`$.

We used the prescription of Stecker & de Jager (1998), together with the hybrid model discussed in the previous section, to recalculate the optical depth of intergalactic space $`\tau (E,z)`$ for redshifts between 0.03 (near the $`z`$ values of Mkn 501 and 421) and 0.3, assuming that the EBL is in place by a time corresponding to $`z=0.3`$. Above $`z0.3`$, evolutionary effects have to be taken into account (Salamon & Stecker 1998). Figure 2 shows $`\tau (E,z)`$ calculated using the baseline and fast evolution models for $`\gamma `$-ray energies down to $`50`$ GeV, which is the approximate threshold energy for meaningful image analyses in next generation ground based $`\gamma `$-ray telescopes such as MAGIC, H.E.S.S. and VERITAS. Where our results overlap, they are in good agreement with the metallicity corrected results of Salamon & Stecker (1998). We note that the Salamon & Stecker (1998) results extend to both lower energies and higher redshifts.

As a quantitative example of the uncertainty in the predicted value of $`\tau (E)`$ at the Mkn 501 redshift of 0.031 produced by the uncertainty in the EBL flux, we note that the 19% difference between our two model curves at the mid-IR wavelength of 10 $`\mu `$m (see Figure 1) maps into a difference of 21% in the optical depth at 10 TeV (see Figure 2). This, in turn yields a difference of 34% in the predicted Mkn 501 intrinsic $`\gamma `$-ray flux at 10 TeV. The uncertainty in the energy above which the optical depth is greater than 1 is indicated in Figure 2. We find that $`\tau (E)1`$ for $`E2.3`$ TeV assuming the “fast evolution” model for the EBL and $`\tau (E)1`$ for $`E3.2`$ TeV for the baseline model EBL.

As in Stecker & de Jager (1998), we also obtained parametric expressions for $`\tau (E,z)`$ (as shown in Figure 2), but the accuracy was improved by increasing the order of the polinomials. The expressions are of the form (for $`z<0.3`$ and $`0.1<E_{\mathrm{TeV}}<50`$):

$$log_{10}[\tau (E_{\mathrm{TeV}},z)]\underset{i=0}{\overset{4}{}}a_i(z)(\mathrm{log}_{10}E_{\mathrm{TeV}}+2)^i$$
(1)
where the z-dependent coefficients are given by

$$a_i(z)=\underset{j=0}{\overset{3}{}}a_{ij}(\mathrm{log}_{10}z)^j.$$
(2)
Tables 1 (baseline model) and 2 (fast evolution, numbers in brackets) give the numerical values for $`a_{ij}`$, with $`j=0,1,2,3`$, and $`i=0,1,2,3,4`$. Equation (2) approximates $`\tau (E,z)`$ correctly within 5% for all bounded values of $`z`$ and $`E`$ considered.

## 4 The Observed Spectrum of Mkn 501 during the 1997 Flare.

Recently the HEGRA group (Aharonian et al. 2001a) rederived the energy spectrum of Mkn 501 during the 1997 flare. The validity of the high energy tail of the spectrum out to 20 TeV in their earlier analysis (Aharonian et al. 1999) was questioned because of spillover effects between energy channels. A reanalysis of the events above 3 TeV with better energy resolution gave a flux at $``$ 21 TeV which was a factor of 3.6 lower than that obtained in their previous analysis. It should be stressed that their revised spectrum is consistent with the original spectrum below 20 TeV. The final spectrum between 0.56 and 21 TeV is shown in Figure 3 and marked “OBSERVED”.

Contemporaneous observations of Mkn 501 by the Whipple group (Krennrich et al. 1999) are consistent (within errors) with the HEGRA spectrum in the overlapping energy range between 0.5 and 10 TeV. This allows us to add the two Whipple points below 0.5 TeV to the HEGRA spectrum, giving an observed spectrum between 0.26 and 22 TeV for the Mkn 501 flare. The result is a single spectrum extending over two decades of energy as shown in Figure 3.

## 5 The Intrinsic VHE Spectrum of Mkn 501 during the 1997 Flare.

By using both the Whipple and HEGRA spectra of Mkn 501 and correcting for absorption by multiplying by $`e^{\tau (E)}`$ evaluated at $`z=0.034`$ with the newly derived values for the opacity calculated as discussed in the previous section, we have derived the intrinsic spectrum of Mkn 501 over two decades of energy. This is given by the data points and two curves marked INTRINSIC in Figure 3. The upper of these curves corresponds to the fast evolution case; the lower curve corresponds to the baseline model case.

Fossati et al. (2000) suggested a parameterization to describe smoothly curving blazar spectra. This parameterization is of the form

$$\frac{dN}{dE}=KE^{\mathrm{\Gamma }_1}\left(1+(\frac{E}{E_B})^f\right)^{(\mathrm{\Gamma }_1\mathrm{\Gamma }_2)/f}$$
(3)
A spectrum of this form changes gradually from a spectral index of $`\mathrm{\Gamma }_1`$ to an index of $`\mathrm{\Gamma }_2`$ when the energy $`E`$ increases through the break energy $`E_B`$. The parameter $`f`$ describes the rapidity (“fastness”) of the change in spectral index over energy. This model was used by Fossati et al. (2000) to describe the curvature in the Mkn 421 spectrum. We have fitted our derived intrinsic Mkn 501 spectrum to the form given by equation (3).

We applied the formalism of Fossati et al. (2000) and found the best-fit parameters $`K`$, $`E_B`$, $`f`$, $`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_2`$ after correcting the observed spectrum for intergalactic absorption. A total of 27 data points between 0.26 ad 21 TeV were used for each of the two intrinsic spectra in Figure 3, and we fitted to the Fossati et al. parameterization of these spectra, keeping all parameters free. Whereas the low energy spectral index $`\mathrm{\Gamma }_1`$ appears to be well constrained, we find that the higher energy index $`\mathrm{\Gamma }_2`$ is unconstrained. Table 3 shows the parameter values for each fit, assuming two arbitrary values for $`\mathrm{\Gamma }_2`$ (2.5 and 3.0). The parameter $`E_M`$ is derived from the fit, and corresponds to the $`\gamma `$-ray energy where the SED peaks in $`\nu F_\nu `$. The fits are acceptable in terms of the $`\chi ^2`$ per degree of freedom.

It is important to note that the intrinsic spectrum which we derive does not have an upturn near 20 TeV as obtained previously by some authors (e,g. Protheroe & Meyer 2000; Dwek 2001) working with the earlier analysis of Aharonian, et al. (1999) as well as the 60 $`\mu `$m flux earlier reported by Finkbeiner et al. (2000 - see footnote 1). The upturn obtained by Dwek (2001) can be effectively eliminated by using the revised HEGRA $`\gamma `$-ray data with the observed flux reduced by a factor of 4 at $``$20 TeV; the stronger upturn obtained by Protheroe & Meyer (2001) is eliminated both by using the revised HEGRA data and a more realistic IR-SED than the one these authors used which had no mid-IR valley and an unreasonably high mid-IR flux.

Aharonian, Timokhin & Plyasheshnikov (2001b) have also found an upturn in their analysis, particularly if they use an SED constructed using the flux claimed by Lagache et al. (2000) at 100 $`\mu `$m, the only point in Figure 1 which is inconsistent with our SEDs at the 2$`\sigma `$ level. Assuming that this upturn is correct, these authors suggested avoiding a “ TeV-IR crisis” by invoking a cold jet near the central engine with a bulk Lorentz factor greater than $`10^7`$, with the upturn resulting from comptonization.

In this context, we note that the COBE-DIRBE group has argued that a real flux derived from the COBE-DIRBE data at 100 $`\mu `$m, as claimed by Lagache et al. (2000), is untenable because isotropy in the residuals (after foreground subtractions) could not be proven. Dwek et al. (1998) have concluded that only a conservative lower limit of 5 m-2sr-1 could be inferred at 100 $`\mu `$m. For Model II of Aharonian et al. (2001b) which neglects the 100 $`\mu `$m point, these authors derive a source spectrum which is quite flat up to an energy of $``$ 20 TeV, as we do in this work. This obviates any need to invoke extreme bulk Lorentz factors in modelling the Mkn 501 source spectrum.

Our intrinsic SED is slightly convex, rather than a single power law, consistent with an SSC hypotheses for the origin of the TeV radiation. We find that this SED peaks at $`E_M89`$ TeV (independent of the unconstrained $`\mathrm{\Gamma }_2`$) in the case where the fast evolution EBL is assumed and that $`E_M5`$ TeV if the baseline EBL is assumed (see Figure 3 and Table 3.)

## 6 A Simple Interpretation of the Intrinsic Mkn 501 Flare Spectrum within the Context of the SSC Paradigm

If one observationally determines a flaring variability time $`\mathrm{\Delta }`$t, the optical depth for $`\gamma `$$`\gamma `$ absorption from pair production is then given by (Dermer & Schlickheiser 1994; Aharonian, et al. 1999)

$$\tau _{\gamma \gamma }0.055\varphi _{10}\mathrm{\Delta }t_{\mathrm{day}}^1\delta _{10}^6H_{65}^2E_{\mathrm{TeV}}$$
(4)
using the normalizations $`\delta _{10}\delta /10`$ and $`\varphi _{10}\varphi _\mathrm{r}/10^{10}`$ erg cm-2s-1 with $`\varphi (ϵ)ϵ^2(dN/dϵ)`$ being the observed energy flux of the low energy radiation. The flux $`\varphi (ϵ)`$ is determined at $`ϵ100\delta _{10}^2E_{\mathrm{TeV}}^1`$ eV which is the energy range of the low energy photons responsible for the intrinsic optical depth of high energy $`\gamma `$-rays in the jet. The value of the relevant low energy radiation flux $`\varphi _{10}`$ is approximately 0.5 in the optical–UV range. $`H_{65}`$ is the Hubble constant is in units of 65 km s-1 Mpc-1.

The time variation of the TeV $`\gamma `$-ray flux scales as the square of the time variation of the X-ray flux. This scaling favours a synchrotron-self Compton origin for the Mkn 501 source with the synchrotron peak in the X-ray range and the Compton peak in the TeV $`\gamma `$-ray range (Krawczynski, et al. 2000). The X-ray SED of Mkn 501 appears to have a maximum in the $``$ 50-100 keV range in the synchrotron component (Catanese et al. 1997; Pian et al. 1998; Petry, et al. 2000). This shape should be reflected in the high energy $`\gamma `$-ray domain, but with different fastness ($`f`$) and slopes ($`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_2`$) owing to the energy dependence of the cross section for inverse Compton scattering as the possibility of having more than one soft target photon component for such scattering (Stecker, De Jager & Salamon 1996).

Since Mkn 501 is a giant elliptical galaxy with little dust, it is also reasonable to assume that the galaxy itself does not produce enough infrared radiation to provide a significant opacity to high energy $`\gamma `$-rays. Also, if our EBL is approximately correct, it is reasonable to assume that the dominant absorption process is intergalactic and that pair-production in the jet in negligible. This hypothesis is supported by the fact that the high energy $`\gamma `$-ray SED did not steepen during the flare. This implies that the optical depth given by eq. (4) is less than unity out to the highest observed energy $`E20`$ TeV. The pair-production opacity constraints then limit the Doppler factor to

$$\delta _{10}1.2\left(\frac{\varphi _{10}}{0.5}\right)^{1/6}\left(\frac{5\mathrm{hours}}{\mathrm{\Delta }t}\right)^{1/6}H_{65}^{1/3}\left(\frac{E_\gamma }{20\mathrm{TeV}}\right)^{1/6}$$
(5)
which is consistent with the lower limit of $`\delta _{10}=1`$ set by Aharonian et al. (1999).

The energy of the synchrotron peak in the X-ray region is given by

$$E_{X,M}(\mathrm{keV})67B_{0.1}\delta _{10}[E_e(\mathrm{TeV})]^2$$
(6)
where $`E_e`$ is the energy of the electrons which are radiating into the synchrotron peak and $`B_{0.1}`$ is in units of 0.1 G. The energy of the Compton peak in the Klein-Nishina range is given by

$$E_M=\delta \kappa E_e$$
(7)
where the inelasticity factor, $`\kappa 1`$ in the Klein-Nishina limit.

From equations (4) and (5), we obtain the SSC expression for $`E_M`$,

$$E_M(\mathrm{TeV})1.2\kappa \delta _{10}^{1/2}B_{0.1}^{1/2}E_{X,M}(\mathrm{keV})^{1/2}.$$
(8)
Note that this expression also holds if SEC is important. Krawczynski et al. (2000) have estimated the $`B`$ field from the observed lag in the cooling time between the 3 keV and 12 keV synchrotron radiation in the jet to be $``$ 0.025 G.

Taking $`\delta _{10}1`$, $`B_{0.1}=1`$, $`\kappa 1`$ and 50 keV $`E_{X,M}`$ 100 keV (Pian, et al. 1998) then eq. (8) gives a prediction of $`E_{TeV,M}`$ 9 – 12 TeV, in good agreement with the results of the $`\chi ^2`$ fit to our intrinsic SED for the Mkn 501 flare of $``$ 5 – 9 TeV. (Note that this value only depends on the square roots of the parameters involved.)

It should be noted that Tavecchio, et al. (1998) neglected intergalactic absorption and used the observed spectrum of Mkn 501 with an assumed Compton peak at sub-TeV energies to determine constraints on SSC jet parameters. The fact that they failed to derive a self-consistent model shows the importance of taking intergalactic absorption into account in theoretical work involving TeV blazars.

Our interpretation is based on time averages of flare spectra measured during the 1997 high state. However, our analysis is validated by the analyses of Petry et al. (2000). These authors also compared the TeV data with the soft- and hard X-ray data after averaging the X-ray data over several flares during the same high state. Our comparison uses both the soft X-ray spectral index and the hard X-ray peak taken from the Petry et al. (2000) averages.

## 7 Conclusions

The IR SEDs which Malkan & Stecker (2001) have derived and argued are consistent with all of the reliable data on the IR EBL connect smoothly with the EBL derived from converging galaxy counts in the near-IR to UV. Furthermore, at the longest wavelengths, they also connect smoothly with the COBE-FIRAS results. This allows us to construct hybrid EBL SEDs which were used to recalculate the optical depth versus redshift for high energy $`\gamma `$-rays between 50 GeV and 100 TeV, the energy above which the Universe becomes highly opaque, even at low redshifts, owing to interactions with photons of the 2.7K cosmic background radiation (Stecker 1969).

We find that our derived intrinsic flaring SED for Mkn 501 exhibits a broad maximum in the $``$ 5-10 TeV range, which is consistent with a value of $``$ 10 TeV predicted from several consistent physical arguments and independent observational constraints for the SSC model (see previous section). The Mkn 501 SED derived in this paper does not show the upturn above $``$ 20 TeV found in some other studies. We agrue here that such a feature is an artifact produced partly by spillover in the last data bin in the original HEGRA analysis and partly by assuming too high an IR background SED.

Our differential photon spectral index, $`\mathrm{\Gamma }_1`$, for the energy range $`E_\gamma E_\mathrm{M}`$ is 1.6 to 1.7 (Table 3), which is close to the differential photon spectral index of 1.76 at soft X-ray energies (Petry et al. 2000). The similarity of these spectral indices implies that $`\gamma `$-rays below $`1`$ TeV are produced in the Thomson regime by scattering off soft photons with energies in the optical-IR range. On the other hand, $`\gamma `$-rays near the $``$ 7 $`\pm `$ 2 TeV Compton peak are the result of scattering in the Klein-Nishina range (see previous section).

Our derived intrinsic Mkn 501 SED is quite flat in the multi-TeV range as shown in Figure 3. This is in marked contrast to the dramatic turnover in its observed SED. We argue that this is strong evidence that the observed spectrum shows just the absorption effect predicted. Among other consequences of this conclusion, the observed absorption feature allows one to put strong constraints on the breaking of Lorentz invariance (Stecker & Glashow 2001).

###### Acknowledgements.

Acknowledgment: We would like to thank Felix Aharonian, Eli Dwek and Karl Mannheim for useful discussions.

## Figure Captions

Figure 1. The SED of the EBL (see text for references and descriptions). All error bars given at the $`\pm 2\sigma `$ level. The Malkan & Stecker (2001) fast evolution model is shown by the upper, thick solid curve between $`\mathrm{log}_{10}\nu =11.8`$ to 13.8) and their baseline model is shown by the thick dashed line over the same frequency range. Convergent Hubble Deep Field galaxy counts (Madau & Pozzetti 2001): open circles; Ground–based galaxy counts limits at $`2.2\mu `$m (see text): thick vertical bar marked “$`2.2\mu `$”; COBE-DIRBE photometric sky residuals (Wright & Reese 2000): small open squares; TeV $`\gamma `$-ray-based upper limits: thick box marked “TeV” (see text for references); ISOCAM lower limit at $`15\mu `$m: upward arrow marked $`15\mu `$ (Elbaz et al. 1999); COBE-DIRBE far-IR sky residuals (Hauser et al. 1998): large open squares at 140 $`\mu `$m and 240 $`\mu `$m; COBE-DIRBE data recalibrated using the COBE- FIRAS calibration (see text): large solid squares without error bars; ISOCAM 170 $`\mu `$m flux (Kiss et al. 2001): solid circle; COBE-FIRAS far-IR sky residuals (Fixsen, et al. 1998): thick dot-dash curve with $`95\%`$ confidence band (thin dot–dash band); $`100\mu `$m COBE-DIRBE point (Lagache et al. 2000): small solid square; flux at 3.5 $`\mu `$m from COBE-DIRBE (Dwek & Arendt 1998): solid diamond.

Figure 2. The optical depth for $`\gamma `$-rays above 50 GeV given for redshifts between 0.03 and 0.3 (as labelled) calculated using the medium (dashed lines) and fast evolution (solid lines) SED of Malkan & Stecker (2001).

Figure 3. The observed spectrum and derived intrinsic spectrum of Mkn 501. The observed spectral data are as measured by HEGRA (solid circles) and Whipple (solid squares). The upper points are the absorption corrected data points (marked “INTRINSIC”) using our fast evolution hybrid EBL (upper data set and solid curve fit) and baseline hybrid EBL (lower data set with dashed curve fit). The fit parameters are given in Table 1.


format: text

text:
# Thermal Comptonization and Disk Thermal Reprocessing in NGC 3516

## 1 Introduction

A ubiquitous feature of the optical/UV emission from Type 1 Seyfert galaxies (Sy1s) is the highly correlated variability of the continuum flux across a wide range of wavelengths (e.g., Peterson et al. 1998). When this behavior was first observed during the 1989–1990 broad line region reverberation mapping campaign for NGC 5548 (Peterson et al. 1991; Clavel et al. 1991), it was immediately recognized that the upper limits for the relative lags between the optical and UV wavebands were too small for the variability to have been mediated by hydrodynamic processes in the putative accretion disk — the differences in the characteristic radii for the emission in the respective wavebands were too large. In fact, the upper limits for the interband lags of $`2`$ days required signal speeds of $`0.1c`$ which could not have occurred within the disk itself. This fact and the observation that the optical/UV continuum became bluer as the flux increased led to the conclusion that the continuum variations were due to thermal reprocessing in the disk of higher energy radiation produced by a single varying source, such as the X-ray source posited to exist at the disk’s inner regions, near the central black hole (Krolik et al. 1991; Courvoisier & Clavel 1991; Collin-Souffrin 1991). Subsequent observations of other Sy1s showed that highly correlated broad band optical/UV variability is generic to these objects (e.g., Courvoisier & Clavel 1991).

Recently, the thermal reprocessing picture has been challenged by simultaneous optical, UV, and X-ray observations of Sy1s. The first object so observed was NGC 7469 by the IUE and RXTE satellites over $`30`$ days in 1997 (Nandra et al. 1998). Although the UV and X-ray light curves showed some similar variability characteristics both in amplitude and time scales, these light curves were not nearly as well correlated as would be expected if the UV emission was due to thermal reprocessing. Naively, the variations in the UV should follow those in the X-rays with some characteristic time delay, but the maxima in the UV light curves preceded similarly shaped maxima in the X-ray light curve by $`4`$ days, while minima in both light curves occurred nearly simultaneously.

The 1998 observations of NGC 3516 by HST, ASCA, and RXTE seem to be equally troubling for the thermal reprocessing model (Edelson et al. 2000). These observations consisted of intensive monitoring in optical, UV, and X-ray bands over 3 days and so were able to probe correlations on far shorter time scales than any previous set of observations. The X-rays showed very strong variations ($`60`$% peak-to-peak), while the changes in the optical continuum were much smaller ($`3`$%). Because of the accurate relative photometry afforded by the HST/STIS CCDs, the measurement uncertainties were sufficiently small to reveal significant variability that was strongly correlated across the optical bands. Hence, the original motivation for the thermal reprocessing model was still present, but as for NGC 7469, the X-ray and optical flux variations were not strongly correlated on any relevant time scale.

More recently, Nandra et al. (2000) have re-examined the 1997 observations of NGC 7469 and found that while the relatively narrow band 2–10 keV X-ray flux wasn’t well-correlated with the UV flux the spectral index of the X-ray emission was. This led them to suggest that the UV emission is actually correlated with the bolometric X-ray flux. In this Letter, we expand upon that idea and describe what is required for the multiwaveband phenomenology exhibited by NGC 3516 to be due to disk thermal reprocessing in the context of a thermal Comptonization model for the X-ray emission.

## 2 X-ray Thermal Reprocessing in the Sphere+Disk Geometry

Most of the important details of the model are described in Zdziarski, Lubiński, & Smith (1999; hereafter ZLS) which the interested reader is encouraged to consult. We will proceed by summarizing the significant features and results from that work. This model was proposed by ZLS to explain the apparent correlation between the X-ray spectral index, $`\mathrm{\Gamma }`$, and the relative magnitude, $`R`$, of the so-called Compton reflection hump seen in the hard X-ray spectra of Sy1s and X-ray binaries. ZLS argued that this correlation results from the reprocessing and Comptonizing interactions between the disk and the hot plasma believed to produce the X-ray continuum. Expanding upon the work of previous authors (e.g., Poutanen, Krolik, & Ryde 1997; Zdziarski et al. 1998), ZLS noted that the $`R`$-$`\mathrm{\Gamma }`$ correlation could arise in a geometry in which the Comptonizing plasma occupies a spherical region of radius $`r_s`$ centered on the black hole. The accretion disk has an inner radius $`r_{\mathrm{min}}`$ which penetrates this sphere at its equator by varying amounts. This varying penetration could be due to changes in $`r_s`$, $`r_{\mathrm{min}}`$, or both. Smaller values of $`r_{\mathrm{min}}/r_s`$ result in a larger fraction of the X-ray photons from the sphere being intercepted by the disk. The covering factor, $`R`$, for Compton reflection and thermal reprocessing is larger, and a larger fraction of the disk soft photons re-enter the plasma precipitating additional cooling via Compton losses and thereby reducing the temperature of the plasma. This results in a softer X-ray spectrum and a larger value for the photon index $`\mathrm{\Gamma }`$. Enforcing energy balance under the assumption that the disk emits only reprocessed radiation, ZLS derived a $`R`$-$`\mathrm{\Gamma }`$ relation which looks very similar to the observed correlations.

We follow ZLS by assuming that the region of Comptonizing plasma is homogeneous and optically thin and that we can neglect internal viscous dissipation in the disk, at least initially. Reprocessing of incident X-rays by a razor-thin disk then gives a radial effective temperature distribution

$$T(r)=\left[\frac{3}{16\pi ^2\sigma _B}\frac{L_x}{r_s^2}h\left(\frac{r}{r_s}\right)(1a)\right]^{1/4}.$$
(1)
Here $`L_x`$ is the total luminosity of the Comptonizing plasma, $`\sigma _B`$ is the Stefan-Boltzmann constant, and $`a`$ is the disk albedo. The dimensionless function $`h(r/r_s)`$ is given by Eq. 1 of ZLS.

Assuming local blackbody emission and that the optical portion of the observed spectrum originates at large radii, we may use the asymptotic form $`h(r/r_s1)(\pi /4)(r_s/r)^3`$ in Eq. 1 to obtain the disk spectrum

$$F_\nu \nu ^3_{r_{\mathrm{min}}}^{\mathrm{}}r𝑑r\frac{1}{\mathrm{exp}\left(E_\nu r^{3/4}/k_B𝒜\right)1},$$
(2)
where $`T=𝒜r^{3/4}`$, $`𝒜[3L_xr_s(1a)/64\pi \sigma _B]^{1/4}`$, and $`E_\nu `$ is the photon energy. Neglecting variations in the disk albedo and assuming that $`r_{\mathrm{min}}`$ is fixed on the time scales of interest, a constant optical flux implies a constant value of $`𝒜`$. In the presence of variations of the X-ray luminosity, this requires $`L_xr_s^1`$.

For $`rr_s`$, i.e., at higher frequencies, it turns out that the large radius limit of $`h(r/r_s)`$ is a surprisingly good approximation for evaluating the disk spectrum since it enters as the 1/4-power in $`T`$, and at $`r/r_s=1`$, the function takes the value $`h(1)=4/3`$ (ZLS), and we have $`(\pi /4)^{1/4}(4/3)^{1/4}1`$. Hence, for $`r_sr_{\mathrm{min}}`$, adopting the asymptotic form for $`h(r/r_s)`$ throughout and setting $`L_xr_s^1`$, all $`r_s`$-dependence drops out of Eq. 2, and in this approximation, the resulting disk spectrum has a constant shape and magnitude for different values of $`r_s`$. The exact expression for $`h(r/r_s)`$ flattens for $`r<r_s`$, so when $`r_{\mathrm{min}}<r_s`$, the disk spectrum predicted by the full calculation will be somewhat softer, but not appreciably at the wavelengths which have actually been observed.

Given the X-ray luminosity $`L_x`$, we estimate the X-ray spectrum produced by thermal Comptonization in the sphere. For this we use a second result from ZLS which relates the amplification factor for Comptonization, $`A`$, to the ratio $`r_{\mathrm{min}}/r_s`$ (ZLS, Eq. 2). Following ZLS, we use the empirical result of Beloborodov (1999a) which gives the X-ray spectral index as a function of $`A`$: $`\mathrm{\Gamma }2.33(1A)^{1/10}`$. To determine the overall magnitude of the X-ray spectrum, we use another empirical result of Beloborodov (1999a) to estimate the temperature of the Comptonizing plasma, $`\mathrm{\Gamma }(9/4)y^{2/9}`$, where the Compton parameter is given by $`y=4(\theta _e+4\theta _e^2)\tau (\tau +1)`$ and $`\theta _e=k_BT_e/m_ec^2`$. We assume a typical value for the Thomson depth, $`\tau =1`$, and solve the preceding equations to find the electron temperature $`T_e`$. We model the X-ray continuum as a power-law with an exponential roll-over, $`L_E(E/E_c)^{1\mathrm{\Gamma }}\mathrm{exp}(E/E_c)`$, where the cut-off energy is taken to be $`E_c=2k_BT_e`$ appropriate for a Thomson depth $`\tau 1`$. The normalization, $`L_0`$, of the X-ray spectrum is given by

$$L_x=L_0_{E_{\mathrm{min}}/E_c}^{\mathrm{}}(E/E_c)^{1\mathrm{\Gamma }}\mathrm{exp}(E/E_c)d(E/E_c).$$
(3)
For $`\mathrm{\Gamma }<2`$, the precise value of the minimum photon energy is not very important so we have set it at a nominal value of $`E_{\mathrm{min}}=0.01`$ keV.

## 3 Comparisons to Observations

We have computed thermal Comptonization and disk reprocessing spectra in this model for NGC 3516. The X-ray luminosity $`L_x`$ is found by fitting the disk spectrum to the mean optical fluxes measured by HST (Edelson et al. 2000). We fix the disk inner radius at $`r_{\mathrm{min}}=3\times 10^{13}`$ cm and plot the resulting spectral energy distributions (SEDs) in Figure Thermal Comptonization and Disk Thermal Reprocessing in NGC 3516 for values of the sphere radius in the range $`r_s=1.5`$$`5\times 10^{13}`$ cm (solid curves). As claimed, the calculated disk spectrum is strikingly constant, with only relatively small variations in the bluer parts of the spectrum. For comparison, we plot the mean fluxes measured by HST. Given the crudeness of this model, the quality of the fit is satisfactory even for the UV flux.111This datum was not included in the fit for $`L_x`$ since the UV flux was only monitored over a significantly shorter period than the optical data due to SAA passages. In addition, the UV data point, at $`1350`$Å, is in a part of the spectrum which may be contaminated by the broad wings of the Siiv emission line (see Edelson et al. 2000, Fig. 3). Unfortunately, similar contamination problems may exist for the optical data as well (see Edelson et al. 2000, Fig. 4).

Since we have thus far neglected local viscous dissipation in the disk, the black hole mass does not enter into the calculation of the disk emission. However, the value of the disk inner radius can, in principle, be constrained by the optical and UV data. The dashed curve in Figure Thermal Comptonization and Disk Thermal Reprocessing in NGC 3516 is the SED obtained for $`r_{\mathrm{min}}=r_s=6\times 10^{13}`$ cm. This shows that while the disk spectrum is relatively insensitive to factor $`3`$ changes in $`r_s`$ for a given optical flux, comparably sized changes in $`r_{\mathrm{min}}`$ produce strong color variations in the UV. Therefore, in order to account for the UV flux, the disk inner radius cannot be too much larger than $`r_{\mathrm{min}}=3\times 10^{13}`$ cm. For a central black hole mass of $`3\times 10^7`$ M(e.g., Edelson & Nandra 1999), this value for the disk inner radius is roughly consistent with innermost radius implied by ASCA observations of the broad Fe K$`\alpha `$ emission line (Nandra et al. 1999).

Despite the fact that the optical/UV continuum barely changes as a function of $`r_s`$, the X-ray spectrum varies drastically. In the 2–10 keV band in particular, the flux changes by factors of $`2`$ for the spectra shown. To relate this result more concretely to the X-ray observations of NGC 3516, we plot in Figure Thermal Comptonization and Disk Thermal Reprocessing in NGC 3516 the X-ray spectral index versus 2–10 keV flux for a family of curves corresponding to different optical fluxes with the position on any given curve being parameterized by the value of $`r_s`$. We propose that the HST/ASCA/RXTE monitoring observations of NGC 3516 can be explained in this model by changes in $`r_s`$ and accompanying changes in $`L_x`$ so that $`L_xr_s^1`$ is approximately maintained, and the source takes a path in the $`\mathrm{\Gamma }`$-$`F_{210}`$ plane which lies nearly along a curve of constant optical flux. Any significant variations in the optical/UV continuum will correspond to deviations from the constant flux curve and therefore from the $`L_xr_s^1`$ relation. Furthermore, the changes associated with such departures will be broad band and will show the same color variations that have been previously associated with thermal reprocessing.

For comparison, we have plotted in Figure Thermal Comptonization and Disk Thermal Reprocessing in NGC 3516 values of $`\mathrm{\Gamma }`$ and $`F_{210}`$ which we have obtained from fitting the publicly archived RXTE NGC 3516 data for the observations in question. Although this simple model overpredicts the X-ray flux by a factor $`4`$ for the observed range of spectral indices, the correlation of spectral index and flux given by the model is strikingly similar to that of the data. The flux discrepancy could be reduced somewhat by selecting a smaller value for the Thomson depth $`\tau `$ since for a given spectral index (and $`y`$-parameter) that will imply a larger value of $`T_e`$ and more of the X-ray luminosity can be “hidden” at higher energies. Alternatively, a large reduction can be obtained by including contributions to the disk emission from local viscous dissipation. This corresponds to adding a term $`(GM\dot{M}/r^3)(1(r_I/r)^{1/2})`$ inside the square brackets in the expression for the disk temperature (Eq. 1) where $`r_I=r_g`$ for a maximal Kerr metric and $`r_I=6r_g`$ for a Schwarzschild metric. Since the viscous dissipation term is independent of $`r_s`$, it acts essentially as a constant offset for the required value of the X-ray luminosity, and $`L_xr_s^1`$ is still needed to maintain a constant optical flux. In order to match the observed X-ray fluxes, roughly 1/4 of the total disk luminosity must be due to thermal reprocessing with the remaining 3/4 being due to viscous dissipation. Assuming $`M=3\times 10^7`$ M and $`r_I=r_g`$, we have solved for the energy balance condition and find that an accretion rate of $`\dot{M}3.6\times 10^3`$ M yr-1 is required to match the mean X-ray flux. This solution is plotted as the thick solid curve in Figure Thermal Comptonization and Disk Thermal Reprocessing in NGC 3516.

Although the energy balance condition we have used only restricts the properties of the Comptonizing plasma, global energy balance is roughly satisfied by this solution. The bolometric X-ray luminosity which gives the observed mean X-ray flux and spectral index is $`L_x7.0\times 10^{43}`$ erg s-1, while the above accretion rate implies a disk luminosity of $`L_{\mathrm{disk}}3.3\times 10^{43}`$ erg s-1 due solely to viscous dissipation. For a disk extending all the way to $`r_I=r_g`$ (rather than being truncated at $`r_s`$), the total implied luminosity is $`L_{\mathrm{tot}}1.0\times 10^{44}`$erg s-1, and we have

$$L_{\mathrm{tot}}L_x+L_{\mathrm{disk}}.$$
(4)
This would be expected if the accretion luminosity that would otherwise emerge as local black body emission from the inner regions of the untruncated disk were instead partitioned to the heating of the Comptonizing plasma. We note that we have arbitrarily set the Thomson depth to $`\tau =1`$ in computing this solution. By fine tuning the values of $`\tau `$ and $`\dot{M}`$, we find that we can satisfy the global energy balance condition arbitrarily well while simultaneously matching the observed mean X-ray flux, spectral index and optical continuum levels.

## 4 Discussion and Conclusions

The fact that an approximately constant optical/UV continuum can in principle result from thermal reprocessing of a highly variable X-ray flux is connected with the fact that the radial distribution of the flux incident upon the disk is well approximated by a power-law, even very close to the innermost disk radius. As we showed above, the constant flux constraint, in conjunction with the assumed constancy of the disk inner radius and disk albedo, causes the radius of the Comptonizing sphere to drop out of the integral describing the disk spectrum. However, this constraint does not demand that the shape of the incident flux distribution be a power-law. The only requirement is that the relevant physical parameters affecting the incident flux can be considered independently from the expression which describes the shape of the distribution. For example, consider the model proposed by Beloborodov (1999b) in which the Comptonizing medium consists of a blob of plasma which is at some height $`r_h`$ above the disk. The blob moves with mildly relativistic velocities either vertically away from or towards the underlying accretion disk, and the effective reflection fraction is controlled by relativistic aberration of the flux from the blob. Using the formulae in Beloborodov (1999b), the radial distribution of the flux incident upon the disk is given by

$$F_{\mathrm{inc}}\frac{L_x}{r_h^2\gamma ^4}\left[\frac{(\stackrel{~}{r}^2+1)^{1/2}}{\beta +(\stackrel{~}{r}^2+1)^2}\right],$$
(5)
where $`\beta `$ is the blob velocity in units of $`c`$, $`\gamma =(1\beta ^2)^{1/2}`$, and $`\stackrel{~}{r}=r/r_h`$ is the disk radius in units of the blob height. For the range of blob velocities required by this model to reproduce the observed $`R`$-$`\mathrm{\Gamma }`$ correlation, $`\beta 0.3`$–0.7, the factor in the square brackets which gives the shape of the flux distribution is largely insensitive to the exact value of $`\beta `$. Therefore, since the lower limit of the radial flux integral can be taken to be zero (cf. Eq. 2), a nearly constant optical flux in this model can be obtained if $`L_xr_h^2\gamma ^4`$.

Although neither of the aforementioned models offer any reason for the luminosity scalings which seem to be required by the NGC 3516 observations, any plausible explanation would still have to contend with why NGC 7469 and probably other Sy1s violate it. Therefore, we will not attempt to posit a physically motivated rationale for the specific relationship between X-ray luminosity and geometry required for any particular object, but a cursory examination of the important time scales may provide some insights. For NGC 3516, the typical variability time scale exhibited by the 2–10 keV light curve is $`20`$–30 ks (Edelson et al. 2000). Assuming a radial Thomson depth $`\tau =1`$ and $`r_sr_{\mathrm{min}}`$, the cooling time is very short: $`t_{\mathrm{cool}}N_ek_BT_e/L_x10`$ s, where $`N_e=(4\pi r_s^3/3)(\tau /r_s\sigma _T)`$ is the total number of electrons in the sphere with mean temperature $`T_e10^9`$K. Unfortunately, the thermal expansion time for a sphere of this size is fairly long, approximately $`t_{\mathrm{exp}}100`$ ks. This assumes that the protons and the electrons are at the same temperature. If we have a two-temperature plasma in which the protons possess a virial temperature of $`10^{12}`$ K, then the expansion time scale is reduced substantially. Alternatively, this could also be achieved if we have a pair plasma instead.

We have assumed for simplicity that the thermal emission from the disk is locally blackbody. Detailed modeling of bare disk atmospheres appropriate for Sy1s indicates that this is a moderately reasonable approximation in the optical, but significant deviations can occur at higher frequencies (e.g., Hubeny et al. 2000, 2001). This appears to remain true for disks powered exclusively by external X-ray illumination (Sincell & Krolik 1997). Our treatment of the reprocessed flux in the optical is therefore probably fairly good, but the full seed photon spectrum illuminating the Comptonizing plasma will not be very accurate. More detailed modeling of the disk spectrum will likely affect the specific values we deduce for the disk inner radius, X-ray luminosity, and Comptonizing sphere radius, but it will not change our general conclusions.

The existing simultaneous optical, UV, and X-ray data plus the global energy balance condition (Eq. 4) provide just enough information to infer unique values for the five model parameters, $`r_{\mathrm{min}}`$, $`r_s`$, $`T_e`$, $`\tau `$, and $`\dot{M}`$, required to describe the entire spectral energy distribution. Therefore, a crucial test of this general scheme for linking thermal Comptonization and thermal reprocessing will be provided by simultaneous optical, UV, and broad band X-ray monitoring observations which include reliable measurements of the thermal roll-over at energies $`100`$ keV. These latter measurements are necessary for making estimates of the X-ray luminosity $`L_x`$ which can then be compared with the luminosities required by the model. Forthcoming INTEGRAL observations of Sy1s will provide some of the necessary spectral coverage to address this issue, but what is really required is a hard X-ray/soft gamma-ray telescope which can measure the X-ray spectrum and thermal roll-over on time scales comparable to the shortest variability time scales seen from these objects.

We thank Andzrej Zdziarski and the anonymous referee for helpful comments. JC was partially supported by NASA ATP grant NAG 5-7723 during the course of this work. OB acknowledges support from NASA grant NAG 5-7075 and NSF grant PHY 99-07949. This research has made use of data obtained from the High Energy Astrophysics Science Archive Research Center (HEASARC), provided by NASA’s Goddard Space Flight Center.

Model SEDs which have been fit to the HST optical fluxes for NGC 3516. For the solid curves, the disk inner radius is set at $`r_{\mathrm{min}}=3\times 10^{13}`$ cm, the Thomson depth of the Comptonizing sphere is $`\tau =1`$, and the radius of the sphere takes values $`r_s=1.5`$$`5\times 10^{13}`$ cm. The smaller values of $`r_s`$ correspond to the harder X-ray spectra. The dashed curve is the SED corresponding to $`r_{\mathrm{min}}=r_s=6\times 10^{13}`$ cm. The dotted lines indicate the 2–10 keV band.

X-ray spectral index versus 2–10 keV flux for families of curves each of which produces a different constant optical continuum level. The rightmost curve represents the model calculations shown in Figure Thermal Comptonization and Disk Thermal Reprocessing in NGC 3516 and the next four curves to its left have been computed for 87.5, 75, 62.5, and 50% of the mean optical flux level for NGC 3516. The X-ray emission can vary significantly and still produce a nearly constant optical/UV continuum if the source follows a path along one of these curves such as the one labeled “NGC 3516”. The thick solid curve has been fit to the full optical flux, but emission due to internal viscous dissipation in the disk has been included to reduce the required X-ray luminosity. The neighboring dashed curves bound the region of $`\pm 1.5`$% optical variations. The data points with error bars are from spectral fits to the archived 1998 RXTE data for NGC 3516.


text:
# Universal Conductance Distributions in the Crossover between Diffusive and Localization Regimes

## Abstract

The full distribution of the conductance $`P(G)`$ in quasi-one-dimensional wires with rough surfaces is analyzed from the diffusive to the localization regime. In the crossover region, where the statistics is dominated by only one or two eigenchannels, the numerically obtained $`P(G)`$ is found to be independent of the details of the system with the average conductance $`G`$ as the only scaling parameter. For $`G<e^2/h`$, $`P(G)`$ is given by an essentially “one-sided” log-normal distribution. In contrast, for $`e^2/h<G2e^2/h`$, the shape of $`P(G)`$ remarkable agrees with those predicted by random matrix theory for two fluctuating transmission eigenchannels.

The large conductance fluctuations occurring in mesoscopic samples have attracted large attention during the last decade. The magnitude of the fluctuations ($`e^2/h`$) is almost independent of the mean value of the conductance $`G`$ and the system size (‘universal conductance fluctuations’). This holds for wires or waveguides in the diffusive regime where the length of the wire $`L`$ is much larger than the mean free path $`\mathrm{}`$ but still much smaller than the localization length $`\xi `$. As the length of wire $`L`$ approaches the localization length $`\xi `$, the conductance fluctuations become of the same order as the averaged conductance. The averaged values are then not enough to describe the transport properties and the knowledge of the conductance probability distribution is of primary interest.

While for a single-mode wire, where $`\xi \mathrm{}`$, the probability distribution $`P(G)`$ can be computed for any system length, it is remarkably difficult to extend this result to the N-mode case . Based on numerical simulations as well as on perturbation theory, $`P(G)`$ is expected to evolve from a Gaussian (deep in the diffusive regime) to a lognormal distribution (deep in the localization regime). However, the behavior of $`P(G)`$ in the crossover between diffusive and localized regimes is not well known. Recently, both experimental and numerical  studies have reported that the conductance distribution at an integer quantum hall transition presents a striking shape in the form of a broad distribution. Related to that, recent numerical studies suggested a universal distribution at the metal-insulator transition. Based on the Dorokhov-Mello-Pereyra-Kumar (DMPK) equation  for quasi-1D wires in absence of time reversal symmetry it has been found that the crossover region between metallic and insulating regimes is highly non-trivial, and shows one-sided log-normal distribution for the conductance at the transition. It was suggested  that these distributions might not be neither a private issue of the quantum Hall transition nor a product of the removal of the time reversal symmetry due to the presence of magnetic fields. Moreover, numerical calculations in 2D and 3D suggest that the form of $`P(G)`$ at the critical point is independent of the dimensionality of the system and of the model.

Here we present the results of extensive numerical calculations of the conductance distribution $`P(G)`$ and transport eigenchannels all the way from the diffusive to the localization regime, paying special attention to the regime crossover. We show that close to the crossover regime, where the averaged conductance $`G`$ is between 1/2 and 1 (in units of $`2e^2/h`$), the general characteristics of $`P(G)`$ are dominated just by two fluctuating eigenchannels. The diffusion-localization transition is fully characterized by a transition from two to only one eigenchannel, occurring at a critical value $`G_c1/2`$. Right at the critical value $`G_c`$, the distribution is almost perfectly flat and almost identical to that obtained for 2D systems in the quantum Hall regime. In the insulating regime ($`G<G_c`$) our results are consistent with the one-side log-normal distributions predicted by the DMPK equation .

Most of the previous work on transport properties of disordered wires has been focused on the influence of multiple scattering with bulk defects. In this work we are going to deal with surface corrugated wires (SCW) as a model system. The interest on these systems is twofold. On one hand, for very thin wires, the main source of multiple scattering comes only from the surface roughness. On the other hand, in absence of bulk defects, the transport through SCW, like in the DMPK approach, does not contain the effects of wave function correlations in the transverse direction. Many statistical properties of transport through SCW  have been found to be in good agreement with those predicted by random matrix theory (RMT) for bulk defects. In particular, the analytical distributions of transmission coefficients  obtained from DMPK equation are nicely reproduced in numerical simulations of SCW . Our model system consists of a surface corrugated two-dimensional (2D) wire (see inset in Fig. 1). The corrugated part of the wire, of total length $`L`$ and perfectly reflecting walls, is composed of slices of length $`l`$. The width of each slice has random values uniformly distributed between $`W_0\delta `$ and $`W_0+\delta `$ about a mean value $`W_0`$. The transmission matrix $`t`$ is exactly calculated for each realization by solving the 2D wave equation by mode matching at each slice, together with a generalized scattering matrix technique. The conductance is given by $`G=\text{trace}\{tt^{}\}`$. To obtain the mean values ($`\mathrm{}`$) and the probability distributions we have used 10000 independent realizations of the disordered guide.

We have performed calculations for three different values of $`W_0/\lambda `$: 4.9, 2.6 and 1.8, allowing $`N=`$ 9, 5 and 3 propagating modes respectively ($`\lambda `$ would be the Fermi wavelength for electron transport or the wavelength of the diffusive incoming waves in the case of transport of classical waves). We have considered $`l/\delta =3/2`$ and different ratios $`W_0/\delta =7`$ (for $`N=5,3`$) and $`W_0/\delta =13.25`$ (for $`N=9`$). The behavior of different mean values of transport coefficients as a function of the length $`L`$ has been discussed previously in detail . In the diffusive regime, the averaged resistance follows a typical ohmic behavior ($`1/G1/N+L/\xi `$) while, in the localization regime, $`\mathrm{ln}GL/\xi `$. In Figure 1 we plot our results both for $`G`$ (Fig. 1a) and for $`\mathrm{ln}G`$ (Fig. 1b) versus $`sNL/\xi `$ for $`N=9`$ (dashed line), $`N=5`$ (solid line) and $`N=3`$ (long-dashed line). The corresponding localization lengths are $`\xi =187W_0`$ ($`N=9`$), $`\xi =34.4W_0`$ ($`N=5`$) and $`\xi =51W_0`$ ($`N=3`$).

The analysis of the conductance distributions $`P(G)`$ along the transition from the diffusive regime $`L<\xi `$ to the localization regime $`L>\xi `$ (see Fig. 2) shows that the distributions are far from evolve smoothly from the Gaussian to the lognormal distribution as one could expect from the behavior of the transmission coefficients . Our results for $`P(G)`$ are shown in Figure 2 for some selected values of $`G`$ around the transition from the diffusive to the localization regime ($`G`$ 1 (a), 4/5 (b), 1/2 (c) and 1/3 (d); those values correspond to the symbols in Fig. 1). The results for $`N=9`$ ($`\mathrm{}`$), $`N=5`$ ($``$) and $`N=3`$ ($`\mathrm{}`$) are almost identical (within our numerical accuracy) suggesting that the only parameter determining the shape of $`P(G)`$ is the average conductance $`G`$ as expected from scaling arguments .

It is remarkable that the almost perfectly flat distribution obtained for $`G=1/2`$ (Fig. 2c) is in full agreement with those obtained for 2D systems in the quantum Hall regime at the critical point : a flat distribution except for a small dip at $`G0`$  and an exponential cutoff for $`G>1`$. As the averaged conductance goes slightly below 1/2, $`P(G)`$ becomes a “one-sided” log-normal distribution with a sharp cut off at $`G>1`$ (Fig. 3) in agreement with both numerical  and analytical  results. In Figure 3 we have plotted our results for $`P(\mathrm{ln}G)`$ for different values of $`G`$ together with the best fit to a one-sided log-normal distribution (thick solid lines). Although the agreement is very good for all values of $`G<1/2`$, the distribution at the critical value of 1/2 is clearly better described by $`P(\mathrm{ln}G,G=1/2)=\mathrm{exp}(\mathrm{ln}G)`$ (thin line in Fig. 3a) which corresponds to a flat distribution ($`P(G,G=1/2)=1`$ , $`0G1`$). Although there is no phase transition in quasi-one-dimensional systems (1D), our results suggest a special generic behavior at a critical $`G_c=1/2`$ which is certainly not consistent with a one-sided log-normal distribution.

In order to get more insight on the exhibited behavior of the conductance distributions we have made an extensive analysis in terms of the transmission eigenchannels $`\{\tau _i\}`$. These eigenchannels are the eigenvectors of the matrix $`tt^{}`$ and they form the natural basis to analyze the properties of the conductance ($`G=\text{trace}(tt^{})=_i\tau _i`$). Following Imry’s work, we arrange the eigenvalues is descending order, and perform the averaging only taking those eigenvalues of the same level into account (e.g. $`\tau _1`$ corresponds to the average of the highest eigenvalue of each realization). In Fig. 4 we show the evolution of the average value of the eigenchannels as a function of $`G`$. In the diffusive regime ($`G>1`$), in agreement with the behavior discussed by Imry, almost all the eigenchannels are either fully open or closed and only a few of them fluctuate giving rise to UCF. However, their exact behavior depends on the particular parameters of the wire. Close to the crossover to localization, when all the eigenchannels are closed and only two (or one in the localization threshold) of them are partially open, the average eigenchannel transmissions (Fig. 4) as well as the conductance distributions (Fig. 2) do not depend on the size of the wire or the defect details. In this sense, our results suggest an universal behavior of the conductance with $`G`$ as the only scaling parameter.

In contrast with the insulating regime $`G<1/2`$, where the conductance is known to be well described by one-sided log-normal distributions , in the crossover regime $`G>1/2`$, where the statistics is dominated by one or two fluctuating channels, there is no analytical result available for the conductance distribution. As a first analytical approach to this problem, given our limited knowledge of the statistical correlations between different eigenchannels, a possible choice of the statistical ensemble is that which maximizes the information entropy subject to the known constrains of flux conservation and time-reversal invariance. In the RMT context, this leads to the circular orthogonal ensemble. In our case, we also know that, close to the onset of localization (for $`1/2G1`$), only two eigenchannels actually contribute to the conductance, independently of the initial number of channels $`N`$. From the polar decomposition of the scattering matrix  it is possible to write the joint probability distribution of $`n`$ transmission eigenchannels (with $`\tau _i=0`$, $`i=n+1,\mathrm{},N`$) as

$$P(\{\tau _i\})\underset{i<j}{\overset{n}{}}|\tau _i\tau _j|\underset{k}{\overset{n}{}}\tau _k^\alpha $$
(1)
where $`\alpha =(1+n)(Gn/2)/(nG)`$ . Within this approach, the conductance distribution is fully described by the number of fluctuating channels $`n`$ and the mean value of the conductance $`G`$ . Simple closed expressions for $`P(G;G,n)`$ can be obtained for $`n=1,2`$,

$$P(G;G,1)G^{\frac{12G}{1G}},\text{ for }\text{0 ¡ G ¡ 1},$$
(2)
and

$$P(G,G,2)\{\begin{array}{cc}\left(\frac{G}{2}\right)^{\frac{2(12G)}{G2}},\hfill & \text{if }0<G1\hfill \\ \left(\frac{G}{2}\right)^{\frac{2(12G)}{G2}}\left(G1\right)^{\frac{12G}{G2}},\hfill & \text{if }1G<2\hfill \end{array}$$
(3)
In Fig. 2 we have plotted the RMT distributions (continuous thick line) for different $`G`$ ($`n=1`$ for $`G1/2`$ and $`n=2`$ for $`1/2G`$), together with our numerical results. As it can be seen, there is a very good agreement between the analytical results and the numerical calculations. In particular, our RMT approach captures some important features of the distributions such as the almost perfectly flat distribution at $`G=1/2`$ and the cusp point in $`P(G,G)`$ at $`G=1`$. Since the RMT results are quite general, we expect qualitative similar conductance distributions in higher dimensions close to the critical regime.

In summary, we have analyzed the evolution of the conductance distribution all the way from the diffusive to the localization regime. Close to the crossover regime we have shown that the distributions are independent of the details of the system with the averaged conductance $`G`$ as the only scaling parameter. For $`G`$ between 1/2 and 1, the statistics is dominated by one or two fluctuating eigenchannels and the numerical conductance distributions are surprisingly well described by RMT results At a critical value $`G_c=1/2`$, the distribution is almost perfectly flat. In the insulating regime ($`G<1/2`$) our results are consistent with one-sided log-normal distributions. The similarities between our results and those obtained in very different situations suggest that the conductance distributions exhibit a universal behavior at the crossover regime.

We specially acknowledge P.A. Mello, M. Nieto-Vesperinas, J.-L. Pichard and P. Wölfle for their fruitful and stimulating suggestions. We are grateful to C. Soukoulis for discussions and bringing the results of the quantum Hall transitions to our attention. We also like to thank A. Cano and L.S. Froufe for interesting discussions. This work has been supported by the Comunidad Autónoma de Madrid and the DGICyT through Grants 07T/0024/1998 and No. PB98-0464.


text:
# Field Induced Partial Order in the Spin-Ice Dysprosium Titanate.

## 1 Introduction

Ho2Ti2O7 and Dy2Ti2O7 are regarded as good examples of spin-ice materials prl1 ; newprl ; art . The magnetic rare earth ions occupy a pyrochlore lattice, a cubic array of corner-linked tetrahedra. In the near neighbour spin ice model the lattice is occupied by $`111`$ Ising spins that are coupled ferromagnetically to nearest neighbours jpcm . The result is a disordered, macroscopically degenerate, ground state that is determined only by the rule that two spins must point into and two out of each tetrahedron. This “two in, two out” rule is analogous to the “ice rules” that control the proton arrangement in water ice, and so the near neighbour spin ice model maps exactly onto Pauling’s model of the proton disorder in water ice jpcm ; pauling . The near neighbour spin ice model qualitatively describes much of the behaviour of the real materials, but for an accurate description the dipolar spin ice model has been developed byron . In this model, the near neighbour ferromagnetic coupling between $`111`$ spins is dipolar in origin and it has been shown that the long range part of the dipolar coupling maintains a spin-ice like disordered ground state even in the presence of an antiferromagnetic near neighbour super-exchange byron . The dipolar spin ice model has been shown to give an accurate description of the zero field magnetic neutron scattering of Ho2Ti2O7 newprl  in the spin ice regime ($`T<2`$ K) as well as the specific heat of Ho2Ti2O7 and Dy2Ti2O7 throughout the low temperature range (50 mK - 20 K) newprl ; byron . The dipolar coupling in Ho2Ti2O7 and Dy2Ti2O7 is of similar magnitude ($`D_{nn}2.35`$K for both Dy and Ho) while the antiferromagnetic exchange is weak in both cases (-1.2K for Dy and -0.52K for Ho) newprl ; byron . The main difference between these two materials appears to arise from the nuclear spins. Ho2Ti2O7 has a single isotope (165Ho) and strong hyperfine coupling that complicates the spin ice freezing process below $`0.8`$ K newprl . Dy2Ti2O7, on the other hand, has much weaker hyperfine coupling that does not affect the dominant spin ice behaviour in the temperature range of interest. It might therefore be considered the simpler spin ice material. However, the natural mixture of Dy isotopes leads to significant neutron absorption (994 barn for 2200ms-1 neutrons). Thus, although neutron scattering should be the ideal microscopic probe of the spin correlations and dynamics of this system, the most significant experimental observations so far are bulk measurements art . To minimise the absorption problem we have prepared a single crystal of 162Dy2Ti2O7. In this paper we present the initial results of our neutron scattering study of this sample, in both zero and applied magnetic field.

## 2 Experimental

The absorption cross section of natural Dysprosium is $`\sigma _a=994`$ barn, the large value being due mainly to the presence of 164Dy. Isotopically enriched Dy2O3 with composition 156Dy$`<0.01\%`$, 158Dy$`<0.01\%`$,160Dy$`=0.02\%`$,
161Dy$`=0.47\%`$, 162Dy$`=96.8\%`$,163Dy$`=2.21\%`$ and 164Dy$`=0.5\%`$ was prepared by Goss Scientific Instruments Ltd. Enrichment in 162Dy was considered optimum as it has a high abundance and the dominant impurity is 163Dy which has a slightly lower absorption cross section. Enrichment in 163Dy would have lead to 164Dy as the dominant impurity. Through enrichment, the absorption cross section was reduced by a factor of five to $`208\pm 13`$ barn.

A single crystal of Dy2Ti2O7 was prepared from the isotopically enriched Dy2O3 and TiO2 by the floating zone technique don . Neutron scattering was carried out at the ISIS facility on the indirect geometry spectrometer PRISMA. It was configured in the diffraction mode so that sixteen 3He tube detectors are used simultaneously. Rotation of the crystal allows a rapid mapping of a large section of reciprocal space, making PRISMA ideal for observing magnetic diffuse scattering.

In the first experiment the crystal was cooled by a 3He sorption refrigerator for measurements in zero applied magnetic field. It was aligned with $`[1\overline{1}0]`$ vertical, such that the scattering plane contained $`(hhl)`$ wavevectors. In the second experiment an Oxford instruments 7T vertical field cryomagnet, with dilution refrigerator insert, was used. For the applied field measurements two field orientations were studied: $`[1\overline{1}0]`$ as above and $`[100]`$ corresponding to a $`(0kl)`$ scattering plane. The data were normalised to monitor counts and vanadium to remove the characteristic flux profile of the spallation source. Since the sample is still significantly absorbing a correction for this was also applied.

## 3 Results

With $`[100]`$ vertical, a map of reciprocal space was made at the base temperature ($`70`$ mK) in zero field. The diffuse scattering maxima observed agree with those predicted in the $`(hhl)`$ scattering plane for the dipolar spin ice model. For example there is a diffuse feature at $`3,0,0`$ in $`(hk0)`$ and $`0,0,3`$ in $`(hhl)`$ newprl . On application of a field the diffuse scattering disappeared and was replaced by magnetic Bragg peaks at the $`Q=0`$ positions. A field of $`0.7`$T was sufficient to saturate these peaks.

This behaviour is readily explained as the applied field breaks the degeneracy of the six “two in, two out” spin configurations of the elementary tetrahedron. The pyrochlore lattice can be described as a face centred cubic lattice with a tetrahedral basis, and the degeneracy breaking means that every tetrahedron adopts the same “two in, two out” state with a net moment in the direction of the applied field, $`[100]`$. This non-collinear ferromagnetic structure allows the observed magnetic Bragg peaks at the $`Q=0`$ positions such as $`2,0,0`$. Interestingly the experimental magnetization did not develop smoothly, but in a series of steps. Hysteresis was observed on cycling the field.

The $`[1\overline{1}0]`$ direction is a hard direction of magnetization Cornelius . In zero field, the magnetic scattering was observed to be diffuse, and characteristic of the disordered low temperature state of dipolar spin ice newprl , see Figure 1. Application of a small field ($`0.1`$ T) in this direction again caused the appearance of the $`Q=0`$ Bragg peaks; however, unlike the $`[100]`$ direction, the diffuse scattering features did not disappear. As the field was raised up to $`1.5`$T the diffuse scattering sharpened around the $`Q=X`$ positions such as $`0,0,1`$ without becoming resolution limited. This is seen clearly in figure 2. Similar features have been observed in the neutron scattering of Ho2Ti2O7 prl1 ; dten , adding strength to the idea that the true ground state of these systems is a $`Q=X`$ structure that is dynamically inhibited from being accessed on experimental timescales Michel . Again the magnetization versus field curve was observed to have several sharp steps and plateaus (Figure 3).

The formation of the $`Q=X`$ structure is consistent with the spin ice rules prl1 . In this field orientation, assuming perfect $`111`$ spins, only two of the spins of the tetrahedral basis have a component along the field direction. These form “in-out” $`[1\overline{1}0]`$ chains parallel to the field forcing the remaining two spins per tetrahedra into “in-out” $`[110]`$ chains perpendicular to the field, as illustrated in Figure 4. The perpendicular chains are not coupled by the spin ice rules which, in the absence of any further neighbour coupling, would lead to two dimensional $`[110]`$ Bragg sheets of scattering extended along $`[00l]`$ in the scattering plane. The diffuse features in the experimental pattern are indeed extended along this direction, but the sharp build up of intensity around the $`Q=X`$ points indicates a strong tendency to prefer $`Q=X`$ short range ordering of the perpendicular rods (see Figure 2). Neglecting interference between the scattering from the two spin sets (perpendicular and parallel to the field) one arrives at the conclusion that the $`Q=0`$ Bragg scattering arises from the parallel rods and the $`Q=X`$ diffuse scattering arises from the perpendicular rods. However, it may be a crude approximation to separate the scattering in this way.

In conclusion, we have obtained accurate neutron scattering data for Dy2Ti2O7 that is in qualitative agreement with theoretical expectations for a spin ice material prl1 ; newprl ; jpcm ; byron . It is of interest that, even in a relatively strong field along $`[1\overline{1}0]`$, the system remains only partially ordered. It is also noteworthy that the magnetic hysteresis loop shows several steps and plateaus. An understanding of these effects awaits a detailed study of the static and dynamic properties of the dipolar spin ice model byron  in an applied magnetic field.

## 4 Acknowledgements

We acknowledge the EPSRC for funding of beam time and studentships (TF and JDMC).


text:
# Enumerating Foldings and Unfoldings between Polygons and Polytopes

## 1 Introduction

We explore the process of folding a simple polygon by gluing its perimeter shut to form a convex polyhedron, and its reverse, cutting a convex polyhedron open and flattening its surface to a simple polygon. We restrict attention to convex polyhedra (henceforth, polytopes), and to simple (i.e., nonself-intersecting, nonoverlapping) polygons (henceforth, polygons). The restriction to nonoverlapping polygons is natural, as this is important to manufacturing applications \[O’R00\]. The restriction to convex polyhedra is made primarily to reduce the scope of the problem. See \[BDD+98\] and \[BDEK99\] for a start on unfolding nonconvex polyhedra.

We enumerate foldings and unfoldings based on two criteria of indistinguishability: geometric congruence, and combinatoric equivalence. The latter especially will need further specification to become precise, but to presage our results crudely, we show that both the number of foldings and the number of unfoldings can be exponential in the number of vertices $`n`$ of the polygon/polytope. Similarly, we show that polygons may fold and polytopes unfold to an infinite number of incongruent polytopes/polygons. We obtain sharper results when attention is restricted to convex polygons. Proofs and details not provided here may be found in \[DDLO00\].

We will use $`P`$ throughout the paper for a polygon, and $`Q`$ for a polytope, $`P`$ and $`Q`$ respectively for their boundaries, and $`n`$ for the number of their vertices.

## 2 Aleksandrov’s Theorem

A key tool in our work is a far-reaching generalization of Cauchy’s rigidity theorem proved by Aleksandrov \[Ale58\] that gives simple conditions for any folding to a polytope. A gluing maps $`P`$ to $`P`$ in a length-preserving manner, as follows. $`P`$ is partitioned by a finite number of distinct points into a collection of open intervals whose closure covers $`P`$. Each interval is mapped one-to-one (i.e., glued) to another interval of equal length. Corresponding endpoints of glued intervals are glued together (i.e., identified). Finally, gluing is considered transitive: if points $`a`$ and $`b`$ glue to point $`c`$, then $`a`$ glues to $`b`$. Aleksandrov proved that any gluing that satisfies these two conditions corresponds to a unique polytope:

1. No more than $`2\pi `$ total face angle is glued together at any point; and
2. The complex resulting from the gluing is homeomorphic to a sphere.

We call a gluing that satisfies these conditions an Aleksandrov gluing. Although an Aleksandrov gluing of a polygon forms a unique polytope, it is an open problem to compute the three-dimensional structure of the polytope \[O’R00\]. Henceforth we will say a polygon folds to a polytope whenever it has an Aleksandrov gluing.

We should mention two features of Aleksandrov’s theorem. First, the polytope whose existence is guaranteed may be flat, that is, a doubly-covered convex polygon. We use the term “polytope” to include flat polyhedra. Second, condition (2) specifies a face angle $`2\pi `$. The case of equality with $`2\pi `$ leads to a point on the polytope at which there is no curvature, i.e., a nonvertex. We make explicit what counts as a vertex below.

## 3 Geometrical Congruence

In this section we address these two natural questions:

1. How many geometrically different polytopes may be folded from one polygon?
2. How many geometrically different polygons may be unfolded from one polytope?

Here “geometrically different” means incongruent. Although we mentioned the rough answer to both question is ‘infinite,’ there are several nuances in the details. For example, the answer to the first question is: ‘sometimes infinite,’ whereas the answer to the second is: ‘always infinite.’

### 3.1 Congruence: Folding

We start with a natural and easily proved claim:

###### Lemma 1

Some polygons cannot be folded to any polytope.

An example is shown in Fig. 2.

It is natural to wonder what the chances are that a random polygon could fold to a polytope. This is difficult to answer without a precise definition of “random,” but we feel any reasonable definition would lead to the same answer: zero. We provide support for this conjecture in \[DDLO00\]. Despite this evidence for the rare ability to fold, convex polygons are fertile in their folding, as we now demonstrate.

Let $`xP`$ be an arbitrary point on the boundary of $`P`$, and let $`yP`$ be the midpoint of perimeter $`L`$ around $`P`$ measured from $`x`$. Let $`(x,y)`$ be the open interval of $`P`$ counterclockwise from $`x`$ to $`y`$. Define a perimeter-halving gluing as one which glues $`(x,y)`$ to $`(y,x)`$. A consequence of Aleksandrov’s theorem is:

See Fig. 2 for an example.

###### Lemma 2

Every convex polygon folds to a polytope via perimeter halving for every $`xP`$.

This result can be strengthened:

###### Theorem 3.1

Any convex polygon $`P`$ folds, via perimeter halving, to a nondenumerably infinite number of noncongruent polytopes.

Using a different type of folding, we can show that every rectangle folds to a continuum of tetrahedra. See Fig. 3.

### 3.2 Congruence: Unfolding

Although it is a long-standing open problem to determine whether every polytope may be cut along polytope edges and unfolded to a polygon, without the “along edges” restriction it is easy to see that every polytope may be cut open to a continuum of noncongruent polygons. To avoid trivial zigzaging of the cuts, it makes sense to restrict the cuts to be geodesics, which unfold (or “develop”) to straight lines, or restrict even further to shortest paths, geodesics which are in addition shortest paths between their endpoints. Still this holds:

###### Lemma 3

Every polytope $`Q`$ may be cut via shortest paths to unfold to a nondenumerably infinite number of noncongruent polygons.

Proof: (Sketch) This may be accomplished via the star-unfolding \[AAOS97\], which cuts along the shortest paths from a source point $`s`$ to every vertex of $`Q`$. For any point $`p`$ in the interior of a “ridge-free region” of $`Q`$, every $`s`$ in a neighborhood of $`p`$ yields a distinct star-unfolding. $`\mathrm{}`$

## 4 Combinatorial Equivalence

Although a natural counterpart to our geometric enumerations would count combinatorially distinct polygons and polytopes, the former class is uninteresting and the latter class seems difficult to capture.111 Some results for convex unfoldings were obtained by Shephard \[She75\].  Instead we focus on the process of folding and unfolding, and ask:

1. How many combinatorially different foldings of a polygon lead to a polytope?
2. How many combinatorially different cuttings of a polytope lead to polygon unfoldings?

It requires some care to define an appropriate notion of “combinatorially different” for both questions.

### 4.1 Combinatorics: Folding

We capture the combinatorics of a polygon folding via its “gluing tree.” Let a polygon $`P`$ have vertices $`v_1,\mathrm{},v_n`$, labeled counterclockwise, and edge $`e_i`$, $`i=1,\mathrm{},n`$ the open segment of $`P`$ after $`v_i`$. The combinatorial gluing tree $`T_G`$ is a labeled tree representing the identification of $`P`$ with itself.

Any point of $`P`$ that is identified with more or less than one other distinct point of $`P`$ becomes a node of $`T_G`$, as well as any point to which a vertex is glued. (Note that this means there may be nodes of degree $`2`$.) So every vertex of $`P`$ maps to a node of $`T_G`$; each node is labeled with the set of all the elements (vertices or edges) that are glued together there. A point of $`P`$ in the interior of a polygon edge that glues only to itself, i.e., where a crease folds the edge in two, we call a fold point. Points $`x`$ and $`y`$ in Fig. 2 are fold points. A fold-point correspond to a leaf of $`T_G`$, and is labeled by the edge label only. Every nonleaf node has at least one vertex label, and at most one edge label. An example is shown in Fig. 4. The polygon shown folds to a tetrahedron by creasing as illustrated in (a). All four tetrahedron vertices are fold points. The corresponding gluing tree is shown in (b) of the figure. The two interior nodes of $`T_G`$ have labels $`\{v_1,v_6,e_1\}`$ and $`\{v_2,v_5,e_5\}`$.

We start with a characterization of the structure of gluing trees, which will form the basis of our enumeration results. Several combinatorial tree structures play a special role, and to which we assign symbols:

1. |: a path.
2. Y: a tree with a single degree-$`3`$ node.
3. I: a tree with two degree-$`3`$ nodes connected by an edge (e.g., Fig. 4b).
4. +: a tree with single degree-$`4`$ node.

Next, define a belt in a gluing tree to be a path between two leaf fold points; for example, between the $`e_1`$ and $`e_3`$ fold points in Fig. 4b. A belt is a rolling belt if there is a nonzero-length interval $`Ie`$ such that for every $`xI`$, the belt folded at $`x`$ is an Aleksandrov gluing. (A belt could instead have a finite number of distinct gluings, perhaps just one.)

Our characterization shows that gluing trees are fundamentally discrete structures, with one or two rolling belts, and two such belts only in very special circumstances:

###### Theorem 4.1

Gluing trees satisfy these properties:

1. At any gluing tree point of degree $`d2`$, at most one point of $`P`$ in the interior of an edge may be glued, i.e., at most one nonvertex may be glued there.
2. At most four leaves of the gluing tree can be fold points, i.e., points in the interior of an edge of $`P`$. The case of four fold-point leaves is only possible when the tree has exactly four leaves, with the combinatorial structure ‘+’ or ‘I’.
3. A gluing tree can have at most two rolling belts.
4. A gluing tree with two rolling belts must have the structure ‘I’, and result from folding a polygon that can be viewed as a quadrilateral with two of its opposite edges replaced by complimentary polygonal paths.

Thus a generic gluing tree has one rolling belt, with trees hanging off it, and one of those trees having a fold-point leaf, as depicted in Fig. 5.

We use this characterization to prove bounds on the number of gluings. First, a lower bound:

###### Theorem 4.2

For any even $`n`$, there is a polygon $`P`$ of $`n`$ vertices that has $`2^{\mathrm{\Omega }(n)}`$ combinatorially distinct Aleksandrov gluings.

Proof: (Sketch) The polygon $`P`$ is illustrated in Fig. 6(a). It is a centrally symmetric star, with $`m`$ vertices, $`m`$ even, with a small convex angle $`\alpha 0`$, alternating with $`m`$ vertices with large reflex angle $`\beta <2\pi `$. All edges have the same (say, unit) length. We call this an $`m`$-star. We choose $`\alpha `$ small enough so that $`m`$ copies of $`\alpha `$ can join with one of $`\beta `$ and still be less than $`2\pi `$. Now we add two vertices $`x`$ and $`y`$ at the midpoints of edges, symmetrically placed so that $`y`$ is half the perimeter around $`P`$ from $`x`$. Let $`n=n^{}+2`$ be the total number of vertices of $`P`$.

The “base” gluing tree is illustrated in Fig. 6(b). $`x`$ and $`y`$ are fold vertices of the gluing. Otherwise, each $`\alpha `$ is matched with a $`\beta `$. Because all edge lengths are the same, and because $`\alpha +\beta <2\pi `$, this path is an Aleksandrov gluing. We label it $`T_{00\mathrm{}0,00\mathrm{}0}`$, where $`m/2`$ zeros $`00\mathrm{}0`$ represent the top chain, and another $`m/2`$ zeros represent the bottom chain.

The other gluing trees are obtained via “contractions” of the base tree. A contraction makes any particular $`\beta `$-vertex not adjacent to $`x`$ or $`y`$ a leaf of the tree by gluing its two adjacent $`\alpha `$-vertices together. Label a $`\beta `$-vertex $`0`$ or $`1`$ depending on whether it is uncontracted or contracted respectively. Then a series of contractions can be identified with a binary string. For example, Fig. 6(c) displays the tree $`T_{010100\mathrm{},00110\mathrm{}0}`$.

We can bound the number of Aleksandrov gluings resulting from these contractions by $`\mathrm{\Omega }(2^{m/21})=2^{\mathrm{\Omega }(n)}`$. $`\mathrm{}`$

It may not be surprising that some polygons have many foldings, but it is perhaps less intuitive that even simple polygonal shapes have many gluings. For example, our enumeration program finds that an equilateral triangle has 19 gluings and a square 43 gluings. Of the latter, 10 foldings are distinct when symmetries are removed: several flat shapes, four tetrahedra, a hexahedron, and a continuum of octahedra. Hirata \[Hir00\] has shown that the Latin square, whose study we initiated in \[LO96\], has 85 distinct gluings. These lead to 21 distinct shapes, including several flat quadrangles, tetrahedra, hexahedra (including a cube), octahedra, and a pentahedron.222 http://daisy.uwaterloo.ca/~eddemain/aleksandrov/cross/ .  Fig. 7 shows crease patterns for two gluings.

We may obtain an upper bound in terms of the number of leaves $`\lambda `$ of the gluing tree:

###### Theorem 4.3

The number of gluing trees with $`\lambda `$ leaves for a polygon $`P`$ with $`n`$ vertices is $`O(n^{2\lambda 2})`$.

This bound is useful when the number of leaves is bounded, which is the case, for example, with convex polygons. The characterization of Theorem 4.1 can be tightened in this case:

###### Lemma 4

For convex polygons, the gluing tree $`T_G`$ has one of these combinatorial structures: when $`n4`$, either ‘|’ or ‘Y’; when $`n=4`$, in addition ‘I’ and ‘+’ are possible.333 This lemma is largely due to Shephard \[She75\].

This leads to a tighter bound for convex polygons:

###### Theorem 4.4

A convex polygon $`P`$ of $`n`$ vertices folds to at most $`O(n^3)`$ different gluing trees. Some convex polygons have $`\mathrm{\Omega }(n^2)`$ gluings.

We leave open the task of closing the gap between quadratic and cubic.

### 4.2 Combinatorics: Unfolding

Finding the “right” way to count unfoldings is more delicate. We start by defining cut trees, which then form the basis of our enumerations. It will be useful to distinguish between a geometric tree $`𝒯`$ composed of a union of line segments, and the more familiar combinatorial tree $`T`$ of nodes and arcs. A geometric cut tree $`𝒯_C`$ for a polytope $`Q`$ is a tree drawn on $`Q`$, with each arc a polygonal path, which leads to a polygon unfolding when the surface is cut along $`𝒯`$, i.e., flattening $`Q𝒯`$ to a plane.

###### Lemma 5

If a polygon $`P`$ folds to a polytope $`Q`$, $`P`$ maps to a tree $`𝒯_CQ`$, the geometric cut tree, with the following properties:

1. $`𝒯_C`$ is a tree.
2. $`𝒯_C`$ spans the vertices of $`Q`$.
3. Every leaf of $`𝒯_C`$ is at a vertex of $`Q`$.
4. A point of $`𝒯_C`$ of degree $`d`$ (i.e., one with $`d`$ incident segments) corresponds to exactly $`d`$ points of $`P`$. Thus a leaf corresponds to a unique point of $`P`$.
5. Each arc of $`T_C`$ is a polygonal path on $`Q`$.

There are several options in defining the combinatorial tree $`T_C`$ for a geometric $`𝒯_C`$:

1. Make every segment of $`𝒯_C`$ an arc of $`T_C`$. Although this is natural, allowing an arbitrarily complicated polygonal path between any two polytope vertices leads to an infinite number of different cut trees for any polytope.
2. Make every point where a path of $`𝒯_C`$ crosses an edge of the polytope a node of $`T_C`$. This again leads to trivially infinite numbers of cut trees when a path of $`𝒯_C`$ zigzags back and forth over an edge of $`Q`$.
3. Exclude this possibility by forcing the paths between polytope vertices to be geodesics, and again make polytope edge crossings nodes of $`T_C`$. This excludes many interesting cut trees and destroys symmetry between $`T_C`$ and $`T_G`$.
4. Make every maximal path of $`𝒯_C`$ consisting only of degree-$`2`$ points a single arc of $`T_C`$. A consequence is that polytope vertices in the interior of such a path disappear from $`T_C`$.

Threading between these possibilities, we define the combinatorial cut tree $`T_C`$ corresponding to a geometric cut tree $`𝒯_C`$ as the labeled graph with a node (not necessarily labeled) for each point of $`𝒯_C`$ with degree not equal to $`2`$, and a labeled node for each point of $`𝒯_C`$ that corresponds to a vertex of $`Q`$ (labeled by the vertex label); arcs are determined by the polygonal paths of $`𝒯_C`$ connecting these nodes. An example is shown in Fig. 8. Note that not every node of the tree is labeled, but every polytope vertex label is used at some node. All degree-$`2`$ nodes are labeled.

This definition has the consequence that, if all degree-$`2`$ nodes are removed by contraction, $`T_C`$ is isomorphic to the corresponding gluing tree $`T_G`$. Although the definition avoids some of the listed pitfalls, it does have the undesirable consequence of counting different geodesics on $`Q`$ between two polytope vertices as the same arc of $`T_C`$ even if one spirals around the polytope twice and the other once (or not at all).

Before turning to enumeration bounds, we make this straightforward observation: Every polytope admits at least the $`n`$ cut trees provided by the star-unfolding \[AO92\], one with each vertex as source. So in particular, every polytope has at least one unfolding to a simple polygon, in contrast to the corresponding open question for edge-unfoldings (p.3.2).

Our main result here is that some polytopes have an exponential number of unfoldings:

###### Theorem 4.5

There is a polytope $`Q`$ of $`n`$ vertices that may be cut open with exponentially many ($`2^{\mathrm{\Omega }(n)}`$) combinatorially distinct cut trees, which unfold to exponentially many geometrically distinct simple polygons.

Proof: (Sketch) $`Q`$ is a truncated cone: the hull of two regular $`n`$-gons of different radii, lying in parallel planes and similarly oriented. Two different cuttings are illustrated in Figs. 10 and 10.

The “base” cut tree, which we notate as $`T_{0000000}`$, unfolds $`Q`$ as shown in Fig. 10. Define a cut tree $`T_{m_{(n1)/2}\mathrm{}m_2m_1m_0}`$, where $`m_i`$ are the digits of a binary number of $`n/21`$ bits, as an alteration of the base tree $`T_{0\mathrm{}0}`$ illustrated by $`T_{1001101}`$ shown in Fig. 10. There are $`2^{n/21}=2^{\mathrm{\Omega }(n)}`$ cut trees, and it is not difficult to show that each leads to a distinct simple polygon unfolding. $`\mathrm{}`$

We conjecture that there is a polytope with an exponential number of convex unfoldings, i.e., those that result in convex polygons.

Our upper bound relies on bounds on the number of spanning trees of triangulated planar graphs:

###### Theorem 4.6

The maximum number of combinatorially distinct edge-unfolding cut trees of a polytope of $`n`$ vertices is $`2^{O(n)}`$, and the maximum number of (arbitrary) combinatorially distinct cut trees is $`2^{O(n^2)}`$.

## 5 Open Problems

Some of the most interesting open questions in this area are algorithmic:

1. Given an Aleksandrov gluing, compute the 3D structure of the polytope. This is an algorithmic version of Aleksandrov’s theorem, for which only a “finite” algorithm is known.444 I. K. Sabitov, Oberwolfach presentation, May 2000. See also http://www.cms.math.ca/CMS/Events/winter98/w98-abs/node46.html.  This problem is closely related to the following problem, which may be easier because of the additional information:
2. Given an Aleksandrov gluing and the unique crease pattern for the folding (Cf. Fig. 7), compute the 3D structure of the polytope. This can be viewed as an algorithmic version of Cauchy’s rigidity theorem.
3. How difficult is it to determine whether a given polytope has a convex unfolding? (Cf. Lemma 3.)
4. How difficult is it to determine whether a given polygon may be folded to a polytope? (Cf. Fig. 2.) We have an exponential-time algorithm, but a polynomial time algorithm is known only for edge-to-edge gluings \[LO96\].