Extending our recent studies of two-dimensional stellar convection to 3D, we compare three-dimensional hydrodynamic simulations to identically set-up two-dimensional simulations, for a realistic pre-main sequence star. We compare statistical quantities related to convective flows including: average velocity, vorticity, local enstrophy, and penetration depth beneath a convection zone. These statistics are produced during stationary, steady-state compressible convection in the star's convection zone. Our simulations with the MUSIC code confirm the common result that two-dimensional simulations of stellar convection have a higher magnitude of velocity on average than three-dimensional simulations. Boundary conditions and the extent of the spherical shell can affect the magnitude and variability of convective velocities. The difference between 2D and 3D velocities is dependent on these background points; in our simulations this can have an effect as large as the difference resulting from the dimensionality of the simulation. Nevertheless, radial velocities near the convective boundary are comparable in our 2D and 3D simulations. The average local enstrophy of the flow is lower for two-dimensional simulations than for three-dimensional simulations, indicating a different shape and structuring of 3D stellar convection. We perform a statistical analysis of the depth of convective penetration below the convection zone, using the model proposed in our recent study (Pratt et al. 2017). Here we analyze the convective penetration in three dimensional simulations, and compare the results to identically set-up 2D simulations. In 3D the penetration depth is as large as the penetration depth calculated from 2D simulations.

This review aims at providing an up-to-date status and a general introduction to the subject of the numerical study of energetic particle acceleration and transport in turbulent astrophysical flows. The subject is also complemented by a short overview of recent progresses obtained in the domain of laser plasma experiments. We review the main physical processes at the heart of the production of a non-thermal distribution in both Newtonian and relativistic astrophysical flows, namely the first and second order Fermi acceleration processes. We also discuss shock drift and surfing acceleration, two processes important in the context of particle injection in shock acceleration. We analyze with some details the particle-in-cell (PIC) approach used to describe particle kinetics. We review the main results obtained with PIC simulations in the recent years concerning particle acceleration at shocks and in reconnection events. The review discusses the solution of Fokker-Planck problems with application to the study of particle acceleration at shocks but also in hot coronal plasmas surrounding compact objects. We continue by considering large scale physics. We describe recent developments in magnetohydrodynamic (MHD) simulations. We give a special emphasize on the way energetic particle dynamics can be coupled to MHD solutions either using a multi-fluid calculation or directly coupling kinetic and fluid calculations. This aspect is mandatory to investigate the acceleration of particles in the deep relativistic regimes to explain the highest Cosmic Ray energies.

We present the results from a particle-in-cell (PIC) simulation that models the interaction between a spatially localized electron-positron cloud and an electron-ion plasma. The latter is permeated by a magnetic field that is initially spatially uniform and aligned with the mean velocity vector of the pair cloud. The pair cloud expels the magnetic field and piles it up into an electromagnetic piston. Its electromagnetic field is strong enough to separate the pair cloud from the ambient plasma in the direction that is perpendicular to the cloud propagation direction. The piston propagates away from the spine of the injected pair cloud and it accelerates the protons to a high nonrelativistic speed. The accelerated protons form an outer cocoon that will eventually become separated from the unperturbed ambient plasma by a fast magnetosonic shock. No electromagnetic piston forms at the front of the cloud and a shock is mediated here by the filamentation instability. The final plasma distribution resembles that of a hydrodynamic jet. Collisionless plasma jets may form in the coronal plasma of accreting black holes and the interaction between the strong magnetic field of the piston and the hot pair cloud may contribute to radio emissions by such objects.

Characterizing stellar convection in multiple dimensions is a topic at the forefront of stellar astrophysics. Numerical simulations are an essential tool for this task. We present an extension of the existing numerical tool-kit A-MaZe that enables such simulations of stratified flows in a gravitational field. The finite-volume based, cell-centered, and time-explicit hydrodynamics solver of A-MaZe was extended such that the scheme is now well-balanced in both momentum and energy. The algorithm maintains an initially static balance between gravity and pressure to machine precision. Quasi-stationary convection in slab-geometry preserves gas energy (internal plus kinetic) on average despite strong local up- and down-drafts. By contrast, a more standard numerical scheme is demonstrated to result in substantial gains of energy within a short time on purely numerical grounds. The test is further used to point out the role of dimensionality, viscosity, and Rayleigh number for compressible convection. Applications to a young sun in 2D and 3D, covering a part of the inner radiative zone as well as the outer convective zone, demonstrate that the scheme meets its initial design goal. Comparison with results obtained for a physically identical setup with a time-implicit code show qualitative agreement.

We study the effect a guiding magnetic field has on the formation and structure of a pair jet that propagates through a collisionless electron-proton plasma at rest. We model with a PIC simulation a pair cloud with the temperature 400 keV and mean speed 0.9c. The cloud propagates through a spatially uniform, magnetized and cool ambient electron-proton plasma that is at rest. Its mean velocity vector is aligned with the background magnetic field. A jet forms in time. Its outer cocoon consists of jet-accelerated ambient plasma and is separated from the inner cocoon by an electromagnetic piston with a thickness that is comparable to the thermal gyroradius of jet particles. A supercritical fast magnetosonic shock will form between the pristine ambient plasma and the jet-accelerated one on a time scale that exceeds our simulation time by an order of magnitude. The inner cocoon is pair plasma that lost its directed flow energy while it swept out the background magnetic field. A beam of electrons and positrons moves along the jet spine at its initial speed. Its electrons are slowed down and some positrons are accelerated as they cross the jet's head. The latter escape upstream along the magnetic field, which yields an excess of MeV positrons ahead of the jet. Some of the protons, which were located behind the electromagnetic piston at the time it formed, are accelerated to MeV energies

By modelling the expansion of a cloud of electrons and positrons with the temperature 400 keV that propagates at the mean speed 0.9c ($c:$ speed of light) through an initially unmagnetized electron-proton plasma with a particle-in-cell (PIC) simulation, we find a mechanism that collimates the pair cloud into a jet. A filamentation instability develops between the protons at rest and the moving positrons. Its magnetic field collimates the positrons and drives an electrostatic shock into the electron-proton plasma. The magnetic field acts as a discontinuity that separates the protons of the shocked ambient plasma, known as the outer cocoon, from the jet's interior region. The outer cocoon expands at the speed 0.15c along the jet axis and at 0.03c perpendicularly to it. The filamentation instability converts the jet's directed flow energy into magnetic energy in the inner cocoon. The magnetic discontinuity cannot separate the ambient electrons from the jet electrons. Both species rapidly mix and become indistinguishable. The spatial distribution of the positive charge carriers is in agreement with the distributions of the ambient material and the jet material predicted by a hydrodynamic model apart from a dilute positronic outflow that is accelerated by the electromagnetic field at the jet's head.

We study lithium depletion in low-mass and solar-like stars as a function of time, using a new diffusion coefficient describing extra-mixing taking place at the bottom of a convective envelope. This new form is motivated by multi-dimensional fully compressible, time implicit hydrodynamic simulations performed with the MUSIC code. Intermittent convective mixing at the convective boundary in a star can be modeled using extreme value theory, a statistical analysis frequently used for finance, meteorology, and environmental science. In this letter, we implement this statistical diffusion coefficient in a one-dimensional stellar evolution code, using parameters calibrated from multi-dimensional hydrodynamic simulations of a young low-mass star. We propose a new scenario that can explain observations of the surface abundance of lithium in the Sun and in clusters covering a wide range of ages, from $\sim$ 50 Myr to $\sim$ 4 Gyr. Because it relies on our physical model of convective penetration, this scenario has a limited number of assumptions. It can explain the observed trend between rotation and depletion, based on a single additional assumption, namely that rotation affects the mixing efficiency at the convective boundary. We suggest the existence of a threshold in stellar rotation rate above which rotation strongly prevents the vertical penetration of plumes and below which rotation has small effects. In addition to providing a possible explanation for the long standing problem of lithium depletion in pre-main sequence and main sequence stars, the strength of our scenario is that its basic assumptions can be tested by future hydrodynamic simulations.

We examine a penetration layer formed between a central radiative zone and a large convection zone in the deep interior of a young low-mass star. Using the Multidimensional Stellar Implicit Code (MUSIC) to simulate two-dimensional compressible stellar convection in a spherical geometry over long times, we produce statistics that characterize the extent and impact of convective penetration in this layer. We apply extreme value theory to the maximal extent of convective penetration at any time. We compare statistical results from simulations which treat non-local convection, throughout a large portion of the stellar radius, with simulations designed to treat local convection in a small region surrounding the penetration layer. For each of these situations, we compare simulations of different resolution, which have different velocity magnitudes. We also compare statistical results between simulations that radiate energy at a constant rate to those that allow energy to radiate from the stellar surface according to the local surface temperature. Based on the frequency and depth of penetrating convective structures, we observe two distinct layers that form between the convection zone and the stable radiative zone. We show that the probability density function of the maximal depth of convective penetration at any time corresponds closely in space with the radial position where internal waves are excited. We find that the maximal penetration depth can be modeled by a Weibull distribution with a small shape parameter. Using these results, and building on established scalings for diffusion enhanced by large-scale convective motions, we propose a new form for the diffusion coefficient that may be used for one-dimensional stellar evolution calculations in the large P\'eclet number regime. These results should contribute to the 321D link.

We present the results of a numerical benchmark study for the MUlti-dimensional Stellar Implicit Code (MUSIC) based on widely applicable two- and three-dimensional compressible hydrodynamics problems relevant to stellar interiors. MUSIC is an implicit large eddy simulation code that uses implicit time integration, implemented as a Jacobian-free Newton Krylov method. A physics based preconditioning technique which can be adjusted to target varying physics is used to improve the performance of the solver. The problems used for this benchmark study include the Rayleigh-Taylor and Kelvin-Helmholtz instabilities, and the decay of the Taylor-Green vortex. Additionally we show a test of hydrostatic equilibrium, in a stellar environment which is dominated by radiative effects. In this setting the flexibility of the preconditioning technique is demonstrated. This work aims to bridge the gap between the hydrodynamic test problems typically used during development of numerical methods and the complex flows of stellar interiors. A series of multi-dimensional tests are performed and analysed. Each of these test cases is analysed with a simple, scalar diagnostic, with the aim of enabling direct code comparisons. As the tests performed do not have analytic solutions we verify MUSIC by comparing to established codes including ATHENA and the PENCIL code. MUSIC is able to both reproduce behaviour from established and widely-used codes as well as results expected from theoretical predictions. This benchmarking study concludes a series of papers describing the development of the MUSIC code and provides confidence in the future applications.

Context: We study the impact of two-dimensional spherical shells on compressible convection. Realistic profiles for density and temperature from a one-dimensional stellar evolution code are used to produce a model of a large stellar convection zone representative of a young low-mass star. Methods: We perform hydrodynamic implicit large-eddy simulations of compressible convection using the MUltidimensional Stellar Implicit Code (MUSIC). Because MUSIC has been designed to use realistic stellar models produced from one-dimensional stellar evolution calculations, MUSIC simulations are capable of seamlessly modeling a whole star. Simulations in two-dimensional spherical shells that have different radial extents are performed over hundreds of convective turnover times, permitting the collection of well-converged statistics. Results: We evaluate basic statistics of the convective turnover time, the convective velocity, and the overshooting layer. These quantities are selected for their relevance to one-dimensional stellar evolution calculations, so that our results are focused toward the 321D link. The inclusion in the spherical shell of the boundary between the radiative and convection zones decreases the amplitude of convective velocities in the convection zone. The inclusion of near-surface layers in the spherical shell can increase the amplitude of convective velocities, although the radial structure of the velocity profile established by deep convection is unchanged. The impact from including the near-surface layers depends on the speed and structure of small-scale convection in the near-surface layers. Larger convective velocities in the convection zone result in a commensurate increase in the overshooting layer width and decrease in the convective turnover time. These results provide support for non-local aspects of convection.

This work is the first attempt to describe the multi-dimensional structure of accreting young stars based on fully compressible time implicit multi-dimensional hydrodynamics simulations. One major motivation is to analyse the validity of accretion treatment used in previous 1D stellar evolution studies. We analyse the effect of accretion on the structure of a realistic stellar model of the young Sun. Our work is inspired by the numerical work of Kley \& Lin (1996, ApJ, 461, 933) devoted to the structure of the boundary layer in accretion disks. We analyse the redistribution of accreted material with a range of values of specific entropy relative to the bulk specific entropy of the material in the accreting object's convective envelope. A primary goal is to understand whether and how accreted energy deposited onto a stellar surface is redistributed in the interior. This study focusses on the high accretion rates characteristic of FU Ori systems. We find that the highest entropy cases produce a distinctive behaviour in the mass redistribution, rms velocities, and enthalpy flux in the convective envelope. This change in behaviour is characterised by the formation of a hot layer on the surface of the accreting object, which tends to suppress convection in the envelope. We analyse the long-term effect of such a hot buffer zone on the structure and evolution of the accreting object with 1D stellar evolution calculations. We study the relevance of the assumption of redistribution of accreted energy into the stellar interior used in the literature. One conclusion is that, for a given amount of accreted energy transferred to the accreting object, a treatment assuming accretion energy redistribution throughout the stellar interior could significantly overestimate the effects on the stellar structure, in particular, on the resulting expansion.

A model, which is widely used for inertial rang statistics of supersonic turbulence in the context of molecular clouds and star formation, expresses (measurable) relative scaling exponents Z_p of two-point velocity statistics as a function of two parameters, beta and Delta. The model relates them to the dimension D of the most dissipative structures, D=3-Delta/(1-beta). While this description has proved most successful for incompressible turbulence (beta=Delta=2/3, and D=1), its applicability in the highly compressible regime remains debated. For this regime, theoretical arguments suggest D=2 and Delta=2/3, or Delta=1. Best estimates based on 3D periodic box simulations of supersonic isothermal turbulence yield Delta=0.71 and D=1.9, with uncertainty ranges of Delta in [0.67, 0.78] and D in [2.04,1.60]. With these 5-10\% uncertainty ranges just marginally including the theoretical values of Delta=2/3 and D=2, doubts remain whether the model indeed applies and, if it applies, for what values of beta and Delta. We use a Monte Carlo approach to mimic actual simulation data and examine what factors are most relevant for the fit quality. We estimate that 0.1% (0.05%) accurate Z_p, with p=1...5, should allow for 2% (1%) accurate estimates of beta and Delta in the highly compressible regime, but not in the mildly compressible regime. We argue that simulation-based Z_p with such accuracy are within reach of today's computer resources. If this kind of data does not allow for the expected high quality fit of beta and Delta, then this may indicate the inapplicability of the model for the simulation data. In fact, other models than the one we examine here have been suggested.

This work is a continuation of our efforts to develop an efficient implicit solver for multidimensional hydrodynamics for the purpose of studying important physical processes in stellar interiors, such as turbulent convection and overshooting. We present an implicit solver that results from the combination of a Jacobian-Free Newton-Krylov method and a preconditioning technique tailored to the inviscid, compressible equations of stellar hydrodynamics. We assess the accuracy and performance of the solver for both 2D and 3D problems for Mach numbers down to $10^{-6}$. Although our applications concern flows in stellar interiors, the method can be applied to general advection and/or diffusion-dominated flows. The method presented in this paper opens up new avenues in 3D modeling of realistic stellar interiors allowing the study of important problems in stellar structure and evolution.

The thin-shell instability has been named as one process, which can generate entangled structures in astrophysical plasma on collisional (fluid) scales. It is driven by a spatially varying imbalance between the ram pressure of the inflowing upstream plasma and the downstream's thermal pressure at a non-planar shock. Here we show by means of a particle-in-cell (PIC) simulation that an analogue process can destabilize a thin shell formed by two interpenetrating, unmagnetized and collisionless plasma clouds. The amplitude of the shell's spatial modulation grows and saturates after about ten inverse proton plasma frequencies, when the shell consists of connected piecewise linear patches.

Collisionless magnetic reconnection is a prime candidate to account for flare-like or steady emission, outflow launching, or plasma heating, in a variety of high-energy astrophysical objects, including ones with relativistic ion-electron plasmas. But the fate of the initial magnetic energy in a reconnection event remains poorly known: what is the amount given to kinetic energy, the ion/electron repartition, and the hardness of the particle distributions? We explore these questions with 2D particle-in-cell simulations of ion-electron plasmas. We find that 45 to 75% of the total initial magnetic energy ends up in kinetic energy, this fraction increasing with the inflow magnetization. Depending on the guide field strength, ions get from 30 to 60% of the total kinetic energy. Particles can be separated into two populations that only weakly mix: (i) particles initially in the current sheet, heated by its initial tearing and subsequent contraction of the islands; and (ii) particles from the background plasma that primarily gain energy via the reconnection electric field when passing near the X-point. Particles (ii) tend to form a power-law with an index $p=-d\log n(\gamma)/d\log\gamma$, that depends mostly on the inflow Alfv\'en speed $V_A$ and magnetization $\sigma_s$ of species $s$, with for electrons $p=5$ to $1.2$ for increasing $\sigma_e$. The highest particle Lorentz factor, for ions or electrons, increases roughly linearly with time for all the relativistic simulations. This is faster, and the spectra can be harder, than for collisionless shock acceleration. We discuss applications to microquasar and AGN coronae, to extragalactic jets, and to radio lobes. We point out situations where effects such as Compton drag or pair creation are important.

We present first results of a long-term project which aims at multi-scale, multi-physics simulations of wind accretion in microquasars and high-mass X-ray binaries. The 3D hydrodynamical simulations cover all scales, from the circum-binary environment down to the immediate vicinity of the black hole. We first introduce the numerical method and parallelization strategy of the AMR A-MAZE code. We then discuss some preliminary results of how, and on what scales, an accretion disk is formed around the black hole. We finally present some characteristics of this disk, which is far from Keplerian. We emphasize that on all scales shocks play a decisive role for the accretion process and the process of structure formation -- for the formation of the large scale, nearly coherent structure of the disk, but also for the formation of turbulent fluctuations.

Magnetic reconnection is a leading mechanism for magnetic energy conversion and high-energy non-thermal particle production in a variety of high-energy astrophysical objects, including ones with relativistic ion-electron plasmas (e.g., microquasars or AGNs) - a regime where first principle studies are scarce. We present 2D particle-in-cell (PIC) simulations of low $\beta$ ion-electron plasmas under relativistic conditions, i.e., with inflow magnetic energy exceeding the plasma rest-mass energy. We identify outstanding properties: (i) For relativistic inflow magnetizations (here $10 < \sigma_e < 360$), the reconnection outflows are dominated by thermal agitation instead of bulk kinetic energy. (ii) At large inflow electron magnetization ($\sigma_e > 80$), the reconnection electric field is sustained more by bulk inertia than by thermal inertia. It challenges the thermal-inertia-paradigm and its implications. (iii) The inflows feature sharp transitions at the entrance of the diffusion zones. These are not shocks but results from particle ballistic motions, all bouncing at the same location, provided that the thermal velocity in the inflow is far smaller than the inflow E cross B bulk velocity. (iv) Island centers are magnetically isolated from the rest of the flow, and can present a density depletion at their center. (v) The reconnection rates are slightly larger than in non-relativistic studies. They are best normalized by the inflow relativistic Alfv\'en speed projected in the outflow direction, which then leads to rates in a close range (0.14-0.25) thus allowing for an easy estimation of the reconnection electric field.

Colliding supersonic bulk flows shape observable properties and internal physics of various astrophysical objects, like O-star winds, molecular clouds, galactic sheets, binaries, or gamma-ray bursts. Using numerical simulations, we show that the bulk flows leave a clear imprint on the collision zone, its mean properties and the turbulence it naturally develops. Our model setup consists of 3D head-on colliding isothermal hydrodynamical flows with Mach numbers between 2 and 43. Simulation results are in line with expectations from self-similarity: root mean square Mach numbers (Mrms) scale linearly with upstream Mach numbers, mean densities remain limited to a few times the upstream density. The density PDF is not log-normal. The turbulence is inhomogeneous: weaker in the zone center than close to the confining shocks. It is anisotropic: while Mrms is generally supersonic, Mrms transverse to the upstream flow is always subsonic. We argue that uniform, isothermal, head-on colliding flows generally disfavor isotropic, supersonic turbulence. The anisotropy carries over to other quantities like the density variance - Mach number relation. Structure functions differ depending on whether they are computed along a line-of-sight perpendicular or parallel to the upstream flow. We suggest that such line-of-sight effects should be kept in mind when interpreting turbulence characteristics derived from observations.

The PIC model relies on two building blocks. The first stems from the capability of computers to handle only up to $\sim10^{10}$ particles, while real plasmas contain from $10^4$ to $10^{20}$ particles per Debye sphere: a coarse-graining step must be used, whereby of the order of $p\sim10^{10}$ real particles are represented by a single computer superparticle. The second is field storage on a grid with its subsequent finite superparticle size. We introduce the notion of coarse-graining dependent quantities, i.e. physical quantities depending on the number $p$. They all derive from the plasma parameter $\Lambda$, which we show to be proportional to $1/p$. We explore three examples: the rapid collision- and fluctuation-induced thermalization of plasmas with different temperatures, that scale with the number of superparticles per grid cell and are a factor $p\sim10^{10}$ faster than in real plasmas; the high level of electrostatic fluctuations in a thermal plasma, with corrections due to the finite superparticle sizes; and the blurring of the linear spectrum of the filamentation instability, where the fastest growing modes do not dominate the total energy because of a high level of fluctuations. We stress that the enhanced collisions and correlations of PIC plasmas must be kept negligible toward kinetic physics.

Aims. Observations show nebulae around some massive stars but not around others. If observed, their chemical composition is far from homogeneous. Our goal is to put these observational features into the context of the evolution of massive stars and their circumstellar medium (CSM) and, more generally, to quantify the role of massive stars for the chemical and dynamical evolution of the ISM. Methods. Using the A-MAZE code, we perform 2d-axisymmetric hydrodynamical simulations of the evolution of the CSM, shaped by stellar winds, for a whole grid of massive stellar models from 15 to 120 Msun and following the stellar evolution from the zero-age main-sequence to the time of supernova explosion. In addition to the usual quantities, we also follow five chemical species: H, He, C, N, and O. Results. We show how various quantities evolve as a function of time: size of the bubble, position of the wind termination shock, chemical composition of the bubble, etc. The chemical composition of the bubble changes considerably compared to the initial composition, particularly during the red-supergiant (RSG) and Wolf-Rayet (WR) phases. In some extreme cases, the inner region of the bubble can be completely depleted in hydrogen and nitrogen, and is mainly composed of carbon, helium and oxygen. We argue why the bubble typically expands at a lower rate than predicted by self-similarity theory. In particular, the size of the bubble is very sensitive to the density of the ISM, decreasing by a factor of around 2.5 for each additional dex in ISM density. The bubble size also decreases with the metallicity of the central star, as low-metallicity stars have weaker winds. Our models qualitatively fit the observations of WR ejecta nebulae.

We present the parallel particle-in-cell (PIC) code Apar-T and, more importantly, address the fundamental question of the relations between the PIC model, the Vlasov-Maxwell theory, and real plasmas. First, we present four validation tests: spectra from simulations of thermal plasmas, linear growth rates of the relativistic tearing instability and of the filamentation instability, and non-linear filamentation merging phase. For the filamentation instability we show that the effective growth rates measured on the total energy can differ by more than 50% from the linear cold predictions and from the fastest modes of the simulation. Second, we detail a new method for initial loading of Maxwell-J\"uttner particle distributions with relativistic bulk velocity and relativistic temperature, and explain why the traditional method with individual particle boosting fails. Third, we scrutinize the question of what description of physical plasmas is obtained by PIC models. These models rely on two building blocks: coarse-graining, i.e., grouping of the order of p~10^10 real particles into a single computer superparticle, and field storage on a grid with its subsequent finite superparticle size. We introduce the notion of coarse-graining dependent quantities, i.e., quantities depending on p. They derive from the PIC plasma parameter Lambda^{PIC}, which we show to scale as 1/p. We explore two implications. One is that PIC collision- and fluctuation-induced thermalization times are expected to scale with the number of superparticles per grid cell, and thus to be a factor p~10^10 smaller than in real plasmas. The other is that the level of electric field fluctuations scales as 1/Lambda^{PIC} ~ p. We provide a corresponding exact expression. Fourth, we compare the Vlasov-Maxwell theory, which describes a phase-space fluid with infinite Lambda, to the PIC model and its relatively small Lambda.

Time-implicit schemes are attractive since they allow numerical time steps that are much larger than those permitted by the Courant-Friedrich-Lewy criterion characterizing time-explicit methods. This advantage comes, however, with a cost: the solution of a system of nonlinear equations is required at each time step. In this work, the nonlinear system results from the discretization of the hydrodynamical equations with the Crank-Nicholson scheme. We compare the cost of different methods, based on Newton-Raphson iterations, to solve this nonlinear system, and benchmark their performances against time-explicit schemes. Since our general scientific objective is to model stellar interiors, we use as test cases two realistic models for the convective envelope of a red giant and a young Sun. Focusing on 2D simulations, we show that the best performances are obtained with the quasi-Newton method proposed by Broyden. Another important concern is the accuracy of implicit calculations. Based on the study of an idealized problem, namely the advection of a single vortex by a uniform flow, we show that there are two aspects: i) the nonlinear solver has to be accurate enough to resolve the truncation error of the numerical discretization, and ii) the time step has be small enough to resolve the advection of eddies. We show that with these two conditions fulfilled, our implicit methods exhibit similar accuracy to time-explicit schemes, which have lower values for the time step and higher computational costs. Finally, we discuss in the conclusion the applicability of these methods to fully implicit 3D calculations.

The medium around massive stars is strongly shaped by the stellar winds. Those winds depend on various stellar parameters (effective temperature, luminosity, chemical composition, rotation, ...), which are varying as a function of the time. Using the wind properties obtained by classical stellar evolution code allows for the multi- D hydrodynamical simulation of the circum-stellar medium accounting for the time variations of the wind. We present here the preliminary results of the simulation of the medium around a fast rotating star. Comparing the results of such simulations with the observed properties of the circum-stellar medium will allow to better understand the interactions of the stellar winds with the interstellar medium, and could give hints on the past history of the mass loss around massive stars, providing constraints for stellar models.

Massive stars are crucial building blocks of galaxies and the universe, as production sites of heavy elements and as stirring agents and energy providers through stellar winds and supernovae. The field of magnetic massive stars has seen tremendous progress in recent years. Different perspectives -- ranging from direct field measurements over dynamo theory and stellar evolution to colliding winds and the stellar environment -- fruitfully combine into a most interesting and still evolving overall picture, which we attempt to review here. Zeeman signatures leave no doubt that at least some O- and early B-type stars have a surface magnetic field. Indirect evidence, especially non-thermal radio emission from colliding winds, suggests many more. The emerging picture for massive stars shows similarities with results from intermediate mass stars, for which much more data are available. Observations are often compatible with a dipole or low order multi-pole field of about 1 kG (O-stars) or 300 G to 30 kG (Ap / Bp stars). Weak and unordered fields have been detected in the O-star zeta-Ori A and in Vega, the first normal A-type star with a magnetic field. Theory offers essentially two explanations for the origin of the observed surface fields: fossil fields, particularly for strong and ordered fields, or different dynamo mechanisms, preferentially for less ordered fields. Numerical simulations yield the first concrete stable (fossil) field configuration, but give contradictory results as to whether dynamo action in the radiative envelope of massive main sequence stars is possible. Internal magnetic fields, which may not even show up at the stellar surface, affect stellar evolution as they lead to a more uniform rotation, with more slowly rotating cores and faster surface rotation. Surface metallicities may become enhanced, thus affecting the mass-loss rates.

This paper describes the first steps of development of a new multidimensional time implicit code devoted to the study of hydrodynamical processes in stellar interiors. The code solves the hydrodynamical equations in spherical geometry and is based on the finite volume method. Radiation transport is taken into account within the diffusion approximation. Realistic equation of state and opacities are implemented, allowing the study of a wide range of problems characteristic of stellar interiors. We describe in details the numerical method and various standard tests performed to validate the method. We present preliminary results devoted to the description of stellar convection. We first perform a local simulation of convection in the surface layers of a A-type star model. This simulation is used to test the ability of the code to address stellar conditions and to validate our results, since they can be compared to similar previous simulations based on explicit codes. We then present a global simulation of turbulent convective motions in a cold giant envelope, covering 80% in radius of the stellar structure. Although our implicit scheme is unconditionally stable, we show that in practice there is a limitation on the time step which prevent the flow to move over several cells during a time step. Nevertheless, in the cold giant model we reach a hydro CFL number of 100. We also show that we are able to address flows with a wide range of Mach numbers (10^-3 < Ms< 0.5), which is impossible with an anelastic approach. Our first developments are meant to demonstrate that the use of an implicit scheme applied to a stellar evolution context is perfectly thinkable and to provide useful guidelines to optimise the development of an implicit multi-D hydrodynamical code.

We report on the development of an implicit multi-D hydrodynamic code for stellar evolution. We present two test-cases relevant for the first scientific goal of the code: the simulation of convection in pulsating stars. First results on a realistic stellar model are also presented.

Shock bound interaction zones (SBIZs) are ubiquitous in astrophysics. We present numerical results for 2D and 3D, plane-parallel, infinitely extended SBIZs. Isothermal settings and parameterized cooling are considered. We highlight and compare characteristic of such zones. We emphasize the mutual coupling between the turbulence within the SBIZ and the confining shocks, point out potential differences to 3D periodic box studies of supersonic turbulence, and contemplate on possible effects on the X-ray emission of such zones.

We present 3D hydrodynamical simulations of the separated binary RS Ophiuchi (RS Oph), a recurrent nova and potential progenitor of a SN Ia. RS Oph is composed of a red giant (RG) and a white dwarf (WD) whose mass is close to the Chandrasekhar limit. In an isothermal scenario, the WD accrets about 10% of a 20 km/s RG wind by a non-Keplerian accretion disk with strong spiral shocks, and about 2% of a 60 km/s RG wind by what we term a 'turbulent accretion ball'. A significantly larger impact have the thermodynamics. In an adiabatic scenario only about 0.7% of the 20 km/s RG wind is accreted. The rate of change of the system separation due to mass and angular momentum loss out of the system is negative in all three cases studied, but is ten times smaller for a fast RG wind (60 km/s) than for a slow RG wind (20 km/s). The results demonstrate that existing nova models and observed recurrence times fit well together with 3D wind accretion and that RS Oph is one of the most promising systems to become an SN Ia.

When the surface angular velocity is above about 70% of the critical angular velocity, many interesting features appear which may be tested by interferometric observations, like significant deformation of stars, variation of the effective temperature with the latitude. Also polar winds become important and equatorial disks may appear. Near the critical limit, convection is also favored in the outer layers. In the present paper, we emphasize the need for a proper estimate of the critical velocity since this is the ratio of the actual velocity of the star to that critical velocity which determines the amplitude of the above effects. We recall the existence of two critical velocities. The first one, also called the classical critical velocity is the one to consider when the star has an Eddington factor inferior to 0.639, while the second one is the one to be considered when the Eddington factor is above 0.639. The features of the star at these two critical limits may be very different.

We discuss the properties of the progenitors of core collapse SNe, as they can be deduced from rotating stellar models of single stars. The type of the SN progenitor was determined from the surface abundances at the pre-SN stage. The type of the SN event was obtained from the masses of hydrogen and helium ejected at the time of the core-collapse SN event. We find that the minimum amount of helium ejected by a core-collapse SN (of whatever type) is around 0.3 MSun. There is no difference between the WC and WO stars in the ejected masses of helium, CNO elements, and heavy elements. Also no difference is expected between the chemical composition of a WC star resulting from a normal or a homogeneous evolution. The progenitors of type Ib SNe are WNL, WNE, or less massive WC stars. Those of type Ic are WC and WO stars. WO stars are produced in a limited mass range (around 60 MSun) and only at low metallicity (for <= 0.010) as already found. The WO stars are the progenitors of only a small fraction of type Ic. Present stellar models indicate that, at solar metallicity, there is about 1 type Ib SN for 1 type Ic, and this ratio rises to 3 type Ic for 1 type Ib SN at twice solar metallicity. At this metallicity, type Ic's are more frequent than type Ib's because most massive stars that go through a WNE stage evolve further into a WC/WO phase. Current models can account for the observed number ratios SN Ib/SN II and SN Ic/SN II and for their observed variation with the metallicity. In case no supernova occurs when a black hole is formed, single-star models can still account for more than half of the observed (SN Ib+SN Ic)/SN II ratio for Z >= ZSun. For the GRB rate, our models produce too large a number for such an event, even if we restrict the progenitor to the WO stars.

RS Ophiuchi is a recurrent nova with a period of about 22 years, consisting of a wind accreting binary system with a white dwarf (WD) very close to the Chandrasekhar limit and a red giant star (RG). The system is considered a prime candidate to evolve into an SNIa. We present a 3D hydrodynamic simulation of the quiescent accretion and the subsequent explosive phase. The computed circumstellar mass distribution in the quiescent phase is highly structured with a mass enhancement in the orbital plane of about a factor of 2 as compared to the poleward directions. The simulated nova remnant evolves aspherically, propagating faster toward the poles. The shock velocities derived from the simulations are in agreement with those derived from observations. For v_RG = 20 km/s and for nearly isothermal flows, we derive a mass transfer rate to the WD of 10% of the mass loss of the RG. For an RG mass loss of 10^{-7} solar masses per year, we found the orbit of the system to decay by 3% per million years. With the derived mass transfer rate, multi-cycle nova models provide a qualitatively correct recurrence time, amplitude, and fastness of the nova. Our simulations provide, along with the observations and nova models, the third ingredient for a deeper understanding of the recurrent novae of the RS Oph type. In combination with recent multi-cycle nova models, our results suggests that the WD in RS Oph will increase in mass. Several speculative outcomes then seem plausible. The WD may reach the Chandrasekhar limit and explode as an SN Ia. Alternatively, the mass loss of the RG could result in a smaller Roche volume, a common envelope phase, and a narrow WD+WD system. Angular momentum loss due to graviational wave emission could trigger the merger of the two WDs and - perhaps - an SN Ia via the double degenerate scenario.

We have carried out a hydrodynamical code comparison study of interacting multiphase fluids. The two commonly used techniques of grid and smoothed particle hydrodynamics (SPH) show striking differences in their ability to model processes that are fundamentally important across many areas of astrophysics. Whilst Eulerian grid based methods are able to resolve and treat important dynamical instabilities, such as Kelvin-Helmholtz or Rayleigh-Taylor, these processes are poorly or not at all resolved by existing SPH techniques. We show that the reason for this is that SPH, at least in its standard implementation, introduces spurious pressure forces on particles in regions where there are steep density gradients. This results in a boundary gap of the size of the SPH smoothing kernel over which information is not transferred.

Colliding hypersonic flows play a decisive role in many astrophysical objects. In this paper, we look at the idealized model of a 2D plane parallel isothermal slab (CDL) and at symmetric settings, where both flows have equal parameters. We performed a set of high-resolution simulations with upwind Mach numbers, 5 < M_u < 90. We find that the CDL is irregularly shaped and has a patchy and filamentary interior. The size of these structures increases with l_cdl, the extension of the CDL. On average, but not at each moment, the solution is about self-similar and depends only on M_u. We find the root mean square Mach number to scale as M_rms ~ 0.2 M_u. Independent of M_u is the mean density, rho_m ~ 30 rho_u. The fraction f_eff of the upwind kinetic energy that survives shock passage scales as f_eff= 1 - M_rms^(-0.6). This dependence persists if the upwind flow parameters differ from one side to the other of the CDL, indicating that the turbulence within the CDL and its driving are mutually coupled. In the same direction points the finding that the auto-correlation length of the confining shocks and the characteristic length scale of the turbulence within the CDL are proportional. In summary, larger upstream Mach numbers lead to a faster expanding CDL with more strongly inclined confining interfaces relative to the upstream flows, more efficient driving, and finer interior structure relative to the extension of the CDL.

We present results from an extensive set of one- and two-dimensional radiation-hydrodynamic simulations of the supernova core collapse, bounce, and postbounce phases, and focus on the protoneutron star (PNS) spin periods and rotational profiles as a function of initial iron core angular velocity, degree of differential rotation, and progenitor mass. For the models considered, we find a roughly linear mapping between initial iron core rotation rate and PNS spin. The results indicate that the magnitude of the precollapse iron core angular velocities is the single most important factor in determining the PNS spin. Differences in progenitor mass and degree of differential rotation lead only to small variations in the PNS rotational period and profile. Based on our calculated PNS spins, at ~ 200-300 milliseconds after bounce, and assuming angular momentum conservation, we estimate final neutron star rotation periods. We find periods of one millisecond and shorter for initial central iron core periods of below ~ 10 s. This is appreciably shorter than what previous studies have predicted and is in disagreement with current observational data from pulsar astronomy. After considering possible spindown mechanisms that could lead to longer periods we conclude that there is no mechanism that can robustly spin down a neutron star from ~ 1 ms periods to the "injection" periods of tens to hundreds of milliseconds observed for young pulsars. Our results indicate that, given current knowledge of the limitations of neutron star spindown mechanisms, precollapse iron cores must rotate with periods around 50-100 seconds to form neutron stars with periods generically near those inferred for the radio pulsar population.

Using the 2D multi-group, flux-limited diffusion version of the code VULCAN/2D, that also incorporates rotation, we have calculated the collapse, bounce, shock formation, and early post-bounce evolutionary phases of a core-collapse supernova for a variety of initial rotation rates. This is the first series of such multi-group calculations undertaken in supernova theory with fully multi-D tools. We find that though rotation generates pole-to-equator angular anisotropies in the neutrino radiation fields, the magnitude of the asymmetries is not as large as previously estimated. Moreover, we find that the radiation field is always more spherically symmetric than the matter distribution, with its plumes and convective eddies. We present the dependence of the angular anisotropy of the neutrino fields on neutrino species, neutrino energy, and initial rotation rate. Only for our most rapidly rotating model do we start to see qualitatively different hydrodynamics, but for the lower rates consistent with the pre-collapse rotational profiles derived in the literature the anisotropies, though interesting, are modest. This does not mean that rotation does not play a key role in supernova dynamics. The decrease in the effective gravity due to the centripetal effect can be quite important. Rather, it means that when a realistic mapping between initial and final rotational profiles and 2D multi-group radiation-hydrodynamics are incorporated into collapse simulations the anisotropy of the radiation fields may be only a secondary, not a pivotal factor, in the supernova mechanism.

We review some of the reasons for believing that the generic core-collapse supernova is neutrino-driven, not MHD-jet driven. We include a discussion of the possible role of rotation in supernova blast energetics and morphology, and speculate on the origin of Cas A's and SN1987A's ejecta fields. Two "explosive" phenomena may be associated with most core collapses, the neutrino-driven supernova itself and an underenergetic jet-like ejection that follows. The latter may be a magnetic wind that easily penetrates the debris created by the much more energetic supernova. We predict that many core-collapse supernova remnants should have sub-dominant jet-like features. In Cas A, we associate this sub-dominant collimated wind with the "jet/counter-jet" structure observed. We suggest that the actual Cas A explosion itself is at nearly right angles to this jet, along a rotation axis that coincides with the bulk of the ejecta, the iron lobes, and the putative direction of motion of the point source. It may be that when rotation becomes sufficiently rapid that the strong-neutrino-driven-supernova/weak-jet duality switches to a strong-MHD-jet scenario that might be associated with hypernovae, and in some cases GRBs. Finally, we present a calculation using a new 2D multi-group, flux-limited radiation/hydrodynamics code we have recently developed for the simulation of core-collapse supernovae. We discuss the rotation-induced anisotropy in the neutrino radiation field, neutrino heating, and the neutrino flux vectors and speculate on rotation's possible role in the supernova mechanism and in the overall supernova phenomenon.

We have developed a time-dependent, multi-energy-group, and multi-angle (S$_n$) Boltzmann transport scheme for radiation hydrodynamics simulations, in one and two spatial dimensions. The implicit transport is coupled to both 1D (spherically-symmetric) and 2D (axially-symmetric) versions of the explicit Newtonian hydrodynamics code VULCAN. The 2D variant, VULCAN/2D, can be operated in general structured or unstructured grids and though the code can address many problems in astrophysics it was constructed specifically to study the core-collapse supernova problem. Furthermore, VULCAN/2D can simulate the radiation/hydrodynamic evolution of differentially rotating bodies. We summarize the equations solved and methods incorporated into the algorithm and present results of a time-dependent 2D test calculation. A more complete description of the algorithm is postponed to another paper. We highlight a 2D test run that follows for 22 milliseconds the immediate post-bounce evolution of a collapsed core. We present the relationship between the anisotropies of the overturning matter field and the distribution of the corresponding flux vectors, as a function of energy group. This is the first 2D multi-group, multi-angle, time-dependent radiation/hydro calculation ever performed in core collapse studies. Though the transport module of the code is not gray and does not use flux limiters (however, there is a flux-limited variant of VULCAN/2D), it still does not include energy redistribution and most velocity-dependent terms.

Observations of molecular clouds show the existence of starless, dense cores, threaded by magnetic fields. Observed line widths indicate these dense condensates to be embedded in a supersonically turbulent environment. Under these conditions, the generation of magnetic waves is inevitable. In this paper, we study the structure and support of a 1D plane-parallel, self-gravitating slab, as a monochromatic, circularly polarized Alfven wave is injected in its central plane. Dimensional analysis shows that the solution must depend on three dimensionless parameters. To study the nonlinear, turbulent evolution of such a slab, we use 1D high resolution numerical simulations. For a parameter range inspired by molecular cloud observations, we find the following. 1) A single source of energy injection is sufficient to force persistent supersonic turbulence over several hydrostatic scale heights. 2) The time averaged spatial extension of the slab is comparable to the extension of the stationary, analytical WKB solution. Deviations, as well as the density substructure of the slab, depend on the wave-length of the injected wave. 3) Energy losses are dominated by loss of Poynting-flux and increase with increasing plasma beta. 4) Good spatial resolution is mandatory, making similar simulations in 3D currently prohibitively expensive.

We have carried out an extensive set of two-dimensional, axisymmetric, purely-hydrodynamic calculations of rotational stellar core collapse with a realistic, finite-temperature nuclear equation of state and realistic massive star progenitor models. For each of the total number of 72 different simulations we performed, the gravitational wave signature was extracted via the quadrupole formula in the slow-motion, weak-field approximation. We investigate the consequences of variation in the initial ratio of rotational kinetic energy to gravitational potential energy and in the initial degree of differential rotation. Furthermore, we include in our model suite progenitors from recent evolutionary calculations that take into account the effects of rotation and magnetic torques. For each model, we calculate gravitational radiation wave forms, characteristic wave strain spectra, energy spectra, final rotational profiles, and total radiated energy. In addition, we compare our model signals with the anticipated sensitivities of the 1st- and 2nd-generation LIGO detectors coming on line. We find that most of our models are detectable by LIGO from anywhere in the Milky Way.

We present 2D and 3D hydrodynamical simulations of the colliding wind interaction zone in WR+O binaries. It is shown that 3D effects can basically explain certain observed, orbit dependent flux variations. Possible connections between the interior structure of the interaction zone and dust formation are outlined and its stability is re-investigated.

First attempts are made to derive astrophysical implications of the collision of clumped stellar winds from order of magnitude estimates and preliminary numerical simulations. Compared to colliding smooth winds, we find that the most significant differences occur in widely separated systems like WR 140. Clumped winds de-stabilize the wind-wind interaction zone of such systems. Highly compressed, cold knots of WR-wind material can form. Hydrogen rich material is likely to be mixed into these knots by the excited turbulence. Such knots, therefore, are good candidates to form dust. We briefly discuss to what degree our results can be applied to other systems and look at different possibilities for the origin and nature of the inhomogeneities in hot star winds.

Aspects of the wind-dynamics in symbiotic binaries, colliding winds and accretion, are reviewed. Inconsistencies between theory and observations of the hot star wind are discussed. If the hot star wind were governed by CAK theory, nearly all symbiotics would be colliding wind binaries. For the case of colliding winds, 3D hydrodynamical simulations reveal that the matter distribution is spirally shaped. Shock confined high-density shells as well as huge voids are found even in the immediate neighborhood of the stars. Synthetic spectra computed on the basis of different 3D hydrodynamical models suggest observational discrimination between them to be possible. Colliding wind models also provide a link between symbiotics and planetary nebulae. Accretion during some time is a necessary condition for symbiotics to exist. However, there is no proof of whether currently accreting systems show the symbiotic phenomenon. Existing accretion models are inconsistent amongst each other, predicting either extended disks or small, high-density accretion wakes. Synthetic spectra allowing to discriminate between two models do not yet exist.

The colliding winds interaction zone in WR+O binaries is a highly complex environment. In this review we summarize the progress made towards its theoretical understanding during the last years. We review the effect of different physical processes on the interaction zone, among them geometry and orbital motion, radiative forces, thermal conduction, instabilities and turbulence, ionizing radiation, dust formation, clumped winds, magnetic fields, and particle acceleration. Implications with regard to observations are discussed. Subsequently, we proceed to the important question of mutual interaction amongst these processes. Because of the wealth of physical processes involved, numerical simulations are usually mandatory. Finally, we turn to the combined role these processes play for the thermal and ionization properties of the colliding winds interaction zone in WR+O binaries.

Radiative shock waves show a strong cooling instability at temperatures above approximately 2 times 10^5 K. We numerically investigate this instability by simulating different astronomical objects in which colliding flows play an outstanding role: Wind bubbles, supernova remnants, and colliding winds. Computing the flow of each object over a large part of its evolutionary time and resolving all physically relevant scales, we find several phenomenologically different types of this instability. If two smooth flows collide, the instability follows a periodic limit cycle with two modes being important. The connection between the radiative loss function and the mode and type of the overstability is discussed. The collision of non-smooth flows can temporarily result in an aperiodic evolution of the system. After a characteristic relaxation time the instability then becomes periodic again. Such disturbances as well as violent types of the instability can excite oscillations of the thin layer of cold compressed gas downstream of the shock, which in turn can influence the stability of the radiative shock.