- Email: [email protected]
Numerical simulations of AGNs
Ady. Space Res. VoL 8. No. 2-3. pp. (2)119—(2)126. 1988 Printed in Great Britain. All rights reserved.
0273—1177/88 $0.00 + .50 Copyright © COSPAR
NUMERICAL SIMULATIONS OF AGNS John F. Hawley Theoretical Astrophysics 130—33, Caltech, Pasadena, CA 91125, U.S.A.
ABSTRACT The role that numerical simulations play in the study of active galactic nuclei is reviewed. Some of the advantages and limitations of numerical work are discussed and illustrated with applications from jet and disk studies. THE ROLE OF NUMERICAL SIMULATIONS As several of the contributors to this meeting have noted, we do not actually have observational proof that blac} holes are tile central engine of AGN. Nevertheless, there are strong theoretical grounds to believe that this is so, and I won’t review that here. The point nevertheless illustrates the gap between what we observe and what we theorize they do not connect even at this most basic level! -
Given this fundamental lack, we may well ask just what theorists are working on, particularly those who would do “detailed numerical simulations” of AGNs. Current theory consists of AGN paradigms, that is, pictures of possible black hole accretion mechanisms. These are based on general physical properties such as conservation laws, and they form a framework upon which more detailed work can be built. Numerical simulations can play a valuable role by investigating the details of these paradigms, details that depend on complex nonlinear physics, e.g., multidimensional compressible flow, nonlinear dynamics, and time dependence. The physics learned in such experiments can feed back into the paradigm; the numerical results may even lead to a new analytic result. By numerical simulations, I will mean the time-dependent solution of partial differential equations in one, two, or three spatial dimensions. This excludes valuable work in the solution of ordinary differential equations, but the general availability of good integration algorithms and codes, and the formal error bounds on can place on these solutions blurs the distinction between “analytic” and “numerical” for these types of problems. Since my own experience is in finite-differencing, my remarks will tend to be biased in that sense, particularly when describing the limitations of computer modeling. Below I present a very limited review of some recent and ongoing computations dealing with various aspects of active galaxies. My goal is not so much to discuss the physics of AGN as to provide an introduction to the role of computer simulations in investigating that physics. EVALUATING NUMERICAL SIMULATIONS Honest evaluation of simulations requires understanding what is involved in a large numerical
(2)119
(2)120
J. F. Hawley
project. Computational work has a kinship with experimental science; in both, the task is to set up a meaningful experiment, run it, and then make sense of the data generated. To begin, we must ask what questions can be profitably addressed with large scale numerical simulation? The answer is those problems that are solvable within the limitations of a typical multidimensional code: Resoluiion: The goal of a numerical simulation is to study a physical phenomenon that occurs on some characteristic length-scale. On the computer, a continuous function is reduced to a finite number of data points. For example, a wave will be represented by some discrete number of grid zones per wavelength. The longest wavelength that can be represented is length to the entire grid; the shortest wavelength is 2i~x.however, there is a crucial difference between mere representation and adequale representation. A feature represented by too few data points will be underresolved, and numerical errors may well dominate. Finite difference errors have their greatest effect on the short wavelength features. Diffusion and dispersion errors are typically the most serious. Diffusion is the preferential damping of short wavelength components: sharp features are spread out over many zones. Dispersion occurs because the finite difference scheme transports different wavelength components at different velocities: a traveling pulse will leave a wake, for instance. The physical length-scale of interest must be adequately resolved so that these numerical limitations do not dominate. Timescales: The issue is resolution in time. Numerical techniques that evolve forward in time by using the values of variables at the current timestep are termed time-explicit. Such schemes are restricted by the “Courant condition”, i.e., i~t must be sufficiently small so that no signal can propagate more than one grid zone per timestep. This means that an explicit hydrodynamics scheme is well suited to study the propagation and effects of shocks, pressure waves, and high Mach number flaws (i.e., tile dynamical timescate), but not for evolving long term effects. In AGN accretion disks there is a dynamical (orbital) timescale, a longer viscous disk evolution timescale, and finally the central engine fueling timescale over which the AGN as a whole evolves. Choosing one timescale means ignoring others; the consequences of such neglect are not always apparent. For example, modeling the long term behavior of an AGN central engine requires a description of the accretion disk that is nearly steady-state. This was one of the motivations behind the paradigm of the thick accretion torus. However, as the work on nonaxisymmetric disk instabilities has shown (discussed below), these “steady-state” configurations can be dynamically unstable. Obviously tile short timescale behavior of an unstable disk can not be neglected. Physics: The physics of a numerical simulation is incorporated in the equations that are being solved. The results must be evaluated in light of both what is and is noi included. The ideal would be to solve the system as completely as possible. This is difficult since, in many cases, the fundamental equations are unknown. When the relevant physics is not known, an approximation might be employed (e.g., parameterized star formation rate, parameterized viscosity), but then one is in danger of having the results tell more about the approximation than about the actual physical system. On the other hand, the physics left out of a simulation may invalidate the results if the inclusion of same would produce significant changes. Once one has determined a problem that is both relevant to theory, and within the capabilities of one’s hardware and software, there remains the decision as to which specific models to compute. Numerical models are discrete entities; a single model does not have much predictive power. Models must be linked through the parameier space of the problem, the set of all solutions to the equations, spanned by a set of fundamental parameters. These parameters must be chosen for their physical significance. Computing a set of models while varying one parameter and fixing the others provides a means for systematizing the numerical results. This approach allows a finite number of models to represent the complete solution space. Small changes in the input parameters should not produce large changes in the resulting model. Predictive power depends on covering parameter space with enough models. Once again it is a question of resolution!
Numerical Simulations
(2)12 1
Below, I describe several numerical simulations that address questions related to active galactic nuclei. The descriptions will be necessarily brief; those interested in further information should refer to the cited literature. JET SIMULATIONS During the 70’s, many theories were developed to explain the extended radio emission observed in certain active galaxies and quasars. The paradigm that eventually came to dominate was tile “twin exhaust” model, i.e., bi-directional collimated beams or jets featuring a terminal working surface associated with radio lobe hot spots /1/. Simply developing the idea of beamed exhaust did not resolve all the questions however. For example, is such a beam dynamically stable? Can one explain the observed knots of emission? Questions such as these are well suited to numerical attack, and simulations of idealized fluid jets provide some of the answers. Norman, Smarr, and Winkler /2-3/ investigated the behavior of an axisymmetric, hydrodynamic, pressure matched jet passing through an ambient medium. There are two fundamental parameters in this simplified system: the Mach number of the jet, and ij, the density ratio of the jet fluid to that of the ambient medium. By computing representative models from this parameter space, they were able to delineate general classes of jet behavior. For example, only high Mach number, light jets produce substantial backflowing cocoons. In their investigation of jet stability they were able to distinguish two distinct regions of nonlinear behavior. Jets with Mach numbers greater than 1 + ..J~jsaturated in the nonlinear regime as a series of reflecting shock waves. Such jets remain well collimated, and the reflecting shocks provide a possible mechanism for producing the observed radio emission knots. Jets with lower Mach numbers ~re disrupted by tile “Ordinary Mode” instability; the jet’s kinetic energy is dissipated in a series of shocks perpendicular to the beam direction. Results such as these have helped to establish the viability of the jet model. ACCRETION SIMULATIONS Theoretical studies of accretion disks pose legions of problems amenable to numerical simulation. Time dependent hydrodynamic codes can be useful for investigating any fluid flow process that occurs on dynamic or freefahl timescales. Examples include the effect of angular momentum on accretion flows, shock formation, multidimensional flow patterns, sonic-point studies, and disk formation. Another useful application is to follow the nonlinear evolution of instabilities. In. two dimensions the usual dynamic instabilities can be studied, e.g., Rayleigh-Taylor overturn in a thick disk. In three dimensions the nonaxisymmetric instability first described by Papaloizou and Pringle /4/ is of current interest. Here I discuss some axisymmetric simulations, and, in the following section, the nonaxisymmetric instability work. Ideal general relativistic fluid dynamics has been used to extend the spherical accretion problem to the infall of fluid with angular momentum /5-6/. A simple parameterization of tile problem is to consider freefalling gas with a range of constant angular momentum values. Simulations reveal that even angular momentum as small as the “marginally stable” value, i.e., tile angular momentum of the last stable black hole orbit, forms shocks in tile flow as the fluid converges towards the equator, creating an evacuated “funnel” along the fluid’s spin axis. Increasing the angular momentum above this value produces the most interesting results. The presence of a centrifugal barrier leads to shock formation and “splash back” from tile black hole. Where there was once supersonic inflowinggas, there is now a hot thick disk. The properties of this thick disk can be compared with the analytic thick disk model. Tile two are remarkably similar, considering the dynamic genesis of the numerical model. Changing the angular momentum distribution of the inflowing gas from constant to a value proportional to sin2 0 drastically changes the the inflow and the structure of tile thick disk that is formed /6/. This has the unfortunate implication that the characteristics of a dynamically
(2)122
i. F. Hawley
formed disk depend strongly on the distribution of angular momentum in the gas. Neither the analytic nor the numerical models provide a self-consistent prescription, except in the limited cases of constant or Keplerian angular momentum distributions. It is also clear that tile limited number of numerical models that have been calculated to date undersample the parameter space. The ultimate goal is to simulate a dynamic accretion disk that includes viscosity, radiation transport, and magnetohydrodynamics. Some of the difficulties associated with this goal are the lack of a good theory for viscosity, the complexity of the radiation transport equations, and the disparate viscous and dynamic timescales in the disk. Eggum, Coroniti, and Katz /7/ computed a viscous disk using flux-limited diffusion for the radiation transport and a parameterized form for the viscosity. Their mildly super-Eddington models display several novel features, including large scale convection cells in the center of the disk, and a radiation driven wind emerging from the disk surface. Despite the super-Eddington accretion rates, only 20 per cent of the gravitational energy released in the flow escapes to infinity. As with tile pure hydrodynamic accretion flows described above, the Eggum ci. al. models suggest ways in which our fundamental accretion paradigms might be modified. Further simulations will be required to determine the extent to which the observed effects depend on the details of the assumptions employed. NONAXISYMMETRIC DISK INSTABILITIES A popular model for accretion disks around black holes is the thick accretion disk, or accretion torus. The evacuated “funnel” of such a torus has been proposed as a possible jet acceleration and collimation site. However, nonaxisymmetric linear perturbation analysis of thick disks demonstrates that in many general circumstances, these tori are unstable. This poses several questions: Do these instabilities saturate at low levels or do they disrupt the disk? Do these nonaxisymmetnc instabilities transport significant angular momentum? Are thick accretion disks and narrow funnels viable? The answers to these questions depend on understanding the full nature of the instability. Linear perturbation analysis /4,8-10/ gives tile fundamental parameters of this problem, the angular momentum distribution through the disk, and tile radial width of the torus, from inner to outer edge. Typically, the angular momentum distribution is further simplified by assuming a radial dependence of the form ~(2~) (for a Newtonian potential). Linear theory provides tile form of the unstable modes and their growth rates, but we must rely on time-dependent hydrodynamic simulations to investigate the nonlinear evolution. The first such simulation was performed by Zurek and Benz /11/ with a three-dimensional smooth particle hydrodynamics code. They ran several models and found that an initially constant angular momentum torus experiences a rapid readjustment of its averaged angular momentum. Angular momentum is transported outwards, and the averaged angular momentum parameter q decreases from q = 2 (constant angular momentum) to a value around 1.7 (q = 1.5 corresponds to a Keplerian distribution). A similar three-dimensional calculation with a finite difference code could be performed, but would either tax the capacity of the current supercomputers, or would suffer from lack of resolution. (Indeed, Zurek and Benz’s code used only 1000 particles to represent the torus and must also be considered underresolved, although the relationship of number of partvicles to finite difference grid zones is not clear.) However, as shown by the linear perturbation analysis by Goldreich, ci. al. /8/ the principle eigenmode of the instability is nearly independent of z for a slender torus (i.e., the width of the disk is much less than the radius from the central mass). It thus becomes possible to study these modes using a two-dimensional simulation in the equatorial plane using (r, q5) coordinates. A series of such models of slender tori /12/ has faithfully reproduced the eigenmodes and growth rates predicted by the linear theory, and has shown that at nonlinear saturation, the instability causes the torus to break into a series of m orbiting fluid blobs, designated “planets”, where m is the azimuthal wavenumber of the fastest growing mode. An example is shown in Figure 1. An example of tile symbiotic relationship between numerical and analytic work is that tile discovery of these fluid “planets” has led to the development of an
Numerical Simulations
(2)123
O~,,’
Figure 1: Density plot showing the evolution of a slender 2D disk (on the left) to a nonaxisymmetric highly distorted planet. The change occurs in 9.6 orbits. analytic planet solution /13/.
It is a bit more problematic to extend the same two dimensional simulation techniques to the case of tlte radially wide torus. A full 3D linear perturbation analysis by Frank /14/ indicates that the unstable eigenrnodes in wide, thick tori are not similarly independent of height z. however, extending the 2D numerical simulations over the complete range of disk widths allows a detailed comparison with linear theory /10/. Further, such work is a clear prerequisite to performing more complex 3D simulations. The 2D work indicates that tile planet mode is a special case for tile slender torus. Wider disks saturate with a low amplitude dcnsity enhancement at the disk pressure maximum that trails a spiral wave through the disk (Fig. 2). In the outer regions of the disk, the wavelength of the eigenmode becomes increasingly short. The numerical simulations have shown that inaccurate growth rates result from the failure to adequately resolve those waves. This means that accurate three dimensional simulations may require more resolution than is currently feasible. These numerical simulations add to our understanding of the nonlinear behavior of tile instability, but their applicability to astrophysical accretion disks is more problematic. The simulations themselves demonstrate that the initial conditions used, e.g., constant angular momentum, are implausible. Tile nonaxisymmetnic structures that develop are themselves dynamically evolving, albeit on a slightly longer timescale. The models suggest new paradigms, but do not establish them. Further work will be required before we can determine what constitutes a viable steadystate thick disk model. MIlD PROSPECTS The observation that some jets may be overpressured and nevertheless confined suggests that the hydrodynamic model may be insufficient. Large scale coherent magnetic fields provide an additional confining and collimating force. As in the hydrodynamic case, numerical simulations will play an important role in refining the jet paradigm, but the work is only now beginning. Tile addition of the magnetohiydrodynamical terms to the jet equations adds considerable complexity. The ratio of magnetic to fluid pressure is one new parameter to consider; also one must model both Loroidal and poloidal fields. Some of the preliminary work along these lines is described
(2)124
J. F. Hawley
S Figure 2: Evolution to a nonaxisymnletric spiral wave. A wide constant angular momentum model forms a spiral wave trailing from a density enhancement, as this plot of logarithmic density shows. by Kevin Lind in his talk (see also /15/). In addition, first attempts have been made to add magnetic fields to Norman’s jet code /16/. Even these early results have uncovered interesting new jet phenomenology. Another area of interest is tile study of relativistic MIlD near the black hole. Sterl Phinney (this volume) has reviewed some aspects of black hole MIlD and its role in powering and collimating jets. Numerical simulations of dynamic MIlD flows around black holes were attempted by James Wilson over a decade ago /17/. The focus of current research is to develop and apply up-to-date finite difference techniques to tile general relativistic MIlD equations /18/. CONCLUSION I hope that this brief presentation has conveyed some sense of what numerical simulations can and have contributed to tile study of active galactic nuclei and jets. At the risk of being labeled overly cautious, I will again emphasize my belief that it is essential to apply tile most conservative and rigorous of standards to numerical simulations, from the initial development stage to the final conclusions. This is analogous to experimental science in every sense of the word; we must understand our apparatus, use the best technique, and honestly evaluate the data, always keeping the error bars in mind. If we do this, large scale computer modeling will undoubtably have a significant impact in advancing our understanding of complicated systems such as AGN. This work was supported in part by NSF grants P11Y85-11426 and AST86-15325. A portion of tile computational work described here was performed at the National Center for Supercomputer Applications.
Numerical Simulations
(2)125
REFERENCES 1. R. D. Blandford, and M. J. Rees, “A Twin Exhaust Model for Double Radio Sources,” Mon. Not. R. astr. Soc., 169, 395-415 (1974) 2. M. L. Norman, L. L. Smarr, and K.-H. A. Winkler, “Shocks, Interfaces, and Patterns in Supersonic Jets,” Physica 12D, 83-106 (1984) 3. M. L. Norman, and K.-H. A. Winkier, “Supersonic Jets,” Los Alamos Science, No. 12, 38-71 (1985) 4. J. C. B. Papaioizou, and J. E. Pringle, “The Dynamical Stability of Differentially Rotating Discs with Constant Specific Angular Momentum,” Mon. Not. R. astr. Soc., 208, 721-750 (1984) 5. J. F. Hawley and L. L. Smarr, “New Paradigms for Black hole Accretion from High Resolution Supercomputer Experiments,” in Magnetospheric Phenomena in Astrophysics, ed. R. Epstein and W. Feldman (New York: AlP), p. 263 (1986). 6. J. F. Hawley, “Hydrodynamics Near tile Central Engine,” in Radiation hydrodynamics in Stars and Compact Objects, ed. D. Mihalas and K.-H. A. Winkler (New York: Springer Verlag), p. 369 (1986) 7. G. E. Eggum, F. V. Coroniti, and J. I. Katz, “Jet Production in Super-Eddington Accretion Disks,” Astrophys. J., 298, L41 (1985) 8. P. Goldreich, J. Goodman, and R. Narayan, “The Stability of Accretion Tori - I. Longwavelength Modd of Slender Tori,” Mon. Not. R. astr. Soc., 221, 339 (1986) 9. W. Glatzel, “On the Stability of Compressible Differentially Rotating Cylinders,” Mon. Not. R. astr. Soc., 225, 227 (1987) 10. 0. M. Blaes, and J. F. Hawley, “Nonaxisymmetric Disk Instabilities: A Linear and Nonlinear Synthesis,” Astrophys. J., in press (1987) 11. W. H. Zurek, and W. Benz, “Redistribution of Angular Momentum by Nonaxisymmetric Instabilities in a Thick Accretion Disk,” Astrophys. J., 308, 123-133 (1986) 12. J. F. Hawley, “Non-linear Evolution of a Non-axisymmetric Disc Instability,” Mon. Not. fl. a.str. Soc., 225, 677 (1987) 13. 3. Goodman, R. Narayan, and P. Goldreich, “The Stability of Accretion Tori— II. Nonlinear Evolution to Discrete Planets,” Mon. Not. R. astr. Soc., 225, 695 (1987) 14. 3. Frank, talk presented at the Institute for Theoretical Physics Workshop on Active Galactic Nuclei, January. 1987, Santa Barbara. 15. K. Lind, “Observations and Gas Dynamics of Extragalactic Radio Jets,” thesis, California Institute of Technology (1987) 16. D. A. Clarke, M. L. Norman, and J. 0. Burns, “Numerical Simulations of Magnetically Confined Jet,” Astrophys. J., 311, L63 (1986) 17. J. R. Wilson, “Some Magnetic Effects in Stellar Collapse and Accretion,” Ann. N.Y. Acad. Sci., 262, 123 (1975)
(2)126
J. F. Hawley
18. C. R. Evans, and J.~F. Hawley, “Simulation of General Relativistic Magnetohiydrodynamic Flows: A Constrained Transport Method,” preprint (1987)
0273—1177/88 $0.00 + .50 Copyright © COSPAR
NUMERICAL SIMULATIONS OF AGNS John F. Hawley Theoretical Astrophysics 130—33, Caltech, Pasadena, CA 91125, U.S.A.
ABSTRACT The role that numerical simulations play in the study of active galactic nuclei is reviewed. Some of the advantages and limitations of numerical work are discussed and illustrated with applications from jet and disk studies. THE ROLE OF NUMERICAL SIMULATIONS As several of the contributors to this meeting have noted, we do not actually have observational proof that blac} holes are tile central engine of AGN. Nevertheless, there are strong theoretical grounds to believe that this is so, and I won’t review that here. The point nevertheless illustrates the gap between what we observe and what we theorize they do not connect even at this most basic level! -
Given this fundamental lack, we may well ask just what theorists are working on, particularly those who would do “detailed numerical simulations” of AGNs. Current theory consists of AGN paradigms, that is, pictures of possible black hole accretion mechanisms. These are based on general physical properties such as conservation laws, and they form a framework upon which more detailed work can be built. Numerical simulations can play a valuable role by investigating the details of these paradigms, details that depend on complex nonlinear physics, e.g., multidimensional compressible flow, nonlinear dynamics, and time dependence. The physics learned in such experiments can feed back into the paradigm; the numerical results may even lead to a new analytic result. By numerical simulations, I will mean the time-dependent solution of partial differential equations in one, two, or three spatial dimensions. This excludes valuable work in the solution of ordinary differential equations, but the general availability of good integration algorithms and codes, and the formal error bounds on can place on these solutions blurs the distinction between “analytic” and “numerical” for these types of problems. Since my own experience is in finite-differencing, my remarks will tend to be biased in that sense, particularly when describing the limitations of computer modeling. Below I present a very limited review of some recent and ongoing computations dealing with various aspects of active galaxies. My goal is not so much to discuss the physics of AGN as to provide an introduction to the role of computer simulations in investigating that physics. EVALUATING NUMERICAL SIMULATIONS Honest evaluation of simulations requires understanding what is involved in a large numerical
(2)119
(2)120
J. F. Hawley
project. Computational work has a kinship with experimental science; in both, the task is to set up a meaningful experiment, run it, and then make sense of the data generated. To begin, we must ask what questions can be profitably addressed with large scale numerical simulation? The answer is those problems that are solvable within the limitations of a typical multidimensional code: Resoluiion: The goal of a numerical simulation is to study a physical phenomenon that occurs on some characteristic length-scale. On the computer, a continuous function is reduced to a finite number of data points. For example, a wave will be represented by some discrete number of grid zones per wavelength. The longest wavelength that can be represented is length to the entire grid; the shortest wavelength is 2i~x.however, there is a crucial difference between mere representation and adequale representation. A feature represented by too few data points will be underresolved, and numerical errors may well dominate. Finite difference errors have their greatest effect on the short wavelength features. Diffusion and dispersion errors are typically the most serious. Diffusion is the preferential damping of short wavelength components: sharp features are spread out over many zones. Dispersion occurs because the finite difference scheme transports different wavelength components at different velocities: a traveling pulse will leave a wake, for instance. The physical length-scale of interest must be adequately resolved so that these numerical limitations do not dominate. Timescales: The issue is resolution in time. Numerical techniques that evolve forward in time by using the values of variables at the current timestep are termed time-explicit. Such schemes are restricted by the “Courant condition”, i.e., i~t must be sufficiently small so that no signal can propagate more than one grid zone per timestep. This means that an explicit hydrodynamics scheme is well suited to study the propagation and effects of shocks, pressure waves, and high Mach number flaws (i.e., tile dynamical timescate), but not for evolving long term effects. In AGN accretion disks there is a dynamical (orbital) timescale, a longer viscous disk evolution timescale, and finally the central engine fueling timescale over which the AGN as a whole evolves. Choosing one timescale means ignoring others; the consequences of such neglect are not always apparent. For example, modeling the long term behavior of an AGN central engine requires a description of the accretion disk that is nearly steady-state. This was one of the motivations behind the paradigm of the thick accretion torus. However, as the work on nonaxisymmetric disk instabilities has shown (discussed below), these “steady-state” configurations can be dynamically unstable. Obviously tile short timescale behavior of an unstable disk can not be neglected. Physics: The physics of a numerical simulation is incorporated in the equations that are being solved. The results must be evaluated in light of both what is and is noi included. The ideal would be to solve the system as completely as possible. This is difficult since, in many cases, the fundamental equations are unknown. When the relevant physics is not known, an approximation might be employed (e.g., parameterized star formation rate, parameterized viscosity), but then one is in danger of having the results tell more about the approximation than about the actual physical system. On the other hand, the physics left out of a simulation may invalidate the results if the inclusion of same would produce significant changes. Once one has determined a problem that is both relevant to theory, and within the capabilities of one’s hardware and software, there remains the decision as to which specific models to compute. Numerical models are discrete entities; a single model does not have much predictive power. Models must be linked through the parameier space of the problem, the set of all solutions to the equations, spanned by a set of fundamental parameters. These parameters must be chosen for their physical significance. Computing a set of models while varying one parameter and fixing the others provides a means for systematizing the numerical results. This approach allows a finite number of models to represent the complete solution space. Small changes in the input parameters should not produce large changes in the resulting model. Predictive power depends on covering parameter space with enough models. Once again it is a question of resolution!
Numerical Simulations
(2)12 1
Below, I describe several numerical simulations that address questions related to active galactic nuclei. The descriptions will be necessarily brief; those interested in further information should refer to the cited literature. JET SIMULATIONS During the 70’s, many theories were developed to explain the extended radio emission observed in certain active galaxies and quasars. The paradigm that eventually came to dominate was tile “twin exhaust” model, i.e., bi-directional collimated beams or jets featuring a terminal working surface associated with radio lobe hot spots /1/. Simply developing the idea of beamed exhaust did not resolve all the questions however. For example, is such a beam dynamically stable? Can one explain the observed knots of emission? Questions such as these are well suited to numerical attack, and simulations of idealized fluid jets provide some of the answers. Norman, Smarr, and Winkler /2-3/ investigated the behavior of an axisymmetric, hydrodynamic, pressure matched jet passing through an ambient medium. There are two fundamental parameters in this simplified system: the Mach number of the jet, and ij, the density ratio of the jet fluid to that of the ambient medium. By computing representative models from this parameter space, they were able to delineate general classes of jet behavior. For example, only high Mach number, light jets produce substantial backflowing cocoons. In their investigation of jet stability they were able to distinguish two distinct regions of nonlinear behavior. Jets with Mach numbers greater than 1 + ..J~jsaturated in the nonlinear regime as a series of reflecting shock waves. Such jets remain well collimated, and the reflecting shocks provide a possible mechanism for producing the observed radio emission knots. Jets with lower Mach numbers ~re disrupted by tile “Ordinary Mode” instability; the jet’s kinetic energy is dissipated in a series of shocks perpendicular to the beam direction. Results such as these have helped to establish the viability of the jet model. ACCRETION SIMULATIONS Theoretical studies of accretion disks pose legions of problems amenable to numerical simulation. Time dependent hydrodynamic codes can be useful for investigating any fluid flow process that occurs on dynamic or freefahl timescales. Examples include the effect of angular momentum on accretion flows, shock formation, multidimensional flow patterns, sonic-point studies, and disk formation. Another useful application is to follow the nonlinear evolution of instabilities. In. two dimensions the usual dynamic instabilities can be studied, e.g., Rayleigh-Taylor overturn in a thick disk. In three dimensions the nonaxisymmetric instability first described by Papaloizou and Pringle /4/ is of current interest. Here I discuss some axisymmetric simulations, and, in the following section, the nonaxisymmetric instability work. Ideal general relativistic fluid dynamics has been used to extend the spherical accretion problem to the infall of fluid with angular momentum /5-6/. A simple parameterization of tile problem is to consider freefalling gas with a range of constant angular momentum values. Simulations reveal that even angular momentum as small as the “marginally stable” value, i.e., tile angular momentum of the last stable black hole orbit, forms shocks in tile flow as the fluid converges towards the equator, creating an evacuated “funnel” along the fluid’s spin axis. Increasing the angular momentum above this value produces the most interesting results. The presence of a centrifugal barrier leads to shock formation and “splash back” from tile black hole. Where there was once supersonic inflowinggas, there is now a hot thick disk. The properties of this thick disk can be compared with the analytic thick disk model. Tile two are remarkably similar, considering the dynamic genesis of the numerical model. Changing the angular momentum distribution of the inflowing gas from constant to a value proportional to sin2 0 drastically changes the the inflow and the structure of tile thick disk that is formed /6/. This has the unfortunate implication that the characteristics of a dynamically
(2)122
i. F. Hawley
formed disk depend strongly on the distribution of angular momentum in the gas. Neither the analytic nor the numerical models provide a self-consistent prescription, except in the limited cases of constant or Keplerian angular momentum distributions. It is also clear that tile limited number of numerical models that have been calculated to date undersample the parameter space. The ultimate goal is to simulate a dynamic accretion disk that includes viscosity, radiation transport, and magnetohydrodynamics. Some of the difficulties associated with this goal are the lack of a good theory for viscosity, the complexity of the radiation transport equations, and the disparate viscous and dynamic timescales in the disk. Eggum, Coroniti, and Katz /7/ computed a viscous disk using flux-limited diffusion for the radiation transport and a parameterized form for the viscosity. Their mildly super-Eddington models display several novel features, including large scale convection cells in the center of the disk, and a radiation driven wind emerging from the disk surface. Despite the super-Eddington accretion rates, only 20 per cent of the gravitational energy released in the flow escapes to infinity. As with tile pure hydrodynamic accretion flows described above, the Eggum ci. al. models suggest ways in which our fundamental accretion paradigms might be modified. Further simulations will be required to determine the extent to which the observed effects depend on the details of the assumptions employed. NONAXISYMMETRIC DISK INSTABILITIES A popular model for accretion disks around black holes is the thick accretion disk, or accretion torus. The evacuated “funnel” of such a torus has been proposed as a possible jet acceleration and collimation site. However, nonaxisymmetric linear perturbation analysis of thick disks demonstrates that in many general circumstances, these tori are unstable. This poses several questions: Do these instabilities saturate at low levels or do they disrupt the disk? Do these nonaxisymmetnc instabilities transport significant angular momentum? Are thick accretion disks and narrow funnels viable? The answers to these questions depend on understanding the full nature of the instability. Linear perturbation analysis /4,8-10/ gives tile fundamental parameters of this problem, the angular momentum distribution through the disk, and tile radial width of the torus, from inner to outer edge. Typically, the angular momentum distribution is further simplified by assuming a radial dependence of the form ~(2~) (for a Newtonian potential). Linear theory provides tile form of the unstable modes and their growth rates, but we must rely on time-dependent hydrodynamic simulations to investigate the nonlinear evolution. The first such simulation was performed by Zurek and Benz /11/ with a three-dimensional smooth particle hydrodynamics code. They ran several models and found that an initially constant angular momentum torus experiences a rapid readjustment of its averaged angular momentum. Angular momentum is transported outwards, and the averaged angular momentum parameter q decreases from q = 2 (constant angular momentum) to a value around 1.7 (q = 1.5 corresponds to a Keplerian distribution). A similar three-dimensional calculation with a finite difference code could be performed, but would either tax the capacity of the current supercomputers, or would suffer from lack of resolution. (Indeed, Zurek and Benz’s code used only 1000 particles to represent the torus and must also be considered underresolved, although the relationship of number of partvicles to finite difference grid zones is not clear.) However, as shown by the linear perturbation analysis by Goldreich, ci. al. /8/ the principle eigenmode of the instability is nearly independent of z for a slender torus (i.e., the width of the disk is much less than the radius from the central mass). It thus becomes possible to study these modes using a two-dimensional simulation in the equatorial plane using (r, q5) coordinates. A series of such models of slender tori /12/ has faithfully reproduced the eigenmodes and growth rates predicted by the linear theory, and has shown that at nonlinear saturation, the instability causes the torus to break into a series of m orbiting fluid blobs, designated “planets”, where m is the azimuthal wavenumber of the fastest growing mode. An example is shown in Figure 1. An example of tile symbiotic relationship between numerical and analytic work is that tile discovery of these fluid “planets” has led to the development of an
Numerical Simulations
(2)123
O~,,’
Figure 1: Density plot showing the evolution of a slender 2D disk (on the left) to a nonaxisymmetric highly distorted planet. The change occurs in 9.6 orbits. analytic planet solution /13/.
It is a bit more problematic to extend the same two dimensional simulation techniques to the case of tlte radially wide torus. A full 3D linear perturbation analysis by Frank /14/ indicates that the unstable eigenrnodes in wide, thick tori are not similarly independent of height z. however, extending the 2D numerical simulations over the complete range of disk widths allows a detailed comparison with linear theory /10/. Further, such work is a clear prerequisite to performing more complex 3D simulations. The 2D work indicates that tile planet mode is a special case for tile slender torus. Wider disks saturate with a low amplitude dcnsity enhancement at the disk pressure maximum that trails a spiral wave through the disk (Fig. 2). In the outer regions of the disk, the wavelength of the eigenmode becomes increasingly short. The numerical simulations have shown that inaccurate growth rates result from the failure to adequately resolve those waves. This means that accurate three dimensional simulations may require more resolution than is currently feasible. These numerical simulations add to our understanding of the nonlinear behavior of tile instability, but their applicability to astrophysical accretion disks is more problematic. The simulations themselves demonstrate that the initial conditions used, e.g., constant angular momentum, are implausible. Tile nonaxisymmetnic structures that develop are themselves dynamically evolving, albeit on a slightly longer timescale. The models suggest new paradigms, but do not establish them. Further work will be required before we can determine what constitutes a viable steadystate thick disk model. MIlD PROSPECTS The observation that some jets may be overpressured and nevertheless confined suggests that the hydrodynamic model may be insufficient. Large scale coherent magnetic fields provide an additional confining and collimating force. As in the hydrodynamic case, numerical simulations will play an important role in refining the jet paradigm, but the work is only now beginning. Tile addition of the magnetohiydrodynamical terms to the jet equations adds considerable complexity. The ratio of magnetic to fluid pressure is one new parameter to consider; also one must model both Loroidal and poloidal fields. Some of the preliminary work along these lines is described
(2)124
J. F. Hawley
S Figure 2: Evolution to a nonaxisymnletric spiral wave. A wide constant angular momentum model forms a spiral wave trailing from a density enhancement, as this plot of logarithmic density shows. by Kevin Lind in his talk (see also /15/). In addition, first attempts have been made to add magnetic fields to Norman’s jet code /16/. Even these early results have uncovered interesting new jet phenomenology. Another area of interest is tile study of relativistic MIlD near the black hole. Sterl Phinney (this volume) has reviewed some aspects of black hole MIlD and its role in powering and collimating jets. Numerical simulations of dynamic MIlD flows around black holes were attempted by James Wilson over a decade ago /17/. The focus of current research is to develop and apply up-to-date finite difference techniques to tile general relativistic MIlD equations /18/. CONCLUSION I hope that this brief presentation has conveyed some sense of what numerical simulations can and have contributed to tile study of active galactic nuclei and jets. At the risk of being labeled overly cautious, I will again emphasize my belief that it is essential to apply tile most conservative and rigorous of standards to numerical simulations, from the initial development stage to the final conclusions. This is analogous to experimental science in every sense of the word; we must understand our apparatus, use the best technique, and honestly evaluate the data, always keeping the error bars in mind. If we do this, large scale computer modeling will undoubtably have a significant impact in advancing our understanding of complicated systems such as AGN. This work was supported in part by NSF grants P11Y85-11426 and AST86-15325. A portion of tile computational work described here was performed at the National Center for Supercomputer Applications.
Numerical Simulations
(2)125
REFERENCES 1. R. D. Blandford, and M. J. Rees, “A Twin Exhaust Model for Double Radio Sources,” Mon. Not. R. astr. Soc., 169, 395-415 (1974) 2. M. L. Norman, L. L. Smarr, and K.-H. A. Winkler, “Shocks, Interfaces, and Patterns in Supersonic Jets,” Physica 12D, 83-106 (1984) 3. M. L. Norman, and K.-H. A. Winkier, “Supersonic Jets,” Los Alamos Science, No. 12, 38-71 (1985) 4. J. C. B. Papaioizou, and J. E. Pringle, “The Dynamical Stability of Differentially Rotating Discs with Constant Specific Angular Momentum,” Mon. Not. R. astr. Soc., 208, 721-750 (1984) 5. J. F. Hawley and L. L. Smarr, “New Paradigms for Black hole Accretion from High Resolution Supercomputer Experiments,” in Magnetospheric Phenomena in Astrophysics, ed. R. Epstein and W. Feldman (New York: AlP), p. 263 (1986). 6. J. F. Hawley, “Hydrodynamics Near tile Central Engine,” in Radiation hydrodynamics in Stars and Compact Objects, ed. D. Mihalas and K.-H. A. Winkler (New York: Springer Verlag), p. 369 (1986) 7. G. E. Eggum, F. V. Coroniti, and J. I. Katz, “Jet Production in Super-Eddington Accretion Disks,” Astrophys. J., 298, L41 (1985) 8. P. Goldreich, J. Goodman, and R. Narayan, “The Stability of Accretion Tori - I. Longwavelength Modd of Slender Tori,” Mon. Not. R. astr. Soc., 221, 339 (1986) 9. W. Glatzel, “On the Stability of Compressible Differentially Rotating Cylinders,” Mon. Not. R. astr. Soc., 225, 227 (1987) 10. 0. M. Blaes, and J. F. Hawley, “Nonaxisymmetric Disk Instabilities: A Linear and Nonlinear Synthesis,” Astrophys. J., in press (1987) 11. W. H. Zurek, and W. Benz, “Redistribution of Angular Momentum by Nonaxisymmetric Instabilities in a Thick Accretion Disk,” Astrophys. J., 308, 123-133 (1986) 12. J. F. Hawley, “Non-linear Evolution of a Non-axisymmetric Disc Instability,” Mon. Not. fl. a.str. Soc., 225, 677 (1987) 13. 3. Goodman, R. Narayan, and P. Goldreich, “The Stability of Accretion Tori— II. Nonlinear Evolution to Discrete Planets,” Mon. Not. R. astr. Soc., 225, 695 (1987) 14. 3. Frank, talk presented at the Institute for Theoretical Physics Workshop on Active Galactic Nuclei, January. 1987, Santa Barbara. 15. K. Lind, “Observations and Gas Dynamics of Extragalactic Radio Jets,” thesis, California Institute of Technology (1987) 16. D. A. Clarke, M. L. Norman, and J. 0. Burns, “Numerical Simulations of Magnetically Confined Jet,” Astrophys. J., 311, L63 (1986) 17. J. R. Wilson, “Some Magnetic Effects in Stellar Collapse and Accretion,” Ann. N.Y. Acad. Sci., 262, 123 (1975)
(2)126
J. F. Hawley
18. C. R. Evans, and J.~F. Hawley, “Simulation of General Relativistic Magnetohiydrodynamic Flows: A Constrained Transport Method,” preprint (1987)