Skip to main content
European Commission logo print header

CAMAP: Computer Aided Modeling for Astrophysical Plasmas

Final Report Summary - CAMAP (CAMAP: Computer Aided Modeling for Astrophysical Plasmas)

Astrophysical plasma is probably the most exciting, complex, physically entangled state of matter in Nature. This is particularly true if the magnetic fields providing the coupling for the plasma constituents become dynamically important. Additional challenges for the physical modelling of astrophysical plasma pop up if its constituents are relativistic or if the plasma moves at bulk speeds very close to the speed of light in vacuum, c, giving rise to a complete subfield of Astrophysics in it own right: Relativistic Plasma Astrophysics. Currently, Relativistic Plasma Astrophysics has become a subject of great interest because it pushes towards the limit our basic theoretical understanding of Nature: speeds very close to c, strong gravitational fields in the vicinity of compact objects, uncertain equations of state for highly compressed mater beyond nuclear matter density, large dynamical ranges in order to couple regimes from Newtonian to ultrarelativistic, the solution of constrained hyperbolic systems of equations in General Relativity, etc. Furthermore, the recent detection of gravitational waves by the Ligo interferometers on September 14h 2015, 100 years after the prediction of Einstein, has been a major breakthrough in Relativistic Astrophysics.

In our project we have obtained an almost integral view of the process of formation of relativistic flows resulting in long gamma-ray bursts (GRBs). Our main focus has been to begin from state-of-the-art stellar evolution models, let them collapse and obtain the consistent evolution of the remnant and any potential outflow. In agreement with prior expectations, the formation of the central engine of a long GRB is like a theater play in three acts lasting a few seconds. On the first act, a compact remnant resulting from the stellar core implosion forms. This act takes a fraction of a second. During the second act, the central remnant contracts as it frees part of its entropy. A small fraction of the remnant entropy is transferred to the shells of matter surrounding the remnant. Whether or not the transferred entropy is enough to halt the accretion of matter onto the central compact object depends on a delicate interplay between neutrino physics, the rotational profile in the collapsing star and the action of huge magnetic fields build up in the course of the core collapse (and afterwards). In the final and definitive act, the central object, a hypermassive neutron star (also known as a proto-neutron star; PNS hereafter), finally undergoes a further dramatic collapse to form a black hole (BH) or the remnant cools down secularly forming an extremely magnetized PNS. Adding to the main action of the former theater play, we have observed that the system releases, on the one hand, large-amplitude gravitational waves as well as an outflow consisting of ejecta moving at different speeds from subrelativistic (the supernova ejecta) to ultrarelativistic (the jets which may release the GRB). We have studied the conditions that massive stars with low metallicity must possess in order to produce either BHs or hypermassive PNSs. Whether the remnant is one or the other has a profound implication on the mechanism to produce GRBs. We find that sufficiently large angular momentum per unit mass in the stellar core added to “moderate” or “large” magnetic fields result in proto-magnetars (i.e. extremely magnetized and massive neutron stars with rotational periods of the order of 1 millisecond and magnetic field strengths beyond 100 million MGauss –equivalent to 10 billion times the Earth surface magnetic field). On the other extreme, either smaller initial rotational rates in the progenitor or smaller initial magnetic fields result into collapsars (i.e. fast rotating BHs surrounded by thick accretion discs in the stellar core).

Regardless of the outcome of the stellar core collapse, the most direct evidence of the type of central engine feeding a GRB comes from the analysis of the gravitational wave signature left by the system dynamics. Thus, we have thoroughly analyzed our numerical models to predict the gravitational waveforms. Our simulations show that the phase during which the PNS survives before collapsing to a BH is optimal for gravitational wave emission. The emitted high-amplitude waves last for several seconds and show a remarkable quasi-periodicity associated with the violent PNS dynamics, namely produced by episodes of convection and the subsequent nonlinear development of the standing-accretion shock instability (SASI). By analysing the spectrogram of our simulations we are able to identify the frequencies associated with the presence of gravity modes (g-modes) and with the SASI motions at the PNS surface. We note that the gravitational waves emitted reach large enough amplitudes to be detected with third-generation of gravitational wave detectors such as the Einstein Telescope within the Virgo Cluster volume at rates 0.1 events per year.

Complementary to the gravitational wave emission and acting as a proxy for its search in gravitational wave detectors, the central engines of long GRBs generate outflows whose synthetic characteristic signature we have tried to characterize. For that, we need to wait until the outflow, generated at scales of ~100 km, becomes optically thin, so that radiation may escape from its interior. This happens at huge distances from the central engine, namely, beyond scales of about a billion kilometers (roughly, the radius of our Sun). At these distances the outflow has already broken out of the surface of the progenitor star, whose outer layers remain unaware of the fatidic fate the stellar core has just undergone. Estimating such a signal helps us to compare with existing observations, especially the ones that have been attributed to the breakout of such outflows through the surface of the stellar progenitor, or as the outflows interact with a complex circumstellar medium. To be more precise, we have argued that a large fraction of long-duration GRBs should take place in massive stellar clusters where the circumburst medium is much more complicated than commonly assumed. These media possess a non-trivial structure as a result of the interaction of stellar winds arising from close by massive stars, some of which may be the progenitors of GRB jets. We have indeed estimated the light curve and spectra resulting in such scenario for various degrees of collimation of the ultrarelativistic ejecta generated in the stellar progenitor. For uncollimated ejecta the afterglow light curve shows a flattening followed by a shallow break on a time-scale from hours up to a week after the burst, which results from their propagation through the shocked wind region. If the ejecta is collimated, the light curve break may become very pronounced with the post-break temporal decline as steep as t^{-5}. We estimate that the inverse Compton up scattering of ultraviolet photons from the nearby star off energetic electrons in the relativistic ejecta leads to a bright GeV afterglow flare that may be detectable by the Fermi satellite.

Besides in stellar progenitors of GRBs, magnetic fields may have an influence on the broader class of massive stars, which undergo core-collapse and supernova explosions. Our group has studied the amplification of magnetic fields in cores during and after collapse in order to determine whether ordinary seed magnetic fields can grow sufficiently to affect the dynamics of the stellar core and potentially trigger an explosion (or any other kind of ejecta). In some of our models, which do not include rotation, field amplification is mostly due to compression and the turbulent flow created by hydrodynamic instabilities. If the initial field is already strong enough, we find that the magnetic field is able to modify these flows and contribute to launching an explosion. In other models, where rotation of the stellar core is included, we find that the magneto-rotational instability (MRI) is the main amplification agent of the preexisting magnetic field in the imploded core. The growth of the magnetic field is saturated by the action of Kelvin-Helmhotz parasitic instabilities feeding off the channel structures developed in the course of the magnetic field amplification. Remarkably, our detailed numerical models suggest that the MRI channel modes provide only a fraction of the amplification that magnetic fields need to reach magnetar-like strengths. The study of this saturation phase is extremely challenging from the numerical simulation point of view and has only been possible thanks to the extremely high accuracy of our numerical tools. These tools have been developed in our group, calibrated and cross validated with the codes of other groups in such a way that we can ensure that our results correspond to sound physics rather than to numerical artifacts. Actually, we have developed a set of test-beds for Eulerian MHD codes, which may probably set the standards for the community to cross-compare the accuracy of their own tools.

The story of the imploded core of a massive star, if it does not finally result into a BH, does not finish right after collapse. Ten to hundred seconds after the formation of a PNS it cools down and contracts, settling into an extremely compact neutron star. These neutron stars are observed typically because they become strong radio-to-X-ray emitters forming a variety of fascinating objects from pulsars to magnetars. The observed emission is certainly related to the existence of gigantic magnetic fields inside of the neutron star and in its surrounding magnetosphere. Whether the observed magnetic fields in cold neutron stars correspond to the ones build up during the collapse of the progenitor stellar core or there are ulterior episodes of matter accretion which modulate the strength of the magnetospheric fields has become a matter of active scientific debate. Actually, there are observations of several neutron stars with significantly lower values of the magnetic field than the average radio-pulsar population. This intriguing conundrum has motivated our group to study the so-called hidden magnetic field scenario, where the accretion of the fallback of the supernova debris on to the neutron star is responsible for the submergence of the field and its apparently low value. The main conclusion of our work is that typical magnetic fields of a few thousands of billions of Gauss can be buried by accreting less than one hundredth of a solar mass over the surface of the neutron star. Thus, we have predicted that anomalously weak magnetic fields should be common in “very young” neutron stars, i.e. in neutron stars created less than a few millennia ago.

Beyond the questions raised by anomalously low magnetic fields in some young neutron stars, it is also found that in a fraction of the observed strongly magnetized neutron stars there exist quasi-periodic oscillations (QPOs), which happen during the giant flares of soft gamma repeaters (a special class of magnetars). The interest of these QPOs resides in the fact that their frequency depends on the equation of state of nuclear matter, on the mass, on the magnetic field, on the superfluidity properties of the neutron star core, etc. The added complexity of evolved magnetars (compared to our magneto-fluid models of PNSs resulting from the stellar core collapse) is that they develop a solid crust and an interior with superfluid neutrons and all of these ingredients must be carefully accounted for with a General Relativistic treatment. We have modelled these QPOs as torsional oscillations of the structure of magnetars (accompanied by the corresponding oscillations of their magnetic field topology). We have concluded that the oscillations of the magnetar crust cannot be responsible for the observed QPOs, regardless of the structure of the magnetic field. Going further, we have suggested that there is, at least, one superfluid species in magnetar cores. In this realm, we have also investigated the coupling of magnetar core oscillations with their surrounding magnetosphere, where the QPO emission processes take place as a result of the resonant cyclo-synchrotron scattering of thermal photons with magnetospheric currents. Indeed, we have computed the structure of twisted neutron star magnetospheres and identified the maximum twist that they can support. At this maximum twist instabilities are expected to appear, which, in our opinion, may be linked to observed flare activity on magnetars.

Furthermore, in this project we have been able to find the missing theoretical link between the process of merging of neutron stars and the formation of the central engine of gamma-ray bursts, one of the most luminous events in the Universe since the Big-Bang. For this discovery, it has also been crucial to understand very well how the process of amplification of the pre-existing magnetic fields in neutron stars works. The action of matter compression and the growth of hydro-magnetic instabilities makes it possible that “small” seed fields, with strength of the order of a hundred thousand times the strongest neodymium-iron-boron magnets produced by humans, can be amplified by, at least, a million times more, to reach strengths beyond 1 billion MGauss. The processes involved in the field amplification are so quick (they may happen in a fraction of a second), and so non-linear that only state-of-the-art numerical simulations can provide satisfactory answers to the problems at hand.

Our group has also developed one of the most advanced electron transport schemes for relativistic flows, which is implemented in the SPEV code. A special version of the SPEV code has been developed for the purpose of studying blazar variability and spectra. In a standard picture the jet in blazars is pointing almost directly towards the Earth and the variability is due to the collision of different parts of the jet due to the velocity variations. We model the collisions as interactions of cylindrical shells made of plasma. We first determined the amount of energy dissipated at the shock fronts caused by the collision, and then compute the time- and frequency-dependent emission by the electrons accelerated by the shocks. Even in such a simplified picture the calculation is extremely demanding and requires the use of supercomputers. The reason is that we need to compute the radiation at frequencies ranging from infrared to gamma-rays at TeV energies, and we must consider different radiative processes responsible for radiation at different wavelengths. The goal of our study was to determine the signature of the presence of magnetized plasma on blazar variability. We concluded that a significant part of the blazar phenomenology ("blazar sequence") could be explained by varying the magnetization and the velocity differences between different parts (shells) in the blazar jet. Going beyond the standard “shock-in-jet” model for blazars, we have found that more realistic electron energy distributions arising from shock acceleration yield an imprint in on the observed spectra. To account for the emission of these more realistic electron energy distributions extending to very low energies, it has been necessary to go beyond the standard synchrotron approximation, valid in the ultrarelativistic regime. Currently, our group possesses a state-of-the-art numerical code to compute the magneto-bremsstrahlung emission from electrons of any energy (also very low energy) and radiating at arbitrary frequencies (also including very low radio frequencies). The imprint that a broad thermal tail of the electron energy distribution has on the blazar spectrum turns out to be notable, especially regarding the shape, frequency extension and flux density reached at the inverse Compton peak characteristic of blazars. We foresee that a careful comparison of our synthetic spectra with actual blazar observations may reveal the thermal content of the plasma blown in the relativistic jets feeding blazars.

The complex machinery developed in the SPEV code has also been applied to a field, which is very akin to that of GRB afterglow, in many aspects, namely the recently observed tidal disruption events (TDEs). Indeed, as we model them, TDEs are oversized versions of GRB afterglows, happening at much smaller Lorentz factors. Indeed, as in the case of GRB afterglows, which are affected by the structure of the circumstellar medium of their stellar progenitors, in TDEs the circumnuclear medium (CNM) surrounding the supermassive BH where they are produced also shapes their light curves and spectra. We have been able of quantifying the link between the CNM density, jet energy and the radio peak flux and time, opening the door to constraining both the jet energies and the CNM properties using follow-up radio observations on the scale of months or years after the tidal disruption.

Loyal to our tradition of confronting theoretical and numerical models to actual observations, we have broadly collaborated with worldwide groups of astronomers. In particular, we have helped to understand the mechanism of a new class of gamma-ray burst characterized by an unusual thermal component (blackbody-dominated GRBs; BB-GRBs), being our results published in Nature and in several other mass media. It has turned out that the peculiarities of BB-GRBs are likely linked to progenitor systems formed by the merger of a neutron star and a Helium star. Strikingly, the central engine that results from this kind of mergers resembles very much that of the collapsar model, namely, a fast spinning BH surrounded by copious amounts of matter that can be accreted over hundreds to thousands of seconds. Beyond the subtleties of the progenitor system, what turns out to be most decisive to suppress a standard afterglow in these events is the complex and high-density circumstellar medium that results from the ejection of the outer envelope of the secondary star of the binary system. Likewise, we have performed key numerical simulations to unveil the controversial nature of the observations of the supernova remnant associated to Vela Jr., which point to the possibility that the supernova that produced this remnant happened in the 13th century, and would have been seen in its whole glory from the ancient Great Zimbabwe.

Also, we have helped to understand the radio observations of the supermassive BH in the center of our own Milky Way. These observations indicate the existence of a non-thermal population of electrons, whose origin is controversial. In collaboration with Dr. Ian Christie we have proposed a model to explain the elusive origin of the non-thermal electrons: a pulsar transits through the accretion disk of our Galaxy’s supermassive BH and relativistic charged particles generated by the resulting shock are injected into the disk. These particles sustain their emission until they are accreted by the BH. Analytical and numerical calculations we have done show that periodic transitions of a young pulsar with a spin up luminosity ten times larger than that of our Sun within
0.1 pc of the BH can provide an adequate supply of non-thermal charged particles that emit in the radio band through synchrotron cooling.

Finally, we remark that our group has developed world-leading numerical codes that incorporate the most advanced techniques in the simulation of astrophysical scenarios where the simultaneous action of General Relativistic gravity, magnetohydrodynamics, microphysics and transport of energy challenge our current understanding of the physics of the most extreme objects of the Universe. The generation of new algorithms to deal with physical and technical problems intrinsic to the numerical modelling of the astrophysical scenarios at hand has been necessary for the success of our project. We outline the fact that these new algorithms have resulted from the multidisciplinary character of our team, which combines astrophysicists, applied mathematicians and computer scientists. A couple of examples of these new algorithms are a new solver for systems of partial differential equations of elliptic type and the new schemes for the stable integration of hyperbolic systems of partial differential equations which may include stiff source terms. The solver for systems of partial differential equations results from a generalization of the Jacobi method, which is so simple that allows for a straightforward massive parallelization, making it perfectly suited for exaescale supercomputers as well as for GPU based high-performance computing systems.