In this section the results from a range of three dimensional simulations are presented. The simulations are listed in the summary Table which gives a model designation, the grid resolution used, the gravitational potential, the initial torus angular momentum distribution q, the initial magnetic field topology, and the time to which the simulation was run.
|CK1||90 × 80 × 24||Cylindrical Newton||1.5||Bz ~ sin(R)||188|
|CK2||90 × 80 × 24||Cylindrical Newton||1.5||Toroidal||285|
|CT1||90 × 80 × 24||Cylindrical Newton||2||Vertical||102|
|CT2||90 × 80 × 24||Cylindrical Newton||2||Toroidal||403|
|CT3||128 × 64 × 32||Cylindrical Pseudo-Newton||2||Vertical||420|
|GT1||128 × 64 × 128||Global Pseudo-Newton||2||Poloidal loops||780|
|GT2||120 × 64 × 90||Global Pseudo-Newton||2||Poloidal loops||727|
|GT3||128 × 128 × 128||Global Pseudo-Newton||2||Toroidal||727|
|GT4||128 × 128 × 128||Global Pseudo-Newton||1.68||Poloidal loops||1283|
The first set of simulations make use of a cylindrical gravitational potential ~ 1/R (the cylindrical disk limit). This is a significant simplification, as it permits the use of periodic boundary conditions in the vertical direction. Further, with this potential there is only one important vertical lengthscale, namely that of the MRI; this reduces the number of z grid zones required. Finally, in the cylindrical limit one can study the evolution of models with vertical fields without the stringent Courant limits due to the high Alfven speeds associated with strong fields passing through the low density region above the disk.
The first two cylindrical simulations assume a Keplerian initial disk. Keplerian simulations provide an immediate comparison with the local shearing box results which also, for the most part, assumed Keplerian flow. Consider first the run labeled CK1 in the summary Table; this is a cylindrical computational domain (R, phi,z) running from 1 to 4 in R, 0 to 1 in z, and 0 to pi/2 in phi. The grid has 90 × 80 × 24 zones. An outflow boundary condition is used at both the outer and inner radial boundaries, and periodic boundary conditions are used for phi and z. The initial disk has a constant mass density from r=1.5 to the outer boundary. The adiabatic sound speed is constant and equal to 5% of the orbital velocity at the inner edge of the disk. The initial magnetic field is vertical and proportional to sin(R)/R between R=1.5 and 3.5, with a maximum strength corresponding to beta=400. The strength of the field, and hence the Alfven speed vA, was chosen so that the characteristic wavelength of the MRI was lambdaMRI = 2 pi vA/ Omega = 0.38 at the location of the field strength maximum. This ensures that the z domain size and the vertical grid resolution will be adequate to resolve the fastest growing modes. Random nonaxisymmetric pressure fluctuations are added to create a full range of low amplitude initial perturbations.
As the evolution proceeds, the MRI sets in, field energy is amplified, and soon the characteristic radial streaming structures (referred to as the channel solution) of the vertical field instability appear, much as they do in the local shearing box models. In the present simulation, these structures develop first at the inner part of the disk where the rotation frequency is the highest. The amplitude of the MRI becomes nonlinear by 3 orbits at the center of the grid (Porb = 24.8), and filaments of strong magnetic field are carried inward and outward by fluid elements well out of Keplerian balance. These reach the outer part of the disk even before the local MRI in that region becomes fully nonlinear. Thus there are two immediate global effects not seen in local simulations: linear growth rates that vary strongly with radius (omegaMRI ~ Omega ~ R-3/2, and extended radial motions of significantly non-Keplerian fluid.
In local simulations with initial vertical fields that vary sinusoidally with radius, the initial phase of the instability promptly breaks down into MHD turbulence. This happens in the global simulation as well. After the onset of turbulence, the disk displays many tightly wrapped (i.e., large radial wavenumber) trailing spiral features with low azimuthal wavenumber, m. The magnetic field energy saturates at about beta=4, with the toroidal component dominant: BR2/Bphi2 = 0.075 and Bz2/BR2 = 0.34. Similar ratios were found in the shearing box simulations with zero net initial magnetic field (HGB2).
The MHD turbulence produces rapid angular momentum transport and mass accretion, with alphapeaking at 0.21 at t=90, and declining beyond this point, dropping to 0.06 at the end of the simulation. The Maxwell stress displays a high correlation with the total magnetic pressure. After an early peak at 0.68, alpham declines with time; at the end of the simulation alpham = 0.35. The ratio Maxwell to Reynolds stress is about 3 at the time of the largest total stress and rises to about 9 at the end of the simulation. The overall angular momentum distribution remains nearly Keplerian throughout the simulation, becoming slightly sub-Keplerian outside of the original inner disk edge, and slightly super-Keplerian in the region between this point and the inner grid boundary.
The mass flux varies both in space and time. After 8 orbits of evolution, over half of the initial disk mass has been lost through the boundaries, particularly the outer boundary. Since the initial Keplerian disk abuts the outer boundary, and the mass increases with radius for a constant initial rho, it is relatively easy for a considerable amount of the mass to leave the grid, driven in part by a pressure gradient created by a rarefaction wave.
This simulation is similar to that of Armitage (1998) except that it uses a smaller domain in the angular direction, and a constant initial density in the disk. Qualitatively, the two runs exhibit similar behavior, and both are consistent with results from local simulations. One difference noted by Armitage, is that alpha ~ 0.1 in the global simulation; this is larger than the value of alpha ~ 0.01 seen in local simulations with zero net vertical field. The larger alpha is, however, consistent with what is seen in local simulations with net vertical field. The global simulations have much larger radial and azimuthal wavelengths available to the MRI. Further, the radial scale over which the vertical field sums to zero in the global simulation is much larger than in the local shearing box. Here Bz initially varies between r=1.5 and 3.5; this range has a large variation in Omega and hence a large range in growth rates of the MRI. Finally, the presence of outflow radial boundaries means the net field need not remain zero over the simulation. In the local simulation, periodic boundary ensure that the net field cannot change, and by definition, the vertical field sums to zero over a Delta R which is small. Reconnection in the local box should be much more efficient, reducing the total magnetic pressure and hence the total stress.
The next simulation, labeled CK2, begins with the same Keplerian disk, but with an initial toroidal field of constant strength, beta=4, lying between r=2 and 3.5. Generally speaking, toroidal field simulations are more challenging than vertical field simulations. Although the presence of the linear instability is independent of the orientation of the background magnetic field, the fastest growing toroidal field modes are associated with large vertical and radial wavenumbers. Further, the critical azimuthal wavenumber is large, unless the magnetic energy is comparable to the rotation energy. Here mcrit = 2 pi R/lambdaMRI ~ 40. Since there are 80 azimuthal grid zones over pi/2 in angle, the critical wavenumbers of the instability will be resolved for this choice of field strength. But it is clear that it takes a relatively strong toroidal field to do so for a reasonable number of grid zones. This illustrates one of the difficulties with global models: it is necessary to reconcile the practical limitations of grid resolution with the length scales associated with the weakest, and hence most interesting, magnetic fields.
With a toroidal initial field, the disk evolves at a slower rate that CK1. The instability grows over the first few orbits, with the fastest rate of growth associated with the innermost radius. The growing modes of the instability have the same appearance as seen in the local simulations; high m structures appear first, building to lower m with time. The vertical and radial structure also features high wavenumbers. The rapid growth phase of the instability ceases after about 8 orbits at the grid center. At this point the field exhibits a low m, tightly wrapped structure, with rapid variations in R and z. This behavior is consistent with the linear analysis (Balbus & Hawley 1992). Turbulence and accretion begin, and by 12 orbits over half of the disk mass has been lost from the grid (much of it through the outer boundary); what remains is piled up near R=1.7 and is slowly accreting. The magnetic field is subthermal and dominated by the toroidal component. The component magnetic energies at late time are quite similar to those seen in the vertical field run CK1.
Angular momentum transport is again mainly due to the Maxwell stress. Compared with CK1, the toroidal field model has lower overall stress. This is partly because the vertical field model has a very vigorous initial phase associated with the saturation of the linear instability. At the end of both CK1 and CK2, alpham = 0.3-0.4. The mass accretion rate, magnetic field strength, and alpha fluctuate considerably both in time and in space. The total alpha value rises to 0.12 and then declines to 0.06 by the end of CK2. The average ratio of the Maxwell to Reynolds stress is 6 at late times. The value of alpha never rises as high as in local vertical field models, including CK1. It is comparable, however, to the typical value found for weak field toroidal field models in the local shearing box simulations that began with comparable initial toroidal field strengths (HGB).
These Keplerian disk models suffer from significant mass loss through the outer boundary. In contrast, tori can be completely and self-consistently contained initially within the computational domain. We next turn to cylindrical models of constant angular momentum (q=2) tori, beginning with a torus with Rin = 2.0 and Rkep = 2.5. The torus outer boundary is at R=3.3. Two simulations are done, model CT1 with an initial vertical field, and CT2 with an initial toroidal magnetic field; both use the same computational grid as above.
Model CT1 begins with a vertical field with constant beta = 100 from r=2.1 to 3.1. This gives lambdaMRI = 0.25 at the pressure maximum. As the evolution proceeds, the magnetorotational instability rapidly develops, again with the characteristic channel. Early on, the beginnings of the Papaloizou-Pringle principal mode can also be seen as kz = 0 oscillations at the edges of the torus. But long before this global instability can develop, the torus is dramatically altered by the local MRI. The torus does not endure as a torus: it expands rapidly outward as the angular momentum distribution shifts from constant toward Keplerian. After 4 orbits at the initial torus pressure maximum the system has evolved to a nearly Keplerian disk that fills the computational domain.
The magnetic field is amplified to an overall average value of beta=2. Toroidal field dominates: Br2/Bphi2 = 0.1, and Bz2/Br2 = 0.3. As always, the Maxwell stress exceeds the Reynolds stress by a factor of several. The value of alpham peaks at 0.7 and declines to 0.4 at the end of the simulation (t=102). The overall alpha value varies throughout the disk; at t=85 it varies between 0.3 and 0.4.
Next take the same torus and apply an initial toroidal field with beta=100 (Model CT2). With this strength field, the critical azimuthal wavenumber at the pressure maximum is mcrit= 63, so the fastest growing modes are underresolved on this grid. As with the Keplerian simulation CK2, the initial toroidal field model becomes turbulent at a later time compared with an initial vertical field model. The total poloidal magnetic field amplification is also considerably smaller than seen in the vertical field model. Despite this, the qualitative outcome of the instability for the torus is largely the same. The slower onset of the MRI allows the principal mode of the Papaloizou-Pringle instability to appear, but soon the transport of angular momentum brings this to a halt. The disk spreads outward with the bulk of the mass slowly moving inward. The overall angular momentum distribution changes from constant to Keplerian from the inner boundary to R=2.5, and sub-Keplerian but increasing beyond this point. As in the previous simulations, the disk exhibits large local fluctuations in density, stress and field strength, and many low azimuthal, high radial wavenumber features. Similar time and space variations were observed in the local models as well; they are characteristic of the MHD turbulence.
At late time the magnetic energy has risen to beta=15, with the toroidal field dominant by a large factor: Br2/Bphi2= 0.04. The ratio of the vertical to radial field energy is 0.16. The value of alphamrises to ~0.3, close to, but below, the value in the Keplerian toroidal field simulation above. The total stress corresponds to an alpha = 0.02-0.03 at late time. Again the Maxwell stress exceeds the Reynolds component by a factor of 3 to 4. These properties are consistent with toroidal field shearing box simulations (HGB).
Finally, consider a constant angular momentum torus, CT3, with an initial field constructed by setting the azimuthal component of the vector potential equal to the density in the torus, Aphi = rho(R). The resulting field is normalized to an average beta of 100. In the cylindrical limit, the density depends only on radius R and the resulting field is vertical with zero net value integrated over the disk. This model also adopts the pseudo-Newtonian potential sim 1/(R-Rg) with Rg=1. The computation domain runs from R=1.5 to 13.5, phi =0 to pi/2 and z=0 to 2. The grid resolution is 128 × 64 × 32. The torus has an inner edge at r=4.5 and a pressure maximum at rkep = 6.5 where Porb = 88; the orbital period at the outer boundary is 289. A new feature in this model is the presence of the marginally stable orbit at Rms=3. Inside of this radius, matter will rapidly accelerate inward.
The MRI grows rapidly in the inner regions of the disk, again with the characteristic radial channel appearance. Accretion through the inner boundary begins at about t=100. The magnetic energy rises to peak at beta = 8 at t=150. The magnetic energy grows more slowly after that point; additional small peaks are observed due to the growth of the MRI in the outer parts of the disk. Mass loss through the outer boundary begins at t=200, after which the initial linear growth phase is over and the disk is fully turbulent. Between t=200 and the end of the simulation, beta decreases from 10 to 6. As always, the toroidal field exceeds the radial field, with Br2 / Bphi2 = 0.16. The radial field is, in turn, greater than the vertical field, Bz2/Br2 = 0.47.
The overall Maxwell stress in the MHD turbulence reaches a peak at t=150, drops off, and then slowly climbs again. At late time alpham = 0.5, and the total stress parameter alpha varies with radius inside the disk, from 0.1 near the initial pressure maximum, to 0.2 at the outer regions. The angular momentum distribution rapidly evolves away from constant, toward Keplerian. By t=420 it is nearly Keplerian from the pressure maximum inward, and sub-Keplerian out beyond this point. Inside of Rms=3 the Maxwell and Reynolds stresses decline with radius, although not as fast as the pressure. This makes alpha rise sharply.
By the end of the simulation (t=420) over one quarter of the disk mass has been lost, and 70% of it has gone in past the marginally stable orbit. An examination of the various vertically- and azimuthally-averaged velocities at the end of the run (Fig. 2) reveals that the inflow velocity vr accelerates rapidly inside of R=4. By the time the inner grid boundary is encountered, the inflow speed strongly supersonic and super-Alfvenic. Inside the disk the radial speed an order of magnitude smaller than the sound speed, and smaller than the toroidal or poloidal Alfven speeds, which themselves remain subthermal. Going inward from the marginally stable orbit the toroidal Alfven speed rises more slowly that the poloidal Alfven speed. The radial field strength is roughly constant, but the fluid density drops inside of Rms as the flow accelerates inward. The specific angular momentum drops 3% from R=3 to the inner boundary at R=1.5 indicating that there is still some net stress inside Rms.
FIGURE 2: Vertically- and azimuthally-averaged velocities as a function of radius in simulation CT3 at time t=420 (4.8 orbits at the initial pressure maximum). The curves trace the toroidal speed vphi, the adiabatic sound speed cs, the toroidal and poloidal Alfven speeds vA phi and vAp, and the radial speed vr. The vertical dashed line indicates the location of the marginally stable orbit in the pseudo-Newtonian potential. The short-dashed line corresponds to the Keplerian velocity.
Axisymmetry: Simulations in the (R,z) plane
Axisymmetry: Simulations in the (R,z) plane
The cylindrical disk limit represents a useful way of simplifying the full global problem. Another potentially useful simplification is the axisymmetric limit. In this series of simulations, the torus evolution problem is considered in the axisymmetric (R,z) plane. These simulations use a pseudo-Newtonian potential, and begin with a pressure-supported torus that is fully contained on the grid. The initial magnetic field is chosen to satisfy two requirements: it must have a poloidal component to allow for the development of the MRI, and, as a practical matter, it should be contained completely within the torus to avoid the Courant limitations caused by high Alfven speeds due to strong fields in low density regions. A suitable initial configuration consists of magnetic field loops lying along equidensity surfaces in the torus. This initial setup will develop strong toroidal field due to shearing of the initial radial field. However, experience has shown that strong toroidal fields are the outcome of all initial field choices, so this should not represent too great an idealization. All of the two-dimensional simulations considered here will also be run in three dimensions.
The first simulation has a radial grid running from R=1.5 to an outer boundary at R=11.5 and a vertical grid running from z=-5 to 5. The grid resolution is 128 × 128. A constant-l torus is placed on the grid with an inner boundary at the marginally stable orbit, Rin=3 and a pressure maximum at Rkep = 4.7. The orbital period at the pressure maximum is 50. The initial field is obtained from an azimuthal vector potential Aphi = rho(R,z), and is normalized to beta=100 using the total integrated magnetic and thermal energies.
Density grayscale plots from this simulation are presented in Figure 3. The initial period of evolution is dominated by the growth of toroidal magnetic field due to shear. As the magnetic pressure increases, the torus expands, particularly at the inner edge. At first, low density material is driven outward perpendicular to the torus surface; subsequently, it flows radially out and around the torus. Because of the initial reflection symmetry across the equator, the toroidal field changes sign there and a strong current sheet forms. This current sheet proves to be unstable, and oscillates around the equator. This is an important part of the evolution, and indicates that it is necessary to simulate the full (R,z) plane rather than adopt the equator as an explicit boundary.
FIGURE 3: Plots of log density in a two-dimensional axisymmetric simulation of a constant angular momentum torus containing poloidal field loops. Each image is labeled by time. At t=50 the torus has expanded due to shear amplification of the toroidal field. At $t=180$ the poloidal field MRI has set in. A period of turbulence follows (t=250, 350) which dies out by the end of the simulation (t=850, 17 orbits at the initial pressure maxiumum).
As the toroidal field pressure increases, it drives inflow through the marginally stable orbit; this initial accretion flow peaks at t=50 then begins to decline. In the meantime, the poloidal field MRI grows within the torus, and begins to manifest itself visibly in the typical form of radial channels by t=180 (3.6 orbits). There follows a period of violent readjustment within the disk, featuring strong mass inflow punctuated by episodic accretion events. This phase lasts until t~ 350 (7 orbits), beyond which the poloidal magnetic energy declines, and with it the level of turbulence in the disk, and the accretion rate. At the end of the simulation (t=850, 17 orbits) about 60% of the initial disk mass has been lost. Most of this mass is lost by t=500; after this time the inflow accretion rate is very small.
Thus there are three distinct phases to the two-dimensional torus evolution: expansion due to the shear-amplified toroidal magnetic pressure, strong nonlinear evolution of the poloidal field MRI, and finally a more quiescent turbulent state with declining poloidal magnetic field strengths. Angular momentum transport occurs in all three phases at different rates. In the first phase there is a growing Maxwell stress from the amplified Bphi field, mainly in the inner region of the disk where the orbital frequency is highest. The Reynolds stress is negligible. The initially constant angular momentum distribution is unchanged within the disk except in this inner region. With the onset of the MRI, angular momentum transport occurs everywhere, and the specific angular momentum begins to increase with radius. During the middle of this second phase, angular momentum transport peaks, with alpha ranging between 0.2 and 0.5 through the disk. As the turbulence subsides, the rate of transport declines with alpha varying between .02 and .1 within the main part of the disk. Since this is an axisymmetric simulation, the strength of the poloidal magnetic field is limited topologically by the antidynamo theorem, and the MHD turbulence must eventually die out.
The next model is a torus located initially at a larger radius, which increases the amount of time that it can evolve prior to reaching the marginally stable orbit. The grid runs from R=1.5 to R=13.5 and from -4.5 to 4.5 in z. The grid resolution is 120 × 90. The constant-l torus has an inner edge at Rin=4.5 and a pressure maximum at Rkep = 6.5. The orbital period at the pressure maximum is 88. Again, magnetic field loops are placed along the equidensity surfaces within the torus, with a total magnetic energy corresponding to beta = 100.
As before, this torus undergoes three phases of evolution: toroidal field amplification, development of the nonlinear poloidal MRI, and subsequent turbulence. The MRI saturates around t=300 (3.4 orbits) with total magnetic energy beta = 2. The phase of strong MRI turbulence is over by t=450 (5.1 orbits) when the average beta rises to 6; fairly steady accretion follows for the remainder of the simulation which runs to time t=830 (9.5 orbits).
In the third simulation, the initial torus has an angular momentum distribution closer to Keplerian, specifically q=1.68. The torus inner boundary is at 4 and its pressure maximum at Rkep=10 (orbital period =179). The computational domain runs from Rin=1.5 to 21.5, and from -10 to 10 in z; the grid resolution is 128 × 128.
Although it began with a much different initial angular momentum distribution, the evolution is quite similar to the two previous cases. In the early stage, toroidal field is amplified by shear. The poloidal MRI soon comes into play, producing the characteristic radial streams. The total magnetic energy peaks at about t=800 (4.5 orbits). After this the poloidal field energy declines steadily, as does the Maxwell stress and the mass accretion through the inner radial boundary. Toroidal field dominates, with the ratio Br2/Bphi2 = 0.017. Beyond t=2000, the toroidal field energy is essentially constant with a beta= 1.3.
Angular momentum transport begins almost immediately. Initially it is confined to the inner regions of the torus where shear is strongest, and dominated by the Maxwell stress associated with the growing toroidal field. However, the growth of the MRI produces stress throughout the torus. By t=1000 the averaged alpha value is between 0.1 and 0.2. At the end of the simulation at t=2840, alpha ranges between 0.01 and 0.04 within the torus. The ratio of the Maxwell stress to the magnetic pressure, alpham, begins small, rises to a value of 0.1 during the initial saturation of the MRI, then steadily declines to 0.03 by the end of the simulation as poloidal field is preferentially destroyed. Since this torus began with an angular momentum distribution similar to the end state of the initially q=2 tori, there is only a slight change in the slope of the average angular momentum. Once the torus expands, this slope is extended in range all the way from Rmin to the outer radial boundary.
Three Dimensional Tori: Full Global Simulations
Three Dimensional Tori: Full Global Simulations
Now we turn to the evolution of fully three dimensional tori using a pseudo-Newtonian potential. The simulations described in this section are three-dimensional versions of the tori considered above in the axisymmetric limit. The first, model GT1, is the constant angular momentum torus with rkep=4.7 and in=3. The orbital period at the pressure maximum is torb = 50. The computational domain runs from 1.5 to 11.5 in radius, from -5 to 5 in z, and from phi=0 to pi/2. The initial magnetic field is constructed by setting the toroidal component of the vector potential equal to the density inside the disk, Aphi = rho (R,z), for all rho greater than a minimum value. The total magnetic energy is then normalized to a value of beta=100, using the total integrated gas pressure of the torus. Different regions within the torus will, of course, have larger or smaller values of beta. The strongest initial fields are found between the pressure maximum and the inner edge of the torus. The simulation is run to time t=780 (15.6 orbits).
Animation:Click on the image to see an MPEG
animation of run GT1 compared with its 2D axisymmetric counterpart.
The evolution of this torus is illustrated with a series of grayscale plots in log density (Fig. 4) of vertical (R,z) and equatorial (R, phi) slices. As with the axisymmetric torus (Fig.~3), the initial phase of evolution is controlled by the shear amplification of the toroidal field. However, in three-dimensions this toroidal field is itself unstable to the MRI. The characteristic large-wavenumber structures associated with the instability quickly appear in the inner regions of the disk (by t=80). This creates turbulence, and momentarily increases the mass accretion rate over that seen in the axisymmetric simulation (Fig.~5). By t=200 the disk contains tightly wrapped, low m trailing spiral waves. The presence of this nonaxisymmetric structure affects the development of the poloidal field instability, eliminating the organized radial streaming flows seen in axisymmetry. The flow is more turbulent and the mass accretion rate is steadier, without the large impulsive spikes in accretion rate associated with the radial streaming in two dimensions. This overall turbulence declines after saturation of the MRI, but continues through the end of the run. This is the stage when the accretion rate in the two-dimensional simulation drops to very small levels. While in three dimensions one might hope that the disk could achieve a nearly steady state mass accretion rate, there is, of course, only a finite amount of disk mass available to accrete. At the end of the run 85% of the initial total mass has been lost from the grid, most of it through the innermost stable orbit. The amount of mass ejected off the outer boundary is about 30% of that accreted through the inner boundary.
FIGURE 4: Plots of log density in simulation GT1, initially a constant angular momentum torus containing poloidal field loops. Each image consists of a side view and an equatorial view, and is labeled by time. At t=80 the torus has expanded due to shear amplification of the toroidal field which has become visibly unstable in the inner regions. Full turbulence sets in by t=200 (4 orbits at the initial pressure maximum) and continues for the remainder of the simulation. The total disk mass drops steadily due to accretion. These images should be compared with the axisymmetric model in Figure 3.
FIGURE 5: Mass flux through the inner grid radius Rmin as a function of time for global models (solid lines) and their axisymmetric counterparts where applicable (dashed lines). In each case, the mass accretion rate is normalized by the initial mass of the torus. In models GT1, GT2 and GT4, the initial accretion is driven mainly by the growth of the toroidal field due to shear. Subsequently the poloidal MRI drives strong accretion. The coherent channel solution in axisymmetry produces particularly strong accretion events. At late times, accretion in the axisymmetric models dies down while in the global models accretion continues at a reduced, but significant, level. Note that despite their quite different initial field strengths and topologies, models GT2 and GT3 have similar accretion rates at late times.
Although the total volume-averaged magnetic energy drops with time beyond t=400, most of this decline is due to loss from accretion. Average beta after t=200 is ~4 and remains fairly constant thereafter. The magnetic energy is dominated by the toroidal component; the poloidal field beta value has an average value of ~30, and also remains nearly constant after t=200. In the axisymmetric simulation the poloidal beta value attains a value of 15 at t=280, but rises steadily to beta=75 by t=800, indicating a poloidal magnetic energy loss rate exceeding that due to accretion alone.
Angular momentum transport results in rapid restructuring of the disk. The angular momentum distribution in the inner half of the torus steepens to a nearly Keplerian value by t=50 (Fig.~6). The angular momentum in the outer part of the torus increases over the remainder of the simulation. At the end time the angular velocity Omega ~ R-q has a slope near q=1.55. This is essentially the same as in the axisymmetric run above (represented by the short dashed line in Fig.~6). This slope adequately characterizes the disk since the angular momentum, and hence angular velocity, is nearly constant on cylinders. After turbulence has developed, the average alpha values in the main part of the disk are between 0.1 and 0.2. The ratio of the Maxwell to Reynolds stress varies, but lies between 1 and 4 throughout the main part of the disk. The time- and space-averaged value of alpham after t=200 is 0.40.
FIGURE 6: Mass distribution as a function of radius (top), and averaged angular momentum as a function of radius (bottom) at several different times in model GT1. The angular momentum curves are labeled by time; the corresponding mass curves run from top to bottom with increasing time. The orbital time at the initial pressure maximum is 50. The long dashed line corresponds to the Keplerian value in the pseudo-Newtonian potential. The short dashed line is the angular momentum distribution at time t=420 in the axisymmetric version of this torus.
Figure 6 shows that the specific angular momentum continues to decline even inside the marginally stable orbit. This indicates that even here there remains a significant net stress. In fact, inside the marginally stable orbit alpha rapidly increases because the gas pressure drops while the Maxwell stress remains roughly constant. The presence of this stress means there is no maximum in the epicyclic frequency, appa2 = 2 Omega/r dl/dr. In the pseudo-Newtonian potential (and, of course in the relativistic potential that the pseudo-Newtonian potential was designed to imitate), the Keplerian value of kappa2 has a maximum that occurs at about 3.7 rg. More generally, the epicyclic frequency will go to zero at a stress-free inner edge of a disk where dl/dr = 0, ensuring a kappa maximum somewhere in the disk. Here, however, the inner boundary of the disk is not characterized by a zero-stress condition, the slope of dl/dr is nearly constant, and the epicyclic frequency does not turn over but continues to rise. It has been proposed that potentially observable oscillatory modes might be trapped near the disk inner edge where the epicyclic frequency turns over (e.g., Nowak & Wagoner 1993). The present result must be regarded as preliminary, but if significant stress inside the marginally stable orbit proves to be a generic property of magnetic turbulence, trapped disk oscillations may not be present.
In some respects, the two- and three-dimensional simulations are similar. They both have shear amplification of toroidal field, they both evolve due to the resulting increase in toroidal magnetic pressure, and they both develop turbulence due to the poloidal field MRI. Both rapidly evolve from constant to nearly Keplerian specific angular momentum distributions (Fig.~6). The three dimensional simulation, however, permits the development of the nonaxisymmetric MRI which increases and sustains turbulence and mass accretion. In two dimensions the organized poloidal field channel solution produces an impulsive accretion rate greater than that seen in three dimensions during the initial saturation of the poloidal MRI. But axisymmetry causes the two-dimensional poloidal field to decline, and with it the Maxwell stress. The contrast between alpham in the three- and the two-dimensional simulations is instructive (Fig.~7). In three dimensions alpham remains relatively constant and near the value typically seen in the local shearing box simulations. In two dimensions alpham declines steadily with time following the saturation of the MRI. Maintenance of the poloidal field through dynamo action is possible only in three dimensions.
FIGURE 7: Time evolution of volume-integrated value of the magnetic alpha value in the global torus models (solid lines) and the axisymmetric version (dashed lines) where appropriate. Individual graphs are labeled by global torus model number. In three dimensions, poloidal field amplitudes are maintained relative to the toroidal field, and the Maxwell stress remains appreciable compared to the total magnetic pressure. In two dimensions the poloidal fields die out and the stress drops.
Although two-dimensional simulations cannot capture these essential features of global evolution, they do have one clear advantage: they are considerably easier to compute. Two dimensional simulations are useful for searching a wide range of initial conditions in support of the more challenging three-dimensional models. To test the idea of using an evolved two-dimensional simulation as an initial condition, consider next a constant-l torus with an initial inner edge at rin= 4.5 and a pressure maximum at r=6.5 where Porb = 88. This same initial torus was run both in the cylindrical limit, and in the axisymmetric (R,z) limit above. Here the computational domain runs from R=1.5 to 13.5, phi=0 to pi/2 and z=-4.5 to 4.5. The grid resolution is 120 × 64 × 90. The three dimensional torus GT2 is initialized by taking the output from the axisymmetric simulation at time t=233 and expanding it out in the phi direction. This time corresponds to shortly before the nonlinear aspects of the MRI first manifest themselves noticeably in the density plots. Random pressure perturbations are added to break the axisymmetry.
While this initialization procedure reduces the total computational time required, it doesn't allow the toroidal field instability the opportunity to grow during the initial phase. This means that the strong coherent poloidal field instability develops as it does in two dimensions. Nevertheless, significant nonaxisymmetric structure appears by t=400, and the overall turbulence is increased over that seen in axisymmetry. Accretion inflow at the inner boundary is about twice what it is in two dimensions after this point in time. After t=400 the average magnetic field strength is beta=4.7, and the poloidal field strength betap = 24. At late time alpha ~ 0.1, and alpham = 0.37, a level similar to that seen in GT1. In two dimensions the Maxwell stress drops and the poloidal beta increase. A the end of the simulation, alpham = 0.17 and betap = 55. The toroidal field beta is fairly constant with time in both two- and three-dimensions runs after t=450.
This run can also be compared with the cylindrical run CT3. Figure 8 is a plot of the radial dependence of the azimuthally and vertically averaged speeds. These curves are very similar to those in Figure 2 from CT3. GT3 shows a stronger initial field amplification phase and a stronger MRI saturation (created by initializing from the axisymmetric run). However the magnetic field amplitudes and stress levels are comparable at late times in both simulations. These similarities indicate that the cylindrical limit provides a useful approximation for investigating the nonlinear evolution of MHD turbulence near the midplane of a disk.
FIGURE 8: Vertically- and azimuthally-averaged velocities as a function of radius at time t=727 (8.3 orbits at the initial pressure maximum) in simulation GT2. The curves trace the toroidal speed vphi, the adiabatic sound speed cs, the toroidal and poloidal Alfv\'en speeds vAphi and vAp, and the radial speed vr. The vertical dashed line indicates the location of the marginally stable orbit in the pseudo-Newtonian potential. The short-dashed line corresponds to the Keplerian velocity.
To summarize, GT2 shows that fully three dimensional turbulence can be rapidly produced and sustained in a simulation initialized from the output of an axisymmetric simulation. The late-time properties of such a simulation, essentially one with complicated nonlinear initial perturbations, are quite similar to those obtained from a simulation evolved from a formal equilibrium and linear perturbations. The two-dimensional simulations, therefore, serve as just another type of initial condition, albeit more complicated than usually adopted.
One characteristic of all these simulations is that at late times the magnetic field is predominantly toroidal. The next simulation considers a torus that begins with only a toroidal field. Model GT3 is of a q=2 torus with an inner boundary at 4.5 and a pressure maximum at Rkep=6.5 (orbital period =88). The initial toroidal field has beta=4 everywhere within the torus. The MRI that will result is nonaxisymmetric; an axisymmetric (R,z) version of this model does not develop turbulence or transport angular momentum. The computational domain runs from Rin=1.5 to 13.5, z from -6 to 6, and over the full 2 pi in angle; the grid resolution is 128 × 128 × 128.
No attempt is made to keep the initial torus in pressure balance; the initial magnetic field simply adds to the equilibrium hydrodynamic pressure. As a consequence, the torus undergoes an axisymmetric readjustment due to this additional magnetic pressure. The toroidal field MRI develops rapidly at the inner edge of the disk where Omega is largest, before spreading throughout the disk. The poloidal field energy grows steadily with time until about t=250 when it reaches a value of betap= 34. After t=400, betap increases with time to 65 at the end of the run. Grayscale images of GT3 are given in Figure 9.
Animation:Click on the image to see an MPEG
animation of run GT3.
FIGURE 9: Plots of log density in simulation GT3, a torus with q=2 and beta=4 toroidal field initially. Each frame consists of a side and an equatorial view, and is labeled by time. At t=120 the MRI has appeared at the inner edge of the torus. Turbulence is soon established throughout the torus.
Models GT2 and GT3 began with the same hydrodynamic torus; they differ in their initial magnetic fields, in the size of the phi domain, and in resolution. Despite these differences, the late time evolution in these two runs is very similar. They have comparable mass accretion rates through Rmin (Fig. 5). Both have alpham = 0.36 at t=720. Figure 10 illustrates the redistribution of mass and angular momentum as a function of time. The slope of the angular momentum distribution at late time is the same in GT3 as in both GT1 and GT2.
FIGURE 10: Time evolution of the vertically- and azimuthally-averaged mass (top) and specific angular momentum (bottom) as a function of radius in simulation GT3. The different lines correspond to different times throughout the evolution: $t=0$, 117, 209, 298, 382, 470, 550, 637, 727. The orbital period at the initial torus pressure maximum is 88. The long dashed line in the angular moment plot corresponds to the Keplerian value. The short dashed line is the curve from the final time in model GT1.
Figure 11 shows the radial run of density Sigma, velocity vR, and mass flux at time t=720. Density is normalized by the maximum Sigma at t=0. Within the torus the radial drift velocity is small compared to the orbital velocity; vR reverses outside of R=8 beyond which there is net outflow. Inside Rms the density decreases rapidly as the inflow accelerates inward. The accretion rate increases slightly from R=7 inward. Figure 12 shows gas and magnetic pressures at this same end time (normalized by the initial pressure maximum value), as well as the Maxwell stress. The magnetic energy remains subthermal throughout, although it drops less rapidly than the gas pressure inside of Rms. The ratio of the Maxwell stress to the gas pressure has a minimum of 0.03 between R=4 and 5, rising sharply inside Rms and more slowly toward the outer boundary.
FIGURE 11: From top to bottom, the vertically averaged density, radial velocity, and mass flux as a function of radius in model GT3 at t=720 (8.2 orbits at the initial pressure maximum). The dashed line in the radial velocity plot is the orbital velocity.
FIGURE 12: Vertically averaged gas pressure P, toroidal magnetic field pressure, radial field pressure, and Maxwell stress as a function of radius at t=720 in model GT3. The pressures are normalized to the initial pressure maximum value.
The final global simulation, GT4, is the q=1.68 torus with an inner boundary at 4 and a pressure maximum at Rkep=10. This angular momentum distribution yields a torus that extends over a large radial distance without becoming too thick in the vertical direction. The computational domain runs from Rin=1.5 to 21.5, from -10 to 10 in z, and over the full 2 pi in angle; the grid resolution is 128 × 128 × 128. The simulation is run to time t=1280. This is longer than the other simulations, although it amounts to fewer orbits at the pressure maximum: t=1280 is 7.2 orbits at Rkep. This evolution time is sufficient to observe all of the stages seen in the other three-dimensional simulations, although not long enough for the torus to have settled into a quasi-steady state.
Animation:Click on the image to see an MPEG
animation of run GT4.
Density plots from GT4 are presented in Figure 13. At the beginning, the toroidal field grows by shearing out the radial field, but as it does so, the toroidal MRI sets in. This soon leads to turbulence. The poloidal MRI develops rapidly in the inner regions of the disk and subsequently spreads throughout. The inner edge of the disk moves slowly inward until it reaches the marginally stable orbit at t=145 after which gas plunges inward. Initially the inflowing fluid is confined largely to the equatorial plane. As time passes, however, this region fills with gas and becomes thicker.
FIGURE 13: Plots of log density in simulation GT4, initially a torus with q=1.68 and poloidal field loops. Each frame consists of a side and an equatorial view, and is labeled by time. The orbital period at the initial pressure maximum is 179. At t=320 the torus has expanded due to shear amplification of the toroidal field which has, in turn, become unstable with m=4 dominating in the equatorial plane. Beyond this time the poloidal field MRI begins and turbulence follows.
As with the previous global simulations, the mass accretion rate is larger in three dimensions compared to two. Near the end of the simulation, the accretion rate at the inner radial boundary is approximately 2.5 times that seen in axisymmetry. The value of alpham is rising with time indicating that the poloidal field strength is still increasing at the end of the simulation.
Figure 14 shows the radial mass and angular momentum distributions at the initial and final times in GT4. Also shown are curves from two times in the equivalent axisymmetric calculation. The average angular velocity parameter q decreases by a very small amount over the course of the evolution. The torus mass, on the other hand, has been substantially redistributed.
FIGURE 14: Vertically and azimuthally averaged mass (top) and specific angular momentum (bottom) as a function of radius for the initial and final time, t=1283 in model GT4. The slope has steepened very slightly over the course of the evolution. The long dashed line in the angular momentum plot corresponds to the Keplerian value. The short dashed lines are from times t=1559 and t=2840 in the axisymmetric simulation of this torus.
Back to title page