Beyond moments: relativistic Lattice-Boltzmann methods for radiative transport in computational astrophysics
L. R. Weih, A. Gabbana, D. Simeoni, L. Rezzolla, S. Succi, R. Tripiccione
MMNRAS , 000–000 (2020) Preprint 6 October 2020 Compiled using MNRAS L A TEX style file v3.0
Beyond moments: relativistic Lattice-Boltzmann methods forradiative transport in computational astrophysics
L. R. Weih ★ , A. Gabbana , D. Simeoni , , , L. Rezzolla , , , S. Succi , , , R. Tripiccione Institut für Theoretische Physik, Max-von-Laue-Str. 1, 60438 Frankfurt, Germany Università di Ferrara and INFN-Ferrara, I-44122 Ferrara, Italy Bergische Universität Wuppertal, D-42119 Wuppertal, Germany University of Cyprus, CY-1678 Nicosia, Cyprus School of Mathematics, Trinity College, Dublin 2, Ireland Helmholtz Research Academy Hesse for FAIR, Max-von-Laue-Str. 12, 60438 Frankfurt, Germany Center for Life Nano Science @ La Sapienza, Italian Institute of Technology, Viale Regina Elena 295, I-00161 Roma, Italy Istituto Applicazioni del Calcolo, National Research Council of Italy, Via dei Taurini 19, I-00185 Roma, Italy Harvard-SEAS, Oxford Street 29, 02130 Cambridge, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We present a new method for the numerical solution of the radiative-transfer equation (RTE)in multidimensional scenarios commonly encountered in computational astrophysics. Themethod is based on the direct solution of the Boltzmann equation via an extension of theLattice Boltzmann (LB) equation and allows to model the evolution of the radiation fieldas it interacts with a background fluid, via absorption, emission, and scattering. As a firstapplication of this method, we restrict our attention to a frequency independent (“grey”)formulation within a special-relativistic framework, which can be employed also for classicalcomputational astrophysics. For a number of standard tests that consider the performance of themethod in optically thin, optically thick and intermediate regimes with a static fluid, we showthe ability of the LB method to produce accurate and convergent results matching the analyticsolutions. We also contrast the LB method with commonly employed moment-based schemesfor the solution of the RTE, such as the M1 scheme. In this way, we are able to highlightthat the LB method provides the correct solution for both non-trivial free-streaming scenariosand the intermediate optical-depth regime, for which the M1 method either fails or providesinaccurate solutions. When coupling to a dynamical fluid, on the other hand, we present thefirst self-consistent solution of the RTE with LB methods within a relativistic-hydrodynamicscenario. Finally, we show that besides providing more accurate results in all regimes, the LBmethod features smaller or comparable computational costs compared to the M1 scheme.
Key words: radiative transfer – radiation: dynamics – neutrinos – scattering – methods:numerical
The proper treatment of the dynamics of radiation, as it interactswith a matter fluid, is a fundamental problem in essentially allastrophysical phenomena and requires the solution of the radiative-transport equation (RTE). Given the complexity of the RTE andthe nonlinear regimes normally encountered in astrophysical andcosmological scenarios, the inclusion of radiative effects inevitablycommands the use of advanced numerical methods to solve theRTE (see, e.g., Paardekooper et al. 2010, for the solution on anunstructured dynamic grid).An important and representative example is a binary system of ★ [email protected] merging neutron stars (see Baiotti & Rezzolla (2017); Paschalidis(2017) for an overview), where the radiative transport of neutri-nos can alter significantly the chemical composition of the ejectedmatter, the efficiency of its ejection, as well as the stability of thepost-merger object. Furthermore, radiative-transport effects are ex-pected to play a fundamental role in the kilonova signal that isproduced (Rosswog et al. 2014; Dietrich & Ujevic 2017; Siegel &Ciolfi 2016; Bovard et al. 2017; Perego et al. 2017; Fujibayashiet al. 2018; Siegel & Metzger 2017; Fernández et al. 2019) andmight even be relevant for producing a short gamma-ray burst as-sociated with such a merger. An equally important astrophysicalscenario where radiation transport plays a fundamental role is theone explored in simulations of neutrino-driven core-collapse super-novae (Mezzacappa et al. 2001; O’Connor 2015; Just et al. 2015; © a r X i v : . [ phy s i c s . c o m p - ph ] O c t Weih et al.
Kuroda et al. 2016), for which Colgate & White (1966) first showedthat the radiation in form of neutrinos is essential for the explosionmechanism. This is because the explosion mechanism may requirea fine balance, of the order of a few percent, between the energydeposited in the stalled accretion shock and the release of potentialgravitational energy by the collapsing matter (see also Janka et al.(2007) for an overview). Finally, one more classical astrophysicalscenario where radiative-transfer effects cannot be ignored is thestudy of accretion flows around black holes (Zanotti et al. 2011;Fragile et al. 2012; Roedig et al. 2012; Sa¸dowski et al. 2013; McK-inney et al. 2014), for which a broad array of techniques has beendeveloped over the years to compare with the observations (EventHorizon Telescope Collaboration et al. 2019).Unfortunately, the computational cost associated with the so-lution of the RTE in numerical astrophysics also represents a sig-nificant obstacle to the inclusion of radiative effects in numericalsimulations. This is due to the properties of the fundamental equa-tion behind the RTE, i.e., the Boltzmann equation for masslessparticles, which lives in a seven-dimensional space of time (onedimension), configuration space (three dimensions) and momen-tum space (three more dimensions). As a result, the solution ofthe RTE for typical astrophysical scenarios as the ones mentionedabove exceeds the current capacities of supercomputers and thusapproximate methods have to be used .A low-order approximation first used in the context of super-nova explosions by Livio et al. (1980) and commonly employedin binary neutron-star simulations is the leakage scheme (Ruffertet al. 1996; Rosswog & Liebendörfer 2003; Galeazzi et al. 2013;Perego et al. 2014; Most et al. 2019), which only allows for cool-ing via the emission of neutrinos and is therefore not useful forcore-collapse supernovae simulations, where heating is essential forreviving the shock. These effects can be included by an approxima-tion of similar simplicity, the flux-limited diffusion approximationfirst implemented for astrophysical simulations of core-collapse su-pernovae by Bruenn (1975) (see also Pomraning 1981; Levermore& Pomraning 1981); a recent implementation of this scheme hasbeen presented by Rahman et al. (2019). In this method, the zerothmoment of the radiation distribution function, i.e., the radiationenergy-density, is evolved together with the fluid quantities. Sincethe zeroth moment does not provide any information about the di-rection of the radiation fluxes, its implementation is useful only forsystems with clear underlying symmetries. These symmetries arenot present in the case of binary neutron-star simulations, so theflux-limited diffusion does not offer but a crude approximation ofthe radiative effects. This is also true for the M0 scheme developedby Radice et al. (2016), which also evolves the lowest moment ofthe distribution function, but in the free-streaming limit.A considerably better approximation is the so-called truncatedmoment-based scheme developed by Thorne (1981) and first im-plemented in general relativity by Rezzolla & Miller (1994) in onedimension and by Shibata et al. (2011); Cardall et al. (2013) inthree dimensions. Within these schemes, the lowest moments upto order 𝑁 of the distribution function are evolved, and the flux-limited diffusion method is then the limiting case for 𝑁 = 𝑁 = In some cases the direct solution of the Boltzmann equation is indeedfeasible by exploiting symmetries and reducing the spatial dimensionalityof the problem.
Such a scheme is known in the literature as “M1 scheme” and is in-deed one of the most commonly used methods for radiative transportthroughout many different applications of computational relativisticastrophysics (Rezzolla & Miller 1994; Roedig et al. 2012; Sa¸dowskiet al. 2013; Fragile et al. 2014; McKinney et al. 2014; O’Connor2015; Foucart et al. 2015; Skinner et al. 2019; Melon Fuksman &Mignone 2019; Weih et al. 2020b). With this method, it is possibleto track the average direction of the radiation momentum, providinga significant improvement over the previously mentioned schemes,but can still lead to rather unphysical results such as those thatemerge when radiation beams interact and cross (see e.g., Fragileet al. 2014; Foucart et al. 2015; Weih et al. 2020b). Furthermore, asis typical in moment-based schemes, in the M1 approach the set ofevolution equations for the zeroth and first moment depend on thesecond moment, which is not known within this hierarchy schemeand has to be approximated in the form of a closure relation. Thisresults in the M1 scheme only being exact in the limits of highand/or low optical depths, but not in the intermediate regime. Thesituation does not improve when going to methods with
𝑁 > 𝑁 + Lattice Boltzmann method (LB method) (Krüger et al. 2017; Succi2018), which is commonly used in computational fluid dynamics asan alternate scheme to direct hydrodynamic solvers.The application of the LB method as a solver for radiationtransfer problems is relatively new (Asinari et al. 2010), and mostmodels proposed so far apply only to the analysis of steady-stateradiation-transport problems in one and two dimensions (Bindra& Patil 2012; Mishra et al. 2014; McCulloch & Bindra 2016; Yiet al. 2016), with very few studies carried out in three dimensions(McHardy et al. 2016; Wang et al. 2019). More recent develop-ments (Mink et al. 2020) have shown in detail that LB offers anaccurate and efficient tool in the diffusive regime of radiation trans-port, though struggling in the transition towards ballistic conditions.Finally, we should stress that all these previous studies deal withradiative transfer in non-relativistic regimes.We introduce a new LB solver for studying the time dependentevolution of radiation that interacts via emission, absorption andscattering with a (dynamic) background fluid. We make use ofhigh-order spherical quadrature rules, which, thanks to their high-
MNRAS000
MNRAS000 , 000–000 (2020) attice-Boltzmann methods for radiative transport in computational astrophysics order isotropy, allow to significantly extend the applicability of themethod to a wider range of kinetic regimes. We work in a special-relativistic framework and present, to the best of our knowledge, thefirst self-consistent coupled simulation of an LB solver for radiativetransport with a dynamically evolving fluid background.Our paper is structured as follows: in Sec. 2 we present a shortsummary and introduction of the classical LB method and discussits advantages, which make it ideal for radiative-transport problemsin computational astrophysics. In Sec. 3 we illustrate the details ofthe new LB method for radiative transport in special relativity andwithin a “grey” (i.e., energy averaged) approximation and show howto implement such a scheme suitably for simulations. We verify thisimplementation by a number of standard tests in Sec. 4. In Sec. 5we couple our new LB-code to a hydrodynamics code, which isrepresentative of the many relativistic-hydrodynamics codes usedin the field of numerical astrophysics, and present a simulation of arelativistic jet; we show that its dynamics changes qualitatively dueto the back-reaction of the produced radiation. Finally, we compareits accuracy and computational cost to the commonly used M1scheme in Sec. 6. We conclude and discuss future prospects of ournew method in Sec. 7.Throughout this paper, we use units with 𝑐 = In this section we provide a brief overview of the LB method.The reader already familiar with the topic may safely jump directlyto Sec. 3. We should stress that, for ease of presentation, in thissection we will summarise the conceptual steps of the derivationand algorithmic structure of LB in a non-relativistic framework.For details on the derivation of the method in special-relativity, theinterested reader is referred to a recent review (Gabbana et al. 2020).The LB method has emerged in the past decades as a compu-tationally efficient numerical tool for the simulation of the dynamicof fluids in classical hydrodynamics. Its origin can be traced back tothe pioneering work on discrete velocity models (Broadwell 1964)in the 1960s and later on to the work done on Lattice Gas CellularAutomata (Hardy et al. 1973; Frisch et al. 1986) in the late 1980s.Since then, the method has evolved as an independent and efficientalternative to direct Navier-Stokes solvers in the field of classicalcomputational fluid dynamics (McNamara & Zanetti 1988; Higueraet al. 1989; Benzi et al. 1992) and allied disciplines, primarily softmatter (Succi 2015; Dünweg & Ladd 2009).Contrary to direct hydrodynamic solvers, the LB method relieson the underlying microscopic dynamics of the fluid constituents –be them molecules or photons – and therefore the natural theoreticalframework to start approaching the method is kinetic theory andits mathematical cornerstone, the Boltzmann equation (here takenwithout source terms): (cid:18) 𝜕𝜕𝑡 + 𝒗 · ∇ (cid:19) 𝑓 ( 𝒓 , 𝒗 , 𝑡 ) = C( 𝒓 , 𝒗 , 𝑡 ) . (1)The distribution function 𝑓 ( 𝒓 , 𝒗 , 𝑡 ) refers to the number of par-ticles with velocity 𝒗 at position 𝒓 at time 𝑡 , while the collisionoperator C( 𝒓 , 𝒗 , 𝑡 ) accounts for collisions between point-particlesin the fluid and in Boltzmann’s theory takes the form of a nonlocal integral in momentum space. It is customary to replace the full collisional operator with a simplified model, such as the well-known Bhatnagar-Gross-Krook (BGK) relaxation time approxima-tion (Bhatnagar et al. 1954), encompassing the natural tendency ofthe system to relax towards an equilibrium, i.e., the tendency of 𝑓 to reach a distribution function 𝑓 eq describing a local equilibriumstate C( 𝒓 , 𝒗 , 𝑡 ) = − 𝜏 ( 𝑓 ( 𝒓 , 𝒗 , 𝑡 ) − 𝑓 eq ( 𝒓 , 𝒗 , 𝑡 )) . (2)In the above, 𝜏 represents the typical timescale needed to reach theequilibrium, a parameter which controls the hydrodynamic transportcoefficients, hence dissipative phenomena within the fluid.For a classical fluid with massive constituents, the equilibriumfunction is represented by the Maxwell-Boltzmann distribution (see,e.g., Rezzolla & Zanotti 2013) 𝑓 eq ( 𝒓 , 𝒗 , 𝑡 ) = 𝜌 ( 𝒓 , 𝑡 ) (cid:18) 𝑚 𝜋𝑘 B 𝑇 (cid:19) 𝑑 exp (cid:20) − 𝑚𝑘 B 𝑇 ( 𝒗 − 𝒖 ( 𝒓 , 𝑡 )) (cid:21) , (3)where 𝑚 is the mass of constituent particles and 𝑘 B is the Boltzmannconstant. The rest-mass density 𝜌 ( 𝒓 , 𝑡 ) and the velocity field 𝒖 ( 𝒓 , 𝑡 ) can be computed as the zeroth and first moment of the distributionfunction, respectively, 𝜌 ( 𝒓 , 𝑡 ) : = 𝑚 ∫ 𝑓 ( 𝒓 , 𝒗 , 𝑡 ) 𝑑 𝒗 , 𝒖 ( 𝒓 , 𝑡 ) : = 𝑚𝜌 ( 𝒓 , 𝑡 ) ∫ 𝒗 𝑓 ( 𝒓 , 𝒗 , 𝑡 ) 𝑑 𝒗 . (4)Historically, the LB method was devised as a noise-free (pre-averaged) version of its lattice-gas cellular automaton (LGCA) an-cestor (McNamara & Zanetti 1988). This represented a major con-ceptual leap, but left all other LGCA shortcomings untouched, pri-marily the exponential complexity barrier associated with Booleancollision operators and the ensuing low collisional rates which pre-vented LGCA from accessing high-Reynolds regimes (turbulentflows). All of the above barriers were lifted just months later in ashort sequence of papers (Higuera et al. 1989; Higuera & Jiménez1989; Higuera & Succi 1989; Benzi et al. 1992), which placed LBon the map of computational fluid dynamics.For the sake of simplicity and continuity with continuum ki-netic theory, it proves expedient to derive LB from the expansion ofthe equilibrium distribution in a series of orthogonal Hermite poly-nomials 𝑯 ( 𝑘 ) ( 𝒗 ) ; the advantage of using Hermite as the expansionbasis is that the expansion coefficients 𝒂 ( 𝑘 ) coincide with the mo-ments of the distribution function (Grad 1949a,b). The expansionis then truncated to the desired order 𝑆 , high enough to recover themacroscopic observables of interest (Shan & He 1998; Shan et al.2006).In 𝑑 dimensions, this results in the following local equilibriumdistribution function: 𝑓 eq = (cid:18) 𝜋 (cid:19) 𝑑 exp (cid:18) − 𝒗 (cid:19) 𝑆 ∑︁ 𝑘 = 𝒂 ( 𝑘 ) ( 𝜌, 𝒖 ) 𝑯 ( 𝑘 ) ( 𝒗 ) , (5)where all physical quantities have been made dimensionless byappropriate scaling with respect to a characteristic velocity ˜ 𝑐 : = √︁ 𝑘 B 𝑇 / 𝑚 , a characteristic temperature 𝑇 , a mass unit 𝑚 , and alength scale 𝐿 . Furthermore, the (microscopic) velocity vector ofphase-space is discretized using a set of 𝑁 pop distinct populations 𝒗 𝑖 with 𝑖 ∈ [ 𝑁 pop ] . As a consequence, the distribution 𝑓 itself be-comes a set of 𝑁 pop functions 𝑓 𝑖 ( 𝒓 , 𝑡 ) = 𝑓 ( 𝒓 , 𝒗 𝑖 , 𝑡 ) , each accountingfor the particles moving along the discrete direction 𝒗 𝑖 .The choice of the discrete velocities lies at the heart of the LBmethod, the informing criterion being of reproducing exactly the setof kinetic moments which describe the low-Knudsen hydrodynamic MNRAS , 000–000 (2020)
Weih et al. regime, namely the mass density (scalar), the flow current (vector)and the momentum flux tensor (second-order tensor). In the above,"exactly" means that no error altogether is incurred by replacing theintegrals in continuum velocity space with the corresponding sum-mations over the discrete velocities which characterise the LB rep-resentation. Formally, this can be linked to a Gauss-Hermite quadra-ture rule, where one defines 𝒗 𝑖 , and the corresponding weights 𝑤 𝑖 ,and requires exact preservation of the relevant hydrodynamic fields.In equations: 𝜌 ( 𝒓 , 𝑡 ) = 𝑁 pop ∑︁ 𝑖 = 𝑓 𝑖 ( 𝒓 , 𝑡 ) , 𝒖 ( 𝒓 , 𝑡 ) = 𝜌 ( 𝒓 , 𝑡 ) 𝑁 pop ∑︁ 𝑖 = 𝒗 𝑖 𝑓 𝑖 ( 𝒓 , 𝑡 ) , (6)with the truncated equilibrium distribution Eq. (5) given by 𝑓 eq 𝑖 = 𝑤 𝑖 𝑆 ∑︁ 𝑘 = 𝒂 ( 𝑘 ) ( 𝜌, 𝒖 ) 𝑯 ( 𝑘 ) ( 𝒗 𝑖 ) . (7)The combination of the velocity discretization with explicit time-marching finally delivers the lattice Boltzmann equation: 𝑓 𝑖 ( 𝒓 + 𝒗 𝑖 Δ 𝑡, 𝑡 + Δ 𝑡 ) = 𝑓 𝑖 ( 𝒓 , 𝑡 ) + Δ 𝑡 C 𝑖 ( 𝒓 , 𝑡 ) , (8)with Δ 𝑡 the time-step, Δ 𝑥 = 𝒗 𝑖 Δ 𝑡 the characteristic mesh spacing,and C 𝑖 = 𝜏 ( 𝑓 𝑖 − 𝑓 eq 𝑖 ) . (9)The evolution of Eq. (8) follows the so-called “stream-and-collide”paradigm, where in the “collide step” each population 𝑓 𝑖 ( 𝒓 , 𝑡 ) isupdated by receiving a local collisional contribution: 𝑓 ∗ 𝑖 ( 𝒓 , 𝑡 ) = 𝑓 𝑖 ( 𝒓 , 𝑡 ) + Δ 𝑡 C 𝑖 ( 𝒓 , 𝑡 ) . (10)In the streaming step, instead, the post collision populations 𝑓 ∗ 𝑖 ( 𝒓 , 𝑡 ) stream along their associated direction 𝒗 𝑖 , landing on the corre-sponding neighbouring lattice site (no particle can fly off-grid): 𝑓 𝑖 ( 𝒓 + 𝒗 𝑖 Δ 𝑡, 𝑡 + Δ 𝑡 ) = 𝑓 ∗ 𝑖 ( 𝒓 , 𝑡 ) . (11)Two major assets associated to the stream-collide paradigmare worth highlighting. First, the non-local operator (streaming) islinear and the nonlinear one (collision) is local, meaning that, atvariance with the hydrodynamic representation, non-linearity andnon-locality are disentangled. This is because information alwaystravels along constant characteristics, the discrete velocities, regard-less of the spacetime complexity of the emergent hydrodynamics. Bycontrast, in the fluid representation, information travels along space-time dependent material lines, defined by the local flow speed. Thisis a major advantage also for the handling of complex boundaryconditions and parallel computing.Second, since dissipation emerges from enslaving of the Boltz-mann distribution to local equilibrium, there is no need for secondorder spatial derivatives. This is a significant advantage for the cal-culation of the stress tensor, especially near solid boundaries. Inaddition, since the collision operator is conservative to machine-accuracy, the LB method usually offers better accuracy than mostgrid-based discretisations of the Laplace operator.The standard LB method described so far is suitable for thedescription of hydrodynamic systems, where the molecular meanfree path is much shorter than the shortest hydrodynamic lengthscale, the ratio of the two being the Knudsen number of the fluid.On the other hand, for a fluid of radiation particles, e.g., photons orneutrinos, one should bear in mind that the radiation constituentsdo not interact among themselves, but only with the background fluid, which effectively produces, destroys, and scatters the radia-tion particles. Hence, when applying the LB method to radiationlocal conservation laws must be revisited and complemented withsuitable source (sink) terms, accounting for the above processes.In particular, energy and momentum are conserved only over thecombined system of radiation and fluid.This still fits within the LB framework, which has been used fordecades to study transport phenomena, such as advection-diffusion-reaction equations, (Massaioli et al. 1993; He et al. 1998; Penget al. 2003; Karlin et al. 2013). Based on this idea, several LBmodels have been proposed to study radiative transport. Initiallymost of these efforts have focused on studying steady-state problemsin one and two spatial dimensions (Asinari et al. 2010; Bindra& Patil 2012; Mishra et al. 2014; McCulloch & Bindra 2016; Yiet al. 2016), considering isotropic as well as anisotropic scattering(Vernekar & Mishra 2014). So far, however, only very few authorshave considered the three-dimensional case (McHardy et al. 2016;Mink et al. 2020).All the models mentioned above discretize the velocity spaceby means of standard space-filling lattices, typically used in LBmethods. As a consequence, information travels in a single timestepto nodes located at different distances, corresponding to a differ-ent magnitude of the discrete velocities, a typical example beingthe two-dimensional nine-velocity lattice comprising a rest particlewith zero speed, four particles connecting to the nearest neighbours(speed 1) and four connecting the diagonals (speed √ The RTE is the master equation describing how radiation propa-gates through a medium that scatters, absorbs and emits radiationparticles. As such, the RTE follows from the Boltzmann equation as-suming massless particles (i.e., photons or neutrinos). All radiativefields are expressed in terms of the distribution function 𝑓 𝜈 ( 𝒙 , ˆ 𝒏 , 𝜈 ) of neutrinos or photons within a given frequency band d 𝜈 at position 𝒙 and velocity within a solid angle d Ω in the direction ˆ 𝒏 . MNRAS000
Weih et al. regime, namely the mass density (scalar), the flow current (vector)and the momentum flux tensor (second-order tensor). In the above,"exactly" means that no error altogether is incurred by replacing theintegrals in continuum velocity space with the corresponding sum-mations over the discrete velocities which characterise the LB rep-resentation. Formally, this can be linked to a Gauss-Hermite quadra-ture rule, where one defines 𝒗 𝑖 , and the corresponding weights 𝑤 𝑖 ,and requires exact preservation of the relevant hydrodynamic fields.In equations: 𝜌 ( 𝒓 , 𝑡 ) = 𝑁 pop ∑︁ 𝑖 = 𝑓 𝑖 ( 𝒓 , 𝑡 ) , 𝒖 ( 𝒓 , 𝑡 ) = 𝜌 ( 𝒓 , 𝑡 ) 𝑁 pop ∑︁ 𝑖 = 𝒗 𝑖 𝑓 𝑖 ( 𝒓 , 𝑡 ) , (6)with the truncated equilibrium distribution Eq. (5) given by 𝑓 eq 𝑖 = 𝑤 𝑖 𝑆 ∑︁ 𝑘 = 𝒂 ( 𝑘 ) ( 𝜌, 𝒖 ) 𝑯 ( 𝑘 ) ( 𝒗 𝑖 ) . (7)The combination of the velocity discretization with explicit time-marching finally delivers the lattice Boltzmann equation: 𝑓 𝑖 ( 𝒓 + 𝒗 𝑖 Δ 𝑡, 𝑡 + Δ 𝑡 ) = 𝑓 𝑖 ( 𝒓 , 𝑡 ) + Δ 𝑡 C 𝑖 ( 𝒓 , 𝑡 ) , (8)with Δ 𝑡 the time-step, Δ 𝑥 = 𝒗 𝑖 Δ 𝑡 the characteristic mesh spacing,and C 𝑖 = 𝜏 ( 𝑓 𝑖 − 𝑓 eq 𝑖 ) . (9)The evolution of Eq. (8) follows the so-called “stream-and-collide”paradigm, where in the “collide step” each population 𝑓 𝑖 ( 𝒓 , 𝑡 ) isupdated by receiving a local collisional contribution: 𝑓 ∗ 𝑖 ( 𝒓 , 𝑡 ) = 𝑓 𝑖 ( 𝒓 , 𝑡 ) + Δ 𝑡 C 𝑖 ( 𝒓 , 𝑡 ) . (10)In the streaming step, instead, the post collision populations 𝑓 ∗ 𝑖 ( 𝒓 , 𝑡 ) stream along their associated direction 𝒗 𝑖 , landing on the corre-sponding neighbouring lattice site (no particle can fly off-grid): 𝑓 𝑖 ( 𝒓 + 𝒗 𝑖 Δ 𝑡, 𝑡 + Δ 𝑡 ) = 𝑓 ∗ 𝑖 ( 𝒓 , 𝑡 ) . (11)Two major assets associated to the stream-collide paradigmare worth highlighting. First, the non-local operator (streaming) islinear and the nonlinear one (collision) is local, meaning that, atvariance with the hydrodynamic representation, non-linearity andnon-locality are disentangled. This is because information alwaystravels along constant characteristics, the discrete velocities, regard-less of the spacetime complexity of the emergent hydrodynamics. Bycontrast, in the fluid representation, information travels along space-time dependent material lines, defined by the local flow speed. Thisis a major advantage also for the handling of complex boundaryconditions and parallel computing.Second, since dissipation emerges from enslaving of the Boltz-mann distribution to local equilibrium, there is no need for secondorder spatial derivatives. This is a significant advantage for the cal-culation of the stress tensor, especially near solid boundaries. Inaddition, since the collision operator is conservative to machine-accuracy, the LB method usually offers better accuracy than mostgrid-based discretisations of the Laplace operator.The standard LB method described so far is suitable for thedescription of hydrodynamic systems, where the molecular meanfree path is much shorter than the shortest hydrodynamic lengthscale, the ratio of the two being the Knudsen number of the fluid.On the other hand, for a fluid of radiation particles, e.g., photons orneutrinos, one should bear in mind that the radiation constituentsdo not interact among themselves, but only with the background fluid, which effectively produces, destroys, and scatters the radia-tion particles. Hence, when applying the LB method to radiationlocal conservation laws must be revisited and complemented withsuitable source (sink) terms, accounting for the above processes.In particular, energy and momentum are conserved only over thecombined system of radiation and fluid.This still fits within the LB framework, which has been used fordecades to study transport phenomena, such as advection-diffusion-reaction equations, (Massaioli et al. 1993; He et al. 1998; Penget al. 2003; Karlin et al. 2013). Based on this idea, several LBmodels have been proposed to study radiative transport. Initiallymost of these efforts have focused on studying steady-state problemsin one and two spatial dimensions (Asinari et al. 2010; Bindra& Patil 2012; Mishra et al. 2014; McCulloch & Bindra 2016; Yiet al. 2016), considering isotropic as well as anisotropic scattering(Vernekar & Mishra 2014). So far, however, only very few authorshave considered the three-dimensional case (McHardy et al. 2016;Mink et al. 2020).All the models mentioned above discretize the velocity spaceby means of standard space-filling lattices, typically used in LBmethods. As a consequence, information travels in a single timestepto nodes located at different distances, corresponding to a differ-ent magnitude of the discrete velocities, a typical example beingthe two-dimensional nine-velocity lattice comprising a rest particlewith zero speed, four particles connecting to the nearest neighbours(speed 1) and four connecting the diagonals (speed √ The RTE is the master equation describing how radiation propa-gates through a medium that scatters, absorbs and emits radiationparticles. As such, the RTE follows from the Boltzmann equation as-suming massless particles (i.e., photons or neutrinos). All radiativefields are expressed in terms of the distribution function 𝑓 𝜈 ( 𝒙 , ˆ 𝒏 , 𝜈 ) of neutrinos or photons within a given frequency band d 𝜈 at position 𝒙 and velocity within a solid angle d Ω in the direction ˆ 𝒏 . MNRAS000 , 000–000 (2020) attice-Boltzmann methods for radiative transport in computational astrophysics The subscript 𝜈 highlights the dependency on frequency, todistinguish from the corresponding frequency-independent quanti-ties to be introduced later on. The RTE describes the evolution ofthe radiation distribution function in the direction ˆ 𝒏 𝑐 𝜕 𝑓 𝜈 𝜕𝑡 + ˆ 𝒏 · ∇ 𝑓 𝜈 = − 𝜅 𝑎,𝜈 𝑓 𝜈 + 𝜂 𝜈 + C scat = : C rad . (12)This expression is equivalent to Eq. (1), but for massless particlesand with an explicit expression for the collisional operator C rad on the right-hand side. This operator splits into three terms: theabsorption, the emission, and the scattering term, respectively. Theabsorption term is proportional to the absorption coefficient 𝜅 𝑎,𝜈 ,the emission term to the emissivity 𝜂 𝜈 , and the scattering term C scat is written in its most general form as (Bruenn 1985; Rampp 1997) C scat = ∫ ∞ 𝜈 (cid:48) d 𝜈 (cid:48) ∫ 𝜋 𝑓 (cid:48) 𝜈 (cid:48) ( − 𝑓 𝜈 ) 𝑅 in − 𝑓 𝜈 ( − 𝑓 (cid:48) 𝜈 (cid:48) ) 𝑅 out d Ω (cid:48) , (13)where 𝑅 in ( 𝜈, 𝜈 (cid:48) ) and 𝑅 out ( 𝜈, 𝜈 (cid:48) ) are the incoming and outgoingscattering kernels, respectively.These kernels depend on the underlying physical process andgenerally do not allow for an analytic solution of the scatteringintegral. In order to simplify the above integral, the incoming andoutgoing scattering kernels are typically expanded as a Legendreseries and truncated to the first two terms (Bruenn 1985; Rampp1997; Shibata et al. 2011). This leads to 𝑅 in / out ( 𝜈, 𝜈 (cid:48) ) ≈ Φ in / out0 ( 𝜈, 𝜈 (cid:48) ) + Φ in / out1 ( 𝜈, 𝜈 (cid:48) ) cos 𝜃 , (14)where Φ in / out ℓ = , are the ℓ -th coefficients of the Legendre expansionand 𝜃 is the angle between the incoming and outgoing particle,cos 𝜃 = ˆ 𝒏 · ˆ 𝒏 (cid:48) .Hereafter, we will consider only iso-energetic scatterings,i.e., we assume that the energy of the radiation particles is leftunchanged by the scattering with the constituents of the underlyingfluid.It follows that: Φ in ℓ ( 𝜈, 𝜈 (cid:48) ) ≡ Φ out ℓ ( 𝜈, 𝜈 (cid:48) ) (cid:67) Φ ℓ ( 𝜈 ) 𝛿 ( 𝜈 − 𝜈 (cid:48) ) .In this way, inserting Eq. (14) in Eq. (13), the scattering termin Eq. (12) reads as: C scat ≈ − 𝜅 ,𝜈 𝑓 𝜈 + 𝜅 ,𝜈 𝐸 𝜈 + 𝜅 ,𝜈 ˆ 𝒏 · 𝑭 𝜈 , (15)where we have used the definition of the zeroth and first moment ofthe distribution function 𝐸 𝜈 (cid:66) 𝜋 ∫ 𝜋 𝑓 𝜈 d Ω , 𝑭 𝜈 (cid:66) 𝜋 ∫ 𝜋 ˆ 𝒏 𝑓 𝜈 d Ω , (16)and defined the energy-dependent opacities 𝜅 ℓ,𝜈 = 𝜋𝜈 Φ ℓ . Notethat the explicit form of these opacities depends on the type ofradiation, namely, whether one is considering photons or neutrinos(see e.g., Eqs. (A.47) and (A.48) in Rampp (1997) for the case ofneutrinos or Rybicki & Lightman (1986) for photons).The frequency-integrated version of the RTE – often referredto as “grey” approximation – is obtained via multiplication by 𝜈 and integration over 𝜈 , i.e.,1 𝑐 𝜕𝐼𝜕𝑡 + ˆ 𝒏 · ∇ 𝐼 = − 𝜅 𝑎 𝐼 + 𝜂 + 𝜅 ( 𝐸 − 𝐼 ) + 𝜅 ˆ 𝒏 · 𝑭 , (17)where we used the definition of the frequency-integrated specificintensity 𝐼 (cid:66) ∫ ∞ 𝐼 𝜈 d 𝜈 (cid:66) ∫ ∞ 𝜈 𝑓 𝜈 d 𝜈 , (18) and the frequency-integrated moments 𝐸 (cid:66) ∫ ∞ 𝜈 𝐸 𝜈 d 𝜈 = 𝜋 ∫ 𝜋 𝐼 d Ω , (19) 𝑭 (cid:66) ∫ ∞ 𝜈 𝑭 𝜈 d 𝜈 = 𝜋 ∫ 𝜋 ˆ 𝒏 𝐼 d Ω . (20)These moments can be interpreted as the radiation energy densityand momentum density, respectively, and arguably represent themost important properties of the radiation field.Finally, the energy-averaged opacities 𝜅 𝑎 , 𝜅 , 𝜅 are given by 𝜅 ∗ (cid:66) ∫ ∞ 𝜅 ∗ ,𝜈 𝐼 𝜈 d 𝜈 ∫ ∞ 𝐼 𝜈 d 𝜈 , (21)where ∗ = 𝑎, , 𝜂 is givenby 𝜂 (cid:66) ∫ ∞ 𝜈 𝜂 𝜈 d 𝜈 . (22) We next describe the steps needed to derive an LB-inspired dis-cretization of the RTE within the grey approximation. We start byrecasting Eq. (17) in a BGK-like form:1 𝑐 𝜕𝐼𝜕𝑡 + ˆ 𝒏 · ∇ 𝐼 = − 𝜅 (cid:0) 𝐼 − 𝐼 eq (cid:1) + 𝑆 , (23)where the source term 𝑆 , accounting for emission and absorption,is given by 𝑆 = − 𝜅 𝑎 𝐼 + 𝜂 , (24)while the scattering term has been rearranged by introducing theequilibrium radiation intensity 𝐼 eq : = 𝐸 + 𝜆 ˆ 𝒏 · 𝑭 , 𝜆 = 𝜅 𝜅 . (25)We then perform a discretization of the velocity space, namely,the directions ˆ 𝒏 along which radiation propagates at the speedof light. In essence, we replace ˆ 𝒏 with a discrete set of 𝑁 pop directions ˆ 𝒏 𝑖 , which define the corresponding discrete intensities 𝐼 𝑖 ( 𝒙 , 𝑡 ) = 𝐼 ( 𝒙 , ˆ 𝒏 𝑖 , 𝑡 ) and their associated weights 𝑤 𝑖 (a more de-tailed description on how to choose the 𝑁 pop directions is given inSec 3.2.1). Consequently, Eq. (23) splits into a set of 𝑁 pop equa-tions, each describing the evolution of a specific intensity 𝐼 𝑖 via the 𝑖 -th expression of the RTE1 𝑐 𝜕𝐼 𝑖 𝜕𝑡 + ˆ 𝒏 𝒊 · ∇ 𝐼 𝑖 = − 𝜅 (cid:16) 𝐼 𝑖 − 𝐼 eq 𝑖 (cid:17) + 𝑆 𝑖 . (26)The discrete counterparts of the integrals in Eq. (19) and Eq. (20)then take the following form: 𝐸 ≈ 𝑁 pop ∑︁ 𝑖 = 𝐼 𝑖 , 𝑭 ≈ 𝑁 pop ∑︁ 𝑖 = ˆ 𝒏 𝑖 𝐼 𝑖 , (27)while higher moments can be computed accordingly as 𝑀 𝑗 ... 𝑗 𝑚 ≈ 𝑁 pop ∑︁ 𝑖 = 𝑛 𝑗 𝑖 . . . 𝑛 𝑗 𝑚 𝑖 𝐼 𝑖 . (28)Finally, a first-order discretization in time with a time-step Δ 𝑡 leads to the radiative Lattice Boltzmann equation: 𝐼 𝑖 ( 𝒓 + 𝑐 ˆ 𝒏 𝑖 Δ 𝑡, 𝑡 + Δ 𝑡 ) = 𝐼 𝑖 ( 𝒓 , 𝑡 ) − 𝑐𝜅 Δ 𝑡 (cid:16) 𝐼 𝑖 ( 𝒓 , 𝑡 ) − 𝐼 eq 𝑖 ( 𝒓 , 𝑡 ) (cid:17) + 𝑐 Δ 𝑡𝑆 𝑖 ( 𝒓 , 𝑡 ) , (29) MNRAS , 000–000 (2020)
Weih et al. where 𝐼 eq 𝑖 : = 𝑤 𝑖 ( 𝐸 ( 𝒓 , 𝑡 ) + 𝜆 ˆ 𝒏 𝒊 · 𝑭 ( 𝒓 , 𝑡 )) , 𝑆 𝑖 : = 𝑤 𝑖 𝜂 − 𝜅 𝑎 𝐼 𝑖 ( 𝒓 , 𝑡 ) . (30) The key of the LB method lies within the velocity discretization,that is, the mapping of the continuum velocity space in terms of afinite set of discrete velocities { 𝒗 𝑖 } , the so-called discrete velocitystencil.Indeed, the definition of the discrete velocity set, and its asso-ciated weights, plays a crucial role in the LB method, and severalapproaches have been developed over the years. The one based onthe Gauss-Hermite quadrature is arguably the most systematic (Shanet al. 2006; Philippi et al. 2006; Shan 2016). Another way to go,which has been used in the early days of LB, is to construct velocitysets as 𝑑 -dimensional projections from known ( 𝑑 + ) velocity setsd ' Humières et al. (1986).Yet, another possible approach consists in defining generalconditions that should be satisfied by the velocity set, in terms ofsymmetry and conservation, typically mass, momentum and mo-mentum flux (isotropy).As mentioned above, since we cannot rely on conservationlaws that would allow us to derive quadrature rules in a systematicway, we need to define the velocity stencil on the basis of symmetryconsiderations and isotropy conditions, privileging those stencilsthat exhibit a sufficiently high order of isotropy, so as to handlethe different kinematic regimes that are typically encountered inastrophysical scenarios.Formally, we define 𝑘 -rank tensors 𝑇 𝛼 ...𝛼 𝑘 , that are con-structed by combining all the products between the different direc-tions 𝑛 𝑖 forming the velocity stencil (appropriately weighted withthe weights 𝑤 𝑖 ) 𝑇 𝛼 ...𝛼 𝑘 (cid:66) ∑︁ 𝑖 𝑤 𝑖 𝑛 𝛼 𝑖 . . . 𝑛 𝛼 𝑘 𝑖 . (31)The microscopic relations that have to be satisfied by the latticein order to ensure n-th order isotropy are given by (Rivet & Boon2001) 𝑇 𝛼 ...𝛼 𝑘 = , ∝ (cid:205) perm (cid:0) 𝛿 𝛼 𝛼 . . . 𝛿 𝛼 𝑘 − 𝛼 𝑘 (cid:1) k even , (32)which need to be satisfied for all 𝑘 ≤ 𝑛 .One additional condition on the definition of the velocity sten-cil is that all the (pseudo)-particles travel at the same speed, i.e., thespeed of light, so 𝒗 𝑖 = 𝑐 ˆ 𝒏 𝑖 . It follows that the discrete directions ˆ 𝒏 𝑖 must have the same magnitude, hence span the surface of a spherein three dimensions (3D), or a circle in two dimensions (2D).The adoption of spherical velocity stencils is not new to LB: ithas been used, for example, in LB models for the simulation of ultra-relativistic hydrodynamics. In such frameworks, since the interestis restricted to the hydrodynamic picture, it is still possible to definestencils which live on the intersection between a Cartesian grid anda sphere of fixed radius (Mendoza et al. 2013; Gabbana et al. 2018),thus preserving a very desirable property of LB: exact-streaming.On the other hand, going beyond hydrodynamics generally requiresa much larger number of discrete directions, making the definitionof on-lattice quadratures impractical. In these cases, it is thereforenecessary to take into consideration off-lattice schemes (Coelhoet al. 2018; Ambruş & Blaga 2018). Since we need to properly model free-streaming regimes, we adopt this latter approach andwork with off-lattice stencils spanning a unit sphere.Moreover, the choice of the velocity set should be such tomaximise the accuracy in the calculation of the moments of thespecific intensity 𝐼 ( 𝒙 , ˆ 𝒏 , 𝑡 ) , such as the energy density 𝐸 ( ˆ 𝒙 , 𝑡 ) andthe momentum density 𝑭 ( ˆ 𝒙 , 𝑡 ) .In practice, we request the discrete sums in Eq. (27) to correctlyreproduce their continuous counterparts, i.e., Eq. (19) and Eq. (20).These are spherical integrals of the form 𝑄 ( ℎ ) = 𝜋 ∫ 𝜋 ℎ ( ˆ 𝒏 ) 𝑑 Ω ≈ 𝑁 pop ∑︁ 𝑖 = 𝑤 𝑖 ℎ ( ˆ 𝒏 𝒊 ) , (33)where ℎ ( ˆ 𝒏 ) represents a generic function of the direction ˆ 𝒏 . To beginwith, we recall that any square integrable function can be expandedon the unit sphere as a series of orthogonal spherical harmonics(Atkinson & Han 2012) ℎ ( ˆ 𝒏 ) = ℎ ( 𝜃, 𝜙 ) = +∞ ∑︁ ℓ = ℓ ∑︁ 𝑚 = − ℓ 𝑐 ℓ𝑚 𝑌 𝑚ℓ ( 𝜃, 𝜙 ) , (34)where the convergence rate depends on the coefficients 𝑐 ℓ𝑚 . Thespherical quadrature rule determining the weights 𝑤 𝑖 and discretedirections ˆ 𝒏 𝒊 = ( 𝜃 𝑖 , 𝜙 𝑖 ) is then said to be of order 𝑝 if it integratesexactly all the spherical harmonics 𝑌 𝑚ℓ ( 𝜃, 𝜙 ) up to the degree ℓ = 𝑝 𝜋 ∫ 𝜋 𝑌 𝑚ℓ ( 𝜃, 𝜙 ) 𝑑 Ω = 𝑁 pop − ∑︁ 𝑘 = 𝑤 𝑘 𝑌 𝑚ℓ ( 𝜃 𝑘 , 𝜙 𝑘 ) ∀ ℓ ≤ 𝑝 . (35)It follows that a spherical quadrature of order 𝑝 integrates exactlyall integrals 𝑄 ( ℎ ) of functions ℎ ( 𝜃, 𝜙 ) described by linear com-binations of the first 𝑝 harmonics. On the other hand, functions ℎ ( 𝜃, 𝜙 ) that contain in their series harmonics of higher order areonly approximated by the quadrature, with errors that depend onthe smoothness of the function itself. Note that it is possible toprove that quadratures of order 𝑝 satisfy Eq. (32) up to the level 𝑝 ,enabling in this way to evaluate the quality of the stencils with theuse of a single parameter, i.e., 𝑝 that we call the quadrature order.Determining the quadrature satisfying the conditions discussedabove is trivial in the case of two dimensions: the integral in Eq.(33) simply has to be computed on the unit circle. Furthermore, theexpansion in Eq. (34) reduces to a Fourier series since the three-dimensional spherical harmonics reduces to circular functions (sineand cosines of the only angular coordinate 𝜙 ) ℎ ( ˆ 𝒏 ) = ℎ ( 𝜙 ) = ∞ ∑︁ 𝑚 = −∞ 𝑐 𝑚 𝑒 𝑖𝑚𝜙 , (36)and one has to satisfy exactly the relations12 𝜋 ∫ 𝜋 𝑒 𝑖𝑚𝜙 𝑑𝜙 = 𝑁 pop − ∑︁ 𝑘 = 𝑤 𝑘 𝑒 𝑖𝑚𝜙 𝑘 ∀ | 𝑚 | ≤ 𝑝 . (37)As a result, to obtain a 2D quadrature of order 𝑝 , i.e., in orderto satisfy Eq. (37), it is sufficient to consider 𝑁 pop = 𝑝 + 𝜙 𝑖 = 𝜋 𝑘 / 𝑁 pop , thus ˆ 𝒏 𝑖 = [ cos ( 𝜙 𝑖 ) , sin ( 𝜙 𝑖 )] 𝑇 , all having equalweights 𝑤 𝑘 = / 𝑁 pop . A representation of such a velocity stencilwith its discrete directions ˆ 𝒏 𝑖 is presented in Fig. 1.The problem becomes considerably more involved for thethree-dimensional case, which is clearly the most relevant one interms of astrophysical applications. Indeed, the definition of quadra-ture rules on the surface of a sphere (i.e., a 2-sphere) is still an ac-tive area of research, with several different approaches coming with MNRAS000
Weih et al. where 𝐼 eq 𝑖 : = 𝑤 𝑖 ( 𝐸 ( 𝒓 , 𝑡 ) + 𝜆 ˆ 𝒏 𝒊 · 𝑭 ( 𝒓 , 𝑡 )) , 𝑆 𝑖 : = 𝑤 𝑖 𝜂 − 𝜅 𝑎 𝐼 𝑖 ( 𝒓 , 𝑡 ) . (30) The key of the LB method lies within the velocity discretization,that is, the mapping of the continuum velocity space in terms of afinite set of discrete velocities { 𝒗 𝑖 } , the so-called discrete velocitystencil.Indeed, the definition of the discrete velocity set, and its asso-ciated weights, plays a crucial role in the LB method, and severalapproaches have been developed over the years. The one based onthe Gauss-Hermite quadrature is arguably the most systematic (Shanet al. 2006; Philippi et al. 2006; Shan 2016). Another way to go,which has been used in the early days of LB, is to construct velocitysets as 𝑑 -dimensional projections from known ( 𝑑 + ) velocity setsd ' Humières et al. (1986).Yet, another possible approach consists in defining generalconditions that should be satisfied by the velocity set, in terms ofsymmetry and conservation, typically mass, momentum and mo-mentum flux (isotropy).As mentioned above, since we cannot rely on conservationlaws that would allow us to derive quadrature rules in a systematicway, we need to define the velocity stencil on the basis of symmetryconsiderations and isotropy conditions, privileging those stencilsthat exhibit a sufficiently high order of isotropy, so as to handlethe different kinematic regimes that are typically encountered inastrophysical scenarios.Formally, we define 𝑘 -rank tensors 𝑇 𝛼 ...𝛼 𝑘 , that are con-structed by combining all the products between the different direc-tions 𝑛 𝑖 forming the velocity stencil (appropriately weighted withthe weights 𝑤 𝑖 ) 𝑇 𝛼 ...𝛼 𝑘 (cid:66) ∑︁ 𝑖 𝑤 𝑖 𝑛 𝛼 𝑖 . . . 𝑛 𝛼 𝑘 𝑖 . (31)The microscopic relations that have to be satisfied by the latticein order to ensure n-th order isotropy are given by (Rivet & Boon2001) 𝑇 𝛼 ...𝛼 𝑘 = , ∝ (cid:205) perm (cid:0) 𝛿 𝛼 𝛼 . . . 𝛿 𝛼 𝑘 − 𝛼 𝑘 (cid:1) k even , (32)which need to be satisfied for all 𝑘 ≤ 𝑛 .One additional condition on the definition of the velocity sten-cil is that all the (pseudo)-particles travel at the same speed, i.e., thespeed of light, so 𝒗 𝑖 = 𝑐 ˆ 𝒏 𝑖 . It follows that the discrete directions ˆ 𝒏 𝑖 must have the same magnitude, hence span the surface of a spherein three dimensions (3D), or a circle in two dimensions (2D).The adoption of spherical velocity stencils is not new to LB: ithas been used, for example, in LB models for the simulation of ultra-relativistic hydrodynamics. In such frameworks, since the interestis restricted to the hydrodynamic picture, it is still possible to definestencils which live on the intersection between a Cartesian grid anda sphere of fixed radius (Mendoza et al. 2013; Gabbana et al. 2018),thus preserving a very desirable property of LB: exact-streaming.On the other hand, going beyond hydrodynamics generally requiresa much larger number of discrete directions, making the definitionof on-lattice quadratures impractical. In these cases, it is thereforenecessary to take into consideration off-lattice schemes (Coelhoet al. 2018; Ambruş & Blaga 2018). Since we need to properly model free-streaming regimes, we adopt this latter approach andwork with off-lattice stencils spanning a unit sphere.Moreover, the choice of the velocity set should be such tomaximise the accuracy in the calculation of the moments of thespecific intensity 𝐼 ( 𝒙 , ˆ 𝒏 , 𝑡 ) , such as the energy density 𝐸 ( ˆ 𝒙 , 𝑡 ) andthe momentum density 𝑭 ( ˆ 𝒙 , 𝑡 ) .In practice, we request the discrete sums in Eq. (27) to correctlyreproduce their continuous counterparts, i.e., Eq. (19) and Eq. (20).These are spherical integrals of the form 𝑄 ( ℎ ) = 𝜋 ∫ 𝜋 ℎ ( ˆ 𝒏 ) 𝑑 Ω ≈ 𝑁 pop ∑︁ 𝑖 = 𝑤 𝑖 ℎ ( ˆ 𝒏 𝒊 ) , (33)where ℎ ( ˆ 𝒏 ) represents a generic function of the direction ˆ 𝒏 . To beginwith, we recall that any square integrable function can be expandedon the unit sphere as a series of orthogonal spherical harmonics(Atkinson & Han 2012) ℎ ( ˆ 𝒏 ) = ℎ ( 𝜃, 𝜙 ) = +∞ ∑︁ ℓ = ℓ ∑︁ 𝑚 = − ℓ 𝑐 ℓ𝑚 𝑌 𝑚ℓ ( 𝜃, 𝜙 ) , (34)where the convergence rate depends on the coefficients 𝑐 ℓ𝑚 . Thespherical quadrature rule determining the weights 𝑤 𝑖 and discretedirections ˆ 𝒏 𝒊 = ( 𝜃 𝑖 , 𝜙 𝑖 ) is then said to be of order 𝑝 if it integratesexactly all the spherical harmonics 𝑌 𝑚ℓ ( 𝜃, 𝜙 ) up to the degree ℓ = 𝑝 𝜋 ∫ 𝜋 𝑌 𝑚ℓ ( 𝜃, 𝜙 ) 𝑑 Ω = 𝑁 pop − ∑︁ 𝑘 = 𝑤 𝑘 𝑌 𝑚ℓ ( 𝜃 𝑘 , 𝜙 𝑘 ) ∀ ℓ ≤ 𝑝 . (35)It follows that a spherical quadrature of order 𝑝 integrates exactlyall integrals 𝑄 ( ℎ ) of functions ℎ ( 𝜃, 𝜙 ) described by linear com-binations of the first 𝑝 harmonics. On the other hand, functions ℎ ( 𝜃, 𝜙 ) that contain in their series harmonics of higher order areonly approximated by the quadrature, with errors that depend onthe smoothness of the function itself. Note that it is possible toprove that quadratures of order 𝑝 satisfy Eq. (32) up to the level 𝑝 ,enabling in this way to evaluate the quality of the stencils with theuse of a single parameter, i.e., 𝑝 that we call the quadrature order.Determining the quadrature satisfying the conditions discussedabove is trivial in the case of two dimensions: the integral in Eq.(33) simply has to be computed on the unit circle. Furthermore, theexpansion in Eq. (34) reduces to a Fourier series since the three-dimensional spherical harmonics reduces to circular functions (sineand cosines of the only angular coordinate 𝜙 ) ℎ ( ˆ 𝒏 ) = ℎ ( 𝜙 ) = ∞ ∑︁ 𝑚 = −∞ 𝑐 𝑚 𝑒 𝑖𝑚𝜙 , (36)and one has to satisfy exactly the relations12 𝜋 ∫ 𝜋 𝑒 𝑖𝑚𝜙 𝑑𝜙 = 𝑁 pop − ∑︁ 𝑘 = 𝑤 𝑘 𝑒 𝑖𝑚𝜙 𝑘 ∀ | 𝑚 | ≤ 𝑝 . (37)As a result, to obtain a 2D quadrature of order 𝑝 , i.e., in orderto satisfy Eq. (37), it is sufficient to consider 𝑁 pop = 𝑝 + 𝜙 𝑖 = 𝜋 𝑘 / 𝑁 pop , thus ˆ 𝒏 𝑖 = [ cos ( 𝜙 𝑖 ) , sin ( 𝜙 𝑖 )] 𝑇 , all having equalweights 𝑤 𝑘 = / 𝑁 pop . A representation of such a velocity stencilwith its discrete directions ˆ 𝒏 𝑖 is presented in Fig. 1.The problem becomes considerably more involved for thethree-dimensional case, which is clearly the most relevant one interms of astrophysical applications. Indeed, the definition of quadra-ture rules on the surface of a sphere (i.e., a 2-sphere) is still an ac-tive area of research, with several different approaches coming with MNRAS000 , 000–000 (2020) attice-Boltzmann methods for radiative transport in computational astrophysics ∆ t = ∆ x ˆ n ˆ n ˆ n ˆ n ˆ n ˆ n ˆ n ˆ n ˆ n ˆ n ˆ n ˆ n Figure 1.
Example of two-dimensional velocity stencils showing 𝑁 pop = 𝑐 Δ 𝑡 , where Δ 𝑡 is assumed to be in units of thegrid-spacing. corresponding advantages and drawbacks, depending on the targetapplication (Beentjes 2015; Gamba et al. 2017; Gross & Atzberger2018; Lutsko & Lam 2018; Stepán, Jirí et al. 2020).Here, we consider three different types of spherical quadratureschemes:(i) Gauss-Legendre quadrature(ii) Lebedev quadrature(iii) Spherical designThe application of the Gauss-Legendre quadrature to a 2-sphere can be obtained by making use of the product quadraturerule. Exploiting the separability of spherical harmonics, the inte-grals in Eq. (35) can be expressed as the product of a circularfunction 𝑒 𝑖𝑚𝜙 and Legendre polynomials: ∫ 𝜋 𝑌 𝑚ℓ ( 𝜃, 𝜙 ) 𝑑 Ω ∝ (cid:18)∫ 𝜋 𝑃 𝑚ℓ ( cos 𝜃 ) sin 𝜃𝑑𝜃 (cid:19) (cid:18)∫ 𝜋 𝑒 𝑖𝑚𝜙 𝑑𝜙 (cid:19) , (38)where 𝑃 𝑚ℓ ( cos ( 𝜃 )) are the associated Legendre polynomials. At thispoint, the two integrals can be evaluated using two one-dimensionalquadrature rules, where the integral in the direction 𝜃 is performedwith a one-dimensional Gauss-Legendre quadrature (Hildebrand1956), while the integral in 𝜙 can be evaluated, for example, withthe trapezoidal rule. As it is apparent from the left panel of Fig.2, the Gaussian-Legendre approach generates an accumulation ofpoints near the north and south poles of the sphere, and it is lessefficient than the other two quadratures, in the sense that it requiresa larger number of points to achieve the same order of precision.The Lebedev quadrature (central panel of Fig. 2), on the otherhand, follows a different approach: instead of considering the prod-uct of two single quadratures, the integrals in Eq. (35) are used tobuild a system of nonlinear equations in the variables { 𝑤 𝑘 , 𝜃 𝑘 , 𝜙 𝑘 } .The central intuition (due to Sobolev (1962)) is that the numberof nonlinear equations can be greatly reduced by considering only the spherical harmonics of degree ≤ 𝑝 that exhibit invariance un-der all transformations that belong to a pre-determined group 𝐺 .This procedure generates quadratures that are invariant under 𝐺 ,that is, these stencils evaluate exactly both 𝑄 ( ℎ ) and 𝑄 ( 𝑔 ( ℎ )) forall elements 𝑔 ∈ 𝐺 ,and still retain the same order of integration 𝑝 . In its original definition Lebedev’s quadrature is by constructioninvariant under the octahedral group (Lebedev 1975, 1976, 1977),although quadratures based on different symmetry groups are alsopresent in the literature (see, e.g., (Ahrens & Beylkin 2009) for aquadrature based on the icosahedral group). From the central panelof Fig. 2 one appreciates that points in Lebedev’s quadratures offera more homogeneous distribution over the sphere with respect tothe Gauss product rule. Nevertheless, a few points can be seen to bealmost overlapping.Lastly, we consider the spherical-design quadrature rules, firstintroduced by Delsarte et al. (1977). This type of quadrature requiresthat the weights associated to all nodes be equal. The task is thento define a minimum set of points which integrates correctly allthe spherical harmonics up to order 𝑝 . There is no known rule forthe definition of a generic order 𝑝 spherical-design quadrature, andthe topic is indeed still object of ongoing research. Nevertheless,numerical results leading to the definition of quadrature rules up tovery high order are available online. In this work we refer to the setof stencils presented by Womersley (2018), for which we providean example in the right panel of Fig. 2. Having presented the equations to be solved and their discretiza-tion, we proceed to summarise the steps required to evolve theLB-discretized RTE.More specifically, Eq. (29) can be solved following the standardstream-and-collide approach, where the streaming step is performedfirst, yielding provisional values of the 𝑁 pop intensities, i.e., 𝐼 ∗ 𝑖 ( 𝒓 , 𝑡 ) = 𝐼 𝑖 ( 𝒓 − ˆ 𝒏 𝑖 Δ 𝑡, 𝑡 ) . (39)However, since the discrete velocities fall off-grid for the abovediscussed stencils, an interpolation is required. For simplicity weconsider a trilinear interpolation scheme, from which follows 𝐼 ∗ 𝑖 ( 𝒓 − ˆ 𝒏 𝑖 Δ 𝑡, 𝑡 ) = Δ 𝑥 Δ 𝑦 Δ 𝑧 × (cid:110) 𝐼 𝑖 ( 𝒓 − ˆ 𝒙 − ˆ 𝒚 − ˆ 𝒛 , 𝑡 ) (cid:16) Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑥𝑖 (cid:12)(cid:12)(cid:17) (cid:16) Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑦𝑖 (cid:12)(cid:12)(cid:17) (cid:16) Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑧𝑖 (cid:12)(cid:12)(cid:17) 𝐼 𝑖 ( 𝒓 − ˆ 𝒚 − ˆ 𝒛 , 𝑡 ) (cid:16) Δ 𝑥 − Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑥𝑖 (cid:12)(cid:12)(cid:17) (cid:16) Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑦𝑖 (cid:12)(cid:12)(cid:17) (cid:16) Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑧𝑖 (cid:12)(cid:12)(cid:17) 𝐼 𝑖 ( 𝒓 − ˆ 𝒙 − ˆ 𝒛 , 𝑡 ) (cid:16) Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑥𝑖 (cid:12)(cid:12)(cid:17) (cid:16) Δ 𝑦 − Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑦𝑖 (cid:12)(cid:12)(cid:17) (cid:16) Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑧𝑖 (cid:12)(cid:12)(cid:17) 𝐼 𝑖 ( 𝒓 − ˆ 𝒙 − ˆ 𝒚 , 𝑡 ) (cid:16) Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑥𝑖 (cid:12)(cid:12)(cid:17) (cid:16) Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑦𝑖 (cid:12)(cid:12)(cid:17) (cid:16) Δ 𝑧 − Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑧𝑖 (cid:12)(cid:12)(cid:17) 𝐼 𝑖 ( 𝒓 − ˆ 𝒛 , 𝑡 ) (cid:16) Δ 𝑥 − Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑥𝑖 (cid:12)(cid:12)(cid:17) (cid:16) Δ 𝑦 − Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑦𝑖 (cid:12)(cid:12)(cid:17) (cid:16) Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑧𝑖 (cid:12)(cid:12)(cid:17) 𝐼 𝑖 ( 𝒓 − ˆ 𝒚 , 𝑡 ) (cid:16) Δ 𝑥 − Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑥𝑖 (cid:12)(cid:12)(cid:17) (cid:16) Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑦𝑖 (cid:12)(cid:12)(cid:17) (cid:16) Δ 𝑧 − Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑧𝑖 (cid:12)(cid:12)(cid:17) 𝐼 𝑖 ( 𝒓 − ˆ 𝒙 , 𝑡 ) (cid:16) Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑥𝑖 (cid:12)(cid:12)(cid:17) (cid:16) Δ 𝑦 − Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑦𝑖 (cid:12)(cid:12)(cid:17) (cid:16) Δ 𝑧 − Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑧𝑖 (cid:12)(cid:12)(cid:17) 𝐼 𝑖 ( 𝒓 , 𝑡 ) (cid:16) Δ 𝑥 − Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑥𝑖 (cid:12)(cid:12)(cid:17) (cid:16) Δ 𝑦 − Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑦𝑖 (cid:12)(cid:12)(cid:17) (cid:16) Δ 𝑧 − Δ 𝑡 (cid:12)(cid:12) ˆ 𝑛 𝑧𝑖 (cid:12)(cid:12)(cid:17)(cid:111) , (40) MNRAS , 000–000 (2020)
Weih et al.
Figure 2.
Examples of three-dimensional velocity stencils, comparing the distribution of the velocity directions (purple dots) for three different quadraturetypologies; from left to right we respectively show an example for Gauss-Legendre, Lebedev, and spherical design. Analogous to Fig. 1 the purple dots arelocated on a sphere of radius 𝑐 Δ 𝑡 and the green squares denote the cell centres. with ˆ 𝒙 = (cid:2) sgn (cid:0) ˆ 𝑛 𝑥𝑖 (cid:1) Δ 𝑥, , (cid:3) 𝑇 (41) ˆ 𝒚 = (cid:2) , sgn (cid:0) ˆ 𝑛 𝑦𝑖 (cid:1) Δ 𝑦, (cid:3) 𝑇 (42) ˆ 𝒛 = (cid:2) , , sgn (cid:0) ˆ 𝑛 𝑧𝑖 (cid:1) Δ 𝑧 (cid:3) 𝑇 . (43)For the two dimensional case a bilinear interpolation is used equiv-alently.Next, the macroscopic moments are computed according toEq. (27), from which the source term and the collisional operatorare computed using Eq. (30).The collision step is then performed as follows: 𝐼 𝑖 ( 𝒓 , 𝑡 + Δ 𝑡 ) = 𝐼 ∗ 𝑖 ( 𝒓 , 𝑡 ) − 𝜅 Δ 𝑡 (cid:2) 𝐼 ∗ 𝑖 ( 𝒓 , 𝑡 ) − 𝐼 eq 𝑖 ( 𝒓 , 𝑡 ) (cid:3) + Δ 𝑡𝑆 𝑖 ( 𝒓 , 𝑡 ) . (44)Note that since in most astrophysical applications the emissivitiesand opacities vary over several orders of magnitude, these sourceterms may take values much larger than the evolved variable 𝐼 𝑖 , thusmaking the RTE a stiff equation that requires special treatment tobe solved efficiently (see, e.g., Weih et al. 2020b).Thus, rather than solving the explicit Eq. (44), we solve theimplicit form, i.e., 𝐼 𝑖 ( 𝒓 + ˆ 𝒏 𝑖 Δ 𝑡, 𝑡 + Δ 𝑡 ) = 𝐼 ∗ 𝑖 ( 𝒓 , 𝑡 )− 𝜅 Δ 𝑡 (cid:2) 𝐼 𝑖 ( 𝒓 , 𝑡 + Δ 𝑡 ) − 𝐼 eq 𝑖 ( 𝒓 , 𝑡 + Δ 𝑡 ) (cid:3) + Δ 𝑡𝑆 𝑖 ( 𝒓 , 𝑡 + Δ 𝑡 ) , (45)Because of the presence of the first two moments in the definition ofthe source term, we obtain a system of 𝑁 pop linear equations with 𝑁 pop unknowns. This could be solved via the inversion of an 𝑁 pop × 𝑁 pop -matrix at every grid-point, which is, however, computationallyunfeasible, especially in 3D, when a large number 𝑁 pop of discretevelocity directions is required. An alternate and computationallymuch more viable method is known as Lambda iteration (Rampp1997), which works as follows:1. Take as initial guess for 𝐸 ( 𝒓 , 𝑡 + Δ 𝑡 ) and 𝑭 ( 𝒓 , 𝑡 + Δ 𝑡 ) themoments computed from 𝐼 𝑖 at time 𝑡 .2. Eq. (45) constitutes a system of 𝑁 pop linear decoupled equa-tions that can be solved analytically and independently, to obtain afirst estimate of 𝐼 𝑖 at time 𝑡 + Δ 𝑡 . 3. Use this estimate of 𝐼 𝑖 ( 𝒓 , 𝑡 + Δ 𝑡 ) to formulate an improvedguess for 𝐸 ( 𝒓 , 𝑡 + Δ 𝑡 ) and 𝑭 ( 𝒓 , 𝑡 + Δ 𝑡 ) .4. Cycle back to step 2. and repeat until all 𝐼 𝑖 ( 𝒓 , 𝑡 + Δ 𝑡 ) converge,with an error below a prescribed threshold.To summarise, our scheme consists of two steps; first the streaming-step according to Eq. (40) and then the collision-step according toEq. (44), which is solved following the iterative procedure describedabove. If the radiation is coupled to a fluid, the radiative four-forcethat enters the standard equations of relativistic (magneto-) hydro-dynamics (RMHD) can be computed after every timestep from themoments 𝐸 and 𝑭 , which themselves are computed approximatelyfrom 𝐼 𝑖 according to Eq. (27). At the beginning of the next timestep,the coefficients 𝜂 , 𝜅 𝑎 , 𝜅 and 𝜅 can then be computed from theupdated fluid variables.As a result, the coupling to a standard RMHD code is exactly the same as for the commonly used M1 scheme (see Weih et al.2020b, for a detailed description of this coupling). Finally, we notethat while we restrict ourselves for simplicity to a first-order time-stepper, it is conceptually straightforward to extend the evolution tohigher orders using for example implicit-explicit (IMEX) schemes(Pareschi & Russo 2005), where the streaming-step is treated ex-plicitly and the collision-step implicitly. In astrophysical simulations, the ordinary fluid interacting with theradiation features optical depths that vary considerably, rangingfrom the optically thin regime – where radiation is in free streaming– to the optically thick regime – where radiation is coupled with thefluid and propagates by diffusion. While in several studies only oneof these regimes is considered, (see e.g., Fragile et al. 2012; Roediget al. 2012), it is our aim to provide a numerical method that canhandle both limits, as well as the intermediate regime. The latteris particularly difficult to be described accurately by moment-basedschemes, mostly because the closure relation – which is essentialand inevitable in these schemes – is normally defined in either the Special attention has to be paid to the case of a moving background fluid(see Sec. 5 and appendix B). MNRAS , 000–000 (2020) attice-Boltzmann methods for radiative transport in computational astrophysics − . . . y | F | =1 − .
25 0 .
00 0 . x − . . . y − .
25 0 .
00 0 . x . . . . . . E/E max
Figure 3.
Left:
Straight beam of radiation with a CFL-number of 1 . . 𝑡 = . Right:
The sameas the left panel, but for a beam propagating diagonally. The region with 𝑥, 𝑦 < − .
25 is frozen via a boundary condition in order to continuouslyshoot the beam into the grid from the bottom left. optically thin or the optically thick regime and is then interpolatedbetween these two limits (Weih et al. 2020b).In the series of tests presented below, we discuss the perfor-mance of the LB method in these different limits, starting in Sec.4.1 with the optically thin one – which represents the most difficultchallenge. The optically thick regime is tested in Sec. 4.2 – wherethe LB method performs extremely well. Finally, in Sec. 4.2.2, wepresent an example of the intermediate regime, where LB is shownto provide very accurate results, at variance with the commonlyused M1 scheme.
Since the LB method is designed for collisional fluids, we beginthe verification of the LB method and its implementation with anumber of code tests in the most difficult regime, namely, the one inwhich the radiation is actually freely streaming, as is the case when 𝜂 = 𝜅 𝑎 = 𝜅 = 𝜅 = We start with a “classical” beam test, namely, the propagation ofa well defined beam of radiation injected from the boundary of − . . . y | F | = 1 LB − . − . . . . x − . . . y M1 . . . . . . E/E max
Figure 4.
Comparison of the performance for the beam-crossing problembetween LB (top) and M1 (bottom). the computational domain. Having set 𝜂 = 𝜅 𝑎 = 𝜅 = 𝜅 = 𝐼 𝑖 is propagated from one grid-cell to thenext following the underlying stencil. The top-left panel of Fig. 3shows an example of a freely streaming beam on a grid of size − . < 𝑥 < . , − . < 𝑦 < .
5, which we cover with 100 equal-size grid cells. This test is performed in 2D using the stencil shownin Fig. 1 with 𝑁 pop = 𝐼 𝑖 = (cid:40) 𝑖 = 𝑖 ≠ , (46)at 𝑥 = − . | 𝑦 | < . 𝑥 -axis fromleft to right at the speed of light. In Fig. 3 we show with a colorcodethe radiation energy density and with arrows the momentum density,as computed according to Eq. (27) at time 𝑡 = .
7, and where wehave used a timestep Δ 𝑡 = Δ 𝑥 = Δ 𝑦 . This timestep is a standardchoice in classical LB methods and leads to perfect streaming,i.e., no diffusion along the beam’s direction of propagation. Thisis due to the fact that the streaming step, i.e., Eqs. (39) and (40),reduces to 𝐼 ∗ 𝑖 ( 𝒓 , 𝑡 + Δ 𝑡 ) = 𝐼 𝑖 ( 𝒓 − ˆ 𝒏 𝑖 Δ 𝑡, 𝑡 ) = | ˆ 𝑛 𝑥𝑖 | Δ 𝑥 × 𝐼 𝑖 ( 𝒓 − ˆ 𝒙 , 𝑡 ) . (47)Considering that | ˆ 𝑛 𝑥𝑖 |/ Δ 𝑥 =
1, where | ˆ 𝑛 𝑥𝑖 | is the 𝑥 -th component ofthe 𝑖 -th velocity vector, we find that the intensity is simply propa-gated from one cell to its right neighbour during each iteration.Most of the astrophysical codes solving the equations ofRMHD employ either finite-volume or finite-differencing schemesand require a timestep Δ 𝑡 = CFL Δ 𝑥 for numerical stability, whereCFL ≤ MNRAS , 000–000 (2020) Weih et al. cal systems, the LB method needs to be coupled to such a RMHDcode, we also perform the above beam tests with Δ 𝑡 = . Δ 𝑥 , whereCFL = . 𝐼 , which correspondsto ˆ 𝒏 = [ cos ( 𝜋 / ) , sin ( 𝜋 / )] , rather than 𝐼 to a nonzero value; ofcourse 𝐼 𝑖 ≠ = Δ 𝑡 = Δ 𝑥 (top panel) and Δ 𝑡 = . Δ 𝑥 (bot-tom panel). In both cases, the beam now diffuses much more, aneffect that can also be observed for the commonly used moment-schemes (Weih et al. 2020b). It is worth remarking that in all ofthe above tests, the beams propagate along one of the discrete ve-locity directions. However, radiation beams can propagate in anydirection in the continuum, not necessarily along a discretized ve-locity direction. In this case, the LB scheme would inevitably incurincreasing errors as we approach the optically thin regime. Thesecan be tamed by developing high-order phase-space interpolators,possibly involving non-local neighbours both in configuration andvelocity space. Clearly, this exposes a tension between accuracy andefficiency which still needs to be explored and resolved in full. In thispaper, we rely upon trilinear and nearest-neighbour interpolation inconfiguration and velocity space, respectively.Next, we consider the performance of our code for another“classical” and yet fundamental free-streaming test: two crossingbeams. Assuming the radiation to consist of photons or neutrinosof the same flavour, one would expect the two beams to cross eachother without interacting. The M1 scheme is known to perform verypoorly under these conditions (Fragile et al. 2014; McKinney et al.2014; Foucart et al. 2015; Rivera-Paleo & Guzmán 2019; Weih et al.2020b), due to the fact that it retains only the two lowest moments,hence only the average direction of propagation.This information could in principle be obtained by employinghigher moments, at the cost of increased computational costs. On theother hand, in an LB method the various directions of propagationare evolved separately and thus crossing beams can be correctlyevolved in time. This is shown in Fig. 4, where the results from anLB simulation (top panel) are compared to the ones from an M1code (bottom panel). In both cases, we perform the simulation on agrid of size − . < 𝑥 < . − . < 𝑦 < .
25 with a resolutionof 200 ×
100 and choose again a CFL coefficient of 0 .
2. For the M1solution, we initialise the momentum densities of the two beamsto point towards the center, while for LB we simply initialise theintensities for the corresponding directions.It can be seen clearly from Fig. 4 that in the LB simulationthe two radiation beams cross as expected, while they merge toan averaged beam when using M1. The proper description of thisbehaviour is of crucial importance in astrophysical simulations,where beams of radiation – such as those emitted from the torus ofa binary neutron-star merger remnant or in a supernova explosion –are expected to meet and interact. Hence, the successful outcome of this test provides encouraging evidence that the present LB methodcan be used in the place of moment-based schemes.
In order to evaluate if radiation is propagated isotropically by thenumerical scheme, we show here the results obtained for the case ofa spherically symmetric propagation of a radiation wave. The waveis expected to expand at the speed of light in all directions, with itsenergy density decreasing over time following an inverse-square lawfor the distance travelled by the wavefront. We run all simulationson a 3D grid with 200 uniform sized grid cells with Δ 𝑡 = .
2, andinitialise the intensities by setting 𝐼 𝑖 = 𝑅 = Δ 𝑥 grid-cells and 𝐼 𝑖 = ( 𝑥, 𝑦 ) and ( 𝑦, 𝑧 ) planes. It is clear that a low number of discrete velocities 𝑁 pop in the stencils does not allow the radiation to stream isotropically.However, upon increasing 𝑁 pop , the isotropy of the system rapidlyimproves up to very satisfactory levels. Also to be noted, none ofthe three methods emerges as a neat winner, all yielding results ofsimilar quality. This is true even for the Lebedev quadrature, whichfeatures a consistently higher quadrature order as the other two typesfor the same 𝑁 pop .Finally, we show that the wave’s maximum energy decays pro-portionally to 𝑟 − , which is expected for this 3D test, where theenergy is conserved and spreads over a spherical shell of area 𝜋𝑟 in time. Figure 6 shows this behaviour for the reference case of theLebedev quadrature of order 𝑝 =
23 (middle panel in Fig. 5). Thefigure reports the profile of the radiation energy density along adiagonal cut at different times. Note that the maxima of these cutsalign well with an inverse square law (red-dashed line).
Having tested the LB method for the solution of the RTE in theoptically thin regime, we now focus on the optically thick regime.We recall that an optically thick medium either absorbs (in the caseof a high value of the absorption opacity 𝜅 𝑎 ) or scatters radiation(in the case of a high scattering opacities 𝜅 and 𝜅 ). As a result, wepresent tests of the absorption scenario in Secs. 4.2.1 and 4.2.2 andthe scattering scenario, i.e., the diffusion limit, in Sec. 4.2.3. We start again with a beam using a similar setup as for the straightbeam in Sec. 4.1.1, but we place an optically thick obstacle onthe beam path. “Thick obstacle” means that we set the absorptionopacity to a large value in the region occupied by the obstacle.From a physical point of view, the obstacle can then be viewed as ablack-body absorbing the incoming radiation without re-emitting it( 𝜂 = − . < 𝑥 < . , − . < 𝑦, 𝑧 < .
25 and covered with 150 × ×
50 uniformsized grid cells, we set 𝜅 𝑎 = within a sphere of radius 𝑅 = . 𝑥 -axis. All MNRAS000
50 uniformsized grid cells, we set 𝜅 𝑎 = within a sphere of radius 𝑅 = . 𝑥 -axis. All MNRAS000 , 000–000 (2020) attice-Boltzmann methods for radiative transport in computational astrophysics − . − . . . . y N pop = 50 p = 9 G a u ss i a n N pop ≈ N pop = 200 p = 17 N pop ≈ N pop = 392 p = 27 N pop ≈ − . − . . . . y N pop = 50 p = 11 L e b e d e v N pop = 194 p = 23 N pop = 434 p = 35 − .
25 0 .
00 0 . x − . − . . . . y N pop = 48 p = 9 Sph e r i c a l − .
25 0 .
00 0 . z − .
25 0 .
00 0 . x N pop = 192 p = 19 − .
25 0 .
00 0 . z − .
25 0 .
00 0 . x N pop = 393 p = 27 − .
25 0 .
00 0 . z . . . . . . E / E m a x Figure 5.
Spherical freely-streaming wave for different types of stencils. From top to bottom we show stencils derived from (i) Gaussian product quadrature,(ii) Lebedev quadrature and (iii) spherical design. Each column compares these three stencils for a comparable number of discrete velocities. Snapshots showthe radiation energy-density in the ( 𝑥, 𝑦 ) (left) and ( 𝑦, 𝑧 ) (right) plane after 200 iterations. In every panel we also report 𝑁 pop and the quadrature order 𝑝 . − .
25 0 .
00 0 . x − . . . y t/t max = 0 . − .
25 0 .
00 0 . x t/t max = 0 . − .
25 0 .
00 0 . x t/t max = 1 . . . . . . x . . . E / E m a x ( t = ) t/t max = 0 . t/t max = 0 . t/t max = 0 . t/t max = 0 . t/t max = 0 . t/t max = 1 . ∝ r − max( E ( t )) . . . . . . E/E max
Figure 6.
Top:
Radiation energy density (colour coded) and momentumdensity (red arrows) at three representative times in the (x,y) plane for theradiation wave test with the Lebedev quadrature of order 𝑝 = Bottom:
Diagonal profiles of the energy density at different times (blue to green) andinverse square law (red-dashed) fitted to the maxima of each time (red dots). other directions do not matter, since the radiation only propagatesparallel to the 𝑥 -axis in this test. Figure 7 shows the beam in the ( 𝑥, 𝑦 ) plane after it has propagated across the whole grid and passedthe obstacle marked as a cyan-coloured circle.It can be seen that, as expected, the beam is blocked by theobstacle and only a negligible amount of radiation diffuses insidethe optically thick region. This is due to the finite resolution of thegrid and is progressively suppressed as the grid is refined.We should note that, in contrast with all the tests performedso far, the set of equations that we solve here is very stiff because − . − . − . . . . . x − . . . y κ a = 10 | F | =1 − − − − − log ( E/E max ) Figure 7.
Radiation beam hitting an optically thick obstacle (cyan circle)with 𝜅 𝑎 = . The colour encodes the radiation energy-density and the redarrows show magnitude and direction of the energy momentum-density. of the high numerical value of 𝜅 𝑎 (as compared to the evolvedvariables). Obtaining a numerically stable solution with a reasonabletimestep (we here use again Δ 𝑡 = . Δ 𝑥 ) is then only possible whenusing an implicit solver in time. As described in Sec. 3.3, we haveimplemented the Lambda iteration -method. It is found that, at leastfor this simple test, it converges in at most three cycles at everytimestep, within a tolerance of Δ 𝐸 / 𝐸 = − . We next consider a test involving the emission of radiation as itoccurs for a background static fluid in thermodynamic equilibriumwith the radiation, that is, when the emissivity 𝜂 is equal to theabsorption opacity 𝜅 𝑎 . More specifically, we simulate a system inwhich 𝜂 = 𝜅 𝑎 = const. within a sphere of radius 𝑅 and zero outsidethe sphere. This is known as the homogeneous-/radiative-/emitting-sphere test, as first proposed by Smit et al. (1997). From a physicalpoint of view, this system can be thought of as a dense sphere MNRAS , 000–000 (2020) Weih et al. − − r/R . . . . . . E κ a = η = 10 /Rκ a = η = 10 /Rκ a = η = 1 /R LB analytic r/R M1 − − r/R . . . . . . F LB r/R M1 Figure 8.
Diagonal cuts for the radiation energy-density (left) and the magnitude of the momentum-density (right) of the equilibrium for the radiating-spheretest. We compare the numerical solution obtained with LB (crosses; left halves) and M1 (circles; right halves) with the analytic solution (solid lines.) for a low(green), moderate (red) and high (blue) value of 𝜅 𝑎 = 𝜂 . with a sharp boundary to vacuum, as is the case of a neutron star,that constantly emits radiation from its surface to the surroundingvacuum.The system eventually finds a steady state, for which the analyt-ical solution of the distribution function in terms of radial distance 𝑟 and azimuthal angle 𝜃 reads as follows: 𝑓 ( 𝑟, 𝜇 ) = 𝑏 ( − 𝑒 − 𝜅 𝑎 𝑠 ( 𝑟,𝜇 ) ) , (48)where 𝜇 (cid:66) cos 𝜃 , 𝑏 = 𝜅 𝑎 / 𝜂 = 𝑠 : = (cid:40) 𝑟 𝜇 + 𝑅𝑔 ( 𝑟, 𝜇 ) 𝑟 < 𝑅 and − < 𝜇 < , 𝑅𝑔 ( 𝑟, 𝜇 ) 𝑟 ≥ 𝑅 and √︁ − 𝑅 / 𝑟 < 𝜇 < , (49)with 𝑔 ( 𝑟, 𝜇 ) : = √︄ − 𝑟 𝑅 ( − 𝜇 ) . (50)Since the final equilibrium is spherically symmetric, the zerothmoment 𝐸 can be obtained via the integration of the distributionfunction over 𝜇 , i.e., 𝐸 ( 𝑟 ) = ∫ − 𝑑𝜇 𝑓 ( 𝑟, 𝜇 ) . (51)Likewise, we obtain the first moment 𝐹 as 𝐹 ( 𝑟 ) = ∫ − 𝑑𝜇 𝜇 𝑓 ( 𝑟, 𝜇 ) . (52)We run this test using a 3D grid with 128 uniform sizedgrid cells and a spherical-design stencil of order 20 ( 𝑁 pop = ) for three different values of 𝜅 𝑎 , i.e., 𝜅 𝑎 = 𝑅 − , 𝑅 − , 𝑅 − ,and where 𝑅 = /( 𝑛 𝑥 ) is the sphere’s radius.The discrete intensities are initialised as 𝐼 𝑖 ( 𝑟, 𝑡 = ) : = 𝑤 𝑖 (cid:40) 𝑟 < 𝑅 ,𝑟 − 𝑟 ≥ 𝑅 , (53)where 𝑤 𝑖 are the weights associated with the underlying velocitystencil.The results of these simulations are reported in Fig. 8 in termsof the radiation energy-density (left panel) and of the momentum-density (right panel). As already remarked by Radice et al. (2013) and Weih et al. (2020b),this test requires three dimensions on a Cartesian grid, since fluxes alsopropagate across grid-cells in the angular directions.
As one can appreciate, the LB method works very well in allof the three cases. As already discussed for the previous shadowtest, numerical stability in the case of 𝜅 = / 𝑅 is only possiblethanks to the implicit time-stepper.Figure 8 also offers a comparison with the corresponding re-sults obtained with an M1 scheme (these are shown in the rightportions of the two panels in Fig. 8). It is clear that for small valuesof 𝜅 𝑎 , the LB method performs significantly better than M1. Asalready reported by Weih et al. (2020b), the M1 code fails in thisspecific case due to the lack of the second moment, i.e., the correctpressure tensor. While the pressure tensor is exact in the limit ofinfinite optical depth (see red and blue curves in Fig. 8), it is not sofor the intermediate regime between optically thick and thin media.Indeed, the case of 𝜅 𝑎 = 𝑅 − (green curves) falls exactly in thisregime, for which the pressure tensor is interpolated inaccurately(see also Murchikova et al. (2017) for a detailed analysis of thistest for the M1 scheme with different closures). The failure of theM1 scheme in this test is particularly evident upon looking at theenergy density, which is systematically above the correct analyticalsolution inside the sphere (note that the tail, which is within the free-streaming regime is reproduced correctly). The same – albeit lessvisible – is true for the corresponding momentum density, which iseverywhere below the analytic solution. As a final test of the solution of the RTE with the LB method, weconsider the effects of scattering by simulating a Gaussian distri-bution that diffuses over a static background fluid. The same testhas been used also in the validation of various M1 codes (see,e.g., Pons et al. 2000; O’Connor 2015; Weih et al. 2020b). Theinitial conditions are given by: 𝐸 ( 𝑟, 𝑡 = ) = 𝐴 exp (cid:32) − ( 𝑟 − 𝑟 ) 𝜎 (cid:33) , (54)where 𝐴 is the energy-density amplitude and 𝜎 the width of theGaussian at 𝑡 = 𝑟 . By neglecting emission and ab-sorption ( 𝜂 = 𝜅 𝑎 = 𝐸 ( 𝑟, 𝑡 ) = 𝐴 𝜎 𝜎 + 𝜎 𝐷 exp (cid:32) − ( 𝑟 − 𝑟 ) ( 𝜎 + 𝜎 𝐷 ) (cid:33) , (55) MNRAS000
As one can appreciate, the LB method works very well in allof the three cases. As already discussed for the previous shadowtest, numerical stability in the case of 𝜅 = / 𝑅 is only possiblethanks to the implicit time-stepper.Figure 8 also offers a comparison with the corresponding re-sults obtained with an M1 scheme (these are shown in the rightportions of the two panels in Fig. 8). It is clear that for small valuesof 𝜅 𝑎 , the LB method performs significantly better than M1. Asalready reported by Weih et al. (2020b), the M1 code fails in thisspecific case due to the lack of the second moment, i.e., the correctpressure tensor. While the pressure tensor is exact in the limit ofinfinite optical depth (see red and blue curves in Fig. 8), it is not sofor the intermediate regime between optically thick and thin media.Indeed, the case of 𝜅 𝑎 = 𝑅 − (green curves) falls exactly in thisregime, for which the pressure tensor is interpolated inaccurately(see also Murchikova et al. (2017) for a detailed analysis of thistest for the M1 scheme with different closures). The failure of theM1 scheme in this test is particularly evident upon looking at theenergy density, which is systematically above the correct analyticalsolution inside the sphere (note that the tail, which is within the free-streaming regime is reproduced correctly). The same – albeit lessvisible – is true for the corresponding momentum density, which iseverywhere below the analytic solution. As a final test of the solution of the RTE with the LB method, weconsider the effects of scattering by simulating a Gaussian distri-bution that diffuses over a static background fluid. The same testhas been used also in the validation of various M1 codes (see,e.g., Pons et al. 2000; O’Connor 2015; Weih et al. 2020b). Theinitial conditions are given by: 𝐸 ( 𝑟, 𝑡 = ) = 𝐴 exp (cid:32) − ( 𝑟 − 𝑟 ) 𝜎 (cid:33) , (54)where 𝐴 is the energy-density amplitude and 𝜎 the width of theGaussian at 𝑡 = 𝑟 . By neglecting emission and ab-sorption ( 𝜂 = 𝜅 𝑎 = 𝐸 ( 𝑟, 𝑡 ) = 𝐴 𝜎 𝜎 + 𝜎 𝐷 exp (cid:32) − ( 𝑟 − 𝑟 ) ( 𝜎 + 𝜎 𝐷 ) (cid:33) , (55) MNRAS000 , 000–000 (2020) attice-Boltzmann methods for radiative transport in computational astrophysics [h!] . . . . r/r max . . . . . . E / E m a x ( t = ) κ ∆ x = 1 analyticLB . . . . r/r max κ ∆ x = 10 t = 0 t = 20 t = 40 t = 100 t = 200 t = 300 t = 500 t = 1000 Figure 9.
Gaussian wave packet diffusing in a static fluid. The analyticsolution (solid lines) is compared with results of simulations (dots) obtainedon a two-dimensional grid of size 100 , with Δ 𝑡 = .
1, and 𝜆 = .
5. Thetwo panels show results at low and high Peclet numbers. where 𝜎 𝐷 = √ 𝐷𝑡 . The connection between the microscopic pa-rameters 𝜅 and 𝜅 in Eq. (17) and the diffusion coefficient 𝐷 isdiscussed in Appendix A by an asymptotic analysis.We perform two different simulations at different values of the“Peclet number” (Pe : = 𝜅 Δ 𝑥 ), tracking the history of the system atdifferent time intervals. In Fig. 9, we compare the analytical solution(lines) with the numerical results (dots) obtained running in 2D, on a100 grid, with 𝜅 Δ 𝑥 = 𝜅 Δ 𝑥 = (right panel).In both cases, we set 𝜆 (cid:66) 𝜅 / 𝜅 = .
5. The same test can beperformed in three dimensions, leading to results of similarly goodquality. Finally, in Fig. 10 we show that, in stark contrast with theresults shown for the free-streaming regime (Fig. 5), the evolutionis well captured even by quadratures of comparatively low orders,featuring a small number of discrete directions. This highlights thefact that the LB method performs extremely well in the diffusionlimit, in fact the one it was born for. In general, we have observedthat in this regime quadratures with order 𝑝 = After having shown that the proposed radiative LB method performswell in all the regimes of the solution of the RTE on a static fluid,we next move on to test the method for a dynamical fluid, i.e.one coupled to radiation via an RMHD simulation of a realisticastrophysical scenario. As a result, the test presented here is ofgreat importance as it allows us to explore our novel approachunder conditions typical of the astrophysical scenarios for which ithas been developed in the first place.We should also note that, strictly speaking, this is not a genuinetest, as the explored scenario does not have an analytic solution tobe compared with. Moreover, simulations of this type have beenperformed before only by Rivera-Paleo & Guzmán (2019) usingan M1 scheme for evolving the radiation but with rather differentprescriptions for the properties of the radiation field.Hence, to gain confidence on the reliability of our results andcontrast them with those obtained when the RTE is solved usingan M1 scheme, we carry out additional simulations of the identicalphysical scenario but making use of the M1 code
FRAC (Weih et al.2020b). While this does not prove the correctness of our results –both codes could be incorrect despite the many tests passed – it does provide us with the confidence necessary to implement the LBmethod for even more realistic astrophysical scenarios.More specifically, we simulate the evolution of a relativisticjet as it propagates through the interstellar medium. Simulations ofthis type have a long history since such jets are of major importancein the study of active galactic nuclei (see Perucho 2019, for a recentreview), where they are produced through the accretion processonto a supermassive black hole (see Porth et al. 2019, for a recentcomparative study).Highly energetic relativistic jets are known to accompany thephenomenology of short gamma-ray bursts and are associated withthe merger of two magnetised neutron stars (Rezzolla et al. 2011).Here, we simulate this problem by coupling our LB code tothe
Black hole accretion code ( BHAC ) (Porth et al. 2017).
BHAC is afinite-volume code that solves the equations of general-relativisticMHD in a fixed and curved spacetime. The results presented here,however, are restricted to a flat background spacetime.In essence, we perform the coupling between
BHAC and the LBcode following the strategy indicated below.1. At every iteration, we pass the fluid rest-mass density 𝜌 , tem-perature 𝑇 and three-velocity 𝑣 𝑖 to the LB code, from which we thencompute the emissivity and opacities (see also below).2. While BHAC advances the conservative RMHD variables intime, the LB code does the same for the 𝑁 pop populations of theradiation specific intensity, 𝐼 𝑖 .3. After performing the streaming and collision step, the LB codecomputes the zeroth ( 𝐸 ), first ( 𝐹 𝑖 ) and second moment ( 𝑃 𝑖 𝑗 ) of theradiation distribution function according to Eqs. (27) and (28).4. From these moments, we compute the radiative source terms S = 𝑊 ( ˜ 𝜅 𝑎 𝐽 − ˜ 𝜂 ) + ˜ 𝜅𝐻 and S 𝑗 = 𝑊 ( ˜ 𝜅 𝑎 𝐽 − ˜ 𝜂 ) 𝑣 𝑗 + ˜ 𝜅𝐻 𝑗 , which wereturn to BHAC . Here 𝑊 is the Lorentz-factor, ˜ 𝜅 (cid:66) ˜ 𝜅 𝑎 + ˜ 𝜅 𝑠 (cid:66) ˜ 𝜅 𝑎 +( ˜ 𝜅 − / 𝜅 ) , and 𝐽 and 𝐻 𝑗 are the radiation energy and momentumdensity in the comoving fluid frame, to which we transform via 𝐽 = 𝑊 (cid:0) 𝐸 − 𝐹 𝑖 𝑣 𝑖 + 𝑃 𝑖 𝑗 𝑣 𝑖 𝑣 𝑗 (cid:1) (56) 𝐻 𝑗 = 𝑊 (cid:0) 𝐹 𝑖 𝑣 𝑖 − 𝐸 (cid:1) 𝑣 𝑗 + 𝑊 ℎ 𝑗𝑖 𝐹 𝑖 − 𝑊 ℎ 𝑗𝑖 𝑣 𝑘 𝑃 𝑖𝑘 , (57)where ℎ 𝑖 𝑗 = 𝑊 𝑣 𝑖 𝑣 𝑗 + 𝛿 𝑖 𝑗 is the projection operator orthogonal tothe fluid velocity. 𝐻 can then be computed from 𝐻 𝜇 𝑢 𝜇 = BHAC has updated its variables and before cycling to1), we add the sources S and S 𝑗 to the energy and momentumequation, respectively, which are solved in conservative form within BHAC .As it is the case in most simulations of high-energy astrophys-ical phenomena, the fluid moves at relativistic speeds. Eq. (17), iswritten in an Eulerian (lab) frame so that the opacities and emissivi-ties 𝜂 , 𝜅 𝑎 , 𝜅 , 𝜅 it employs are to be evaluated in the same Eulerianframe. However, the microphysics used to derive such quantities iswell defined only in the fluid’s rest-frame (i.e., the frame co-movingwith the fluid) where the corresponding quantities ˜ 𝜂 , ˜ 𝜅 𝑎 , ˜ 𝜅 , ˜ 𝜅 areisotropic and can be written in a compact way. Therefore, care needsto be taken in transforming the opacities and emissivities betweenthe two frames, as we discuss in detail in Appendix B.For the setup of our simulation, we follow Martí et al. (1997),who have extensively analysed relativistic jets in a purely hydrody-namical context. The jet is simply injected through a circular nozzleat the lower edge of the computational domain and propagating par-allel to the coordinate 𝑧 -axis in a Cartesian grid. The simulation isthen characterized by four parameters: the Newtonian Mach num-ber M : = 𝑣 jet / 𝑐 𝑠 , with 𝑐 𝑠 the local sound speed, the jet Lorentzfactor 𝑊 : = ( − 𝑣 ) − / with 𝑣 jet the jet propagation velocity, MNRAS , 000–000 (2020) Weih et al. − . − . . . . y N pop = 8 p = 3 G a u ss i a n N pop = 18 p = 5 N pop = 72 p = 11 − . − . . . . y N pop = 6 p = 3 L e b e d e v N pop = 14 p = 5 N pop = 50 p = 11 − . . . x − . − . . . . y N pop = 8 p = 3 Sph e r i c a l − . . . z − . . . x N pop = 18 p = 5 − . . . z − . . . x N pop = 72 p = 11 − . . . z . . . . . . E / E m a x Figure 10.
Scattering wave in 3D. The simulations are performed on a grid of size 128 , with Δ 𝑡 = . 𝜅 Δ 𝑥 = 𝜆 = − .
25. From top to bottom we showthe results obtained using stencils derived from (i) Gaussian product quadrature, (ii) Lebedev quadrature and (iii) spherical design. Each column comparesthese three stencils for a comparable number of discrete velocities. Snapshots show the radiation energy-density in the ( 𝑥, 𝑦 ) (left) and ( 𝑦, 𝑧 ) (right) planeafter 100 / Δ 𝑡 iterations. Contour lines are used to compare the analytical solution (continuous red lines) with the numerical results (white dotted lines). In eachpanel, we also report 𝑁 pop and the quadrature order 𝑝 . the ratio of the jet rest-mass density to that of the ambient medium R : = 𝜌 jet / 𝜌 amb and the pressure ratio K : = 𝑝 jet / 𝑝 amb . We simulatea pressure matched jet, i.e., K =
1, with R = . 𝑊 = M =
42 on a grid of size − . < 𝑥, 𝑦 < . , < 𝑧 <
60 covered by160 × ×
640 grid cells, where the jet is injected at 𝑧 = 𝑟 jet = 𝑁 pop =
154 discrete velocity directions.We initialise the populations at zero and let the radiation evolveself-consistently during the simulation. To this purpose, we set˜ 𝜅 𝑎 = 𝜌 𝑇 − . and ˜ 𝜂 = 𝜎 SB / 𝜋 ˜ 𝜅 𝑎 𝑇 , where the fluid temperature 𝑇 is computed from an ideal gas equation of state. This absorp-tion opacity is motivated by the Rosseland mean opacity for ther-mal bremsstrahlung (Rybicki & Lightman 1986) and the emissivitysimply follows from Kirchhoff’s law upon assuming black-bodyradiation.We should note that in contrast to the pure-RMHD simulation,which is scale invariant, the coupled RMHD-radiation simulationfixes the length scale via the value of the Stefan-Boltzmann constant 𝜎 SB and the density assumed for the ambient medium 𝜌 amb . Here,we simply use 𝜌 amb = 𝜎 SB = .
1, which does not lead toa physically realistic setup, but ensures that a moderate amountof radiation is produced, that neither dominates the fluid nor isdominated by it.Finally, we also add scattering using ˜ 𝜅 = − 𝜌 and ˜ 𝜅 =
0. This choice is motivated by the microphysical process of Thomsonscattering, which is proportional to the number of scatters in themedium (hence, the choice for ˜ 𝜅 ) and has no preferred direction(hence, the choice for ˜ 𝜅 ).We compare the results of the pure-RMHD and the RTE-coupled simulations in Fig. 11, whose left panel refers to the pure-RMHD jet, the central panel to the RTE-coupled solution obtainedwith the LB method, and the right panel to the corresponding evo-lution when the RTE-coupled solution is obtained with the M1scheme. A quick comparison of the pure-RMHD jet morphologyin the left panel shows that it is in good agreement with the onepresented by Martí et al. (1997) and Aloy et al. (1999), where thistype of jets has been studied extensively. On the other hand, thesolutions employing a coupling with the radiation are considerablydifferent. In particular, the RTE-coupled simulation with the LBmethod shows that the jet propagates more slowly, i.e., the Lorentzfactor is ∼
15% smaller than for the pure-RMHD case. This is due tothe fact that the fluid making up the jet loses energy via the emissionof radiation. By stark contrast, the RTE-coupled solution obtainedwith the M1 scheme shows that the jet propagates more rapidly.While this behaviour is similar to the one reported by Rivera-Paleo& Guzmán (2019) – who, in addition, continuously injected energyin the radiation field – we believe it is actually incorrect.The origin of this substantially different dynamical behaviouris due to the poor handling by the M1 scheme of radiation interactingwith itself. We find this to cause significantly different distributionsof the energy-density. In the case of the M1-evolution, the radiationenergy density accumulates mostly in a narrow region along the 𝑧 -axis, while it is more evenly spread in the case of LB-evolution. Theactual cause of this radiation focussing is the same as at the originof the poor performance of the M1 scheme in the beam-crossingtest reported in Fig. 4. Also in this case, in fact, different beams ofradiation originating from the recollimation shock produced near the MNRAS000
15% smaller than for the pure-RMHD case. This is due tothe fact that the fluid making up the jet loses energy via the emissionof radiation. By stark contrast, the RTE-coupled solution obtainedwith the M1 scheme shows that the jet propagates more rapidly.While this behaviour is similar to the one reported by Rivera-Paleo& Guzmán (2019) – who, in addition, continuously injected energyin the radiation field – we believe it is actually incorrect.The origin of this substantially different dynamical behaviouris due to the poor handling by the M1 scheme of radiation interactingwith itself. We find this to cause significantly different distributionsof the energy-density. In the case of the M1-evolution, the radiationenergy density accumulates mostly in a narrow region along the 𝑧 -axis, while it is more evenly spread in the case of LB-evolution. Theactual cause of this radiation focussing is the same as at the originof the poor performance of the M1 scheme in the beam-crossingtest reported in Fig. 4. Also in this case, in fact, different beams ofradiation originating from the recollimation shock produced near the MNRAS000 , 000–000 (2020) attice-Boltzmann methods for radiative transport in computational astrophysics − x [ r jet ] z [ r j e t ] pure hydro − x [ r jet ] hydro + LB − x [ r jet ] hydro + M1 − . − . − . − . − . − . . . l og ( ρ / ρ a m b ) Figure 11.
Cut through the ( 𝑥, 𝑧 ) plane for the relativistic jet after 𝑡 = 𝑟 jet . Shown is the rest-mass density for the pure-hydro (left) and thecoupled hydro-radiation using LB (middle) and M1 (right). injection region of the jet (Mizuno et al. 2015), intersect along the 𝑧 -axis and lead to an inconsistent combination of radiation fluxes.This is illustrated in the top panels of Fig. 12, which reportsa cut through the ( 𝑥, 𝑧 ) plane of the radiation energy density 𝐸 atan early time during the jet evolution, namely at 𝑡 = 𝑟 jet , usingeither the LB method (left) or the M1 scheme (right). Note thatin both cases there is a triangular region of small 𝐸 originatingfrom the recollimation shock, and ultimately due to the contactdiscontinuity between propagating jet and ambient medium (seeAloy & Rezzolla 2006, for a discussion of the role of the contactdiscontinuity in accelerating the jet). This region of small 𝐸 issurrounded on both sides by high- 𝐸 streams, which appear as whitein the colorcode used. In the case of the M1-evolution (right panel),these two streams merge at the top of the triangular low- 𝐸 regionleading to the same pathology discussed in Sec. 4.1.1 (cf., Fig. 4)and to the artificial accumulation of radiation along the 𝑧 -axis. Thelatter then provides additional momentum, pushing the jet forwardand overcompensating the linear momentum lost in the productionof the radiation. In the LB-evolution, on the other hand, the twobeams do not merge but cross correctly. As a result, the radiationenergy density is not artificially focused and 𝐸 spreads out over alarger region, leading to a broader bow shock ahead of the jet, whichis also propagating more slowly.Note that these problems in the M1-evolution affect the dy-namics of the jet only downstream of the first recollimation shock.This is very clearly illustrated in the bottom panels of Fig. 12, whichrefer instead to the rest-mass densities in the two cases and to thesame time as the top panels, i.e., 𝑡 = 𝑟 jet . As one would expect,the differences are in this case very small since the radiative effectshave not yet been able to play a role, which they will instead do for 𝑡 (cid:38) 𝑟 jet . z [ r j e t ] hydro + LB hydro + M1 − . . . x [ r jet ] z [ r j e t ] − . . . x [ r jet ] − − − − − − l og ( E / E m a x ) − . − . − . − . − . − . . l og ( ρ / ρ a m b ) Figure 12.
Top panel:
Cut through the ( 𝑥, 𝑧 ) plane of the radiation energydensity at 𝑡 = 𝑟 jet for the coupled RMHD-radiation simulations using theLB method (left) and the M1 scheme (right). Bottom panel:
The same as inthe top panel, but for the rest-mass density.
In summary, these results go well-beyond our intention of pro-viding proof-of-principle evidence for a correct implementation ofthe LB method in a fully coupled relativistic-fluid configuration. Inparticular, they clearly show the ability of the LB method to han-dle correctly scenarios with physical conditions that are very closeto those encountered in relativistic astrophysical phenomena. Moreimportantly, they highlight that under these very same conditions,the M1 approach commonly employed – even by us (Weih et al.2020a) – may suffer from artefacts that may affect the dynamics ofrelativistic jets, such as those expected to be generated in a binarysystem of merging neutron stars, and whose accurate description isessential to gain insight in the launching of relativistic jets in shortgamma-ray bursts.
MNRAS , 000–000 (2020) Weih et al.
Novel computational methods are expected to i) feature good scal-ability on parallel computers, ii) perform better than existing onesin accuracy, efficiency or both.In what follows we show that the LB method proposed herepossesses both these qualities, i.e. parallel scalability and efficiency.Indeed, one key ingredient in the success and widespread adop-tion of LB is parallel efficiency. Thanks to the synchronous algorith-mic flow stemming from the stream-collide paradigm, LB schemesexhibit high amenability to parallel coding, making them naturaltargets for efficient implementations on modern computer architec-tures (Godenschwager et al. 2013; Bernaschi et al. 2009; Mazzeo &Coveney 2008; Bernaschi et al. 2010; Succi et al. 2019).In this section, we provide a brief performance evaluation of thenumerical scheme designed in this paper. With respect to classicalLB models, our scheme presents two main differences: i) we employoff-lattice quadratures, which require interpolation to implement thestreaming phase; ii) the number of velocity components forming thestencil is significantly larger with respect to standard LB schemes.We have implemented our numerical scheme using directive-based programming environments, such as OpenMP and
OpenACC ,to expose parallelism (Calore et al. 2016a). The advantage of thisapproach is that the code is portable and can therefore be compiledand executed on diverse architectures, from commodity CPU-basedprocessors to GPU accelerators. We test the code on an Intel Sky-lake 20-cores processor, a commonly adopted architecture in theHPC-market. We measure performances in terms of Million LatticeUpdates per Second (MLUPS):MLUPS = 𝐿 𝑁 iter 𝑡 exe , (58)where 𝑡 exe measures the execution time (in seconds) required tosimulate 𝑁 iter timesteps on a grid of 𝐿 points. In essence, for aproblem of fixed size and iterations 𝐿 and 𝑁 iter , a parallel imple-mentation should lead to an execution time that decreases linearlywith the number of cores and hence a linearly growing MLUPS.This is shown in Fig. 13, where we assess the scalability of the codeon a single Intel Skylake board. We solve the emitting-sphere bench-mark with an explicit Euler stepper, on an 𝐿 = grid, and witha spherical-design quadrature of order 𝑝 =
20 with 𝑁 pop =
200 dis-crete components. The code scales up to 20 threads with a parallelefficiency above 70%.These figures are already quite good, but could be further im-proved. A limiting performance factor is given by inefficient mem-ory accesses, which in turn leads to a sub-optimal use of the cacheand of the vector unit of the processor. Omitting details, we simplymention here that a careful optimisation of the data-layout used tostore the grid in memory should be taken into consideration. Thetwo most common data layouts used in several stencil applicationsare the so-called array of structures (AoS) and structure of arrays(SoA) schemes; in the AoS layout, all populations associated to onelattice site are stored in contiguous memory locations. Conversely,in the SoA scheme all the populations having the same index 𝑖 arestored contiguously, while populations belonging to the same latticesite are stored far from each other at non unit-stride addresses. Re-cently, more sophisticated data-layouts have been designed explic-itly for LB applications (Shet et al. 2013) yielding strong benefits,especially when targeting multi-core architectures with wide vectorunits and large cache memories (Calore et al. 2019). The applicationof these optimisation and a more in-depth analysis will be reportedelsewhere.In Fig. 14, we show the performance obtained on a dual- Number of Threads M L U P S . . . . . . P a r a ll e l e ffi c i e n c y Figure 13.
Analysis of performance scalability on a single Intel Skylakeboard. The results refer to simulations of the emitting-sphere benchmarkwith an explicit Euler stepper, on an 𝐿 = grid, using a stencil with 𝑁 pop =
200 components. The figures reported correspond to the best out of10 trials for each set of parameters.
Skylake processor as a function of 𝑁 pop . We show results for anexplicit and also for the implicit stepper described in Sec. 3.3, whichextends the stability of the method also to the case of very stiff terms.To simplify the comparison, we have executed the Lambda iteration loop exactly 5 times at each grid-cell. Since the execution time ofthe current implementation is completely bounded by memory ac-cesses, the difference between the implicit and the explicit schemeis almost negligible.The figure also presents a comparison with the performanceof the special-relativistic version of the M1 code
FRAC (Weih et al.2020b), featuring a similar level of optimisation. For M1 the majorbottleneck is the closure relation, which requires a complex root-finding algorithm for obtaining the pressure tensor. In principle, theperformance of
FRAC can be arbitrarily good or bad depending onthe desired accuracy of the root-finding step. Here we consider astandard setup (i.e., a relative accuracy of 10 − for root-finding andcomputing the closure on cell-centres as well as cell-faces). As canbe seen from the dashed line in Fig. 14, the LB method outperformsthe M1 code for 𝑁 pop (cid:46) − ) and lower (10 − )relative accuracy for the root-finding, which is shown as a blue bandin Fig. 14.Finally, we comment that the LB method is particularly wellsuited for GPU implementations (Bernaschi et al. 2010; Calore et al.2016b). For this reason we also test a GPU-optimized version of ourcode on an NVIDIA V100 GPU. We observe the performances tobe systematically one order of magnitude higher than those reportedon the CPU, both for stencils with 10 −
100 components, suitablefor simulations in which scattering terms are relevant, as well as forstencils with >
200 components, which instead allow to extend theapplicability of the solver to free-streaming regimes. The presentversion of the code allows for efficient memory accesses on GPUsarchitectures, in turn exposing the higher costs of the implicit solver,which for large values of 𝑁 pop is found to be 2 − . × slower thanthe explicit Euler method. MNRAS000
200 components, which instead allow to extend theapplicability of the solver to free-streaming regimes. The presentversion of the code allows for efficient memory accesses on GPUsarchitectures, in turn exposing the higher costs of the implicit solver,which for large values of 𝑁 pop is found to be 2 − . × slower thanthe explicit Euler method. MNRAS000 , 000–000 (2020) attice-Boltzmann methods for radiative transport in computational astrophysics N pop M L U P S M1 LB − Explicit stepperLB − Implicit stepper
Figure 14.
Comparison of the performance (measured MLUPS) achieved ona dual-Skylake processor (40 threads) as a function of 𝑁 pop , the number ofthe discrete components forming the stencil. The results refer to simulationsof the emitting-sphere benchmark, on an 𝐿 = grid. Dots correspond tosimulations based on an explicit Euler stepper, while squares correspond tosimulations using an implicit stepper. The dashed line shows the performanceof the M1 code FRAC with the blue-shaded band corresponding to differentaccuracies for the M1 closure. The figures reported are the best out of 10trials for each set of parameters.
We have presented an extension of the LB method, which is com-monly used in classical fluid dynamics, to the solution of the RTEin special relativity, thus making it applicable to three-dimensionalsimulations of high-energy astrophysical phenomena.After implementing the new method in flat spacetime and underthe grey approximation, we analysed its performance in a numberof code tests. In this way, we have shown that while the LB methodperforms extremely well in the diffusion limit, the free-streamingregime represents the most difficult one to treat accurately. In thisregime, a large number of discrete velocity directions is requiredsince the propagation of radiation is restricted to these directionsonly. Nevertheless, even in this optically thin limit, the LB methodproves superior to the commonly used M1 method, as witnessed bythe fact that radiation beams cross correctly.This feature is important for astrophysical systems with anaccretion disk and torus, from which radiation is expected to crossand lead to the phenomenology observed in short gamma-ray bursts.Besides this important advantage, the new LB method also outper-forms the M1 method in the intermediate regime, between diffusionand free-streaming, mostly because it does not need to rely on aclosure relation to compute higher moments, as it is the case for theM1 scheme, where the pressure tensor is simply interpolated via theclosure.The accuracy in the calculation of the moments also dependson the underlying quadrature, for which we have compared threepossibilities. While in our code tests only minor differences couldbe seen among these quadratures, it is mathematically clear that theLebedev and spherical-design quadratures are more accurate thanthe Gauss-Legendre quadrature commonly used in direct Boltzmannsolvers (see e.g., Nagakura et al. 2018).We have also coupled our new radiation code to the GRMHD-code
BHAC and simulated a relativistic jet, where the radiation back-reacts dynamically onto the fluid. Not only does this represent thefirst such simulation using an LB scheme, but also proves thatour new method is indeed applicable to high-energy astrophysics. Furthermore, when comparing the corresponding results obtainedwith the more standard M1 code
FRAC , that employs a moment-based scheme, we have shown that the LB-method solution doesnot suffer from the inaccuracies that plague the M1 method.Finally, we have also shown that the LB method is faster thanthe M1 method for a number of discrete directions 𝑁 pop (cid:46) 𝑁 pop ≈ −
800 might be necessary. Inthis case the LB method is slightly more expensive than the M1scheme. Considering, however, the high amenability of LB to GPUimplementations, a major speed-up can certainly be achieved alongthis line.While this paper is meant to provide a presentation of the LBmethod for the solution of the radiative-transfer equation in com-putational astrophysics, a number of improvements are possible,both in terms of astrophysical applications, and in terms of mathe-matical and numerical developments. The former involves a moredetailed and realistic investigation of the role played by radiationin impacting the dynamics and imaging of astrophysical relativisticjets. The latter necessarily involves the extension of the method to ageneral-relativistic framework, which requires the extension of thestreaming-step to curved spacetimes.Finally, the numerical developments will include the possibilityof using numerical grids with various forms of static and dynamicmesh-refinements. Also in this case, an adjustment of the streaming-step will be needed to allow for the streaming from coarse to finegrid cells and vice-versa. Work along these lines is in progress.
ACKNOWLEDGMENTS
LRW acknowledges support from HGS-HIRe. AG acknowledgesfunding by "Contributo 5 per mille assegnato all’Università degliStudi di Ferrara-dichiarazione dei redditi dell’anno 2017". DS hasbeen supported by the European Union’s Horizon 2020 research andinnovation programme under the Marie Sklodowska-Curie grantagreement No. 765048. SS acknowledges funding from the Euro-pean Research Council under the European Union’s Horizon 2020framework programme (No. P/2014-2020)/ERC Grant AgreementNo. 739964 (COPMAT). He also wishes to thank A. Ferrara, P.Mocz, D. Spergel and J. Stone for illuminating discussions. Sup-port also comes in part from “PHAROS”, COST Action CA16214;LOEWE-Program in HIC for FAIR; the ERC Synergy Grant “Black-HoleCam: Imaging the Event Horizon of Black Holes” (Grant No.610058). The simulations were performed on the SuperMUC andSuperMUC-NG clusters at the LRZ in Garching, on the LOEWEcluster in CSC in Frankfurt, on the HazelHen cluster at the HLRSin Stuttgart, and on the COKA computing cluster at Università diFerrara.
DATA AVAILABILITY
The data underlying this article will be shared on reasonable requestto the corresponding author.
MNRAS , 000–000 (2020) Weih et al.
REFERENCES
Ahrens C., Beylkin G., 2009, Proceedings of the Royal Society A: Mathe-matical, Physical and Engineering Sciences, 465, 3103Aloy M. A., Rezzolla L., 2006, Astrophys. J., 640, L115Aloy M. A., Ibáñez J. M. ., Martí J. M. ., Gómez J.-L., Müller E., 1999,Astrophys. J., 523, L125Ambruş V. E., Blaga R., 2018, Phys. Rev. C, 98, 035201Asinari P., Mishra S., Borchiellini R., 2010, Numerical Heat Transfer, PartB: Fundamentals, 57, 126Atkinson K., Han W., 2012, Spherical Harmonics and Approxima-tions on the Unit Sphere: An Introduction. Springer Berlin Hei-delberg, doi:10.1007/978-3-642-25983-8, http://dx.doi.org/10.1007/978-3-642-25983-8
Baiotti L., Rezzolla L., 2017, Rept. Prog. Phys., 80, 096901Beentjes C. H. L., 2015, International Conference on Mathematics andComputational Methods Applied to Nuclear Science and Engineering,pp 1–15Benzi R., Succi S., Vergassola M., 1992, Phys. Rep., 222, 145Bernaschi M., Melchionna S., Succi S., Fyta M., Kaxiras E., Sircar J., 2009,Computer Physics Communications, 180, 1495Bernaschi M., Fatica M., Melchionna S., Succi S., Kaxiras E., 2010, Con-currency and Computation: Practice and Experience, 22, 1Bhatnagar P. L., Gross E. P., Krook M., 1954, Phys. Rev., 94, 511Bindra H., Patil D. V., 2012, Phys. Rev. E, 86, 016706Bovard L., Martin D., Guercilena F., Arcones A., Rezzolla L., Korobkin O.,2017, Phys. Rev. D, 96, 124005Broadwell J. E., 1964, The Physics of Fluids, 7, 1243Bruenn S. W., 1975, in Bergman P. G., Fenyves E. J., Motz L., eds, TexasSymposium on Relativistic Astrophysics Vol. 262, Seventh Texas Sym-posium on Relativistic Astrophysics. pp 80–94, doi:10.1111/j.1749-6632.1975.tb31422.xBruenn S. W., 1985, Astrophys. J., Supp., 58, 771Calore E., Gabbana A., Kraus J., Schifano S. F., Tripiccione R., 2016a,Concurrency and Computation: Practice and Experience, 28, 3485Calore E., Gabbana A., Kraus J., Pellegrini E., Schifano S. F., TripiccioneR., 2016b, Parallel Computing, 58, 1Calore E., Gabbana A., Schifano S., Tripiccione R., 2019, The InternationalJournal of High Performance Computing Applications, 33, 124Cardall C. Y., Endeve E., Mezzacappa A., 2013, Phys. Rev. D, 87, 103004Chapman S., Cowling T. G., 1970, The Mathematical Theoryof Non-Uniform Gases, 3rd ed. Cambridge University Press,doi:10.1119/1.1942035Coelho R. C., Mendoza M., Doria M. M., Herrmann H. J., 2018, Computers& Fluids, 172, 318Colgate S. A., White R. H., 1966, Astrophys. J., 143, 626Delsarte P., Goethals J. M., Seidel J. J., 1977, Geometriae Dedicata, 6, 363Dietrich T., Ujevic M., 2017, Classical and Quantum Gravity, 34, 105014Dünweg B., Ladd A. J. C., 2009, Lattice Boltzmann Simulations of SoftMatter Systems. Springer Berlin Heidelberg, Berlin, Heidelberg, pp89–166, doi:10.1007/978-3-540-87706-6_2Event Horizon Telescope Collaboration et al., 2019, Astrophys. J. Lett., 875,L1Fernández R., Tchekhovskoy A., Quataert E., Foucart F., Kasen D., 2019,Mon. Not. R. Astron. Soc., 482, 3373Foucart F., 2018, Mon. Not. R. Astron. Soc., 475, 4186Foucart F., et al., 2015, Phys. Rev. D, 91, 124021Fragile P. C., Gillespie A., Monahan T., Rodriguez M., Anninos P., 2012,Astrophys. J., Supp., 201, 9Fragile P. C., Olejar A., Anninos P., 2014, Astrophys. J., 796, 22Frisch U., Hasslacher B., Pomeau Y., 1986, Phys. Rev. Lett., 56, 1505Fromm C. M., Perucho M., Porth O., Younsi Z., Ros E., Mizuno Y., ZensusJ. A., Rezzolla L., 2018, Astron. Astrophys., 609, A80Fujibayashi S., Kiuchi K., Nishimura N., Sekiguchi Y., Shibata M., 2018,Astrophys. J., 860, 64Gabbana A., Mendoza M., Succi S., Tripiccione R., 2018, Computers &Fluids, 172, 644 Gabbana A., Simeoni D., Succi S., Tripiccione R., 2020, Physics Reports,863, 1Gairola A., Bindra H., 2017, Annals of Nuclear Energy, 99, 151Galeazzi F., Kastaun W., Rezzolla L., Font J. A., 2013, Phys. Rev. D, 88,064009Gamba I. M., Haack J. R., Hauck C. D., Hu J., 2017, SIAM Journal onScientific Computing, 39, B658Godenschwager C., Schornbaum F., Bauer M., Köstler H., Rüde U., 2013, inProceedings of the International Conference on High Performance Com-puting, Networking, Storage and Analysis. SC ’13. Association for Com-puting Machinery, New York, NY, USA, doi:10.1145/2503210.2503273Grad H., 1949a, Communications on Pure and Applied Mathematics, 2, 325Grad H., 1949b, Communications on Pure and Applied Mathematics, 2, 331Gross B., Atzberger P., 2018, Journal of Computational Physics, 371, 663Hardy J., Pomeau Y., de Pazzis O., 1973, Phys. Rev. Lett., 31, 276He X., Chen S., Doolen G. D., 1998, Journal of Computational Physics, 146,282Higuera F. J., Jiménez J., 1989, EPL (Europhysics Letters), 9, 663Higuera F. J., Succi S., 1989, Europhysics Letters (EPL), 8, 517Higuera F. J., Succi S., Benzi R., 1989, EPL (Europhysics Letters), 9, 345Hildebrand F. B., 1956, Introduction to numerical analysis. New York :McGraw-Hill, https://trove.nla.gov.au/work/3478395
Janka H.-T., Langanke K., Marek A., Martínez-Pinedo G., Müller B., 2007,Physics Reports, 442, 38Just O., Obergaulinger M., Janka H. T., 2015, Mon. Not. R. Astron. Soc.,453, 3386Karlin I. V., Sichau D., Chikatamarla S. S., 2013, Phys. Rev. E, 88, 063310Krüger T., Kusumaatmaja H., Kuzmin A., Shardt O., Silva G., ViggenE. M., 2017, The Lattice Boltzmann Method. Springer InternationalPublishing, doi:10.1007/978-3-319-44649-3, http://dx.doi.org/10.1007/978-3-319-44649-3
Kuroda T., Takiwaki T., Kotake K., 2016, Astrophys. J., Supp., 222, 20Lebedev V., 1975, USSR Computational Mathematics and MathematicalPhysics, 15, 44Lebedev V., 1976, USSR Computational Mathematics and MathematicalPhysics, 16, 10Lebedev V. I., 1977, Siberian Mathematical Journal, 18, 99Levermore C. D., Pomraning G. C., 1981, Astrophys. J., 248, 321Livio M., Buchler J. R., Colgate S. A., 1980, Astrophys. J. Lett., 238, L139Lutsko J. F., Lam J., 2018, Phys. Rev. E, 98, 012604Martí J. M., Müller E., Font J. A., Ibáñez J. M., Marquina A., 1997, Astro-phys. J., 479, 151Massaioli F., Benzi R., Succi S., 1993, Europhysics Letters (EPL), 21, 305Mazzeo M., Coveney P., 2008, Computer Physics Communications, 178,894McCulloch R., Bindra H., 2016, Computers & Fluids, 124, 261McHardy C., Horneber T., Rauh C., 2016, Optics Express, 24, 16999McKinney J. C., Tchekhovskoy A., Sadowski A., Narayan R., 2014, Mon.Not. R. Astron. Soc., 441, 3177McNamara G. R., Zanetti G., 1988, Phys. Rev. Lett., 61, 2332Melon Fuksman J. D., Mignone A., 2019, Astrophys. J., Supp., 242, 20Mendoza M., Karlin I., Succi S., Herrmann H. J., 2013, Phys. Rev. D, 87,065027Mezzacappa A., Liebendörfer M., Messer O. E. B., Hix W. R., ThilemannF.-K., Bruenn S. W., 2001, Phys. Rev. Lett., 86, 1935Mihalas D., Auer L. H., 2001, Journal of Quantitative Spectroscopy andRadiative Transfer, 71, 61Miller J. M., Ryan B. R., Dolence J. C., 2019, Astrophys. J., Supp., 241, 30Mink A., McHardy C., Bressel L., Rauh C., Krause M. J., 2020, Journal ofQuantitative Spectroscopy and Radiative Transfer, 243, 106810Mishra S. C., Poonia H., Vernekar R. R., Das A. K., 2014, Heat TransferEngineering, 35, 1267Mizuno Y., Gómez J. L., Nishikawa K.-I., Meli A., Hardee P. E., RezzollaL., 2015, Astrophys. J., 809, 38Most E. R., Papenfort L. J., Rezzolla L., 2019, Mon. Not. R. Astron. Soc.,490, 3588Murchikova E. M., Abdikamalov E., Urbatsch T., 2017, Mon. Not. R. Astron.Soc., 469, 1725 MNRAS000
Kuroda T., Takiwaki T., Kotake K., 2016, Astrophys. J., Supp., 222, 20Lebedev V., 1975, USSR Computational Mathematics and MathematicalPhysics, 15, 44Lebedev V., 1976, USSR Computational Mathematics and MathematicalPhysics, 16, 10Lebedev V. I., 1977, Siberian Mathematical Journal, 18, 99Levermore C. D., Pomraning G. C., 1981, Astrophys. J., 248, 321Livio M., Buchler J. R., Colgate S. A., 1980, Astrophys. J. Lett., 238, L139Lutsko J. F., Lam J., 2018, Phys. Rev. E, 98, 012604Martí J. M., Müller E., Font J. A., Ibáñez J. M., Marquina A., 1997, Astro-phys. J., 479, 151Massaioli F., Benzi R., Succi S., 1993, Europhysics Letters (EPL), 21, 305Mazzeo M., Coveney P., 2008, Computer Physics Communications, 178,894McCulloch R., Bindra H., 2016, Computers & Fluids, 124, 261McHardy C., Horneber T., Rauh C., 2016, Optics Express, 24, 16999McKinney J. C., Tchekhovskoy A., Sadowski A., Narayan R., 2014, Mon.Not. R. Astron. Soc., 441, 3177McNamara G. R., Zanetti G., 1988, Phys. Rev. Lett., 61, 2332Melon Fuksman J. D., Mignone A., 2019, Astrophys. J., Supp., 242, 20Mendoza M., Karlin I., Succi S., Herrmann H. J., 2013, Phys. Rev. D, 87,065027Mezzacappa A., Liebendörfer M., Messer O. E. B., Hix W. R., ThilemannF.-K., Bruenn S. W., 2001, Phys. Rev. Lett., 86, 1935Mihalas D., Auer L. H., 2001, Journal of Quantitative Spectroscopy andRadiative Transfer, 71, 61Miller J. M., Ryan B. R., Dolence J. C., 2019, Astrophys. J., Supp., 241, 30Mink A., McHardy C., Bressel L., Rauh C., Krause M. J., 2020, Journal ofQuantitative Spectroscopy and Radiative Transfer, 243, 106810Mishra S. C., Poonia H., Vernekar R. R., Das A. K., 2014, Heat TransferEngineering, 35, 1267Mizuno Y., Gómez J. L., Nishikawa K.-I., Meli A., Hardee P. E., RezzollaL., 2015, Astrophys. J., 809, 38Most E. R., Papenfort L. J., Rezzolla L., 2019, Mon. Not. R. Astron. Soc.,490, 3588Murchikova E. M., Abdikamalov E., Urbatsch T., 2017, Mon. Not. R. Astron.Soc., 469, 1725 MNRAS000 , 000–000 (2020) attice-Boltzmann methods for radiative transport in computational astrophysics Nagakura H., et al., 2018, ApJ, 854, 136O’Connor E., 2015, Astrophys. J., Supp., 219, 24Paardekooper J. P., Kruip C. J. H., Icke V., 2010, A&A, 515, A79Pareschi L., Russo G., 2005, Journal of Scientific Computing, 25, 129Paschalidis V., 2017, Classical and Quantum Gravity, 34, 084002Peng Y., Shu C., Chew Y. T., 2003, Phys. Rev. E, 68, 026701Perego A., Rosswog S., Cabezón R. M., Korobkin O., Käppeli R., ArconesA., Liebendörfer M., 2014, Mon. Not. R. Astron. Soc., 443, 3134Perego A., Radice D., Bernuzzi S., 2017, Astrophys. J. Lett., 850, L37Perucho M., 2019, Galaxies, 7, 70Philippi P. C., Hegele L. A., dos Santos L. O. E., Surmas R., 2006, Phys.Rev. E, 73, 056702Pomraning G. C., 1981, J. Quant. Spectrosc. Radiative Transfer, 26, 385Pons J. A., Ibáñez J. M., Miralles J. A., 2000, Mon. Not. R. Astron. Soc.,317, 550Porth O., Olivares H., Mizuno Y., Younsi Z., Rezzolla L., Moscibrodzka M.,Falcke H., Kramer M., 2017, Computational Astrophysics and Cosmol-ogy, 4, 1Porth O., Chatterjee K., Event Horizon Telescope Collaboration 2019, As-trophys. J. Supp., 243, 26Radice D., Abdikamalov E., Rezzolla L., Ott C. D., 2013, Journal of Com-putational Physics, 242, 648Radice D., Galeazzi F., Lippuner J., Roberts L. F., Ott C. D., Rezzolla L.,2016, Mon. Not. R. Astron. Soc., 460, 3255Rahman N., Just O., Janka H. T., 2019, Mon. Not. R. Astron. Soc., 490,3545Rampp M., 1997, PhD thesis, -Rezzolla L., Miller J. C., 1994, Class. Quantum Grav., 11, 1815Rezzolla L., Zanotti O., 2013, Relativistic Hydro-dynamics. Oxford University Press, Oxford, UK,doi:10.1093/acprof:oso/9780198528906.001.0001Rezzolla L., Giacomazzo B., Baiotti L., Granot J., Kouveliotou C., AloyM. A., 2011, Astrophys. J. Letters, 732, L6Rivera-Paleo F. J., Guzmán F. S., 2019, Astrophys. J., Supp., 241, 28Rivet J.-P., Boon J. P., 2001, Lattice Gas Hydrodynamics. Cam-bridge Nonlinear Science Series, Cambridge University Press,doi:10.1017/CBO9780511524707Roedig C., Zanotti O., Alic D., 2012, Mon. Not. R. Astron. Soc., 426, 1613Rosswog S., Liebendörfer M., 2003, Mon. Not. R. Astron. Soc., 342, 673Rosswog S., Korobkin O., Arcones A., Thielemann F.-K., Piran T., 2014,Mon. Not. R. Astron. Soc., 439, 744Ruffert M., Janka H.-T., Schaefer G., 1996, Astron. Astrophys., 311, 532Rybicki G. B., Lightman A. P., 1986, Radiative Processes in Astrophysics.Wiley-VCHSa¸dowski A., Narayan R., Tchekhovskoy A., Zhu Y., 2013, Mon. Not. R.Astron. Soc., 429, 3533Shan X., 2016, Journal of Computational Science, 17, 475Shan X., He X., 1998, Phys. Rev. Lett., 80, 65Shan X., Yuan X.-F., Chen H., 2006, Journal of Fluid Mechanics, 550, 413Shet A. G., Sorathiya S. H., Krithivasan S., Deshpande A. M., Kaul B.,Sherlekar S. D., Ansumali S., 2013, Phys. Rev. E, 88, 013314Shibata M., Kiuchi K., Sekiguchi Y., Suwa Y., 2011, Progress of TheoreticalPhysics, 125, 1255Siegel D. M., Ciolfi R., 2016, Astrophys. J., 819, 14Siegel D. M., Metzger B. D., 2017, Physical Review Letters, 119, 231102Skinner M. A., Dolence J. C., Burrows A., Radice D., Vartanyan D., 2019,Astrophys. J., Supp., 241, 7Smit J. M., Cernohorsky J., Dullemond C. P., 1997, Astron. Astrophys., 325,203Sobolev S., 1962, Soviet Mathematics Doklady, 3, 1307Stepán, Jirí Bestard, Jaume Jaume Trujillo Bueno, Javier 2020, A&A, 636,A24Succi S., 2015, EPL (Europhysics Letters), 109, 50001Succi S., 2018, The Lattice Boltzmann Equation: ForComplex States of Flowing Matter. OUP Oxford,doi:10.1093/oso/9780199592357.001.0001Succi S., Amati G., Bernaschi M., Falcucci G., Lauricella M., MontessoriA., 2019, Computers & Fluids, 181, 107 Thorne K. S., 1981, Mon. Not. R. Astron. Soc., 194, 439Vernekar R. R., Mishra S. C., 2014, International Journal of Heat and MassTransfer, 77, 218Wang Y., Ma Y., Xie M., 2019, Progress in Nuclear Energy, 110, 341Weih L. R., Hanauske M., Rezzolla L., 2020a, Phys. Rev. Lett., 124, 171103Weih L. R., Olivares H., Rezzolla L., 2020b, Mon. Not. R. Astron. Soc.,495, 2285Womersley R. S., 2018, in , Contemporary Computational Mathematics -A Celebration of the 80th Birthday of Ian Sloan. Springer InternationalPublishing, pp 1243–1285, doi:10.1007/978-3-319-72456-0_57Yi H.-L., Yao F.-J., Tan H.-P., 2016, Phys. Rev. E, 94, 023312Zanotti O., Roedig C., Rezzolla L., Del Zanna L., 2011, Mon. Not. R. Astron.Soc., 417, 2899d ' Humières D., Lallemand P., Frisch U., 1986, Europhysics Letters (EPL),2, 291
APPENDIX A: LINK BETWEEN MICROSCOPIC ANDMACROSCOPIC PARAMETERS
In this section of the appendix we perform an asymptotic analysisto link the microscopic parameters with the macroscopic ones. Inparticular we consider the specific limit of zero emission and ab-sorption ( 𝜅 𝑎 = 𝜂 = 𝜕 𝑡 𝐸 = 𝐷 Δ 𝐸 . (A1)In Appendix A1 the link between the diffusive coefficient 𝐷 and the scattering parameters ( 𝜅 , 𝜅 ) is established analyticallythrough a Chapman-Enskog expansion (Chapman & Cowling 1970).In Appendix A2 the analytic expressions are extended to accountfor extra numerical corrections by fitting results from numericalsimulations. A1 Chapman-Enskog Analysis
We start from Eq. (29) assuming zero absorption and emission, i.e., 𝐼 𝑖 ( 𝒓 + 𝑐 ˆ 𝒏 𝑖 Δ 𝑡, 𝑡 + Δ 𝑡 ) − 𝐼 𝑖 ( 𝒓 , 𝑡 ) = − 𝑐𝜅 Δ 𝑡 (cid:16) 𝐼 𝑖 ( 𝒓 , 𝑡 ) − 𝐼 eq 𝑖 ( 𝒓 , 𝑡 ) (cid:17) . (A2)with 𝐼 eq 𝑖 ( 𝒓 , 𝑡 ) given by Eq. (30). Note that throughout this section wewill write 𝑐 explicitly and write vector components with Greek in-dices rather than using boldface vectors. We also adopt the Einsteinconvention of summing over repeated indices.Taking a Taylor expansion of the left-hand side of Eq. (A2),and including terms up to the second order, gives: 𝐼 𝑖 ( 𝒙 + 𝒏 𝑖 Δ 𝑡, 𝑡 + Δ 𝑡 ) − 𝐼 𝑖 ( 𝒙 , 𝑡 ) = Δ 𝑡 (cid:0) 𝜕 𝑡 + 𝑐𝑛 𝛼𝑖 𝜕 𝛼 (cid:1) 𝐼 𝑖 + Δ 𝑡 (cid:0) 𝜕 𝑡 + 𝑐𝑛 𝛼𝑖 𝜕 𝛼 (cid:1) 𝐼 𝑖 + O( Δ 𝑡 ) . (A3)We also expand the differential operator with respect to time 𝜕 𝑡 = 𝜖 𝜕 ( ) 𝑡 + 𝜖 𝜕 ( ) 𝑡 + O( 𝜖 ) , (A4)and space: 𝜕 𝛼 = 𝜖 𝜕 ( ) 𝛼 + O( 𝜖 ) , (A5)where 𝜖 (cid:28) 𝐼 𝑖 = 𝐼 ( ) 𝑖 + 𝜖 𝐼 ( ) 𝑖 + 𝜖 𝐼 ( ) 𝑖 + O( 𝜖 ) , (A6)where 𝐼 ( ) 𝑖 ≡ 𝐼 eq 𝑖 . We also recall the definition of the first and secondmoment of the distribution: 𝐸 = ∑︁ 𝑖 𝐼 𝑖 , (A7) MNRAS , 000–000 (2020) Weih et al. 𝐹 𝛼 = ∑︁ 𝑖 𝑛 𝛼𝑖 𝐼 𝑖 . (A8)Assuming the most basic level of isotropy for the stencil used in thenumerical method we have ∑︁ 𝑖 𝑤 𝑖 = , ∑︁ 𝑖 𝑤 𝑖 𝑛 𝛼𝑖 = , ∑︁ 𝑖 𝑤 𝑖 𝑛 𝛼𝑖 𝑛 𝛽𝑖 = 𝑑 𝛿 𝛼𝛽 , ∑︁ 𝑖 𝑤 𝑖 𝑛 𝛼𝑖 𝑛 𝛽𝑖 𝑛 𝛾𝑖 = , (A9)where 𝛼, 𝛽 , and 𝛾 run over the spatial indexes in 𝑑 dimensions.By integrating Eq. (30), in combination with Eq. (A9), weget the following definitions for the moments of the equilibriumdistribution: ∑︁ 𝑖 𝐼 eq 𝑖 = 𝐸 , (A10) ∑︁ 𝑖 𝐼 eq 𝑖 𝑛 𝛽𝑖 = 𝜆𝑑 𝐹 𝛽 , (A11) ∑︁ 𝑖 𝐼 eq 𝑖 𝑛 𝛽𝑖 𝑛 𝛾𝑖 = 𝑑 𝛿 𝛽𝛾 𝐸 . (A12)Since we are neglecting absorption and emission and consider thediffusion limit, where the radiation is in thermodynamic equilibriumwith the underlying fluid, we find due to conservation ∫ ( 𝐼 − 𝐼 eq ) 𝑑 Ω = ⇒ ∫ 𝐼𝑑 Ω = ∫ 𝐼 eq 𝑑 Ω , (A13)which in the discretized form leads to the condition ∑︁ 𝑖 𝐼 𝑖 = ∑︁ 𝑖 𝐼 eq 𝑖 = 𝐸 , ∑︁ 𝑖 𝐼 ( 𝑘 ) 𝑖 = ∀ 𝑘 ≥ . (A14)We do not write down any condition for the first moment, since itspreservation depends on the choice of 𝜆 .Replacing the left-hand side of Eq. (A2) with its second orderTaylor expansion, i.e., Eq. (A3), we get Δ 𝑡 (cid:0) 𝜕 𝑡 + 𝑐𝑛 𝛼𝑖 𝜕 𝛼 (cid:1) 𝐼 𝑖 ( 𝒙 , 𝑡 ) + Δ 𝑡 (cid:0) 𝜕 𝑡 + 𝑐𝑛 𝛼𝑖 𝜕 𝛼 (cid:1) 𝐼 𝑖 ( 𝒙 , 𝑡 ) = − 𝑐 Δ 𝑡𝜅 (cid:16) 𝐼 𝑖 ( 𝒙 , 𝑡 ) − 𝐼 eq 𝑖 ( 𝒙 , 𝑡 ) (cid:17) . (A15)We now plug in the above Eqs. (A4), (A5), and (A6) and performa multi-scale expansion in which we keep track separately of termsup to order 𝜖 and 𝜖 . The resulting equations are O( 𝜖 ) : (cid:16) 𝜕 ( ) 𝑡 + 𝑐𝑛 𝛼𝑖 𝜕 ( ) 𝛼 (cid:17) 𝐼 ( ) 𝑖 ≈ − 𝑐𝜅 𝐼 ( ) 𝑖 , (A16) O( 𝜖 ) : (cid:18) − Δ 𝑡 𝑐𝜅 (cid:19) (cid:16) 𝜕 ( ) 𝑡 + 𝑐𝑛 𝛼𝑖 𝜕 ( ) 𝛼 (cid:17) 𝐼 ( ) 𝑖 + 𝜕 ( ) 𝑡 𝐼 ( ) 𝑖 ≈ − 𝑐𝜅 𝐼 ( ) 𝑖 . (A17)We then take into consideration Eq. (A16) and integrate (i.e.sum over all 𝑁 pop populations), getting ∑︁ 𝑖 (cid:16) 𝜕 ( ) 𝑡 𝐼 ( ) 𝑖 + 𝑐𝑛 𝛼𝑖 𝜕 ( ) 𝛼 𝐼 ( ) 𝑖 (cid:17) ≈ − 𝑐𝜅 ∑︁ 𝑖 𝐼 ( ) 𝑖 , (A18)which using Eq. (A14) and the definition of the first order momentleads to 𝜕 ( ) 𝑡 𝐸 + 𝑐 𝜆𝑑 𝜕 ( ) 𝛼 𝐹 𝛼 = . (A19)Next, starting again from Eq. (A16), we multiply by 𝑛 𝛽𝑖 andintegrate, which yields 𝜆𝑑 𝜕 ( ) 𝑡 𝐹 𝛽 + 𝑑 𝑐𝜕 𝛽 𝐸 = − 𝑐𝜅 ∑︁ 𝑖 𝑛 𝛽𝑖 𝐼 ( ) 𝑖 . (A20) Note that the RHS of the above equation vanishes when the conser-vation of the first moment is ensured; for the moment we leave it ina general form and evaluate this term later on.We now integrate Eq. (A17) and obtain ∑︁ 𝑖 (cid:18) − Δ 𝑡 𝑐𝜅 (cid:19) (cid:16) 𝜕 ( ) 𝑡 + 𝑐𝑛 𝛼𝑖 𝜕 ( ) 𝛼 (cid:17) 𝐼 ( ) 𝑖 + ∑︁ 𝑖 𝜕 ( ) 𝑡 𝐼 ( ) 𝑖 ≈ − 𝜅 ∑︁ 𝑖 𝐼 ( ) 𝑖 , (A21) (cid:18) − Δ 𝑡 𝑐𝜅 (cid:19) 𝑐𝜕 ( ) 𝛼 ∑︁ 𝑖 𝑛 𝛼𝑖 𝐼 ( ) 𝑖 + 𝜕 ( ) 𝑡 𝐸 = , (A22) 𝜕 ( ) 𝑡 𝐸 = − (cid:18) − Δ 𝑡 𝑐𝜅 (cid:19) 𝑐𝜕 ( ) 𝛼 ∑︁ 𝑖 𝑛 𝛼𝑖 𝐼 ( ) 𝑖 . (A23)The RHS can be derived from Eq. (A20), leading to 𝜕 ( ) 𝑡 𝐸 = 𝑐𝜅 (cid:18) − Δ 𝑡 𝑐𝜅 (cid:19) 𝑐𝜕 ( ) 𝛼 (cid:18) 𝜆𝑑 𝜕 ( ) 𝑡 𝐹 𝛼 + 𝑑 𝑐𝜕 𝛼 𝐸 (cid:19) , (A24)which can be re-arranged as 𝜕 ( ) 𝑡 𝐸 = 𝐷 (cid:18) Δ ( ) 𝐸 + 𝜆 𝑐 𝜕 ( ) 𝛼 𝜕 ( ) 𝑡 𝐹 𝛼 (cid:19) (A25)with 𝐷 = 𝑐 𝑑𝑐𝜅 (cid:18) − Δ 𝑡 𝑐𝜅 (cid:19) . (A26)We note that Eq. (A25) is the diffusion equation with 𝐷 its diffusioncoefficient plus an extra error term. This error term is the same onethat arises for the equations of the M1 system, when considering theoptically thick limit and neglecting absorption and emission. It isremarkable that despite the M1 system being a set of macroscopicequations, where the connection between 𝐷 and the system’s scat-tering coefficient is directly evident, in contrast to the microscopicequation of our LB method, we arrive at the exact same macroscopicequation for the diffusion limit. A2 Numerical Fit
In the previous section we have established a link between the mi-croscopic parameters and the diffusion coefficient by performingthe Chapman-Enskog expansion. In principle one would need toextend the expansion to include corrections coming from higherorder terms. Besides, one should also account for extra dissipa-tive effects coming from the fact that the LB method described inthe present work is not based on a space-filling Cartesian latticeand consequently requires interpolation. Since the overall analysiswould become rather tedious from an analytical point of view, inthis section we numerically evaluate the corrections to Eq. (A26)that so far have been neglected.We start by taking into consideration the effect of varying theparameter 𝜆 . We consider the same numerical setup discussed inSec. 4.2.3, working in 2D, on a 100 ×
100 grid, with Δ 𝑡 = Δ 𝑥 and 𝜅 Δ 𝑥 =
1. In Fig. A1 we show the results of numerical simulationsfor a few selected values of 𝜆 ; the results clearly show how 𝜆 impactsthe diffusion speed. Since Eq. (A26) does not depend on 𝜆 , but onlyon 𝜅 , we extend it by assuming a dependency on the parameter 𝜒 = 𝜅 (cid:18) + 𝛼 𝑑 𝜆 (cid:19) , (A27)i.e., a linear combination of 𝜅 and 𝜅 , with 𝛼 a coefficient left tobe determined. MNRAS000
1. In Fig. A1 we show the results of numerical simulationsfor a few selected values of 𝜆 ; the results clearly show how 𝜆 impactsthe diffusion speed. Since Eq. (A26) does not depend on 𝜆 , but onlyon 𝜅 , we extend it by assuming a dependency on the parameter 𝜒 = 𝜅 (cid:18) + 𝛼 𝑑 𝜆 (cid:19) , (A27)i.e., a linear combination of 𝜅 and 𝜅 , with 𝛼 a coefficient left tobe determined. MNRAS000 , 000–000 (2020) attice-Boltzmann methods for radiative transport in computational astrophysics − . − . . . . r/r max . . . . . E / E m a x ( t = ) λ = − λ = 0 λ = t/ ∆ t = 50 t/ ∆ t = 150 t/ ∆ t = 300 Figure A1.
Effect of the parameter 𝜆 on the diffusion speed of a Gaussianpulse. The simulations are performed in a bidimensional grid of size 𝐿 = Δ 𝑡 = 𝑘 Δ 𝑥 = In Sec. 4.1 we have discussed the artificial diffusivity intro-duced by the interpolation scheme used to implement the propa-gation step. In order to try to capture these effects we propose thefollowing extension for Eq. (A26): 𝐷 = 𝑐𝑑𝜒 (cid:20) − (cid:18) + 𝛼 (cid:19) Δ 𝑡𝑐𝜒 (cid:21) + Δ 𝑥 Δ 𝑡 𝛼 . (A28)In the above, 𝛼 introduces a correction to the leading term of theTaylor expansion of the propagation step, which is also present inother off-grid LB schemes (see e.g. Coelho et al. (2018)). This cor-rection becomes vanishingly small as one takes smaller timesteps.We attribute this correction to the specific interpolation schemeemployed. Besides, we also introduce an extra coefficient 𝛼 whichintroduces a background diffusivity which depends on the specificstencil taken into consideration.Our task consists now in determining 𝛼 , 𝛼 , 𝛼 . To this aimwe perform several numerical simulations, in which we once againreproduce the benchmark described in Sec. 4.2.3, varying 𝜅 Δ 𝑥 , Δ 𝑡 and 𝜆 . We follow the evolution of each simulation up to 𝑡 = Δ 𝑡 . Next, we estimate the numerical diffusion coefficient ofeach simulation. To do so we compare the numerical results withthe analytic solution (Eq. (55) ), and fit a numerical value of thediffusion coefficient which leaves the L2-norm of the relative errorbelow 1% over several finite time-steps ( 𝑡 / Δ 𝑡 = , , , . . . 𝜆 = . 𝛼 ≈ − in both 2 and 3 dimensions. Theparameters 𝛼 and 𝛼 , however, depend both on the dimensionalityand on the specific stencil taken into consideration, although bothtend to stabilise when considering stencils formed by a sufficientlylarge number of components. To give an example, conducting theanalysis in 2D with 𝑁 pop =
120 we get 𝛼 ≈ . 𝛼 ≈ .
65. In3D, instead, using a spherical-design quadrature with order 𝑝 = 𝑁 pop = 𝛼 ≈ . 𝛼 ≈ . .
50 0 .
75 1 .
00 1 .
25 1 .
50 1 .
75 2 . κ ∆ x . . . . . . . D ∆ t / ∆ x ∆ t = 0 . t = 0 . t = 0 . t = 0 . t = 1 . Figure A2.
Example of the numerical analysis performed to fit the expressionfor the diffusion coefficient as in Eq. (A28). The results are shown forthe specific case 𝜆 = .
25. Crosses represent estimates of the numericaldiffusion coefficient which ensures that the L2-norm of the relative erroris within 1% with respect to the analytic solution. Dashed lines show thepredicted diffusion coefficient using Eq. (A28) with 𝛼 = / 𝛼 = . 𝛼 = . APPENDIX B: EMISSIVITY AND OPACITIES IN THELAB FRAME
For radiative transfer simulations on a moving fluid background it issimpler to express the microscopic quantities describing emission,absorption and scattering processes in the comoving fluid-frame,in which the fluid is at rest. Since our LB scheme is designed inthe lab frame, it is then necessary to transform these fluid-framequantities to their lab-frame counterparts. For doing so we followthe derivation in Mihalas & Auer (2001), which for completenesswe summarise here for the grey approximation.For a fluid moving with three-velocity 𝑣 𝑖 and Lorentz factor 𝑊 ,the frequency 𝜈 in the lab frame of a radiation particle propagatingin direction ˆ 𝑛 𝑖 is transformed to the comoving fluid frame via˜ 𝜈 = 𝑊𝜈 ( − 𝑣 𝑖 ˆ 𝑛 𝑖 ) . (B1)It was then shown that the frequency-dependent quantities, whichare isotropic in the fluid frame, transform like 𝜂 𝜈 = 𝜈 ˜ 𝜈 ˜ 𝜂 ˜ 𝜈 = ˜ 𝜂 ˜ 𝜈 𝑊 ( − 𝑣 𝑖 ˆ 𝑛 𝑖 ) (B2) 𝜅 𝑎,𝜈 = ˜ 𝜈𝜈 ˜ 𝜅 𝑎, ˜ 𝜈 = 𝑊 ( − 𝑣 𝑖 ˆ 𝑛 𝑖 ) ˜ 𝜅 𝑎, ˜ 𝜈 . (B3)The frequency-averaged absorption opacity defined by Eq. (21) thensimply follows as 𝜅 𝑎 = ∫ ∞ 𝑊 ( − 𝑣 𝑖 ˆ 𝑛 𝑖 ) ˜ 𝜅 𝑎, ˜ 𝜈 𝐼 𝜈 𝑑𝜈 ∫ ∞ 𝐼 𝜈 𝑑𝜈 = 𝑊 ( − 𝑣 𝑖 ˆ 𝑛 𝑖 ) ˜ 𝜅 𝑎 , (B4)where ˜ 𝜅 𝑎 is the frequency-integrated fluid-frame absorption opacity.The transformation of the frequency-integrated emissivity is 𝜂 = ∫ ∞ 𝜈 𝜂 𝜈 𝑑𝜈 = ∫ ∞ ˜ 𝜈 ˜ 𝜂 ˜ 𝜈 𝑑𝜈𝑊 ( − 𝑣 𝑖 ˆ 𝑛 𝑖 ) = ∫ ∞ ˜ 𝜈 ˜ 𝜂 ˜ 𝜈 𝑑 ˜ 𝜈𝑊 ( − 𝑣 𝑖 ˆ 𝑛 𝑖 ) (B5) = ˜ 𝜂𝑊 ( − 𝑣 𝑖 ˆ 𝑛 𝑖 ) , where ˜ 𝜂 is the frequency-integrated fluid-frame emissivity. MNRAS , 000–000 (2020) Weih et al.
The scattering opacities are more complicated to transformand we here only consider iso-energetic isotropic scattering likethe Thomson scattering process, which we use in Sec. 5. We thenrecognize that for 𝜅 we have two terms on the RHS of Eq. (17). Theone proportional to 𝐼 acts like the absorption term and transformscorrespondingly, i.e., 𝜅 𝐼 = 𝑊 ( − 𝑣 𝑖 ˆ 𝑛 𝑖 ) ˜ 𝜅 𝐼 . (B6)The second term proportional to 𝐸 acts like an emission term, forwhich we use the same transformation that we used in Eq. (B5). Wethen find 𝜅 𝐸 = ˜ 𝜅 𝐽𝑊 ( − 𝑣 𝑖 ˆ 𝑛 𝑖 ) , (B7)where 𝐽 is the radiation energy density in the fluid frame, which canbe computed from the lab-frame moments according to Eq. (56).Taking all of the above transformations together means that forthe simulation of the relativistic jet in Sec. 5, we solve the mixed-frame equation1 𝑐 𝜕𝐼𝜕𝑡 + ˆ 𝒏 · ∇ 𝐼 = − 𝑊 ( − 𝑣 𝑖 ˆ 𝑛 𝑖 )( ˜ 𝜅 𝑎 + ˜ 𝜅 ) 𝐼 + ˜ 𝜂 + ˜ 𝜅 𝐽𝑊 ( − 𝑣 𝑖 ˆ 𝑛 𝑖 ) (B8)instead Eq. (17). MNRAS000