text:
# Possible contribution of olefins and heteroatoms to the Unidentified Infrared Bands

## 1 Introduction

A large number of high-quality mid-IR spectra are now available in emission from HII regions (e.g. Crété et al. cre99 (1999)), PDRs (PhotoDissociation Regions; see Bregman et al. breg94 (1994)), PNe (Planetary Nebulae; e.g. Allamandola et al. all89 (1989), Molster et al. mol96 (1996), Witteborn et al. wit89 (1989)), transition objects (e.g. Buss et al. bus93 (1993), Geballe et al. geb92 (1992), Justtanont et al. jus96 (1996), Kwok et al. kwo01 (2001)), M supergiants (Sylvester et al. syl94 (1994)), RNe (Reflection Nebulae; see Moutou et al. mou99 (1999), Sloan et al. slo99 (1999), Uchida et al. uch00 (2000)) and, in absorption, from protostars (e.g. Bregman et al. breg00 (2000)) and the GC (Galactic Center; see Chiar et al. chi02 (2002), Tielens et al. tiel96 (1996)). Several authors have assembled spectra of different origins (e.g. Boulanger and Cox boul98 (1998), Peeters et al. pee02 (2002)) and given “generic” (typical) spectra (e.g. Hony hon01 (2001)) of the UIBs (Unidentified IR Bands).

In the vast majority of these spectra, the “11.3” feature peaks near 890 cm-1 (11.25 $`\mu `$m) and emerges conspicuously from a plateau which extends from $``$11 to $``$ 13.5 $`\mu `$m. A second, clearly separate and typically weaker, feature also sits upon this plateau and peaks at 12.7 $`\mu `$m (790 cm-1). When it is possible to delineate an underlying continuum for the plateau, the intensity of the latter is of the same order as that of the 12.7 feature, typically much weaker than that of the 11.3 feature.

While the 12.7 feature remains unexplained, the 11.3 feature has been traditionally assigned to oop (out-of-plane) bending of C-H bonds at the periphery of aromatic clusters and in the absence of other, immediately adjacent, such bonds (solo bending). In laboratory spectra of PAHs (Polycyclic Aromatic Hydrocarbons) of different sizes and shapes, and various substitutions for peripheral hydrogen atoms, this type of vibration is found to occur indeed in the range 11-12 $`\mu `$m but, to my knowledge, no single known PAH displays a solitary, prominent, line at the peak position of the celestial feature, which is quite stable and well defined. When a feature is present in a PAH, in the 11-12 $`\mu `$m range, it is generally accompanied by a few other bands of comparable intensity, near 12 and 13 $`\mu `$m; these are attributed to duo, trio and quartet bendings, i.e. *coordinated* oop vibrations of, respectively, 2, 3 and 4 *adjacent* C-H bonds (see Colthup et al. col90 (1990), Smith smi79 (1979)). Present models are, therefore, unable to reproduce typical UIB spectra. This has recently prompted reassessments of previous assignments, leading to various new proposals: essentially ionized PAHs (cations) and large, compact, neutral clusters. The present research addresses the same issue and suggests new assignments for the 11.3 and 12.7 features, based on an unbiased perusal of IR data bases, extended to non-aromatics and to two heteroatoms among the most abundant in the universe: nitrogen and sulfur.

## 2 Discussion of the current aromatic assignments of the 11.3 feature

As implied by its designation, the solo C-H bending vibration is that of a C-H bond isolated from the next one by at least a bare C site. One conceivable way of obtaining this vibration without exciting duos, trios, etc. too is to take a linear chain of orthocondensed benzenic rings and bend it into a closed annulus: here, all the oop C-H bendings are indeed of the solo type. Unfortunately, such a structure is unknown in the laboratory and its solo vibration is not guarantied to fall at 11.25 $`\mu `$m.

Hony et al. (hon01 (2001)) proposed a solution consisting of a large diamond-shaped graphene (planar sheet made of compactly condensed hexagonal carbon rings). In this case, all C-H oop bends are again of the solo type, except at the four corners, where duos occur. Both objections raised above against the annular chain of benzene rings also apply here. Moreover, large graphene sheets are known to curl up and fold over into fullerenes or nanotubes (Robertson et al. rob92 (1992)), whose IR features were never observed in the sky. Finally, even if such large graphenes existed, they would necessarily display a strong continuum because of their large electronic conductivity, which would intolerably reduce the UIB feature contrast.

For any known isolated PAH, enough dehydrogenation may conceivably eliminate duos, trios, etc. But, for one thing, there is no reason for dehydrogenation to remain nicely distributed around the cluster, so as to avoid all adjacencies. More importantly, the relatively high density of atomic hydrogen in UIB sources in general, and the high measured rate of H2 formation on dust (R$``$3 10-17cm3sec-1 in diffuse clouds) make it highly improbable that many dangling C bonds will remain simultaneously unoccupied for any long period of time.

For want of an ideal, *single*, aromatic hydrocarbon which would mimic the single 11.3 feature, let us now explore the possibility that aromatic hydrocarbons *as a family* could mimic the whole, or part, of the 11-14 $`\mu `$m spectrum. For this purpose, Papoular et al. (rp89 (1989)) displayed the average of the spectra of the 18 smaller PAHs up to coronene. As expected, this composite spectrum exhibits a “forest” of features, between 11 and 14 $`\mu `$m, which blend into an asymmetrical band peaking at $``$13.7 $`\mu `$m, in sharp contrast with the UIB spectrum.

Later on, Schlemmer et al. (schl94 (1994)) studied the emission spectra of “hot” gaseous molecules of naphtalene and pyrene between 3 and 7.5 $`\mu `$m. They noted the presence of strong features between 5 and 6 $`\mu `$m (not observed in the sky) and a ratio of mid- to near-IR feature intensities much lower than in UIB spectra. They concluded that small, neutral, PAHs cannot be the carriers of the UIBs. Nevertheless, Allamandola and co-workers (see references in Allamandola et al. all99 (1999)) extended the study of neutrals to larger sizes, as far as dicoronylene C43H20, isolated in low-temperature rare gas matrices. Spectra of individual, and mixtures of, neutral PAHs are given in Allamandola et al. (all99 (1999)), Hudgins et al. (hud99 (1999)), Sloan et al. (slo99 (1999)) and Hony et al. (hon01 (2001)). These results confirm the expectation, in the spectra of neutral PAHs, of 4 dominant oop C-H bending bands, peaking around 11.3 (solo), 12 (duo), 13 (trio) and 13.5 $`\mu `$m (quartet). All intensities are of the same order, so no mixture could be proposed to mimic the prominent, isolated, 11.3 feature of generic UIB spectra.

The same can be said of the random mixtures of small and intermediate PAHs which are known to constitute the backbone of bulk kerogens (Papoular rp01 (2001)) and coals (Papoular et al. rp89 (1989)). *For these materials, the oop C-H bending bands form a plateau extending from $``$11 to $``$13.5 $`\mu `$m, with 3 or 4 low-contrast peaks corresponding to the various adjacencies*. Guillois et al. (gui96 (1996)) showed that, as a consequence, these materials correctly mimick the rippled 12-$`\mu `$m plateau of type-B sources (see Tokunaga tok97 (1997) and Kwok et al. kwo01 (2001)). This is not the case with type-A objects, where the plateau is also present, but dominated by a lone and narrower 11.25-$`\mu `$m feature in which we are primarily interested here.

Theorists in the early ’90s started *ab initio* spectral calculations of PAHs, including charged ones, in the hope of achieving a better fit to the observed relative feature intensities. Langhoff (lang96 (1996)) then produced high-quality theoretical spectra of PAH neutrals, cations and anions, using density functional theory. This demonstrated once again the inadequacy of neutrals as UIB carriers, but showed that cations do a far better job regarding relative feature intensities, although the frequency agreement with observations was worsened to some extent. The 11.3 feature was tentatively assigned to the duo C-H oop mode, and the 12.7 to the trio mode of *cations*.

These results fuelled the interest of experimentalists and, a few years later, a rich laboratory data base was produced (see Sloan et al. slo99 (1999) and Hony et al. hon01 (2001) for references and typical spectra). From these rather concurrent theoretical and experimental results, it appears that, in the spectral range of interest here,

–the solo oop C-H bending of neutrals fits the celestial 11.3 feature in frequency far better than does that of cations; but it is too strong;

–the solo oop C-H bending of cations falls near 11$`\mu `$m, and so may be the carrier of the weak feature often detected on the blue wing of the 11.3 feature; but if this were the case, and if cations were ubiquitous, the 11 $`\mu `$m feature should be much stronger than observed;

–both neutrals and cations also exhibit equally strong duos, trios, etc. as far as 13.5 $`\mu `$m; so neither can explain the dominance of the 11.3 feature;

–both in theory (see Langhoff lang96 (1996); Bakes et al. ba01 (2001)) and in experiment (see Hony et al. hon01 (2001)), the whole spectrum of PAHs depends heavily on their degree of ionization; one would therefore expect very different spectra from regions of different UV field intensities as measured by their G values. In fact, the observed UIB spectra vary only slightly from G=1 to G=104 (see Boulanger and Cox boul98 (1998); Sellgren sel01 (2001)); it seems therefore that, in the framework of the PAH model of the UIBs, this remakable stability can hardly be accounted for, unless one is prepared to invoke very stringent and finely tuned electrical, optical and thermal properties for these molecules, as well as very well defined mixtures thereof, depending on environmental conditions, the reality or realizability of which have yet to be proven (see Verstraete et al. vers01 (2001); Bakes et al. ba01 (2001)). We have therefore to conclude with Boulanger and Cox that our problem is not likely to be solved by ionized PAHs.

Even the loose assignment of the 11.3 feature to the general class of aromatic bonds is called into question by angular scans of the sky, indicating the non-coincidence of the intensity peak of this feature with that of the definitely aromatic C-H stretch feature at 3.28 $`\mu `$m (Bregman et al. breg94 (1994); Crété et al. cre99 (1999)).

## 3 Non-aromatic assignment of the 11.3-$`\mu `$m feature

The slight variations in band shape and peak wavenumber of the UIBs from different sources suggest that these bands are mostly carried by *functional groups* rather than particular molecules with definite spectra. These groups are made of small numbers of atoms, always bonded together in a configuration whose geometry changes only slightly according to the main structure to which it is attached. Each group is characterized by a small number of vibration modes. *Correlation charts have been drawn, which graphically display the vibration strength and the range covered by the changing peak wave number of each mode of every group as its attachment changes*. These charts are not limited to aromatic hydrocarbon groups, but include olefinic and aliphatic bondings as well as higher mass atoms. No group is expected to carry all the UIBs. We have surveyed such charts in Colthup et al. (col90 (1990)), Smith (smi79 (1979)) and others, in search of alternative assignments to the 11.3 feature, then to the 12.7 feature (see Sec. 4).

The best fit to the former turns out to be the (olefinic) *vinylidene* group R2C=CH2, where R represents an alkane substituent (saturated carbon; see fig. 1). This has a *strong* CH2 wag in the range 885-895 cm-1 (11.2-11.3 $`\mu `$m), in which both H atoms oscillate in phase and oop. The other IR-active modes are: a CH2 deformation (weak) near 1415 cm-1 (7.07 $`\mu `$m); a C=C stretch (medium strong) at 1639-1661 cm-1 (6.02-6.1 $`\mu `$m); a weak overtone of the CH2 wag at 1775-1792 cm-1 (5.58-5.63 $`\mu `$m); a weak symmetric (in phase) stretch near 2985 cm-1 (3.35 $`\mu `$m) and a medium strong asymmetric (out-of-phase) C-H stretch at 3070-3100 cm-1 (3.23-3.26 $`\mu `$m). An example of known molecule including this group is dimethylethylene, where R=CH3. The vinylidene group may also be attached to a 5- or 6-membered carbon ring without noticeable spectral changes, as in methylene cyclohexane (Colthup et al. col90 (1990)).

Arguments in favour of this assignment are: a) there are no misplaced strong features and no added continuum; b) two secondary features fall in the vicinity of weak UIBs observed by Peeters et al. (pee02 (2002)) at 6 and 7 $`\mu `$m; c) the shift of the asymmetric CH stretch relative to the 3.28-$`\mu `$m feature may help explaining the large variations observed in the ratio of the 3.3 to 11.3 UIBs (see Bregman et al. breg94 (1994)) and the lack of spatial coincidence of the peak intensities of these bands (Crété et al. cre99 (1999)); d) the non-aromaticity of the group may explain part of the large fluctuations of the intensity ratio of the 6.2 to 11.3 UIBs (see Hony et al. hon01 (2001)); e) this assignment is consisitent with the established presence in space of carbon double bonds as in cumulene carbenes (see Thaddeus et al. tha98 (1998)).

In the similar *vinyl* group, one of the R’s of vinylidene is replaced by an H (Colthup et al. col90 (1990)). This is interesting in that the strong CH2 wag is then shifted to 905-910 cm-1 ($``$11 $`\mu `$m), near the weak feature observed in the wing of the 11.25-$`\mu `$m UIB (see Sloan et al. slo99 (1999)). Vinyl also has a (weaker) *trans* C-H wag at 985-995 cm-1 ($``$10.1 $`\mu `$m), in which the 2 opposite C-H bonds oscillate oop. If vinyl groups are present in space, this feature should be detectable in the UIB spectrum whenever the 11.1 feature is observed, as is perhaps the case in Sloan et al. (slo99 (1999)).

## 4 The 12.7-$`\mu `$m feature ($``$787 cm-1)

In space, the 12.7 feature is not directly correlated with the 11.3 UIB (see Hony et al. hon01 (2001)) and, indeed the candidates considered above do not exhibit this feature. An independent search is therefore required for it. This feature falls in the gap between the duo and trio aromatic C-H oop peaks. Exceptionally, naphtalene does exhibit a strong feature at 780 cm-1 (see Schrader schr89 (1989)); numerical simulation of this molecule (Papoular rp01b (2001)) shows that this frequency corresponds to all 8 CH’s wagging in phase oop, the two opposite edge pairs being most active.

A somewhat similar geometry is provided by para-xylene, a benzene ring where CH3’s are substituted for 2 opposite H’s , which carries a very strong feature at 793 cm-1 (Schrader schr89 (1989)). Several other para-substituted benzenes exhibit a very strong feature in the range covered by the 12.7 UIB. 4-picoline C6H7N has a closely related symmetry (Schrader schr89 (1989)): one C is substituted with N and the opposite H is replaced by a methyl (CH3); the resonance at 799 cm-1 again corresponds to 2 opposite, in phase, duo C-H wags and is very strong. An interesting feature of these molecules is their strong vibrations near 21 $`\mu `$m. On the other hand, they also have undesirable, albeit weaker, bands near 7 $`\mu `$m.

It seems, therefore, preferable to turn to primary and secondary *organic nitrites* R-O-N=O, of which there are two isomers, according to the relative orientations of the R-O and N=O bonds (Lin-Vien et al. lin90 (1990)). The desired vibration is the *strong* N-O stretch at 775-814 cm-1 (12.3-12.9 $`\mu `$m; fig. 1). There is only one other main band, due to N=O stretching, at 1613-1681 cm-1 and likely peaking near 6$`\mu `$m. A subsidiary band occurs at 565-625 cm-1, due to group bending and is reminiscent of the weak 16.4 $`\mu `$m-band observed by Moutou et al. (mou00 (2000)).

In aliphatic *amines*, R-NH2, the NH2 wagging and twisting vibrations give rise to broad, strong, usually multiple bands around 800 cm-1 (Colthup et al. col90 (1990)), and so may also be considered as a candidate. A test of its presence could be the detection of its N-H stretch, which is solitary around 3400 cm-1 (fig. 1), but hard to observe because of telluric absorption.

## 5 Inclusion in condensed dust

As a by-product of the present correlation search, *sulfones*, SO2 turned out to have 2 strong stretch bands near the 7.7- and 8.6-$`\mu `$m UIBs, as well as weaker deformation around 17.4 and wag around 20 $`\mu `$m (fig. 1). Since O, N and S are among the most abundant elements, and since molecular combinations thereof have been detected in space, the classes of functional groups in fig.1 may all be considered as possible contributors to the UIBs. The shape and variability of the latter speak in favour of the groups being included randomly in condensed matter rather than molecules. A possible model of that matter is provided by the 3-D models of natural coals and kerogens, which invoke an amorphous hydrocarbon, only partly aromatic, interspersed with substitutions for hydrogen and carbon atoms with precisely the heteroatoms considered above: O, N, S (see Papoular et al. rp89 (1989) and Papoular rp01 (2001)), and “a complex collection of alkane and alkene side chains” (Kwok et al. kwo01 (2001)). In this framework, the coal/kerogen-like material, in its less evolved forms, provides models for the broad, mostly aliphatic, features of B-type objects.

On the other hand, the large range of its more mature (aromatic) forms, augmented with adequate concentrations of functional groups such as those proposed above, will help tailoring the great variety of relative intensities, bandwidths and contrast of aliphatic and aromatic features in celestial spectra (see Peeters et al. pee02 (2002)), with no need to invoke dehydrogenation, ionization or size effects. By analogy with laboratory heat treatment (Papoular et al. rp89 (1989)), this continuum of maturation states (tending towards type-A objects) may result, in space, from successive gas shocks or radiative heating or, still, photochemistry (Kwok et al. kwo01 (2001)).


format: text

text:
# Roughness of stylolites: a stress-induced instability with non local interactions

## Abstract

We study the roughness of stylolite surfaces (i.e. natural pressure-dissolution surfaces in sedimentary rocks) from profiler measurements at laboratory scales. The roughness is shown to be nicely described by a self-affine scaling invariance. At large scales, the roughness exponent is $`\zeta _10.5`$ and very different from that at small scales where $`\zeta _21.1`$. A cross-over length scale at around $`\lambda _c=1`$mm is well characterized and interpreted as a possible fossil stress measurement if related to the Asaro-Tiller-Grinfeld stress-induced instability. Measurements are consistent with a Langevin equation that describes the growth of stylolite surfaces in a quenched disordered material with long range elastic correlations.

Stylolites are geological patterns that are very common in polished marbles, a material largely used for floors and walls of buildings and monuments. They exist in many sedimentary rocks such as limestones, sandstones or evaporites Bathurst71 . They are rock-rock interfaces, that are formed during diagenesis and result from combined stress-induced dissolution and precipitation processes Dewers90 ; Ortoleva94 ; Park68 . They exist on a very large range of scales, from micro-meters to tens of meters. They are often observed as thin irregular interfaces that look like printed lines on rock cuts which explain their name. At larger scales, they are planar structures that are typically perpendicular to the tectonic load (i.e. lithostatic pressure). They might form complex network with or without connections.

However, despite their abundance, stylolites are as mentioned by Gal et al. Gal98  “among the least well-explained of all pressure-solution phenomena”. First they are complex 3D structures that are often only described from 2D cross-sections. Second, they develop in various petrological and tectonic contexts which lead to very different geometries. Third they are often transformed and deformed because of post-processes like diagenesis and metamorphism that develop after their initiation.

In this letter we propose a new mechanism for stylolites growth that is consistent with recent 3D roughness measurements of a stylolite interface. The first part of the letter deals with the topography measurement and its description in terms of scaling invariance, namely self-affinity. The second part is devoted to a physical model of the stylolites roughening. It is based on a Langevin equation that accounts for a stress-induced instability in a quenched disorder with capillary effects and long range elastic interactions.

The roughness measurement has been performed on a stylolite interface included in a limestone sample from Juras Mountains in France (see Fig. 1). The sample has been collected in a newly open quarry, thus preserved from late breakage and chemical erosion. The opening procedure was possible for this sample because of the accumulation of undissolved minerals like clays that form a weak layer along the stylolite interface. The concentration of these minerals provides an estimate of the cumulative strain that the sample has undergone Renard03 . As shown in Fig. 1, peaks along the interface are randomly distributed in space and of various sizes (up to several millimeters). Locally slopes along the topography might be very important which makes the roughness measurement difficult.

We used two different profilers to sample the stylolite roughness. First, with a mechanical profiler Schmittbuhl95 ; Lopez98  we extracted four profiles of 1030 points each with a horizontal step of $`\mathrm{\Delta }x=30\mu `$m. The specificity of the mechanical profiler is to measure the surface height from the contact of a needle onto the surface. The radius of the needle tip is 25 $`\mu `$m. The vertical resolution is $`3\mu `$m over the available range of 5 $`cm`$. One of this profile is shown in Fig. 2. We compare the mechanical measurement to an optical profiling Meheust02 . This technique is based on a laser triangulation of the surface without any contact with the surface. The laser beam is $`30\mu `$m wide. Horizontal steps between measurement points were $`\mathrm{\Delta }x=\mathrm{\Delta }y=50\mu m`$. The main advantage of this technique is the acquisition speed that can be significantly larger compared to the mechanical profiler since there is no vertical move. However, a successful comparison with mechanical measurement is necessary to ensure that optical fluctuations are height fluctuations and not material property fluctuations.

Fig. 3 shows a gray level map of the stylolite surface heights as a high resolution image ($`600\times 600`$ pixels) that covers over a surface area of 3$`\times `$3 cm2. The vertical resolution is $`2\mu `$m over the $`6`$ mm vertical range of the setup.

We analyzed the height distribution in terms of self-affinity Barabasi95  which states that the surface remains statistically unchanged for the transform: $`\mathrm{\Delta }x\lambda \mathrm{\Delta }x`$, $`\mathrm{\Delta }y\lambda \mathrm{\Delta }y`$, $`\mathrm{\Delta }z\lambda ^\zeta \mathrm{\Delta }z`$, where $`\lambda `$ can take any real value. The exponent $`\zeta `$ is the so-called roughness exponent. A 1D Fourier spectrum of a self-affine profile is shown Barabasi95  to behave as a power law with a slope $`12\zeta `$, and provides an estimate of the roughness exponent $`\zeta `$. We applied this technique whose reliability has been tested Schmittbuhl95b , to profiles extracted from the optical map shown in Fig. 3, either along $`x`$ or $`y`$ directions. Results are plotted in Fig. 4. The figure shows with a thick solid line the $`y`$-direction average spectrum of 600 profiles extracted along the $`x`$-direction. The spectrum clearly exhibits two regimes. At small wavenumbers (i.e. large length scales), a power law behavior is observed with a slope in the log-log plot of -2 which is consistent with a roughness exponent of $`\zeta _1=0.5`$. At large wavenumber (i.e. small length scales), a second power law behavior is observed with a larger slope (-3.2) in agreement with a roughness exponent $`\zeta _2=1.1`$. The crossover length scale is sharp and defines a characteristic length scale $`\lambda _c1`$mm. The average spectrum of profiles taken along the perpendicular direction ($`y`$-direction) provides a very consistent result and confirms the isotropy of the surface in terms of scaling invariance. We checked in Ref. Renard03  that mechanical profiles show the same properties specially at large wavenumbers since they are sampled at a higher resolution and allow an extended description of the unusual high roughness exponent $`\zeta _2=1.1`$. We also checked that another analysis technique, namely the Average Wavelet Coefficient technique Simonsen98 , was providing very consistent results. We also show in Ref. Renard03  that this behavior can be observed in several other stylolite samples. The difference between samples comes from the crossover length scale that can be larger or smaller, changing drastically the aspect of the surface. However, the two power law behaviors are conserved with the same values of the roughness exponents.

The second part of the letter is devoted to a modeling of the stylolite roughening in order to understand the origin of the self-affine behaviors and the parameter sensitivity of the characteristic length $`\lambda _c`$. As usually done, we reduce the problem to a 2D situation $`z(x)`$, i.e. invariance along the $`y`$-direction.

We propose a phenomenological approach that follows the work reported in Kassner et al. Kassner01  who studied the dynamics of a strained solid in contact with its melt and more specifically the Asaro-Tiller-Grinfeld instability. In the case of a solid/fluid interface, the chemical potential between the solid and the fluid at the boundary can be written as Asaro72 ; Gal98 :

$$\mathrm{\Delta }\mu =\mathrm{\Omega }(u_e+\gamma \kappa )$$
(1)
where $`u_e`$ is the elastic energy per unit volume in the solid, $`\gamma `$ is the surface energy, $`\kappa `$ the curvature and $`\mathrm{\Omega }`$ a molecular volume. We have assumed that gravity effects are negligible. If a bulk diffusion holds in the fluid, the evolution of the interface is directly related to the chemical potential: $`v_n=m\mathrm{\Delta }\mu `$ where $`v_n`$ is the normal velocity and $`m`$ is the mobility.

Following the work of Grinfeld and Hazzledine Grinfeld96 , we propose to extend the approach to the case of a solid/solid interface. Moreover, we shall also assume the presence of a fluid in the pores. Porosity of the material is supposed to be sufficiently high for having both a bulk diffusion in the fluid and capillary effects. Doing so, Eq.(1) holds. Since a complete match between both solid surfaces is assumed, the normal velocity of the interface is now driven by the relative dissolution between the two solids.

The difficult task is to expand the local elastic energy along the interface that takes into account the surface corrugations. Using the representation theorem for the displacement field and the Hook’s law, it is possible to include the influence of non local elastic interactions in the limit of small perturbations of the interface through the elastic Green function as Bilby68 ; Gao89 :

$$\sigma (x)=\sigma _0\left(1+\frac{1}{2\pi }PV_{\mathrm{}}^{\mathrm{}}\frac{z(x^{})z(x)}{(x^{}x)^2}𝑑x^{}\right)$$
(2)
where $`\sigma (x)`$ is the local stress along the interface between the solid, $`\sigma _0`$ is the average uniaxial external stress or an effective stress that includes the fluid pressure, and $`PV`$ is the principal value. According to Ref. Kassner01 , the elastic energy can be approximated for plain strain by:

$$u_e\frac{(1\nu ^2)}{E}\sigma ^2(x)$$
(3)
where $`E`$ is an effective Young’s modulus and $`\nu `$ an effective Poisson coefficient.

The last aspect of the modeling concerns a description of the noise that exists in the chemical potential because of material fluctuations. We assume that this noise fluctuates on scales significantly smaller than the scales where elasticity and capillary effects are considered as for a Langevin approach Barabasi95 . However, unlike classical use of a Langevin equation, we assume here that the noise $`\eta (x,z(x))`$ is quenched to represent spatial material fluctuations, and known to have a strong influence on the scaling properties. For instance, for the Edwards and Wilkinson (EW) problem Edwards82 ; Roux94 , the roughness exponent is 0.5 in the case of an annealed noise and 1.2 for a quenched noise.

In the present framework, the interface evolution is described by the following equation:

$`{\displaystyle \frac{1}{m}}{\displaystyle \frac{dz}{dt}}`$ $`=`$ $`\mathrm{\Omega }{\displaystyle \frac{(1\nu ^2)}{E}}\sigma _0^2\left(1+{\displaystyle \frac{1}{\pi }}PV{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{z(x^{})z(x)}{(x^{}x)^2}}𝑑x^{}\right)`$ (4)
$`+\mathrm{\Omega }\gamma {\displaystyle \frac{d^2z}{dx^2}}+\eta (x,z(x))`$
which includes three driving forces: the long range elastic interactions, the local capillary effect which has a stabilizing influence, and the noise fluctuations. A crossover length scale is obtained by the balance of the average elastic and capillary terms in Eq.(4) Kassner01 :

$$\lambda ^{}=\frac{\gamma E}{\sigma _0^2(1\nu ^2)}$$
(5)
Typical values for limestones are $`E=810^{10}`$ N/m2, $`\nu =0.25`$, and $`\gamma =20`$ J/m2. If we consider a stress at 2 km depth ($`\sigma _0=40`$ MPa), we obtain $`\lambda ^{}1`$ mm. We propose that the $`\lambda ^{}`$ length scale is the crossover length scale $`\lambda _c`$ observed in Fig. 4.

Actually, the mechanical regime ($`\lambda \lambda _c`$) and the capillary regime ($`\lambda \lambda _c`$) have both been extensively explored but independently. Indeed, the capillary regime is nothing else than the EW problem in a quenched noise. Roux and Hansen Roux94  have shown in this case that the interface is self-affine with an exponent $`\zeta _21.2`$. In the mechanical regime, Eq.(4) is reduced to the quasistatic propagation of an elastic line in a disordered material. This has been largely studied Schmittbuhl95c ; Ramanathan98 ; Rosso02  and the roughness exponent has been shown to be $`\zeta _10.4`$.

We performed a numerical simulation of Eq.(4) using an even driven algorithm Schmittbuhl95c ; Schmittbuhl99  where there is no time evolution. The periodic interface is discretized in N=2048 cells. At start, the interface is flat. At each step, we are searching for the cell that exhibits the maximum speed $`dz/dt`$ according to Eq.(4). This cell is then advanced by a random amount $`dz`$ uniformly sampled in $`[0,1]`$. The local fluctuation of the chemical potential $`\eta (x,z(x))`$ is updated from a prescribed distribution chosen here to be uniform between $`[0,1]`$. Doing so, we are always dealing with the most unstable cell and let the interface grow there.

After a transient regime, we observed a stationary width of the rough interface that develops. The power spectrum of the interface position $`P(k)`$, averaged over 500 samples, is shown in Fig. 5. It reveals a behavior very similar to the stylolite measurements: A characteristic length scale defined as the cross-over between two power law regimes. The crossover length scale $`\lambda ^{}`$ is controlled by the balance between the magnitude of the mechanical and capillary effects and is compared to the charateristic length scale $`\lambda _c`$. As expected, the mechanical regime with an exponent $`\zeta _10.4`$, is dominating at small wavenumbers (i.e. large length scales). On the contrary, at large wave numbers, the capillary regime dominates with a roughness exponent close to $`\zeta _2=1.2`$. Roughness exponents obtained from the modeling are slightly different from the measurements. This might be explained by the dimension difference: experimental surfaces are full 3D interfaces, on the contrary the model is 2D. A complete 3D numerical modeling is on-going to account for this situation. However, 3D computations are much heavier because of the long delay to reach the stationary regime of the interface growth.

In conclusion, we presented a quantitative description of stylolite interfaces that is consistent with a model of interface growth. The experimental measurement is a high resolution profiling of a 3D stylolite topography. We show that the surface is self-affine but with two regimes. At small scales, the roughness exponent is unexpectedly high, $`\zeta _2=1.1`$, and consistent with a capillary dominated regime. At large scales, the stylolites morphology is controlled by mechanical effects, i.e. long range elastic stress redistribution. In this case the roughening is important with a low roughness exponent $`\zeta _1=0.5`$. The modeling is based on the description of a stress-induced instability previously reported as the Asaro-Tiller-Grinfeld instability. Such framework provides a prediction of a characteristic length $`\lambda _c`$ which is the crossover between the two scaling regimes. It is important for geological implications to note that the characteristic length $`\lambda _c`$ is very sensitive to the average stress $`\sigma _0`$. Indeed, a measurement of $`\lambda _c`$ from today roughness profiling could provide an estimate of the stress magnitude during the stylolite growth, that is, in the past. Thus stylolites could be considered as stress fossils.

We acknowledge D. Rothman, J. Rice, A. Lobkovsky, B. Evans, Y. Bernabé, B. Goffé, P. Meakin, and E. Merino for very fruitful discussions. This work was partly supported by the ACI “Jeunes Chercheurs” of the French Ministry of Education (JS), and the ATI of the CNRS (FR).


text:
# X-ray imaging and spectroscopy of (radio-quiet) AGN: highlights from Chandra and XMM-Newton

## 1 Introduction

### 1.1 X-ray probes of the black hole region of AGN

X-ray emission originates from the immediate vicinity of the black hole. The detection of luminous, hard, power-law-like X-ray emission, rapid variability, and the recent discovery of relativistic effects in the iron-K line profile provided excellent evidence for the presence of supermassive black holes (SMBHs) in active galaxies. X-ray observations currently constitute the most powerful means to explore the black hole region of AGN.

X-rays at the centers of AGN arise in the accretion-disk – corona system. On larger scales, but still within the central region, X-rays might be emitted by a hot intercloud medium at distances of the broad or narrow-line region. The X-rays which originate from the accretion-disk region are reprocessed in form of absorption and partial re-emission (e.g., George & Fabian 1991, Netzer 1993, Krolik & Kriss 1995, Collin-Souffrin et al. 1996, Komossa & Fink 1997) as they make their way out of the nucleus. The reprocessing bears the disadvantage of veiling the intrinsic X-ray spectral shape, and the spectral disentanglement of many different potentially contributing components is not always easy. However, reprocessing also offers the unique chance to study the physical conditions and dynamical states of the reprocessing material, including the outer parts of the accretion disk; the ionized absorber; the torus, which plays an important role in AGN unification schemes (Antonucci 1993); and the BLR and NLR. Detailed modeling of the reprocessor(s) is also indispensable to recover the shape and properties of the intrinsic X-ray spectrum.

### 1.2 Chandra and XMM-Newton

There are two key abilities of the new X-ray observatories XMM-Newton (e.g., Jansen et al. 2001) and Chandra (e.g., Tananbaum & Weisskopf 2001, Weisskopf et al. 2002). One is the superb spatial resolution of Chandra, providing images rich in details on the 1 arcsec-level for the first time. The other is the increase in spectral resolution of XMM and Chandra, providing X-ray spectra with such a wealth of features including many narrow absorption and emission lines, detected in X-rays for the first time.

Both missions have imaging detectors (energy bandpass $``$(0.1–10) keV) and grating spectrometers aboard. Below, a short review of results from these observatories is given. Due to space limitations it will necessarily be incomplete, both in topics covered and in citations given. My apologies in advance.

## 2 X-ray imaging spectroscopy of AGN

### 2.1 Nearby Seyfert galaxies

Among the first published AGN X-ray images taken with Chandra was that of NGC 1068 (Young et al. 2001). The image demonstrated the outstanding spatial resolution power of the X-ray observatory. While it had been noted before that the X-ray emission of NGC 1068 was extended, with Chandra for the first time a wealth of details was seen. One could actually immediately realize one was looking at a galaxy just by inspecting the X-ray image (Fig. 1). The following details could be recognized (Young et al. 2001): a bright compact core (165 pc), emission extending to the north-east coincident with the NLR, and large-scale emission partly tracing the spiral arms. Similarly spectacular is the X-ray spectrum of NGC 1068, addressed below (Sect. 3.1).

The nearby Seyfert galaxy NGC 4151 is among the best studied AGN across the electromagnetic spectrum. Again, its extended X-ray emission was noted previously (e.g., Elvis et al. 1983). Details remained unknown, though. Chandra resolved the extended emission (Ogle et al. 2000), showed that it is spatially coincident with the NLR and of high temperature.

### 2.2 Ultra-luminous infrared galaxies

Ultra-luminous infrared galaxies (ULIRGs) are characterized by their huge luminosity output in the infrared, powered by central starburst and/or AGN activity. X-rays provide an excellent tool to search for deeply obscured AGN at the centers of ULIRGs, and also allow to trace starburst-superwind activity. With Chandra, the ULIRGs Mrk 273 (Xia et al. 2002), Mrk 231 (Gallagher et al. 2002b) and Arp 220 (Clements et al. 2002) were resolved in X-rays, revealing many details. Since there is a separate review on ULIRGs (Xia, these proceedings), the rest of this section will concentrate on NGC 6240, a recent Chandra observation of which revealed some surprises: Despite decades of searches and debates whether an AGN is present at all in this source (for an overview see Sect. 2 of Komossa & Schulz 1999), a clear detection remained elusive prior to the registration of hard X-ray emission beyond 2-10 keV with ROSAT, ASCA, and particularly BeppoSAX and RXTE. Both, Beppo-Sax and RXTE did not posses spatial resolution across the field of view, but the most likely counterpart of the hard X-ray emission was argued to be NGC 6240. The unprecedented spatial resolution of Chandra, coupled with simultaneous energy information, let to the detection of luminous, hard X-ray emission and strong neutral iron lines from both cores of the galaxy (Fig. 1,2). These properties are the characteristic features of the presence of AGN, implying both cores of NGC 6240 are active (Komossa et al. 2003). The two supermassive black holes are presently separated by $``$3000 ly. The final merging of these black holes is expected to produce a strong gravitational wave signal of the kind detectable with the future space-borne observatory LISA (e.g., Danzmann 1996, Centrella 2003).

In addition to the double black hole, NGC 6240 exhibits luminous extended X-ray emission, first detected with the instruments aboard ROSAT (e.g., Schulz et al. 1998), and spatially and spectroscopically resolved with Chandra (Komossa et al. 2003). It is very likely related to starburst-superwind activity and changes its rich structure in dependence of energy. While the more extended starburst loops have a lower X-ray temperature, there is also a more compact component surrounding the two nuclei which is of higher temperature and shows a tendency for higher metal abundances than the more extended component.

Chandra observations of more ULIRGs are expected to provide a sharp and detailed look at the very centers of these galaxies in the coming years.

## 3 High-resolution X-ray spectroscopy of AGN

For the first time, high-resolution spectroscopy of AGN can be performed in X-rays. With the exception of a few high-quality ASCA spectra of nearby AGN, the number of X-ray features typically seen in AGN spectra was 2-3. This strongly increased now; typically $``$20-30 features are detected in the spectra of nearby AGN and even more in the brightest, best-studied objects. A special role is played by warm absorbers. Recognized earlier (see Komossa 1999 for a pre-Chandra review), they have now become the most important X-ray diagnostic of the central region of AGN, imprinting many absorption lines on soft X-ray spectra of Seyfert 1 galaxies. There have even been suggestions that the emission-lines detected in nearby Seyfert 2 galaxies are also related to ionized absorbers, seen from the side.

Basically, the trend emerged that Seyfert 1 X-ray spectra are dominated by absorption, while Seyfert 2 X-ray spectra are dominated by emission-lines. There are some exceptions, like NGC 4151, a classical Seyfert 1 galaxy of type 1.5 (e.g., Osterbrock & Koski 1976), which in X-rays looks more like a Seyfert 2 (e.g., Kahn et al. 2001). Only very few AGN X-ray spectra turned out to be completely featureless (e.g., Marshall et al. 2002).

### 3.1 Seyfert 2 galaxies: emission-line spectroscopy

I would like to start, again, with NGC 1068. Both, Chandra and XMM-Newton had a deep look at this galaxy (e.g., Kahn et al. 2001, Kinkhabwala et al. 2002, Brinkman et al. 2002). Its soft X-ray spectrum is extremely rich in emission lines and is dominated by H-like and He-like ions of low $`Z`$, and by Fe-L shell ions (Fig. 3). Kinkhabwala et al. (2002) inferred that the emission-line spectrum is photoionized by the nuclear continuum, and pointed out that the inferred column densities of the ions in the X-ray gas match those derived for warm absorbers detected by absorption-spectroscopy of Seyfert 1 galaxies. The relatively cold X-ray gas spatially coincides with the optical ionization cones of NGC 1068. The inferred broad distribution of ionization parameters necessary to explain the X-ray spectrum, requires the presence of a distribution of densities (several orders of magnitude) at each radius in the ionization cone (Brinkman et al. 2002, Kinkhabwala et al. 2002). It is interesting to note, that photoionization modelling of the optical radiation cone spectrum of NGC 4151 (Schulz et al. 1993) and of the optical spectra of Seyfert 2 galaxies (Komossa & Schulz 1997) reached the same conclusion (i.e., the requirement of a range in densities at fixed radius in order to reproduce the observed emission lines).

While most attention now is focussed on employing the new observatories in orbit, a number of surprises and exciting results may still linger in the archives of ROSAT, ASCA and BeppoSAX data. Making use of archival RXTE, BeppoSAX and ASCA data, Colbert et al. (2002) reported the detection of a flare in the 6.7 keV FeK emission-line of NGC 1068. The line is interpreted as arising in warm reflecting material in the center ($`<`$ 0.2 pc from the AGN core) of the galaxy.

The X-ray emission lines detected in the spectra of several more Seyfert 2 galaxies (e.g., Mrk 3 and the Circinus galaxy) and in NGC 4151, contain important information on the physical conditions in the line-emitting medium, like temperature, density, and the main gas excitation/ionization mechanism; photoionization or collisional ionization. Of particular importance in determining the main power mechanism of the lines are the Helium-like triplets, the widths of the radiative recombination continua, and the strengths of the Fe-L complexes.

The emission lines in the X-ray spectrum of NGC 4151 (Ogle et al. 2000) are narrow and exhibit properties characteristic for the narrow-line region. This is the first X-ray NLR detected.111Ultimately, the goal will be to combine IR, optical, UV and X-ray observations in a renewed attempt to finally derive the parameters the NLR of NGC 4151 and to infer the unobserved EUV-SED of this galaxy; many attempts were done during the last decade, demonstrating that several different solutions exist (e.g., Contini et al. 2002 and ref. therein).

The general trend emerged, that the extended gas in Seyfert 2 galaxies appears to be photoionized out to relatively large distances from the nucleus.

### 3.2 Seyfert 1 galaxies: absorption spectroscopy

Neutral (‘cold’) or ionized (‘warm’) gaseous material is ubiquitous in the AGN/SMBH environment, and therefore of utmost importance in understanding the AGN phenomenon, the evolution of active galaxies, their link with starburst galaxies and ULIRGs, and the X-ray background. X-ray absorption and emission features provide valuable diagnostics of the physical conditions in the X-ray gas.

Chandra spectroscopy of NGC 5548, performed with the Low Energy Transmission Grating Spectrometer, LETGS, revealed for the first time a system of deep narrow absorption lines in the X-ray spectrum of an AGN (Kaastra et al. 2000, 2002). These lines arise from ionized material in the nucleus of NGC 5548 and show a range of ionization states.

While all AGN in deep fields are exposed for megaseconds, the record for the deepest grating observation of a nearby galaxy is held by NGC 3783. The warm absorber of this source was discovered in pioneering work by Turner et al. (1993). The 900 ks Chandra spectrum (e.g., Fig. 1 of Kaspi et al. 2002) exhibits a wealth of absorption (and emission) features.

Using Chandra and XMM-Newton, absorption-complexes were detected and studied in several more AGN, including IRAS 13349+2438 (Sako et al. 2001), NGC 4051 (Collinge et al. 2001), NGC 3227 (Komossa et al. 2001), MCG $``$6-30-15 (Branduardi-Raymont et al. 2001, Lee et al. 2001), Mrk 509 (Yaqoob et al. 2002), NGC 4593 (Steenbrugge et al. 2002), Ark 564 (Marshall 2002, Matsumoto et al. 2002) and NGC 3516 (Netzer et al. 2002). Generally, the trend emerges that the warm absorbers are complex multi-component structures with a range in ionization parameters.

The detailed atomic physics involved in producing the X-ray spectra of Seyfert galaxies, and the consequences for our understanding of the central region of AGN, are still under scrutiny. The X-ray spectra of the brightest AGN are so rich, that it will take a while until all available information is extracted from them.

### 3.3 Iron lines

The iron lines beyond 6 keV are the most important diagnostics at these energies (see Fabian 2001, O’Brien et al. 2003 for reviews). If the line is formed in the inner parts of the accretion-disk, its profile reflects the general relativistic effect of gravitational redshift and the special relativistic effects of beaming and transverse Doppler effect. The best-studied iron-line case is MCG $``$6-30-15. At certain times, the red wing of MCG$``$6-30-15 is very broad, extending down to very soft energies (Wilms et al. 2001).

XMM and Chandra confirmed and further corroborated that the iron line profiles are complex and that many components of the active nucleus may contribute to iron-line emission, including likely the BLR (NGC 5548), the torus (NGC 3783, Mrk 205), the X-ray ionization cone of NGC 1068, and a contribution from the outer parts of the accretion disk (MCG$``$6-30-15). A particularly interesting case is the iron-line complex of NGC 3516. Turner et al. (2002) reported the presence of several narrow components within the FeK$`\alpha `$ profile of this galaxy. Their origin is still somewhat unclear. More details on the iron-line topic are given in the contributions by J. Wilms and T. Wang (these proceedings).

## 4 (Distant) quasars

### 4.1 BAL quasars

Broad absorption line (BAL) quasars are characterized by broad UV absorption lines. It has been suggested that these lines arise in a flow of gas which rises vertically from a narrow range of radii from the accretion disk. The flow then bends and forms a conical wind moving radially outwards (Elvis 2000). Variants of radiatively-driven disk-winds were explored (e.g., Murray et al. 1995, Proga et al. 2000, Proga 2001, Everett et al. 2002). In some of these models, an X-ray absorber shields the wind downstream from soft X-rays, allowing resonant-line driving to remain effective and accelerate the outflowing BAL wind up to $``$0.1c.

Pre-Chandra/XMM detections of BAL quasars in X-rays were rare. Generally, BAL quasars are X-ray weak, which is usually interpreted in terms of strong excess absorption (e.g., Green et al. 1995, Gallagher et al. 1999, Brinkmann et al. 1999, Brandt et al. 2000, Wang et al. 2000). Chandra provided valuable new constraints on the amount of absorption towards selected BALs (e.g., Sabra & Hamann 2001, Oshima et al. 2001, Gallagher et al. 2002, Sabra et al. 2002), and an XMM observation let to a new identification of the X-ray counterpart of the BAL quasar PHL 5200 (Brinkmann et al. 2002). Almost all data still suffer from low S/N (typically 50 to few hundred detected X-ray photons), though. There are indications that the BAL material is ionized instead of neutral. This is definitely the case for the quasar APM 08279+5255 which has the best-measured X-ray spectrum of any BAL quasar I am aware of. APM 08279+5255 is magnified by a gravitational lense and is among the most luminous objects in the universe (even after correction for lensing). A 100ks XMM-Newton observation (Fig. 4) led to the detection of a strong absorption feature of ionized iron, interpreted as K-edges, arising from a warm absorber of high column density (Hasinger et al. 2002). A recent Chandra spectrum shows similar, but not identical features, in the sense that the high-energy dip in the spectrum is interpreted/modeled in terms of two highly blueshifted (in the quasar rest frame) iron-lines rather than edges (Chartas et al. 2002). Combining the available observations, variability in the BAL-flow is implied.

The inferred abundance of Fe/O $``$ 3 $`\times `$ solar is rather high, given the high redshift of the object (Hasinger et al. 2002). Taken at face value, it implies an efficient iron production mechanism in the early universe and provides independent evidence for the presence of a cosmological constant (see Fig. 4 of Komossa & Hasinger 2003).

The unusual strength of the iron features and the indications for strong variability make APM 08279+5255 an excellent target for simultaneous, deep follow-up observations with Chandra and XMM-Newton.

### 4.2 Absorption in high-redshift quasars

Evidence for excess X-ray absorption was found in high-redshift, mostly radio-loud, quasars (e.g., Wilkes et al. 1992, Elvis et al. 1994, Schartel et al. 1997, Yuan et al. 1998, Vignali et al. 2001). The ionization state of the absorber remained largely unknown. However, there is now growing evidence that these absorbers are ionized, not neutral. As shown by Schartel et al. (1997) the spectrum and spectral changes of the high-redshift quasar PKS 2351-154 ($`z`$=2.67) are well explained by the presence of an ionized absorber of column density $`\mathrm{log}N_\mathrm{w}=22.4`$ which changes its ionization state in response to intrinsic luminosity changes of the quasar. PKS 2351-154 is one of the very few high-$`z`$ quasars which show a variable UV absorption system as well. For several years, this quasar held the record of being the most distant X-ray warm-absorber candidate known, recently exceeded by GB 1428+42 and PMN J0525-33 (Fabian et al. 2001a,b).

Among the open questions related to absorption in (distant) quasars are: What is the origin and nature of the high-z excess absorbers ? Is this material preferentially warm or cold ? Why has it been more abundant in the past ? How does it evolve ? Why is excess absorption mainly seen in high-redshift radio-loud quasars whereas a number of (non-BAL) high-z radio-quiet quasars appear to be absorption-free ? Answers to these questions are crucial for understanding the formation and evolution of AGN. Apart from measuring ionic column densities, a very interesting prospect is to determine element abundances in dust and gas at high redshift.

### 4.3 Quasars beyond redshift 5

Chandra and XMM were used to search for and successfully detect X-rays from the highest redshift quasars, identified in the course of the SDSS (e.g., Brandt et al. 2001, $`z`$=5.80; Mathur et al. 2002, Brandt et al. 2002, Bechtold et al. 2002, Schwartz 2002, $`z`$=5.82, 5.99 and 6.28) and discovered through the Chandra multiwavelength project ChaMP (Silverman et al. 2002; $`z`$=4.93), presently the most distant X-ray selected quasar. Most of the highest redshift quasars tend to be underluminous in X-rays. It is presently being investigated whether this is due to evolutionary effects in the spectrum or the amount of absorption, or something else. The detailed analysis of accreting black holes in the early universe, i.e., the highest redshift quasars, will also be a major goal of future X-ray missions.

The status of the identification of many nearby and distant AGN in deep fields, and the implications, is summarized by G. Hasinger (these proceedings) and Hasinger (2003).

## 5 Future X-ray missions

While current X-ray missions like Chandra and XMM-Newton are still expected to provide many new exciting results, the new generation of X-ray survey missions will constitute a very useful supplement, by repeatedly re-scanning the whole X-ray sky. This will lead to detections of numerous new X-ray sources, and all kinds of unusual variability events among AGN. Among the planned missions are LOBSTER (Fraser 2001), ROSITA (Predehl 2003), and MAXI (Mihara 2001) which will cover different energy bands with different spatial and spectral resolution. In the long run, follow-ups of Chandra and XMM will come: the Japanese missions ASTRO-E II and NeXT (Kunieda 2001), the US mission Constellation-X (White & Tananbaum 2001) and the European mission XEUS (Parmar 2003), the latter two presently planned to be launched in the 2012-2015 time frame. These will combine wide energy coverage with high sensitivity and spectral resolution and with medium spatial resolution.

Acknowledgments
It is a pleasure to thank the Max-Planck Gesellschaft and the Chinese Academy of Science for providing travel funds in the course of the MPG-CAS exchange program between China and Germany, and Wolfgang Voges for reading the manuscript. This and related papers are also available at: http://www.xray.mpe.mpg.de/$``$skomossa/


text:
# Comments on tails in Schwarzschild spacetimes.

## I Introduction

Waves propagating in a curved spacetime usually undergo backscatter; that leads to the emergence of two classes of effects, the so-called quasinormal modes (discovered by Vishveshvara ) and the tails. Price, who investigated the temporal behaviour of tails in 1972, has shown some kind of universality. The late tail behaviour of scalar, electromagnetic and gravitational waves happened to depend on the order of the multipole expansion . The exposition of  is somewhat confusing and that forces us to say more than would normally be needed for the purpose of a paper reporting only numerical results.

Price states If a static l-pole field is present outside the star, prior to the onset of collapse, the field will fall off as $`t^{(2l+2)}`$. If there is no field initially outside the star, but an l-pole perturbation develops during the collapse process, it will fall off as $`t^{(2l+3)}`$ at large time . That requires a comment. Since the evolution equation (2) (see below) is linear and of the form $`_t^2\mathrm{\Psi }+L\mathrm{\Psi }=0`$, where $`L`$ is a linear time-independent operator, its solutions can be superposed. That means that if there is a static solution, then it cannot influence any evolving perturbation. Irrespective of whether there are static solutions or not, the perturbation will decay in the same way. Therefore the quoted statement cannot be taken at the face value. The natural understanding would be the following: A static field  means that there is initially a field ($`\mathrm{\Psi }(t=0)0`$), but $`_t\mathrm{\Psi }=0`$ (hence data are momentarily static). The initial absence of field in turn means that $`\mathrm{\Psi }=0`$ but $`_t\mathrm{\Psi }0`$. This agrees with the standard way of posing the initial value problem in mathematical literature, and with this in mind we inspect the later papers on this subject. The notion of initially static l-pole has been clearly specified in two publications that appeared in late seventies; it becomes an initially stationary multipole of order l (that falls off like $`1/t^{2l+2}`$) . Therein the term initially stationary means momentarily static as defined above, but with $`\mathrm{\Psi }`$ satisfying special asymptotic condition. On the other hand a detailed discussion in  shows that for initially nonstatic, $`_t\mathrm{\Psi }0`$  analytic initial data satisfying certain boundary conditions (at the event horizon of a black hole and at spatial infinity) the decay of a solution agrees with that of the Price’s initially static l-pole field. This might suggest in turn that one can understand the ”static l-pole field” as a specially prepared general initial wave profile. Corresponding initial data (both $`\mathrm{\Psi }`$ and $`_t\mathrm{\Psi }`$ are nonzero) have a noncompact support (the spatial support has infinite extension, as measured in terms of the tortoise variable $`r^{}`$ – see below). The work of Gundlach et al.  numerically confirmed the analytically derived $`1/t^{2l+2}`$ decay of late tails. Let us remark an internal inconsistency – the contention must be now that the fall off actually depends on the profile of initial data, in contrast to the primary expectation that the ”tail” is a universal phenomenon due to the backscatter of the wave signal off the curvature of the geometry.

Leaving aside the area of initially static fields, let us notice that in the case of data of compact support the predictions of  are clear – solutions should have a falloff $`t^{(2l+3)}`$. That estimate applies to solutions generated by generic initial data ($`\mathrm{\Psi }0`$, $`_t\mathrm{\Psi }0`$) . From a naive inspection of the formal solution $`\mathrm{\Psi }=_tG\mathrm{\Psi }(t=0)+G_t\mathrm{\Psi }(t=0)`$ (where $`G`$ is the Green function), one expects that in the case of momentarily static initial data the fall off is faster (by $`1/t`$) than that corresponding to general initial data. It is interesting that no numerical investigation has been performed in the case of moment of time symmetry initial data of compact support. It should be noticed also that the status of the analytic investigation is still not quite satisfactory, since it consists either of heuristic analysis  or rather formal investigation of the Green function formulae (, ) . It is only recently that the (generic case) falloff $`1/t^3`$ has been rigorously proven by Dafermos and Rodnianski for the monopole mode of the real scalar field . The same conclusion is being derived, in a different (spectral) method, by Machedon and Stalker .

The main focus of our numerical work will be on finding the falloff of tails generated by initial data of compact support in either of the aforementioned two fundamental cases. We study the evolution of the scalar monopole ($`l=0`$), the electromagnetic dipole ($`l=1`$) and the gravitational axial quadrupole ($`l=2`$). The obtained results reveal a falloff $`1/t^{2l+3}`$ for the data $`_t\mathrm{\Psi }_l0,\mathrm{\Psi }_l=0`$, in agreement with , and $`1/t^{2l+4}`$ for the initially static initial data . This latter result on the behaviour of late tails generated by initial data $`_t\mathrm{\Psi }_l=0,\mathrm{\Psi }_l0`$ is new in the numerical literature. Thus the generic initial data have the falloff $`1/t^{2l+3}`$, which agrees with Leaver , Ching et al.  and Poisson .

Our results support the view that the l-th moment of any of the waves (scalar, electromagnetic or gravitational) will decay like $`1/t^{2l+3}`$ for general initial data (in accordance with ,  and ) and $`1/t^{2l+4}`$ for the l-th moments evolving from initially static data (in agreement with the formal analysis of  and ). These results should be of practical significance for numerical relativity. The determination of the late tail behaviour is a nontrivial but at the same time feasible numerical task that can serve as a useful test for the accuracy and stability of numerical codes.

## II Definitions and equations

The spacetime geometry is defined by the line element

$$ds^2=\eta _Rdt^2+\frac{dR^2}{\eta _R}+R^2d\mathrm{\Omega }^2,$$
(1)
where $`t`$ is a time coordinate, $`R`$ is the radial areal coordinate, $`\eta _R=1\frac{2m}{R}`$ and $`d\mathrm{\Omega }^2=d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2`$ is the line element on the unit sphere, $`0\varphi <2\pi `$ and $`0\theta \pi `$. Throughout this paper the Newtonian constant $`G`$ and the velocity of light $`c`$ are put equal to 1.

The propagation of the scalar, the dipole electromagnetic and the axial (quadrupole) gravitational waves is given by

$$(_t^2+_r^{}^2)\mathrm{\Psi }=V\mathrm{\Psi }.$$
(2)
Here $`r^{}(R)=R+2m\mathrm{ln}\left(\frac{R}{2m}1\right)`$ is the tortoise coordinate while the potential term reads: for the $`l=0`$ mode of the scalar field

$$V(R)=2m\frac{\eta _R}{R^3},$$
(3)
for the $`l=1`$ (dipole) electromagnetic mode

$$V(R)=2\frac{\eta _R}{R^2}$$
(4)
and for the quadrupole axial mode of GW

$$V(R)=6\frac{\eta _R^2}{R^2}(1\frac{m}{R}).$$
(5)
## III Numerical schemes

We will choose following classes of initial data:

i) $`\mathrm{\Psi }(R,t=0)=\mathrm{sin}^4\left(\pi \frac{r^{}r^{}(a)}{4m}\right)`$ and $`_t\mathrm{\Psi }|_{t=0}=_r^{}\mathrm{\Psi }_{t=0}`$ for $`r^{}(r^{}(a)=r^{}(3m),r^{}(a)+4m)`$, and $`\mathrm{\Psi }(r,0)=_t\mathrm{\Psi }=0`$ for $`r^{}<r^{}(a)`$ or $`r^{}>r^{}(a)+4m`$. The expected asymptotic tail falloff should be the same as for generic initial data.

ii) Symmetric initial data with $`\mathrm{\Psi }(R,t=0)=\mathrm{sin}^4\left(\pi \frac{r^{}r^{}(a)}{4m}\right)`$ for $`r^{}(r^{}(a)=r^{}(3m),r^{}(a)+4m)`$ and $`\mathrm{\Psi }(r,0)=0`$ for $`r^{}<r^{}(a)`$ or $`r^{}>r^{}(a)+4m`$, and $`_t\mathrm{\Psi }_{t=0}=0`$ everywhere. It is shown below that the corresponding tail decays faster than that of i).

iii) $`_t\mathrm{\Psi }=\mathrm{sin}^4\left(\pi \frac{r^{}r^{}(a)}{4m}\right)`$ for $`r^{}<r^{}(a)=r^{}(3m)`$ or $`r^{}>r^{}(a)+4m`$ and $`_t\mathrm{\Psi }(r,0)=0`$ for $`r^{}<r^{}(a)`$ or $`r^{}>r^{}(a)+4m`$, while $`\mathrm{\Psi }|_{t=0}=0`$ everywhere; the tail behaviour happens to be like in i).

One finds convenient to split the wave equation (2), into the pair of first order differential equations

$$(_t+_r^{})\mathrm{\Psi }=\mathrm{\Phi }$$
(6)
$$(_t_r^{})\mathrm{\Phi }=V(r^{})\mathrm{\Psi }.$$
(7)
We use the following mesh scheme (see Fig. 1)

and two differencing schemes: the one-step implicit alghoritm and the multistep Adams-like algorithm. In the first approach we approximate the equations (7) by

$$\mathrm{\Psi }_{ab}\mathrm{\Psi }_{a1}=(\mathrm{\Phi }_{ab}+\mathrm{\Phi }_{a1})\frac{dt}{2}$$
(8)
$$\mathrm{\Phi }_{ab}\mathrm{\Phi }_{b1}=(V_{ab}\mathrm{\Psi }_{ab}+V_{b1}\mathrm{\Psi }_{b1})\frac{dt}{2}.$$
(9)
In the next step these equations are solved for $`\mathrm{\Psi }_{ab},\mathrm{\Phi }_{ab}`$,

$$\mathrm{\Psi }_{ab}=\mathrm{\Psi }_{ab}(\mathrm{\Psi }_{a1},\mathrm{\Psi }_{b1},\mathrm{\Phi }_{a1},\mathrm{\Phi }_{b1})$$
(10)
$$\mathrm{\Phi }_{ab}=\mathrm{\Phi }_{ab}(\mathrm{\Psi }_{a1},\mathrm{\Psi }_{b1},\mathrm{\Phi }_{a1},\mathrm{\Phi }_{b1}).$$
(11)
The second alghoritm is based on the Adams-Moulton corrector formula for ordinary differential equations

$$y_{n+1}=y_n+\frac{dt}{12}(5y_{n+1}^{^{}}+8y_n^{^{}}y_{n1}^{^{}})+O(dt^4)$$
(12)
Applying this formula to each of the equations (7) we get approximatelly

$$\mathrm{\Psi }_{ab}=\mathrm{\Psi }_{a1}+\frac{dt}{12}(5\mathrm{\Phi }_{ab}+8\mathrm{\Phi }_{a1}\mathrm{\Phi }_{a2})$$
(13)
and

$$\mathrm{\Phi }_{ab}=\mathrm{\Phi }_{b1}\frac{dt}{12}(5V_{ab}\mathrm{\Psi }_{ab}+8V_{b1}\mathrm{\Psi }_{b1}V_{b2}\mathrm{\Psi }_{b2})$$
(14)
Again this linear set of equations can be solved for $`\mathrm{\Psi }_{ab},\mathrm{\Phi }_{ab}`$

$`\mathrm{\Psi }_{ab}`$ $`=`$ $`\mathrm{\Psi }_{ab}(\mathrm{\Psi }_{a1},\mathrm{\Psi }_{b1},\mathrm{\Psi }_{a2},\mathrm{\Psi }_{b2},\mathrm{\Phi }_{a1},\mathrm{\Phi }_{b1},\mathrm{\Phi }_{a2},\mathrm{\Phi }_{b2})`$ (16)
$`\mathrm{\Phi }_{ab}=\mathrm{\Phi }_{ab}(\mathrm{\Psi }_{a1},\mathrm{\Psi }_{b1},\mathrm{\Psi }_{a2},\mathrm{\Psi }_{b2},\mathrm{\Phi }_{a1},\mathrm{\Phi }_{b1},\mathrm{\Phi }_{a2},\mathrm{\Phi }_{b2}).`$
Let us stress that these two alghoritms give almost the same results, hinting at the numerical stability of our methods.

## IV Numerical results

We expect the field $`\mathrm{\Psi }_l(r^{},t)`$ to behave (for a fixed $`r^{}`$) like

$$\underset{t\mathrm{}}{lim}\mathrm{\Psi }(r^{},t)=Ct^\alpha $$
(18)
Our goal is to get a value of $`\alpha `$ with a reasonable accuracy. Therefore we calculate the function

$$f(t)=\frac{d\mathrm{log}(\mathrm{\Psi }_l(r^{},t))}{d\mathrm{log}(t)}=t\frac{d\mathrm{log}(\mathrm{\Psi }_l(r^{},t))}{dt}$$
(19)
which asymptotically should be equal to $`\alpha `$.

In the first instance we analysed the case with waves, which belongs to the class i) of Sec. 3, and for the initial data of class iii) of the preceding section. The results are presented in Figs. 2 – 4 which show the temporal behaviour of the funtion $`f(t)`$ as seen at $`r^{}=r^{}(3m)+4m`$. In all figures (2 – 7) below the time is put on the x-axis, which is scaled in units of $`m`$.

The numerical time ranged from 5000 m (for the gravitational and electromagnetic waves) up to 40000 m (for the scalar waves). It turned out that the precise shape of this function, without numerical noise, can be calculated using numbers with 64 significant digits (gravitational and electromagnetic cases) and the ordinary double precision numbers (scalar case). In order to achieve this aim we have used the freely distributable ”qd” (quad precision) and ”arprec” (arbitrary precision) numerical libraries .

The high accuracy calculations are time consuming and therefore the achieved calculational time for the scalar case exceeds by a factor of 8 the evolution time for the remaining cases. The numerically achieved values of the exponent , for the general initial data i), range from 2.999 for the scalar field, through 4.99 for the electromagnetic field up to 6.99 for the gravitational waves (figures 2 – 4, solid line). Both sets of exponents is in a good agreement with 3, 5 and 7, respectively, obtained by Ching et al. . Similar results (broken line Fig. 3) have been obtained for the initial data of class iii) of the preceding section, in agreement with .

The momentarily static initial data (class ii) from Sec. 3) produced profiles $`f(t)`$ shown in Figs. 5 – 7. In this case the numerical exponents have been determined at $`r^{}(b)=r^{}(3m)+4m`$ and also at another obervation point, $`r^{}=100m+r^{}(b)`$. In the first case we obtained 4, 6 and 8 (up to the fourth digit number) while in the latter case we arrived at $`3.99`$, $`5.88`$ and $`7.84`$, for the scalar, electromagnetic and gravitational waves, respectively. These exponents seem to be stabilize at 4, 6 and 8 for the scalar, electromagnetic and gravitational fields, respectively; notice that the results detected at $`r^{}=r^{}(3m)+4m`$ converge quicker to the prospective limiting value than those taken at the point located farther. The asymptotic tail regime is evidently achieved quicker at a region closer to the horizon, which is a new feature, unknown in the existing literature. But in both detection points, the asymptotic values of $`\alpha `$ are similar and very close in the case of the scalar waves, see Fig. 5. The calculation of these exponents constitutes the main result of this paper.

Acknowledgments. This work has been supported in part by the KBN grant 2 PO3B 006 23. One of us (EM) acknowledges a discussion with Eric Poisson.


text:
# Tidal effects cannot be absent in a vacuum

## 1 Introduction

Non-conformally flat space-times for which the metric is an exact solution of the Einstein field equations

$$G_{ab}R_{ab}\frac{1}{2}Rg_{ab}+\mathrm{\Lambda }g_{ab}=T_{ab}$$
(1)
and in which there exists a family of observers with 4-velocity $`u^a`$ ($`u_au^a=1`$) such that the gravito-electric (or tidal) part of the Weyl-tensor vanishes,

$$E_{ac}C_{abcd}u^bu^d=0,$$
(2)
are called purely gravito-magnetic space-times. They are remarkable as the remaining gravito-magnetic part of the Weyl-tensor,

$$H_{ac}C_{abcd}^{}u^bu^d,$$
(3)
does not appear in the equation of geodesic deviation, which implies that in a purely gravito-magnetic vacuum a congruence of observers would exist for which the geodesic deviation would be identically zero:

$$\frac{D^2\xi }{d\tau ^2}𝐄.\xi =0$$
(4)
It has been conjectured that purely gravito-magnetic vacuum space-times simply do not exist , but so far a complete proof has not been given. A partial proof exists for the special cases where the Petrov type is D , or where the timelike congruence $`𝐮`$ is shear-free  or normal . The latter results extend earlier work on normal and shear-free congruences  and were generalised recently also to space-times in which there are less stringent restrictions on the shear and vorticity tensors . A clear indication that the field equations for a purely gravito-magnetic vacuum probably are not consistent when the congruence is geodesic, was given in , where it was shown that for a dust filled universe (hence $`\dot{𝐮}=0`$) a complicated chain of integrability conditions has to be satisfied: although the analysis was done for vanishing vorticity only, the reasoning suggested a possible way of attack. It is the purpose of the present paper to demonstrate explicitly that the equations are indeed inconsistent, at least for vacuum, with or without cosmological constant.

## 2 Relevant equations

I present below the relevant dynamical equations for a purely gravito-magnetic vacuum space-time in which the timelike congruence $`u^a`$ is geodesic. As in  I will follow the notations and conventions of the orthonormal tetrad formalism , with the coefficients $`n_{aa}`$ being redefined as follows:

$$n_{11}=(n_2+n_3)/2,n_{22}=(n_3+n_1)/2,n_{33}=(n_1+n_2)/2$$
(5)
and with the tetrad being specified as an eigenframe of $`H_{ab}`$. The system of equations being SO(3)-invariant, each triplet of equations will be represented by a single equation. The vanishing of the gravito-electric part of the Weyl-tensor can then be expressed by the 9 equations

$`E_{11}_0\theta _1\theta _1^2\sigma _{12}^2\sigma _{13}^2+\omega _2^2+\omega _3^2+2\sigma _{12}\mathrm{\Omega }_32\sigma _{13}\mathrm{\Omega }_2+{\displaystyle \frac{1}{3}}\mathrm{\Lambda }=0`$ (6)
$`E_{12}_0(\sigma _{12}+\omega _3)(\theta _1+\theta _2)(\sigma _{12}+\omega _3)(\sigma _{13}\omega _2)(\sigma _{23}\omega _1)`$
$`+\mathrm{\Omega }_1(\sigma _{13}\omega _2)\mathrm{\Omega }_2(\sigma _{23}\omega _1)+\mathrm{\Omega }_3(\theta _2\theta _1)=0`$ (7)
$`E_{21}_0(\sigma _{12}\omega _3)(\theta _1+\theta _2)(\sigma _{12}\omega _3)(\sigma _{23}+\omega _1)(\sigma _{13}+\omega _2)`$
$`+\mathrm{\Omega }_1(\sigma _{13}+\omega _2)\mathrm{\Omega }_2(\sigma _{23}+\omega _1)+\mathrm{\Omega }_3(\theta _2\theta _1)=0`$ (8)
The vanishing of the off-diagonal components of $`H_{ab}`$ on the other hand leads to

$`H_{12}_0(n_{12}+a_3)_1\mathrm{\Omega }_2\theta _1(n_{12}+a_3)(n_{23}a_1)(\sigma _{13}+\omega _2)+{\displaystyle \frac{1}{2}}n_2(\sigma _{12}\omega _3)`$
$`+\mathrm{\Omega }_1(n_{13}a_2)\mathrm{\Omega }_2(n_{23}a_1)+{\displaystyle \frac{1}{2}}\mathrm{\Omega }_3(n_1n_2)=0`$ (9)
$`H_{21}_0(n_{12}a_3)_2\mathrm{\Omega }_1\theta _2(n_{12}a_3)(n_{13}+a_2)(\sigma _{23}\omega _1)+{\displaystyle \frac{1}{2}}n_1(\sigma _{12}+\omega _3)`$
$`+\mathrm{\Omega }_1(n_{13}+a_2)\mathrm{\Omega }_2(n_{23}+a_1)+{\displaystyle \frac{1}{2}}\mathrm{\Omega }_3(n_1n_2)=0`$ (10)
Together with the Jacobi-identities (which guarantee the symmetry of $`E_{ab}`$ and $`H_{ab}`$), we obtain from these equations the evolution for $`\theta _a`$, $`\sigma _{ab}`$, $`\omega _a`$, $`a_a`$, $`n_{ab}`$ and $`n_a`$:

$`_0n_1=2_1(\omega _1+\mathrm{\Omega }_1)+2_2\sigma _{13}2_3\sigma _{12}+4(\omega _2+\mathrm{\Omega }_2)n_{13}4(\omega _3+\mathrm{\Omega }_3)n_{12}`$
$`n_1\theta _1n_2(\theta _1\theta _3)n_3(\theta _1\theta _2)+4n_{23}\sigma _{23}`$ (11)
From (2) one can eliminate the curl of the shear by using the diagonal components of $`H_{ab}`$:

$`H_{11}=_2\sigma _{13}_3\sigma _{12}+_1\omega _1{\displaystyle \frac{1}{2}}(\theta _1(n_2+n_3)\theta _2n_3\theta _3n_2)+2n_{23}\sigma _{23}`$
$`n_{12}(\sigma _{12}+\omega _3)a_2(\sigma _{13}+\omega _2)n_{13}(\sigma _{13}\omega _2)+a_3(\sigma _{12}\omega _3)`$ (12)
The remaining Jacobi identities contain spatial gradients of the kinematical scalars only and can be used to simplify the integrability conditions which result by considering the following commutators:

$$[_0,_1]\theta _2[_0,_2](\sigma _{12}\omega _3)\text{and}[_0,_1]\theta _3[_0,_3](\sigma _{13}+\omega _2)$$
(13)
One obtains then three pairs of equations,

$`(\sigma _{23}\omega _1+\mathrm{\Omega }_1)H_{11}(\sigma _{23}\omega _12\mathrm{\Omega }_1)H_{22}=0`$ (14)
$`(2\sigma _{23}+2\omega _1+\mathrm{\Omega }_1)H_{11}+(\sigma _{23}+\omega _1+2\mathrm{\Omega }_1)H_{22}=0`$ (15)
in which one recognizes the familiar relation  between $`\sigma `$, $`\omega `$ and $`𝐇`$,

$$\sigma \times 𝐇=3\omega 𝐇,$$
(16)
together with a relation between $`\sigma `$, $`\omega `$ and the rotation rate $`𝛀`$ of the $`𝐇`$-eigenframe with respect to a Fermi-propagated triad:

$$\omega _1^2+2\omega _1\mathrm{\Omega }_1=\sigma _{23}^2$$
(17)
(+ cyclic permutations). In fact the latter equations hold also in the general vacuum case, when the acceleration is non-zero! The last bit of information we need is the time evolution of the curvature, in the form

$$_0H_{11}=\theta _2(H_{33}H_{11})+\theta _3(H_{22}H_{11}),$$
(18)
which can be obtained directly by considering the $`[_0,_2](\sigma _{13}\omega _2)[_0,_3](\sigma _{12}+\omega _3)`$ commutators. Using these one can simplify the equations which result by substituting the expressions obtained from (16, 17) for $`\sigma `$ and $`𝛀`$, namely

$`\sigma _{12}=3\omega _3{\displaystyle \frac{H_{11}+H_{22}}{H_{11}H_{22}}}`$ (19)
$`\mathrm{\Omega }_1=2\omega _1{\displaystyle \frac{(H_{11}H_{22})(H_{11}H_{33})}{(H_{22}H_{33})^2}}`$ (20)
in the evolution equations for the shear. Note that $`𝐇`$ is not allowed to have equal eigenvalues . One finds then the following algebraic relations between $`\omega `$ and $`𝐇`$:

$`2\omega _3(H_{22}H_{33})(H_{11}H_{33})[\theta _1(H_{22}H_{33})+\theta _2(H_{11}H_{33})+3\theta _3(H_{11}H_{22})]`$
$`+3\omega _1\omega _2(H_{11}H_{22})(5H_{11}^2+8H_{11}H_{22}+5H_{22}^2)=0`$ (21)
Eliminating $`\theta _a`$ from the latter equation and its cyclic permutations results in

$`\left(2H_{11}^{}{}_{}{}^{2}+2H_{11}H_{22}+5H_{22}^{}{}_{}{}^{2}\right)\left(H_{11}H_{33}\right)^4\omega _{1}^{}{}_{}{}^{2}\omega _{3}^{}{}_{}{}^{2}`$
$`+\left(2H_{22}^{}{}_{}{}^{2}+2H_{22}H_{33}+5H_{33}^{}{}_{}{}^{2}\right)\left(H_{22}H_{11}\right)^4\omega _{2}^{}{}_{}{}^{2}\omega _{1}^{}{}_{}{}^{2}`$
$`+\left(2H_{33}^{}{}_{}{}^{2}+2H_{11}H_{33}+5H_{11}^{}{}_{}{}^{2}\right)\left(H_{22}H_{33}\right)^4\omega _{3}^{}{}_{}{}^{2}\omega _{2}^{}{}_{}{}^{2}=0`$ (22)
All the coefficients in the above expression are strictly positive (as a degenerate $`𝐇`$ is not allowed), finishing the proof that a purely magnetic vacuum is inconsistent with the assumption of a geodesic congruence.
Note that the above analysis breaks down when one or two components of the vorticity vanish. Taking into account however the evolution equations for the vorticity, it follows that at least two components must vanish, say $`\omega _2=\omega _3=0`$, such that in stead of (19) and its cyclic permutations one obtains only a single relation (19). Substituting this in the evolution for the shear leads then to a single algebraic relation between the $`H_{ab}`$ and $`\theta _a`$, a further time derivative of which is needed to obtain an inconsistency between the signs of the involved $`\omega _1^2`$ and curvature terms.

## References