RUNKO is a new open-source plasma simulation framework implemented in C++ and PYTHON. It is designed to function as an easy-to-extend general toolbox for simulating astrophysical plasmas with different theoretical and numerical models. Computationally intensive low-level kernels are written in modern C++ taking advantage of polymorphic classes, multiple inheritance, and template metaprogramming. High-level functionality is operated with PYTHON scripts. The hybrid program design ensures good code performance together with ease of use. The framework has a modular object-oriented design that allows the user to easily add new numerical algorithms to the system. The code can be run on various computing platforms ranging from laptops (shared-memory systems) to massively parallel supercomputer architectures (distributed-memory systems). The framework supports heterogeneous multiphysics simulations in which different physical solvers can be combined and run simultaneously. Here, we showcase the framework’s relativistic particle-in-cell (PIC) module by presenting (i) 1D simulations of relativistic Weibel instability, (ii) 2D simulations of relativistic kinetic turbulence in a suddenly stirred magnetically-dominated pair plasma, and (iii) 3D simulations of collisionless shocks in an unmagnetized medium.

We perform 2D and 3D kinetic simulations of reconnection-mediated turbulent flares in a magnetized electron-positron plasma, with weak and strong radiative cooling. Such flares can be generated around neutron stars and accreting black holes. We focus on the magnetically dominated regime where tension of the background magnetic field lines exceeds the plasma rest-mass density by a factor σ₀ > 1. In the simulations, turbulence is excited on a macroscopic scale l0, and we observe that it develops by forming thin, dynamic current sheets on various scales. The deposited macroscopic energy dissipates by energizing thermal and nonthermal particles. The particle energy distribution is shaped by impulsive acceleration in reconnecting current sheets, gradual stochastic acceleration, and radiative losses. We parameterize radiative cooling by the ratio Α of light-crossing time l₀/c to a cooling timescale, and study the effect of increasing Α on the flare. When radiative losses are sufficiently weak, Α>σ₀^-1, the produced emission is dominated by stochastically accelerated particles, and the radiative power depends logarithmically on Α. The resulting radiation spectrum of the flare is broad and anisotropic. In the strong-cooling regime, Α>σ₀^-1, stochastic acceleration is suppressed, while impulsive acceleration in the current sheets continues to operate. As Α increases further, the emission becomes dominated by thermal particles. Our simulations offer a new tool to study particle acceleration by turbulence, especially at high energies, where cooling competes with acceleration. We find that the particle distribution is influenced by strong intermittency of dissipation, and stochastic acceleration cannot be described by a universal diffusion coefficient.

The theory governing the strong nuclear force-quantum chromodynamics-predicts that at sufficiently high energy densities, hadronic nuclear matter undergoes a deconfinement transition to a new phase of quarks and gluons(1). Although this has been observed in ultrarelativistic heavy-ion collisions(2,3), it is currently an open question whether quark matter exists inside neutron stars(4). By combining astrophysical observations and theoretical ab initio calculations in a model-independent way, we find that the inferred properties of matter in the cores of neutron stars with mass corresponding to 1.4 solar masses (M-circle dot) are compatible with nuclear model calculations. However, the matter in the interior of maximally massive stable neutron stars exhibits characteristics of the deconfined phase, which we interpret as evidence for the presence of quark-matter cores. For the heaviest reliably observed neutron stars(5,6) with mass M approximate to 2M(circle dot), the presence of quark matter is found to be linked to the behaviour of the speed of sound c(s) in strongly interacting matter. If the conformal bound cs2 <= 1/3 (ref. (7)) is not strongly violated, massive neutron stars are predicted to have sizable quark-matter cores. This finding has important implications for the phenomenology of neutron stars and affects the dynamics of neutron star mergers with at least one sufficiently massive participant.

When the accretion disc around a weakly magnetised neutron star (NS) meets the stellar surface, it should brake down to match the rotation of the NS, forming a boundary layer. As the mechanisms potentially responsible for this braking are apparently inefficient, it is reasonable to consider this layer as a spreading layer (SL) with negligible radial extent and structure. We perform hydrodynamical 2D spectral simulations of an SL, considering the disc as a source of matter and angular momentum. Interaction of new, rapidly rotating matter with the pre-existing, relatively slow material co-rotating with the star leads to instabilities capable of transferring angular momentum and creating variability on dynamical timescales. For small accretion rates, we find that the SL is unstable for heating instability that disrupts the initial latitudinal symmetry and produces large deviations between the two hemispheres. This instability also results in breaking of the axial symmetry as coherent flow structures are formed and escape from the SL intermittently. At enhanced accretion rates, the SL is prone to shearing instability and acts as a source of oblique waves that propagate towards the poles, leading to patterns that again break the axial symmetry. We compute artificial light curves of an SL viewed at different inclination angles. Most of the simulated light curves show oscillations at frequencies close to 1 kHz. We interpret these oscillations as inertial modes excited by shear instabilities near the boundary of the SL. Their frequencies, dependence on flux, and amplitude variations can explain the high-frequency pair quasi-periodic oscillations observed in many low-mass X-ray binaries.

We computed accurate atmosphere models of rotation-powered millisecond pulsars in which the polar caps of a neutron star (NS) are externally heated by magnetospheric return currents. The external ram pressure, energy losses, and stopping depth of the penetrating charged particles were computed self-consistently with the atmosphere model, instead of assuming a simplified deep-heated atmosphere in radiative equilibrium. We used exact Compton scattering formalism to model the properties of the emergent X-ray radiation. The deep-heating approximation was found to be valid only if most of the heat originates from ultra-relativistic bombarding particles with Lorentz factors of γ ≳ 100. In the opposite regime, the atmosphere attains a distinct two-layer structure with an overheated optically thin skin on top of an optically thick cool plasma. The overheated skin strongly modifies the emergent radiation: It produces a Compton-upscattered high-energy tail in the spectrum and alters the radiation beaming pattern from limb darkening to limb brightening for emitted hard X-rays. This kind of drastic change in the emission properties can have a significant impact on the inferred NS pulse profile parameters as performed, for example, by Neutron star Interior Composition ExploreR. Finally, the connection between the energy distribution of the return current particles and the atmosphere emission properties offers a new tool to probe the exact physics of pulsar magnetospheres.

In this White Paper we present the potential of the Enhanced X-ray Timing and Polarimetry (eXTP) mission for determining the nature of dense matter; neutron star cores host an extreme density regime which cannot be replicated in a terrestrial laboratory. The tightest statistical constraints on the dense matter equation of state will come from pulse profile modelling of accretion-powered pulsars, burst oscillation sources, and rotation-powered pulsars. Additional constraints will derive from spin measurements, burst spectra, and properties of the accretion flows in the vicinity of the neutron star. Under development by an international Consortium led by the Institute of High Energy Physics of the Chinese Academy of Sciences, the eXTP mission is expected to be launched in the mid 2020s.

We use ensemble machine learning algorithms to study the evolution of magnetic fields in magnetohydrodynamic (MHD) turbulence that is helically forced. We perform direct numerical simulations of helically forced turbulence using mean field formalism, with electromotive force (EMF) modeled both as a linear and non-linear function of the mean magnetic field and current density. The form of the EMF is determined using regularized linear regression and random forests. We also compare various analytical models to the data using Bayesian inference with Markov chain Monte Carlo (MCMC) sampling. Our results demonstrate that linear regression is largely successful at predicting the EMF and the use of more sophisticated algorithms (random forests, MCMC) do not lead to significant improvement in the fits. We conclude that the data we are looking at is effectively low dimensional and essentially linear. Finally, to encourage further exploration by the community, we provide all of our simulation data and analysis scripts as open source IPYTHON notebooks.

Transitional millisecond pulsars provide a unique set of observational data for understanding accretion at low rates onto magnetized neutron stars. In particular, PSR.J1023+0038 exhibits a remarkable bimodality of the X-ray luminosity (low and high modes), pulsations extending from the X-ray to the optical band, GeV emission, and occasional X-ray flares. We discuss a scenario for the pulsar interaction with the accretion disk capable of explaining the observed behavior. We suggest that during the high mode the disk is truncated outside the light cylinder, allowing the pulsar wind to develop near the equatorial plane and strike the disk. The dissipative wind-disk collision energizes the disk particles and generates synchrotron emission, which peaks in the X-ray band and extends down to the optical band. The emission is modulated by the pulsar wind rotation, resulting in a pulse profile with two peaks 180 degrees apart. This picture explains the high mode luminosity, spectrum, and pulse profile (X-ray and optical) of PSR.J1023+0038. It may also explain the X-ray flares as events of sudden increase in the effective disk cross section intercepting the wind. In contrast to previously proposed models, we suggest that the disk penetrates the light cylinder only during the low X-ray mode. This penetration suppresses the dissipation caused by the pulsar wind-disk collision, and the system enters the propeller regime. The small duty cycle of the propeller explains the low spin-down rate of the pulsar.

We present a Bayesian method to constrain the masses and radii of neutron stars (NSs) using the information encoded in the X-ray pulse profiles of accreting millisecond pulsars. We model the shape of the pulses using oblate Schwarzschild approximation, which takes into account the deformed shape of the star together with the special and general relativistic corrections to the photon trajectories and angles. The spectrum of the radiation is obtained from an empirical model of Comptonization in a hot slab in which a fraction of seed black-body photons is scattered into a power-law component. By using an affine-invariant Markov chain Monte Carlo ensemble sampling method, we obtain posterior probability distributions for the different model parameters, especially for the mass and the radius. To test the robustness of our method, we first analysed self-generated synthetic data with known model parameters Similar analysis was then applied for the observations of SAX J1808.4-3658 by the Rossi X-ray Timing Explorer (RXTE). The results show that our method can reproduce the model parameters of the synthetic data, and that accurate constraints for the radius can be obtained using the RXTE pulse profile observations if the mass is a priori known. For a mass in the range 1.5-1.8 M-circle dot, the radius of the NS in SAX J1808.4-3658 is constrained between 9 and 13 km. If the mass is accurately known, the radius can be determined with an accuracy of 5% (68% credibility). For example, for the mass of 1.7 M-circle dot the equatorial radius is R-eq = 11.9(-0.4)(+0.5) km. Finally, we show that further improvements can be obtained when the X-ray polarization data from the Imaging X-ray Polarimeter Explorer will become available.

Thermonuclear X-ray bursts on the surface of neutron stars (NSs) can enrich the photosphere with metals, which may imprint photoionization edges on the burst spectra. We report here the discovery of absorption edges in the spectra of the type I X-ray burst from the NS low-mass X-ray binary GRS 1747-312 in Terzan 6 during observations by the Rossi X-ray Timing Explorer. We find that the edge energy evolves from 9.45 +/- 0.51 to similar to 6 keV and then back to 9.44 +/- 0.40 keV during the photospheric radius expansion phase and remains at 8.06 +/- 0.66 keV in the cooling tail. The photoionization absorption edges of hydrogen-like Ni, Fe, or an Fe/Ni mixture and the bound-bound transitions of metals may be responsible for the observed spectral features. The ratio of the measured absorption edge energy in the cooling tail to the laboratory value of the hydrogen-like Ni(Fe) edge energy allows us to estimate the gravitational redshift factor 1 + z = 1.34 +/- 0.11(1 + z = 1.15 +/- 0.09). The evolution of the spectral parameters during the cooling tail are well described by metal-rich atmosphere models. The combined constraints on the NS mass and radius from the direct cooling method and the tidal deformability strongly suggest very high atmospheric abundance of the iron group elements and limit the distance to the source to 11 +/- 1 kpc.