Observing supermassive black holes in virtual reality

We present a full 360 degree (i.e., 4$\pi$ steradian) general-relativistic ray-tracing and radiative transfer calculations of accreting supermassive black holes. We perform state-of-the-art three-dimensional general relativistic magnetohydrodynamical simulations using the BHAC code, subsequently post-processing this data with the radiative transfer code RAPTOR. All relativistic and general-relativistic effects, such as Doppler boosting and gravitational redshift, as well as geometrical effects due to the local gravitational field and the observer's changing position and state of motion, are therefore calculated self-consistently. Synthetic images at four astronomically-relevant observing frequencies are generated from the perspective of an observer with a full 360-degree view inside the accretion flow, who is advected with the flow as it evolves. As an example, we calculated images based on recent best-fit models of observations of Sagittarius A*. These images are combined to generate a complete 360-degree Virtual Reality movie of the surrounding environment of the black hole and its event horizon. Our approach also enables the calculation of the local luminosity received at a given fluid element in the accretion flow, providing important applications in, e.g., radiation feedback calculations onto black hole accretion flows. In addition to scientific applications, the 360-degree Virtual Reality movies we present also represent a new medium through which to communicate black hole physics to a wider audience, serving as a powerful educational tool.


Introduction
Active Galactic Nuclei (AGN) are strong sources of electromagnetic radiation from the radio up to γ-rays. Their source properties can be explained in terms of a galaxy hosting an accreting supermassive black hole (SMBH) in its core. The Milky Way also harbours a candidate SMBH, Sagittarius A* (Sgr A*), which is subject to intensive Very-Long-Baseline Interferometric (VLBI) studies [1][2][3][4][5][6]. Sgr A* is one of the primary targets of the Event Horizon Telescope Collaboration (EHTC), which aims to image for the very first time the "shadow" of a black hole [7]. Theoretical calculations predict this shadow to manifest as a darkening of the inner accretion flow image anticipated to be observed due to the presence of a black hole event horizon, representing the region within which no radiation can escape [7][8][9]. The apparent size on the sky of this shadow is constrained by Einstein's General Theory of Relativity (GR) [9][10][11][12][13][14][15][16], and observational measurements of the black hole shadow size and shape can in principle provide a strong test of the validity of GR in the strong-field regime [7,9,15,17].
The theoretical aspects of the observational study of Sgr A* require the generation of general-relativistic magnetohydrodynamical (GRMHD) simulation data of the accretion flow onto a black hole, which is subsequently used to calculate synthetic observational data for physically-motivated plasma models which can be compared to actual observational data. In the past, synthetic observational data was generated by ray-tracing radiative transfer codes which calculate the emission originating from the accreting black hole and measured by a far away observer by solving the equations of radiative transfer along geodesics, i.e., the paths of photons (or particles) as they propagate around the black hole in either static spacetimes [e.g. [18][19][20][21][22][23][24][25][26][27][28][29][30] or dynamical spacetimes [31,32].
These models vary only in the dynamics of the black hole accretion flow, with the observer remaining stationary through the calculations. In this work, we consider the most general case of an observer who can vary arbitrarily in both their position (with respect to the black hole) and their state of motion. In particular, the observer is chosen to follow the flow of the accreting plasma in a physically-meaningful manner through advection, and therefore all dynamical effects introduced by the motion of the observer around the black hole are also correctly included in the imaging calculation.
With recent developments in Graphical Processor Units (GPUs) and Virtual Reality (VR) rendering, it has become possible to visualise these astrophysical objects at high resolutions in a 360 • (i.e., 4π steradian) format that covers the entire celestial sphere of an observer, enabling the study of the surroundings of an accreting black hole from within the accretion flow itself. Virtual Reality is a broad concept that encompasses different techniques, such as immersive visualisation, stereographic rendering, and interactive visualisations. In this work, we explore the first of these three, by rendering the full celestial sphere of the observer along a trajectory. The viewer can then look in any direction during the animation; this is also known as 360 • VR. Another important feature of VR, stereographic rendering, presents different images to each eye, so that the viewer experiences stereoscopic depth. For our application, however, this technique is not relevant, since the physical distance between the eyes of the observer is much smaller than the typical length scale of a supermassive black hole (which is 6.645 × 10 11 cm for Sagittarius A*), and therefore we would not see any depth in the image (just as we do not see stereoscopic depth when looking at the Moon). Interactive visualisations, where the viewer also has the freedom to change his or her position, would require real-time rendering of the environment, which is beyond the reach of current computational resources.
Our new way of visualising black holes enables the study of accretion from the point of view of an observer close to the black hole event horizon, with the freedom to image in all directions, as opposed to the perspective of an observer far away from the source with a fixed position and narrow field of view. In the case of a distant observer, the source appears projected onto the celestial sphere (thus appearing two-dimensional). Since one cannot easily distinguish three-dimensional structures within the accretion flow, placing the observer inside the flow itself opens a new window in understanding the geometrical structure and dynamical properties of such systems. Several researchers have previously considered an observer moving around, or falling into a black hole, e.g., (1) falling through the event horizon as illustrated through the gravitational lensing distortions of different regions (e.g., the ergo-region and event horizon), represented as chequerboard patterns projected onto an observer's image plane [33], (2) a flight through a simulation of a non-rotating black hole [34], (3) a flight through an accretion disk of a black hole using an observer with a narrow field of view camera [35], (4) a 360 • VR movie of an observer falling into a black hole surrounded by vacuum with illumination provided exclusively by background starlight, i.e., without an accretion flow [36], (5) a 360 • VR movie of a hotspot orbiting a SMBH [37], and (6) a 360 • VR movie of an N-body/hydrodynamical simulation of the central parsec of the Galactic center [38]. In this study, we consider a self-consistent three-dimensional GRMHD simulation of the accretion flow onto a spinning (Kerr) black hole, determining its time evolution and what an observer would see in full 360 • VR as they move through the dynamically evolving flow. To image accreting black holes in VR, we use the generalrelativistic radiative-transfer (GRRT) code RAPTOR [30]. The code incorporates all important general-relativistic effects, such as Doppler boosting and gravitational lensing in curved spacetimes, and can be compiled and run on both Central Processing units (CPU's) and GPU's by using NVIDIA's OpenACC framework.
In this work, we investigate the environment of accreting black holes from within the accretion flow itself with a virtual camera. As an example astrophysical case we model the supermassive black hole Sgr A*, although the methods presented in this work are generally applicable to any black hole as long as the radiation field's feedback onto the accreting plasma has a negligible effect on the plasma's magnetohydrodynamical properties, which is the case for Low Luminosity AGNs or low/hard state X-ray binaries.
The trajectory of this camera consists of two phases: a hovering trajectory, where the observer moves with a pre-defined velocity, and a particle trajectory, where the observer's instantaneous velocity is given by a trajectory of a tracer particle computed with a seperate axisymetric GRMHD simulation. The tracer particle follows the local plasma velocity (specifically, it is obtained by interpolating the plasma velocity of the GRMHD simulation cells to the camera's location) . We present a 360 • VR simulation of Sgr A*, demonstrating the applications of VR for studying not just accreting black holes but also for education, public outreach and data visualisation and interpretation amongst the wider scientific community. In section 2 we describe the camera setup, present several black hole shadow lensing tests, describe the camera trajectories and outline the radiative transfer calculation. In section 3 we present our 360 • VR movie of an accreting black hole. In section 4 we discuss our results and outlook.

Methods
In this section, we introduce the virtual camera setup, present black hole shadow vacuum lensing tests using both stationary and free-falling observers at different radial positions, discuss the different camera trajectories used in the VR movie shown later in this article and introduce the GRMHD plasma model that is used as an input for the geometry of the accretion flow onto the black hole.
The original RAPTOR code [30] initialises rays (i.e., photon geodesics) using impact parameters determined form coordinate locations on the observer's image plane [39]. This method is not suitable for VR since it only applies to distant observers where geometrical distortions in the image which arise from the strong gravitational field (i.e., spacetime curvature) of the black hole are negligible. To generate full 360 • images as seen by an observer close to the black hole, we have extended the procedure of [19] to use an orthonormal tetrad basis for the construction of initial photon wave vectors, distributing them uniformly as a function of θ ∈ [0, π] and φ ∈ [0, 2π] over a unit sphere.
The advantage of this approach is that all geometrical, relativistic, and generalrelativistic effects on the observed emission are naturally and self-consistently folded into the imaging calculation, providing a complete and physically-accurate depiction of what would really be seen from an observer's perspective.
The first step in building the tetrad basis is using a set of trial vectors (specifically, 4-vectors), t µ (a) , to find the tetrad basis vectors, e µ (a) . Herein, parenthesised lowercase Roman letters correspond to tetrad frame indices while Greek letters correspond to coordinate frame indices. Unless stated otherwise, all indices are taken to vary over 0-3, with 0 denoting the temporal component and 1-3 denoting the spatial components of a given 4-vector. Given a set of {θ, φ} pairs (typically on a uniform grid), the corresponding wave vector components in the tetrad frame, k (a) , are given by: where it is trivial to verify that this wave vector satisfies k (a) k (a) = 0, as expected for null geodesics.
In order to determine the wave vector defined in eqs. (1)-(4) in the coordinate frame, k α , it is necessary to first construct the tetrad vectors explicitly. The first trial vector we use is the four-velocity of the observer, t µ (0) = u µ obs . This vector is, by virtue of sensible initial conditions and preservation of the norm during integration, normalised. Using the four-velocity as an initial trial vector also ensures that Doppler effects due to the motion of the camera is included correctly. It is then possible to build a set of orthonormal basis vectors e µ (a) by using the Gram-Schmidt orthonormalisation procedure. The required trial vectors for this procedure are given by: = (0, 0, 1, 0) , This set of trial vectors is chosen such that the observer always looks towards the black hole in a right-handed basis. Any other initialisation, e.g., along with the velocity vector, could cause discomfort when used in VR due to high azimuthal velocities. The wave vector may now be found by taking the inner product of the tetrad basis vectors and the wave vector in the observer's frame as: The observer's camera is then initialised at a position X µ cam and uniformly-spaced rays are launched in all directions from this point. This method is fully covariant and is therefore valid in any coordinate system.

Black holes and gravitational lensing
In this work, we adopt geometrical units, G = M = c = 1, such that length and time scales are dimensionless. Hereafter M denotes the black hole mass, and setting M = 1 is equivalent to rescaling the length scale to units of the gravitational radius, r g := GM/c 2 , and the time scale to units of r g /c = GM/c 3 . To rescale lengths and times to physical units, one simply scales r g and r g /c using the appropriate black hole mass. For Sgr A* these scalings are given by r g = 5.906 × 10 11 cm and r g /c = 19.7 seconds, respectively.
The line element in GR determines the separation between events in space-time, and is defined as: where g µν is the metric tensor and dx µ an infinitesimal displacement vector. The metric is a geometrical object that contains all the information concerning the spacetime under consideration (in this study a rotating Kerr black hole) and is used to raise and lower tensor indices, e.g., g αµ A µν1ν2...νn = A ν1ν2...νn α , where the Einstein summation convention is implicitly assumed. The line element for a rotating black hole is given by the Kerr metric [40], which is written in Boyer-Lindquist coordinates x µ = (t, r, θ, φ) as: where and a is the dimensionless spin parameter of the black hole.
In the above form, the Kerr metric has a coordinate singularity at the outer (and inner) event horizon, which presents difficulties for both the numerical GRMHD evolution and the GRRT calculations. This also prohibits the observer's camera from passing smoothly through this region. To avoid this we transform (10) from x µ into horizon-penetrating Kerr-Schild coordinatesx µ = t ,r,θ,φ as: where In eq. (14) the outer horizon is given by r out ≡ 1 + √ 1 − a 2 , and the inner horizon by r in ≡ 1 − √ 1 − a 2 . Hereafter the coordinate system employed in this study is the modified Kerr-Schild (MKS) system, denoted by X µ , which is related to the aforementioned Kerr-Schild coordinates,x µ , as: To visualise the effect of a moving camera compared to a stationary camera, we calculate light rays originating from both a stationary observer and a free-falling observer. This calculation is performed at two different positions, which in MKS coordinates are given by: X µ 1 = (0, ln 10, 0, 0) and X µ 2 = (0, ln 3, 0, 0) .
Consequently, the observer positions 1 and 2 correspond to radial distances of 10 r g and 3 r g , respectively. An observer at rest has a four-velocity where α := (−g tt ) −1/2 is the lapse function. At the positions X µ 1 and X µ 2 the free-falling observer has the following corresponding four-velocity components: The free-falling velocities were obtained by numerically integrating the geodesic equation for a free-falling massive particle.
To visualise the effect of the observer's motion on the observed field of view, we place a sphere around both the observer and the black hole, which is centred on the black hole. This is what we subsequently refer to as the "celestial sphere". The black hole spin is taken to be a = 0.9375, the exact value of the spin parameter for Sgr A* is unknown, the chosen value was the best fit of a parameter survey [48]. The observer is positioned in the equatorial plane of the black hole (i.e., θ = 90 • ), where the effects of gravitational lensing are most significant and asymmetry in the shadow shape due to the rotational frame dragging arising from the spin of the black hole is most pronounced.
Each quadrant of the celestial sphere is then painted with a distinct colour and lines of constant longitude and latitude are included to aid in the interpretation of the angular size and distortion of the resulting images. The celestial sphere in Minkowski spacetime, where we used cartesian coordinates to integrate the geodesics, as seen by an observer positioned at 10 r g can be seen in Figure 1. The number of coloured patches in the θ and φ directions is (n θ , n φ ) = (8,16). Therefore, excluding the black lines of constant latitude and longitude (both 1.08 • in width), each coloured patch subtends an angle of 22.5 • in both directions. We also calculated 25 light rays for each of these observers, distributing them equally over (θ, φ) in the frame of the observer (see bottom rows of Figs. 2 & 3) in order to interpret the geometrical lensing structure of the images in terms of their constituent light rays. Figure 2 presents black hole shadow images and background lensing patterns for the Kerr black hole as seen by both a stationary observer (top panel) and a radially infalling observer (middle panel) located at a distance of 10 r g . The angular size of the shadow is larger for the stationary observer. This observer, being in an inertial frame, is essentially accelerating such that the local gravitational acceleration of the black hole is precisely counteracted by the acceleration of their reference frame. This gives rise to a force on the observer directed away from the black hole itself, reducing the angular momentum of photons oriented towards the black hole (seen as the innermost four rays being bent around the horizon), effectively increasing the black hole's capture cross-section and producing a larger shadow. Strong gravitational lensing of the image due to the presence of the compact mass of the black hole is evident in the warping of the grid lines.
In Figure 3 the observers are now placed at 3 r g , i.e., very close to the black hole. For the stationary observer, all photons within a field of view centred on the black hole of > 180 • in the horizontal direction and over the entire vertical direction, are captured by the black hole. Such an observer looking at the black hole would see nothing but the darkness of the black hole shadow in all directions. This is clear in the corresponding bottom-left plot of photon trajectories. As the observer approaches the event horizon the entire celestial sphere begins to focus into an ever shrinking point adjacent to the observer. For the infalling observer, the lensed image is far less extreme. Whilst the shadow presents a larger size in the observer's field of view, this is mostly geometrical, i.e., due to the observer's proximity to the black hole. There is also visible magnification of regions of the celestial sphere behind the observer. These results clearly follow from the photon trajectories in the bottom-right panel.
In all images of the shadow, repeated patches of decreasingly small area and identical colours are visible. In particular, multiple blue and yellow patches whose photons begin from behind the observer are visible near the shadow. These are a consequence of rays which perform one or more orbits of the black hole before reaching the observer, thereby appearing to originate from in front of the observer.

Camera trajectories
As described in Section 1, we consider two distinct phases for the camera trajectory. The first phase assumes a hovering observer positioned either at a fixed point or on a hovering trajectory around the black hole (i.e., the camera's motion is unaffected by the plasma motion and is effectively in an inertial frame). For the second phase of the trajectory, the observer's four-velocity is determined from an axisymmetric GRMHD simulation which includes tracer particles that follow the local plasma velocity. The choice to perform a separate tracer-particle simulation that is axisymmetric, in contrast to the 3D plasma simulation, was made to omit turbulent features in the φ direction which can be nauseating to watch in VR environments. This makes the movie scientifically less accurate, but is necessary to prevent viewers from experiencing motion sickness. Since the methods presented in this paper are not dependent on the dimensionality of the tracer particle simulation, they can be used for full 3D tracer particle simulations as well. In the following subsections, these two camera trajectories are described in detail.

Hovering trajectory
In the first phase of the trajectory, the observer starts in a vacuum, with only the light from the distant background stars being considered in the calculation. The observer is initially at a radius of 400 r g and moves inward to 40 r g . After this, the observer rotates around the black hole, which we term the "initialisation scene", and comprises 1600 frames. Each frame is separated by a time interval of 1 r g /c. The first phase of the movie, which includes the time-evolving accretion flow, consists of 2000 frames from the perspective of an observer at a radius of 40 r g and an inclination of 60 • with respect to the spin axis of the black hole. We refer to this first phase as "Scene 1". We then subsequently rotate around the black hole whilst simultaneously moving inward to a radius of 20 r g over a span of 1000 frames, which we refer to as "Scene 2". Within Scene 2, after the first 500 frames the observer then starts to decelerate until stationary once more.

Particle trajectory
For the second phase of the trajectory, the observer moves along a path that is calculated from an axisymmetric GRMHD simulation which includes tracer particles.   Fig. 2, now with the observer located at r = 3 rg. Differences between the shadow size and shape as seen by the two observers are now significant. See corresponding text for further discussion.
The tracer particles act like test masses: their velocity is found by interpolating the local plasma four-velocity (which is stored in a grid-based data structure) to the position of the particle. A first-order Euler integration scheme is then employed to update the position of each particle. For the camera, we are concerned with particles which are initially located within the accretion disk, begin to accrete towards the black hole, and then subsequently leave the simulation domain via the jet. To identify particles which satisfy all of these conditions we create a large sample of particle trajectories. The number of injected particles, N inj , within a grid cell with index {i, j} is set by two parameters: the plasma density, ρ, of the bounding cell, and the total mass, M tot , within the simulation domain. The number of injected particles is then calculated as where the weight factor ensures that only a predefined number of particles, N tot , after appropriate weighting, are then injected into a given simulation cell of volume V cell = √ −g dx 1 dx 2 dx 3 , where g is the determinant of the metric tensor. The code then randomly distributes these particles inside the simulation cell. The particles in the resulting movie. The blue square represents the initial position of the tracer particle used for the camera. The blue curve shows the trajectory corresponding to this tracer particle.
are initially in Keplerian orbits and co-rotate with the accretion disk. The disk then quickly becomes turbulent due to the growth of the magneto-rotational instability (MRI). As the particles are advected with the flow they can be classified into three different types: (1) accreted particles which leave the simulation at the inner radius (i.e., plunge into the event horizon) and remain gravitationally bound, (2) wind particles which become gravitationally unbound, travel through weakly magnetised regions and then exit the simulation at the outer boundary, (3) accelerated jet particles which are similar to wind particles but additionally undergo rapid acceleration within the highly-magnetised jet sheath. To discriminate between these three types of particle, several key hydrodynamical and magnetohydrodynamical criteria are examined. The first criterion is that the hydrodynamical Bernoulli parameter of the particle satisfies Bern = −hu t > 1.02, where h is the (specific) enthalpy of the accretion flow and u t is the covariant time component of the four-velocity. When this condition is satisfied the particle is, by definition, unbound. The boundary transition between bound and unbound happens at Bern = −hu t > 1.00, but we take a slightly larger value to select the part of the outflow that has a substantial relativisitc velocity. A similar value for the Bernoulli parameter was used in e.g. [52,57]. The second criterion is that the particle resides in high magnetisation regions where σ = B 2 /ρ > 0.1, where B := b µ b µ is the magnetic field strength and b µ is the magnetic field 4-vector. Satisfying this second criterion ensures that the particle ends up inside the jet sheath. The third criterion is that the particle's radial position is at a substantial distance from the black hole, typically r 300 r g , at the end of the simulation.
We simulate the particles with the axisymmetric GRMHD code HARM2D [41]. The simulation begins with N tot = 10 5 particles, a simulation domain size of r out = 1000 r g , and is evolved until t final = 4000 r g /c. The spacetime is that of a Kerr black hole, and the dimensionless spin parameter is set to be a = 0.9375. For this value of the spin, the black hole (outer) event horizon radius is r h = 1.344 r g and the simulation inner boundary lies within r h (i.e., we can track particles inside the event horizon). The specific particle used to initialise the camera trajectory is shown in Fig. 4 (blue square and curve). The full particle trajectory and velocity profile for all components u µ are shown in Fig. 5. Rapid variations in the azimuthal 4velocity, u 3 , as well as the angular velocity, Ω := u 3 /u 0 , in the right panel of Fig. 5 are consistent with the tightly wound trajectory in the left panel. This trajectory, which we term "Scene 3", begins immediately after Scene 2 (i.e. after frame 4600), and comprises 4000 frames, ending at frame 8599.

Radiative-transfer calculations and background images
To create images of an accreting black hole, it is necessary to compute the trajectories of light rays from the radiating plasma to the observer. For imaging applications, such as the present case, it is most computationally efficient to start the light rays at the observer instead -one for each pixel in the image the observer sees -and then trace them backward in time. Given a ray's trajectory, the radiative-transfer equation is solved along that trajectory, in order to compute the intensity seen by the observer. The radiative-transfer code RAPTOR uses a fourth-order Runge-Kutta method to integrate the equations of motion for the light rays (i.e., the geodesic equation). It simultaneously solves the radiative-transfer equation using a semianalytic scheme (for a more detailed description of RAPTOR, see Bronzwaer et al. [30]). The same methodology is applied here in order to create images of the black hole accretion disk, with one small addition. When accretion disks, which tend to be roughly toroidal in shape, are filmed against a perfectly black background, the Figure 5 Left panel: the trajectory of the tracer particle that is used to initialise the camera trajectory. Right panel: the velocity profile of the tracer particle. The velocity peaks when the particle is closest to the black hole, where the angular velocity is high. The time shown on the x-axis is the time range of the frames used for Scene 3. resulting animations fail to convey a natural sense of motion and scale for the observer as they orbit the black hole. In order to increase the immersiveness of the observer and provide a physically-realistic sense of scale and motion, the present work expands on the aforementioned radiative-transfer calculations by including an additional source of radiation in the form of a background star field that is projected onto the celestial sphere surrounding the black hole and observer. This is achieved by expressing the intensity received by the observer in Lorentzinvariant form and integrating this intensity from the camera to its point of origin within the plasma, i.e., eq. (37) in [30]. This can then be expressed in integral form (upon including a term for the background radiation) as where the optical depth along the ray is calculated as Here, I ν describes a ray's specific intensity, ν its frequency, and j ν and α ν refer respectively to the plasma emission and absorption coefficients evaluated along the ray, which is itself parametrised by the affine parameter, λ. The subscript "∞" denotes quantities evaluated at the outer integration boundary (i.e., far from the black hole), while the subscript "obs" refers to the observer's current location. The background radiation is encoded in the term I ν,∞ /ν 3 ∞ . The first term on the righthand-side of eq. (20) is constant and represents the intensity of the background radiation, weighted by the local optical depth. The second term on the right-handside of eq. (20) is evaluated at a given observer position, λ obs , and specifies the accumulated intensity of emitted radiation after taking into account the local emissivity and absorptivity of the accreting plasma. See [42], [23], [30] for further details.
A physical description of the radiation is needed for I ν,∞ /ν 3 ∞ . Since this quantity is projected onto the celestial sphere, it is a function of two coordinates (θ,φ). Note that for the ray coordinates, in the limit r → ∞, both θ →θ and φ →φ, i.e., spacetime is asymptotically flat. We also note that only rays which exit the simulation volume (as opposed to rays which plunge towards the horizon) are assigned a nonzero background intensity after integration. In order to evaluate I ν for a given ray, we therefore take the ray's (θ, φ) coordinates after the ray leaves the simulation volume, and use them as the coordinates (θ,φ) on the celestial sphere. Finally, we transform these coordinates into pixel coordinates (x, y) of a PNG image in order to evaluate the intensity. The transformation from celestial coordinates to pixel coordinates is given by where z ≡ floor(z) is the floor function (which outputs the greatest integer ≤ z), and W and H are the width and height (in pixels) of the background image, respectively.
Using the scheme described above, it is possible to fold the background radiation field directly into the radiative transfer calculations of the accretion disk plasma. A second approach is to render separate movies for both the background and for the plasma, create a composite image for all corresponding time frames between the two movies in post-processing, and then create the new composite movie from the composite images. We adopt the second approach in all results shown in this paper.
We have chosen a background that is obtained from real astronomical star data from the Tycho 2 catalogue which are not in the Galactic Plane. The original equirectangular RGB 3K image was generated by [43] and converted to a greyscale 2K image.

Plasma and radiation models
In this work, we seek to model the SMBH Sgr A*. To this end we use a black hole mass of M BH = 4.0 × 10 6 M [44], and a dimensionless spin parameter of a = 0.9375, consistent with the particle simulation. The plasma flow was simulated with the GRMHD code BHAC [45]. The simulation domain had an outer radius of r outer = 1000 r g . The simulation is initialised with a Fishbone-Moncrief torus [46] with an inner radius of r inner = 6 r g , and with a pressure maximum at r max = 12 r g . Magnetic fields were inserted as poloidal loops that follow iso-contours of density, and the initial magnetisation was low, i.e., β = P gas /B 2 = 100, where P gas is the gas pressure of the plasma. The simulation was performed in three dimensions, with a resolution of 256, 128, 128 cells in the r, θ and φ directions, respectively. We simulated the flow up to t = 7000 r g /c.
In this work we use, an electron model by [52] where the electrons are cold inside the accretion disk and hot inside the highly magnetized outflows. For the electron distribution function, we adopt a thermal distribution, where [57] showed that this model accurately describes the quiescent state of Sgr A*. The used model [52] is capable of recovering the observational parameters of Sgr A*, such as radio fluxes and intrinsic source sizes [2,4,5,14].
We calculated the synthetic images at four different radio frequencies: 22 GHz (1.2 cm), 43 GHz (7 mm), 86 GHz (3 mm), and 230 GHz (1.3 mm). These frequencies were chosen since they correspond to the frequencies at which, e.g., the Very Long Baseline Array (VLBA) (1.2 mm, 7 mm, 3 mm), Global mm-VLBI Array (GMVA) (3 mm) and the Event Horizon Telescope (EHT) (1.3 mm) operate. After ray-tracing these frequencies were converted into separate PNG image files, where distinct colourmaps were chosen for each of the four frequencies. In postprocessing, these images were then combined into a single image by averaging over the RGB channels of the four different input images. A star-field background was also included to serve as a reference point for the observer during their motion. This star-field background was rendered separately from the radio images, although the opacity at 22 GHz was used to obscure stars located behind the accretion disk. This background was then also averaged together with the radio images using the same RGB channel averaging. The four separate frequencies, the star-field background, and the resulting combined image are presented in Fig. 6.

VR movie
The resulting VR movie contains 8600 frames at a resolution of 2000 × 1000 pixels. As a proof of concept, this resolution was chosen to balance image quality and computational resources. Current VR headsets also upscale the provided resolution with interpolation routines. We tested the resolution with the Oculus VR headset, which turned out to be sufficient. Since the provided methods are not limited by the resolution, a larger resolution can in principle be achieved. The movie is available on Youtube VR [56]. In this section, we discuss several snapshots from this movie.
The first set of snapshots is shown in 7. In Fig. 7 we show a set of snapshots from Scene 1, (1600, 2300, 3000), matter starts to accrete onto the black hole and the jet is launched. The jet then propagates through the ambient medium of the simulation, forming a collimated funnel that is mainly visible at lower frequencies. Since the accretion rate peaks at this point in the simulation (see Fig. 8), the black-hole shadow is barely visible.
In Fig. 9 we show snapshots from Scene 2 (3700, 4050, 4400), the jet propagates outward to the boundary of our simulation domain, the accretion rate settles and the black hole shadow becomes visible.
In Fig. 10 we show snapshot from Scene 3. When the observer moves along with the flow in Scene 3 (5100, 5800, 6150), small hot blobs of plasma orbiting the black hole are distinguishable. At closest approach (around 6 r g , frame 6150), the scene changes rapidly. This is due not only to rapid rotation of the black hole but also to the rapid decrease of observed flux. It is hard to distinguish individual stars and the only observable emission is at 230 GHz. At the end of Scene 3 (7200, 7900, 8599) the observer exits the accretion disk via the jet, whereafter a rapid increase in radial velocity is clearly seen. To obtain a better quantitative understanding of the movie we also calculate the total bolometric luminosity as received by the observer's camera. This is shown in the top panel of Fig. 11. At 6150 a decrease in luminosity is evident at the three lowest frequencies, which corresponds to where the observer is closest to the black hole event horizon and has entered the optically-thick accretion disk. A magnified version of this Figure in the optically-thick part is shown in the bottom panel of Fig. 11. A frame corresponding to this particular moment is shown in Fig. 10, panel 6150. At closest approach, the total luminosity detected at 230 GHz peaks, and the observer is exposed to ≈ 25L .

Discussion and conclusion
In this work, we have detailed our methods for visualising the surroundings of accreting black holes in virtual reality. We presented a visualisation of a three-dimensional fully-general-relativistic accreting black hole simulation in a full 360 • VR movie with radiative models based on physically-realistic GRMHD plasma simulations. In order to produce representative images, the radiative-transfer capabilities of our code RAPTOR were extended to include background starlight and an observer in an arbitrary state of motion. To model the emission emerging from the vicinity of a black hole we coupled the GRMHD simulation with our radiative-transfer code to produce a VR movie based on our recent models for Sgr A* [52,57]. These methods can be applied to accreting black holes of any size, so long as radiation feedback onto the accretion flow has a negligible impact on the flow's magnetohydrodynamical properties.
The trajectory of the camera consisted of two phases: a hovering observer and an advected observer. For this second phase, we used an axisymmetric GRMHD simulation, in contrast to the plasma simulation used to calculate the radiation, which was fully-three-dimensional. This choice, whilst scientifically less accurate, was intentional and somewhat necessary. Turbulent features in the φ direction were omitted since they can be nauseating to watch in VR environments and commonly lead to motion sickness. A composition of starfield and accretion flow images at four frequencies was then used to create a movie, consisting of 8600 frames, which is freely available on YouTube.
This movie couples GRMHD simulations with GRRT post-processing in VR. Since we do not make any strong a-priori assumptions regarding the field-of-view of the observer, we can calculate the full radiation field measured at a specific point in the accretion disk, where we include all GR effects. This enabled us to calculate light curves of the total measured luminosity at multiple frequency bands at the position of a particle being advected in the flow. This way of calculating the full self-irradiation of the disk is of potential interest in, e.g., studies of X-ray reflection models in AGN, or coupling to GRMHD simulation to calculate the proper radiative feedback onto an emitting, absorbing (and even scattering) plasma in GR in a selfconsistent way.
Finally, beyond the aforementioned scientific applications, VR represents a new medium for scientific visualisation which can be used, as demonstrated in this work, to investigate the emission that an observer would measure from inside the accretion flow. It is natural, and of contemporary interest even in the film industry [see e.g. has not yet begun, which can be seen as the faint, stationary equilibrium accretion torus configuration in the centre of the image. By frame 2300 accretion has begun (see also Fig. 8) and the dim jet (upper half of image) and dimmer counter jet (lower half of image) propagate outwards through the ambient medium. At frame 3000 the jet has propagated further outwards, and angular momentum transport has shifted torus material outward, as can be seen by the increased angular size of the inner accretion flow. The black hole shadow is not visible since the accretion rate has yet to reach a quasi-stationary state. Figure 8 Simulation accretion rate as a function of time (in code units). At t=2500 the MRI start to saturate. The time shown on the x-axis is the time of the frames used for "Scene 2" and "Scene 3". 58, 59] to ask the question as to what an observer would see if they were in the immediate vicinity of a black hole. In this work, we have sought to address this question directly, by using state-of-the-art numerical techniques and astrophysical models in a physically-self-consistent manner. Given the EHTC is anticipated to obtain images of the black hole shadows in Sgr A* and M87 in the near future, the calculations we have presented are timely. The VR movies presented in this work also provide an intuitive and interactive way to communicate black hole physics to wider audiences, serving as a useful educational tool.

Availability of Data and Materials
Please contact author for data requests.

Competing interests
The authors declare that they have no competing interests.

Figure 9
Movie snapshots from Scene 2. By frame 3700 the MRI has begun to saturate and the accretion rate reaches a quasi-stationary state. At frame 4050 the jet and counter-jet have propagated further away from the black hole and reached the boundary of our simulation domain. Due to the steadier accretion rate, by frame 4400 the central region surrounding the event horizon becomes cooler and more optically thin. The upper-half of the black-hole shadow is now visible.

Figure 10
Movie snapshots from Scene 3. The observer now begins their journey through the accretion flow (panels with frames 5100-6150), before being advected away from the black hole via the large-scale jet (panels with frames 7200-8599). At frame 6150 the observer is at their point of closest approach to the black hole, where the incident flux is as high as ≈ 25L . This region is highly optically thick, completely obscuring the observer's view of the black hole shadow.
As the observer is advected further away, by frame 8599 the angular size of the black hole and the surrounding accretion flow is greatly reduced and appears almost point-like.