We present the results of simulations done with the code HARM-COOL developed in the CTP PAS Warsaw research group over the years 2017-2019. It is based in the original GR MHD scheme proposed by Gammie et al. (2003) for the simulation of Active Galactic Nucleus, but now it has been suited for the engine of a short Gamma Ray Burst event. We compute time-dependent evolution of a black hole accretion disk, in two-dimensional, axisymmetric scheme. The code includes neutrino cooling and accounts for nuclear structure of dense, degenerate matter. Free protons, neutrons, and electron-positron pairs form a neutron-rich, magnetically driven outflow that provides site for subsequent r-process nucleosynthesis. Here the heavy elements up to the Uranium and Gold are synthesized and may contribute to the chemical enrichment of the circum-burst medium. Their radio-active decay will give signal in lower energies in a timescale of weeks-months after the GRB prompt phase. In addition, the magnetic fields are responsible for the launching of ultra-relativistic jets along the rotation axis of the central black hole, according to the well-known Blandford-Znajek mechanism. These jets are sites of variable high energy emission in gamma rays. We find that the magnetic field and the black hole spin account for the observed variability timescales and jet energetics.

Short gamma ray bursts are presumably results of binary neutron star mergers, which lead to the formation of a stellar mass black hole, surrounded by a remnant matter. The strong magnetic fields help collimate jets of plasma, launched along the axis of the black hole rotation. We study the structure and evolution of the accreting plasma in the short GRBs and we model the formation of the base of a relativistic, Poynting-dominated jets. Our numerical models are based on the general relativistic MHD, axisymmetric simulations. We discuss the origin of variability in the GRB jet emission, which timescales are related to the action of the magneto-rotational instability in the accreting plasma. We also estimate the value of a maximum achievable Lorentz factor in the jets produced by our simulations, and reached at the large distances, where the gamma ray emission is produced.

We model the structure and evolution of black hole accretion disks, and their neighboring regions, using numerical simulations. The numerics is governed by the equations of general relativistic magneto-hydrodynamics (GRMHD). In particular, such disks and outflows can be found at the base of very energetic ultra-relativistic jets produced by cosmic explosions, so called gamma-ray bursts (GRBs). Another, more persistent type of the jet phenomena, are blazars, emitted from the centers of galaxies. Long-lasting, detailed computations are essential to properly determine the physics of these explosions, and confront the theoretical models with any potential observables. From the point of view of numerical methods and computational techniques, three ingredients need to be considered. First, the numerical scheme must work in a conservative manner, which is achieved by solving a set of non-linear equations at each time-step, to advance the conserved quantities from one time step to the next. Second, the efficiency of computations intrinsically depends on the code parallelization methods, which may use various techniques. Third, the analysis of results is possible via the post-processing of the computed time-dependent physical quantities, and visualization of the flow properties. This is done via implementing various packages and libraries that are standardized in the field of computational astrophysics and supported by community developers. In the present paper, we discuss the physical picture of the cosmic sources which are modeled using numerical framework. We also describe several technical issues, in the particular context of our own experience with the performance of the GRMHD code which we develop. We also present a suite of performance tests, done on the High-Performance Computer cluster (HPC) in the Center for Mathematical Modeling of the Warsaw University.

The collapsar model was proposed to explain the long-duration gamma-ray bursts (GRBs), while the short GRBs are associated with the mergers of compact objects. In the first case, mainly the energetics of the events is consistent with the proposed progenitor models, while the duration, time variability, as well as the afterglow emission may shed some light on the detailed properties of the collapsing massive stars. In the latter case, the recent discovery of the binary neutron star (NS-NS) merger in the gravitational wave observation made by LIGO (GW170817), and the detection of associated electromagnetic counterparts, for the first time gave a direct proof of the NS-NS merger being a progenitor of a short GRB. In general, all GRBs are believed to be powered by accretion through a rotationally supported torus, or by fast rotation of a compact object. For long ones, the rotation of the progenitor star is a key property in order to support accretion over relatively long activity periods, and also to sustain the rotation of the black hole itself. The latter is responsible for ejection of the relativistic jets, which are powered due to the extraction of the BH rotational energy, mitigated by the accretion torus and magnetic fields. The jets must break through the stellar envelope though, which poses a question on the efficiency of this process. Similar mechanisms of powering the jet ejection may act in short GRBs, which in this case may freely propagate through the interstellar medium. The power of the jets launched from the rotating black hole is at first associated mostly with the magnetic Poynting flux, and then at large distances it is transferred to the kinetic and finally radiative energy of the expanding shells.

We study the structure and evolution of the accreting plasma in gamma ray burst central engines. The models are based on the general relativistic MHD simulations. The nuclear equation of state adequate for dense and degenerate plasma, is incorporated to the numerical scheme instead of a simple polytropic EOS. Plasma is cooled by neutrinos and energy is extracted from the rotating BH by magnetic fields. We discuss our results in the frame of the observable GRBs, and speculate about the origin of variability in the GRB jet emission. We also discuss the possible contribution from the heavy elements synthesized in winds from GRB engines to the kilonova emission.

We construct and analyze a model of the relativistic steady-state magnetohydrodynamic (MHD) rarefaction that is induced when a planar symmetric flow (with one ignorable Cartesian coordinate) propagates under a steep drop of the external pressure profile. Using the method of self-similarity we derive a system of ordinary differential equations that describe the flow dynamics. In the specific limit of an initially homogeneous flow we also provide analytical results and accurate scaling laws. We consider that limit as a generalization of the previous Newtonian and hydrodynamic solutions already present in the literature. The model includes magnetic field and bulk flow speed having all components, whose role is explored with a parametric study.

Relativistic jets associated with long/soft gamma-ray bursts are formed and initially propagate in the interior of the progenitor star. Because of the subsequent loss of their external pressure support after they cross the stellar surface, these flows can be modeled as moving around a corner. A strong steady-state rarefaction wave is formed, and the sideways expansion is accompanied by a rarefaction acceleration. We investigate the efficiency and the general characteristics of this mechanism by integrating the steady-state, special relativistic, magnetohydrodynamic equations, using a special set of partial exact solutions in planar geometry (r self-similar with respect to the "corner"). We also derive analytical approximate scalings in the ultrarelativistic cold/magnetized, and hydrodynamic limits. The mechanism is more effective in magnetized than in purely hydrodynamic flows. It substantially increases the Lorentz factor without much affecting the opening of the jet; the resulting values of their product can be much grater than unity, allowing for possible breaks in the afterglow light curves. These findings are similar to the ones from numerical simulations of axisymmetric jets by Komissarov et al and Tchekhovskoy et al, although in our approach we describe the rarefaction as a steady-state simple wave and self-consistently calculate the opening of the jet that corresponds to zero external pressure.