text:
# Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis

## 1 INTRODUCTION

The synchrotron emission from extragalactic radio sources is a function of both the local magnetic fields and the relativistic particle populations residing within. These two components are important to the energy budget of such objects. So, pinning down their relative contributions is crucial to understanding their overall behavior. Unfortunately, the optically-thin synchrotron emission alone cannot be used to extract the individual particle and field components. However, it has long been known that in principle radio synchrotron observations can be combined with X-ray observations of inverse-Compton scattered cosmic microwave background photons (hereafter IC/3K) to extract information about particles and fields in emission regions (e.g., Harris & Romanishin, 1974; Cooke et al., 1978; Fabbiano et al., 1979; Harris & Grindlay, 1979). The advent of the Chandra and XMM-Newton observatories has made this kind of analysis possible for a large number of classes and objects.

X-ray emission has now been detected in the jets, hotspots, and lobes of numerous radio galaxies and attributed to a host of different physical processes ranging from synchrotron emission (e.g., Wilson et al. (2001)) to inverse Compton scattered emission off of one or more of several photon fields. Brunetti et al. (1999) and Setti et al. (2002) reported the detection of X-ray emission in the lobes of 3C 219 and in 3C 215 and 3C 334 respectively, that possibly arises from inverse-Compton scattering of IR photons from the quasar nucleus. IC/3K lobe emission has been reported in several radio galaxies, including Fornax A (Kaneda et al., 1995; Feigelson et al., 1995), Centaurus B (Tashiro et al., 1998), Abell85 0038-096 (Bagchi et al., 1998), 3C 295 (Brunetti et al., 2001), 3C 330 (Hardcastle et al., 2002), and 3C 263 (Crawford & Fabian, 2003). Synchrotron self-Compton (SSC) emission was reported in the hotspots of several powerful galaxies, including Cygnus A (Harris et al., 1994), 3C 295 (Harris et al., 2000), 3C 273 (Röser et al., 2000), 3C 123 (Hardcastle et al., 2001), 3C 207 (Brunetti et al., 2002), and 3C 263 and 3C 330 (Hardcastle et al., 2002). Nonthermal X-rays of uncertain origin have recently been detected in the hotspots of 3C 280 and 3C 254 (Donahue et al., 2003).

When combined with radio observations, X-ray detections allow one to infer the magnetic field strength in an emitting region. Yet, in practice the derived fields vary significantly, depending on the assumed mechanism for X-ray emission. So, interpretation is often uncertain. Indeed, closeness to the equipartition value is sometimes used as a validation criterion for a magnetic field estimated by any other means. That seems unacceptable, since there is no convincing theoretical argument for equipartition between radio-emitting electrons and magnetic fields, and empirical evidence argues that not all sources are in equipartition (e.g., Centaurus B (Tashiro et al., 1998), PKS 0637-752 (Schwartz et al., 2000)).

There are obvious complications with all of the field measures. For example, deconvolution of the particle and field information also requires assumptions about the particle and field filling factors, as these values are impossible to extract from observations. For simplicity, a uniform magnetic field distributed through the emitting region is customarily assumed. As explained below, caution is needed, particularly when assuming magnetic field uniformity. It is important to be aware that the radio and X-ray emissions may be dominated by physically different regions in the source, and so one must be careful that the same particle population is sampled by the radio and X-ray observations. This paper is based on the expectation that synthetic observations of numerically simulated radio galaxies may help us understand these issues in a way that makes minimal use of convenient, simplifying mathematical assumptions, while retaining the benefit of complete knowledge of the actual physical conditions being observed.

Tregillis et al. (2001b) (Hereafter TJR01) recently carried out three dimensional time dependent MHD simulations of radio-jet flows that included nonthermal relativistic electron transport in space and momentum, enabling them to create the first synthetic radio observations from nonthermal electron distributions that were consistently evolved within the plasma flows. Those synthetic observations were used in conjunction with a detailed dynamical analysis to study how the dynamical and nonthermal particle transport processes lead to observable radio synchrotron surface brightness and spectral index patterns. Preliminary results of the synthetic observations were also reported in Tregillis et al. (2001a), Tregillis et al. (2002a), and Tregillis et al. (2002b).

Here we look anew at the simulation data presented in TJR01, shifting our attention from a dynamical analysis to the extraction of physical properties of the simulated objects using standard radio and X-ray observational analyses and then comparing results to the actual physical properties of the simulated objects. Through this effort we aim for new insights into how well standard analysis assumptions and techniques work to capture the true physical nature of a source. In addition, we attempt to identify crucial issues for successful extraction of physical properties of a given object. We emphasize that our purpose here is not to simulate specific real world radio galaxies. Neither do we intend to reevaluate the observations and analyses of particular sources. Additional details and discussion can be found in Tregillis (2002).

The paper is structured as follows. § 2 provides a brief exposition of our numerical methods, a review of the model parameters from TJR01, and an overview of our synthetic observation techniques. These techniques are then applied to the models in §§ 3-5. In § 3 we look at what can be learned from the correlations between X-ray and radio surface brightness. § 4 is devoted to an examination of the magnetic field inferred from synchrotron and inverse-Compton surface brightness, $`B_{ic}`$ and the minimum-energy field, $`B_{me}`$. Combined use of $`B_{ic}`$ with the minimum-energy field $`B_{me}`$ as a possible tool to estimate local nonthermal energy content is discussed in § 5. The key findings are then summarized in § 6.

## 2 METHODS

### 2.1 Numerics

Our numerical methods were detailed in Jones et al. (1999) and TJR01. In short, we follow the bulk flow through a 3D TVD Eulerian ideal MHD scheme and evolve a passive relativistic electron population on the Eulerian grid through the standard particle kinetic equation.

Our MHD code is based on an extension of the “Total Variation Diminishing” scheme (Harten, 1983), as detailed in Ryu & Jones (1995) and Ryu et al. (1995). The code preserves $`B=0`$ at each time step using a constrained transport scheme (Dai & Woodward, 1998; Ryu et al., 1998). We include a passive “mass fraction” or “color tracer”, $`C_j`$, to distinguish material entering the grid through the jet orifice ($`C_j=1`$) from ambient plasma ($`C_j=0`$).

Our electron transport scheme takes practical advantage of the mismatch between bulk flow and diffusive transport scales for GeV electrons of relevance to radio and X-ray emissions within radio galaxies. In short, the lengths and times appropriate to the dynamics are orders of magnitude larger than those for electrons at energies relevant to radio synchrotron radiation. The electron momentum (energy) distribution, $`f(p)`$, is then sufficiently broad that it can be adequately represented by a piecewise power-law form within a few broad bins in momentum space. At shocks, rapid diffusive acceleration for $`10`$ GeV electrons ensures that they will emerge “instantaneously” from shocks with power-law momentum distributions. Subsequent, downstream cooling can be treated in a straightforward way. We therefore divide the momentum domain into a small number $`N`$ of logarithmically spaced bins, and estimate particle fluxes across momentum bin boundaries by representing $`f(p)p^{q(p)}`$ within bins, where $`q(p)`$ varies in a regular way (Jones et al., 1999).

In the test-particle limit for diffusive shock acceleration used here, electrons emerge from shocks with a power-law spectral index $`q=3\sigma /(\sigma 1)`$, where $`\sigma `$ is the shock compression ratio. The magnetic field is sufficiently weak in these simulations that all strong shocks behave dynamically essentially as hydrodynamic shocks. In accord with current understandings of collisionless shocks (e.g., Kang et al. (2002)) we also inject and accelerate electrons from the thermal plasma passing through shocks. Injection physics at shocks is still poorly understood, so we simply assume that a small fraction, $`ϵ`$, of the thermal electron flux through a shock becomes extended via diffusive shock acceleration into the aforementioned power-law beginning at momenta just above the postshock electron thermal values.

### 2.2 Simulated Jet Properties

Here we give a short rundown of the physical parameters of the three simulations first introduced and analyzed in TJR01. The models are summarized in Table 1, which is reproduced from TJR01.

The simulated MHD flows of TJR01 are dynamically identical. Each jet entered the grid with a simple “top hat” velocity profile with $`M_j=u_j/c_a=80`$, where $`c_a`$ is the sound speed in the uniform ambient medium. The jets entered the grid at $`x=0`$ in gas pressure balance with the ambient medium and with a density contrast $`\eta =\rho _j/\rho _a=10^2`$, giving a jet-based Mach number $`M_j=8`$. The initial jet core radius was $`r_j=15`$ zones, while the entire $`576\times 192\times 192`$ uniformly zoned grid is $`38\frac{2}{5}r_j\times 12\frac{2}{5}r_j\times 12\frac{2}{5}r_j`$. In units of initial jet radius ($`r_j=1`$) and ambient sound speed ($`c_a=(\gamma P_a/\rho _a)^{1/2}=1`$, with $`\gamma =5/3`$) the simulations stopped at $`\tau _{end}=5.4`$ time units when the bow shock reached the boundary at $`x=38\frac{2}{5}`$. In physical units, the jet radius $`r_j=2`$ kpc and the inflow velocity $`u_j=0.05c`$. This leads to a physical time unit $`10^7`$ years and a computational grid length $`77`$ kpc. Open boundary conditions were used everywhere except at the jet orifice. The initial axial and background magnetic field ($`B_x=B_{x0};B_y=B_z=0`$) (also termed “fiducial” below) exerted a magnetic pressure 1% of the ambient gas pressure ($`\beta =10^2`$). In addition to the axial component, $`B_{x0}`$, the in-flowing jet also carried a toroidal magnetic field component derived from a uniform axial current and a return current on the jet surface; i.e., $`B_\varphi =2\times B_{x0}(r/r_j)`$ for $`rr_j`$. To break cylindrical symmetry, we added a modest wobble to the in-flowing jet velocity; that is, it was slowly precessed on a cone of opening angle $`5\mathrm{°}`$ with five periods during the run.

TJR01 presented three idealized examples of electron transport, designed to isolate individual transport behaviors. Those models are also summarized in Table 1. Here, briefly, are other vital details. Electrons were transported explicitly over the momentum range $`p_0<p<p_N`$ with $`p_0=10`$ and $`p_N1.63\times 10^5`$ (with $`p`$ in units of $`m_ec`$) for all models. Below $`p_0`$ the distribution function $`f(p)`$ was continued as a power-law. At $`p_N`$ the gradient in the distribution function slope, $`dq(p_N)/d\mathrm{ln}p`$ was continued. Eight momentum bins ($`N=8`$) were used for each simulation. Since the simulated nonthermal cosmic-ray electrons were passive, all results can be scaled for other choices of $`p_0`$ ($`p_N/p_0`$ is fixed). All three transport models included adiabatic cooling and diffusive shock acceleration, although second-order Fermi acceleration and Coulomb losses were neglected. In each model, the jet nonthermal electron population entered with a momentum index $`q=4.4`$, representing a synchrotron spectral index, $`\alpha =(q3)/2=0.7`$, where $`S_\nu \nu ^\alpha `$.

As listed in Table 1 we excluded local shock electron injection in Models 1 (hereafter called the “Control Model”) and 3 (hereafter called the “Cooling Model”), setting $`ϵ=0`$. Those models isolated evolutionary behaviors of pre-existing electron populations in the jet. In Model 2 the injection parameter $`ϵ=10^4`$, so we label Model 2 the “Injection Model”. Details of these $`ϵ`$ choices are given in Jones et al. (1999) and TJR01. We note that since the electron population in the simulation was passive, the Injection Model results could be simply rescaled for an alternate $`ϵ`$. The in-flowing jet nonthermal electron population in this model was made much smaller than the other two models, so that local enhancements of the electron populations at shocks could be isolated. The Control Model and the Injection Model are “adiabatic” in the sense that electrons experience negligible synchrotron aging. That feature eliminates convex spectral curvature, although some concave curvature can result from spatial mixing of dissimilar power-law populations. Model 3, on the other hand, includes significant radiative aging from synchrotron and inverse-Compton processes, thus its Cooling Model label.

In order to parameterize radiative particle aging, we defined a characteristic synchrotron cooling time, $`\tau _{s0}`$, for electrons with $`p=10^4`$ in the fiducial magnetic field. For the adiabatic Models 1 and 2, we set $`\tau _{s0}=1.6\times 10^3`$ (compared to $`\tau _{end}=5.4`$) by setting $`B_{x0}=0.39\mu `$G and ignored inverse-Compton losses, as well, to ensure negligible radiative aging for electrons of interest. For Model 3 we made $`\tau _{s0}=\tau _{end}=5.4`$ by setting $`B_{x0}=5.7\mu `$G. Inverse-Compton losses from the CMB were taken into account by including a term with $`B_{cmb}=3.2\mu `$G corresponding to the current epoch. Once again, the physical unit for $`\tau `$ in these simulations is approximately 10 Myr.

### 2.3 Synthetic Observation Techniques

Our synthetic observation technique is straightforward. We combine vector magnetic field and nonthermal electron distribution data from our simulations to calculate self-consistent radio and nonthermal X-ray volume emissivities in every zone of the computational grid. In order to simplify the analysis we restrict emissions to zones dominated by bulk plasma originating at the jet orifice, by requiring $`C_j0.99`$. A raytracing procedure integrates along lines of sight in the optically thin limit to project emissions from the simulated objects on the sky for an arbitrary orientation. We write the output data in fits format, and analyze it using conventional observational packages (miriad and karma (Gooch, 1995)).

Since the simulations are high resolution, the synthetic maps produced are of much higher resolution and dynamic range than typical real observations. To see the influence of resolution on our analysis we also have convolved the synthesized brightness distributions with circular Gaussian beams using the miriad task CONVOL to several lower resolutions. For convenience, we placed all the objects at a fixed luminosity distance of 100 Mpc, although the choice has no influence on our conclusions. For this distance, the nominal unconvolved resolution is 0.28 arcseconds, and the projected jet length is about 110 arcseconds. We present results for convolved resolutions of 3.0 and 22.0 arcseconds as well, chosen to correspond to roughly 37 and 5 telescope beams along the jet, respectively. These choices, while arbitrary, match qualitatively what is commonly achieved in many real source observations. All synthetic observation images are set to 512 $`\times `$ 512 pixels (not every pixel contains a nonzero brightness value, however).

As in TJR01 we confine our discussion to one representative point in time; namely, $`t=4.0`$ in simulation units or about $`4\times 10^7`$ yr. All synthetic observations in this chapter refer to the same source orientation on the sky as that used in TJR01, in which the jet axis is about 45 degrees from the plane of the sky. Our conclusions are unaffected by these choices. The bulk flow is nonrelativistic, so no Doppler corrections have been applied.

#### 2.3.1 Radio Synchrotron Emissivity

In each spatial zone we compute a synchrotron emissivity based on the local vector magnetic field, B, and nonthermal electron distribution, $`f(p)`$, as evolved by our transport scheme. As given by Jones et al. (1974), the emissivity is

$$j_s(\nu )=j_{\alpha 0}\frac{4\pi e^2}{c}f(p_s)p_s^q\left(\frac{\nu _B_{}}{\nu }\right)^\alpha \nu _B_{}.$$
(1)
The spectral index $`\alpha `$ is related to the local electron momentum index $`q`$ via $`\alpha =(q3)/2`$, $`\nu _B_{}=eB\mathrm{sin}\mathrm{\Omega }/(2\pi m_ec)`$, where $`\mathrm{\Omega }`$ projects the local field onto the sky, and $`j_{\alpha 0}`$ is an order-unity dimensionless constant, defined in Jones et al. (1974). For a selected observing frequency, $`\nu `$, the distribution, $`f(p_s)`$, and the index, $`q`$, are determined for each point on the grid by establishing the relevant electron momentum from the relation $`p_s=[2\nu /(3\nu _B_{})]^{(1/2)}`$, with $`p_s`$ in units $`m_ec`$. We note for the magnetic fields in our simulated objects and radio frequencies in the GHz band that typical $`p_s10^410^5`$.

#### 2.3.2 X-ray Emissivities

We compute an X-ray emissivity, $`j_X`$, including inverse-Compton contributions from the CMB (Hereafter IC/3K) and synchrotron self-Compton (SSC); viz, $`j_X=j_{3K}+j_{SSC}`$. We ignore inverse-Compton scattered AGN photons, as well as X-ray synchrotron emission. Although we include SSC emission, since its analysis does provide some useful insights, we mention that it is typically several orders of magnitude less intense in our simulated objects than IC/3K emission.

The X-ray IC/3K emissivity, $`j_{3K}(\nu _X)`$, can be simply expressed at a selected X-ray frequency, $`\nu _X`$, in terms of the synchrotron emissivity, $`j_s(\nu _X)`$, in equation (1) extrapolated to $`\nu _X`$; namely Jones et al. (1974),

$$j_{3K}(\nu _X)=e_{\alpha 0}^{3K}\frac{\sigma _Tc}{e^2}\frac{cu_\mu \nu _\mu ^{\alpha 1}}{4\pi \nu _B_{}^{1+\alpha }}(1+z)^{3+\alpha }j_s(\nu _X),$$
(2)
where $`u_\mu =aT_0^4`$ and $`\nu _\mu =kT_0/h`$ are the energy density and characteristic frequency, respectively, at the current epoch of the CMB, while $`e_{\alpha 0}^{3K}`$ is another order-unit constant that can be obtained from Jones et al. (1974). We note in equation (2) that $`j_s(\nu _X)`$ is normalized and $`\alpha `$ is determined at $`p_\mu =(\nu _X/\nu _\mu )^{1/2}`$. For $`h\nu _X=1.27.5`$keV, considered below, $`p_\mu 2\times 10^35\times 10^3`$, which is substantially less than momenta responsible for the GHz synchrotron emission.

The SSC emissivity, $`j_{SSC}`$, depends upon the synchrotron intensity distribution incident upon each zone. To keep the calculation manageable, we adopt the common approximation of an isotropic incident intensity, so that $`j_{SSC}`$ can be expressed in terms of the omnidirectional incident flux, $`\mathrm{\Phi }_\nu ^S`$. Except when the radiation field is dominated by a very intense anisotropic local source, this approximation should be good, so adequate for our present purposes. Then, using a convenient expression in terms of the IC/3K X-ray emissivity from the same electron population we have (Jones et al., 1974),

$$j_{SSC}(\nu _X)=e_{\alpha 0}^{SSC3K}\left(\frac{\mathrm{\Phi }_{\nu _k}^S\nu _k^{\alpha _k}}{cu_\mu \nu _\mu ^{\alpha 1}}\right)\mathrm{ln}\left[\frac{2}{3}p_N^4\frac{\nu _B_{}}{\nu _X}\right]j_{3K}(\nu _X).$$
(3)
Here $`\alpha _k`$ is the synchrotron spectral index at the low-frequency synchrotron cutoff $`\nu _k`$, and $`\alpha `$ is the spectral index determined for equation (2). The constant $`e_{\alpha 0}^{SSC3K}`$ is order-unity and obtainable from Jones et al. (1974). We note that since the dominant synchrotron emission occurs at much lower frequencies than $`\nu _\mu `$, electrons responsible for SSC emission are generally of higher energy than those producing the IC/3K emission.

We use a simple and fast FFT-based scheme to estimate the omnidirectional synchrotron flux $`\mathrm{\Phi }_\nu ^S`$. Details of the calculation can be found in Appendix A.

## 3 RESULTS: RADIO AND X-RAY BRIGHTNESS CORRELATIONS

We begin with a brief, qualitative comparison between the synthetic radio and X-ray brightness distributions. Radio surface brightness maps for the three models were shown in TJR01. Here, Figures Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis through Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis show the corresponding IC/3K and SSC X-ray surface brightness images. From the images it is apparent that dominant dynamical features such as the jet and the radio hotspots are visible in both bands, but also that many details differ. For instance, X-ray brightness distributions are smoother, with much less dynamic range than the radio distributions. This simply reflects the fact that the inverse-Compton brightness represents only the column density of relativistic electrons over a narrow range of electron energies. The synchrotron brightness also depends strongly on the (intermittent) magnetic field distribution.

In §§4 and 5 we will combine the synthetic radio and X-ray data to explore their effectiveness for inferring source magnetic field properties. First, however, we demonstrate that some basic dynamical relationships connecting particles and fields can be extracted from correlations between the various intensity distributions, which are shown in Figures Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis, Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis and Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis. As we examine our three simulated objects, we should keep in mind that the magnetic field structures are identical in all three models, except for magnitude; only the nonthermal electron distributions differ. Also, we emphasize, once again, that our goal is to examine ways of extracting reliable information about the objects we have in hand; that is, how we extract from observations the physical properties that are actually present. Because we know their properties exactly, the simulated objects are uniquely valuable testing grounds, independent of how closely those detailed properties match astrophysical objects.

Figures Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis and Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis are scatter plots of 7.5 keV IC/3K and SSC intensities versus 2.9 GHz radio intensity, respectively, for the three electron transport models (Table 1). Similar distributions would be obtained for other observing frequency choices. Every other nonzero pixel is represented, for just under $`2.5\times 10^4`$ points in each plot. In each panel the X-ray and radio brightness distributions are broadly correlated. It is possible to extract from these trends insights about the dynamical behaviors of the source magnetic fields, as we illustrate. We emphasize, however, that at a given radio (X-ray) brightness, the X-ray (radio) intensities generally spread over more than an order of magnitude, so it would be unrealistic to try to predict the inverse-Compton brightness distribution from the radio distribution or vice versa. This is especially so in the lobe structure, where the spread in synchrotron brightness is enormous, due to the wide range of magnetic field strengths.

Comparison between Figures Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis and Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis shows that the SSC brightness distributions of our simulated objects are generally several orders of magnitude less than the IC/3K intensities. The difference is more than 10 orders of magnitude for both the Control Model and Injection Model, where the magnetic field magnitude was set very small to reduce synchrotron energy losses. Even for the Cooling Model, with its stronger field, the SSC intensities are mostly three to four orders of magnitude less than the IC/3K intensities. This behavior comes from the fact that the synchrotron omnidirectional flux is generally much smaller than the CMB omnidirectional flux in our simulated objects. That mirrors the situation in real radio galaxies, as well, since only the most intense X-ray hotspots have been found to be SSC dominated (e.g., Harris et al., 2000). Despite this, it is instructive to include the SSC brightness distributions in our analysis, since they reveal some useful insights.

The IC/3K vs radio trends in Figure Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis are well described for all three models by the proportionality $`I_{IC/3K}I_S^m`$, with $`m1/31/2`$. To understand this, suppose the local magnetic field scales with the electron density as $`Bn_e^b`$. Since $`I_{IC/3K}n_e𝒟`$ and $`I_Sn_e𝒟B^{1+\alpha }`$, where $`𝒟`$ is the path length, we have $`I_{IC/3K}I_S^{1/(1+b(1+\alpha ))}I_S^m`$.

On the other hand, by standard arguments $`B\rho ^{2/3}`$ for a disordered field in compression-dominated flows. Assuming $`n_e\rho `$, this gives $`m=3/(5+2\alpha )`$. For $`\alpha 0.51`$, typical of the synthetic radio sources, this would give $`m0.50.43`$. Alternatively, at perpendicular shocks one expects $`B\rho `$, resulting in $`m=1/(2+\alpha )`$, or $`m0.40.33`$. Thus, the observed trends are consistent with field evolution dominated by these dynamical processes. On the other hand, field evolution controlled by stretching of flux tubes would satisfy $`B\mathrm{}`$, where $`\mathrm{}`$ is the length of the flux tube (e.g., Gregori et al., 2000). There is no explicit interdependence between density and magnetic field. We know that the magnetic fields in the simulated object are filamented, so flux tube stretching certainly takes place. The above exercise brings out the fact that the global field evolution in the simulated object is dominated by compression, however. This matches our conclusions about the simulated field behaviors in TJR01. Sheared field evolution, while clearly involved and locally important, predominantly adds scatter to the correlations in Figure Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis.

With the notable exception of the dominant hotspot in the Injection Model, the SSC and radio intensity distributions follow similar correlations as the IC/3K-radio correlation; namely, $`I_{SSC}I_S^m`$, with $`m1/31/2`$. From this match between SSC and IC/3K trends one might expect that $`I_{SSC}I_{IC/3K}`$, and this is nicely confirmed for the Control Model in Figure Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis. This proportionality is also evident in the other two models, excepting the Injection Model hotspot, although there is a lot more scatter.

For the general distribution we note from equation 3 that one can write $`I_{SSC}\mathrm{\Phi }^SI_{IC/3K}`$. The observed correlations then imply, especially for the Control Model, that $`\mathrm{\Phi }^Sconstant`$ within the radio lobes. That represents the fact that the average omnidirectional flux for the source is a good approximation to the local flux at any point in the source. This follows if the omnidirectional synchrotron flux at each point is dominated by the collective contributions of distant emissions within the source. In that case the average synchrotron emissivity becomes an effective value to use in estimating the SSC emission from the object.

The situation is much different in the dominant Injection Model hotspot. The SSC/radio intensity correlation is described by $`m11.5`$, whereas the SSC and IC/3K intensities follow $`I_{SSC}I_{IC/3K}^{2.53}`$. On the other hand, the correlation between IC/3K and radio surface brightness is similar to what is seen elsewhere in the source and in the other two models. Previously, we interpreted the IC/3K vs radio correlation as evidence for magnetic fields evolving according to the scaling $`Bn_e^{2/3}`$, which, therefore, seems applicable in the hotspot, as well. It is the SSC behavior that is different in the Injection Model hotspot, and the distinctive correlations there can be understood by supposing that $`\mathrm{\Phi }^SI_S`$; that is, the synchrotron radiation field is dominated by local emissions in the hotspot, as we might expect. Using again the emissivity relations in §2.3 we can write in this case that $`I_{SSC}I_S^{(8+2\alpha )/(5+2\alpha )}`$ and $`I_{SSC}I_{IC/3K}^{(8+2\alpha )/3}`$. Taking $`\alpha 0.65`$ for the Injection Model hotspot radio spectral index we would predict $`I_{SSC}I_S^{1.47}`$ and $`I_{SSC}I_{IC/3K}^{3.1}`$, very close to what is observed.

Finally, we note that on the whole SSC X-ray emission is no more tightly correlated with the radio intensity than is the IC/3K intensity, despite the fact that it is a consequence of the radio emission itself. That comes from the fact that the SSC emissivity reflects radio emissions throughout the source, not just the local radio emission. Small regions of high radio brightness embedded in more diffuse emission tend to wash out the SSC/radio correlation, because they act like internal point sources that anisotropically illuminate their neighborhoods. If the surrounding emission is particularly diffuse, such regions may contribute flux at considerable distances throughout the source. That effect is particularly striking in the Injection Model, where the very bright radio hotspot is compact and surrounded by more diffuse radio emission.

## 4 RESULTS: ESTIMATED MAGNETIC FIELD STRENGTHS FROM OBSERVED INTENSITIES

As mentioned in the introduction, total synchrotron intensities are commonly applied in two ways to estimate magnetic field strengths in optically thin synchrotron sources. This is not straightforward, of course, since synchrotron intensity from an optically thin source alone cannot determine the strength of the source magnetic field. Some additional information must be supplied. The more traditional and still most-common method minimizes the total energy (or sometimes pressure) for the relativistic particles and magnetic field, or what is close to the same thing, an equipartition is assumed between them. A second, and increasingly common method combines information carried in the synchrotron intensity with observed inverse-Compton scattered intensities presumably produced by the same electron population. There are serious uncertainties in the application of both methods, ranging from the absence of convincing theoretical arguments for the minimum energy condition in an unrelaxed system to the likelihood of inhomogeneous magnetic field and particle distributions. In this section we will apply standard methods for these two approaches to synthetic observations of our simulated objects in order to see how they compare and to understand better how the inferred field properties match actual conditions for objects with complex, quasi-realistic structures.

### 4.1 The Analysis Procedure

For this exercise we “observed” synchrotron emission at $`\nu =2.9`$ GHz and IC/3K emission at 1.2 keV in all three models. To explore the connections between observationally inferred magnetic fields and actual magnetic field distributions we selected six lines of sight (LOS) for close analysis. Those LOS are marked in Figure Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis. Two of the LOS pass through lobe structures (LOS 1, LOS 4). Three penetrate the jet (LOS 3, LOS 5, LOS 6) and one centers on the dominant hotspot (LOS 2). LOS 4 allows us to look at the influence of particle aging on the analyses, since the spectrum there is strongly convex in the Cooling Model. Similarly, LOS 5 shows a concave spectrum in the Injection Model, due to strong blending of dissimilar power-law electron distributions introduced at shocks. (These spectra along LOS 4 and LOS 5 are shown in Figure Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis.) LOS 6 passes through the jet base, where magnetic field and particle populations are particularly simple, and known from the jet inflow conditions.

Magnetic field strength can be calculated very simply from synchrotron and IC/3K intensities for a uniform medium when the electron distribution is a power-law. We used the following expression for this estimated field, $`B_{ic}`$ (e.g., Jones et al., 1974; Harris & Romanishin, 1974; Harris & Grindlay, 1979):

$`B_{ic}^{1+\alpha }=(1.06\times 10^{11})(2.09\times 10^4)^{\alpha 1}(1+z)^{3+\alpha }`$
$`\times \left({\displaystyle \frac{\nu _r}{\nu _X}}\right)^\alpha \left({\displaystyle \frac{j_{\alpha 0}^{3K}}{j_{\alpha 0}}}\right){\displaystyle \frac{I_S(\nu _r)}{I_{IC/3K}(\nu _X)}}\mu \mathrm{Gauss}.`$ (4)
For uniform particle and field distributions, this expression is exactly equivalent to inverting the analytic calculation of synchrotron and IC/3K surface brightnesses in our synthetic observations.

Similarly, the standard expression used to compute the minimum-energy magnetic field, $`B_{me}`$, is (Miley (1980))

$`B_{me}=5.69\times 10^5\left[\left({\displaystyle \frac{1+k}{\eta }}\right){\displaystyle \frac{F(\nu _r)}{\nu _r^\alpha }}{\displaystyle \frac{(1+z)^{3+\alpha }}{\theta _x\theta _y\mathrm{}\mathrm{sin}^{3/2}\vartheta }}\right]^{2/7}`$
$`\times \left[{\displaystyle \frac{\nu _2^{1/2\alpha }\nu _1^{1/2\alpha }}{(1/2\alpha )}}\right]^{2/7}\mathrm{Gauss}.`$ (5)
Here $`F(\nu _r)`$ is the observed radio flux density within an observing beam, $`kU_{proton}/U_E`$, $`\eta `$ is a magnetic field volume filling factor, $`\vartheta `$ is the angle between the magnetic field and the line of sight, $`\theta _x`$ and $`\theta _y`$ are the semimajor and semiminor axes of the observing beam in arcseconds, $`\mathrm{}`$ is the path length through the source in kpc, and $`\nu _1`$ and $`\nu _2`$ are the fixed lower and upper synchrotron cutoff frequencies in the source frame, expressed in GHz. Note that while $`B_{me}`$ is expressed in terms of synchrotron flux density, it really depends on the mean $`I_S(\nu _r)/\mathrm{}`$; that is, the mean synchrotron emissivity. We assume below for simplicity that $`k=\eta =\mathrm{sin}\vartheta =1`$. The magnetic field is so tangled and intermittent that assuming $`\mathrm{sin}\vartheta =1`$ introduces errors of only a few percent into our analysis.

In order to incorporate the full nonthermal electron distributions as well as possible into $`B_{me}`$, the frequency limits needed in equation (5) correspond in each model to the characteristic synchrotron frequencies of the lowest and highest energy electrons and the fiducial magnetic field for each model (see Table 1). Those turn out to be $`\nu _1=100`$ Hz and $`\nu _2=30`$ GHz for the Control Model and the Injection Model; for the Cooling Model $`\nu _1=1500`$ Hz and $`\nu _2=450`$ GHz. Our conclusions do not depend on any of these parameter choices.

We note that since both empirical field measures depend only on intensities, they are independent of distance assumptions. Nevertheless, for consistency with our discussions in TJR01 we place the objects at a distance of 100 Mpc, giving 0.28 arcseconds for the projected size of a numerical resolution element.

Summaries of our analysis are given in Table 2 and Figure Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis. At high resolution, $`B_{ic}`$ generally falls within a factor of two of $`B_{rms}`$ in the sampled volume. $`B_{me}`$, on the other hand, can show considerably larger scatter, over an order of magnitude in some cases, depending on the specifics of the particle transport. Both estimates are sensitive to spectral curvature. Estimates for jet structures tend to show less scatter than for lobe structures, which typically embody a wider range of physical conditions. Results for 3.0 and 22.0 arcsecond Gaussian beams are also included in Table 2. The influences of bigger beams are predictable in that derived fields resemble the analogous properties of the larger regions surveyed.

Figure Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis compares $`B_{ic}`$ and $`B_{me}`$ to $`B_{rms}`$. It is evident that both empirical magnetic field estimates correlate roughly with $`B_{rms}`$ along the selected lines of sight. The bolometric synchrotron intensity does depend on the rms magnetic field along the line of sight, of course. The spectral emissivity $`j_s(\nu )B^{1+\alpha }`$, is similar, since in our sources $`\alpha 0.71`$ are common. Still it is not obvious a priori how well such a simple measure as $`B_{rms}`$ should compare to inferred values, so the experiment is valuable. Two-thirds of the $`B_{ic}`$ ($`B_{me}`$) points are within a factor of two (three) of $`B_{rms}`$. On the other hand, while the $`B_{ic}`$ values are approximately randomly distributed with respect to $`B_{rms}`$, there are obvious biases in $`B_{me}`$ that are dependent on the electron transport model in the simulation. Those biases correctly reflect actual deviations from the minimum energy condition. In particular, our simulated objects are not in electron/magnetic field energy equipartition. Nor, as we will discuss in §5, is there physics in the simulations expected to produce this kind of equipartition.

Magnetic field profiles for each of the LOS are shown in Figure Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis for the Control Model and Injection Model. The field distributions in the Cooling Model are identical to the other two sources except for being everywhere a factor 14.8 times larger. However, since the particle populations evolve differently in the Cooling Model, the inferred magnetic field behaviors are not necessarily the same as in either of the other two models. Points where $`C_j<0.99`$ are set to zero in this analysis; that is, we filtered out plasma that did not originate in the jet in order to match the assumptions made in computing nonthermal emissions. In this regard we comment that the contact discontinuity between the “jet plasma” and the “ambient plasma”; i.e., the boundary between $`C_j=1`$ and $`C_j=0`$, is not simple. It is actually quite convoluted, so that distinct fingers of ambient plasma penetrate into the cocoon. Consequently, there are “dropouts” in the emission and associated physical variables according to the prescription we have followed; that is, these regions do not contribute to $`B_{ic}`$, $`B_{me}`$, nor to $`B_{peak}`$ or $`B_{rms}`$.

In §5 we will explore how effectively $`B_{ic}`$ and $`B_{me}`$ can be used to reveal information about the partitioning of energy between magnetic fields and nonthermal electrons. Here we focus on the two field estimators separately. We summarize the important features in the next subsection.

### 4.2 Field Properties on Selected Lines of Sight

In this section we look directly at the properties of the magnetic field along the individual selected LOS as revealed in Figure Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis.

LOS 1 and LOS 4 pass through the diffuse lobe of the sources. Both lines exhibit very complex magnetic fields, and considerable entrainment of ambient plasma into the lobe structure is evident from the large number of “dropouts” in Figure Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis. The strongest field regions are revealed as the relatively bright filaments in Figures Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis-Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis and Figure Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis (see also Fig. 5 in TJR01). The nonthermal electron density also exhibits much fine structure (not shown), as one would expect from the discussion in §3. While the IC/3K-radio analysis (§3) indicates the field and particle structures are correlated, it also shows very wide scatter, especially in the lobes, where flows are chaotic. That does impact on $`B_{ic}`$ as a quantitative field measure.

Except for one case (LOS 4 in the Cooling Model) values of $`B_{ic}`$ along these lobe LOS all fall below $`B_{rms}`$. This bias results from the presence of substantial electron populations in weak field regions that contribute little to the radio emission, but that do produce X-rays. In effect the IC/3K intensity provides an overestimate of the number of radio emitting electrons, so that under the uniform source hypothesis the field apparently needed to account for the radio emissions is weakened. On the other hand, a comparison with Figure Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis shows that $`B_{ic}`$ for the Control Model and Injection Model are representative of the field values being sampled along both LOS, so it does give a “meaningful” result, if not a simply defined quantitative one.

The relatively higher value of $`B_{ic}`$ for the Cooling Model on LOS 4 is a consequence of spectral aging; that is, the spectrum is convex. Intentionally, no correction was made for this effect, in order to expose its potential influence. Recall that the 1.2 keV IC/3K emission used to compute $`B_{ic}`$ comes from 1 GeV electrons, while the 2.9 GHz emission is produced by electrons with energies about an order of magnitude higher, even in the strongest field regions along this LOS. With a convex electron spectrum the computed $`B_{ic}`$ will be artificially increased, as simple arguments can show. Suppose, for example, we measured the bolometric synchrotron intensity, $`I_S`$, from electrons of energy $`\gamma _S`$ and the bolometric IC/3K intensity, $`I_{3K}`$, from electrons of energy $`\gamma _{3K}`$. If the object were homogeneous, and we fixed the observed radio band, it is simple to show that $`B(I_S/I_{3K})N_{\gamma _{3K}}/N_{\gamma _S}`$, where $`N_{\gamma _{3K}}`$ and $`N_{\gamma _S}`$ represent the number of electrons required to produce the observed intensities. If we assumed that $`N_{\gamma _{3K}}`$ and $`N_{\gamma _S}`$ were connected by a powerlaw in $`\gamma `$ with the powerlaw index determined by the high energy electrons responsible for the radio band, but the distribution is actually convex, we would always overestimate $`B`$, since we would overestimate $`N_{\gamma _{3K}}/N_{\gamma _S}`$. This same influence has an even more striking impact on the value of $`B_{me}`$ for the Cooling Model on LOS 4. Under the equipartition that accompanies the minimum energy assumption, the effective overestimate of the electron population also exaggerates the estimated magnetic field energy, leading to an estimated magnetic field about two orders of magnitude greater than $`B_{rms}`$ and an order of magnitude greater than $`B_{peak}`$ along this LOS. Thus, it is very important when analyzing sources that exhibit “aged” spectra to account for the curvature when calculating the electron energy (I. L. Tregillis and L. Biggs, in preparation).

The values of $`B_{me}`$ found in the Injection Model along these LOS, by contrast are almost an order of magnitude below the RMS field values. This correctly reflects the fact that the magnetic field energy in this model actually does greatly exceed the nonthermal electron energy in most locations that generate synchrotron emission.

LOS 3, LOS 5 and LOS 6 intersect the jet. $`B_{ic}`$ values are close to $`B_{rms}`$ in each model when the highest resolution is used in the observations. For the Control Model and Cooling Model the inferred field is almost independent of angular resolution. The Injection Model jet is relatively less dominant in accounting for emissions, since the radiating electron population is small there. At high resolution the inferred field is also close to $`B_{rms}`$. The concave property of the synchrotron spectrum on LOS 5 in the Injection Model has at most a modest depressing influence on $`B_{ic}`$. As beam size increases on LOS 3 and LOS 5, $`B_{ic}`$ in the Injection Model increases markedly, because emissions become influenced by the nearby hotspot, where the electron population increases dramatically. That influence is greater in the Injection Model since its hotspot is much more intense than in the other two models (see Fig. Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis).

$`B_{me}`$ values follow a similar pattern to those seen in $`B_{ic}`$ for lobe LOS; that is, the relation between $`B_{me}`$ and $`B_{rms}`$ follows the actual relationship between nonthermal particle and magnetic field energies in the predominant emission regions.

LOS 6 passes through the jet near its origin. In Figure Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis the jet is contained roughly between position coordinates 245 - 295, measured in computation zones. The magnetic field structure is simple and approximates that introduced at the computational boundary. Indeed, for all three models $`B_{ic}`$ is very close to $`B_{rms}`$ as well as the projected axial base jet magnetic field in the simulation.

LOS 2 passes through a hotspot resulting from a shock produced as the jet impinges on the near boundary of its cocoon. The magnetic field is relatively compressed and organized there, accounting for the dominant peak in the field distribution near line coordinate 200. Virtually all of the radio emission in this direction originates in the hotspot, so that the $`B_{ic}`$ and $`B_{me}`$ values for all three models lie close to $`B_{rms}`$, irrespective of angular resolution. This is the only LOS we sampled that returns a $`B_{me}`$ estimate in the Injection Model that is close to $`B_{rms}`$. It is, in fact, the only LOS we sampled that has emissions in that electron transport model predominantly from regions close to electron/magnetic field equipartition.

## 5 RESULTS: OBSERVATIONAL ESTIMATION OF PARTICLE/FIELD ENERGY PARTITIONING

In the previous section we found that a synchrotron/inverse-Compton analysis provides reasonable estimates of characteristic magnetic fields, in particular $`B_{rms}`$, in our simulated objects. There are biases that can degrade the estimates when the fields are significantly more intermittent than the nonthermal electrons or when the electron energy spectra are strongly aged. However, the $`B_{ic}`$ values we derived were mostly within a factor of two of the RMS fields. On the other hand the minimum energy magnetic fields were close to the RMS fields (or any other simple statistical measure) only when the emission was dominated by relativistic plasma close to equipartition. That is what we should hope for and expect, of course. At the same time, we noted that the biases in $`B_{me}`$ correctly reflected the degree of deviation from equipartition of the dominant plasma.

This suggests that combining $`B_{ic}`$ and $`B_{me}`$ might reveal a meaningful estimate of the ratio of nonthermal electron and magnetic field energy and, consequently, provide an improved estimate of the total nonthermal energy compared to the minimum energy. We test that possibility in this section for our simulated objects, where we not only can perform radio and X-ray observations, but also know the internal energetics directly.

### 5.1 The Analysis Procedure

The obvious parameter for comparing magnetic and nonthermal particle energies is the ratio of their local densities; namely,

$$d\frac{U_B}{U_E}.$$
(6)
Standard expressions for the minimum-energy magnetic field, $`B_{me}`$, (i.e., assuming a uniform source and a power-law $`f(p)`$, and a fixed frequency band) (see eq. ), lead to the simple relation

$$d=\frac{3}{4}(1+k)\left(\frac{B}{B_{me}}\right)^{7/2},$$
(7)
where $`B`$ is the actual field strength. The exact value of $`k`$ makes no substantial difference to our conclusions, so for convenience we will still apply the commonly used value $`k=1`$. The total energy density in nonthermal particles and magnetic fields, $`U_T`$, is then obviously

$`U_T=U_B(1+{\displaystyle \frac{1+k}{d}})`$
$`={\displaystyle \frac{3}{7}}U_{me}\left[\left({\displaystyle \frac{B}{B_{me}}}\right)^2+{\displaystyle \frac{4}{3}}\left({\displaystyle \frac{B_{me}}{B}}\right)^2\right],`$ (8)
where $`U_{me}=(7/3)(B_{me}^2/(8\pi ))`$ is the combined minimum energy density.

It is convenient to define $`d_{min}\frac{3}{4}(1+k)`$ corresponding to $`B=B_{me}`$. For $`k=1`$ this gives $`d_{min}=1.5`$. Using $`d_{min}`$ we can write simply

$$d=d_{min}\left(\frac{B}{B_{me}}\right)^{7/2},$$
(9)
or

$$\frac{B}{B_{me}}=\left(\frac{d}{d_{min}}\right)^{2/7}.$$
(10)
Note for the discussion below, when $`d1`$ ($`d1`$), that $`U_Td^{4/7}U_{me}`$ ($`U_Td^{4/7}U_{me}`$).

To set the stage for our discussion we note that the jet inflow boundary conditions discussed in §2.2 give $`d_0=0.16`$ for the Control Model and the Cooling Model. Thus, these two jets enter the grid relatively close to electron/field equipartition. Accordingly, from equation (10) one finds at the jet orifice that $`B0.5B_{me}`$. In contrast, the Injection Model jet enters very far from equipartition, with $`d_0=1.6\times 10^3`$, so $`B7.3B_{me}`$.

How should we expect $`d`$ to change from these input values as the flows evolve? Recall that in our simulations particle energies are dominated by the thermal population; the nonthermal electrons are “test particles” that respond to the underlying flow dynamics, which is based on otherwise fully consistent MHD. The expected variation of $`d`$ under adiabatic compression in those flows is straightforward to estimate. Recall from §3 that compression of an isotropic magnetic field varies $`B\rho ^{2/3}`$, so that $`U_B\rho ^{4/3}`$. We found in §3 that this behavior seems to account for much of the field structure in the lobe. At the same time, the nonthermal electron gas is relativistic, so $`U_E\rho ^{4/3}`$, as well. Thus, we may expect adiabatic expansion or compression to produce relatively little change in $`d`$.

On the other hand, in sheared flows $`B`$ may have little or no correlation with $`\rho `$. In particular, an incompressible stretched flux tube has a field strength depending only on the length of the tube (see §3), which could produce very large increases in $`d`$ inside flux tubes.

Variations of $`d`$ across shocks are also difficult to quantify, although qualitatively we can expect $`d`$ to decrease there, perhaps by a large factor. At most $`U_B\rho ^2`$ (in a perpendicular shock), but $`U_E`$ should increase more strongly than this in response to diffusive shock acceleration. As shock strength increases the jump in $`U_B`$ increases asymptotically, but the jump in $`U_E`$ does not, both due to the diffusive reacceleration of the incident particle flux and also to any fresh injection due to “thermal leakage” at the shock.

### 5.2 Comparison Between Inferred and Actual Energetics

Because of initial conditions, most regions within the Control Model and Cooling Model coincidentally lie within an order of magnitude of this equipartition condition, with nonthermal electron energy being typically somewhat greater than magnetic field energy (see Figs. Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis,Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis). In contrast for much of the source volume in the Injection Model, magnetic field energy greatly exceeds the nonthermal electron energy density (Fig. Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis). That contrasting condition simply reflects the small population of nonthermal electrons entering with the jet in this model. It conveniently provides an excellent opportunity to study how different energy balances between particles and fields reveal themselves through their emission properties.

We again use the selected LOS 1 - 6 to explore how energy partitioning inferred from observations compares to actual physical conditions within the sources. Table 3 lists the RMS and mean values, $`d_{rms}`$ and $`d`$, along each line of sight, as well as the value of $`d_{obs}`$ inferred from the ratio $`B_{ic}/B_{me}`$ used in equation (9), with $`k=\eta =\mathrm{sin}\vartheta =1`$. The direct values of $`d`$ were computed from the magnetic field strength in each numerical bin and the nonthermal electron energy density integrated between $`p_0=10`$ and $`p_N=1.63\times 10^5`$, corresponding to the synchrotron cutoffs assumed in calculating $`B_{me}`$. Figure Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis shows graphically the relationship between the $`d_{obs}`$ and $`d`$ data in Table 3. For comparison, Figures Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis-Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis display distributions of $`d`$ along LOS 1-6 using the actual, full nonthermal electron distribution and the magnetic field properties from the simulations. Average values, $`d`$, are indicated in those plots by the dotted lines and $`d_{obs}`$ by dashed lines.

First looking at Figure Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis we see a rough correspondence between $`d_{obs}`$ and $`d`$ for most of the data. On the other hand agreement is not generally better than about an order of magnitude. From equation (8) that would correspond to an uncertainty in the total energy of about a factor of four, if $`U_{me}`$ were known precisely. We should not expect exact correspondence, since $`d_{obs}`$ is based on two indirect measures weighted to different points along each line of sight. We note below a couple of evident patterns that are useful to examine in order to obtain better insights to the limits of information contained in $`d_{obs}`$.

The two lobe LOS (LOS 1 and 4) in the Cooling Model give $`d_{obs}`$ values two and four orders of magnitude, respectively, below $`d`$. This strong bias is once again due to the influence of spectral aging in the lobes of this model, which have not been corrected for in this simple experiment. That is, if the electron spectrum is convex, one must be careful to count properly the low energy electrons, or $`B_{me}`$ will be seriously over estimated, which has an even larger impact on $`d_{obs}`$. Most of the other LOS in the Control Model and the Cooling Model scatter reasonably around the $`d_{obs}=d`$ line within about a factor of 4, and examination of Figures Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis and Synthetic Observations of Simulated Radio Galaxies I: Radio and X-ray Analysis shows that $`d_{obs}`$ is usually “representative”. In that regard it is worth noting along the jet LOS 3, 5 and 6 that $`d_{obs}`$ agrees well with physical values of $`d`$ in the jet itself, which is the dominant emission source. On the other hand $`d_{obs}`$ for the Control Model hotspot (LOS 2) misses the actual energy partitioning in the hotspot by about a factor of seven. Here the field geometry is more complex, so that $`B_{me}`$ is an under estimate of $`B`$ in the hotspot, even though it exceeds $`B_{rms}`$ along the LOS.

All of the $`d_{obs}`$ values in the Injection Model fall significantly below $`d`$. The closest match comes from LOS 6, which is dominated by the (simple) jet base. Otherwise, in this model there are very large fluctuations in $`d`$ that do not correspond well to the strongest emission regions. So, in this model $`d`$ is biased upwards compared to values representative of the emitting plasma and that determine $`d_{obs}`$.

The estimation of the local energy partitioning $`d`$ based on $`B_{ic}`$ and $`B_{me}`$ works best in regions that are not strongly influenced by radiative aging and where the partitioning is not expected to vary strongly over short distances, such as inside jet structures. There is increasing evidence that the X-ray emission observed from some jets is doppler-boosted IC/3K emission (e.g., Tavecchio et al., 2000; Sambruna et al., 2001; Celotti et al., 2001). In that case this kind of analysis could be useful still, but it would require taking beaming and boosting effects into account, which we have not done because the flows here are nonrelativistic. That obviously adds another level of uncertainty, since those corrections can be large. Our analysis here was directed at the simpler case, but one that still contains many real world features that must be accepted in any complex synchrotron source.

## 6 CONCLUSIONS

We have performed extensive synthetic radio and X-ray observational analyses of the numerically simulated radio galaxies introduced in TJR01. These are the first synthetic observations with sufficient detail to allow the application of standard observational techniques to numerical simulations of radio galaxies. Standard observational techniques were applied to the simulation data in order to understand better how these techniques recover physical properties of observed objects. The simulated objects have the advantage of known physical properties and evolve with many fewer simplifying assumptions than required in analytic studies of this kind. We emphasize that our goal was to compare observed properties with known physical properties, not to present the simulated objects as direct models for real objects. Our models were intentionally idealized to isolate important nonthermal particle transport behaviors. We concentrated in this paper on magnetic field strength and source energy content calculations derived from synchrotron and inverse-Compton intensities. An analysis of information extracted through polarimetry will be presented in a forthcoming companion paper.

We enumerate some of the practical messages from our analysis:

1. The synchrotron to inverse-Compton intensity ratio provides a reasonable tool for estimating magnetic field strengths in complex radio sources when the inverse-Compton photons come predominantly from the CMB. In our synthetically-observed simulated radio galaxies the standard radio/X-ray analysis returned magnetic field values that fell within about a factor two of the RMS field in a sampled volume, unless the electron spectrum was strongly convex; i.e., in the absence of strong aging. Strongly aged spectra return field values that are too high, so need to be corrected for that effect. The effectiveness of this tool seems to apply even in regions where the electron distribution and the field structure are spatially intermittent; i.e., where they have small filling factors, such as in the lobes of the simulated radio galaxies. It is largely independent of the relative partitioning of energy between electrons and magnetic field.
2. The standard synchrotron emission minimum energy analysis returns magnetic field values consistent with the RMS fields only when the actual energy partitioning between electrons and fields is close to equipartition. Otherwise the inferred fields are biased in directions that correctly reflect the actual deviations from equipartition in the regions being sampled. In addition, our analysis demonstrated that minimum energy estimates based on strongly aged spectra can seriously overestimate the actual mimimum energy magnetic field, unless the spectral curvature is properly accounted for.
3. The energy partitioning and the total source energy can be roughly estimated utilizing the ratio of the minimum energy magnetic field to that inferred from the relative synchrotron and inverse-Compton intensities. In our analysis the actual energy content was recovered to within a little better than an order of magnitude, once again in the absence of strong spectral aging. Without correction, however, total energy contents were overestimated by much more than an order of magnitude in regions with emissions dominated by strongly aged electron spectra.
4. If sources are well resolved, it may be practical to examine correlations between the radio and X-ray intensity distributions as a probe of dynamical relationships between the particles and the magnetic field. In our synthetic observation analysis of these distributions in the simulated objects we were able to recover correctly the physical correlation between magnetic field and plasma density present in the objects.
5. In regions with very large intensity contrasts, such as near a bright hotspot, smoothing at low resolution naturally leads to biases of the inferred properties in the direction of those physically in the dominant emission regions.

Finally, a note on relativistic effects. Our surface brightness computations would have to be modified in order to obtain meaningful numbers for the case of relativistic flows. However, our main purpose has not been to study the details of the surface brightness distributions so much as it has been to try to understand how well we recover meaningful information from the observations if we model them in simple but appropriate ways. It is difficult and probably of only limited usefulness to attempt a direct comparison to the relativistic case without genuine calculations. A valuable extension of this work would be to apply the same kind of analysis to appropriately-modeled emission (i.e., including beaming and boosting effects) from relativistic flows.

The work by I. L. T. and T. W. J. was supported by the NSF under grants AST96-16964 and AST00-71176 and by the University of Minnesota Supercomputing Institute. The work by D.R. was supported in part by KOSEF through grant R01-1999-00023. We gratefully acknowledge Larry Rudnick for many helpful comments and discussions, the assistance of William Ryan in developing data analysis tools, and the anonymous referee for numerous valuable suggestions. The Los Alamos National Laboratory strongly supports academic freedom and a researcher’s right to publish; therefore, the Laboratory as an institution does not endorse the viewpoint of a publication or guarantee its technical correctness.

## Appendix A CALCULATION OF OMNIDIRECTIONAL SYNCHROTRON FLUX

Assuming an isotropic emissivity, $`j_{s,\nu }j_s(\nu )`$, and an optically thin intervening plasma, the omnidirectional synchrotron flux at a given position, $`\mathrm{\Phi }_\nu ^S(𝐱)`$, is estimated as

$$\mathrm{\Phi }_\nu ^S(𝐱)=\frac{1}{4\pi }\frac{j_{s,\nu }(𝐱^{})}{(𝐱𝐱^{})^2}d^3𝐱^{}.$$
(A1)
On a uniform grid of zone size $`\mathrm{\Delta }`$, the above equation is written as

$$\mathrm{\Phi }_\nu ^S(l,m,n)=\frac{\mathrm{\Delta }}{4\pi }\underset{l^{}=0}{\overset{L1}{}}\underset{m^{}=0}{\overset{M1}{}}\underset{n^{}=0}{\overset{N1}{}}G(ll^{},mm^{},nn^{})j_{s,\nu }(l^{},m^{},n^{}),$$
(A2)
where

$$G(ll^{},mm^{},nn^{})\frac{1}{(ll^{})^2+(mm^{})^2+(nn^{})^2}.$$
(A3)
Here, $`L`$, $`M`$, and $`N`$ are the numbers of grid zones along the $`x`$, $`y`$, and $`z`$-directions, respectively. Then, from the convolution theorem, the omnidirectional flux can be computed as

$$\mathrm{\Phi }_\nu ^S(l,m,n)=\frac{\mathrm{\Delta }}{4\pi }^1\left[\widehat{G}^k(p,q,r)\widehat{j}_{s,\nu }^k(p,q,r)\right],$$
(A4)
where $`\widehat{G}^k(p,q,r)`$ and $`\widehat{j}_{s,\nu }^k(p,q,r)`$ are the Fourier transform of $`G(l,m,n)`$ and $`j_{s,\nu }(l,m,n)`$, respectively, and $`^1`$ denotes the inverse Fourier transform. For $`\widehat{G}^k(p,q,r)`$, we use

$$\widehat{G}^k(p,q,r)=\frac{2\pi ^2}{\left[\mathrm{sin}^2(2\pi p/L)+\mathrm{sin}^2(2\pi q/M)+\mathrm{sin}^2(2\pi r/N)\right]^{1/2}}.$$
(A5)
We note that the discrete Fourier transform can be used only with a periodic grid, while ours is not. However, as described, for instance, in Binney & Tremaine (1987), we can make $`G(l,m,n)`$ and $`j_{s,\nu }(l,m,n)`$ periodic by doubling our grid along each direction (that is, by extending the range of the summation into $`L`$, $`M`$, and $`N`$ in equation A2). However, since the synchrotron emissions are generally confined well inside the computational grid, aliasing is not significant. We have established empirically that there was little difference between using the doubled grid and using the original grid with periodic extensions.


format: text

text:
# Relativistic Collapse of Rotating Supermassive Stars to Supermassive Black Holes

## Introduction

There is substantial evidence that supermassive black holes (SMBHs) of mass $`10^610^{10}M_{}`$ exist and are the engines that power active galactic nuclei (AGNs) and quasars . There is also ample evidence that SMBHs reside at the centers of many, and perhaps most, galaxies , including the Milky Way .

Since quasars have been discovered out to redshift $`z>6`$ , the first SMBHs must have formed by $`z_{\mathrm{BH}}>6`$, or within $`t_{\mathrm{BH}}<10^9`$ yrs after the Big Bang. However, the cosmological origin of SMBHs is not known. This issue remains one of the crucial, unresolved aspects of structure formation in the early universe. Gravitationally, black holes are strong-field objects whose properties are governed by Einstein’s theory of relativistic gravitation — general relativity. General relativistic simulations of gravitational collapse to black holes therefore may help reveal how, when and where SMBHs, or their seeds, form in the universe. Simulating plausible paths by which the very first seed black holes may have arisen comprises a timely computational challenge (see Fig. 1). It is an area in which the tools of numerical relativity may be exploited to address a fundamental question in cosmology. Performing such simulations may also help identify plausible astrophysical scenarios and sites for promising gravitational wave sources involving black holes.

We are actively developing new algorithms and new computer codes to solve Einstein’s field equations of general relativity, coupled to the equations of relativistic hydrodynamics, in three spatial dimensions plus time (3+1). As our codes have come online, we have applied them to explore different astrophysical scenarios involving strong gravitational fields and the generation of gravitational waves. These have included the inspiral and coalescence of binary neutron stars and binary black holes, the growth of instabilities in rotating stars, the nonlinear evolution of unstable stars, and the collapse and collision of rotating stars and rotating clusters of collisionless matter, to name a few. \[For recent reviews of simulations of compact binary stars and references, see  and \]. We have also performed a wide range of simulations in recent years to study alternative scenarios leading to the formation of supermassive black holes, or their seeds. \[For a review of these calculations and references, see \] Here we will focus on one specific scenario, namely, the collapse of a rotating, supermassive star (SMS) to a supermassive black hole.

## The Onset of Dynamical Instability

The formation of SMBHs through the growth of black hole seeds by gas accretion is supported by the consistency between the total energy density in QSO light and the SMBH mass density in local galaxies, adopting a reasonable accretion rest-mass–to–energy conversion efficiency . But SMBHs must be present by $`z_{\mathrm{BH}}>6`$ to power quasars. It has been argued  that if they grew by accretion from smaller seeds, substantial seeds of mass $`>10^5M_{}`$ must already be present at $`z9`$ to have had sufficient time to build up to a typical quasar black hole mass of $`10^9M_{}`$. A likely progenitor is a very massive object (e.g., an SMS) supported by radiation pressure.

SMSs ($`10^3<M/M_{}<10^{13}`$) may form when contracting or colliding primordial gas builds up sufficient radiation pressure to inhibit fragmentation and prevent star formation (see, e.g., ). SMSs supported by radiation pressure will evolve in a quasi-stationary manner to the point of onset of dynamical collapse due to general relativity  Unstable SMSs with $`M>10^5M_{}`$ and metallicity $`Z<0.005`$ do not disrupt due to thermonuclear explosions during collapse . In fact, recent Newtonian simulations suggest that evolved zero-metallicity (Pop III) stars $`>300M_{}`$ do not disrupt but collapse with negligible mass loss . This finding could be important since the first generation of stars may form in the range $`10^210^3M_{}`$ . A combination of turbulent viscosity and magnetic fields likely will keep a spinning SMS in uniform rotation (; but see  and the final section below for an alternative). As they cool and contract, uniformly rotating SMSs reach the maximally rotating mass-shedding limit and subsequently evolve in a quasi-stationary manner along a mass-shedding sequence until reaching the instability point. At mass-shedding, the matter at the equator moves in a circular geodesic with a velocity equal to the local Kepler velocity .

It is straightforward to understand the radial instability induced by general relativity in a SMS by using an energy variational principle . Let $`E=E(\rho _c)`$ be the total energy of a momentarily static, spherical fluid configuration characterized by central mass density $`\rho _c`$. The condition that $`E(\rho _c)`$ be an extremum for variations that keep the total rest mass and specific entropy distribution fixed is equivalent to the condition of hydrostatic equilibrium and establishes the relation between the equilibrium mass and central density:

$$\frac{E}{\rho _c}=0M_{\mathrm{eq}}=M_{\mathrm{eq}}(\rho _c)(\mathrm{equilibrium}).$$
(1)
The condition that the second variation of $`E(\rho _c)`$ be zero is the criterion for the onset of dynamical instability. This criterion shows that the turning point on a curve of equilibrium mass vs. central density marks the transition from stability to instability:

$$\frac{^2E}{\rho _{c}^{}{}_{}{}^{2}}=0\frac{M_{\mathrm{eq}}}{\rho _c}=0(\mathrm{onset}\mathrm{of}\mathrm{instability}).$$
(2)
Consider the simplest case of a spherical Newtonian SMS supported solely by radiation pressure and endowed with zero rotation. This is an $`n=3`$, ($`\mathrm{\Gamma }=1+1/n=4/3`$) polytrope, with pressure

$$P=P_{\mathrm{rad}}=\frac{1}{3}aT^4=K\rho ^{\frac{4}{3}},$$
(3)
where $`K=K(s_{\mathrm{rad}})`$ is a constant determined by the value of the (constant) specific entropy $`s_{\mathrm{rad}}=\frac{4}{3}aT^3/n`$ in the star. Here $`T`$ is the temperature, $`n`$ is the baryon number density, and $`a`$ is the radiation constant. Consider a sequence of configurations with the same specific entropy but different values of central density. The total energy of each configuration is

$$E(\rho _c)=U_{\mathrm{rad}}+W,$$
(4)
where $`U_{\mathrm{rad}}`$ is the total internal radiation energy and $`W`$ is the gravitational potential energy. Applying the equilibrium condition (1) to this functional yields $`M_{\mathrm{eq}}=M_{\mathrm{eq}}(s_{\mathrm{rad}})`$, i.e. the equilibrium mass depends only on the specific entropy and is independent of central density (see Fig. 2a). Applying the stability condition (2) then shows that all equilibrium models along this sequence are marginally stable to collapse.

Now let us account for the effects of general relativity. If we include the small (de-stabilizing) Post-Newtonian (PN) correction to the gravitational field, we must also include a comparable (stabilizing) correction to the equation of state arising from thermal gas pressure:

$$P=P_{\mathrm{rad}}+P_{\mathrm{gas}}=\frac{1}{3}aT^4+2nkT,$$
(5)
where we have taken the gas to be pure ionized hydrogen. Note that $`P_{\mathrm{gas}}/P_{\mathrm{rad}}=8/(s_{\mathrm{rad}}/k)1`$. The energy functional of a star now becomes

$$E(\rho _c)=U_{\mathrm{rad}}+W+\mathrm{\Delta }U_{\mathrm{gas}}+\mathrm{\Delta }W_{\mathrm{PN}},$$
(6)
where $`\mathrm{\Delta }U_{\mathrm{gas}}`$ is the internal energy perturbation due to thermal gas energy and $`\mathrm{\Delta }W_{\mathrm{PN}}`$ is the PN perturbation to the gravitational potential energy. Applying the equilibrium condition (1) now yields $`M_{\mathrm{eq}}M_{\mathrm{eq}}^{\mathrm{Newt}}`$ times a slowly varying function of $`\rho _c`$ (see Fig. 2b). The turning point on the equilibrium curve marks the onset of radial instability; the marginally stable critical configuration is characterized by

$`\rho _{c,\mathrm{crit}}`$ $`=`$ $`2\times 10^3M_6^{7/2}\mathrm{gm}\mathrm{cm}^3,`$
$`T_{c,\mathrm{crit}}`$ $`=`$ $`(3\times 10^7)M_6^1\mathrm{K},`$ (7)
$`(R/M)_{\mathrm{crit}}`$ $`=`$ $`1.6\times 10^3M_6^{1/2},`$
where $`M_6`$ denotes the mass in units of $`10^6M_{}`$. (Here and throughout we adopt gravitational units and set $`G=1=c`$.)

Finally, let us consider a uniformly rotating SMS spinning at the mass-shedding limit . A centrally condensed object like an $`n=3`$ polytrope can only support a small amount of rotation before matter flys off at the equator. At the mass-shedding limit, the ratio of rotational kinetic to gravitational potential energy is only $`T/|W|=0.899\times 10^21`$. Most of the mass resides in a nearly spherical interior core, while the low-mass (Roche) envelope bulges out in the equator: $`R_{\mathrm{eq}}/R_{\mathrm{pole}}=3/2`$. When we include the contribution of rotational kinetic energy to the energy functional, we must now also include the effects of relativistic gravity to Post-Post-Newtonian (i.e. PPN) order, since both $`T`$ and $`\mathrm{\Delta }W_{\mathrm{PN}}`$ scale with $`\rho _c`$ to the same power. The energy functional becomes

$$E(\rho _c)=U_{\mathrm{rad}}+W+\mathrm{\Delta }U_{\mathrm{gas}}+\mathrm{\Delta }W_{\mathrm{PN}}+\mathrm{\Delta }W_{\mathrm{PPN}}+T$$
(8)
Applying the equilibrium condition (1), holding $`M`$, angular momentum $`J`$ and $`s`$ fixed, now yields $`M_{\mathrm{eq}}M_{\mathrm{eq}}^{\mathrm{Newt}}`$ times a slowly varying function of $`\rho _c`$ (see Fig. 2c). If we restrict our attention to rapidly rotating stars with $`M>10^5M_{}`$ (the typical size of the seed black holes adopted in some recent galaxy merger simulations; see, e..g., ) the influence of thermal gas pressure is unimportant in determining the critical point of instability. The turning point on the equilibrium curve then shifts to higher density and compaction than the critical values for nonrotating stars, reflecting the stabilizing role of rotation:

$`\rho _{c,\mathrm{crit}}`$ $`=`$ $`0.9\times 10^1M_6^2\mathrm{gm}\mathrm{cm}^3,`$
$`T_{c,\mathrm{crit}}`$ $`=`$ $`(9\times 10^7)M_6^{1/2}\mathrm{K},`$ (9)
$`(R_{\mathrm{pole}}/M)_{\mathrm{crit}}`$ $`=`$ $`427,`$
$`(J/M^2)_{\mathrm{crit}}`$ $`=`$ $`0.97.`$
The actual values quoted above for the critical configuration were determined by a careful numerical integration of the general relativistic equilibrium equations for rotating stars ; they are in close agreement with those determined analytically by the variational treatment. The numbers found for the nondimensional critical compaction and angular momentum are quite interesting. First, they are universal ratios that are independent of the mass of the SMS. This means that a single relativistic simulation will suffice to track the collapse of a marginally unstable, maximally rotating SMS of arbitrary mass. Second, the large value of the critical radius shows that a marginally unstable configuration is nearly Newtonian at the onset of collapse. Third, the fact that the angular momentum parameter of the critical configuration $`J/M^2`$ is below unity suggests that, in principle, the entire mass and angular momentum of the configuration could collapse to a rotating black hole in vacuum without violating the Kerr limit for black hole spin (but see the next section below!).

## The Outcome of Collapse

There are several plausible outcomes that one might envision a priori for the dynamical collapse of a uniformly rotating SMS once it reaches the marginally unstable critical point identified above. It could collapse to a clumpy, nearly axisymmetric disk, similar to the one arising in the Newtonian SPH simulations for the isothermal ($`\mathrm{\Gamma }=1`$) implosion of an initially homogeneous, uniformly rotating, low-entropy cloud . Alternatively, the disk might develop a large-scale, nonaxisymmetric bar. After all, the onset of a dynamically unstable bar mode in a spinning equilibrium star occurs when the ratio $`T/|W|0.27`$ (see, e.g.,  and  for Newtonian treatments and  and  for simulations in general relativity). Since $`T/|W|`$ is $`0.899\times 10^2`$ at the onset of collapse and scales roughly as $`R^1`$ during collapse, assuming conservation of mass and angular momentum, this ratio climbs above the dynamical bar instability threshold when the SMS collapses to $`R/M20`$, well before the horizon is reached. The growth of a bar might begin at this point. Indeed, a weak bar forms in simulations of rotating supernova core collapse , but in supernova the equation of state stiffens ($`\mathrm{\Gamma }>4/3`$) at the end of the collapse, triggering a bounce and thereby allowing more time for the bar to develop. A rapidly rotating unstable SMS might not form a disk at all, but instead collapse entirely to a Kerr black hole; not surprisingly, a nonrotating spherical SMS has been shown to collapse to a Schwarzschild black hole . Alternatively, the unstable rotating SMS might collapse to a rotating black hole and an ambient disk.

Recently we have performed two simulations that together resolve the fate of a marginally unstable, maximally rotating SMS of arbitrary mass $`M`$. In  we followed the collapse in full $`3D`$, but assumed PN theory. We tracked the implosion up to the point at which the central spacetime metric begins to deviate appreciably from flat space at the stellar center. We found that the massive core collapses homologously during the Newtonian epoch of collapse, and that axisymmetry is preserved up to the termination of the integrations. This calculation motivated us  to follow the collapse in full general relativity by assuming axisymmetry from the beginning to maximize spatial resolution (see Fig. 3). We found that the final object is a Kerr-like black hole surrounded by a disk of orbiting gaseous debris. The final black hole mass and spin were determined to be $`M_h/M0.9`$ and $`J_h/M_h^20.75`$. The remaining mass goes into the disk of mass $`M_{\mathrm{disk}}/M0.1`$. A disk forms even though the total spin of the progenitor star is safely below the Kerr limit. This outcome results from the fact that the dense inner core collapses homologously to form a central black hole, while the diffuse outer envelope avoids capture because of its high angular momentum. Specifically, in the outermost shells, the angular momentum per unit mass $`j`$, which is strictly conserved on cylinders, exceeds $`j_{\mathrm{ISCO}}`$, the specific angular momentum at the innermost stable circular orbit about the final hole. This fact suggests how the final black hole and disk parameters can be calculated analytically from the initial SMS density and angular momentum distribution . The result applies to the collapse of any marginally unstable $`n=3`$ polytrope at mass-shedding. Maximally rotating stars which are characterized by stiffer equations of state and smaller $`n`$ (higher $`\mathrm{\Gamma }`$) do not form disks, typically, since they are more compact and less centrally condensed at the onset of collapse .

The above calculations show that a SMBH formed from the collapse of a maximally rotating SMS is always born with a “ready-made” accretion disk. This disk might provide a convenient source of fuel to power the central engine. The calculations also show that the SMBH will be born rapidly rotating. This fact is intriguing in light of suggestions that observed SMBHs rotate rapidly (e.g., ).

The implosion will result in a burst of gravitational waves. Current relativistic collapse calculations for this problem break down soon after the apparent horizon of the hole approaches its asymptotic size. Moreover, the need to resolve the centrally condensed inner core precludes extending the numerical grid very far out into the wave zone. These two restrictions prevent a reliable determination of the complete burst waveform at present. Crude estimates can be generated for the frequency and amplitude based on the quadrupole formula, and the results suggest that they are well within the range of detectability for LISA for a reasonable spectrum of masses. More interesting, perhaps, is the possible generation of quasi-periodic gravitational waves from nonaxisymmetric instabilities that might arise in the ambient disk; the characteristics of these waves are likely to be comparable those discussed in the next section and in Eqn. (10) below.

## An Alternative Scenario

The previous analysis assumes that a combination of turbulent viscosity and magnetic fields combine to keep a spinning SMS rotating uniformly prior to the onset of collapse. While this is the most likely possibility , New and Shapiro  have discussed an alternative. In the absence of viscosity or magnetic fields, an evolving SMS will conserve angular momentum on cylinders as it cools and contracts in a quasi-stationary manner. So even if it is rotating uniformly at formation, it will evolve to a state of differential rotation. Our evolution calculations indicate that an axisymmetric configuration that is slowly rotating uniformly and nearly spherical at birth ($`T/|W|1`$) will evolve to a differentially rotating toroid (see Fig. 4). Prior to encountering any relativistic radial instability, the star reaches the mass-shedding limit, at which point $`T/|W|0.27`$ and differential rotation is extreme, $`\mathrm{\Omega }_{\mathrm{pole}}/\mathrm{\Omega }_{\mathrm{eq}}2.5\times 10^3`$.

Toroidal configurations which have high values of $`T/|W|`$ and substantial differential rotation and central mass concentration are subject to nonaxisymmetric dynamical instabilities, especially the one-armed spiral $`m=1`$ instability . Softer equations of state with higher polytropic indicies are very susceptible to this instability (see Fig. 5), so this mode is likely to be triggered in a rotating SMS should it evolve with low viscosity and magnetic field. Such a mode typically induces a secondary m=2 bar mode of smaller amplitude, and the bar mode can excite quasi-periodic gravitational waves (see Fig. 6). Typical wave amplitudes and wave frequencies for a source at a distance $`r`$ are given by

$`h`$ $``$ $`10^{16}(T/|W|)(M/R)M_6r_{\mathrm{GPC}}^110^{20}M_6r_{\mathrm{GPC}}^1,`$ (10)
$`f_{\mathrm{GW}}`$ $``$ $`10^1(T/|W|)^{1/2}(M/R)^{3/2}M_6^1Hz10^4M_6^1Hz.`$
where we have set $`T/|W|0.27`$ and $`R/M10^3`$ in evaluating the above quantities. These values also apply to waves generated by bar modes that may arise in a similar fashion in the ambient disk surrounding the SMS black hole that forms from the collapse of a relativistically unstable, uniformly rotating SMS (see discussion above). Both the amplitudes and frequencies put these waves well within the detectable range for LISA for a realistic distribution of masses. Future simulations that we are planning will explore this possibility in greater detail, and help resolve which scenario best describes the final fate of a rotating SMS.

## acknowledgments

This work was supported in part by NSF Grants PHY-0090310 and PHY-0205155 and NASA Grant NAG5-10781 at the University of Illinois at Urbana-Champaign.


text:
# Ionization of clusters in intense laser pulses through collective electron dynamics

## Abstract

The motion of electrons and ions in medium-sized rare gas clusters ($``$ 1000 atoms) exposed to intense laser pulses is studied microscopically by means of classical molecular dynamics using a hierarchical tree code. Pulse parameters for optimum ionization are found to be wavelength dependent. This resonant behavior is traced back to a collective electron oscillation inside the charged cluster. It is shown that this dynamics can be well described by a driven and damped harmonic oscillator allowing for a clear discrimination against other energy absorption mechanisms.

The interaction of of intense laser radiation with clusters has been of continuing interest mcth+94 ; diti+97 ; ledo+98 ; dizw+99 ; zwdi+99 ; kosc+99 ; po01 , pushed by prominent findings as the emission of keV-photons mcth+94 , highly charged ions diti+97 ; ledo+98 , or fast fragments dizw+99 . Basically, all these phenomena are caused by the exceedingly effective absorption of energy from the laser field into the cluster. This enhanced absorption (when compared to atoms or bulk matter) is possible due the initially *solid-like atomic density in the cluster* in combination with the *rapid expansion of the cluster* on a femtosecond time scale, i. e. typically during the laser pulse. More detailed insight into the mechanism of energy absorption can be gained by pulse length variation. Using pulses of equal energy, i. e. longer pulses have lower intensities, it has been found that there is an optimum pulse length with maximum absorption zwdi+99 ; kosc+99 . This holds true for small rare gas and metal clusters (of the order of some 10 atoms) as well as for large clusters, although the underlying reason for maximum absorption can be very different and is in fact an issue of current debate. For small rare gas clusters under laser pulses of peak intensities in the tunneling regime ($`I10^{15}`$ W/cm2) such an optimum could be attributed to the mechanism of enhanced ionization scro0203  known from molecules zuba95seiv+95 . For small metal clusters exposed to similar laser pulses the existence of optimal absorption was interpreted as a plasmon resonance phenomenon kosc+99  in analogy to the well known dipole resonance of the valence electrons in perturbative photo-absorption br93  or low-intensity laser pulses sure00 . For large clusters (of more than $`10^5`$ atoms) it has been proposed that strong laser pulses create a nanoplasma inside the cluster dido+96 ; mimc+01 . In the course of the expansion of the cluster the electron density and consequently the plasma frequency decreases resulting in strong energy absorption at resonance with the laser.

Separating the different mechanisms from each other requires specific and clear signatures for each process which are difficult to identify in a multi particle system such as a cluster. Surprisingly, dipole resonant absorption dynamics in a cluster can be very well characterized by a simple driven damped harmonic oscillator, which describes the dipole response of the electrons inside the cluster.

In the following we will demonstrate the validity of this simple description with full dynamical microscopic calculations for Xenon clusters ($`10^2`$…103 atoms) in strong optical laser pulses (wavelengths $`\lambda =520\mathrm{}1170`$ nm, intensities $`I10^{14}`$…1016 W/cm2, pulse lengths $`T10\mathrm{}1000`$ fs). Our approach is similar to those used before for intense laser-cluster interaction rosc+97 . However, we have been forced to use a completely new propagation scheme, namely a hierarchical tree code pfgi96 , to handle of the order of $`10^4`$ charged particles ($``$ 1000 ions and $``$ 8000 electrons) with their mutual interactions. Originally developed for gravitational $`N`$-body problems in cosmology bahu86 , the hierarchical tree code allows us to follow the dynamics of all charged particles over a few hundred femtoseconds with typical time steps of attoseconds.

Atoms are initially arranged in so-called Mackay icosahedra ho79  known to be optimal structures of Lennard-Jones clusters. We distinguish between free and bound electrons, whereby the latter ones are not treated explicitly. The condition for creation of an electron is that none of the other electrons has a negative binding energy to the ion under consideration. In this case a new electron is “born” at the position of that ion with a kinetic energy to satisfy the ionization potential. The charge of the ion is increased by one. Free electrons as well as ions are described classically as charged particles in an oscillating field interacting via a softened Coulomb interaction $`W`$. This is defined for a pair of particles with charges $`q_1`$ and $`q_2`$ and a distance of $`r_{12}`$ as $`W=q_1q_2/\sqrt{r_{12}^{}{}_{}{}^{2}+1}`$.

First of all, we present the calculated pulse length dependence for three different laser wavelengths $`\lambda `$. Figure 1 shows the final averaged charge state per atom from Xe561 clusters after laser impact as a function of the pulse length $`T`$ or the peak intensity $`I`$, respectively. The applied laser pulse was linearly polarized with a field

$$F(t)=\{\begin{array}{cc}\sqrt{I/I_0}\mathrm{cos}^2\left(\frac{\pi }{2}\frac{t}{T}\right)\mathrm{cos}(\omega t)\hfill & \text{for }|t|<T\hfill \\ 0\hfill & \text{elsewhere}\hfill \end{array}$$
(1)
where $`I_0=3.5110^{16}`$ W/cm2. In order to keep the energy of the different pulses constant we fixed the product of intensity and pulse length $`I\times T=I_0\times 4`$ fs. For short pulses ($`T\mathrm{\hspace{0.17em}30}`$ fs), where the cluster atoms have not enough time to react on the charging, the final charge state decreases with an increasing pulse length due to the lower intensity. For longer pulses ($`T\mathrm{\hspace{0.17em}50}`$ fs), however, the final charge state increases despite the intensities become smaller. This can only be understood if one considers the expansion of the cluster, for a detailed explanation see below. Finally, for very long pulses ($`T\mathrm{\hspace{0.17em}400}`$ fs) the cluster is already completely fragmented before the laser pulse reaches its peak intensity rendering the ionization similar compared to the case of single atoms. Qualitatively, this behavior is the same for all three frequencies (Fig. 1). However, the shift of the optimal pulse length towards longer pulses for longer wavelengths is characteristic for a resonant ionization mechanism.

To gain insight into the mechanism of ionization we discuss the dynamics of a Xe923 cluster in a laser pulse with $`I=910^{14}`$ W/cm2, cf. Fig. 2. The pulse has short rise and fall times of 20 fs and a long plateau of 160 fs in order to eliminate effects from the time dependence of the laser pulse itself. As can be clearly seen in Fig. 2a, the cluster ionization occurs in two steps and the total cluster charge $`Q`$ (sum of total ionic charge $`Q_{\mathrm{ion}}`$ and charge of those electrons which are inside the cluster volume) as a function of time $`t`$ is well represented by

$$Q(t)=\underset{i=1,2}{}\frac{Q_i}{1+\mathrm{exp}((tt_i)/\delta t_i)}.$$
(2)
In the 1st step during the rising of the pulse ($`t_1\pm \delta t_1=86`$ fs $`\pm `$ 3 fs) electrons are emitted mainly due to field ionization. This process slows down, however, already before the plateau intensity is reached due to the increasing space charge (cluster charge $`Q1000`$ at $`t=80`$ fs). Because of this attractive space charge, one may distinguish between inner ionization, which accounts for excitation from localized electrons to quasi-free electrons moving inside the cluster volume, and outer ionization, which corresponds to the final escape of the quasi-free electrons from the cluster into the continuum. Note, that the space charge can hold an appreciable number of quasi-free electrons which engage into collective motion discussed below, in contrast to the ionization dynamics of molecules or small clusters with almost no quasi-free electrons scro0203 . The charging up of the cluster leads to an expansion as can be seen for $`t60`$ fs in Fig. 2b. During the expansion a 2nd ionization step occurs which lasts for a much longer time and leads to an increase of the average ionic charge from about 1 to more than 6. During this time, the quasi-free electrons in the cluster are driven collectively back and forth along the polarization direction of the laser which is evident from their centre-of-mass (cm) velocity $`v_{\mathrm{CM}}`$ shown in Fig. 2c.

This oscillation can be modelled by a driven and damped classical harmonic oscillator

$$\ddot{X}(t)+2\mathrm{\Gamma }_t\dot{X}(t)+\mathrm{\Omega }_{t}^{}{}_{}{}^{2}X(t)=F(t)$$
(3)
with $`X(t)`$ the cm position of the electron cloud, $`F(t)`$ the driving laser amplitude, and $`\mathrm{\Omega }_t`$ and $`\mathrm{\Gamma }_t`$ the eigenfrequency and damping rate, respectively, which are determined by the cluster. The index $`t`$ indicates that due to ionization and expansion of the cluster, both, $`\mathrm{\Omega }_t`$ and $`\mathrm{\Gamma }_t`$, may depend parametrically on time. Under the assumption of a spherical, uniformly charged cluster with total ionic charge $`Q_{\mathrm{ion}}`$ and radius $`R`$ the potential inside the cluster is harmonic with an eigenfrequency $`\mathrm{\Omega }_t=\sqrt{Q_{\mathrm{ion}}(t)/R(t)^3}`$. The damping is caused by both, internal heating of the quasi-free electrons in the cloud and energy transfer to bound electrons. These two effects are responsible for outer and inner ionization, respectively.

For periodic driving $`F(t)=F_0\mathrm{cos}(\omega t)`$ the dynamics is given by $`X(t)=A_t\mathrm{cos}(\omega t\varphi _t)`$ with lali94

$`A_t`$ $`=`$ $`F_0/\sqrt{\left(\mathrm{\Omega }_{t}^{}{}_{}{}^{2}\omega ^2\right)^2+(2\mathrm{\Gamma }_t\omega )^2},`$ (4a)
$`\varphi _t`$ $`=`$ $`\mathrm{arctan}\left(2\mathrm{\Gamma }_t\omega /(\mathrm{\Omega }_{t}^{}{}_{}{}^{2}\omega ^2)\right).`$ (4b)
The energy balance of the dynamics (4) is characterized, on one hand, by energy loss $`E_{\mathrm{loss}}`$ due to the damping and, on the other hand, by energy gain $`E_{\mathrm{gain}}`$ from the external laser field. The cycle-averaged energy transfer rates read

$`\dot{E}`$ $`=`$ $`\dot{E}_{\mathrm{loss}}+\dot{E}_{\mathrm{gain}}`$ (5)
$`=`$ $`\mathrm{\Gamma }_tA_{t}^{}{}_{}{}^{2}\omega ^2+{\displaystyle \frac{1}{2}}F_0A_t\omega \mathrm{sin}\varphi _t.`$
Obviously and well known lali94 , maximum $`\dot{E}_{\mathrm{gain}}`$ or optimal heating requires $`\varphi _t=\pi /2`$, i. e. resonant behaviour $`\mathrm{\Omega }_t=\omega `$. As shown in Fig. 2d, the phase shift $`\varphi _t`$ changes in time from 0 to $`\pi `$, thereby passing the resonance $`\varphi _t=\pi /2`$. This change is directly connected with the increased ionization of the cluster, cf. Fig. 2a. In particular, the resonance time $`t_{\mathrm{res}}`$ coincides with the time $`t_2`$ of maximal ionization rate, cf. Eq. (2). This applies to other laser wavelengths as well (For $`\lambda =1170`$ nm $`t_2`$ is somewhat smaller due to the early laser switch-off at $`t=80`$ fs.):     $`\lambda _{\mathrm{laser}}`$ 520 nm 780 nm 1170 nm     $`t_2\pm \delta t_2`$ ($``$24$`\pm `$13) fs (13$`\pm `$24) fs (40$`\pm `$26) fs   $`t_{\mathrm{res}}`$ $``$22 fs 19 fs 72 fs

Passing through the resonance is not necessarily connected with large amplitude oscillations, if the damping strength $`\mathrm{\Gamma }_t`$ is comparable to or larger than the eigenfrequency $`\mathrm{\Omega }_t`$, see Fig. 3. Otherwise, the phase shift $`\varphi _t`$ at resonance is independent of the damping (Fig. 3).

In order to validate the applicability of the driven and damped harmonic oscillator model of Eq. (3) we use amplitude $`A_t`$ and phase shift $`\varphi _t`$ of the collective oscillation from our results of Fig. 2 to determine frequency and damping rate according to

$`\mathrm{\Omega }_{t}^{}{}_{}{}^{2}`$ $`=`$ $`\omega ^2+(F_0/A_t)\mathrm{cos}\varphi _t,`$ (6a)
$`\mathrm{\Gamma }_t`$ $`=`$ $`\left(F_0/(2A_t\omega )\right)\mathrm{sin}\varphi _t.`$ (6b)
Figure 4 shows these parameters as a function of time $`t`$ for the same system as in Fig. 2. The calculated eigenfrequency $`\mathrm{\Omega }_t`$ (gray circles in Fig. 4) closely resembles the frequency of a uniformly charged sphere (solid line). Moreover, both frequencies match the laser frequency $`\omega `$ at the same time $`t20`$ fs providing additional support for the collective oscillator model. Along with the decrease of the eigenfrequency $`\mathrm{\Omega }_t`$ the damping term $`\mathrm{\Gamma }_t`$ (white diamonds in Fig. 4) rises for times $`t60\mathrm{}30`$ fs. This accounts for energy transfer to deeper and deeper bound electrons, which does not occur in other theoretical studies where also resonant behaviour was discussed lajo99 ; paco+01 . Either inner ionization was not considered lajo99  or deuterium clusters composed of single electron atoms were discussed paco+01 . However, in order to understand the experimentally observed high charge states diti+97 ; ledo+98 ; zwdi+99 ; kosc+99  it is of utmost importance to take this into consideration. It is just this continuous cycle of effective heating and induced inner ionization at the resonance which causes the increased ionization rate and the high charge states of the fragments. At this point it is interesting to note, that the damping almost completely compensates the heating as can be deduced from the almost negligible increase of the amplitude of $`v_{\mathrm{CM}}`$ before reaching the resonance (Fig. 2c). Once the cluster has passed the resonance, however, the damping falls off rapidly ($`t30`$ fs in Fig. 4). Obviously, the oscillating electron cloud becomes unable to drive further inner ionization. This implies a weaker damping and thus an increase of the oscillation amplitude, cf. Fig. 2c.

In summary, we have shown that the electron emission in medium-sized rare gas clusters ($``$ 103 atoms) is enhanced by resonant energy absorption in agreement with experimental data zwdi+99 ; kosc+99 . (The metallic nature of the clusters used in kosc+99  should be of minor importance for the creation of the high charge states $`8`$ since the delocalized valence electrons are emitted early in the pulse.) Our microscopic calculations of the motion of ions and electrons using a hierarchical tree code reveal a laser-driven collective oscillation of the cloud of quasi-free electrons which are held back inside the cluster volume by the space charge of the cluster. The eigenfrequency of this oscillation is determined by charge and size of the cluster. Electron emission and cluster expansion change these quantities in time and eventually enable the matching of eigenfrequency $`\mathrm{\Omega }_t`$ and laser frequency $`\omega `$ during the pulse. This resonance allows for an effective energy transfer to the collective motion. The changing phase shift between driving field and driven electron cloud clearly indicates the different stages of energy absorption of the cluster electrons from the laser field. The fact that the collective electron dynamics can be well described by a simple damped harmonic oscillator helps to clarify the nature of this dynamics and provides a clear signature of this type of collective dynamics. We expect that the other two mechanisms for effective energy absorption, namely enhanced ionization and nanoplasma excitation will have an almost vanishing amplitude for the electronic cm velocity. Small clusters, which exhibit enhanced ionization, produce only a few quasi-free electrons which cannot create a sizable cm-velocity amplitude. In large clusters, however, where a nanoplasma is formed, electrons are heated resonantly at their critical density leaving their cm position at rest.