text:
# Ionization, Magneto-rotational and Gravitational Instabilities in Thin Accretion Disks Around Supermassive Black Holes

## 1. Introduction

Thin accretion disks around compact objects are locally unstable at central temperatures $`10^4`$ K, due to large opacity changes when hydrogen is partially ionized (Meyer & Meyer-Hofmeister 1981). It is generally believed that this instability drives the outbursts of dwarf novae (DN; Smak 1984; Cannizzo 1993 and references therein) and soft X-ray transients (SXTs; van Paradijs & Verbunt 1984; Mineshige & Wheeler 1989; Dubus et al. 2000). Physically, the disk cannot accrete at a steady rate because a range of radii are thermally and viscously unstable. As a result the disk alternates between two stable branches (see §2): (1) a hot, ionized branch corresponding to large accretion rates (outburst phase) and (2) a cold, neutral branch corresponding to low accretion rates (quiescent phase).

Global time-dependent calculations of the thermal limit cycle show that, in order to reproduce observed DN outbursts, viscosity on the cold, neutral, branch must be much smaller than on the hot, ionized, branch ($`\alpha _{\mathrm{cold}}0.01<<\alpha _{\mathrm{hot}}0.1`$; e.g., Smak 1984). For $`\alpha _{\mathrm{cold}}\alpha _{\mathrm{hot}}`$, the disk only undergoes small amplitude luminosity fluctuations, rather than the large amplitude outbursts observed. Gammie & Menou (1998; hereafter GM) showed that the magnetic Reynolds number on the cold branch in DN disks is $`Re_M3\times 10^3`$. Numerical simulations of the magneto-rotational instability (Balbus & Hawley 1991; 1998) show that, at such low magnetic Reynolds numbers, MHD turbulence and its associated angular momentum transport are significantly reduced (Hawley, Gammie, & Balbus 1996; Fleming, Stone, & Hawley 2000). Menou (2000) used Fleming et al’s calculation of angular momentum transport as a function of $`Re_M`$ in a time-dependent numerical model for DN outbursts. He found that, as proposed by GM, MHD turbulence dies away in quiescence. This interplay between the thermal ionization instability and the magneto-rotational instability provides a physical motivation for the observational inference that $`\alpha _{\mathrm{cold}}\alpha _{\mathrm{hot}}`$; it suggests that either accretion completely stops in quiescence or an additional, less efficient, angular momentum transport mechanism operates (perhaps spiral density waves excited by the companion; Spruit 1987; Menou 2000).

Thin accretion disks around supermassive black holes in active galactic nuclei are also expected to have regions with $`T_c10^4`$ K, and are thus also subject to the thermal ionization instability (e.g., Lin & Shields 1986; Cannizzo 1992; Siemiginowska, Czerny, & Kostyunin 1996). Applications of the disk instability model to AGN have, however, received less attention than those to DN and SXTs; this is no doubt because of the much longer recurrence times, which are observationally inaccessible in individual sources. The thermal ionization instability could, however, significantly influence the luminosity function of AGN, perhaps accounting for the short “duty cycle” of bright Seyferts and quasars (e.g., Mineshige & Shields 1990; Siemiginowska & Elvis 1997). In addition, it could lead to a radially varying accretion rate in the disk, complicating comparisons between theoretical models and observed broad band spectra. Understanding the role of the ionization instability and other sources of large amplitude variability in AGN is thus of considerable importance.

In this paper we generalize GM’s analysis and argue that, in contrast to SXTs and DN, the gas is well coupled to the magnetic field in AGN disks even on the cold branch of the thermal limit cycle (quiescence). This is essentially because the disks in AGN are physically much larger than those in DN; for a fixed resistivity, the magnetic Reynolds number is thus much larger. Our analysis implies that the thermal ionization instability is unlikely to be a significant source of large amplitude variability in AGN.

Accretion disks around supermassive black holes also differ from disks in DN and SXTs in that they may develop both local and global gravitational instabilities, which can lead to efficient angular momentum transport (e.g. Shlosman, Begelman & Frank 1990; Kumar 1999; Gammie 2000). In this paper, we argue that the degree to which the magneto-rotationally and gravitationally unstable portions of the accretion disk overlap is important for understanding time-dependent accretion in AGN.

We proceed as follows. We first calculate the structure of AGN disks in regions subject to the thermal ionization instability (§2). We then estimate the importance of resistivity, ambipolar diffusion, and Hall currents in these regions (§3). In §4 we reassess the role of the ionization instability in AGN. In §5 we discuss the interplay between MHD (un)stable111Throughout this paper, we use “MHD unstable” to mean unstable to developing MHD turbulence via the magneto-rotational instability. and gravitationally (un)stable regions of AGN accretion disks, and in §6 we conclude.

## 2. The Local Thermal Structure of AGN Disks

We consider a gaseous, thin accretion disk around a supermassive black hole (BH) of mass $`M`$. We ignore X-ray irradiation in calculating thermal equilibria (but see §5).222X-ray irradiation can lead to higher temperatures than considered below; this can substantially change the thermal ionization instability (e.g., Dubus et al. 1999). The local thermal equilibria of the disk are then found by equating the local viscous dissipation rate $`Q^+`$ to the radiative cooling rate $`Q^{}`$, where $`Q^{}`$ is a function of the disk opacity. Such equilibria have previously been calculated by several authors for AGN disks (see, e.g., Mineshige & Shields 1990; Cannizzo 1992; Siemiginowska et al. 1996). As for disks in DN and SXTs, these equilibria are thermally and viscously unstable in regions of partial ionization, corresponding to $`T_c10^4`$ K.

Figure 1 shows two thermal equilibrium curves (“S-curves”) for an annulus located at $`R=10^{14}`$ cm around a BH of mass $`M=10^7M_{}`$; for the upper S-curve, $`\alpha _{\mathrm{hot}}=\alpha _{\mathrm{cold}}=0.1`$, while for the lower S-curve, $`\alpha _{\mathrm{hot}}=0.1\alpha _{\mathrm{cold}}=0.01`$. The thermal equilibria have been found by calculating a grid of models for the disk vertical structure, for different values of the surface density $`\mathrm{\Sigma }`$ and the central temperature $`T_c`$, at the radius of interest (Hameury et al. 1998). Equilibria correspond to ($`\mathrm{\Sigma }`$,$`T_c`$) pairs such that $`Q^+=Q^{}`$.

The middle, inverted sections of the S-curves are where the gas in the disk is partially ionized and unstable. Any disk annulus ostensibly in this state instead undergoes a limit cycle in which it time-dependently evolves between the two stable (hotter and colder) thermal equilibria.

Point A in Figure 1 represents the lowest temperature reached by any unstable annulus during the thermal limit cycle. This point is of interest because it is the place in the limit cycle where the gas in the disk is least coupled to magnetic fields (it has the smallest ionization fraction in LTE). As a conservative estimate, we use $`T_c2000`$ K for the central temperature of the unstable annulus at point A. This is appropriate for $`\alpha _{\mathrm{cold}}=0.01`$.333This result can be checked by solving the vertically-averaged equation $`Q^+=(9/8)\nu \mathrm{\Sigma }\mathrm{\Omega }^2=Q^{}=4\sigma T_c^4/(3\mathrm{\Sigma }\kappa )`$ for an optically-thick disk, where $`\nu `$ is the kinematic viscosity and $`\kappa `$ is the opacity (which can be obtained from a standard opacity table for the range of density and temperature of interest). It is also interesting to note that the temperature $`T_A2000`$ K at point A on the S-curve is, for $`\alpha _{\mathrm{cold}}0.01`$, a generic prediction for any black hole mass and at any radius in the disk (see, e.g., GM for an identical result in the case of a DN disk). If, as we argue in §3-4, MHD turbulence is responsible for accretion in AGN disks in the vicinity of point A, $`\alpha _{\mathrm{cold}}\alpha _{\mathrm{hot}}0.1`$ would be more appropriate (see Hawley 2000 for global MHD simulations which give $`\alpha 0.1`$). This only strengthens our conclusions.

## 3. Coupling of the Gas to the Magnetic Field

GM showed that $`Re_M3\times 10^3`$ at point A in DN disks; here we generalize their analysis and assess the coupling of gas to magnetic fields in unstable annuli of accretion disks around compact objects of arbitrary mass.

The importance of ohmic diffusion can be estimated using the magnetic Reynolds number:

$$Re_M=\frac{c_sH}{\eta },$$
(1)
where $`Hc_s/\mathrm{\Omega }`$ is the scale height of the disk, $`c_s`$ is the sound speed $`\eta =c^2m_e\nu _{en}/(4\pi n_ee^2)`$ is the resistivity, $`\nu _{en}8.3\times 10^{10}T^{1/2}n_n`$ s-1 is the electron-neutral collision frequency (Draine, Roberge, & Dalgarno 1983), and $`n_n(n_e)`$ is the number density of neutrals (electrons). For the relatively low density disks of interest in AGN, ambipolar diffusion can be more important than resistive diffusion. Its importance can be estimated using the dimensionless number (e.g., Blaes & Balbus 1994)

$$Re_A=\frac{\nu _{ni}}{\mathrm{\Omega }}.$$
(2)
$`Re_A`$ measures the ratio of the neutral-ion collision frequency, $`\nu _{ni}`$, to the orbital frequency in the disk, $`\mathrm{\Omega }`$; for $`Re_A1`$, the neutrals are well coupled to the magnetic field by collisions with the ions. In equation (2), $`\nu _{ni}=\gamma \rho _i`$, $`\gamma 3.5\times 10^{13}`$ cm3 s-1 g-1 is the drag coefficient for ion-neutral collisions444We take $`m_i=30m_H`$ and $`m_n=2.33m_H`$, where $`m_i`$ is the ion mass, $`m_n`$ is the neutral mass, and $`m_H`$ is the mass of the hydrogen atom. (Draine et al. 1983), and $`\rho _i`$ is the ion mass density.

To estimate $`Re_M`$ and $`Re_A`$, we use the physical conditions in the limit cycle least favorable to MHD turbulence, namely those at point A; in particular, we use $`T_A2000`$ K, a somewhat conservative value (§2). From our numerical calculations of S-curves, we find that the surface density at point A can be approximately fit by (see also Cannizzo 1992; Siemiginowska et al. 1996)

$$\mathrm{\Sigma }_A\mathrm{\Sigma }_0\alpha _{0.1}^{3/4}rm_9^{2/3},$$
(3)
where $`\mathrm{\Sigma }_01000`$ g cm-2, $`\alpha _{0.1}=(\alpha /0.1)`$, $`r`$ is the radius of the unstable annulus in Schwarzschild units, and $`m_9=(M/10^9M_{})`$. Finally, for a given rate at which mass is fed to the accretion flow at large radii ($`\dot{M}`$, in $`M_{}\mathrm{yr}^1`$), the radius around which the ionization instability sets in is given by

$$r_Ar_0\dot{M}^{0.4}m_9^{2/3},$$
(4)
where $`r_0300`$. Equations (3) and (4) can be roughly reproduced analytically by fixing $`T10^4`$ K in a Shakura-Sunyaev disk (thus isolating the portions of the disk where hydrogen starts being partially ionized). It is also interesting to note that equations (3) and (4) are applicable to DN and SXTs as well as AGN, with $`m_9`$ then $`10^9`$, $`\dot{M}10^9M_{}\mathrm{yr}^1`$, and $`r_A10^5`$.

Using the Saha equation for a solar metallicity gas, we find that, in LTE, the ionization fraction at $`T_A2000`$ K in the thermally unstable annulus is $`x_en_e/n_n10^610^7`$ for $`\dot{M}ϵ[1,10^{10}]M_{}\mathrm{yr}^1`$.555The neutral density at point A in the unstable annulus is $`n_n3\times 10^{15}\dot{M}^{0.2}`$ cm-3. From the Saha equation, $`x_en_n^{1/2}\dot{M}^{0.1}`$. Thus, at fixed $`T_c`$, the $`10`$ order of magnitude change in $`\dot{M}`$ from DN to luminous AGN implies a change in $`x_e`$ by a factor of $`10`$. Larger values of $`x_e`$ are expected in high metallicity AGN disks because $`n_e`$ is dominated by contributions from metals with small ionization potentials. Substituting equations (3) and (4) into equations (1) and (2) we then find the smallest values of $`Re_M`$ and $`Re_A`$ in the thermal limit cycle (obtained at point A):

$$Re_M10^8\left(\frac{x_e}{10^7}\right)\left(\frac{r_0}{300}\right)^{3/2}\dot{M}^{0.6}$$
(5)
and

$$Re_A3\times 10^7\alpha _{0.1}^{3/4}\left(\frac{x_e}{10^7}\right)\left(\frac{r_0}{300}\right)\left(\frac{\mathrm{\Sigma }_0}{1000}\right)\dot{M}^{0.4}.$$
(6)
For a given accretion rate, equations (5) and (6) determine the coupling of the gas to the magnetic field in the annulus which is unstable to the ionization instability; this annulus lies at a radius given by equation (4). Note that $`Re_M`$ and $`Re_A`$ are independent of black hole mass at fixed $`\dot{M}`$.

Hall currents can significantly modify the dynamics of low density, weakly-ionized disks (e.g., Wardle 1999; Balbus & Terquem 2000). Their importance can be measured using the ion Hall parameter, $`\beta _i\omega _i/\gamma \rho _n`$, where $`\omega _i`$ is the ion cyclotron frequency. For $`\beta _i1`$ the ambipolar diffusion limit applies, while for $`\beta _i1`$ Hall currents are important. Assuming equipartition magnetic fields, we estimate $`\beta _i0.3\dot{M}^{0.1}(\mathrm{\Sigma }_0/1000)^{1/2}`$ in the unstable annuli of interest. The role of Hall currents in modifying MHD turbulence is determined by $`Re_H\beta _iRe_A`$; since $`\beta _i1`$, $`Re_A`$ is a reasonable proxy for $`Re_H`$: if $`Re_A1`$, then $`Re_H1`$, MHD turbulence is fully developed, and Hall currents are unimportant. By contrast, if $`Re_A1`$ MHD turbulence is modified (and Hall currents must be considered because of their potentially destabilizing effects).

## 4. Implications

MHD turbulence and its associated angular momentum transport are expected to be suppressed for $`Re_M`$ and/or $`Re_A`$ less than critical values $`Re_M^c`$ and $`Re_A^c`$. Numerical simulations with explicit resistivity are consistent with $`Re_M^c10^4`$ (Hawley et al. 1996; Fleming et al. 2000), while simulations of ion-neutral disks give $`Re_A^c100`$ (Hawley & Stone 1998; see also Mac Low et al. 1995; Brandenburg et al. 1995).

Equations (5) and (6) then show that there is an important difference between accretion in DN and SXTs, where $`\dot{M}10^9M_{}\mathrm{yr}^1`$, and accretion in luminous AGN, where $`\dot{M}0.0110M_{}\mathrm{yr}^1`$. As shown by GM, $`Re_M<Re_M^c`$ in DN and SXTs. MHD turbulence is thus suppressed on the “cold” branch of the thermal limit cycle.

By contrast, for reasonable $`\dot{M}`$ in AGN, $`Re_MRe_M^c`$ and $`Re_ARe_A^c`$. This implies that the disk is MHD-turbulent throughout the thermal limit cycle. In contrast to DN and SXTs, the efficiency of angular momentum transport should then be comparable on the hot and cold branches ($`\alpha _{\mathrm{cold}}\alpha _{\mathrm{hot}}`$). In this case global time-dependent calculations show that the thermal ionization instability is probably not a significant source of time-dependent accretion in AGN; it leads to luminosity “flickering,” rather than large amplitude outbursts analogous to DN and SXTs (e.g., Smak 1984; Mineshige & Shields 1990).

The analysis presented above does not, of course, preclude that accretion in AGN is highly time dependent; it does, however, suggest that the thermal ionization instability, as applied to DN and SXTs, is not relevant to AGN. Radiation-pressure induced thermal and viscous instabilities at small radii may be an important source of time dependence in AGN accretion disks (Lightman & Eardley 1974; Piran 1978). Another, less well explored, possibility lies in the interplay between MHD (un)stable and gravitationally (un)stable portions of the accretion disk. We briefly address several of the relevant issues below; much more theoretical work on these problems is needed.

## 5. Magneto-rotationally and Gravitationally Unstable Regions

### 5.1. Magneto-rotational Stability

An obvious shortcoming of the argument in §3 is that it is local. We examined annuli in thin accretion disks which are subject to the thermal ionization instability and found that they are always MHD turbulent so long as $`\dot{M}10^7M_{}\mathrm{yr}^1`$. At radii somewhat larger than the thermally unstable region, however, the disk is everywhere stable, but on the cold branch of the S-curve in Figure 1; if viscous dissipation is the only source of heating, the central temperature in the disk eventually drops to $`1000`$ K. In LTE the ionization fraction then becomes very small and MHD turbulence dies away.

X-rays and cosmic-rays may provide sufficient heating and/or nonthermal ionization to maintain magnetic coupling in the disk at large radii. Blaes & Balbus (1994) estimated $`Re_M1`$ and $`Re_A1`$ from cosmic-ray ionization in the parsec scale circumnuclear disk in our Galaxy. Simple estimates of X-ray ionization yield similar results.

As a concrete example, consider the $`0.10.25`$ pc masing disk around the $`M=3.6\times 10^7M_{}`$ black hole in NGC 4258 (Miyoshi et al. 1995). Neufeld & Maloney (1995; NM) proposed that the masing region is strongly irradiated by the central X-ray source, due to a significant warp in the disk. For a disk illuminated at an angle cos$`{}_{}{}^{1}\mu `$ to the normal direction, the incident X-ray flux on the disk is $`f(\mu L_X)/(4\pi R^2)`$, where $`L_X4\times 10^{40}`$ ergs s-1 for NGC 4258 (Makishima et al. 1994) and $`\mu `$ may be as large as $`0.2`$ (NM; Herrnstein 1997). The volume ionization rate in the bulk of the disk due to the incident X-ray flux is $`\xi \mathrm{min}(1,\tau )f/(HE_i)`$ cm-3 s-1, where $`E_i37`$ eV is the typical energy required to generate one secondary electron (Glassgold, Najita, & Igea 1997), $`\tau =\mu ^1\mathrm{\Sigma }/\mathrm{\Sigma }_c`$ is the shielding optical depth of the disk, and $`\mathrm{\Sigma }_c1`$ g cm-2 is the column density corresponding to order unity optical depth for hard X-rays $`>\mathrm{\hspace{0.33em}10}`$ keV (the factor $`\mathrm{min}(1,\tau )`$ is included to take into account the possibility of a disk which is optically-thin to X-rays). NM find $`n_n3\times 10^7`$ cm-3 in the masing region and so $`\mathrm{\Sigma }0.1`$ g cm-2. Balancing X-ray ionization against dissociative recombination at a rate $`8.7\times 10^6n_e^2T^{1/2}`$ cm-3 s-1 yields $`n_e300`$ cm-3 and thus $`x_e10^5`$. From equations (1) and (2), this corresponds to $`Re_M10^{11}`$ and $`Re_A10^3`$. Somewhat counterintuitively, the disk in NGC 4258 may be both masing and MHD-turbulent! Observational support for this conclusion is provided by the work of Wallin, Watson, & Wyld (1998, 1999); they showed that the maser spectral line features observed in NGC 4258 can be reproduced reasonably well by including a turbulent velocity field in a thin Keplerian disk model.

In order to account for the absence of masing at radii $`<\mathrm{\hspace{0.33em}0.1}`$ pc, NM propose that X-ray irradiation “shuts off” inside $`0.1`$ pc (because the warp decreases with radius, $`\mu 0`$). If correct, MHD turbulence probably ceases to be a viable angular momentum transport mechanism (cosmic-ray ionization alone yields $`Re_M10^7\xi _{17}^{1/2}`$, but $`Re_A0.1\xi _{17}^{1/2}`$, where $`\xi =10^{17}\xi _{17}`$ s-1 is the cosmic-ray interaction rate per hydrogen atom).666The ion Hall parameter is $`\beta _i10^3`$ assuming equipartition fields, so Hall currents should not change this result. Scattered X-rays from the vicinity of the AGN could, however, provide significant ionizing photons; their importance is very uncertain.

This illustrates a more general worry; the role of nonthermal heating and ionization is quite uncertain, depending on, e.g., the geometry of the disk (warps, flaring, etc.), the presence of X-ray scatterers, and the possibility of an anomalously large cosmic-ray flux due to a nearby jet. It is possible that many parsec scale disks in AGN are not magnetically well coupled. In this case gravitational instability is the most promising angular momentum transport mechanism.

### 5.2. Gravitational Stability

Thin disks are locally gravitationally unstable for $`Q<\mathrm{\hspace{0.33em}1}`$, where $`Q=c_s\mathrm{\Omega }/(\pi G\mathrm{\Sigma })`$. The nonlinear outcome of local gravitational instability depends on the ratio of the cooling time to the rotational period of the disk (Shlosman & Begelman 1989; Shlosman et al. 1990; Gammie 2000). For long cooling times, the disk reaches a steady state with efficient angular momentum transport produced by “gravito-turbulence” (Shlosman et al. 1990; Gammie 2000). For parsec-scale disks in AGN, however, the cooling time is generally much shorter than the rotational period of the disk; in this case the disk fragments, and probably undergoes efficient star formation (Shlosman & Begelman 1989; Gammie 2000). It is unlikely that local gravitational instability produces significant accretion in this regime.

For $`M_{\mathrm{disk}}2\pi \mathrm{\Sigma }R^2>M`$, however, the disk is globally gravitationally unstable. Using $`\dot{M}=2\pi \mathrm{\Sigma }R|v_r|`$ and $`|v_r|\alpha c_s^2/v_\varphi `$, where $`v_\varphi `$ is the (Keplerian) rotational velocity, this condition can be rewritten as (Shlosman & Begelman 1989)

$$\dot{M}>\frac{\alpha c_s^2v_\varphi }{G}10^2\alpha _{0.1}T_3v_{100}M_{}\mathrm{yr}^1,$$
(7)
where $`T=10^3T_3`$ K and $`v_\varphi =100v_{100}`$ km s-1.

If nonthermal heating and ionization are insufficient to maintain MHD coupling, disks will become MHD stable if $`T_3<\mathrm{\hspace{0.33em}1}`$. For luminous AGN with $`\dot{M}>\mathrm{\hspace{0.33em}10}^2M_{}\mathrm{yr}^1`$, however, the disk is globally gravitationally unstable in this region (eq. ). It is therefore likely that bar formation in the gaseous disk leads to significant angular momentum transport – and thus accretion – on a dynamical timescale (Shlosman et al. 1990). At smaller radii, where $`T_3>\mathrm{\hspace{0.33em}1}`$, the disk is both MHD and gravitationally unstable; the dynamics in this regime has been poorly explored but may be accessible with numerical simulations. More generally, in luminous AGN the accretion disk probably evolves from gravitationally unstable and MHD stable at large radii to gravitationally unstable, MHD unstable, and finally to gravitationally stable, MHD unstable at the smallest radii. It is rather uncertain whether these transitions occur at constant $`\dot{M}`$.

In low-luminosity AGN with $`\dot{M}<\mathrm{\hspace{0.33em}10}^2M_{}\mathrm{yr}^1`$ equation (7) shows that the disk is globally gravitationally stable where it is MHD stable; for $`\dot{M}<\mathrm{\hspace{0.33em}10}^3M_{}\mathrm{yr}^1`$ the disk has $`Q>\mathrm{\hspace{0.33em}1}`$ and is locally stable as well; this raises the interesting possibility of an MHD and gravitationally stable “dead zone” in which accretion may proceed very inefficiently (or perhaps not at all).777Accretion may proceed in an ionized surface layer analogous to Gammie’s (1996) proposal for T-Tauri disks. Mass will build up in this region until $`M_{\mathrm{disk}}>M`$, until the disk warp evolves so as to provide significant X-ray ionization, or until a global hydrodynamic instability sets in (analogous to GM’s proposal for DN). It is possible that such piling up of mass is presently occurring in many nearby galaxies. In particular, returning again to NGC 4258, NM’s masing model gives $`\dot{M}10^5\alpha _{0.1}M_{}\mathrm{yr}^1`$ and thus a very gravitationally stable disk: $`M_{\mathrm{disk}}100M_{}M`$ and $`Q1001`$ (see Gammie, Narayan, & Blandford 1999 for alternate $`\dot{M}`$ estimates). If X-ray irradiation and MHD turbulence indeed shut off for $`R<\mathrm{\hspace{0.33em}0.1}`$ pc, it is unlikely that significant accretion occurs at these radii. It is then also unlikely that NM’s accretion rate estimate is relevant for the energy-producing portions of the accretion flow at small radii.

## 6. Conclusion

We have studied the combined role of the thermal ionization, magneto-rotational and gravitational instabilities in thin accretion disks around supermassive black holes.

In contrast to dwarf novae and soft X-ray transients, we find that there is no physical motivation for a reduced efficiency of angular momentum transport when AGN disks become neutral via the thermal ionization instability. This is because even the predominantly neutral gas is well coupled to the magnetic field ($`Re_M1`$ and $`Re_A1`$). MHD turbulence due to the magneto-rotational instability should thus be present and similar to that in an ionized plasma. In this case, numerical simulations have shown that the thermal ionization instability leads to small amplitude luminosity “flickering,” rather than large amplitude outbursts analogous to DN and SXTs.

We propose, however, that accretion in AGN can be time-dependent because of the interplay between the magneto-rotational and gravitational instabilities, both of which are viable candidates for transporting angular momentum in AGN accretion disks. The hot inner regions of the disk are preferentially subject to the magneto-rotational instability, while the cold outer regions are prone to developing gravitational instabilities. In high luminosity systems (accretion rate $`\dot{M}>\mathrm{\hspace{0.33em}10}^2M_{}\mathrm{yr}^1`$), the magneto-rotationally and gravitationally unstable regions of the accretion disk overlap, while in lower luminosity systems ($`\dot{M}<\mathrm{\hspace{0.33em}10}^3M_{}\mathrm{yr}^1`$) they may not. In systems where these regions do not overlap, such as nearby low-luminosity galactic nuclei on scales of $`0.11`$ parsec, there is no well understood process for removing angular momentum from accreting gas. Mass may therefore be temporarily building up in a non-accreting “dead-zone.”

## Acknowledgments

We thank Bruce Draine, Pawan Kumar, Sergei Nayakshin, the referee, and especially Charles Gammie for useful comments or discussions. Charles Gammie provided a Mathematica routine to calculate ionization fractions. Support for this work was provided by NASA through Chandra Fellowship grant PF9-10008 (to EQ) and PF9-10006 (to KM) awarded by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-39073.


format: text

text:
# Calculation of minor hysteresis loops under metastable to stable transformations in vortex matter

## I Introduction

Supercooled or metastable states have been reported across first order phase transitions in vortex matter . It has also been established that isothermal field excursions cause the metastable phase to be converted to the stable phase . We have proposed  that the isothermal field variations provide a fluctuation energy that causes the metastable supercooled phase to cross the free-energy barrier and transform to the stable equilibrium phase.

The experimental techniques used to study the magnetic signatures of such transformations in vortex matter are (i) bulk dc measurements using a SQUID or a vibrating sample magnetometer which yields the magnetisation M of the entire sample ; (ii) bulk ac measurements of susceptibility which probe a region near the surface of the sample ; and (iii) local measurements of magnetic induction using magneto-optic or microhall probes that allow mapping the spatial profile B(x) . The last of these three techniques has been used recently to show that different (metastable and equilibrium) phases exist simultaneously in different regions of the sample . Similar inference has also been drawn by studies using the dc magnetisation  and ac susceptibility  techniques. Experiments have thus shown that metastable to stable transformations occur over local regions,and prompted by these developments we have given a formalism to calculate spatially-resolved energy dissipation under an isothermal field variation . In this paper we shall use this formalism to calculate how the metastable to stable transformation progresses inwards from the surface of the sample, under experimentally relevant field excursions. We shall also calculate the sample magnetisation M and the spatial field profile B(x) as the experimental observables.

## II Modelling the Peak Effect

The response of a hard superconductor to external magnetic fields is understood in terms of Bean’s critical state model (CSM) . Bean had assumed that the critical current density JC is independent of field. While detailed agreement with experiments has required introduction of various functional forms of JC(B) (see e.g. ref. ), much of the essential physics is captured even by assuming a field-independent JC.

Experiments have recently been addressing the region below and near the onset of a peak in JC(B) at B=B1, where a first order phase transition is seen in some superconductors. The occurrence of this ”peak-effect” has been known in various superconductors for a very long time, but attempts to have a CSM describing this JC(B) have been made only recently . While our detailed analytical model  could be used for the subsequent calculations, our focus here is to understand whether qualitatively new and anomalous signatures in recent experiments  can be arising from metastable to stable transformations in vortex matter. In this paper we shall use a simple Bean-like assumption for the two phases of vortex matter viz.

$`J_C(B)`$ $`=`$ $`J_1\text{ for }BB_1`$ (1)
$`=`$ $`J_2\text{ for }BB_1`$ (2)
We stress that, at a fixed temperature, our model has only two constant parameters viz. (J2/J1) and B1. We have assumed above that phase 1 (characterised by J1) is the stable phase for B(x)$`B_1`$, and phase 2 (characterised by J2) is the stable phase for B(x)$`B_1`$. We recognise that JC is not a thermodynamic quantity, but is a physical property that changes discontinuously across the phase transition. Phase 2 can exist for B(x) $`B_1`$ as a supercooled metastable phase, and we shall address this possibility in the next section. In this section we shall assume, however, that the free energy barrier surrounding the metastable phase drops very sharply as B(x)=B1 is crossed, and supercooling or superheating does not occur.

To obtain magnetisation-vs-field (or M-H) curves, we consider the sample to be in the form of an infinite slab in parallel field, as this geometry has the simplest algebra amongst the zero demagnetization factor cases of infinite cylinders in parallel field. We shall also continue with Bean’s simplifying assumption of HC1=0 followed usually in the CSM .

We follow standard procedures  to solve the CSM, and show in figure 1 the envelope M-H curves obtained with equation (1), with the parameters B1=1500 mTesla, J1R=4mTesla, and J2R=10 mTesla. Here the slab has surfaces at x=$`\pm `$R, is infinite along the y and z directions, and the magnetic field is applied along the z-axis. (We shall consider only positive values of x in this paper; there is a symmetry about x=0.) Note that the first order transition shows different widths in M vs H when measured along the field-increasing and along the field-decreasing directions. This is because the shielding current density at x is dictated  by the local magnetic induction B(x) through equation (1), and B(x) is different from the applied field as well as different in the field-increasing and field-decreasing cases. In figure 2 we plot B(x) for some values of applied field H corresponding to the field-increasing and field-decreasing cases. We note that phase 1 and 2 exist simultaneously in two different regions of the sample. We emphasize that there is no metastability because the stable phase 1 exists wherever B(x)$`B_1`$ and the stable phase 2 exists wherever B(x)$`B_1`$.

The calculation above is for some fixed temperature T1, and we note that the phase transition field B1 falls as the temperature T1 rises .

## III Supercooling and metastable-to-stable transformations

We now consider that we have applied a field H1 which is smaller than B1(T1). But we apply this field at a much higher temperature T2 such that H1 is much larger than B1(T2). So, B(x) throughout the sample is larger than B1 at that temperature, and the entire sample is in phase 2. We assume further that B(x) is constant at T2. This happens if the critical current density J2 in phase 2 vanishes at T2. One can, however, also achieve a constant B(x) by applying an external field H1+hcos(wt), with (H1-B1(T2)) $`>`$ h $`>`$ J2(T2)R, and then slowly reducing the amplitude h to zero .

We now lower the sample temperature (i.e. field-cool) to T1 such that the sample is supercooled and is metastable in phase 2. As discussed in references , there is a free energy barrier fB(T) that keeps phase 2 metastable, where fB(T) is determined uniquely by B(x) and T. The vortex matter in the neighbourhood of x will transform to phase 1 when the fluctuation energy Pd(x) created by an isothermal field variation is larger than \[fB(T) - kT\].

The field profile B(x) in the field-cooled sample is constant at H1 (see figure 3), and we now start lowering the applied field H with the temperature fixed at T1. Since the sample is in the supercooled phase 2, the shielding currents set up initially will have a magnitude J2. The variations in field will cause a fluctuation energy Pd(x) given by equations (4) and (5) of reference , and the vortex matter in the neighbourhood of x=x0 will transform to the stable phase at Pd(x0)= fB(T) - kT =P0. This transformation is triggered from the surface  and the shielding current magnitude will drop to J1 for x $`>`$ x0. The point x0 moves from (x/R)=1 to (x/R)=0 as the applied field is lowered, and B(x) are shown in figure 3 for representative values of the applied field. We have used H1=1480 mTesla, and P0=18 (mTesla)2. The large (small) slopes of B(x) correspond to large (small) magnitudes of the shielding current density, and thus to vortex matter being in phase 2 (phase 1). From these B(x) one can readily calculate  the sample magnetisation as the field H is lowered. In figure 4 we show the minor hysteresis loop (MHL) obtained as the applied field is lowered after field-cooling. We show also the field-decreasing envelope curve from figure 1. Note that the MHL first shoots out above the envelope curve,and then slowly merges from above. This nature is in qualitative agreement with published data . If we had used the detailed JC(B) of reference  to model the peak-effect, instead of the simple model of equation (1), the peak of the MHL would be less sharp and the merger with the envelope curve would be slower. The simple model used to obtain figure 4 brings out the qualitative behaviour observed and captures the essential underlying origin of anomalous MHLs as being due to phase 2 being supercooled and the transformation from the metastable phase 2 to the stable phase 1 occurring progressively deeper into the sample.

We have assumed that P0 (and thus fB(T1)) is constant at 18 (mTesla)2. fB(T1) is actually dictated by B(x0), and falls monotonically as B(x0) moves farther from the phase transition line B1(T). As is seen in figure 3, B(x0) varies only by less than a few mTesla as the MHL merges with the envelope curve. The assumption of a constant P0 over an MHL is thus justified. If, however, we field-cool to the same temperature T1 at a lower field H2, then fB(T) will be lower . This implies that P0(HFC=H2) will be smaller than P0(HFC=H1). We show, in figure 5(a), the MHL for the case when the sample was field-cooled to H2=1400 mTesla where P0 is taken to be 4.5 (mTesla)2. In figure 5(b) we have taken HFC=H3=1300 mTesla, where fB must be still lower and is taken as P0=2 (mTesla)2. The MHLs again shoot out of the envelope curve, but to peak values progressively smaller than in figure 4. The merger of the MHLs with the envelope curve also occurs over a progressively narrower range of field reduction than in figure 4. This qualitative change in the nature of the MHLs with reduction of HFC is also consistent with published data .

## IV Conclusion

We have used the ideas developed in references  to calculate the isothermal field-cooled MHLs, and the spatial field profiles B(x). The model calculation was done, in the spirit of Bean’s original work , with field-independent critical current densities. The only parameters were (J2/J1), and the onset field B1 at which the peak effect starts in the field increasing case. We used the fact  that fB becomes smaller as B(x) falls below B1, and that metastable to stable transformations occur in local regions of the sample .

The formalism of reference  can similarly be used to calculate MHLs after different thermomagnetic histories. We assert here that the simple model of equation (1) reproduces qualitative features of various observations  of anomalous MHLs. As was stated in the Introduction, more detailed tests of the extent of phase coexistence are possible and calculated B(x) can be compared with field profiles measured with local probes. Our model also predicts the spatial region over which the two phases coexist, and the evolution of these regions under isothermal field variation.

We thank Drs. S.B. Roy and Sujeet Chaudhary for many helpful discussions.


text:
# Optical Variability of the T Tauri Star HH30Based on observations obtained with the 1.2 m Telescope at the F. L. Whipple Observatory

## 1. Introduction

Signatures of hot and cool starspots have been observed in numerous T Tauri stars, primarily in the form of photometric variability, with periodic variability detected in hundreds of systems (e.g., Bouvier et al. 1993; Wichmann et al. 1998; Choi & Herbst 1996; Makidon et al. 1997; Stassun et al. 1999). While cool spots have been mapped through Doppler imaging of weak-lined T Tauri stars (e.g., Hatzes 1995), hot spots tend to be associated with variability in Classical T Tauri stars (e.g., Kenyon et al. 1994; Bouvier et al. 1993; Herbst et al. 1994). The hot spots are thought to be associated with the accretion process, and the currently favored magnetic accretion models naturally provide for hot spots on stellar surfaces (Ghosh & Lamb 1979a, 1979b; Königl 1991; Shu et al. 1994; Ostriker & Shu 1995; Najita 1995).

Hot spots yield specific brightness and polarization variations as they rotate into and out of view (Wood et al. 1996; Mahdavi & Kenyon 1998; Stassun & Wood 1999). Wood & Whitney (1998) investigated the effects of hot starspots on the morphology of scattered light disks. The non-axisymmetric illumination of the disk by hot spots leads to an asymmetric brightening of the disk. Such asymmetric brightening has been detected in HST observations of the edge-on disk system HH30 (Burrows et al. 1996; Stapelfeldt et al. 1999a), hinting at possible agreement with the magnetic accretion model.

The gross consistency between the existing observations of HH30 and the predictions of the magnetic accretion model is tantalizing, but ambiguous due to the large timescales separating the HST observations. Periodic variability due to magnetic accretion is expected to match the rotational periods for Classical T Tauri stars: of the order less than one day to tens of days (Attridge & Herbst 1992; Bouvier et al. 1993; Edwards et al. 1993; Eaton et al. 1995; Stassun et al. 1999; Herbst, et al. 2000). In some cases the hotspots are stable over many rotation periods.

Motivated by the apparent success of the magnetic accretion model at explaining the morphological variations observed in the limited HST imagery, we have undertaken a ground-based photometric monitoring campaign on HH30 to determine the timescale of any variability. In Section 2 we present the results of our observational campaign, in Section 3 we present models for the observed photometric variability, and we summarize our findings in Section 4.

## 2. Observations

Data for HH30 (RA = 04 31 37.5, DEC = 18 12 26.0; $`V19.5`$, $`I17.3`$, Mundt & Fried 1983) were obtained with the 1.2 m telescope at the F. L. Whipple Observatory (FLWO), using the “4Shooter” CCD mosaic (Szentgyorgyi et al. 2000) with four thinned, back-side illuminated, AR coated Loral $`2048^2`$ CCDs. The camera has a pixel scale of 0.335 arcsec/pixel and field of view of roughly 11.5 arcmin for each chip. The data were taken in the $`2\times 2`$ CCD binning mode. To obtain multi-wavelength variability information we used $`VRI`$ filters from the FLWO “Harris Set”. We obtained data (see Figure 1), between September 7th 1999 and February 28th 2000, usually consisting of up to three sets of 300 sec, 180 sec, 180 sec $`VRI`$ exposures. In total we have 24 useful images of HH30 in the $`V`$-band and 27 useful images in each of the $`RI`$-bands, spanning about 174 days. The seeing ranged from $`FWHM=1.3`$ arcsec to about $`4.0`$ arcsec. The $`VRI`$ data can be accessed from anonymous ftp at ftp://cfa-ftp.harvard.edu/pub/kstanek/HH30.

Preliminary processing of the CCD frames was done with the standard routines in the IRAF-CCDPROC package.111IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Associations of Universities for Research in Astronomy, Inc., under cooperative agreement with the NSF Photometric variability was determined using two separate techniques: aperture photometry and point-spread function (PSF) fitting. PSF fitting was done using the photometric pipeline of the project DIRECT (Kaluzny et al. 1998; Stanek et al. 1998), based on the Daophot/Allstar photometry package (Stetson 1987, 1991). The differential light curves obtained by both methods agreed well, which was reassuring considering that the HH30 nebulosity is somewhat resolved in our images. As a further check, we saw no correlation between the derived photometric changes and the seeing or sky brightness.

The resulting $`VRI`$ differential light curves are shown in Figure 1 and are adjusted so the faintest measurement in each band corresponds to $`\mathrm{\Delta }`$mag$`=0`$. The HH30 nebulosity was observed to vary in brightness over timescales of a few days by $`\mathrm{\Delta }V1.4`$mag, $`\mathrm{\Delta }R1.0`$mag, and $`\mathrm{\Delta }I1.1`$mag. The changes in brightness are very well correlated between the three photometric bands. There are hints of periodicity with the rise and fall at the start of our observations and again around HJD 1570 in Figure 1. Independently for each band we searched for periodicity using a variant of the Lafler-Kinman (1965) string-length technique, proposed by Stetson (1996) and also the Lomb (1976) normalized periodogram method. The strongest signals in our periodogram techniques were found at 19.8 days (97% probability) and 11.6 days (95% probability). An 11.6 day period is more typical of T Tauri stars and further monitoring of HH30 will possibly yield a more reliable period. We emphasize that this period detection is very tentative given the present data, and we encourage follow-up observations to confirm or refute our finding.

Photometric variability in excess of 2 mags may be obtained by hot spots on the stellar surface (e.g., Mahdavi & Kenyon 1998). Hot spots produce larger photometric variability at shorter wavelengths due to the increased spot/star luminosity ratio. Our observations show that $`\mathrm{\Delta }V>\mathrm{\Delta }I`$, but $`\mathrm{\Delta }R<\mathrm{\Delta }I`$. The apparent inconsistency of the $`R`$ band variability may be attributed to contamination from strong H$`\alpha `$ and \[S ii\] emission associated with HH30’s jet (Burrows et al. 1996). Bouvier et al. (1999) discounted a hotspot model for AA Tau because they did not detect any color variations ($`\mathrm{\Delta }`$ mag$`1.4`$ in all $`BVRI`$ bands). Our HH30 data do show color variations and in the next section we investigate hot spot models for HH30 and restrict our simulations to the $`V`$ and $`I`$ bands due to the likely contamination of the $`R`$ band observations by jet emission.

## 3. Models

The large variability we observe may occur for a model with two diametrically opposed hot spots when the second spot is occulted (either by the star or the circumstellar disk) and the observed spot is seen directly at maximum light and is obscured by the star at minimum light. However, for the high inclination of HH30 ($`i>80^{}`$, Burrows et al. 1996) we only see a scattered light nebula and do not see the star directly. For a model with two diametrically opposed hot spots, we found that at high inclinations the photometric variability was only around 0.2 mag (Wood & Whitney 1998). As one spot rotates behind the star, the second spot rotates to the front yielding the small amplitude variability. In addition, because we are only detecting scattered light we always see some reflected light from the hot spots yielding small variations in the total intensity.

The fact that we are viewing scattered light makes it impossible for us to apply inversion techniques to determine spot parameters (e.g., Bouvier et al. 1993; Vrba, Herbst, & Booth 1988), since the light curves are sensitive not only to the star and spot parameters, but also the disk structure and circumstellar dust properties. The new WFPC2 observations presented by Stapelfeldt et al. (1999a) and the large variability present in Figure 1 suggests that one hot spot dominates the nebular and photometric variability. We now construct single hot spot models for HH30.

HH30 has been modeled by Burrows et al. (1996) and Wood et al. (1998) with a flared disk geometry,

$$\rho =\rho _0\mathrm{exp}\frac{1}{2}[z/h(\varpi )]^2/\varpi ^\alpha ,$$
(1)
where $`\varpi `$ is the radial coordinate in the disk midplane and the scale height increases with radius, $`h=h_0\left(\varpi /R_{}\right)^\beta `$. Following Burrows et al. (1996) and Wood et al. (1998) we adopt: $`\beta =9/8`$, $`\alpha =15/8`$, $`h_0=0.05R_{}`$ (giving $`h[100\mathrm{A}\mathrm{U}]=15`$AU), and a disk mass of $`2.5\times 10^4M_{}`$. We construct scattered light models using our Monte Carlo radiation transfer code (Whitney & Hartmann 1992; Wood & Whitney 1998; Stassun & Wood 1999) and assume there is one circular hot spot of radius $`\theta _s`$ at latitude $`\varphi _s`$ on the stellar surface. The relative number of photons released from the spot and star is,

$$\frac{N_s}{N_{}}=\frac{1\mathrm{cos}\theta _s}{1+\mathrm{cos}\theta _s}\frac{B_\lambda (T_s)}{B_\lambda (T_{})},$$
(2)
where we assume that the luminosities of the star and spot are Planck functions at the spot ($`T_s`$) and stellar ($`T_{}`$) temperatures.

Figure 2 shows the simulated $`V`$ and $`I`$ band photometric variability for $`\theta _s=20^{}`$, $`\varphi _s=65^{}`$, $`T_s=10^4`$K, $`T_{}=3800`$K (Wood & Whitney 1998; Stapelfeldt et al. 1999). We assumed that the circumstellar dust can be characterized by a Kim, Martin, & Hendry (1994, hereafter KMH) size distribution, typical of grains in the interstellar medium. For KMH grains the $`V`$ and $`I`$ band dust parameters (opacity, albedo, peak polarization, and asymmetry parameter in the Heyney-Greenstein scattering phase function) are displayed in Table 1.

Viewing this model at low inclinations we can easily obtain $`\mathrm{\Delta }V1.5`$ mag. However at $`i=82^{}`$, we cannot reproduce the observed large amplitude variability — we only see scattered light and always detect reflected light from the hot spot, irrespective of rotational phase. Therefore the amplitude of variability is smaller ($`\mathrm{\Delta }V0.5`$ mag, $`\mathrm{\Delta }I0.3`$ mag) than when we see the star directly. One way of increasing the variability is to change the spot size and temperature. Alternatively, if the forward beaming of the scattered light is increased (increasing the asymmetry parameter, $`g`$), we will see less scattered light from the spot at minimum light when it is occulted by the star. Such an increase in the asymmetry parameter, indicative of large grains in the disk, was found by Burrows et al. (1996). Other modeling efforts have also found evidence for larger, more forward throwing grains in protostellar environments (e.g., Lucas & Roche 1998). We have therefore constructed models with larger, more forward throwing grains (see parameters in Table 1). This population of grains was used by Cotera et al. (2000) and Whitney & Wolff (2000) to model multiwavelength HST images of HH30. The dust grains are homogeneous spheres composed of either amorphous carbon (BE1, Rouleau & Martin 1991) or revised astronomical silicate (Weingartner & Draine 2000). The size distribution for each component is specified using a power-law with exponential decay (i.e., $`a^pexp[a/a_c]`$) with a cross-section weighted average radius (summed over both components) of 0.092 µm. The relative numbers of each grain type are such that the dust completely consumes slightly “super-solar” abundances in carbon and silicate: 400 ppm for C/H and 40 ppm for Si/H. We refer the reader to Whitney & Wolff (2000) for a more detailed discussion of the dust model.

With the larger grains the amplitude of the variability is unaffected at low inclinations because the direct starlight is much brighter than the scattered light. For $`i=82^{}`$, we found larger variability ($`\mathrm{\Delta }V1`$mag, $`\mathrm{\Delta }I0.5`$mag) than for the simulation using KMH grains due to the increased beaming of the scattered light by the forward throwing grains.

Increasing the phase function asymmetry through larger grains increases the photometric variability, but it is still less than observed. Therefore, we have investigated other models which have larger spot/star luminosity ratios and different spot parameters. Figure 3 shows a simulation in which we obtain $`\mathrm{\Delta }V1.45`$mags and $`\mathrm{\Delta }I1.05`$mags, in agreement with our observations. This simulation has $`T_s=10^4`$K, $`T_{}=3000`$K, $`\theta _s=20^{}`$, and $`\varphi _s=60^{}`$. We can obtain the same amplitude of variability with a model that has $`T_{}=3800`$K and $`T_s=2\times 10^4`$K. As the central star in the HH30 system is not observed directly, its effective temperature is somewhat uncertain and may be cooler than the $`T_{}=3800`$K, M0 spectral type determined by Kenyon et al. (1998). The large spot/star luminosity, required to match our observations, could result in considerable veiling allowing for a later spectral type (see Kenyon et al. 1998). The spot size and latitude may be constrained through detailed modeling of time series images of the scattered light disk (Stapelfeldt, et al. 1999a).

## 4. Summary

We have presented $`VRI`$ photometric observations of HH30 that show it is highly variable over timescales of a few days ($`\mathrm{\Delta }V1.4`$ mags, $`\mathrm{\Delta }R1.0`$ mags, $`\mathrm{\Delta }I1.1`$ mags), with hints of periodicity at 11.6 or 19.8 days. This contrasts with the recent work of Stapelfeldt et al. (1999b) whose HST observations suggest a characteristic timescale on the order of 15 years, possibly related to inhomogeneities in the disk. The much shorter timescale variability in our data indicates a stellar origin for the photometric variability we observe. High resolution imaging is required to determine whether the variability we observe is related to the morphological variations in the scattered light images.

Comparing our findings to classical T Tauri stars, HH30’s variability is among the largest yet reported. We have modeled the variability in the context of the magnetospheric accretion model with a single hot spot on the stellar surface. The large variability requires a large spot/star luminosity ratio and that the circumstellar grains are larger and more forward throwing than interstellar grains. Variability due to accretion hot spots is often stochastic (e.g., Herbst et al. 1994) and while our modeling has adopted a single stable hotspot, further observations are required to better sample the lightcurve and determine changes in the spot sizes and temperatures. Our modeling of the circumstellar dust is consistent with many other studies that are providing evidence for grain growth in dense environments of T Tauri stars.

While our observations and modeling support the magnetic accretion hypothesis for HH30, alternative sources of the variability may be orbiting dust clouds, inhomogeneities in the disk structure close to the star, a binary star, and stellar flares. Further monitoring of HH30 with sampling on the order of days will yield further insight into this system. If the variability is due to magnetic accretion there will be specific time sequence variations of the phometry, polarimetry, and morphology of the scattered light disk due to illumination by hot spots. Even if the accretion is stochastic there will be specific correlations between the photometry and polarimetry (Stassun & Wood 1999, Fig. 7). Rotational modulation will manifest itself in high resolution imaging studies through asymmetric brightening of scattered light images (Wood & Whitney 1998; Stapelfeldt et al. 1999a). Lower resolution imaging should still be able to detect the asymmetric brightening via shifting of the photometric centroid of the system as pointed out by Stapelfeldt et al. (1999a).

###### Acknowledgements.

We thank Scott Kenyon for discussions relating to this project and the following people who helped us in obtaining observations of HH30 at FLWO: Pauline Barmby, Ann Bragg, Paul Green, Saurabh Jha, Amy Mossman, Jose Munoz, and Rudy Schild. The referee, Bill Herbst, provided comments and suggestions that have improved this paper. We acknowledge financial support from NASA’s Long Term Space Astrophysics Research Program, NAG5 6039 (KW), NAG5 8412 (BW); the National Science Foundation, AST 9909966 (BW and KW); the HST Archival Research Program AR-08367.01-97A, (BW, KW, KS); and NASA contract NAS8-39073 (SJW). Support for KZS was provided by NASA through Hubble Fellowship grant HF-01124.01-99A from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.


text:
# Visualization of Large Multi-Dimensional Datasets

## 1. Introduction

How can huge data sources (gigabytes up to terabytes) be quickly and easily analyzed? There is no off-the-shelf technology for this. There are devastating computational and statistical difficulties; manual analysis of such data sources is now passing from being simply tedious into a new, fundamentally impossible realm where the data sources are just too large to be assimilated by humans. The only alternative is to provide extensive computer support for the process of discovery.

The focus of this workshop is the challenge posed by the next generation of large astronomical sky surveys. Specifically, we concentrate on the Sloan Digital Sky Survey (SDSS), which will create over one terabyte of reduced data over the next 5 years: How does one navigate such a huge, multi-dimensional, dataset? The techniques of information visualization and visual discovery may be extended to such datasets, allowing the scientist to interactively explore and understand her results in real time.

This work is the product of a collaboration at Carnegie Mellon University and the University of Pittsburgh, involving astronomers, experts in traditional scientific visualization and interactive information visualization, and computer and computational scientists. Condensed data representations from the data mining and machine learning community (Nichol et al. in these proceedings) make it possible to explore these huge datasets interactively.

## 2. Visualization Challenges of the Virtual Observatory

Over the next ten years, we will witness a revolution in how astrophysical research is performed. This is primarily due to the large number of new sky surveys presently underway (or completed) that are designed to map the Universe to higher sensitivity and resolution than ever previously envisaged. We are quickly approaching the prospect of a Virtual Observatory, where one can digitally reconstruct the whole sky. These surveys, and the virtual observatory, present scientists with a “gold mine” that the next generation of astrophysicists will spend their whole careers exploring.

Two cornerstones of the Virtual Observatory are the 2 Micron All Sky Survey (2MASS) and the optical Sloan Digital Sky Survey (SDSS). The 2MASS survey (which is 91released) is a near infrared imaging survey covering the full sky in three passbands (from 1- 2.2 microns). The SDSS is an imaging and spectroscopic survey that will cover one quarter of the sky at five different wavelengths. Together these two surveys will detect over 200 million objects (galaxies and stars) and from these detections positions, fluxes, shapes, textures and bitmaps will be extracted. In addition to the scientific information will be bookkeeping information that describes the observations themselves, e.g. whether the sky was cloudy or there were problems with the instrumentation. We must also understand these possible systematic uncertainties present within the data. The total 2MASS and SDSS surveys are expected to acquire 500 GB of cataloged attributes and 1 TB of postage stamp images (cutout images around each detected object) over the next 5 years of operation.

For each object detected hundreds of attributes will be recorded. The size and large dimensionality of these new data sets means that simple visualization and analysis techniques cannot be applied directly, because they do not scale effectively. The questions then become: How do we quickly determine the important dimensions within such a data set? Which dimensions tell us about how galaxies or stars form and how matter is distributed throughout the individual galaxies or, the Universe as a whole? New techniques developed for 2MASS and SDSS will be applicable to the analyses of all large observational data sets (such as the all-sky surveys of GALEX, ROSAT, MAP and PLANCK) and for the visualization of other large physical and biological experiments.

## 3. The Breakdown of Visual Discovery

Traditional astronomy has relied on the study of small numbers of rare objects. With hundreds of millions of objects in the Virtual Observatory, “rare” objects will themselves number in the millions. Statistical methods become necessary to group objects into classes for study.

This situation is opposite to that for which traditional visualization methods are best suited. Large dataset size alone kills interactivity. The record-oriented nature of the dataset makes it difficult to assimilate, because the brain is evolved to understand continuous systems that exist in three dimensions plus time. If too broad a view of the data is taken, important details can be literally too small to see. Too narrow a view is also disastrous. It becomes very easy to be distracted by structures that look important.

Supercomputer simulations routinely produce very large datasets. Standard methods of visualizing such data typically include the generation of one or more animations (assuming the system evolves in time), or interactive visualization of tiny subsets of the data. In most problem domains the researcher’s physical intuition applies, and can help her to assimilate the evolution of the system.

Large computing facilities are developing interactive visualization tools for these terabyte-sized datasets. Direct application of traditional visualization typically requires a supercomputer, and even with one available interactivity is very difficult to attain. Current rendering hardware can draw at most about 10 million polygons per second (though this number is rising), so interactive drawing of hundreds of millions of objects is simply impossible.

## 4. Desirable Capabilities

Here are illustrative questions that astrophysicists will want to ask of the data:

* A range of counting queries such as:

+ How many elliptical galaxies, with a redshift above 0.3, are there?
    + What is the mean and variance of ellipticity among radio galaxies within clusters of galaxies compared to outside clusters of galaxies?
    + How does the distribution of galaxy colors observed on 1 May 2001 compare to those seen on 1 June 2002?
* Sophisticated statistical queries that require the clustering, classification, regression, filtering and newer probabilistic inference techniques.
* Visualization requests to answer questions such as:

+ Give me a smoothed map of all X-ray detected galaxies.
    + Display all emission-line galaxies in detected clusters of galaxies; are they in the center or on the outskirts?

## 5. XGobi, a Simple Tool for Visualizing Multidimensional Data

XGobi, one of the most widely used information visualization tools, was written at Lucent Technologies 111XGobi is freely available from http://www.research.att.com/areas/stat/xgobi/. (Swayne, Cook, & Buja 1998). This tool interactively produces arbitrary 3D projections of n-dimensional data. For example, figure 1 shows the correspondence between data groups which differ in one projection but are cluster together in another. XGobi deals only with scalar data, and is designed for relatively small datasets.

XGobi supports a method called brushing, whereby the user selects data items in one window and the selection is propagated to the corresponding items in all other windows. This has proven to be a useful technique for mentally integrating multiple views.

## 6. Visage, an Information Exploration Tool

At Carnegie Mellon University, in collaboration with Maya Design Group, we have already developed an interactive data exploration system called Visage (Kolojejchick, Roth, & Lucas 1997). It is effective for analyzing high dimensional data, but it uses only 2D visualizations and is limited to small discrete datasets. Below, we describe Visage’s features, while in the next subsection, we discuss our plans to generalize these features and expand the capacity of Visage to handle massive datasets. We call this new system TeraVisage.

Visage presents data as graphemes (visual elements) organized by frames (lightweight nestable windows that impose some visualization discipline on the graphemes presented within them). Within a frame, the data are presented by graphemes such as bars, text labels, marks and gauges. The visual properties of the graphemes (e.g., the color or size of a plot point) encode attributes of the objects they represent (e.g., a galaxy). While each grapheme stands for one object, an object may be represented by many graphemes of varying appearance in different frames.

Some of the basic operations provided by Visage are:

* Drag-and-drop objects from one frame to another,
* Navigate (drill-down) from an object along a relation creates graphemes for the related objects,
* Aggregate (roll-up) a set of selected objects to create a new object with properties computed from its members,
* Brush an object or a set of objects in a choice of colors, as described in Section 5.,
* Dynamically query by an attribute to render invisible all objects within a frame whose values for that attribute fall outside a range selected by a slider widget.

The results of these operations depend only on the underlying data object, and are uniformly applicable in any type of Visage frame.

### 6.1. Going to Large Datasets

In this subsection, we discuss the methods by which we propose to adapt Visage to very large datasets. Interactive visualization requires interactive speeds; the challenge is to maintain these speeds during interrogation of a large (500GB) database. Our plan is to maintain a hierarchy of representations and subsets, with a hierarchy of access speeds.

Figure 2 presents our vision of how this would play out for exploring an astronomical dataset. The data comes from a simulation of the coming merger of our Milky Way galaxy with the Andromeda galaxy222Welling carried out the simulation using conditions specified by John Dubinski and code by Lars Hernquist and others \[Mihos, 1998 #176\]; further information is available from http://www.cita.utoronto.ca/ dubinski/index.html. An animation of the merger can be found at http://www.psc.edu/ welling/big-merger.mpg..

The dataset we examine here corresponds to the last frame of the simulation. In the course of the collision large “tidal tails” of stars, dust, and gas have been drawn from both galaxies. At the moment in question the bulk of the galactic matter has formed a roughly elliptical collision remnant, but two tidal tails remain. The matter in one tail is to be selected and examined.

In the left frame labeled “All Star Groups”, the astronomer has displayed 7 attributes of all the star groups in the two galaxies. The 3D visualization shows spatial position (x, y, and z). The histograms additionally show the distribution of distance from the center of the galaxy, magnitude of velocity, x component of velocity, and original galaxy each star group belonged to. On the galaxy histogram the astronomer has brushed the Milky Way star groups purple and the Andromeda ones yellow. The relative distributions can be seen in the other histograms and the density plot. It is apparent that star groups that originated in the Milky Way dominate the streamers. The astronomer selects a group of stars from one of the streamers using a bounding box. These stars are brightened in the 3D density plot. The brushing and selection operations should take about one second.

The astronomer then creates a new frame with a density plot showing the x, y, and z components of the velocities and drags (copies) the selected star groups there. She labels the new frame “+Z Tail” (right), and copies the histograms as well. Appropriate condensed representations should allow this operation to be carried out in about 10 seconds. Being interested in the relationship between the distance from the center and the magnitude of the velocity, she creates a 2D scatter plot of those variables. By moving the distance slider, she controls which stars remain fully visible and which are grayed out. Feedback from slider changes should occur in roughly 100ms.

It is easy to see how to extend most of the graphemes of Visage to large datasets. Sets of individual marks become density distributions in some space, and are selected by selecting regions in the space and/or using DQ filters. Gauges remain gauges, but now represent statistical summaries of subsets of the data. Text labels remain, but become labels for sub-aggregates found in traversing kd-trees or in more sophisticated models. Many confusing issues remain, however.

## 7.  Interactive Visualization Through Condensed Representations

In addition to facilitating machine learning and database operations, appropriate condensed representations can accelerate such tasks as volume rendering by allowing data items to be grouped and drawn together if they are similar to within specified error bounds. A large dataset can easily have more data elements than there are pixels on the display, and there is no sense rendering at a level of detail which will be invisible to the user. These hierarchical rendering methods are well studied in computer graphics and map well to condensed representations for knowledge discovery. For example, see be Laur & Hanrahan (1991).

## References

Kolojejchick, J. A., Roth, S. F., & Lucas, P. 1997, in Computer Graphics and Applications, Vol. 17 No 4, Information Appliances and Tools in Visage, 32, http://www.cs.cmu.edu/ sage/PDF/Appliances.pdf

Laur & Hanrahan 1991 in Computer Graphics, Vol. 25 No. 4, Hierarchical Splatting: A Progressive Refinement Algorithm for Volume Rendering, 285

Swayne, D. F., Cook, D., & Buja, A. 1998, in Journal of Computational and Graphical Statistics 7 (1), XGobi: Interactive Dynamic Data Visualization in the X Window System


text:
# Branching annihilating random walk on random regular graphs

## ACKNOWLEDGMENTS

Support from the Hungarian National Research Fund (T-23552) is acknowledged.