text:
# A NEW WAY TO DETECT MASSIVE BLACK HOLES IN GALAXIES: The Stellar Remnants of Tidal Disruption

## 1 Tidal Disruption of Giants: The Core Remains

A star of mass $`M_{}`$ and radius $`R_{}`$ can be tidally disrupted by approaching within $`R_t`$ of a black hole (BH).

$$R_t=\left(\eta ^2\frac{M_{bh}}{M_{}}\right)^{\frac{1}{3}}R_{},$$
(1)
$`M_{bh}`$ is the BH’s mass; $`\eta `$ is a parameter of order of unity. Let $`R_S`$ be the BH’s Schwarzschild radius. For $`M_{bh}`$ $`>`$ $`10^8M_{},`$ $`R_t<R_S`$ for most main-sequence stars, and only tidal disruptions (TDs) of giants lead to observable effects. The rate of such TDs, $`\dot{N}_G,`$ depends on $`M_{bh}`$, the kinematics of stars near the galactic nucleus, and the mass and age distribution of these stars (see, e.g., Magorrian & Tremaine 1999). Computed rates are typically in the range $`10^610^4`$ yr-1. A convenient expression is

$$\dot{N}_G10^5\left(\frac{L}{L_{}}\right)^{1.2}yr^1,$$
(2)
with $`L_{}=1.8\times 10^{10}L_{};`$ this was derived by Syer & Ulmer (1999), who also provide cautions about the use of a simple scaling law. Studies of TD have focused on possible associated accretion events, lasting for months or decades, with luminosities as high as $`10^{44}10^{46}`$ erg/s (see, e.g., Hills 1975, Lidskii & Ozernoi 1979, Gurzadyan & Ozernoi 1980, Rees 1988, Loeb & Ulmer 1997, Ulmer, Paczyński, & Goodman 1998, Ulmer 1999).

The Hot Core of a Disrupted Giant:  The disruption leaves an end-product, the giant’s hot dense core, whose presence, influence, and observability are the subjects of this paper. The core remains hot ($`T>10^5`$ K) and bright ($`L>10^2L_{}`$) for $`10^310^6`$ years, thereby providing the longest-lasting signal of a TD.

Detectability:  The soft X-ray sensitivities of Einstein, ROSAT, Chandra, and XMM-Newton have allowed them to detect supersoft X-ray sources (SSSs) in other galaxies. SSSs (Greiner et al. 1991, Rappaport, Di Stefano & Smith 1994, Kahabka & van den Heuvel 1997) have luminosities ($`10^{35}10^{38}`$ erg/s) and temperatures ($`kT10100`$ eV), comparable to those of hot stellar cores. Chandra’s 0.5–1′′ angular resolution allows individual sources to be resolved, even in the dense central regions of nearby galaxies.

The Example of M31:  We have studied 8.8 ksec of data collected by Chandra’s ACIS-I detector. To identify SSSs we applied a hardness ratio test to all point sources in a $`16^{}\times 16^{}`$ region centered on the center of M31, selecting all sources with more than 50% of their photons below 0.7 keV. Three sources satisfied this condition. Remarkably, one of these is coincident with the center of M31. The luminosity is $`10^{37}`$ erg/s; the data are consistent with little or no emission above $`1.5`$ keV. The source appears to be variable; data from earlier missions (Einstein, ROSAT), e.g., are consistent with the source providing zero flux on one occasion. Garcia et al. (2000) have argued that the density of point X-ray sources in this field is so low that a chance coincidence with the center is highly unlikely; this argument is even stronger when applied to SSSs. The observed SSS-like behavior is therefore most likely related to the environment of the center, presumably to the presence of a MBH. The data do not seem consistent with ADAF disk models (see, e.g., Garcia et al. 2000), but are compatible with the signature of a hot stellar core. Whether or not the source is a hot core, the M31 observations establish that it is possible to detect such an object in the centers of nearby galaxies.

## 2 The Nature of the Cores

The nature and characteristics of the core depend on the initial mass, $`M_{}`$ of the stripped star, and on its state of evolution at the time of disruption.

### 2.1 $`M_{}<8M_{}`$

Tidal stripping brings a premature halt to nuclear burning. While any residual nuclear burning continues, the luminosity remains at the value it had when the star was a giant. The radius shrinks, and the effective temperature increases. Depending on the core mass, these systems will stay at $`L>10^2L_{}`$, and $`T>10^5`$ K for $`10^310^5`$ yrs (see, e.g., Blöcker 1995; Gorny, Tylenda & Szczerba 1994).

Nature has provided direct analogs, in the form of the hot cores at the centers of planetary nebulae. A notable example is 1E 0056.8-7154, a planetary nebula system in the Small Magellanic Cloud, with $`L=2\times \mathrm{\hspace{0.17em}10}^{37}`$ erg/s and $`kT=3040`$ eV (Wang & Wu 1992).

The evolution of the hot core is influenced by its mass, its state at the time of disruption (e.g., its thermal pulse-cycle phase), and the mass of the original star (see, e.g., Iben 1995 and references therein). The hot core of a star that was disrupted at some arbitrary time in its evolution, is descended from a more massive star than a core of the same mass emerging during the final stages of an uninterrupted evolution. The core of the disrupted star is therefore younger, and is generally less compact, with a significantly higher core temperature; it cools more slowly than typical post-AGB stars. (See, e.g., Blöcker 1995; Gorny, Tylenda, & Szczerba 1994.) Although these results were arrived at by studying winds, they are expected to be valid for the extreme mass-loss scenario associated with disruption.

The stellar mass function favors stars in this mass range, and they are therefore expected to be disrupted more frequently than stars of higher mass. Nevertheless, we go on to consider high-mass stars, because their disruption may also produce remnants with observable signatures.

### 2.2 $`M_{}>8M_{}`$

The stripped core is a helium star (He star). Two processes (in addition to TD) may produce He stars. First, a very massive star may eject such copious winds that it effectively ejects its own hydrogen envelope. Second, a massive star in a binary may overfill its Roche lobe as it evolves and lose its hydrogen envelope by either transferring it to its companion (if the mass ratio is not too extreme), or losing it during a common envelope episode. The evolution of isolated He stars and He stars in binaries has been well-studied; see, e.g., Maeder (1981), Delgado & Thomas (1981), Habets (1986, 1987), de Loore & De Greve (1992), Reese (1993), Woosley, Langer, & Weaver (1995). The characteristics of each stage of the He star’s evolution depend on many factors, including the mass of the star when its envelope was lost, and the continuing rate of mass loss to winds. A common feature of these evolutions is, however, that for a wide range of system parameters, there is a time interval lasting up to several million years, during which the star is both very luminous ($`10^{38}10^{39}`$ ergs/s) and very hot, with $`kT`$ near $`10`$ eV. Characteristic evolutionary sequences include an epoch of helium burning, a shorter epoch during which heavier elements are burned, and, in some cases, a Type Ib or Type Ic, hydrogen-poor, supernova.

Typical temperatures are lower, and typical luminosities higher than those considered in §2.1. Another distinction is that the core may not be hottest and brightest immediately after the TD. This is because the further evolution of the He star will typically increase its luminosity and alter its temperature; in addition, winds that can block emergent soft X-rays tend to decrease with time.

Wolf-Rayet stars may be the closest well-observed analogs of the He-star remnants of TDs. They tend to eject such significant winds, however, that their appearance may be dominated by the interaction of this surrounding matter with radiation from the He star. Perhaps because of this, no known He star has a “supersoft” signature. It is possible, however, that some of the known SSSs, whose fundamental nature is not yet understood, could be He stars.

He-star remnants of TDs may be more likely than other He stars to be observed as SSSs. This is because the BH continues to act as a sort of vacuum cleaner for mass outside the tidal radius ($`R_t[M_{}/M_{bh}]^{\frac{1}{3}}`$). Furthermore, since self-generated winds tend to decrease with time (see, e.g., Woosley, Langer, & Weaver 1995), there will typically be less and less material within the tidal radius, even as nuclear evolution continues. Thus, He stars near MBHs, especially those that have undergone TD, may provide the cleanest signatures of the underlying properties of the star, distinct from the signatures of its interactions with its environment.

## 3 Observational Signatures

Luminosity: For both white dwarfs and He stars, the luminosity of the stripped core should not exceed that of typical giant stars. Luminosities in excess of $`10^{39}`$ erg/s must be due to some alternative or additional effect.

Spectrum: The maximum temperature expected for hot white dwarfs (He stars) is $`100`$ ($`15`$) eV. Therefore, when $`N_H`$ is low, the X-ray spectrum is dominated by soft photons. There can also be a modest contribution from photons at higher energies. First, the spectrum could differ from a simple thermal shape, e.g. due to opacity effects, as is well known in the case of white dwarf atmospheres (see, e.g., Dupree & Raymond 1983, van Teeseling et al. 1994, Hartmann & Heise 1997); thus, atmospheric effects must be included when fitting the X-ray spectrum. Second, there could be contributions from a corona, an infall disk within the stellar core’s Roche lobe, or shock heating of the ambient medium. Finally, it may not be possible to resolve the hot core and the MBH, which should be accreting matter at a low level, e.g., from stellar winds. Note that any radio emission is more likely to be originating from the region around the MBH than from the hot remnant of TD.

It is important to look for distinctions between the signature of tidally stripped remnants and other possible explanations for SSS signatures. One distinction, e.g., between SSSs in binaries and stripped cores, should be that the distinctive disk signatures of the former (Popham & Di Stefano 1996) should be missing in the latter.

lonization:  Even if the total flux from a stripped core over its lifetime (which, for the purposes of ionization signatures, may be defined to be the time interval during which $`T>50,000`$ K) is smaller than that from the accretion flare event, the number of high-energy ionizing photons it emits may be significantly larger. The distinctive signature is photoionization of high-energy states.

Time Variability:  Intrinsic variability may occur at many time scales. Here we concentrate on variability due to absorption of soft X-radiation by a varying amount of matter lying between the source and our detectors. Figure 1 demonstrates that the decline in count rate with increasing column density can be dramatic.

Figure 1: Chandra ACIS-I detector count rates of a black body radiator with a bolometric luminosity of 1038 erg/s in M 31 for different temperatures (solid lines) and depending on the intervening absorbing column.

If the stripped core orbits the MBH, the flux we receive may decrease as it travels behind gas in the disk around the MBH. If the mass and radius of the star were $`M_{}`$ and $`R_{}`$ before the TD, and if the orbital separation between the core and the MBH is presently equal to the separation at which the star filled its Roche lobe, then

$$P_{orb}=0.32years\eta \left[\frac{R_{}}{100R_{}}\right]^{\frac{3}{2}}\left[\frac{M_{}}{M_{}}\right]^{\frac{1}{2}}$$
(3)
If the core escapes from the vicinity of the MBH, it could also exhibit variability as it passes behind gas clouds. Making a set of simple assumptions, we estimate the amount of time, $`\tau ,`$ it takes for the hot core to pass behind a cloud with column density $`N_H,`$ located a distance $`a`$ from the massive black hole.111 Consider a cloud with mass $`m`$ located a distance $`a`$ from the MBH. If the cloud is corotating, the largest volume it could have is roughly equal to the volume of a sphere whose radius is the effective Roche lobe radius. This allows us to relate the mass of the cloud to its radius $`R`$, its distance $`a`$ from the MBH, and to $`M_{bh}`$. Assume the cloud is composed of hydrogen and that we can approximate an average value of $`N_H`$ by the number of hydrogen atoms divided by the square of the cloud radius. This allows us to relate $`R`$ to $`N_H.`$ Finally, we use $`\tau =R/v`$.

$$\tau =\frac{3.5days}{v_{150}}\xi \left[\frac{N_H}{N_{H,O}}\right]\left[\frac{a}{pc}\right]^3\left[\frac{3\times 10^7M_{}}{M_{bh}}\right],$$
(4)
where $`\xi `$ is of order unity, $`v_{150}`$ is the transverse speed in units of $`150`$ km/s, and $`N_{H,0}=1.9\times 10^{22}`$. The actual duration of obscuration may be shorter or longer, depending on the geometry of the cloud and the track of the hot remnant.

Doppler shift:  The orbital speed of the star in a circular orbit around a $`3\times 10^7M_{}`$ black hole, with orbital separation, $`a,`$ such that a $`100R_{}`$ star fills it Roche lobe, is $`5\%c;`$ the speed of the hot remnant could be greater. Thus, the most direct confirmation that an observed SSS is a tidally stripped remnant, would be a periodic Doppler shift compatible with an orbit small enough to allow TD of the progenitor star. Significant blue or red shifts could possibly be inferred from spectral fits. Even if the remnant does not orbit the MBH, it might be possible to measure its velocity and establish the likelihood that it is escaping from the central gravitational potential.

Spatial Distribution:  If the hot stellar core escapes from the region around the MBH , traveling along a direction perpendicular to our line of sight, it could traverse a considerable distance before it dims and cools. If the average speed is $`10^8`$ cm/s, the core could travel $`1`$ ($`0.1`$) kpc in $`10^6`$ ($`10^5`$) years. If the product of the TD rate and average hot-core lifetimes is larger than unity, we might expect to see several hot cores ringing the MBH, with the oldest cores further away. Any cores located $`1`$ kpc from the galaxy center would likely be He stars. It is interesting to note that the other two SSS candidates near the center of M31 are each several arcminutes from the galaxy’s center, produce count rates similar to that of the central SSS, and are each considerably softer than the central source.

Population Signatures:  (1) Consistency with the model:  If $`N`$ M31’s SSSs are remnants of TD events, each represents a much larger population of giants in the galaxy’s central region at the time of the TD. We must therefore test whether the present-day population could have evolved from one with the requisite number of giants some $`10^310^6`$ years ago. If we determine that a subset of the TD remnants are He stars, it is likely that there were a large enough number of centrally-located high-mass stars active in the recent past to lead to some supernovae. We can ask if the characteristics of the central region (up to $`1`$ kpc from the center) are consistent with the required number of Type II or Type Ib and Ic supernovae. The degree and characteristics of the ionization of the central region must also be consistent with the effects expected from the $`N`$ TD remnants and the requisite total population of high-mass stars. It may be instructive to consider the population of He stars near the center of the Milky Way, which is much better studied. This population is large and includes members with apparently unusual properties. (See, e.g., Tamblyn et al. 1996, Najarro et al. 1994 and references therein.) Unfortunately the large column density toward the center of the Milky Way makes it impossible to study any soft X-ray emission from these stars.

(2) Eliminating other possibilities:  Population studies of the region around the galactic center can lead to estimates of the numbers of post-AGB stars and also of X-ray binaries of various types. The physical processes that give rise to each of these objects also may give rise to SSSs that have nothing to do with TD. We therefore need: (1) observations of the central population, (2) theoretical estimates of the number of SSSs not due to TD expected, and (3) comparisons with studies of populations located away from the central region. These should allow us to estimate the probability that soft sources observed near the center have their genesis through processes that are ordinary and expected parts of stellar and binary evolution.

## 4 Other Galaxies

Radiation from a hot stripped core can be detected in other galaxies if (1) the count rate is high enough, and (2) the source can be resolved.

Count Rate:  The count rate from the soft X-ray source in the center of M31 (Garcia et al. 2000) was $`10`$ counts/ksec. Based on spectral fits to nearby (brighter) sources, the absorption was estimated to be $`2.8\pm 1.0\times 10^{21}`$ cm-2. The number of counts expected from an $`E`$ ksec exposure of a source with similar spectrum, located a distance $`D`$ from us, is approximately

$$N_\gamma =10counts\times E\times \left(\frac{Mpc}{D}\right)^2\times f\left(\frac{N_H}{2.8\times 10^{21}}\right).$$
(5)
$`N_H`$ is the column density between us and the new source. For face-on or early-type galaxies $`N_H`$ could be smaller than $`2.8\times 10^{21}`$, increasing the count rate in a way that depends on the spectrum, but which is likely to be described by a non-linear monotonically decreasing function, $`f(\frac{N_H}{2.8\times 10^{21}}).`$ (See Figure 1.)

A $`50`$ ksec exposure could therefore yield 20 counts from a hot stellar core an in M31-like galaxy as far away as 5 Mpc; for galaxies with less internal absorption, the count rate would be higher. Should the photons’ energies be clustered below 1 keV, we would have convincing evidence for a centrally located soft X-ray source. XMM-Newton, with its higher effective area, could achieve similar numbers of counts for galaxies up to 10 Mpc away. There are about 200 (600) galaxies within 5 (10) Mpc listed in Simbad and the Principal Galaxies Catalog (PGC). The galaxy type or orientation could allow the detection of soft central X-ray sources in roughly $`1/3`$ of these galaxies.

Spatial Resolution In M31 there are $`4`$ other X-ray sources within $`10^{\prime \prime }`$ of the nucleus. If the central density of X-ray sources in other galaxies is similar, we might not be able to resolve the central source in galaxies farther from us than $`23`$ Mpc. Note, however, that, if the soft central source is variable, as the source in M31 appears to be, it may be possible to successfully apply to X-ray searches for central sources, the same difference imaging techniques that have been applied to optical searches for supernovae and microlensing. If the hard X-radiation is steady, or varies in a way well-known from studying other X-ray stars, then we may be able to infer the existence of the soft central source. XMM-Newton, which has $`45`$ times Chandra’s effective area at low energies, will make such differencing studies possible for M31-type galaxies at least as far away as $`5`$ Mpc. This is particularly so if simultaneous observations with Chandra can resolve some sources and help to determine their spectral characteristics.

Statistical Tests:  Statistical questions can be posed when we study a population of galaxies. For example, are galaxies for which there is good evidence (perhaps from dynamical studies) for a MBH more likely to house a centrally located population of soft X-ray sources, than galaxies for which the dynamical studies do not indicate the presence of a MBH?

## 5 Conclusion

When a giant is tidally disrupted, the stripped core may provide a hot and bright signature for $`10^310^6`$ years. Detecting such cores and relating them to TD events is challenging. Fortunately, the present generation of instruments, including Chandra, XMM-Newton, and HST, promise opportunities for progress.

This new method of studying MBHs complements searches for flare events that may be due to TDs. The latter must of course be used for distant galaxies, and there is a growing body of evidence for luminous UV and X-ray events lasting for months to years, which may be consistent with the accretion of the envelopes of tidally-disrupted stars. 222Examples: the UV flare in NGC 4552 (Renzini et al. 1995) the $`5\times 10^{43}`$ erg/s outburst in IC 3599 (Brandt et al. 1995, Grupe et al. 1995), the $`2\times 10^{43}`$ erg/s outburst in NGC 5905 (Bade et al. 1996, Komossa & Bade 1999), the $`>\mathrm{\hspace{0.17em}9}\times 10^{43}`$ erg/s flare in RXJ 1242.6-1119 (Komossa & Greiner 1999), the $`10^{44}`$ erg/s flare in RXJ 1624.9+7554 (Grupe et al. 1999), the $`10^{44}`$ erg/s flare in RXJ 1331.9-3243 (Reiprich & Greiner 2000), or the $`2\times 10^{44}`$ erg/s flare in RX J1420.4+5334 (Greiner et al. 2000). Such events are so rare, however, that we would be lucky to detect one in the near future among the small number of galaxies within $`5`$ or $`10`$ Mpc. In nearby galaxies, however, the long-lived signatures of the stripped cores may allow us to infer that TDs have occurred. The two complementary modes of study should help us to better understand the frequency and consequences of tidal disruption events and to use them as probes of massive black holes in the center of galaxies.


format: text

text:
# 1 Introduction

## 1 Introduction

The Cosmic Microwave Background (CMB) was discovered by Penzias and Wilson  in 1965, and interpreted by Dicke et al. . This radiation comes from the first moments of our Universe, and was first predicted by Gamow, Alpher and Hermann  in 1948, in the context of the Big Bang theory. The Big Bang theory predicted that our Universe was in expansion, which was first observed by Hubble  in 1929, via the redshift of galaxies. As a consequence of the expansion, the Universe also cools down, meaning that at earlier times it was hotter. Going back in the past, our Universe was so hot that matter and photons were tightly coupled with each other. They formed a plasma in thermal equilibrium. As the Universe cooled down, the photon energy became too small to ionise the matter, below 0.1 eV (below 13.6 eV due to the high photon to electron ratio), the mean free path of the photons became larger than the horizon. The photons were free to cross the Universe towards our detectors. The blackbody dsitribution and spatial properties of these photons remain unchanged due to their negligeable cross-section with matter, the blackbody temperature is just cooled to 2.7 K due to the Universe expansion. Through this radiation, we obtain a picture of our Universe 300 000 years after the Big Bang.

## 2 Motivations

After the extremely accurate CMB spectrum measurement by the COBE satellite in 1989, proving that the radiation was a perfect blackbody , and the first detection of deviations  from the perfect homogeneity at the level of 10-5 K, most experiments have concentrated on improving the measurement of the anisotropies. In most theories and with consistent available data, the statistical distribution of the CMB anistropies is Gaussian. The anisotropies are therefore represented by their spherical harmonic power spectrum (non-Gaussianity is nevertheless looked for). The most accurate results are from COBE on large angular scale ($`>`$ 7 degrees) and on small angular scales ($`<`$ 2 degrees) from Boomerang , Maxima , and DASI  (see Figure 1). The intermediate angular scales, between COBE and the other experiments, lack measurements, due to the small sky coverage of these experiments. The main Archeops goal is to link COBE and Boomerang-Maxima-DASI angular scales, with high-sensitivity detectors covering a large fraction of sky (30% of the sky). Archeops is also a testbed for Planck-HFI, using the same type of cryogenic system and detectors; and it measures dust and galactic sources in a polarized frequency band.

## 3 Instrument and Scan-strategy

A complete description of Archeops can be found in  or on http://www.archeops.org, but the main characteristics of Archeops are:

* high angular resolution (10 arcminutes), due to a large primary mirror of 1.5 meters and to our horns, which guide the radiation to our bolometers;
* 30% sky coverage due to a scan-strategy consisting of making large circles on the sky. Associated with the first characteristic, this allows Archeops to measure the CMB power spectrum from angular scales of about 20 degrees (l=30) down to 10 arcminutes (l=800);
* bolometers cooled to 100 mK, by the same type of cryogenic system which will be used by Planck HFI, and does not rely on gravity;
* a frequency band, the 353 GHz band, polarized using OMT, is dedicated to the measurement of Galactic dust and point sources polarization.

The Archeops telescope is an off-axis gregorian telescope (see Figure 2), pointed at 41 degrees elevation, which defines the size of our circles on the sky at about 200 degrees.

To make circle on the sky the gondola spins at constant elevation, and as the earth rotates, the center of the circle describes also a circle on the sky. This describes a typical ring or donut-like sky coverage (see Figure 3). For Archeops (41 degrees elevation and 68 degrees of latitude), the sky coverage is about 30% for a 24 hour flight.

After the secondary mirror, the radiation goes through a polypropylene membrane (see Figure 4), and enters the horns, which define the part of the sky seen by each bolometer; they play a crucial role for the angular resolution of the experiment.

The detectors are spiderweb bolometers (Figure 5), which have a small calorific capacity so a fast response time, and are less sensitive to cosmic ray hits. These bolometers are cooled to 100 mK by a 3He-4He dilution, which is produced in capillary tubes placed around the focal plane (see Figure 5). The focal plane is filled with 22 bolometers, with the following frequency band distribution :

* 8 detectors at 143 GHz, where the CMB is the dominant emission at large galactic latitude;
* 6 detectors at 217 GHz, where the CMB is still dominant but dust contamination is larger;
* 6 detectors at 353 GHz, where the dust and atmospheric emission are dominant; these 6 channels are polarized using OMT : 3 pairs of 2 detectors share the same horn, the signal being separated into 2 orthogonal polarized components;
* 2 detectors at 545 GHz, where dust and atmospheric emission are dominant.

This frequency distribution is used to distinguish the different contributions of astrophysical origin or from parasitic signals, such as atmospheric emission; this is illustrated in Figure 6.

## 4 Flights

The Archeops experiment has already flown twice, once for a technical flight launched in Trapani (Sicily, Italy) crossing the mediterranean sea towards Spain in July 1999, producing 4 hours of night data, but with only 4 detectors. The second flight started from Esrange, a SSC (Swedish Space Corporation) base, near Kiruna (Lapland in the north of Sweden). The gondola was launched by the CNES (Centre National d’Etudes Spatiales***The French National Space Agency) on January, 29th 2001 at 2pm (local time) and landed in Russia at Syktyvkar around midnight.

The flight was shortened due to strong stratospheric winds (about 400km/h), which pushed the gondola too rapidly towards the east. We however obtained 7.5 hours of scientific data (during night and at a float altitude of 31.5 km), covering 22.7% of the sky. The cryostat remained below 100 mK during the entire flight.

The average size of our optical beam is around 8 arcminutes, the average effective beam size is around 12 arcminutes because of the bolometers time response. Some beams are also a bit asymetric, but this only impacts the very smallest angular scales, not attainable by Archeops for this flight (see Figure 8). The quality of the data was excellent, and will allow us to compute an accurate CMB anisotropies spectrum, especially on the low $`\mathrm{}`$ edge of the first acoustic peak.

## 5 Expected performance and future prospects

The data analysis is not yet finished, with various tests to check the power spectrum estimation still to be performed. Using the complete data processing pipeline and the known noise properties of our data, we have performed simulations for 10 detectors, based on the noise power spectrum (in the time domain) of our best channel. The power spectrum estimated from these simulations is shown in Figure 8. It is the average of power spectra given by each detector, weighted by their correlation matrices. In addition we performed an average over several simulations to check that the method was not biased.

The error bars are very small (smaller than available measurement) on large angular scale, because of our large sky coverage, and then increase after the first acoustic peak because of the small signal-to-noise ratio due to the small integration time. This is the motivation for a longer duration flight of at least 24 hours, as the signal-to-noise ratio increases nearly quadratically after the first 10 hours (because the sky coverage does not increase further). Two flights are therefore scheduled for next Winter (December 2001 - January 2002), from Esrange, with the same configuration. The expected performance for a 24 hour flight with 10 bolometers is shown in Figure 9.

The power spectrum estimated this way (Figure 9) is very accurate up to the end of the first peak and gives good constraints on the second one, allowing precise estimation of cosmological parameters.

## 6 Conclusions

Archeops is a balloon-borne experiment, designed to measure the CMB anistropies from large (about 10 degrees) to small angular scales (about 10 arcminutes). Archeops is also a testbed for Planck HFI, and will provide useful information on galactic dust and point source polarisation. The balloon has already flown twice successfully: once during a test flight and another time for its first scientific flight from Kiruna. The first scientific flight yielded 7.5 hours of scientific data, with 14 detectors dedicated to CMB measurements. Results from this flight will come soon, and we showed that it will give good constrains on large angular scale, but less constraining on small angular scale due to the short duration of the flight. Archeops will fly again this Winter, in order to increase the accuracy on large scales and measure the small scales.


text:
# ON ENERGETICS AND STRUCTURE OF SUB-PARSEC SCALE JETS IN QUASARS

## 1 Introduction

Observational data imply that jets in radio-loud quasars can extend over 8 decades of distance, from milliparsecs to hundreds of kiloparsecs. Such jets are launched in the vicinity of super-massive black holes and transport energy at a rate which is sometimes comparable with the accretion luminosity, i.e. $`10^{45}10^{47}`$ erg s<sup>-1</sup> Jets are relativistic and therefore Doppler boosting often makes them observable only on one side of the quasar. However, the hot spots where jets terminate are usually seen on both sides of the central source, and emission on both sides shows similar flux. This suggests that the radiation from hot spots is nearly isotropic, and indicates that jets are “light”, i.e. have lower mass densities than the external medium. The lightness is necessary to explain the nonrelativistic speeds of hot spots and formation of extended radio lobe structures.  Using the asymmetry of radio brightness of jets and counter-jets, Wardle and Aaron  derived the bulk Lorentz factor $`\mathrm{\Gamma }3`$. Larger Lorentz factors are suggested by X–ray observations ; however, because X–ray sensitivity is too low to observe counter-jets, such estimates are not direct and depend on the details of models of the X–ray production.

Quasar jets are imaged in the radio band down to parsec scales. There, they join the central cores, which in turn, due to the synchrotron self-absorption, produce very flat spectra, with the energy spectral index $`\alpha <0.5`$ ($`\alpha :F_\nu \nu ^\alpha `$) (see, e.g., review by Zensus.  In quasars with jets oriented close to the line of sight, the flux of the radio cores strongly dominates over that from radio lobes, and the total radio spectra are flat. Because of this, core-dominated radio quasars are often called FSRQs (flat spectrum radio quasars), as opposed to the lobe dominated radio quasars, which have steep radio spectra and are called SSRQs (steep spectrum radio quasars). Radio cores in FSRQs are often variable on monthly time scales, and variability of the radio flux is often accompanied by an appearance of ejecta propagating with apparently superluminal velocities.

The physics of the parsec scale jets (energetics, magnetic fields, particle distribution) is determined by the use of images, spectra, variability, and polarization properties. In particular, the angle of linear polarization provides information about magnetic field orientation relative to the jet axis. Recent data suggest that such orientation changes from being perpendicular to the jet closer to the central source to the roughly parallel at somewhat larger distances.  This may indicate that at smaller distances, the jets are free (ballistic) and that their energy is dissipated mainly in perpendicular shocks resulting from collisions of the jet inhomogeneities moving with different radial velocities. At larger distances, on the other hand, most dissipation is likely to come from oblique, reconfinment shocks produced due to interaction of a jet with external medium. For a number of objects, there are also circular polarization measurements available.  If circular polarization is produced via the Faraday conversion mechanism, it provides information on the number of electrons at mildly relativistic energies. Combined with the information about the total power of a jet, this can be used to estimate the pair content. It should be emphasized that no information about number of mildly relativistic electrons is available from direct radio flux observations, because the synchrotron radiation by mildly relativistic electrons is self-absorbed.

The parts of jets located close to the cores can be studied also in other spectral bands besides radio, despite the fact that they are by many orders too compact to be resolved. In most FSRQs, typical thermal components of the quasar emission – such as the UV radiation from an accretion disc, the X–ray radiation from the disc corona, and the infrared radiation from dust – are overshined by the Doppler-boosted nonthermal jet radiation. Variability of this radiation is measured on time scales often shorter than those observed in the radio band, and this strongly suggests its parsec/sub-parsec origin. High polarization in the optical and IR bands implies synchrotron mechanism, while the distinct high energy components, with luminosity peaks in the MeV – GeV range,  are presumably products of Comptonization of broad emission lines and IR radiation of hot dust.  If the EC (external-Compton) scenario is indeed responsible for the high energy components, then observations of the X–ray spectra can provide an exceptional opportunity to study the population of the lower energy population of the relativistic electrons directly.  Furthermore, due to narrowness of the electromagnetic spectrum bands covered by dominant broad emission lines and by infrared radiation of hot dust, any features such as breaks in the energy distribution of electrons should be sharply imprinted in the EC radiation components. In particular, as it will be demonstrated below, the shapes of the observed high energy spectra are consistent with the two-power-law injection function of electrons, and the latter can result from a two-step acceleration process of electrons as discussed below.

The two component nonthermal spectra are also observed in BL Lac objects, which differ from FSRQs as having very weak or undetectable emission lines. FSRQs and BL Lac objects form together a class of objects called blazars. It is shown by Fossati et al that broadband spectra of blazars form a sequence which can be parametrized by their total luminosities. In this sequence, FSRQs are the most luminous objects. Both their low and high energy spectral components appear to be least extended to the high energies, and their $`\gamma `$–ray luminosities during flares strongly dominate over synchrotron luminosities. The least luminous blazars are represented by the X–ray selected BL Lac objects. Their synchrotron spectra extend up to hard X–rays, and the $`\gamma `$–ray spectra reach TeV energies. $`\gamma `$–ray luminosities in the TeV-emitting BL Lac objects usually do not dominate over synchrotron luminosities. Low luminosity BL Lac objects are probably associated with nuclei of radio galaxies accreting at a low accretion rate, and $`\gamma `$–rays are very likely produced by the SSC (synchrotron-self-Compton) process. In the more luminous, radio selected BL Lac objects the EC process can be dominant,  in similarity to FSRQs.

In this presentation we focus on jets in quasars, where the interaction of the jet with external radiation is significant and the presence or lack of radiative effects of this interaction can be used to determine physical parameters and constraints on structure of jets on sub-parsec scales. This paper is organized as follows: in §2, we analyze the energetics of intrinsic collisions of material in the jet; in §3, we present the motivation for introduction of the two-step stochastic acceleration of electrons; in §4, we discuss the formation of spectral breaks in the high energy spectral components; in §5, §6, and §7, respectively, we derive the magnetic field intensity, the pair content, and the average power of a jet. In §8, we formulate our predictions regarding the production of soft X–ray precursors via Comptonization of the external radiation by the cold electrons in a jet, and list the main conclusions in §9.

## 2 Energy dissipation

Cooling time scales of ultrarelativistic electrons, which in quasar jets produce synchrotron radiation in the optical band and inverse-Compton radiation in the MeV – GeV range, is much shorter than the dynamical/propagation time scale. Therefore, such electrons must be accelerated in situ, in sites where the jet loses part of its energy and the shocks are formed. The jet energy can be dissipated in internal shocks – formed in collisions between inhomogeneities in a jet  or following reconnection of magnetic fields,  – and/or in external shocks, formed via reconfinment of the jet by external medium  or due to collisions of a jet with external clouds . Perpendicular orientation of magnetic fields to the jet axis, inferred from the polarization measurements in the optical band  and at the high radio frequencies  favors the scenario where at sub-parsec distances, the particles are accelerated in transverse shocks which are likely to be produced via collisions of inhomogeneities flowing down the jet with different radial velocities. A given collision event is then responsible for a given flare, and the time scale of such a collision, as measured in the comoving frame of the shocked plasma, is

$$t_{coll}^{}t_{fl}𝒟,$$
(1)
where $`t_{fl}`$ is the observed time scale of the flare, and

$$𝒟\frac{1}{\mathrm{\Gamma }(1\beta \mathrm{cos}\theta _{obs})}$$
(2)
is the Doppler factor of the shocked plasma. The amount of energy dissipated during the collision is

$$E_{diss}^{}\frac{L_{fl}^{}t_{coll}^{}}{\eta _{rad}\eta _e}\frac{L_{fl}t_{fl}}{f𝒟^3\eta _{rad}\eta _e},$$
(3)
where $`\eta _e`$ is the fraction of dissipated energy used to accelerate electrons, $`\eta _{rad}`$ is the (average) radiative efficiency of electrons, and

$$L_{fl}4\pi \left(\frac{L_{fl}}{\mathrm{\Omega }}\right)=4\pi \left(\frac{L_{fl}^{}}{\mathrm{\Omega }^{}}\right)𝒟^4=L_{fl}^{}f𝒟^4$$
(4)
is the apparent luminosity. All primed quantities are as measured in the source comoving frame. The factor $`f`$ expresses the additional dependence of the apparent luminosity on the observed angle, if the intrinsic radiation is anisotropic. For the EC process, $`f(𝒟/\mathrm{\Gamma })^2`$, while for synchrotron and SSC radiation $`f1`$

In the above relations it was assumed that during a collision, the shocked plasma moves with a constant bulk Lorentz factor. In general this is not true: depending on the initial parameters of inhomogeneities, the shock can accelerate or decelerate, and during the collision the shock structure can evolve from a double one to a single shock (forward or reverse). However, if inhomogeneities are intrinsically identical, then the forward-reverse shock structure is symmetric in the shocked plasma frame (= discontinuity contact surface frame), and the bulk Lorentz factor of the shock is constant. This case is adopted in our quantitative analysis, not only for simplicity, but also because it provides the highest efficiency of energy dissipation.

Assuming that the inertia of inhomogeneities is dominated by protons (i.e. $`n_e/n_pm_p/m_e`$), and noting that

$$E_{diss}\eta _{diss}(E_1+E_2)=\mathrm{\Gamma }E_{diss}^{},$$
(5)
where $`E_1(N_p/2)\mathrm{\Gamma }_1m_pc^2`$ and $`E_2(N_p/2)\mathrm{\Gamma }_2m_pc^2`$ are energies of inhomogeneities, we find that total number of protons involved in the collision is

$$N_p\frac{2\mathrm{\Gamma }E_{diss}^{}}{\eta _{diss}(\mathrm{\Gamma }_1+\mathrm{\Gamma }_2)m_pc^2}\frac{2L_{fl}t_{fl}}{\kappa f𝒟^3\eta _{rad}\eta _em_pc^2},$$
(6)
where

$$\kappa \frac{E_{diss}^{}}{N_pm_pc^2}=\frac{\eta _{diss}(\mathrm{\Gamma }_1+\mathrm{\Gamma }_2)}{2\mathrm{\Gamma }}$$
(7)
is amount of energy dissipated per proton in units of $`m_pc^2`$. For $`\mathrm{\Gamma }_2>\mathrm{\Gamma }_11`$, energy and momentum conservations give

$$\mathrm{\Gamma }\sqrt{\mathrm{\Gamma }_1\mathrm{\Gamma }_2}$$
(8)
and

$$\eta _{diss}\frac{((\mathrm{\Gamma }_2/\mathrm{\Gamma }_1)^{1/2}1)^2}{(\mathrm{\Gamma }_2/\mathrm{\Gamma }_1)+1},$$
(9)
and then

$$\kappa \frac{((\mathrm{\Gamma }_2/\mathrm{\Gamma }_1)^{1/2}1)^2}{2(\mathrm{\Gamma }_2/\mathrm{\Gamma }_1)^{1/2}}.$$
(10)
For $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_1=2.5`$, which corresponds with $`\eta _{diss}0.1`$, $`\kappa 0.11`$.

## 3 Particle acceleration

As it was demonstrated by Sikora and Madejski,  energy flux in powerful jets in quasars cannot be dominated by pair plasma. This is because such jets would produce much larger fluxes of soft X–ray radiation than is observed. Hence, we assume that the inertia of jets is dominated by protons. In this case, the structure of shocks and the structure of the turbulence around them are both determined by protons. Provided that the magnetic fields are amplified up to equipartition with protons, the protons are effectively accelerated both by 1st order and by 2nd order Fermi process.  Those which reach energies $`>10^9`$ GeV, interact efficiently with ambient photons and trigger (mainly via photo-meson process) the synchrotron-supported pair cascades.  However, such a model fails to reproduce very hard X–ray spectra of FSRQs.  This may indicate that too few protons reach sufficiently high energies to power pair cascades.

Alternative possibility is that the high energy radiation in FSRQs is produced by directly accelerated electrons. However, for the acceleration of electrons to occur via the Fermi process, these electrons must be first preheated/preaccelerated up to energies $`\gamma _b`$, at which point, the magnetic rigidity of electrons becomes comparable with magnetic rigidity of thermal protons, i.e. when their momenta are equal:

$$m_e\sqrt{\gamma _b^21}m_p\sqrt{\gamma _{p,th}^21},$$
(11)
where

$$\gamma _{p,th}1=\eta _{p,th}\kappa $$
(12)
is the average thermal proton energy in the shocked plasma, and $`\eta _{p,th}`$ is the fraction of the dissipated energy tapped to heat the protons. For the reasonable assumption that $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_1<10`$, $`\gamma _{p,th}1<1`$, i.e. the thermal proton plasma is at most mildly relativistic.

Noting that the average energy of injected electrons is

$$\overline{\gamma }_{inj}1\frac{n_pm_p}{n_em_e}\eta _e\kappa =\frac{n_pm_p}{n_em_e}\frac{\eta _e}{\eta _{p,th}}(\gamma _{p,th}1),$$
(13)
we find that for $`\gamma _b`$ and $`\overline{\gamma }_{inj}1`$

$$\frac{\gamma _b}{\overline{\gamma }_{inj}}\frac{n_e}{n_p}\frac{\eta _{p,th}}{\eta _e}\left(\frac{\gamma _{p,th}+1}{\gamma _{p,th}1}\right)^{1/2}.$$
(14)
Hence, for non- or mildly relativistic shocks and for $`\eta _{p,th}\eta _e`$, the threshold energy for the Fermi acceleration of electrons significantly exceeds the average energy of electrons, even if $`n_e=n_p`$. That implies that the often-considered bulk preheating process is not able to provide an adequate number of electrons with $`\gamma \gamma _b`$. But, as it was recently demonstrated in numerical PIC (particle-in-cell) simulations, the preheating/preacceleration mechanism has a stochastic character and a significant fraction of electrons can reach such energies.  Furthermore, the fact that in FSRQs, the X–ray spectral index $`\alpha _X`$ is never negative, implies that most of electrons occupy the lowest energies, i.e. the median Lorentz factor is $`\overline{\gamma }_{inj}`$. Therefore, this scenario does not predict a thermal bump in the preaccelerated electron energy distribution. Instead, the preacceleration mechanism, in similarity to the Fermi process, produces a power-law energy distribution of electrons. It should be emphasized here that because the X–rays in FSRQs are produced in the slow cooling regime (see next Section), electrons with low energies are not populated by cooling effect, but are directly injected with such energies.

## 4 Electromagnetic spectra

The basic feature of the high energy spectra in FSRQs — a spectral break between the X–ray and the $`\gamma `$–ray bands — has a natural explanation in terms of the EC model. In this model, X–ray spectra are produced by electrons with radiative cooling time scale $`t_{cool}^{}`$, longer than the collision time scale $`t_{coll}^{}`$ (slow cooling regime), whereas $`\gamma `$–rays are produced by electrons with $`t_{cool}^{}<t_{coll}^{}`$ (fast cooling regime). Noting that the angle averaged cooling rate of electrons, dominated by Comptonization of external radiation, is

$$|\dot{\gamma }|\frac{\sigma _T}{m_ec}u_{ext}^{}\gamma ^2,$$
(15)
we obtain that the angle averaged cooling time scale is

$$t_{cool}^{}\frac{\gamma }{|\dot{\gamma }|}\frac{m_ec}{\sigma _T}\frac{1}{\gamma u_{ext}^{}},$$
(16)
where $`u_{ext}^{}`$ is the energy density of an external radiation field. Then, from $`t_{cool}^{}=t_{coll}^{}`$, where $`t_{coll}^{}`$ is given by Eq. (1), the break in the electron distribution is at energy

$$\gamma _c\frac{m_ec}{\sigma _T}\frac{1}{u_{ext}^{}t_{fl}𝒟}.$$
(17)
For $`\gamma <\gamma _c`$, the slope of the electron distribution is the same as the slope of the injection function; for $`\gamma >\gamma _c`$, the slope of the electron energy distribution is steeper by $`\mathrm{\Delta }s=1`$ ($`s:N_\gamma \gamma ^s`$).

Since time scales of flares in FSRQs are rarely shorter than 1 day, the distances from the central source where they are produced,

$$r_{fl}(r_{fl}/\mathrm{\Delta }r_{coll})ct_{fl}𝒟\mathrm{\Gamma },$$
(18)
are expected to be larger than 0.1 parsec. At such distances, contribution to $`u_{ext}^{}`$ is dominated by the diffuse components of the broad emission line light and of the infrared radiation of hot dust. Noting that

$$u_{ext}^{}=\frac{1}{c}I_{ext}^{}𝑑\mathrm{\Omega }^{}=\frac{1}{c}I_{ext}𝒟_{in}^2𝑑\mathrm{\Omega }u_{diff}\mathrm{\Gamma }^2$$
(19)
where

$$𝒟_{in}=\frac{1}{\mathrm{\Gamma }(1\beta \mathrm{cos}\theta _{in})}$$
(20)
and $`\theta _{in}`$ is the angle between the photon direction and the jet axis, we predict that the break in an electron energy distribution at $`\gamma _c`$ should be imprinted in the electromagnetic spectrum at frequency

$$\nu _{EC,c}𝒟^2\gamma _c^2\nu _{diff}\left(\frac{m_ec^2}{\sigma _T}\right)^2\frac{\nu _{diff}}{u_{diff}^2t_{fl}^2\mathrm{\Gamma }^4}$$
(21)
and that the spectrum should change the slope around $`\nu _{EC,c}`$ by $`\mathrm{\Delta }\alpha _{X,\gamma }0.5`$. For typical energy densities expected in the BELR (broad emission line region) as well as that provided by the hot dust radiation field, the break is predicted to be located in the range of $`10^{20}10^{22}`$ Hz (Sikora et al.), in agreement with observations.

In most FSRQs, during high energy flares the slope of $`\gamma `$–ray spectra measured by EGRET  is $`\alpha _\gamma 1`$, while the X–ray spectra show slopes $`\alpha _X0.60.7`$ The smaller than the predicted value $`0.5`$ of $`\mathrm{\Delta }\alpha _{X,\gamma }`$ can be explained by taking into account contribution to the soft and mid-X–ray bands of the SSC component . More challenging are the so called “MeV blazars.” In these objects, the spectra peak in the $`300`$ keV - $`30`$ MeV range and the $`\gamma `$–ray spectra are very soft ($`\alpha _\gamma >1.4`$), while the X–ray spectra are very hard ($`\alpha _X0.5`$).  The values of $`\mathrm{\Delta }\alpha _{X,\gamma }1.0`$ measured in those sources are much too large to be explained solely by the effect of cooling, and instead, it is necessary to postulate an additional break in the electron energy distribution besides $`\gamma _c`$. It is tempting to identify such a break with the break at $`\gamma _b`$, predicted by the two-step acceleration process discussed in the Section 3. External photons scattered by electrons with energies $`\gamma _b`$, are boosted to energies

$$\nu _{EC,b}𝒟^2\gamma _b^2\nu _{diff}.$$
(22)
For $`2.5<\mathrm{\Gamma }_2/\mathrm{\Gamma }_1<10`$ and $`\eta _{p,th}0.5`$, the break in the electron injection function is $`600<\gamma _b<1700`$ (see Eqs. 10-12 and Table 1), and should appear in the electromagnetic spectrum at around $`1`$ GeV if $`\nu _{diff}=\nu _{BEL}`$, and at lower energies (by a factor of $`1030`$) if $`\nu _{diff}=\nu _{IR}`$. In the former case $`\nu _{EC,b}\nu _{EC,c}`$ and the $`\gamma `$–ray spectra in the EGRET range should have slopes $`\alpha _\gamma 1`$, while in the latter case $`\nu _{EC,b}\nu _{EC,c}30`$ MeV, and the $`\gamma `$–ray spectra in the EGRET range should be soft. Hence, we conclude that for blazars producing most of their radiation in the GeV band (and hereafter called GeV blazars) the $`\gamma `$–ray flux is probably produced at smaller distances (closer to the central source), where the production of $`\gamma `$–rays is dominated by Comptonization of broad emission lines, while the MeV blazar phenomenon is likely to be produced at larger distances, where production of $`\gamma `$–rays is dominated by the Comptonization of infrared radiation from dust.

The above scheme can also explain two other differences between MeV blazars and GeV blazars. One, already mentioned above, is that in MeV blazars the X–ray spectra should be much harder than in GeV blazars. This can result from the lower contribution of the SSC component to the X–ray band in MeV blazars as compared to GeV blazars (Sikora et al., in preparation). The other difference is that in the spectra of MeV blazars, in contrast to GeV blazars, the thermal UV bump is quite prominent.  This difference can be explained noting that magnetic fields are weaker at larger distances, and therefore in MeV blazars the synchrotron spectra are shifted to lower frequencies, revealing the UV bump.

## 5 Magnetic fields

Electrons with energies from the range $`[\gamma ;\gamma +d\gamma ]`$ produce synchrotron radiation at the rate

$$L_{\nu _{syn}^{}}d\nu _{syn}^{}(N_\gamma d\gamma )|\dot{\gamma }_{syn}|m_ec^2(N_\gamma d\gamma )c\sigma _T\gamma ^2u_B^{}$$
(23)
and the angle averaged EC radiation at the rate

$$L_{\nu _{EC}^{}}d\nu _{EC}^{}(N_\gamma d\gamma )|\dot{\gamma }_{EC}|m_ec^2(N_\gamma d\gamma )c\sigma _T\gamma ^2u_{ext}^{}$$
(24)
where

$$\nu _{syn}^{}(2e/3\pi m_ec)\gamma ^2B^{}$$
(25)
$$\nu _{EC}^{}\gamma ^2\nu _{ext}^{}$$
(26)
and $`u_B^{}=B_{}^{}{}_{}{}^{2}/8\pi `$ is the magnetic energy density. Noting that $`L_{\nu _{syn}}d\nu _{syn}=L_{\nu _{syn}^{}}d\nu _{syn}^{}𝒟^4`$ and $`L_{\nu _{EC}}d\nu _{EC}=L_{\nu _{EC}^{}}d\nu _{EC}^{}(𝒟^6/\mathrm{\Gamma }^2)`$ we find that the ratio of synchrotron to EC luminosity, both produced by the same population of electrons in the same energy range, radiating in the optically thin regime, is

$$_{syn/EC}\frac{\nu _{syn}L_{\nu _{syn}}}{\nu _{EC}L_{\nu _{EC}}}\frac{u_B^{}}{u_{ext}^{}}\left(\frac{\mathrm{\Gamma }}{𝒟}\right)^2.$$
(27)
From Eqs. (17) and (21), we have

$$u_{ext}^{}\frac{m_ec}{\sigma _T}\frac{1}{t_{fl}(\nu _{EC,c}/\nu _{diff})^{1/2}}$$
(28)
and inserting this into Eq. (27), we obtain

$$u_B^{}u_{ext}^{}_{syn/EC}(𝒟/\mathrm{\Gamma })^2\frac{m_ec}{\sigma _T}\frac{_{syn/EC}(𝒟/\mathrm{\Gamma })^2}{t_{fl}(\nu _{EC,c}/\nu _{diff})^{1/2}}.$$
(29)
For $`t_{fl}3`$ day, $`_{syn/EC}=0.1`$, $`\nu _{EC,c}=10`$ MeV and $`\nu _{diff}=10`$ eV, and the observer located at $`\theta _{obs}=1/\mathrm{\Gamma }`$, Eq. (29) gives $`B^{}\sqrt{8\pi u_B^{}}0.6`$ Gauss.

We can estimate the energy density of the relativistic electrons and determine whether they are in equipartition with magnetic fields. Relations $`E_{rad}^{}=\eta _{rad}\eta _eE_{diss}^{}`$ and $`E_e^{}=(1\eta _{rad})\eta _eE_{diss}^{}`$ give

$$E_e^{}=\frac{1\eta _{rad}}{\eta _{rad}}E_{rad}^{},$$
(30)
where

$$E_{rad}^{}L_{fl}^{}t_{coll}^{}\frac{L_{fl}t_{fl}}{𝒟^3(𝒟/\mathrm{\Gamma })^2}.$$
(31)
Hence, energy density of the relativistic electrons is

$$u_e^{}=\frac{E_e^{}}{\pi a^2\mathrm{\Delta }_{sh}^{}}\frac{1\eta _{rad}}{\eta _{rad}}\frac{L_{fl}}{2\pi c^3(a/a_\mathrm{\Gamma })^2t_{fl}^2g𝒟^6(𝒟/\mathrm{\Gamma })^2},$$
(32)
where $`aa_\mathrm{\Gamma }=ct_{fl}𝒟`$ is the cross-sectional dimension (radius) of a shocked plasma sheet, and $`\mathrm{\Delta }_{sh}^{}=2ct_{fl}𝒟g`$ is the thickness of the sheet. The scaling factor $`g`$, which relates the thickness of the shocked plasma sheet to the time scale of the flare/collision, depends on $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_1`$ and on the adiabatic index. The analytical form of $`g`$ is derived and presented in Moderski et al. (in preparation), and the example values of $`g`$ are shown in the Table 1.

Combining Eqs. (29) and (32), we obtain

$$\frac{u_B^{}}{u_e^{}}\frac{2\pi m_ec^4}{\sigma _T}\frac{\eta _{rad}}{1\eta _{rad}}\frac{g(a/a_\mathrm{\Gamma })^2_{syn/EC}}{(\nu _{EC,c}/\nu _{diff})^{1/2}}\frac{t_{fl}}{L_{fl}}𝒟^6(𝒟/\mathrm{\Gamma })^4.$$
(33)
The equation demonstrates very strong dependence of $`u_B^{}/u_e^{}`$ on $`\mathrm{\Gamma }`$ and on the observed angle $`\theta _{obs}`$. Provided that blazars are observed at $`\theta _{obs}1/\mathrm{\Gamma }`$, the equipartition condition, $`u_B^{}u_e^{}`$, implies

$$\mathrm{\Gamma }\left(\frac{\sigma _T}{2\pi m_ec^4}\frac{(\nu _{EC,c}/\nu _{diff})^{1/2}}{_{syn/EC}g(a/a_\mathrm{\Gamma })^2}\frac{1\eta _{rad}}{\eta _{rad}}\frac{L_{fl}}{t_{fl}}\right)^{1/6}.$$
(34)
For flares with luminosities $`L_{fl}10^{4849}`$ergs s<sup>-1</sup> and time scales $`t_{fl}3`$d, and assuming $`a=a_\mathrm{\Gamma }`$, $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_1=2.5`$, and $`\eta _{rad}=0.5`$, we obtain $`\mathrm{\Gamma }1826`$.

## 6 Pair content

The pair content is

$$\frac{n_+}{n_p}=(\frac{n_e}{n_p}1)/2$$
(35)
where

$$\frac{n_e}{n_p}=\eta _e\frac{m_p\kappa }{m_e(\overline{\gamma }_{inj}1)}$$
(36)
$$\overline{\gamma }_{inj}=\frac{_1^{\gamma _{max}}\gamma Q𝑑\gamma }{_1^{\gamma _{max}}Q𝑑\gamma }$$
(37)
and $`Q`$ is the electron injection function. For a double power-law injection function, $`Q\gamma ^{p_L}`$ for $`1<\gamma <\gamma _b`$, and $`Q\gamma ^{p_H}`$ for $`\gamma _b<\gamma <\gamma _{max}`$, and observationally suggested indices, $`1<p_L<2`$ and $`p_H>2`$,

$$\overline{\gamma }_{inj}\frac{(p_L1)(p_Hp_L)}{(2p_L)(p_H2)}\gamma _b^{2p_L}.$$
(38)
For $`p_L=1.8`$ ($`\alpha _X=(p_L1)/20.45`$ and $`\alpha _\gamma =p_L/20.9`$ in the EGRET band of GeV blazars), and $`p_H=2.8`$ ($`\alpha _\gamma =p_H/2=1.4`$ in the EGRET band of MeV blazars), and $`2.5<\mathrm{\Gamma }_2/\mathrm{\Gamma }_1<10`$, the number of pairs per proton is $`2.4<n_+/n_p<15.5`$ (see Table 1).

## 7 The power of a jet

Energy carried by the cold inhomogeneities prior to the collision is (see Eq. 6)

$$E_i\frac{N_p}{2}m_pc^2\mathrm{\Gamma }_i\frac{L_{fl}t_{fl}}{\kappa f𝒟^3\eta _{rad}\eta _e}\mathrm{\Gamma }_i,$$
(39)
where $`i=1,2`$. The time it takes to eject the inhomogeneities from the central source is

$$t_{ej,i}\frac{\lambda _0}{c\beta _i\mathrm{\Gamma }_i}\frac{𝒟}{\mathrm{\Gamma }_i}g_0t_{fl}$$
(40)
where $`\lambda _0=g_0c𝒟t_{fl}`$ is the proper length of an inhomogeneity (the thickness of a shell segment), and $`g_0`$ is the factor depending on $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_1`$ and on the adiabatic index of the shocked plasma (see Table 1). The time-averaged power of a jet is

$$P_j\mathrm{\Psi }\frac{E_1+E_2}{t_{ej,1}+t_{ej,2}}\mathrm{\Psi }\frac{N_pm_pc^3\mathrm{\Gamma }^2}{\lambda _0}\mathrm{\Psi }\frac{L_{fl}}{2\eta _e\eta _{rad}\kappa g_0}\frac{\mathrm{\Gamma }^2}{f𝒟^4},$$
(41)
where $`\mathrm{\Psi }`$ is the duty cycle factor. For $`L_{fl}10^{48}`$erg s<sup>-1</sup>, $`𝒟\mathrm{\Gamma }15`$, $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_12.5`$, $`\eta _e\eta _{rad}0.5`$ and $`\mathrm{\Psi }0.1`$, Eq. (41) gives $`P_j3\times 10^{46}`$ erg s<sup>-1</sup>.

## 8 X–ray precursors

Prior to the collision, electrons in inhomogeneities are cold, but being pulled by streaming protons they scatter external photons and produce the beamed bulk-Compton radiation.  The spectrum of this radiation peaks around

$$\nu _i𝒟_i\mathrm{\Gamma }_i\nu _{diff}$$
(42)
and has luminosity

$$L_{X,i}\delta N_{e,i}|\dot{ϵ}_{e,i}|𝒟_i^4,$$
(43)
where

$$\delta N_{e,i}=(N_e/2)\mathrm{Min}[1;r_{BC}/\lambda _0𝒟_i],$$
(44)
$$N_e\frac{E_{e,inj}^{}}{\overline{\gamma }_{inj}m_ec^2}\frac{\eta _eE_{diss}^{}}{\overline{\gamma }_{inj}m_ec^2}\frac{L_{fl}t_{fl}}{f𝒟^3\eta _{rad}\overline{\gamma }_{inj}m_ec^2},$$
(45)
$$|\dot{ϵ}_{e,i}|\mathrm{\Gamma }_i^2c\sigma _Tu_{ext},$$
(46)
$$𝒟_i\frac{1}{\mathrm{\Gamma }_i(1\beta _i\mathrm{cos}\theta _{obs})},$$
(47)
$`r_{BC}`$ is the distance within which most of the bulk Compton radiation is produced. The factor Min$`[1;r_{BC}/\lambda _0𝒟_i]`$ in Eq. (44) expresses the fact that in general, due to light travel effects and the finite length of the source, not all electrons contribute to the radiation observed at a given instant.

Bulk-Compton flares from cold inhomogeneities are expected to precede the nonthermal flares by

$$\delta t_i\frac{r_{fl}}{c\mathrm{\Gamma }_i𝒟_i}\frac{𝒟\mathrm{\Gamma }}{𝒟_i\mathrm{\Gamma }_i}t_{fl}\frac{\mathrm{\Delta }r_{coll}}{r_{fl}}.$$
(48)
Above formulae are very approximate, because they do not take into account that the observed bulk Compton radiation at any given instant includes contributions from a finite distance range over which $`u_{ext}`$ can change significantly. However, making the conservative assumption that the jet is developed (accelerated and collimated) not earlier than at a distance where electron energy losses are already dominated by Compton interactions with the broad emission line flux rather than the direct radiation from an accretion disc, we can calculate the minimum bulk Compton luminosity, using an approximation that $`u_{ext}u_{BEL}`$ const for $`r<r_{BEL}`$ and 0 for $`r>r_{BEL}`$. Noting that $`r_{BEL}3\times 10^{17}\sqrt{L_{UV,46}}`$ cm and $`u_{BEL}0.03`$ erg cm<sup>-3</sup> (Peterson ), we find that for $`\theta _{obs}=1/\mathrm{\Gamma }=1/15`$, $`\mathrm{\Gamma }_2/\mathrm{\Gamma }_1=2.5`$, $`L_{fl}=10^{48}`$ erg s<sup>-1</sup> and $`t_{fl}=3`$ day, the precursor produced by the faster-moving inhomogeneity should have luminosity $`L_{X,2}3\times 10^{46}`$ ergs s<sup>-1</sup>; should peak at around h$`\nu _23.2`$ keV; and should precede the nonthermal flare by $`\delta t_20.7t_{fl}`$. The precursor produced by the slower moving inhomogeneity would have luminosity $`L_{X,1}5\times 10^{45}`$ erg s<sup>-1</sup>; would peak around h$`\nu _11.8`$ keV; and would precede the nonthermal flare by $`\delta t_11.8t_{fl}`$. Since typical soft X–ray luminosities are of the order $`10^{46}`$ erg s<sup>-1</sup>, we can infer from the above estimates that jets propagating through much denser quasi-isotropic radiation fields than would be provided by the BELR are predicted to produce much stronger soft/mid X–ray radiation than is observed. This puts severe constrains on the jet structure at $`10^210^3`$ gravitational radii from the BH, where quasi-isotropic radiation is strongly dominated by the very dense direct radiation of an accretion disc and its corona. Hence, the lack of very large soft X–ray excesses suggests that jets near their bases are very wide (meaning that they have relatively large opening angle) and/or are slower than at larger distances.

## 9 Summary

Multi-wavelength observations of blazars provide a valuable tool allowing us to learn about the structure and physical properties of sub-parsec-scale jets in quasars. The observations suggest:

$``$ Rapid flares and transverse magnetic fields inferred from polarization data $``$ Streams of matter in small-scale jets are unsteady, and collisions between inhomogeneities propagating down the jet with different radial velocities can easily form transverse shocks;

$``$ Large powers of jets and relatively low soft/mid X–ray luminosities $``$ Inertia of jets is dominated by protons and the number of pairs per proton is $`<15`$;

$``$ Very hard X–ray spectra $``$ X–rays and $`\gamma `$–rays are produced by electrons accelerated directly, just as is the case for the synchrotron radiation, rather than being injected following the proton-induced pair cascades;

$``$ Differences between properties of GeV and MeV blazars $``$ EC model can explain these differences assuming a two-power-law electron injection function. The break in the injection function can be related to the threshold energy for the Fermi acceleration of electrons;

$``$ Luminous, short term $`\gamma `$–ray flares in FSRQs $``$ magnetic fields are much below equipartition with the relativistic plasma, unless $`\mathrm{\Gamma }1826`$;

$``$ No prominent spectral signatures of bulk Compton process in the soft/mid X–ray bands $``$ Collimation and/or acceleration of jets occurs beyond $`10^210^3`$ gravitational radii from the central source.

## Acknowledgments

This work was partially supported by Polish KBN grant 5 P03D 002 21 and NASA Chandra observing grant to Stanford University. M.S. thanks the Stanford Linear Accelerator Center for its hospitality during the completion of this work, and the Korean Institute for Advanced Study for hosting this interesting meeting.

## References


text:
# Toward a structural understanding of turbulent drag reduction: nonlinear coherent states in viscoelastic shear flows

## Abstract

Nontrivial steady flows have recently been found that capture the main structures of the turbulent buffer layer. We study the effects of polymer addition on these “exact coherent states” (ECS) in plane Couette flow. Despite the simplicity of the ECS flows, these effects closely mirror those observed experimentally: structures shift to larger length scales, wall-normal fluctuations are suppressed while streamwise ones are enhanced, and drag is reduced. The mechanism underlying these effects is elucidated. These results suggest that the ECS are closely related to buffer layer turbulence.

Rheological drag reduction, the suppression by additives of skin friction in turbulent flow, has received much attention since its discovery in 1947 Virk (1975); Lumley (1969); McComb (1990). For many polymer-solvent systems, the pressure drop measured in the pipe flow of the solution can be $`3050\%`$ less than for the solvent alone. The central rheological feature of drag-reducing additives is their extensional behavior in solution: for dilute polymer solutions in particular the stresses arising in extensional flow can be orders of magnitude larger than those developed in a shear flow. This fact is well-recognized; nevertheless the mechanism of interaction between polymer stretching and turbulent structure is not well-understood and the goal of the present work is to attempt to shed light on this interaction.

A key structural observation from experiments and direct numerical simulations (DNS) of drag-reducing solutions is the modification of the buffer region near the wall Walker and Tiederman (1990); Donohue et al. (1972); den Toonder et al. (1997); Draad et al. (1998); Escudier et al. (1999); Sureshkumar et al. (1997); Dimitropoulos et al. (2001). It has long been known that the flow in this region is very structured, containing streamwise vortices that lead to streaks in the streamwise velocity Robinson (1991); these structures are thickened in both the wall-normal and spanwise directions during flow of drag reducing solutions Walker and Tiederman (1990); Donohue et al. (1972). Because of its importance in the production and dissipation of turbulent energy Robinson (1991), any effort to mechanistically understand rheological drag reduction should address this region.

To better understand the effect of the polymer on the buffer layer, we wish to study a model flow that has structures similar to those seen in this region but without the full complexities of time-dependent turbulent flows. Fortunately, a family of such flows exists, in the recently-discovered “exact coherent states” (ECS) found by computational bifurcation analysis in plane Couette and plane Poiseuille flows Nagata (1990); Clever and Busse (1997); Waleffe (1998); Eckhardt et al. (1998); Waleffe (2001). These are three-dimensional, traveling wave flows (hence steady in a traveling reference frame) that appear via saddle-node bifurcations 111A saddle-node bifurcation, also known as a turning point, limit point or fold, arises when two steady states appear spontaneously as a parameter increases Strogatz (1994).  at a Reynolds number somewhat below the transition value seen in experiments Daviaud et al. (1992); Bottin et al. (1998). The structure of the ECS captures the counter-rotating staggered streamwise vortices that dominate the structure in the buffer region. From the dynamical point of view, there is evidence that these states form a part of the dynamical skeleton of the turbulent flow: i.e., they are saddle points that underlie the strange attractor of turbulence Jimenez and Simens (2001); Kawahara and Kida (2001). Finally, the nonlinear self-sustaining mechanism underlying these states has been elucidated Waleffe (1997). A perturbation of the base flow in the form of streamwise vortices redistributes the streamwise momentum of the flow. This redistribution creates spanwise fluctuations in the streamwise velocity, the “streaks”. The spanwise inflections in the streamwise velocity profile lead to a three-dimensional instability that develops into staggered nearly-streamwise vortices that regenerate the streaks. Because the ECS capture the structures of the buffer region and are mechanistically well-understood, we believe that they provide an excellent simplified, yet still exact, model flow for studying polymer drag reduction. The leading order effect of viscoelasticity on the ECS is therefore the focus of the present study.

To begin, we briefly describe a general result relating polymer stretch to flow kinematics. For a trajectory in a flow field, the Liapunov exponents give the Lagrangian time-averaged rate of stretch of material lines. If the largest Liapunov exponent, $`\sigma _{max}`$, is positive, the flow is extensional on average. In particular, for homogeneous turbulence, the expected value of $`\sigma _{max}`$ is positive Girimaji and Pope (1990) and we show below that this is also the case for the ECS. Now consider the dynamics of a Hookean dumbbell model of a polymer in a flow field. The end-to-end vector $`𝐪`$ of the dumbbell evolves in the flow field, $`𝐯`$, as:

$$\frac{D𝐪}{Dt}=𝐪\mathbf{}𝐯\frac{1}{2\lambda }𝐪+𝝃(t),$$
(1)
where $`D/Dt`$ is the time derivative evaluated on a fluid element, $`\lambda `$ is the stress relaxation time for the dumbbell and $`𝝃(t)`$ is the random Brownian force. Noting that an infinitesimal material line satisfies the same expression but with $`\lambda ^1`$ and $`𝝃`$ set to zero, it is straightforward to show that Hookean dumbbells will stretch indefinitely in a flow if and only if $`\lambda \sigma _{max}\text{We}_\sigma >\frac{1}{2},`$ where $`\text{We}_\sigma `$ is a Weissenberg number based on $`\sigma _{max}`$. This result is a specific statement of an idea that originated with Lumley Lumley (1972) (see also Ryskin (1987); Leal (1990); Chertkov (2000); Balkovsky et al. (2000)). The computations below show the importance of $`\text{We}_\sigma `$ in determining the effect of polymers on coherent structures.

We study here the effect of polymer on the exact coherent states that arise in a variant of plane Couette flow Waleffe (1998). Denoting the streamwise direction as $`x`$, the wall-normal direction as $`y`$, and the spanwise, or vorticity, direction as $`z`$, we consider a flow with boundary conditions $`\frac{v_x}{y}=1,v_y=\frac{v_z}{y}=0\text{ at }y=\pm 1.`$ The characteristic velocity, $`U`$, and the half-height of the channel, $`l`$, have been used to scale the velocity and positions, respectively. These “constant vorticity” boundary conditions provide an advantage over no-slip conditions in that they allow us to model only the buffer region in our domain by eliminating the viscous sublayer. (Exact coherent states found using no-slip BCs Waleffe (2001) show a qualitatively identical vortical structure, only offset from the wall by a small region comprising the viscous sublayer.) Periodic boundary conditions are applied in the streamwise and spanwise directions. For this study, the wavelength of the structures in the streamwise and spanwise directions are fixed at $`\mathrm{}_x=2\pi /0.40`$ and $`\mathrm{}_z=2\pi /1.0`$, respectively (or 165 and 66, if expressed in wall units at a Reynolds number of 110). For this flow, a trivial (Couette) base state exists, $`v_x(y)=y`$; the maximum mean velocity for the ECS is significantly reduced compared to the base state velocity due to the enhanced transport of momentum Waleffe (1998).

In our formulation, time, $`t`$, is scaled with $`l/U`$, and pressure, $`p`$, with $`\rho U^2`$, where $`\rho `$ is the fluid density. The stress due to the polymer, $`𝝉_p`$, is nondimensionalized with the polymer elastic modulus, $`G=\eta _p/\lambda `$, where $`\eta _p`$ is the polymer contribution to the zero-shear viscosity and $`\lambda `$ is the relaxation time for the polymer. The momentum and mass balances are

$$\frac{D𝐯}{Dt}=\mathbf{}p+\beta \frac{1}{Re}^2𝐯+(1\beta )\frac{1}{Re^2}\frac{1}{El}(\mathbf{}𝝉_p),$$
(2)
$$\mathbf{}𝐯=0,$$
(3)
where $`\eta _s`$ is the solvent viscosity, $`El=\lambda (\eta _s+\eta _p)/\rho l^2`$ and $`\beta =\eta _s/(\eta _s+\eta _p)`$. The Reynolds number, $`Re`$, is based on the total viscosity, $`Re=\rho Ul/(\eta _s+\eta _p)`$.

We calculate the polymer stress with the commonly used FENE-P model Bird et al. (1987b), which idealizes the polymer molecules as bead-spring dumbbells with finitely extensible springs. With this model, the non-dimensional structure tensor $`𝜶`$ ($`=\mathrm{𝐪𝐪}`$, where $``$ denotes ensemble average) evolves according to:

$$\frac{𝜶}{1\frac{tr𝜶}{b}}+\text{We}\left(\frac{D𝜶}{Dt}𝜶\mathbf{}𝐯\mathbf{}𝐯^T𝜶\right)=\frac{b𝜹}{b+2},$$
(4)
$$𝝉_p=\frac{b+5}{b}\left(\frac{𝜶}{1\frac{tr𝜶}{b}}\left(1\frac{2}{b+2}\right)𝜹\right),$$
(5)
where $`\text{We}=\frac{\lambda U}{l}`$ is the Weissenberg number based on the wall shear rate and $`b`$ is proportional to the maximum extension of the dumbbell: $`tr𝜶`$ cannot exceed $`b`$. A simple measure of the importance of extensional polymer stress is the magnitude of the parameter $`Ex=\frac{2}{3}\frac{b\eta _p}{\eta _s}.`$ In uniaxial extension with extension rate $`\dot{\epsilon }`$, $`Ex=1`$ implies that $`𝝉_p=𝝉_v\text{as}\dot{\epsilon }\mathrm{}`$ where $`𝝉_v`$ is the solvent contribution to the stress. The polymer can significantly affect the flow only when $`Ex1`$. In the present flow this parameter also represents the maximum ratio of polymer stress to Reynolds shear stress (the flux of $`x`$-momentum due to fluctuations in $`y`$-velocity); the scaling theory of Tabor and de Gennes Tabor and de Gennes (1986); de Gennes (1990); Sreenivasan and White (2000) treats the regime $`Ex1`$.

The governing equations are solved through a Picard iteration. A given velocity field is first used to calculate the polymer stress tensor, $`𝝉_p`$, by time-integrating Eq. 4 until a steady state is attained. For the new $`𝝉_p`$, a steady state of the momentum and continuity equations is found by Newton iteration. The resulting velocity field, is used to compute the new $`𝝉_p`$, and the process is repeated until the velocity field converges. Equations 2-3 are discretized as in Waleffe (1998), using a Fourier-Galerkin formulation with typically a $`7\times 19\times 7`$ grid. Equation 4 is discretized with a Fourier-pseudospectral method, typically with a $`32\times 32\times 32`$ grid, and time-integration performed with an Adams-Bashforth method. To achieve numerical stability, a small diffusive term is added to Eq. 4 (cf. Sureshkumar et al. (1997)) and integrated with a Crank-Nicholson scheme.

Before presenting the effects of the polymer on the ECS, we recall the result that $`We_\sigma >1/2`$ implies large stretch of polymer chains. For the Newtonian ECS at $`Re=110`$ on the lower branch of the bifurcation diagram (see Fig. 1), the velocity field is very nearly ergodic, with $`\sigma _{max}0.030`$. The condition $`\text{We}_\sigma >1/2`$ thus translates into $`\text{We}17`$ for large polymer stretch, and for $`Ex=O(1)`$, will define the onset condition for the polymer to begin to strongly affect the flow field. In DNS of a FENE-P fluid in plane channel flow, Sureshkumar, et al. Sureshkumar et al. (1997), found no drag reduction at $`\text{We}=12.5`$ and significant drag reduction at $`\text{We}=25`$; another recent DNS study Sibilla and Baron (2002) places the onset value at $`\text{We}20`$. This close correspondence between the onset condition predicted from the ECS kinematics and that found by DNS strongly suggests that the ECS model captures the essential structure of the buffer layer.

Fig. 1 shows how the addition of polymer stress affects the bifurcation diagram (i. e. the locus of steady state flows) for $`El=0.20`$ ($`We_\sigma 2/3`$) and $`Ex`$ = 1 to 3. On the $`y`$-axis of the diagram is the maximum value of the streamwise- and spanwise-averaged streamwise velocity $`v_x`$. (The trivial Couette solution is at $`\mathrm{max}v_x=1`$.) When $`Re`$ attains a certain value that we denote $`Re_{sn}`$, two new steady solutions appear in a saddle-node bifurcation. For small values of $`Ex`$, $`Re_{sn}`$ decreases compared to the Newtonian value, but once the polymer stress begins to exceed the viscous ($`Ex1.5`$), $`Re_{sn}`$ increases above the Newtonian value – the presence of the polymer is suppressing the ECS. Qualitatively identical behavior is observed experimentally in the onset Reynolds number for turbulent pipe flow Draad et al. (1998). Note that, for a given velocity $`U`$ and total viscosity $`\eta `$, the increase in critical Reynolds number corresponds to an increase in the characteristic length scale of the coherent structure, again consistent with experimental observationsWalker and Tiederman (1990); Donohue et al. (1972); Sureshkumar et al. (1997). Finally, we see that $`v_x`$ becomes larger for the viscoelastic flows than for the Newtonian – drag reduction occurs.

To examine more closely the effect of the polymer stress on the velocity fields, Fig. 2 shows results at constant $`Re`$ while varying $`El`$, or, equivalently, $`\text{We}_\sigma `$ (based on $`\sigma _{max}=0.030`$). Here we plot $`\mathrm{max}v_y^2`$, where the prime denotes the fluctuating part of a quantity, here wall-normal velocity. At $`Ex=1`$, after an initial increase, $`\mathrm{max}v_y^2`$ decreases below the Newtonian value and eventually saturates, as the polymer stress asymptotes at high $`\text{We}_\sigma `$ to a fixed value relative to the viscous stress. In this case, the polymer stretch becomes nearly uniformly large throughout the domain. The decrease in wall-normal velocity with $`\text{We}_\sigma `$ is even more drastic as the extensibility parameter $`Ex`$ increases. Similar trends are seen in the streamwise enstrophy and Reynolds shear stress. Since $`Ex`$ is related to the extensional viscosity of the viscoelastic solution, these results show the importance of extensional stresses in affecting the ECS. The spatial maximum of $`tr𝜶`$, which is proportional to the square of the polymer extension, is also presented in Fig. 2. Finally, note that the majority of both the polymer stretch and the change in wall-normal velocity occurs in the range $`0.1<\text{We}_\sigma <1.0`$. In contrast to the decreases in wall normal fluctuations, streamwise enstrophy and Reynolds shear stress for $`\text{We}_\sigma 0.1`$, the streamwise fluctuations $`v_x^2`$ are found to increase. All of these trends are observed in DNS and experiments Virk (1975); McComb (1990); Sureshkumar et al. (1997).

Turning now to how, and why, the flow structure is changed by the polymer dynamics, figures 3a-b show (a) the streamwise velocity $`v_x`$ at $`y=0.875`$ (where the maximum in the polymer stress occurs) for the Newtonian lower branch solution at $`Re=110`$ and (b) the difference $`v_{x,VE}v_{x,N}`$ between the viscoelastic (VE) and Newtonian (N) solutions. Here we see the “streak” (white ribbon) and – by adding the pattern of fig. 3b to that of 3a – that this streak is “straightened out” by the viscoelasticity.

Figure 3c shows a contour plot of $`f_x`$, the $`x`$-component of the force $`𝒇`$ exerted by the polymer on the fluid, ($`𝒇=(1\beta )Re^2El^1\mathbf{}𝝉_p`$), at $`y=0.875`$. This force is significantly negative and corresponds spatially to the region where fluid elements are leaving the streak to move into the vorticity-dominated regions. The $`y`$ and $`z`$ components of the force have also been examined; they are smaller than the $`x`$-component but are clearly seen to work directly against the vortex motions: e.g. where $`v_y`$ is highly positive in an upwelling, $`f_y`$ is highly negative. This behavior is also seen in the buffer layer structures in the DNS study of a drag-reducing polymer solution by De Angelis *et al.* de Angelis et al. (2002). The origin for this structure of the polymer force field becomes apparent on examination of polymer stresses along fluid trajectories: polymer molecules stretch in (or moving into) the streak regions, remaining highly stretched until they begin to leave the streak. As molecules move from the streak into and around the vortices, they relax. The spatial gradients in stress accompanying this relaxation work against the vortices, “unwinding” them. This vortex suppression leads to collapse of the self-sustained process – or more precisely to a shift of the process to larger scales – and thus to drag reduction.

To summarize, we list several points of agreement between our results and observations from DNS and experiments in fully turbulent flow, namely: (1) the ECS bear a strong similarity to the structures observed or educed from structural studies of the buffer layer and apparently underlie its dynamics; (2) for $`O(1)`$ values of $`Ex`$, the onset Weissenberg number for drag reduction predicted from the ECS kinematics agrees closely with DNS results; (3) the effects of viscoelasticity on the Couette flow ECS are very similar to those observed in the buffer layer: (a) wall normal velocity fluctuations are suppressed and streamwise ones enhanced, (b) Reynolds shear stress decreases, (c) streamwise vorticity decreases, (d) the velocity fluctuations and polymer force are anticorrelated, and (e) drag is reduced. Finally, at fixed $`U`$ and $`\eta `$ the upward shift in the onset Reynolds number corresponds to an increase in length scale for the structures, again consistent with experiments. These successes show that studying the ECS holds promise for capturing the essential physics of drag reduction. Indirectly, they also validate the view that the ECS underlie Newtonian turbulence, because the effects of polymers on the ECS so closely mirror their effects on full turbulence.

The authors gratefully acknowledge support from NSF and the donors of the Petroleum Research Fund, administered by the American Chemical Society.


text:
# Broken Ergodicity in a Stochastic Model with Condensation

\[

## Abstract

We introduce a variant of the asymmetric random average process with continuous state variables where the maximal transport is restricted by a cutoff. For periodic boundary conditions, we show the existence of a phase transition between a pure high flow phase and a mixed phase, whereby the latter consists of a homogeneous high flow and a condensed low flow substate without translation invariance. The finite system alternates between these substates which both have diverging lifetimes in the thermodynamic limit, so ergodicity is broken in the infinite system. However, the scaling behaviour of the lifetimes in dependence of the system size is different due to different underlying flipping mechanisms.

\]

The asymmetric random average process (ARAP)  represents a simple fundamental model for numerous applications, especially in the field of interdisciplinary questions. It has a close relationship to the q model describing force distributions in granular materials , aggregation and fragmentation processes (e.g. mass transport problems)  and traffic flow theory, in particular a continuous Nagel-Schreckenberg model .

The truncated ARAP (TARAP) studied in this letter, although equipped with rather simple homogeneous short-range interaction and periodic boundary conditions, features the surprising behaviour of broken ergodicity, i.e. a state space decomposition into dynamically unconnected subsets. In addition, these substates are totally different: the system may reside either in a homogeneous high flow phase or in a condensed low flow phase where a finite fraction of the total mass is located on a randomly chosen site. While examples of infinite aggregation  and ergodicity breaking  are known for 1D nonequilibrium systems, the simultaneous occurrence of these phenomena seems to be undiscussed and relevant for a lot of applications.

The ARAP is defined on a 1D periodic lattice with $`L`$ sites. Each site $`i`$ carries a non-negative continuous mass variable $`m_i_0^+`$. In every discrete time step $`tt+1`$ for each site a random number $`r_i[0,1]`$ is generated from a time-independent probability density function (pdf) $`\varphi `$, sometimes called fraction density, that may depend on the actual configuration $`m=(m_1,\mathrm{},m_L)`$. So the universal form can be written as $`\varphi =\varphi (r_1,\mathrm{},m_1,\mathrm{})`$. The fraction $`r_i`$ determines the amount of mass $`r_im_i`$ transported from site $`i`$ to site $`i+1`$. The transport is completely asymmetric, i.e. no mass moves in the opposite direction $`i+1i`$, and we get

$`m_i(1r_i)m_i+r_{i1}m_{i1}.`$ (1)
These update rules correspond to a parallel dynamics. Due to the conservation of the total mass $`M=_im_i`$ the density $`\rho =\frac{M}{L}`$ is fixed. We would like to add that this so called stick representation is equivalent to a particle picture  not used in this letter.

The simplest version of the ARAP is obtained if we use the state independent uniform distribution defined by the fraction density $`\varphi =1`$. We will refer to this system as the free ARAP. Some results have been derived for this model so far whereby the exactness of product measure for $`L\mathrm{}`$ and its form

$$P(m)=\underset{i}{}P(m_i)\text{with}P(m_i)=\frac{4m_i}{\rho ^2}e^{2m_i/\rho }$$
(2)
is the most relevant for our work. Eq. (2) has been obtained in the context of the q model  that shares many properties with the ARAP. Further results for the ARAP can be found in .

In the ARAP the mass $`r_im_i`$ transferred from site $`i`$ to site $`i+1`$ is in principle unbounded. This is different in the truncated ARAP that we introduce here. In the TARAP all transfers of masses $`r_im_i`$ larger than a cutoff $`\mathrm{\Delta }>0`$ are rejected. We like to focus on the simplest model only, the truncated free ARAP. The corresponding fraction density is then given by $`\varphi =_i\varphi (r_i,m_i)`$ with

$$\varphi (r_i,m_i)=\left[1R(m_i)\right]\delta (r_i)+\mathrm{\Theta }\left(R(m_i)r_i\right),$$
(3)
where $`\mathrm{\Theta }`$ is the Heaviside step-function and

$$R(m_i)\mathrm{min}(1,\frac{\mathrm{\Delta }}{m_i})$$
(4)
represents the maximum possible fraction. Note that $`\varphi `$ has become locally state dependent, i.e. the pdf depends on the mass $`m_i`$ at the corresponding site explicitly.

Without loss of generality we set $`\rho =1`$ for the rest of the letter because every TARAP defined by $`(\rho ,\mathrm{\Delta })`$ can be mapped onto a $`(1,\frac{\mathrm{\Delta }}{\rho })`$-system (Fig. 1). Furthermore we introduce the rescaled cutoff

$$\stackrel{~}{\mathrm{\Delta }}=2L^{\frac{1}{2}}\mathrm{\Delta }$$
(5)
which ensures $`L`$-independence of the critical point.

We begin by investigating the relation between the steady state current $`J`$, defined by the average mass transfer per site, and the cutoff parameter $`\mathrm{\Delta }`$ (Fig. 1). For $`\mathrm{\Delta }0`$ the flow vanishes because the transferred mass per site is always smaller than $`\mathrm{\Delta }`$ while for $`\mathrm{\Delta }\mathrm{}`$ the system behaves like a free ARAP and the value of $`J`$ tends to its maximum $`J^{\text{max}}\frac{1}{2}`$.

If we study the process in more detail we discover the system to exist in two phases whereby the transition point $`\stackrel{~}{\mathrm{\Delta }}_c=1`$ has been determined numerically and by analytical approximations.

For rescaled cutoffs $`\stackrel{~}{\mathrm{\Delta }}>\stackrel{~}{\mathrm{\Delta }}_c`$ the steady state mass distribution is nearly identical to the one of the free ARAP and approximately given by (2). The flow is independent of $`\stackrel{~}{\mathrm{\Delta }}`$ and corresponds to the maximum current $`J_{\text{high}}^{\text{max}}J^{\text{max}}`$. Therefore we refer to this parameter range as the high flow phase.

For rescaled cutoffs smaller than $`\stackrel{~}{\mathrm{\Delta }}_c`$ the steady state of the finite system is a composition of two different substates: the system can either exist in a high flow state (with properties as described above) or a state given by a macroscopic condensate, i.e. a finite fraction of the total mass $`M`$ resides on one randomly chosen site (even in the thermodynamic limit). So the mass aggregation is proportional to $`L`$, but not extensive in space. The remaining mass is distributed equally with an algebraically decaying mass distribution. The current of the condensate state depends on $`\stackrel{~}{\mathrm{\Delta }}`$ but does not exceed $`J_{\text{low}}^{\text{max}}\frac{1}{4}`$. Therefore we call this state low flow state and denote the parameter region $`\stackrel{~}{\mathrm{\Delta }}<\stackrel{~}{\mathrm{\Delta }}_c`$ as mixed phase.

In a finite system the transition probabilities between low and high flow states are small but nonzero. The system switches between these states while evolving in time and an alternating current-time relation is obtained (Fig. 2). Note that the switching time between the two states is much smaller than their lifetimes.

This alternating behaviour is very similar to systems with spontaneously broken symmetry . However, in general the flipping is between states of broken symmetry. For example in  an ASEP with two particles ($`+`$) and ($``$) is introduced and in regimes of spontaneously broken symmetry the system switches between states dominated by ($`+`$) or ($``$) particles. In case of the TARAP the flipping occurs between the symmetric high flow state which is translation invariant and a low flow state in which this symmetry is broken.

In the thermodynamic limit the average lifetimes $`\tau _H`$ and $`\tau _L`$ of the high and low flow states diverge (Figs. 3 and 4). This implies that the steady state in the mixed phase is not unique. Ergodicity is broken and the steady state can either be in the high flow state or the low flow state depending on the initial condition. However, the lifetimes of the substates do not scale equally with size $`L`$ (see next paragraph).

Although Monte Carlo simulations are difficult in the mixed phase since the lifetimes of the substates are very large it is safe to say that both the average lifetimes of low and high flow states diverge in the thermodynamic limit. Fig. 3 indicates that $`\tau _H`$ increases exponentially for all $`\stackrel{~}{\mathrm{\Delta }}<\stackrel{~}{\mathrm{\Delta }}_c`$. Assuming $`\tau _H\mathrm{exp}(\alpha (\stackrel{~}{\mathrm{\Delta }})L)`$ we obtain that $`\alpha `$ is an increasing function of $`\stackrel{~}{\mathrm{\Delta }}`$. Note that in the pure high flow phase $`\stackrel{~}{\mathrm{\Delta }}>\stackrel{~}{\mathrm{\Delta }}_c`$ the relation $`\tau _H=\mathrm{}`$ holds for any $`L`$. The corresponding data for $`\tau _L`$ in the low flow state are better fitted algebraically than exponentially (Fig. 4). We see again that the exponents obtained by the assumption $`\tau _LL^{\beta (\stackrel{~}{\mathrm{\Delta }})}`$ increase while moving into the mixed phase, i.e. for $`\stackrel{~}{\mathrm{\Delta }}`$ tending to zero.

Although both transitions are driven by fluctuations, the above scaling behaviours reflect different switching mechanisms. While the high $``$ low transition is based on a collective effect which involves all sites of the lattice, the opposite transition involves only the lattice site where the macroscopic condensate is located. Under these assumptions we have analytically derived approximations that reflect qualitatively the lifetime behaviour .

Note that in spontaneous symmetry breaking both transitions obey the same purely fluctuation driven mechanism yielding the same exponential behaviour of the lifetimes. A flipping mechanism between two alternating symmetric substates that is not purely driven by fluctuations is studied in , resulting in a logarithmic lifetime dependence.

The above results are mainly based on Monte Carlo simulations. We continue now by presenting a few analytical approximations which support the numerical findings.

Starting point is a closer look at the average mass transfer per site (local current)

$$J(m)rm_{\varphi (r,m)}=\{\begin{array}{cc}\frac{1}{2}m\hfill & \mathrm{for}0m<\mathrm{\Delta }\hfill \\ \frac{\mathrm{\Delta }^2}{2m}\hfill & \mathrm{for}\mathrm{\Delta }m<\mathrm{}\hfill \end{array}.$$
(6)
This shows that the average mass shift $`J(m)`$ tends to zero for $`m\mathrm{}`$ and $`m0`$. So high (low) columns shrink (grow) very slowly and accordingly low flow states, resp. one-stick configurations, are very stable. On the other hand, homogeneous configurations maximize the current. We will exemplify this in the following paragraphs.

First we study the low flow state using the approximation

$$J(m)_{P_i}=J\left(m_{P_i}\right).$$
(7)
Here $`P_i`$ denotes the steady state single-site distribution of site $`i`$. Equation (7) holds for stable distributions $`P_i(m)\delta (mm_i)`$ which we assume here. Introducing the nomenclature $`m_im_{P_i}`$ we obtain from the continuity equation in the steady state

$$J_{\text{low}}J(m_i)=J(m_{i+1})\text{for all }i.$$
(8)
Restricting now to the case where only one high column exists at site $`j`$, i.e. $`m_+m_j>\mathrm{\Delta }`$, the remaining columns are all of the same mass (see (8)), i.e. $`m_{}m_i<\mathrm{\Delta }`$ with $`ij`$. So, using (6) and (8), the quantities $`m_\pm `$ are related by $`m_+m_{}=\mathrm{\Delta }^2`$ and the mass conservation law $`\rho L=(L1)m_{}+m_+`$.

Equipped with these relations we are able to compute $`m_\pm `$ and finally the current in the low flow state $`J_{\text{low}}`$. Assuming $`L1`$ our calculations lead to the formula

$$J_{\text{low}}=J_{\text{low}}^{\text{max}}\left\{1\sqrt{1\stackrel{~}{\mathrm{\Delta }}^2}\right\}.$$
(9)
So the current of the low flow state increases from zero to $`J_{\text{low}}^{\text{max}}`$ for $`0<\stackrel{~}{\mathrm{\Delta }}\stackrel{~}{\mathrm{\Delta }}_c`$. For rescaled cutoffs larger than $`\stackrel{~}{\mathrm{\Delta }}_c`$ relation (9) is not defined and the one-stick model is not valid anymore. Thus we have derived an upper bound for the occurrence of the low flow state in accordance with the simulations. As Fig. 5 shows the analytical result agrees very well with Monte Carlo data.

By extending the calculations above to $`N`$ mass aggregates  the validity range in terms of $`\stackrel{~}{\mathrm{\Delta }}`$ shrinks and we obtain smaller critical cutoffs $`\stackrel{~}{\mathrm{\Delta }}_c`$. Furthermore the flux increases. This indicates that the mixed phase is described best by the one-condensate picture, especially in the vicinity of the critical point. This is confirmed by numerical investigations that have shown that configurations with two or more aggregates are not stable. The pre-condensates struggle for masses until only one stick is left. This mechanism of *competition* is different to *coarsening* where aggregates emerge at several sites and then coalesce .

Focussing on the high flow state, resp. phase, every site carries the same mass $`\rho `$ on average. For large $`L`$ we approximate, motivated by the numerical results, the mass distribution by the product measure solution of the free ARAP (2) and calculate by the help of (6) the high flow state current

$$J_{\text{high}}=J_{\text{high}}^{\text{max}}\left\{1(1+L^{\frac{1}{2}}\stackrel{~}{\mathrm{\Delta }})e^{L^{\frac{1}{2}}\stackrel{~}{\mathrm{\Delta }}}\right\}\underset{L1}{}J_{\text{high}}^{\text{max}}$$
(10)
for all rescaled cutoffs. So the high flow state exists both in the mixed phase and in the high flow phase. Equation (10) also predicts that the current in the high flow state differs from $`J_{\text{high}}^{\text{max}}`$ only for small absolute cutoffs $`\mathrm{\Delta }`$ in finite systems. We have verified this by Monte Carlo simulations although it is no longer possible to distinguish the states by their flows which are nearly identical in this parameter range. So we used the appearance of a macroscopic condensate as a criterion.

A model that is strongly related to the TARAP is studied in  where a similar kind of phase separation as in the mixed domain has been observed. The underlying process is the continuum limit (Krauss model) of the Nagel-Schreckenberg cellular automaton model  that is defined by continuous velocities and spatial coordinates. Although this traffic model is given by more complex dynamics than the TARAP, both processes have a common feature: moves may be rejected if a uniformly distributed random variable exceeds a given threshold. So the TARAP can be viewed as a toy model that catches some of the fundamental physics behind the Krauss model. We have also rewritten a simplified version of the Krauss model in terms of fraction densities to point out the similarities with the TARAP explicitly .

In the Krauss model an additional congested phase has been observed . An analogous phase, corresponding to a pure condensed phase, could be expected for the TARAP in the limit of small cutoffs. However, in that regime the lifetime of the low flow state is much larger than the corresponding lifetime of the high flow state, in particular for small systems. Therefore, in Monte Carlo simulations it is difficult to distinguish between a pure condensed phase and the mixed regime. Furthermore, our analytical calculations have not shown any evidence for a transition point separating a (new) pure low flow phase and the mixed phase.

What are the essential ingredients of the TARAP leading to the observed phenomena? First *truncation* and an *unbounded local state space* allow for an unlimited condensation. The formation of infinite aggregates can also be found in a model of aggregation and fragmentation  or a special zero-range process , both equipped with unbounded but discrete state variables. However, in these examples the condensed phase is unique and not associated with a high flow counterpart, i.e. high flow phase and congested phase are separated by a phase transition. The same is true for a discretized version of the TARAP that we have investigated .

Therefore we believe that in case of the TARAP the *continuous state variables* allow for the coexistence of high and low flow state in the mixed phase, yielding nonsymmetric ergodicity breaking in the *thermodynamic limit*. The TARAP is free of an intrinsic mass scale in contrast to models defined on integer state space, equipped with a smallest mass unit $`1`$. Thus, in the TARAP we obtain a trivial dependence on the density $`\rho `$ which is also reflected in the scaling law $`(\rho ,\mathrm{\Delta })(1,\mathrm{\Delta }/\rho )`$. However, in discrete systems physics may change in the low density regime, e.g. an ARAP with total mass $`M=1`$ is nothing else than a one-particle ASEP. This difference may indicate why in the continuous variant the homogeneous high flow state can persist over the whole parameter range, resulting in the coexistence phenomenon.

We would like to thank Joachim Krug, Kavita Jain and Martin Evans for interesting und helpful discussions.