We develop a model of the white dwarf (WD) - red dwarf (RD) binaries AR Sco and AE Aqr as systems in a transient propeller stage of highly asynchronous intermediate polars. The WDs are relatively weakly magnetized with magnetic field of $\sim 10^6$ G. We explain the salient observed features of the systems due to the magnetospheric interaction of two stars. Currently, the WD's spin-down is determined by the mass loading of the WD's magnetosphere from the RD's at a mild rate of $\dot{M}_{WD} \sim 10^{-11} M_\odot $/yr. Typical loading distance is determined by the ionization of the RD's wind by the WD's UV flux. The WD was previously spun up by a period of high accretion rate from the RD via Roch lobe overflow with $\dot{M} \sim 10^{-9} M_\odot $/yr, acting for as short a period as tens of thousands of years. The non-thermal X-ray and optical synchrotron emitting particles originate in reconnection events in the magnetosphere of the WD due to the interaction with the flow from the RD. In the case of AR Sco, the reconnection events produce signals at the WD's rotation and beat periods - this modulation is due to the changing relative orientation of the companions' magnetic moments and resulting variable reconnection conditions. Radio emission is produced in the magnetosphere of the RD, we hypothesize, in a way that it is physically similar to the Io-induced Jovian decametric radiation.

Many important problems in astrophysics, space physics, and geophysics involve flows of (possibly ionized) gases in the vicinity of a spherical object, such as a star or planet. The geometry of such a system naturally favors numerical schemes based on a spherical mesh. Despite its orthogonality property, the polar (latitude-longitude) mesh is ill suited for computation because of the singularity on the polar axis, leading to a highly non-uniform distribution of zone sizes. The consequences are (a) loss of accuracy due to large variations in zone aspect ratios, and (b) poor computational efficiency from a severe limitations on the time stepping. Geodesic meshes, based on a central projection using a Platonic solid as a template, solve the anisotropy problem, but increase the complexity of the resulting computer code. We describe a new finite volume implementation of Euler and MHD systems of equations on a triangular geodesic mesh (TGM) that is accurate up to fourth order in space and time and conserves the divergence of magnetic field to machine precision. The paper discusses in detail the generation of a TGM, the domain decomposition techniques, three-dimensional conservative reconstruction, and time stepping.

Magnetic fields advected along with low angular momentum accretion flows predominantly become toroidal due to the strong azimuthal velocity close to a black hole. We study self-consistently the movements of these flux tubes inside an advective disc and how they dynamically influence the flow. We find that the centrifugal barrier slows down the radial motion of the flux tubes. In this case, the large magnetic flux tubes with a significant drag force escape along the vertical axis due to buoyancy. Magnetic pressure rises close to the black hole and together with the centrifugal force, it combats gravity. The tug-of-war among these forces causes the centrifugal pressure supported shock to oscillate radially. We study the effects of successive injection of flux tubes and find how the flux tube could be trapped inside the disc in regions of highest entropy. Most interestingly, the shock wave remains at its average location and is not destroyed. We show that the toroidal field loops contribute significantly to collimate and accelerate the outflows from the centrifugal barrier and suggest this mechanism to be a way to collimate and accelerate jets.

There is a great need in several areas of astrophysics and space-physics to carry out high order of accuracy, divergence-free MHD simulations on spherical meshes. This requires us to pay careful attention to the interplay between mesh quality and numerical algorithms. Methods have been designed that fundamentally integrate high order isoparametric mappings with the other high accuracy algorithms that are needed for divergence-free MHD simulations on geodesic meshes. The goal of this paper is to document such algorithms that are implemented in the geodesic mesh version of the RIEMANN code. The fluid variables are reconstructed using a special kind of WENO-AO algorithm that integrates the mesh geometry into the reconstruction process from the ground-up. A novel divergence-free reconstruction strategy for the magnetic field that performs efficiently at all orders, even on isoparametrically mapped meshes, is then presented. The MHD equations are evolved in space and time using a novel ADER predictor algorithm that is efficiently adapted to the isoparametrically mapped geometry. The application of one-dimensional and multidimensional Riemann solvers at suitable locations on the mesh then provides the corrector step. The corrector step for the magnetic field uses a Yee-type staggering of magnetic fields. This results in a scheme with divergence-free update for the magnetic field. The use of ADER enables a one-step update which only requires one messaging operation per complete timestep. This is very beneficial for parallel processing. Several accuracy tests are presented as are stringent test problems. PetaScale performance is also demonstrated on the largest available supercomputers.

We obtained high temporal resolution spectroscopy of the unusual binary system AR Sco covering nearly an orbit. The H$\alpha$ emission shows a complex line structure similar to that seen in some polars during quiescence. Such emission is thought to be due to long-lived prominences originating on the red dwarf. A difference between AR Sco and these other systems is that the white dwarf in AR Sco is rapidly spinning relative to the orbital period. "Slingshot" prominences stable at 3 to 5 stellar radii require surface magnetic fields between 100 and 500 G. This is comparable to the estimated WD magnetic field strength near the surface of the secondary. Our time-resolved spectra also show emission fluxes, line equivalent widths, and continuum color varying over the orbit and the beat/spin periods of the system. During much of the orbit, the optical spectral variations are consistent with synchrotron emission with the highest energy electrons cooling between pulses. On the time-scale of the beat/spin period we detect red and blue-shifted H$\alpha$ emission flashes that reach velocities of 700 km/s. Red-shifted Balmer emission flashes are correlated with the bright phases of the continuum beat pulses while blue-shifted flashes appear to prefer the time of minimum in the beat light curve. We propose that much of the energy generated in AR Sco comes from fast magnetic reconnection events occurring near the inward face of the secondary and we show that the energy generated by magnetic reconnection can account for the observed excess luminosity from the system.

The small-scale turbulent dynamo is an important process contributing to the cosmic magnetization. In partially ionized astrophysical plasmas, the dynamo growth of magnetic energy strongly depends on the coupling state between ions and neutrals and the ion-neutral collisional damping effect. A new damping stage of turbulent dynamo in a weakly ionized medium was theoretically predicted by Xu \& Lazarian (2016). By carrying out a 3D two-fluid dynamo simulation, here we for the first time numerically confirmed the physical conditions and the linear-in-time growth of magnetic field strength of the damping stage of dynamo. The dynamo-amplified magnetic field has a characteristic length as the damping scale, which increases with time and can reach the injection scale of turbulence after around eight largest eddy-turnover times given sufficiently low ionization fraction and weak initial magnetic field. Due to the weak coupling between ions and neutrals, most turbulent energy carried by neutrals cannot be converted to the magnetic energy, resulting in a relatively weak magnetic field at the end of dynamo. This result has important implications for the growth of magnetic fields in the partially ionized interstellar medium and shock acceleration of Galactic cosmic rays.

We study time evolution of sub-Keplerian transonic accretion flows onto black holes using a general relativistic numerical simulation code. We perform simulations around the black holes having non-zero rotation. We first compare one-dimensional simulation results with theoretical results and validate the performance of our code. Next, we present results of axisymmetric, two-dimensional simulation of advective flows. In the literature, there is no solution which describes steady shock solutions in two dimensions. However, our simulations produce these centrifugal force supported steady shock waves even in presence of strong dragging of inertial frames. Since the post-shock region could be hot and upscatter photons through Comptonization, these shock would put imprints on the spectra. Thus, our solutions, which represent truly new results, could be useful to measure spins through radiation spectrum of accreting Kerr black holes.

We study time evolution of sub-Keplerian transonic accretion flows onto black holes using a general relativistic numerical simulation code. We perform simulations in Schwarzschild spacetime. We first compare one-dimensional simulation results with theoretical results and validate the performance of our code. Next, we present results of axisymmetric, two-dimensional simulation of advective flows. We find that even in this case, for which no complete theoretical analysis is present in the literature, steady state shock formation is possible.

As computational astrophysics comes under pressure to become a precision science, there is an increasing need to move to high accuracy schemes for computational astrophysics. Hence the need for a specialized review on higher order schemes for computational astrophysics. The focus here is on weighted essentially non-oscillatory (WENO) schemes, discontinuous Galerkin (DG) schemes and PNPM schemes. WENO schemes are higher order extensions of traditional second order finite volume schemes which are already familiar to most computational astrophysicists. DG schemes, on the other hand, evolve all the moments of the solution, with the result that they are more accurate than WENO schemes. PNPM schemes occupy a compromise position between WENO and PNPM schemes. They evolve an Nth order spatial polynomial, while reconstructing higher order terms up to Mth order. As a result, the timestep can be larger. Time-dependent astrophysical codes need to be accurate in space and time. This is realized with the help of SSP-RK (strong stability preserving Runge-Kutta) schemes and ADER (Arbitrary DERivative in space and time) schemes. The most popular approaches to SSP-RK and ADER schemes are also described. The style of this review is to assume that readers have a basic understanding of hyperbolic systems and one-dimensional Riemann solvers. Such an understanding can be acquired from a sequence of prepackaged lectures available from http://www.nd.edu/~dbalsara/Numerical-PDE-Course. We then build on this understanding to give the reader a practical introduction to the schemes described here. The emphasis is on computer-implementable ideas, not necessarily on the underlying theory, because it was felt that this would be most interesting to most computational astrophysicists.

In our prior papers, we considered the non-relativistic linear stability analysis of magnetized jets that do not have current sheet at the boundary. In this paper, we extend our analysis to relativistic jets. In order to find the unstable modes of current sheet-free, magnetized relativistic jets, we linearize full relativistic magnetohydrodynamics equations and solve them numerically. We find the dispersion relation of the pinch and kink mode instabilities. By comparing the dispersion relations of mildly relativistic jet (Lorentz factor 2) with moderately relativistic jet (Lorentz factor 10), we find that the jet with higher Lorentz factor is significantly more stable in both pinch and kink modes. We show that inclusion of the current sheet-free magnetic field in the jet further enhances the stability. Both pinch and kink mode instabilities become progressively more stable with increasing magnetization. We also show a scaling relation between the maximum temporal growth rate of the unstable mode and the Lorentz factor of the jet. The maximum temporal growth rates of the unstable modes are inversely proportion to the Lorentz factors for most of the modes that we study. However, for the fundamental pinch mode it is inversely proportional to the square of the Lorentz factor. This very beneficial scaling relation holds regardless of the presence of a magnetic field.

In the magnetosphere of a rotating black hole, an inner Alfven critical surface (IACS) must be crossed by inflowing plasma. Inside the IACS, Alfven waves are inward directed toward the black hole. The majority of the proper volume of the active region of spacetime (the ergosphere) is inside of the IACS. The charge and the totally transverse momentum flux (the momentum flux transverse to both the wave normal and the unperturbed magnetic field) are both determined exclusively by the Alfven polarization. Thus, it is important for numerical simulations of black hole magnetospheres to minimize the dissipation of Alfven waves. Elements of the dissipated wave emerge in adjacent cells regardless of the IACS, there is no mechanism to prevent Alfvenic information from crossing outward. Thus, numerical dissipation can affect how simulated magnetospheres attain the substantial Goldreich-Julian charge density associated with the rotating magnetic field. In order to help minimize dissipation of Alfven waves in relativistic numerical simulations we have formulated a one-dimensional Riemann solver, called HLLI, which incorporates the Alfven discontinuity and the contact discontinuity. We have also formulated a multidimensional Riemann solver, called MuSIC, that enables low dissipation propagation of Alfven waves in multiple dimensions. The importance of higher order schemes in lowering the numerical dissipation of Alfven waves is also catalogued.

In various astrophysics settings it is common to have a two-fluid relativistic plasma that interacts with the electromagnetic field. While it is common to ignore the displacement current in the ideal, classical magnetohydrodynamic limit, when the flows become relativistic this approximation is less than absolutely well-justified. In such a situation, it is more natural to consider a positively charged fluid made up of positrons or protons interacting with a negatively charged fluid made up of electrons. The two fluids interact collectively with the full set of Maxwell's equations. As a result, a solution strategy for that coupled system of equations is sought and found here. Our strategy extends to higher orders, providing increasing accuracy. Three important innovations are reported here. In our first innovation, the magnetic field within each zone is reconstructed in a divergence-free fashion while the electric field within each zone is reconstructed in a form that is consistent with Gauss' law. In our second innovation, a multidimensionally upwinded strategy is presented which ensures that the magnetic field can be updated via a discrete interpretation of Faraday's law and the electric field can be updated via a discrete interpretation of the generalized Ampere's law. Our third innovation consists of an efficient design for several popular IMEX schemes so that they provide strong coupling between the finite-volume-based fluid solver and the electromagnetic fields at high order. Several accuracy analyses are presented showing that our method meets its design accuracy in the MHD limit as well as in the limit of electromagnetic wave propagation. Several stringent test problems are also presented. We also present a relativistic version of the GEM problem, which shows that our algorithm can successfully adapt to challenging problems in high energy astrophysics.

In a prior paper (Kim et al. 2015) we considered the linear stability of magnetized jets that carry no net electric current and do not have current sheets. In this paper, in addition to physically well-motivated magnetic field structures, we also include the effects of jet shear. The jets we study have finite thermal pressure in addition to having realistic magnetic field structures and velocity shear. We find that shear has a strongly stabilizing effect on various modes of jet instability. Increasing shear stabilizes the fundamental pinch modes at long wavelengths and short wavelengths. Increasing shear also stabilizes the first reflection pinch modes at short wavelengths. Increasing shear has only a very modest stabilizing effect on the fundamental kink modes at long wavelengths; however, increasing shear does have a strong stabilizing effect on the fundamental kink modes at short wavelengths. The first reflection kink modes are strongly stabilized by increasing shear at shorter wavelengths. Overall, we find that the combined effect of magnetic field and shear stabilizes jets more than shear alone. In addition to the results from a formal linear stability analysis, we present a novel way of visualizing and understanding jet stability. This gives us a deeper understanding of the enhanced stability of sheared, magnetized jets. We also emphasize the value of our numerical approach in understanding the linear stability of jets with realistic structure.

The relativistic magnetohydrodynamics (RMHD) set of equations has recently seen increased use in astrophysical computations. Even so, RMHD codes remain fragile. The reconstruction can sometimes yield superluminal velocities in certain parts of the mesh. In this paper we present a reconstruction strategy that overcomes this problem by making a single conservative to primitive transformation per cell followed by higher order WENO reconstruction on a carefully chosen set of primitives that guarantee subluminal reconstruction of the flow variables. For temporal evolution via a predictor step we also present second, third and fourth order accurate ADER methods that keep the velocity subluminal during the predictor step. The RMHD system also requires the magnetic field to be evolved in a divergence-free fashion. In the treatment of classical numerical MHD the analogous issue has seen much recent progress with the advent of multidimensional Riemann solvers. By developing multidimensional Riemann solvers for RMHD, we show that similar advances extend to RMHD. As a result, the face-centered magnetic fields can be evolved much more accurately using the edge-centered electric fields in the corrector step. Those edge-centered electric fields come from a multidimensional Riemann solver for RMHD which we present in this paper. In this paper we also develop several new test problems for RMHD. We show that RMHD vortices can be designed that propagate on the computational mesh as self-preserving structures. These RMHD vortex test problems provide a means to do truly multidimensional accuracy testing for RMHD codes. Several other stringent test problems are presented. We show the importance of resolution in certain test problems. Our tests include a demonstration that RMHD vortices are stable when they interact with shocks.

We present the first high order one-step ADER-WENO finite volume scheme with Adaptive Mesh Refinement (AMR) in multiple space dimensions. High order spatial accuracy is obtained through a WENO reconstruction, while a high order one-step time discretization is achieved using a local space-time discontinuous Galerkin predictor method. Due to the one-step nature of the underlying scheme, the resulting algorithm is particularly well suited for an AMR strategy on space-time adaptive meshes, i.e.with time-accurate local time stepping. The AMR property has been implemented 'cell-by-cell', with a standard tree-type algorithm, while the scheme has been parallelized via the Message Passing Interface (MPI) paradigm. The new scheme has been tested over a wide range of examples for nonlinear systems of hyperbolic conservation laws, including the classical Euler equations of compressible gas dynamics and the equations of magnetohydrodynamics (MHD). High order in space and time have been confirmed via a numerical convergence study and a detailed analysis of the computational speed-up with respect to highly refined uniform meshes is also presented. We also show test problems where the presented high order AMR scheme behaves clearly better than traditional second order AMR methods. The proposed scheme that combines for the first time high order ADER methods with space--time adaptive grids in two and three space dimensions is likely to become a useful tool in several fields of computational physics, applied mathematics and mechanics.

In this paper we consider stability of magnetized jets that carry no net electric current and do not have current sheets (Gourgouliatos et al. 2012). The non-relativistic MHD equations are linearized around the background velocity and the magnetic field structure of the jet. The resulting linear equations are solved numerically inside the jet. We find that introduction of current-sheet-free magnetic field significantly improves jet stability relative to unmagnetized jets or magnetized jets with current sheets at their surface. This particularly applies to the fundamental pinch and kink modes - they become completely suppressed in a wide range of long wavelengths that are known to become most pernicious to jet stability when the evolution enters the non-linear regime. The reflection modes, both for the pinch and kink instability, also become progressively more stable with increased magnetization.

We investigate the nature of the Alfv\'enic turbulence cascade in two fluid MHD simulations in order to determine if turbulence is damped once the ion and neutral species become decoupled at a critical scale called the ambipolar diffusion scale (L$_{AD}$). Using mode decomposition to separate the three classical MHD modes, we study the second order structure functions of the Alfv\'en mode velocity field of both neutrals and ions in the reference frame of the local magnetic field. On scales greater than L$_{AD}$ we confirm that two fluid turbulence strongly resembles single fluid MHD turbulence. Our simulations show that the behavior of two fluid turbulence becomes more complex on scales less than L$_{AD}$. We find that Alfvenic turbulence can exist past L$_{AD}$ when the turbulence is globally super-Alfv\'enic, with the ions and neutrals forming separate cascades once decoupling has taken place. When turbulence is globally sub-Alfvenic and hence strongly anisotropic with a large separation between the parallel and perpendicular decoupling scales, turbulence is damped at L$_{AD}$. We also find that the power spectrum of the kinetic energy in the damped regime is consistent with a $k^{-4}$ scaling (in agreement with the predictions of Lazarian, Vishniac & Cho 2004).

A high order one-step ADER-WENO finite volume scheme with Adaptive Mesh Refinement (AMR) in multiple space dimensions is presented. A high order one-step time discretization is achieved using a local space-time discontinuous Galerkin predictor method, while a high order spatial accuracy is obtained through a WENO reconstruction. Thanks to the one-step nature of the underlying scheme, the resulting algorithm can be efficiently imported within an AMR framework on space-time adaptive meshes. We provide convincing evidence that the presented high order AMR scheme behaves better than traditional second order AMR methods. Tests are shown of the new scheme for nonlinear systems of hyperbolic conservation laws, including the classical Euler equations and the equations of ideal magnetohydrodynamics. The proposed scheme is likely to become a useful tool in several astrophysical scenarios.

We present high resolution simulations of two-fluid (ion-neutral) MHD turbulence with resolutions as large as 512^3. The simulations are supersonic and mildly sub-Alfvenic, in keeping with the conditions present in molecular clouds. Such turbulence is thought to influence star formation processes in molecular clouds because typical cores form on length scales that are comparable to the dissipation scales of this turbulence in the ions. The simulations are motivated by the fact that recent studies of isophotologue lines in molecular clouds have found significant differences in the linewidth-size relationship for neutral and ion species. The goals of this paper are to explain those observations using simulations and analytic theory, present a new set of density-based diagnostics by drawing on similar diagnostics that have been obtained by studying single-fluid turbulence, and show that our two-fluid simulations play a vital role in reconciling alternative models of star formation. The velocity-dependent diagnostics display a complementarity with the density-dependent diagnostics. We find that the linewidth-size relationships show a prominent difference between ions and neutrals when the line of sight is orthogonal to the mean field. This is because the MHD waves in the ions differ from hydrodynamic waves in the neutrals. We also find that the density probability distribution functions (PDFs) show prominent differences between the ions and neutrals when the line of sight is parallel to the mean field. This is because the velocity fluctuations in the ions tend to produce column density fluctuations along field lines. When the magnetic field makes an angle to the line of sight, both observable differences should be visible. These diagnostics should be easy for observers to test. This analysis assumes optically thin lines and a mean magnetic field that is uniform in direction within a cloud.

This report describes a new magnetohydrodynamic numerical model based on a hexagonal spherical geodesic grid. The model is designed to simulate astrophysical flows of partially ionized plasmas around a central compact object, such as a star or a planet with a magnetic field. The geodesic grid, produced by a recursive subdivision of a base platonic solid (an icosahedron), is free from control volume singularities inherent in spherical polar grids. Multiple populations of plasma and neutral particles, coupled via charge-exchange interactions, can be simulated simultaneously with this model. Our numerical scheme uses piecewise linear reconstruction on a surface of a sphere in a local two-dimensional "Cartesian" frame. The code employs HLL-type approximate Riemann solvers and includes facilities to control the divergence of magnetic field and maintain pressure positivity. Several test solutions are discussed, including a problem of an interaction between the solar wind and the local interstellar medium, and a simulation of Earth's magnetosphere.

Very high energy gamma-ray flares from the Crab nebular detected by AGILE and Fermi satellites challenge our understanding of the pulsar wind nebulae. The short duration of the flares, only few days, is particularly puzzling since it is much shorter than the dynamical times scale of the nebular. In this work we investigate analytically and via numerical simulations the electromagnetic signatures expected from the large amplitude low frequency magnetosonic waves generated within the Crab nebular which induce the corrugation perturbations of the termination shock. As a result, the oblique termination shock produces time-dependent, mildly relativistic post-shock flow. Using the relativistic MHD version of the RIEMANN code, we simulate the interaction of the termination shock with downstream perturbations. We demonstrate that mild Doppler boosting of the synchrotron emission in the post-shock flow can produce bright, short time scale flares.

We have studied saturated, MRI-driven turbulence using three-dimensional, isothermal simulations with resolutions that extend from 64 to 192 zones in each direction. The simulations were performed with several higher order Godunov algorithms. A variety of reconstruction strategies as well as a variety of Riemann solvers are tried. We show that the details of the isothermal MRI-driven turbulence depend principally on the Riemann solver and secondarily on the reconstruction strategy. Furthermore, we find that the effective viscosity parameter parameter tends to show progressively smaller decrements with increasing resolution when the best reconstruction strategy (WENO) and the best Riemann solver (linearized)are used. We attribute this result to the more sophisticated dissipation mechanisms that are used in higher-order Godunov schemes. Spectral analysis and transfer functions have been used to quantify the dissipative processes in these higher-order Godunov schemes.

We present a combination of theoretical and simulation-based examinations of the role of two-fluid ambipolar drift on molecular line widths. The dissipation provided by ion-neutral interactions can produce a significant difference between the widths of neutral molecules and the widths of ionic species, comparable to the sound speed. We demonstrate that Alfven waves and certain families of magnetosonic waves become strongly damped on scales comparable to the ambipolar diffusion scale. Using the RIEMANN code, we simulate two-fluid turbulence with ionization fractions ranging from 10^{-2} to 10^{-6}. We show that the wave damping causes the power spectrum of the ion velocity to drop below that of the neutral velocity when measured on a relative basis. Following a set of motivational observations by Li & Houde (2008), we produce synthetic line width-size relations that shows a difference between the ion and neutral line widths, illustrating that two-fluid effects can have an observationally detectable role in modifying the MHD turbulence in the clouds.

We present a series of simulations of turbulent stratified protostellar discs with the goal of characterizing the settling of dust throughout a minimum-mass solar nebula. We compare the evolution of both compact spherical grains, as well as highly fractal grains. Our simulations use a shearing-box formulation to study the evolution of dust grains locally within the disc, and collectively our simulations span the entire extent of a typical accretion disc. The dust is stirred by gas that undergoes MRI-driven turbulence. This establishes a steady state scale height for the dust that is different for dust of different sizes. This sedimentation of dust is an important first step in planet formation and we predict that ALMA should be able to observationally verify its existence. When significant sedimentation occurs, the dust will participate in a streaming instability that significantly enhances the dust density. We show that the streaming instability is pervasive in the outer disc. We characterize the scale heights of dust whose size ranges from a few microns to a few centimeters. We find that for spherical grains, a power-law relationship develops for the scale height with grain size, with a slope that is slightly steeper than -1/2. The sedimentation is strongest in the outer disc and increases for large grains. The results presented here show that direct measurements of grain settling can be made by ALMA and we present favorable conditions for observability. The streaming instability should also be directly observable and we provide conditions for directly observing it. We calculate collision rates and growth rates for the dust grains in our simulations of various sizes colliding with other grains, and find that these rates are significantly enhanced through the density enhancement arising from the streaming instability.

Using a 2.5D time-dependent axisymmetric numerical code we recently developed, we solve the full compressible Navier-Stokes equations (including an alpha-viscosity prescription) to determine the structure of the boundary layer between the white dwarf and the accretion disk in non-magnetic cataclysmic varia ble systems. In this preliminary work, our numerical approach does not include radiation. In the energy equation, we either take the dissipation function (Phi) into account or we assumed that the energy is instantly radiated away (Phi). For a slowly rotating non magnetized accreting white dwarf, the accretion disk e xtends all the way to the stellar surface. There, the matter impacts and spread s towards the poles as new matter continuously piles up behind it. We carried out numerical simulations for different values of the alpha viscosity parameter (alpha), corresponding to different mass accretion rates. In the high viscosity cases (alpha=0.1), the spreading boundary layer sets off a gravity wave in the s urface matter. The accretion flow moves supersonically over the cusp making it s usceptible to the rapid development of gravity wave and/or Kelvin-Helmholtz shea ring instabilities. This BL is optically thick and extends more than 30 degrees to either side of the disk plane after only 3/4 of a Keplerian rotation period (19s). In the low viscosity cases (alpha=0.001), the spreading boundary layer does not set off gravity waves and it is optically thin.

(Abridged) In this paper we have used the RIEMANN code for computational astrophysics to study the interaction of a realistic distribution of dust grains with gas in a vertically stratified protostellar accretion disc. The disc was modeled to have the density and temperature of a minimum mass solar nebula, and was driven to a fully-developed turbulence via the magnetorotational instability (MRI). We find that the inclusion of standard dust to gas ratios does not have any significant effect on the MRI even when the dust sediments to the midplane of the accretion disc. The density distribution of the dust reaches a Gaussian profile, and the scale heights for the dust that we derive are shown to be proportional to the reciprocal of the square root of the dust radius. The largest dust shows a strong tendency to settle to the midplane of the accretion disc, and tends to organize itself into elongated clumps of high density. The dynamics of these clumps is shown to be consistent with a streaming instability. The streaming instability is seen to be very vigorous and persistent once it forms. Each stream of high density dust displays a reduced RMS velocity dispersion, and the densest clumpings of large dust are shown to form where the streams intersect. We have also shown that the mean free path and collision time for the dust that participates in the streaming instability is reduced by almost two orders of magnitude relative to the average mean free paths and collision times. We show that some of the large dust in our 10 au simulations should have a propensity for grain coalescence.

We present a semi-implicit method for isothermal two-fluid ion-neutral ambipolar drift that is second-order accurate in space and time. The method has been implemented in the RIEMANN code for astrophysical fluid dynamics. We present four test problems that show the method works and correctly tracks the propagation of MHD waves and the structure of two-fluid C-shocks. The accurate propagation of MHD waves in the two-fluid approximation is shown to be a stringent test of the algorithm. We demonstrate that highly accurate methods are required in order to properly capture the MHD wave behaviour in the presence of ion-neutral friction.

We present a large number of two and a half dimensional simulations of supernova remnants expanding into interstellar media having a range of densities, temperatures and magnetic field strengths. The volume of hot gas produced is strongly dependent on the inclusion of thermal conduction and magnetic fields. The four-volumes and three-areas of hot gas have been catalogued and their dependence on interstellar parameters documented. Simulated line widths of radioactive species ejected by supernovae have also been catalogued.

Anisotropic thermal conduction plays an important role in various astrophysical systems. One of the most stringent tests of thermal conduction can be found in supernova remnants. In this paper we study anisotropic thermal conduction and examine the physical nature of the flux of thermal conduction in the classical and saturated limits. We also present a temporally second-order accurate implicit-explicit scheme for the time-update of thermal conduction terms within a numerical MHD scheme. Several simulations of supernova remnants are presented for a range of ISM parameters. The role of thermal conduction in such remnants has been studied. We find that thermal conduction produces cooler temperatures and higher densities in the hot gas bubbles that form in the remnants. The effect of thermal conduction in changing the thermal characteristics of the hot gas bubble increases as the remnant propagates through denser ISMs. Remnants evolving in denser ISMs are shown to make a faster transition to a centre-bright x-ray morphology, with the trend emerging earlier in hard x-rays than in the soft x-rays.

We present new observations of the fundamental ro-vibrational CO spectrum of V1647 Ori, the young star whose recent outburst illuminated McNeil's Nebula. Previous spectra, acquired during outburst in 2004 February and July, had shown the CO emission lines to be broad and centrally peaked-similar to the CO spectrum of a typical classical T Tauri star. In this paper, we present CO spectra acquired shortly after the luminosity of the source returned to its pre-outburst level (2006 February) and roughly one year later (2006 December and 2007 February). The spectrum taken in 2006 February revealed blue-shifted CO absorption lines superimposed on the previously observed CO emission lines. The projected velocity, column density, and temperature of this outflowing gas was 30 km/s, 3^{+2}_{-1}E18 cm^{-2$, and 700^{+300}_{-100} K, respectively. The absorption lines were not observed in the 2006 December and 2007 February data, and so their strengths must have decreased in the interim by a factor of 9 or more. We discuss three mechanisms that could give rise to this unusual outflow.

On October 28, 2003 an earthward-directed coronal mass ejection (CME) was observed from SOHO/LASCO imagery in conjunction with an X17 solar flare. The CME, traveling at nearly 2000 km/s, impacted the Earth on October 29, 2003 causing ground-based particle detectors to register a counting rate drop known as a Forbush decrease. In addition to affecting the rate of cosmic rays, the CME was also responsible for causing anisotropies in the direction of incidence. Data from Project GRAND, an array of proportional wire chambers, are presented during the time of this Forbush decrease. A simple model for CME propagation is proposed and we present an argument based on gyroradius that shows that a magnetic field of the radius calculated for the ejecta is sufficient to deflect energetic charged particles of an energy detectable by GRAND.

Since the pioneering observations of Spite & Spite in 1982, the constant lithium abundance of metal-poor ([Fe/H]<-1.3) halo stars near the turnoff has been attributed a cosmological origin. Closer analysis revealed that the observed abundance lies at $\Delta$ Li~0.4 dex below the predictions of Big Bang Nucleosynthesis. The measurements of deuterium abundances on the lines of sight toward quasars and the recent data from the Wilkinson Microwave Anisotropy Probe independently confirmed this gap. We suggest here that part of the discrepancy is explained by the first generation of stars that efficiently depleted lithium. Assuming that the models for lithium evolution in the halo turnoff stars and $\Delta$ Li estimates are correct, we infer that between 1/3 and 1/2 of the baryonic matter of the early halo (~10^9 Mo) was processed through Population III stars. This new paradigm proposes a very economical solution to the lingering difficulty of understanding the properties of the Spite Plateau and its lack of star-to-star scatter down to [Fe/H]=-2.5. It is moreover in agreement both with the absence of lithium detection in the most metal-poor star presently known (HE 1327-2326), and also with new trends of the Plateau suggesting its low metallicity edge may be reached around [Fe/H]=-2.5. We discuss the turbulent mixing associated with enhanced supernovae explosions in the early interstellar medium in this picture. Finally we show that other chemical properties of the extremely metal-poor stars are in agreement with a significant Population III processing in the halo, provided these models include mass-loss and rotationally-induced mixing.

We explore the role of anisotropic thermal conduction on the evolution of supernova remnants through interstellar media with a range of densities via numerical simulations. We find that a remnant expanding in a dense environment can produce centre-bright hard x-ray emission within 20 kyr, and centre-bright soft x-ray emission within 60 kyr of the supernova event. In a more tenuous environment, the appearance of a centre-bright structure in hard x-rays is delayed until about 60 kyr. The soft x-ray emission from such a remnant may not become centre bright during its observable lifetime. This can explain the observations that show that mixed-morphology supernova remnants preferentially occur close to denser, molecular environments. Remnants expanding into denser environments tend to be smaller, making it easier for thermal conduction to make larger changes in the temperatures of their hot gas bubbles. We show that the lower temperatures make it very favorable to use high-stage ions as diagnostics of the hot gas bubbles in SNRs. In particular, the distribution of O VIII transitions from shell-bright at early epochs to centre-bright at later epochs in the evolution of an SNR expanding in a dense ISM when the physics of thermal conduction is included.

We explore the importance of anisotropic thermal conduction in the evolution of supernova remnants via numerical simulations. The mean temperature of the bubble of hot gas is decreased by a factor of ~3 compared to simulations without thermal conduction, together with an increase in the mean density of hot gas by a similar factor. Thus, thermal conduction greatly reduces the volume of hot gas produced over the life of the remnant. This underscores the importance of thermal conduction in estimating the hot gas filling fraction and emissivities in high-stage ions in Galactic and proto-galactic ISMs.

We present high-resolution infrared spectra of four YSOs (T Tau N, T Tau S, RNO 91, and HL Tau). The spectra exhibit narrow absorption lines of 12CO, 13CO, and C18O as well as broad emission lines of gas phase12CO. The narrow absorption lines of CO are shown to originate from the colder circumstellar gas. We find that the line of sight gas column densities resulting from the CO absorption lines are much higher than expected for the measured extinction for each source and suggest the gas to dust ratio is measuring the dust settling and/or grain coagulation in these extended disks. We provide a model of turbulence, dust settling and grain growth to explain the results. The techniques presented here allow us to provide some observationally-motivated bounds on accretion disk alpha in protostellar systems.

For a slowly rotating non-magnetized white dwarf the accretion disk extends all the way to the star. Here the matter impacts and spreads towards the poles as new matter continuously piles up behind it. We have solved the 3d compressible Navier-Stokes equations on an axisymmetric grid to determine the structure of this boundary layer for different viscosities corresponding to different accretion rates. The high viscosity cases show a spreading BL which sets off a gravity wave in the surface matter. The accretion flow moves supersonically over the cusp making it susceptible to the rapid development of gravity wave and/or Kelvin-Helmholtz instabilities. This BL is optically thick and extends more than 30 degrees to either side of the disk plane after 3/4 of a Keplerian rotation period (t=19s). The low viscosity cases also show a spreading BL, but here the accretion flow does not set off gravity waves and it is optically thin.

We describe the results of numerical simulations of the dynamics of the boundary layer (BL) between the accretion disk and the surface of a non-magnetic white dwarf (WD) for different viscosities which corresponds to different stages for dwarf novae burst cycles. The simulations cover the inner part of the accretion disk, the BL, and the upper atmosphere of the star. The high viscosity case, which corresponds to a dwarf nova in outburst, shows a optically thick BL which after one Keplerian rotation period (t_K=19s) extends more than 30 degrees to either side of the disk plane. The BL is optically thick and thus occludes part of the star. The low viscosity case, which corresponds to a dwarf nova in quiescense, also shows a BL, but it is optically thin.

In this paper we further advance the study of magnetic field amplification in the interstellar medium that was started in Balsara et al (2004, Paper I). We show that the flux growth rate is comparable to the rate of magnetic energy growth found in Paper I. We also demonstrate the role of intermittency in field amplification. The density shows a double-peaked PDF, consistent with the cooling curve that was used. The PDF of the magnetic field shows a high-end tail, providing a tell-tale signature of the operation of the small scale dynamo. The magnetic field strength correlates positively with the density. As a result, the field amplification takes place more vigorously in the lower temperature, denser gas. The Lagrangian chaos in the simulated turbulent flows is studied in substantial detail. It is shown that the stretching rate of material lines as well as the Lyapunov exponents can be used to gain important insights into the growth of magnetic field. The cancellation exponent for the small scale supernova-driven dynamo is derived and it is shown that constructive folding of field lines in the dynamo is very inefficient. We also show that our Lagrangian approach can yield actual measures of the turbulent diffusivity in the simulated ISM. The turbulent diffusivity provides insights into the mixing of elements from supernova ejecta on macroscopic scales. The high rates of line stretching in interstellar turbulence suggests that the eventual diffusion of elements at the molecular level is very efficient. Many of the diagnostics of turbulence that are presented here can be used to make direct connections between simulations and observations.

Observations have suggested substantial departures from pressure equilibrium in the interstellar medium (ISM) in the plane of the Galaxy, even on scales under 50 pc. Nevertheless, multi-phase models of the ISM assume at least locally isobaric gas. The pressure then determines the density reached by gas cooling to stable thermal equilibrium. We use numerical models of the magnetized ISM to examine the consequences of supernova driving for interstellar pressures. In this paper we examine a (200 pc)^3 periodic domain threaded by magnetic fields. Individual parcels of gas at different pressures reach widely varying points on the thermal equilibrium curve: no unique set of phases is found, but rather a dynamically-determined continuum of densities and temperatures. A substantial fraction of the gas remains entirely out of thermal equilibrium. Our results appear consistent with observations of interstellar pressures. They also suggest that the high pressures observed in molecular clouds may be due to ram pressures in addition to gravitational forces. Much of the gas in our model lies far from equipartition between thermal and magnetic pressures, with ratios ranging from 0.1 to $10^4$ and ratios of uniform to fluctuating magnetic field of 0.5--1. Our models show broad pressure probability distribution functions with log-normal functional forms produced by both shocks and rarefaction waves, rather than power-law distributions produced by isolated supernova remnants. The width of the distribution can be described quantitatively by a formula derived from the work of Padoan, Nordlund, & Jones (1997).

Several lines of evidence suggest that protogalactic and galactic environments manage to grow magnetic field very rapidly. This makes the theoretical problem of achieving rapid growth of magnetic fields in such environments difficult. Following Balsara, Benjamin and Cox (2001), we show that a multiphase ISM that is dominated by realistic energy inputs via supernova explosions will naturally become a strongly turbulent medium with large positive and negative values of the fluid helicity. The resultant interstellar turbulence seems to have some of the essential characteristics of a fast dynamo. Magnetic field growth in such a turbulent dynamo is studied in this paper. It is shown that the field undergoes a linear phase of evolution where it has a time of growth that is comparable to the eddy turn-over time. Energetics, power spectra, statistics and structures of the turbulent dynamo are studied here. It is shown that despite the fact that the ISM is driven by compressive motions, the kinetic energy in this high Mach number flow is mainly concentrated in solenoidal motions rather than compressive motions.

In this paper we study the problem of divergence-free numerical MHD and show that the work done so far still has four key unresolved issues. We resolve those issues in this paper. The problem of reconstructing MHD flow variables with spatially second order accuracy is also studied. The other goal of this paper is to show that the same well-designed second order accurate schemes can be formulated for more complex geometries such as cylindrical and spherical geometry. Being able to do divergence-free reconstruction in those geometries also resolves the problem of doing AMR in those geometries. The resulting MHD scheme has been implemented in Balsara's RIEMANN framework for parallel, self-adaptive computational astrophysics. The present work also shows that divergence-free reconstruction and the divergence-free time-update can be done for numerical MHD on unstructured meshes. All the schemes designed here are shown to be second order accurate. Several stringent test problems are presented to show that the methods work, including problems involving high velocity flows in low plasma-b magnetospheric environments.

In recent years, several different strategies have emerged for evolving the magnetic field in numerical MHD. Some of these methods can be classified as divergence-cleaning schemes, where one evolves the magnetic field components just like any other variable in a higher order Godunov scheme. The fact that the magnetic field is divergence-free is imposed post-facto via a divergence-cleaning step. Other schemes for evolving the magnetic field rely on a staggered mesh formulation which is inherently divergence-free. The claim has been made that the two approaches are equivalent. In this paper we cross-compare three divergence-cleaning schemes based on scalar and vector divergence-cleaning and a popular divergence-free scheme. All schemes are applied to the same stringent test problem. Several deficiencies in all the divergence-cleaning schemes become clearly apparent with the scalar divergence-cleaning schemes performing worse than the vector divergence-cleaning scheme. The vector divergence-cleaning scheme also shows some deficiencies relative to the staggered mesh divergence-free scheme. The differences can be explained by realizing that all the divergence-cleaning schemes are based on a Poisson solver which introduces a non-locality into the scheme, though other subtler points of difference are also catalogued. By using several diagnostics that are routinely used in the study of turbulence, it is shown that the differences in the schemes produce measurable differences in physical quantities that are of interest in such studies.

We have studied forced turbulence of compressible magnetohydrodynamic (MHD) flows through two-dimensional simulations with different numerical resolutions. First, hydrodynamic turbulence with Mach number $<M_s >_{\rm init} \equiv < v >_{\rm rms}/ c_s = 1$ and density compression ${< \delta\rho / \rho >}_{\rm rms} \simeq 0.45$ was generated by enforcing a random force. Then, initial, uniform magnetic fields of various strengths were added with Alfv\'enic Mach number $<M_A >_{\rm init} \equiv < v >_{\rm rms} / c_{A, {\rm init}} \gg 1$. An isothermal equation of state was employed, and no explicit dissipation was included. After the MHD turbulence is saturated, the resulting flows are categorized as very weak field (VWF), weak field (WF), and strong field (SF) classes, which have $<M_A > \equiv < v >_{\rm rms} / < c_A >_{\rm rms} \gg 1$, $<M_A > > 1$, and $<M_A > \sim 1$, respectively. Not only in the SF regime but also in the WF regime, turbulent transport is suppressed by the magnetic field. In the SF cases, the energy power spectra in the inertial range, although no longer power-law, exhibit a range with slopes close to $\sim1.5$, hinting the Iroshnikov-Kraichnan spectrum. Our simulations were able to produce the SF class behaviors only with high resolution of at least $1024^2$ grid cells. The specific requirements for the simulation of the SF class should depend on the code (and the numerical scheme) as well as the initial setup, but our results do indicate that very high resolution would be required for converged results in simulation studies of MHD turbulence.

The advent of robust, reliable and accurate higher order Godunov schemes for many of the systems of equations of interest in computational astrophysics has made it important to understand how to solve them in multi-scale fashion. This is so because the physics associated with astrophysical phenomena evolves in multi-scale fashion and we wish to arrive at a multi-scale simulational capability to represent the physics. Because astrophysical systems have magnetic fields, multi-scale magnetohydrodynamics (MHD) is of especial interest. In this paper we first discuss general issues in adaptive mesh refinement (AMR). We then focus on the important issues in carrying out divergence-free AMR-MHD and catalogue the progress we have made in that area. We show that AMR methods lend themselves to easy parallelization. We then discuss applications of the RIEMANN framework for AMR-MHD to problems in computational astophysics.

In this paper we present a full-fledged scheme for the second order accurate, divergence-free evolution of vector fields on an adaptive mesh refinement (AMR) hierarchy. We focus here on adaptive mesh MHD. The scheme is based on making a significant advance in the divergence-free reconstruction of vector fields. In that sense, it complements the earlier work of Balsara and Spicer (1999) where we discussed the divergence-free time-update of vector fields which satisfy Stoke's law type evolution equations. Our advance in divergence-free reconstruction of vector fields is such that it reduces to the total variation diminishing (TVD) property for one-dimensional evolution and yet goes beyond it in multiple dimensions. Divergence-free restriction is also discussed. An electric field correction strategy is presented for use on AMR meshes. The electric field correction strategy helps preserve the divergence-free evolution of the magnetic field even when the time steps are sub-cycled on refined meshes. The above-mentioned innovations have been implemented in Balsara's RIEMANN framework for parallel, self-adaptive computational astrophysics which supports both non-relativistic and relativistic MHD. Several rigorous, three dimensional AMR-MHD test problems with strong discontinuities have been run with the RIEMANN framework showing that the strategy works very well.

(Abridged) We present the results of three dimensional calculations for the MHD evolution of an adiabatic supernova remnant in both a uniform and turbulent interstellar medium using the RIEMANN framework of Balsara. In the uniform case, which contains an initially uniform magnetic field, the density structure of the shell remains largely spherical, while the magnetic pressure and synchrotron emissivity are enhanced along the plane perpendicular to the field direction. This produces a bilateral or barrel-type morphology in synchrotron emission for certain viewing angles. We then consider a case with a turbulent external medium as in Balsara & Pouquet, characterized by $v_{A}(rms)/c_{s}=2$. Several important changes are found. First, despite the presence of a uniform field, the overall synchrotron emissivity becomes approximately spherically symmetric, on the whole, but is extremely patchy and time-variable, with flickering on the order of a few computational time steps. We suggest that the time and spatial variability of emission in early phase SNR evolution provides information on the turbulent medium surrounding the remnant. The shock-turbulence interaction is also shown to be a strong source of helicity-generation and, therefore, has important consequences for magnetic field generation. We compare our calculations to the Sedov-phase evolution, and discuss how the emission characteristics of SNR may provide a diagnostic on the nature of turbulence in the pre-supernova environment.

Observations have suggested substantial departures from pressure equilibrium in the interstellar medium (ISM) in the plane of the Galaxy, even on scales under 50 pc. Nevertheless, multi-phase models of the ISM assume at least locally isobaric gas. The pressure then determines the density reached by gas cooling to stable thermal equilibrium. We use two different sets of numerical models of the ISM to examine the consequences of supernova driving for interstellar pressures. The first set of models is hydrodynamical, and uses adaptive mesh refinement to allow computation of a 1 x 1 x 20 kpc section of a stratified galactic disk. The second set of models is magnetohydrodynamical, using an independent code framework, and examines a 200 pc cubed periodic domain threaded by magnetic fields. Both of these models show broad pressure distributions with roughly log-normal functional forms produced by both shocks and rarefaction waves, rather than the power-law distributions predicted by previous work, with rather sharp thermal pressure gradients. The width of the distribution of the logs of pressure in gas with log T < 3.9 is proportional to the rms Mach number in that gas, while the distribution in hotter gas is broader, but not so broad as would be predicted by the Mach numbers in that gas. Individual parcels of gas reach widely varying points on the thermal equilibrium curve: no unique set of phases is found, but rather a dynamically-determined continuum of densities and temperatures. Furthermore, a substantial fraction of the gas remains entirely out of thermal equilibrium. Our results appear consistent with observations of interstellar pressures, and suggest that the pressures observed in molecular clouds may be due to ram pressure rather than gravitational confinement.

We describe the results of a sequence of simulations of gravitational collapse in a turbulent magnetized region. The parameters are chosen to be representative of molecular cloud material. We find that several protostellar cores and filamentary structures of higher than average density form. The filaments inter-connect the high density cores. Furthermore, the magnetic field strengths are found to correlate positively with the density, in agreement with recent observations. We make synthetic channel maps of the simulations and show that material accreting onto the cores is channelled along the magnetized filamentary structures. This is compared with recent observations of S106, and shown to be consistent with these data. We postulate that this mechanism of accretion along filaments may provide a means for molecular cloud cores to grow to the point where they become gravitationally unstable and collapse to form stars.