Numerical methods for radiative and ideal relativistic hydrodynamics applied to the study of gamma-ray bursts

This thesis is devoted to the application of high-resolution numerical methods for relativistic hydrodynamics (RHD) to the study of gamma-ray bursts (GRBs), as well as to the development of new schemes able to describe radiative transfer in relativistic magnetized and unmagnetized flows. On one side, we have performed RHD simulations of relativistic plasma outbursts within the binary-driven hypernova model, developed throughout the last years in the International Center of Relativistic Astrophysics Network (ICRANet). This model is based on the so-called induced gravitational collapse scenario, proposed to explain the observed temporal coincidence of GRBs and supernovae (SN) of type Ic. This scenario considers a carbon-oxigen star (CO core) forming a tight binary system with a companion neutron star (NS). When the collapse of the CO core produces a type Ic SN, part of the ejected material is accreted by the NS, which in turn collapses and forms a black hole (BH). It has been proposed, although the details of this process are a matter of current research, that this collapse creates an optically thick electron-positron plasma around the BH that expands due to its own internal pressure and originates a GRB. Our work in this context has focused on the description of such expanding plasma and its interaction with the surrounding SN ejecta, for which we have followed a hydrodynamical approach using the open-source code PLUTO. This allowed us to study this process in high-density regions that had not been explored thus far, and to perform consistency checks of the model taking into account both theoretical and observational constraints such as the system’s size, the initial plasma energy, the observed timing and the Lorentz factor of the outbursts. Three different scenarios are here considered: (I) the expansion of the plasma in low-density regions, proposed to produce most of the GRB emission in the prompt phase; (II) a model in which X-ray flares are produced due to the breakout of shocks created when the plasma interacts with high-density regions of the SN ejecta; and (III) a model for the emission of secondary bursts due to the creation of reflected waves caused by the same interaction. The second part of this thesis is devoted to the main part of our work, which consists in the development of a numerical code for radiative transfer integrated in PLUTO. Our implementation is able to solve the equations of relativistic radiation magnetohydrodynamics (Rad-RMHD) under the so-called M1 closure, which allows the radiation transport to be handled in both the free-streaming and diffusion limits. Since we use frequency-averaged opacities, this approach is unable to describe frequency-dependent phenomena; instead, the main focus is put on the transport of total energy and momentum. To avoid numerical instabilities arising due to the possibly large timescale disparity caused by the radiation–matter interaction terms, the Rad-RMHD equations are integrated following implicit–explicit (IMEX) schemes. In this way, interaction terms are integrated implicitly, whereas transport and all of the remaining source terms are solved explicitly by means of the same Godunov-type solvers included in PLUTO. Among these, we have introduced a new Harten–Lax–van Leer–contact (HLLC) solver for optically thin radiation transport. The code is suitable for multidimensional computations in Cartesian, spherical, and cylindrical coordinates using either a single processor or parallel architectures. Adaptive grid computations are also made possible by means of the CHOMBO library. We explain in this work the implementation of all of these methods, after which we show the code’s performance in several problems of radiative transfer in magnetized and unmagnetized flows. We pay particular attention to the behavior of the solutions in the free-streaming and diffusion limits, and show the efficiency and scalability properties of the code as compared with its usual nonradiative implementation. Finally, we show an application of this code to the mentioned model for X-ray flares

No self-shadowing instability in 2D radiation-hydrodynamical models of irradiated protoplanetary disks

Theoretical models of protoplanetary disks including stellar irradiation often show a spontaneous amplification of scale height perturbations, produced by the enhanced absorption of starlight in enlarged regions. In turn, such regions cast shadows on adjacent zones that consequently cool down and shrink, eventually leading to an alternating pattern of overheated and shadowed regions. Previous investigations have proposed this to be a real self-sustained process, the so-called self-shadowing or thermal wave instability, which could naturally form frequently observed disk structures such as rings and gaps, and even potentially enhance the formation of planetesimals. All of these, however, have assumed in one way or another vertical hydrostatic equilibrium and instantaneous radiative diffusion throughout the disk. In this work we present the first study of the stability of accretion disks to self-shadowing that relaxes these assumptions, relying instead on radiation-hydrodynamical simulations. We first construct hydrostatic disk configurations by means of an iterative procedure and show that the formation of a pattern of enlarged and shadowed regions is a direct consequence of assuming instantaneous radiative diffusion. We then let these solutions evolve in time, which leads to a fast damping of the initial shadowing features in layers close to the disk surface. These thermally relaxed layers grow towards the midplane until all temperature extrema in the radial direction are erased in the entire disk. Our results suggest that radiative cooling and gas advection at the disk surface prevent a self-shadowing instability from forming, by damping temperature perturbations before these reach lower, optically thick regions.Comment: 26 pages, 19 figures. Accepted for publication in ApJ. Complementary videos can be found in https://youtu.be/RT8IFe8W13

Vertical shear instability in two-moment radiation-hydrodynamical simulations of irradiated protoplanetary disks II. Secondary instabilities and stability regions

The vertical shear instability (VSI) is a hydrodynamical instability likely to produce turbulence in the dead zones of protoplanetary disks. Various aspects of this instability remain to be understood, including the disk regions where it can operate and the physical phenomena leading to its saturation. In this work, we studied the growth and evolution of secondary instabilities parasitic to the VSI, examining their relation with its saturation in axisymmetric radiation-hydrodynamical simulations of protoplanetary disks. We also constructed stability maps for our disk models, considering temperature stratifications enforced by stellar irradiation and radiative cooling and incorporating the effects of dust-gas collisions and molecular line emission. We found that the flow pattern produced by the interplay of the axisymmetric VSI modes and the baroclinic torque forms bands of nearly uniform specific angular momentum. In the high-shear regions in between these bands, the Kelvin-Helmholtz instability (KHI) is triggered. The significant transfer of kinetic energy to small-scale eddies produced by the KHI and possibly even the baroclinic acceleration of eddies limit the maximum energy of the VSI modes, likely leading to the saturation of the VSI. A third instability mechanism, consisting of an amplification of eddies by baroclinic torques, forms meridional vortices with Mach numbers up to $\sim 0.4$. Our stability analysis suggests that protoplanetary disks can be VSI-unstable in surface layers up to tens of au for reasonably high gas emissivities, even in regions where the midplane is stable. This picture is consistent with current observations of disks showing thin midplane millimeter-sized dust layers while appearing vertically extended in optical and near-infrared wavelengths.Comment: Accepted for publication in Astronomy & Astrophysic

Protostellar disks subject to infall: a one-dimensional inviscid model and comparison with ALMA observations

A new one-dimensional, inviscid, and vertically integrated disk model with prescribed infall is presented. The flow is computed using a second-order shock-capturing scheme. Included are vertical infall, radial infall at the outer radial boundary, radiative cooling, stellar irradiation, and heat addition at the disk-surface shock. Simulation parameters are chosen to target the L1527 IRS disk which has been observed using ALMA (Atacama Large Millimeter Array). The results give an outer envelope of radial infall and $u_\phi \propto 1/r$ which encounters a radial shock at $r_\mathrm{shock} \sim 1.5\ \times$ the centrifugal radius ($r_\mathrm{c}$) across which the radial velocity is greatly reduced and the gas temperature rises from a pre-shock value of $\approx 25$ K to $\approx 180$ K over a spatially thin region calculated using a separate shock structure code. At $r_\mathrm{c}$, the azimuthal velocity $u_\phi$ transitions from being $\propto 1/r$ to being nearly Keplerian. These results qualitatively agree with recent ALMA observations which indicate a radial shock where SO is sublimated as well as a transition from a $u_\phi \sim 1/r$ region to a Keplerian inner disk. However, in one set of observations, the position-velocity map of cyclic-C$_3$H$_2$, together with a certain ballistic maximum velocity relation, has suggested that the radial shock coincides with a ballistic centrifugal barrier, which places the shock at $r_\mathrm{shock} = 0.5 r_\mathrm{c}$, i.e, inward of $r_\mathrm{c}$, rather than outward as given by our simulations. It is argued that radial velocity plots from previous magnetic rotating-collapse simulations also indicate that the radial shock is located outward of $r_\mathrm{c}$. The discrepancy with observations is analyzed and discussed, but remains unresolved.Comment: Originally, we incorrectly took Semenov etal. opacities to be m$^2$ per gm of dust rather than gas. Thus our opacities were too low by a factor of 100. Making the correction reduced the temperature across the shock but left velocities and densities nearly unchanged. To account for SO sublimation in L1527 observed by ALMA, we performed a separate 1D shock calculation including non-LTE effect

Vertical shear instability in two-moment radiation-hydrodynamical simulations of irradiated protoplanetary disks I. Angular momentum transport and turbulent heating

We studied the linear and nonlinear evolution of the Vertical Shear Instability (VSI) in axisymmetric models of protoplanetary disks, focusing on the transport of angular momentum, the produced temperature perturbations, and the applicability of local stability conditions. We modeled the gas-dust mixture via high-resolution two-moment (M1) radiation-hydrodynamical simulations including stellar irradiation with frequency-dependent opacities. We found that, given sufficient depletion of small grains (with a dust-to-gas mass ratio of $10\%$ of our nominal value of $10^{-3}$ for $<0.25$ $\mu$m grains), the VSI can operate in surface disk layers while being inactive close to the midplane, resulting in a suppression of the VSI body modes. The VSI reduces the initial vertical shear in bands of approximately uniform specific angular momentum, whose formation is likely favored by the enforced axisymmetry. Similarities with Reynolds stresses and angular momentum distributions in 3D simulations suggest that the VSI-induced angular momentum mixing in the radial direction may be predominantly axisymmetric. The stability regions in our models are well explained by local stability criteria, while the employment of global criteria is still justifiable up to a few scale heights above the midplane, at least as long as VSI modes are radially optically thin. Turbulent heating produces only marginal temperature increases of at most $0.1\%$ and $0.01\%$ in the nominal and dust-depleted models, respectively, peaking at a few (approximately three) scale heights above the midplane. We conclude that it is unlikely that the VSI can, in general, lead to any significant temperature increase since that would either require it to efficiently operate in largely optically thick disk regions or to produce larger levels of turbulence than predicted by models of passive irradiated disks.Comment: Accepted for publication in Astronomy & Astrophysic

Neutrino production from proton-proton interactions in binary-driven hypernovae

We estimate the neutrino emission from the decay chain of the $\pi$-meson and $\mu$-lepton, produced by proton-proton inelastic scattering in energetic ($E_{\rm iso}\gtrsim 10^{52}$~erg) long gamma-ray bursts (GRBs), within the type I binary-driven hypernova (BdHN) model. The BdHN I progenitor is \textcolor{red}{a} binary system composed of a carbon-oxygen star (CO$_{\rm core}$) and a neutron star (NS) companion. The CO$_{\rm core}$ explosion as supernova (SN) triggers a massive accretion process onto the NS. For short orbital periods of few minutes, the NS reaches the critical mass, hence forming a black hole (BH). Recent numerical simulations of the above scenario show that the SN ejecta becomes highly asymmetric, creating a \textit{cavity} around the newborn BH site, due to the NS accretion and gravitational collapse. Therefore, the electron-positron ($e^{\pm}$) plasma created in the BH formation, during its isotropic and self-accelerating expansion, engulfs different amounts of ejecta baryons along different directions, leading to a direction-dependent Lorentz factor. The protons engulfed inside the high-density ($\sim 10^{23}$~particle/cm$^3$) ejecta reach energies in the range $1.24\lesssim E_p\lesssim 6.14$ GeV and interact with the unshocked protons in the ejecta. The protons engulfed from the low density region around the BH reach energies $\sim 1$ TeV and interact with the low-density ($\sim1$~particle/cm$^3$) protons of the interstellar medium (ISM). The above interactions give rise, respectively, to neutrino energies $E_{\nu}\leq 2$ GeV and $10\leq E_{\nu}\leq 10^3$ GeV, and for both cases we calculate the spectra and luminosity.Comment: 19 pages, 22 figures, 2 tables, re-submitted to Physical Review Letters

Thermal instabilities in accretion disks II: Numerical Experiments for the Goldreich-Schubert-Fricke Instability and the Convective Overstability in disks around young stars

The linear stability analysis of a stratified rotating fluid (see paper I) showed that disks with a baroclinic stratification under the influence of thermal relaxation will become unstable to thermal instabilities. One instability is the Goldreich-Schubert-Fricke instability (GSF), which is the local version of the Vertical Shear Instability (VSI) and the other is a thermal overstability, the Convective Overstability (COS). In the present paper we reproduce the analytic predicted growth rates for both instabilities in numerical experiments of small axisymmetric sections of vertically isothermal disks with a radial temperature gradient, especially for cooling times longer than the critical cooling time for VSI. In this cooling time regime our simulations reveal the simultaneous and independent growth of both modes: COS and GSF. We consistently observe that GSF modes exhibit a faster growth rate compared to COS modes. Near the midplane, GSF modes eventually stop growing, while COS modes continue to grow and ultimately dominate the flow pattern. Away from the midplane, we find GSF modes to saturate, when bands of constant angular momentum have formed. In these bands we observe the formation and growth of eddies driven by the baroclinic term, further enhancing the velocity perturbations. In geophysics this effect is known as horizontal convection or sea-breeze instability. Three-dimensional simulations will have to show whether similar effects will occur when axisymmetry is not enforced. Our local simulations help to reveal the numerical resolution requirements to observe thermal instabilities in global simulations of disks around young stars.Comment: ApJ, in press: 27 pages, 18 figures, 4 Movie

GRB 110731A within the IGC paradigm

Bright gamma-ray burst (GRB) 110731A was simultaneously observed by Fermi and Swift observatories, with a follow up optical observation which inferred the redshift of z = 2.83. Thus, available data are spanning from optical to high energy (GeV) emission. We analyze these data within the induced gravitational collapse (IGC) paradigm, recently introduced to explain temporal coincidence of some long GRBs with type Ic supernovae. The case of binary-driven hypcrnova (BdHN) assumes a close system, which starts as an evolved core - neutron star binary. After the core-collapse event, the new NS - black hole system is formed, emitting the GRB in the process. We performed the time-resolved and time-integrated analysis of the Fermi data. Preliminary results gave isotropic energy Eiso = 6.05 × 1053 erg and the total P-GRB energy of Ep–GRB = 3.7 × 1052 erg. At transparency point we found a Lorentz factor Γ ~ 2.17 × 103 laboratory radius of 8.33 x 1013 cm, P-GRB observed temperature of 168 keV and a baryon load B = 4.35 × 10-4. Simulated light-curve and prompt emission spectra showed the average circum burst medium density to be n ~ 0.03 particles per cm3. We reproduced the X-ray light-curve within the rest-frame of the source, finding the common late power-law behavior, with α = –1.22. Considering these results, we interpret GRB 110731A as a member of a BdHNe group