Code PLUTO Research Articles

Simulating the Escaping Atmosphere of GJ 436 b with Two-fluid Magnetohydrodynamic Models

Observations of transmission spectra reveal that hot Jupiters and Neptunes are likely to possess escaping atmospheres driven by stellar radiation. Numerous models predict that magnetic fields may exert significant influences on the atmospheres of hot planets. Generally, the escaping atmospheres are not entirely ionized, and magnetic fields only directly affect the escape of ionized components within them. Considering the chemical reactions between ionized components and neutral atoms, as well as collision processes, magnetic fields indirectly impact the escape of neutral atoms, thereby influencing the detection signals of planetary atmospheres in transmission spectra. In order to simulate this process, we developed a magnetohydrodynamic multi-fluid model based on MHD code PLUTO. As an initial exploration, we investigated the impact of magnetic fields on the decoupling of H+ and H in the escaping atmosphere of the hot Neptune GJ436b. Due to the strong resonant interactions between H and H+, the coupling between them is tight even if the magnetic field is strong. Of course, alternatively, our work also suggests that merging H and H+ into a single flow can be a reasonable assumption in MHD simulations of escaping atmospheres. However, our simulation results indicate that under the influence of magnetic fields, there are noticeable regional differences in the decoupling of H+ and H. With the increase of magnetic field strength, the degree of decoupling also increases. For heavier particles such as O, the decoupling between O and H+ is more pronounced. Our findings provide important insights for future studies on the decoupling processes of heavy atoms in the escaping atmospheres of hot Jupiters and hot Neptunes under the influence of magnetic fields.

FR0 jets and recollimation-induced instabilities

Context. The recently discovered population of faint Fanaroff-Riley type 0 (FR0) radio galaxies has been interpreted as the extension to low power of the classical FRI sources. Their radio emission appears to be concentrated in very compact parsec scale cores, any extended emission is very weak or absent, and very-long-baseline interferometry (VLBI) observations show that jets are already mildly or sub-relativistic at parsec scales. Based on these observational properties, we propose here that the jets of FR0s are strongly decelerated and disturbed at the parsec scale by hydrodynamical instabilities. Aims. With the above scenario in mind, we studied the dynamics of a low-power relativistic jet propagating into a confining external medium, focusing on the effects of entrainment and mixing promoted by the instabilities developing at the jet-environment interface downstream of a recollimation shock. Methods. We performed a 3D relativistic hydrodynamical simulation of a recollimated jet by means of the state-of-the-art code PLUTO. The jet was initially conical, relativistic (with an initial Lorentz factor Γ = 5), cold, and light with respect to the confining medium, whose pressure is assumed to slowly decline with distance. The magnetic field is assumed to be dynamically unimportant. Results. The 3D simulation shows that, after the first recollimation and reflection shock system, a rapidly growing instability develops, as a result of the interplay between recollimation-induced instabilities and Richtmyer-Meshkov modes. In turn, the instabilities promote strong mixing and entrainment that rapidly lead to the deceleration of the jet and spread its momentum to slowly moving, highly turbulent external gas. We argue that this mechanism could account for the peculiarities of the low-power FR0 jets. For outflows with a higher power, Lorentz factor, or magnetic field, we expect that the destabilizing effects are less effective, allowing the survival of the jet up to the kiloparsec scale, as observed in FRIs.

Magnetically threaded accretion disks in resistive magnetohydrodynamic simulations and asymptotic expansion

Aims. A realistic model of magnetic linkage between a central object and its accretion disk is a prerequisite for understanding the spin history of stars and stellar remnants. To this end, we aim to provide an analytic model in agreement with magnetohydrodynamic (MHD) simulations. Methods. For the first time, we wrote a full set of stationary asymptotic expansion equations of a thin magnetic accretion disk, including the induction and energy equations. We also performed a resistive MHD simulation of an accretion disk around a star endowed with a magnetic dipole, using the publicly available code PLUTO. We compared the analytical results with the numerical solutions, and discussed the results in the context of previous solutions of the induction equation describing the star-disk magnetospheric interaction. Results. We found that the magnetic field threading the disk is suppressed by orders of magnitude inside thin disks, so the presence of the stellar magnetic field does not strongly affect the velocity field, nor the density profile inside the disk. Density and velocity fields found in the MHD simulations match the radial and vertical profiles of the analytic solution. Qualitatively, the MHD simulations result in an internal magnetic field similar to the solutions previously obtained by solving the induction equation in the disk alone. However, the magnetic field configuration is quantitatively affected by magnetic field inflation outside the disk; this is reflected in the net torque. The torque on the star is an order of magnitude larger in the magnetic than in the non-magnetic case. Spin-up of the star occurs on a timescale comparable to the accretion timescale in the MHD case, and is an order of magnitude slower in the absence of a stellar magnetic field.

On Earth’s habitability over the Sun’s main-sequence history: joint influence of space weather and Earth’s magnetic field evolution

ABSTRACT The aim of this study is to analyse the Earth habitability with respect to the direct exposition of the Earth atmosphere to the solar wind (SW) along the Sun’s evolution on the main sequence including the realistic evolution of the space weather conditions and the Earth magnetic field. The MHD code PLUTO in spherical coordinates is applied to perform parametric studies with respect to the SW dynamic pressure and the interplanetary magnetic field intensity for different Earth magnetic field configurations. Quiet space weather conditions may not impact the Earth habitability. On the other hand, the impact of interplanetary coronal mass ejections (ICME) could lead to the erosion of the primary Earth atmosphere during the Hadean eon. A dipolar field of 30 μT is strong enough to shield the Earth from the Eo-Archean age as well as 15 and 5 μT dipolar fields from the Meso-Archean and Meso-Proterozoic, respectively. Multipolar weak field period during the Meso-Proterozoic age may not be a threat for ICME-like space weather conditions if the field intensity is at least 15 μT and the ratio between the quadrupolar (Q) and dipolar (D) coefficients is $\frac{Q}{D} \le 0.5$. By contrast, the Earth habitability in the Phanerozoic eon (including the present time) can be hampered during multipolar low field periods with a strength of 5 μT and $\frac{Q}{D} \ge 0.5$ associated with geomagnetic reversals. Consequently, the effect of the SW should be considered as a possible driver of Earth’s habitability.

The Mass Fractionation of Helium in the Escaping Atmosphere of HD 209458b* *Revised Manuscript on 2022 June 7

The absorption signals of metastable He in HD 209458b and several other exoplanets can be explained via an escaping atmosphere model with a subsolar He/H ratio. The low abundance of helium can be a result of planet formation if there is a small amount of helium in their primordial atmosphere. However, another possibility is that the low He/H ratio is caused by the process of mass fractionation of helium in the atmosphere. In order to investigate the effect of fractionation in the hydrogen-helium atmosphere, we developed a self-consistent multi-fluid 1D hydrodynamic model based on the well-known open-source MHD code PLUTO. Our simulations show that a lower He/H ratio can be produced spontaneously in the multi-fluid model. We further modeled the transmission spectra of He 10830 lines for HD 209458b in a broad parameter space. The transmission spectrum of the observation can be fitted in the condition of 1.80 times the X-ray and extreme-ultraviolet flux of the quiet Sun. Meanwhile, the ratio of the escaping flux of helium to hydrogen, F He/F H, is 0.039. Our results indicate that the mass fractionation of helium to hydrogen can naturally interpret the low He/H ratio required by the observation. Thus, in the escaping atmosphere of HD 209458b, decreasing the abundance of helium in the atmosphere is not needed even if its He abundance is similar to that of the Sun. The simulation presented in this work hints that in the escaping atmosphere, mass fractionation can also occur on other exoplanets, which needs to be explored further.

HelioCast: heliospheric forecasting based on white-light observations of the solar corona

We present a new 3D magnetohydrodynamic (MHD) heliospheric model for space-weather forecasting driven by boundary conditions defined from white-light observations of the solar corona. The model is based on the MHD code PLUTO, constrained by an empirical derivation of the solar wind background properties at 0.1 au. This empirical method uses white-light observations to estimate the position of the heliospheric current sheet (HCS). The boundary conditions necessary to run HelioCast are then defined from pre-defined relations between the necessary MHD properties (speed, density, and temperature) and the distance to the current sheet. We assess the accuracy of the model over six Carrington rotations during the first semester of 2018. Using point-by-point metrics and event-based analyses, we evaluate the performances of our model varying the angular width of the slow solar wind layer surrounding the HCS. We also compare our empirical technique with two well-tested models of the corona: Multi-VP and WindPredict-AW. We find that our method is well-suited to reproduce high-speed streams, and does – for well-chosen parameters – better than full MHD models. The model shows, nonetheless, limitations that could worsen for rising and maximum solar activity.

MHD Study of Extreme Space Weather Conditions for Exoplanets With Earth‐Like Magnetospheres: On Habitability Conditions and Radio‐Emission

AbstractThe present study aims at characterizing the habitability conditions of exoplanets with an Earth‐like magnetosphere inside the habitable zone of M and F stars, caused by the direct deposition of the stellar wind on the exoplanet surface. Also, the radio emission generated by exoplanets with a Earth‐like magnetosphere is calculated for different space weather conditions. The study is based on a set of MHD simulations performed by the code PLUTO. Exoplanets hosted by M stars at 0.2 au are protected from the stellar wind during regular and coronal mass ejection (CME)‐like space weather conditions if the star rotation period is slower than 3 days. Exoplanets hosted by a F stars at ≥2.5 au are protected during regular space weather conditions, but a stronger magnetic field compared to the Earth is mandatory if the exoplanet is close to the inner edge of the star habitable zone (2.5 au) during CMEs. The range of radio emission values calculated in the simulations are consistent with the scaling proposed by Zarka (2018,https://doi.org/10.1007/978-3-319-55333-7_22) during regular and common CME‐like space weather conditions. If the radio telescopes measure a relative low radio emission signal with small variability from an exoplanet, that may indicate favorable exoplanet habitability conditions. The radio emission power calculated for exoplanets with an Earth‐like magnetosphere is in the range of 3 × 107to 2 × 1010 W for SW dynamic pressures between 1.5 and 100 nPa and interplanetary magnetic field intensities between 50 and 250 nT, and is below the sensitivity threshold of present radio telescopes at parsec distances.

MRI turbulence in accretion discs at large magnetic Prandtl numbers

ABSTRACT The effect of large magnetic Prandtl number Pm (the ratio of viscosity to resistivity) on the turbulent transport and energetics of the magnetorotational instability (MRI) is poorly understood, despite the realization of this regime in astrophysical environments as disparate as discs from binary neutron star (BNS) mergers, the inner regions of low-mass X-ray binaries and active galactic nuclei, and the interiors of protoneutron stars. We investigate the MRI dynamo and associated turbulence in the regime Pm > 1 by carrying out fully compressible, 3D MHD-shearing box simulations using the finite-volume code pluto, focusing mostly on the case of Keplerian shear relevant to accretion discs. We find that when the magnetic Reynolds number is kept fixed, the turbulent transport (as parameterized by α, the ratio of stress to thermal pressure) scales with the magnetic Prandtl number as α ∼ Pmδ, with δ ∼ 0.5−0.7 up to Pm ∼ 128. However, this scaling weakens as the magnetic Reynolds number is increased. Importantly, compared to previous studies, we find a new effect at very large Pm – the turbulent energy and stress begin to plateau, no longer depending on Pm. To understand these results we have carried out a detailed analysis of the turbulent dynamics in Fourier space, focusing on the effect of increasing Pm on the transverse cascade – a key non-linear process induced by the disc shear flow that is responsible for the sustenance of MRI turbulence. Finally, we find that α–Pm scaling is sensitive to the box vertical-to-radial aspect ratio, as well as to the background shear.

The stress–pressure lag in MRI turbulence and its implications for thermal instability in accretion discs

ABSTRACT The classical alpha-disc model assumes that the turbulent stress scales linearly with – and responds instantaneously to – the pressure. It is likely, however, that the stress possesses a non-negligible relaxation time and will lag behind the pressure on some time-scale. To measure the size of this lag we carry out unstratified 3D magnetohydrodynamic shearing box simulations with zero-net-magnetic-flux using the finite-volume code pluto. We impose thermal oscillations of varying periods via a cooling term, which in turn drives oscillations in the turbulent stress. Our simulations reveal that the stress oscillations lag behind the pressure by ∼5 orbits in cases where the oscillation period is several tens of orbits or more. We discuss the implication of our results for thermal and viscous overstability in discs around compact objects.

Numerical Modeling and Physical Interplay of Stochastic Turbulent Acceleration for Nonthermal Emission Processes

Abstract Particle acceleration is a ubiquitous phenomenon in astrophysical and space plasma. Diffusive shock acceleration (DSA) and stochastic turbulent acceleration (STA) are known to be the possible mechanisms for producing very highly energetic particles, particularly in weakly magnetized regions. An interplay of different acceleration processes along with various radiation losses is typically observed in astrophysical sources. While DSA is a systematic acceleration process that energizes particles in the vicinity of shocks, STA is a random energizing process, where the interaction between cosmic ray particles and electromagnetic fluctuations results in particle acceleration. This process is usually interpreted as a biased random walk in energy space, modeled through a Fokker–Planck equation. In the present work, we describe a novel Eulerian algorithm, adopted to incorporate turbulent acceleration in the presence of DSA and radiative processes like synchrotron and inverse Compton emission. The developed framework extends the hybrid Eulerian−Lagrangian module in a full-fledged relativistic Magneto-hydrodynamic (RMHD) code PLUTO. From our validation tests and case studies, we showcase the competing and complementary nature of both acceleration processes. Axisymmetric simulations of an RMHD jet with this extended hybrid framework clearly demonstrate that emission due to shocks is localized, while that due to turbulent acceleration originates in the backflow and is more diffuse, particularly in the high-energy X-ray band.

Systematic description of wind-driven protoplanetary discs

Aims. Planet-forming discs are believed to be very weakly turbulent in the regions outside of 1 AU. For this reason, it is now believed that magnetised winds could be the dominant mechanism driving accretion in these systems. However, currently, no self-consistent approach can describe discs that are subject to a magnetised wind in a way similar to the α disc model. In this article, I explore the parameter space of wind-driven protoplanetary discs in a systematic manner and present scaling laws that can be used in reduced models in a similar way to α disc models. Methods. I computed a series of self-similar wind solutions, assuming that the disc is dominated by ambipolar and Ohmic diffusion. These solution were obtained by searching for stationary solutions in the finite-volume code PLUTO using a relaxation method and continuation. Results. Self-similar solutions are obtained for values of plasma β ranging from 102 to 108 for several Ohmic and ambipolar diffusion strengths. Mass accretion rates of about 10−8 M⊙ yr−1 are obtained for the poloidal field strength β = O(104) or equivalently, 1 mG at 10 AU. In addition, the ejection efficiency is always close to 1, implying that wind mass-loss rate can be higher than the inner mass accretion rate when the wind-emitting region is large. The resulting magnetic lever arms are typically lower than 2, possibly reaching 1.5 in the weakest field cases. Remarkably, the mean transport properties (accretion rate and mass-loss rate) mostly depend on the field strength and much less on the disc diffusivities or surface density. The disc internal structure is nevertheless strongly affected by Ohmic resistivity, strongly resistive discs being subject to accretion at the surface while ambipolar only models lead to mid-plane accretion. Finally, I provide a complete set of scaling laws and semi-analytical wind solutions, which can be used to fit and interpret observations. Conclusions. Magnetised winds are unavoidable in protoplanetary discs as soon as they are embedded in an ambient poloidal magnetic field. Very detailed disc microphysics are not always needed to describe them, and simplified models such as self-similar solutions can capture most of the physics seen in full 3D simulations. The remaining difficulty to set up a complete theory of wind-driven accretion lies in the transport of the large-scale field, which remains poorly constrained and is not well understood.

Magnetohydrodynamic convection in accretion discs

ABSTRACT Convection has been discussed in the field of accretion discs for several decades, both as a means of angular momentum transport and also because of its role in controlling discs’ vertical structure via heat transport. If the gas is sufficiently ionized and threaded by a weak magnetic field, convection might interact in non-trivial ways with the magnetorotational instability (MRI). Recently, vertically stratified local simulations of the MRI have reported considerable variation in the angular momentum transport, as measured by the stress to thermal pressure ratio α, when convection is thought to be present. Although MRI turbulence can act as a heat source for convection, it is not clear how the two instabilities will interact dynamically. Here, we investigate their interplay in controlled numerical experiments and isolate the generic features of their interactions. We perform vertically stratified, 3D magnetohydrodynamic shearing box simulations with a perfect gas equation of state with the conservative, finite-volume code pluto. We find two characteristic outcomes of the interaction between the two instabilities: (a) straight MRI and (b) MRI/convective cycles, with the latter exhibiting alternating phases of convection-dominated and MRI-dominated flow. During the latter phase, we find that α is enhanced by nearly an order of magnitude, reaching peak values of ∼0.08. In addition, we find that convection in the non-linear phase takes the form of large-scale and oscillatory convective cells. Convection can also help the MRI persist to lower Rm than it would otherwise do. Finally, we discuss how our results help interpret simulations of dwarf novae.

A new ionization network and radiation transport module in pluto

ABSTRACT We introduce a new general-purpose time-dependent ionization network (IN) and a radiation transport (RT) module for the magnetohydrodynamic (MHD) code pluto. Our IN is reliable for temperatures ranging from 5 × 103 to 3 × 108 K and includes all ionization states of H, He, C, N, O, Ne, Mg, Si, S, and Fe, making it suitable for studying a variety of astrophysical scenarios. Radiation loss for each ion–electron pair is calculated using cloudy-17 data on the fly. Photoionization and charge exchange are the main processes contributing to chemical heating. The IN is fully coupled to the RT module over a large range of opacities at different frequencies. The RT module employs a method of short characteristics assuming spherical symmetry. The radiation module requires the assumption of spherical symmetry, while the IN is compatible with full 3D. We also include a simple prescription for dust opacity, grain destruction, and the dust contribution to radiation pressure. We present numerical tests to show the reliability and limitations of the new modules. We also present a post-processing tool to calculate projected column densities and emission spectra.

Modelling of particle acceleration in the pulsar wind nebulae with bow shocks

The observations of the pulsar wind nebulae (PWNe) indicate the presence of efficient acceleration of positrons and electrons in these sources. The Fermi acceleration in the colliding shock flows can explain the observed hard synchrotron emission spectra of PWNe with bow shocks (BSPWNe). This may result in their maximal luminosities in the far ultraviolet range (FUV; 1250 — 2000 Å, ∼ 6 — 10 eV) due to the synchrotron emission of pairs rather than due to the thermal emission of the shocked interstellar matter. Fine spectroscopic observations of sufficiently bright sources with Hubble Space Telescope could be applied to distinguish between these two scenarios. In this paper we simulate BSPWNe flows structure with the relativistic magnetohydrodynamic code PLUTO, consider particle transport and their synchrotron emission for a number of BSPWNe. We calculate the synchrotron FUV luminosities of these BSPWNe and discuss the prospective of their observation in FUV. We also consider possible contribution of PSR J1741-2054 to the positron excess detected by AMS-02 and PAMELA.

How Alfv\xe9n waves induce compressive flows in the neighborhood of a 2.5D magnetic null-point

The aim of the present study is to provide insight on the induced compressive perturbations together with the modifications of the environmental parameters in the course of Alfvén wave interaction with a solar magnetic null-point. The shock-capturing Godunov-type code PLUTO is used to solve the set of ideal magnetohydrodynamic equations. The nonlinear effects connected with an initial Alfvén pulse nearing a magnetic null point induces fast and slow magnetoacoustic waves with anti phase conduct. The induced current density and flows are independent of the local plasma-beta at the reconnection site. The induced inflows and outflows highly depend on the polarization. The inflows have a stronger effect compared to the outflows in both the x and y directions showing its peak in the x-direction. The dominant wave that couples to flows is the fast wave due to the in-phase harmony between perturbations of the compressive parameters and the fast wave. The induced current density possesses a steady orientation at the reconnection site which governs the diffusion or propagation of the waves. Induced perturbations by the nonlinear forces together with their back reaction on the Alfvén wave have a significant role in the current density excitation being responsible for the creation of inflows and outflows that are possible candidates for the creation of solar jets which has a significant contribution towards coronal seismology.

Three-Dimensional Simulations of Solar Wind Preconditioning and the 23 July 2012 Interplanetary Coronal Mass Ejection

Predicting the large-scale eruptions from the solar corona and their propagation through interplanetary space remains an outstanding challenge in solar- and helio-physics research. In this article, we describe three-dimensional magnetohydrodynamic simulations of the inner heliosphere leading up to and including the extreme interplanetary coronal mass ejection (ICME) of 23 July 2012, developed using the code PLUTO. The simulations are driven using the output of coronal models for Carrington rotations 2125 and 2126 and, given the uncertainties in the initial conditions, are able to reproduce an event of comparable magnitude to the 23 July ICME, with similar velocity and density profiles at 1 au. The launch time of this event is then varied with regards to an initial 19 July ICME and the effects of solar wind preconditioning are found to be significant for an event of this magnitude and to decrease over a time-window consistent with the ballistic refilling of the depleted heliospheric sector. These results indicate that the 23 July ICME was mostly unaffected by events prior, but would have traveled even faster had it erupted closer in time to the 19 July event where it would have experienced even lower drag forces. We discuss this systematic study of solar wind preconditioning in the context of space weather forecasting.

Phase-modulated X-Ray Emission from Cepheids due to Pulsation-driven Shocks

Abstract Cepheids are pulsating variable stars with a periodic chromospheric response at UV wavelengths close to their minimum radius phase. Recently, an X-ray variable signature was captured in observations during the maximum radius phase. This X-ray emission came as a surprise and is not understood. In this work, we use the modern astrophysical code PLUTO to investigate the effects of pulsations on Cepheid X-ray emission. We run a number of hydrodynamic numerical simulations with a variety of initial and boundary conditions in order to explore the capability of shocks to produce the observed phase-dependent X-ray behavior. Finally, we use the Simulated Observations of X-ray Sources (SOXS) package to create synthetic spectra for each simulation case and link our simulations to observables. We show that, for certain conditions, we can reproduce observed X-ray fluxes at phases 0.4–0.8 when the Cepheid is at maximum radius. Our results span a wide range of mass-loss rates, 2 × 10−13 M ⊙ yr−1 to 3 × 10−8 M ⊙ yr−1, and peak X-ray luminosities, 5 × 10−17 erg cm−2 s−1 to 1.4 × 10−12 erg cm−2 s−1. We conclude that Cepheids exhibit two-component emission with (a) shock waves being responsible for the phase-dependent variable emission (phases 0.2–0.6) and (b) a separate quiescent mechanism being the dominant emission mechanism for the remaining phases.

Magnetic field configurations in neutron stars from MHD simulations

ABSTRACT We have studied numerically the evolution of magnetic fields in barotropic neutron stars, by performing non-linear magnetohydrodynamical simulations with the code pluto. For both initially predominantly poloidal and toroidal fields, with varying strengths, we find that the field settles down to a mixed poloidal–toroidal configuration, where the toroidal component contributes between ${\rm 10}$ and $20 {{\ \rm per\ cent}}$ of the total magnetic energy. This is, however, not a strict equilibrium, as the instability leads to the development of turbulence, which, in turn, gives rise to an inverse helicity cascade, which determines the final ‘twisted torus’ setup. The final field configuration is thus dictated by the non-linear saturation of the instability, and is not stationary. The average energy of the poloidal and toroidal components, however, is approximately stable in our simulations, and a complex multipolar structure emerges at the surface, while the magnetic field is dipolar at the exterior boundary, outside the star.

Interaction of a supernova remnant with a wind of young massive star: MHD simulations

Fast stellar winds of young massive stars and supernovae play an essential role in the evolution of the interstellar medium in the Milky Way. They dominate the dynamics of the active starforming regions in starburst galaxies and in particularly their non-thermal radiation. The non-thermal emission from starburst galaxies NGC 253, M82, NGC 1068 was detected from radio to gamma-rays. To model the observed radiation, one needs to know the complex structure of the flows and magnetic fields in the starforming regions. In this work a 2.5D magnetohydrodynamic (MHD) simulation (i.e., we consider 3D vector fields, but assume a 2D axisymmetric geometry) of the interaction between a stellar wind and a supernova remnant shock wave is carried out using the MHD module of the code PLUTO. The structure of the flows in the collision region is obtained taking into account the magnetic fields of the rotating source stars. We present the profiles of the bulk plasma velocity, density and magnetic field in the collision region. The amplification of regular magnetic fields in such systems is studied.

A Radiative Transfer Module for Relativistic Magnetohydrodynamics in the PLUTO Code

Abstract We present a numerical implementation for the solution of the relativistic radiation hydrodynamics and magnetohydrodynamics equations, designed as an independent module within the freely available code PLUTO. The radiation transfer equations are solved under the gray approximation and imposing the M1 closure, which allows the radiation transport to be handled in both the free-streaming and diffusion limits. Equations are integrated following an implicit–explicit scheme, where radiation–matter 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 introduce 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. The algorithm performance is demonstrated through a series of numerical benchmarks by investigating various different configurations with a particular emphasis on the behavior of the solutions in the free-streaming and diffusion limits.