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.

VR camera

The original RAPTOR code (Bronzwaer et al. 2018) initialises rays (i.e., photon geodesics) using impact parameters determined form coordinate locations on the observer’s image plane (Bardeen et al. 1972). 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 Noble et al. (2007) to use an orthonormal tetrad basis for the construction of initial photon wave vectors, distributing them uniformly as a function of \(\theta\in [0,\pi]\) and \(\phi\in[0,2\pi]\) over a unit sphere.

The advantage of this approach is that all geometrical, relativistic, and general-relativistic 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)}^{\mu}\), to find the tetrad basis vectors, \(e_{(a)}^{\mu}\). 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 \(\{\theta,\phi \}\) pairs (typically on a uniform grid), the corresponding wave vector components in the tetrad frame, \(k^{(a)}\), are given by:

$$\begin{aligned} &k^{(0)} = +1, \end{aligned}$$ (1)

$$\begin{aligned} &k^{(1)} = -\cos(\phi) \cos(\theta), \end{aligned}$$ (2)

$$\begin{aligned} &k^{(2)} = -\sin(\theta), \end{aligned}$$ (3)

$$\begin{aligned} &k^{(3)} = -\sin(\phi)\cos(\theta), \end{aligned}$$ (4)

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^{\alpha}\), 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)}^{\mu} = u^{\mu}_{\mathrm {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)}^{\mu}\) by using the Gram–Schmidt orthonormalisation procedure. The required trial vectors for this procedure are given by:

$$\begin{aligned} &t_{(1)}^{\mu}= (0, -1, 0, 0), \end{aligned}$$ (5)

$$\begin{aligned} &t_{(2)}^{\mu}=(0, 0, 1, 0), \end{aligned}$$ (6)

$$\begin{aligned} &t_{(3)}^{\mu}= (0, 0, 0, 1). \end{aligned}$$ (7)

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:

$$ k^{\mu}= e_{(a)}^{\mu}k^{(a)}. $$ (8)

The observer’s camera is then initialised at a position \(X_{\mathrm {cam}}^{\mu}\) 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_{\mathrm {g}}:= GM/c^{2}\), and the time scale to units of \(r_{\mathrm {g}}/c = GM/c^{3}\). To rescale lengths and times to physical units, one simply scales \(r_{\mathrm {g}}\) and \(r_{\mathrm {g}}/c\) using the appropriate black hole mass. For Sgr A* these scalings are given by \(r_{\mathrm {g}} = 5.906 \times 10^{11}\) cm and \(r_{\mathrm {g}}/c = 19.7\) seconds, respectively.

The line element in GR determines the separation between events in space-time, and is defined as:

$$ ds^{2} = g_{\mu

u} \,dx^{\mu}dx^{

u}, $$ (9)

where \(g_{\mu

u}\) is the metric tensor and \(dx^{\mu}\) an infinitesimal displacement vector. The metric is a geometrical object that contains all the information concerning the space-time under consideration (in this study a rotating Kerr black hole) and is used to raise and lower tensor indices, e.g., \(g_{\alpha\mu}A^{\mu

u_{1}

u_{2}\ldots

u_{\mathrm {n}}} = A^{

u_{1}

u_{2}\ldots

u_{\mathrm {n}}}_{\alpha}\), where the Einstein summation convention is implicitly assumed. The line element for a rotating black hole is given by the Kerr metric (Kerr 1963), which is written in Boyer–Lindquist coordinates \(x^{\mu}= ( t,r,\theta,\phi )\) as:

$$\begin{aligned} ds^{2}={}&{-} \biggl(1-\frac{2r}{\varSigma} \biggr) \,dt^{2} - \frac{4ar\sin ^{2}\theta}{\varSigma}\,dt \,d\phi +\frac{\varSigma}{\Delta}\,dr^{2} \\ &{}+ \varSigma\, d \theta^{2} \\ &{}+ \biggl(r^{2}+a^{2}+ \frac{2ra^{2}\sin^{2}\theta}{\varSigma} \biggr)\sin ^{2}\theta \,d\phi^{2}, \end{aligned}$$ (10)

where

$$\begin{aligned} &\Delta:= r^{2} - 2r + a^{2}, \end{aligned}$$ (11)

$$\begin{aligned} &\varSigma:=r^{2} +a^{2} \cos^{2}\theta, \end{aligned}$$ (12)

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^{\mu}\) into horizon-penetrating Kerr–Schild coordinates \(\tilde{x}^{\mu}= (\tilde{t}, \tilde{r}, \tilde {\theta}, \tilde{\phi} )\) as:

$$ \begin{aligned} &\tilde{t} = t + \ln\Delta+ 2 \mathcal{R},\qquad \tilde{r} = r, \\ &\tilde{\theta} = \theta,\qquad \tilde{\phi} = \phi+ a \mathcal{R}, \end{aligned}$$ (13)

where

$$ \mathcal{R} \equiv\frac{1}{r_{\mathrm {out}}}-r_{\mathrm {in}}\ln \biggl( \frac{r-r_{\mathrm {out}}}{r-r_{\mathrm {in}}} \biggr). $$ (14)

In Eq. (14) the outer horizon is given by \(r_{\mathrm {out}} \equiv1 + \sqrt{1-a^{2}}\), and the inner horizon by \(r_{\mathrm {in}} \equiv1 - \sqrt{1-a^{2}}\). Hereafter the coordinate system employed in this study is the modified Kerr–Schild (MKS) system, denoted by \(X^{\mu}\), which is related to the aforementioned Kerr–Schild coordinates, \(\tilde{x}^{\mu}\), as:

$$ X^{0} = \tilde{t},\qquad X^{1} = \ln\tilde{r},\qquad X^{2} = \tilde{\theta }/\pi,\qquad X^{3} = \tilde{\phi}. $$ (15)

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^{\mu}_{1} = ( 0, \ln10, 0, 0) \quad\text{and}\quad X^{\mu}_{2} = ( 0, \ln3, 0, 0). $$ (16)

Consequently, the observer positions 1 and 2 correspond to radial distances of \(10~r_{\mathrm { g}}\) and \(3~r_{\mathrm { g}}\), respectively. An observer at rest has a four-velocity

$$ u^{\mu}_{0} = ( \alpha, 0, 0, 0 ), $$ (17)

where \(\alpha:= (-g^{tt} )^{-1/2}\) is the lapse function. At the positions \(X^{\mu}_{1}\) and \(X^{\mu}_{2}\) the free-falling observer has the following corresponding four-velocity components:

$$ \begin{aligned} &u^{\mu}_{1} = ( 1.10, -0.029, 0, -0.0011) \quad\text{and}\\ &u^{\mu}_{2} = ( 1.34, -0.26, 0, -0.034). \end{aligned}$$ (18)

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 (Mościbrodzka et al. 2009). The observer is positioned in the equatorial plane of the black hole (i.e., \(\theta=90^{\circ}\)), 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_{\mathrm { g}}\) can be seen in Fig. 1. The number of coloured patches in the θ and ϕ directions is \((n_{\theta},n_{\phi} )= (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 \((\theta,\phi )\) 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 1 Celestial sphere in Minkowksi spacetime for an observer at \(r=10~r_{\mathrm { g}}\). The different colors represent different quadrants of the sky, with yellow and blue being behind the observer, while red and green are in front of the observer. The black lines represent lines of constant longitude and lattitude Full size image

Figure 2 Celestial sphere and black hole shadow images for an observer located at \(r=10~r_{\mathrm { g}}\). Top panel: celestial sphere and shadow image as seen by a stationary observer. The different colors represent different quadrants of the sky, yellow and blue being behind the observer, while red and green are in front of the observer. The black lines represent lines of constant longitude and lattitude while the black, circular region in the center is the black-hole shadow. Middle panel: as top panel, but seen by a radially in-falling observer. Bottom-left panel: photons originating from a stationary observer’s camera, as used to generate the top panel. Bottom-right panel: photons originating from a radially in-falling observer’s camera, as used to generate the middle panel. The black hole event horizon is shown as the black region in both bottom panels. The shadow sizes are similar in both panels, but differences are clearly visible. See corresponding text for further discussion Full size image

Figure 3 Celestial sphere and black hole shadow images for an observer located at \(r=3~r_{\mathrm { g}}\). As in Fig. 2, now with the observer located at \(r=3~r_{\mathrm { g}}\). Differences between the shadow size and shape as seen by the two observers are now significant. See corresponding text for further discussion Full size image

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_{\mathrm { 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 Fig. 3 the observers are now placed at \(3~r_{\mathrm { 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 Sect. 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_{\mathrm { g}}\) and moves inward to \(40~r_{\mathrm { 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_{\mathrm { 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_{\mathrm { 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_{\mathrm { 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. 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, \({\mathcal{N}}_{\mathrm {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_{\mathrm {tot}}\), within the simulation domain. The number of injected particles is then calculated as

$$ {\mathcal{N}}_{\mathrm {inj}} (i,j ) = N_{\mathrm {tot}} \biggl( \frac {\rho (i,j ) V_{\mathrm {cell}}}{M_{\mathrm {tot}}} \biggr), $$ (19)

where the weight factor ensures that only a predefined number of particles, \(N_{\mathrm {tot}}\), after appropriate weighting, are then injected into a given simulation cell of volume \(V_{\mathrm {cell}} = \sqrt{-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 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 \(\operatorname{ Bern} = - h u_{\mathrm{ t}} > 1.02\), where h is the (specific) enthalpy of the accretion flow and \(u_{\mathrm{ 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 \(\operatorname{ Bern} = - h u_{\mathrm{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. Mościbrodzka et al. (2014); Davelaar et al. (2018b). The second criterion is that the particle resides in high magnetisation regions where \(\sigma= B^{2}/\rho> 0.1\), where \(B:=\sqrt{b_{\mu} b^{\mu}}\) is the magnetic field strength and \(b^{\mu}\) 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 \gtrsim300~r_{\mathrm { g}}\), at the end of the simulation.

We simulate the particles with the axisymmetric GRMHD code HARM2D (Gammie et al. 2003). The simulation begins with \(N_{\mathrm {tot}} = 10^{5}\) particles, a simulation domain size of \(r_{\mathrm {out}} = 1000~r_{\mathrm { g}}\), and is evolved until \(t_{\mathrm {final}}=4000~r_{\mathrm { 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_{\mathrm{ h}}=1.344~r_{\mathrm { g}}\) and the simulation inner boundary lies within \(r_{\mathrm{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^{\mu}\) are shown in Fig. 5. Rapid variations in the azimuthal 4-velocity, \(u^{3}\), as well as the angular velocity, \(\varOmega:=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.

Figure 4 Snapshots of the tracer particles in the advection simulation. Left panel: initial distribution of particles inside the initial torus. Middle panel: snapshot of the advection HARM2D simulation at \(t=2000~r_{\mathrm { g}}/c\). Right panel: later snapshot at \(t=4000~r_{\mathrm { g}}/c\). The two times correspond to the advection simulation time, i.e., frames 4600–7600 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 Full size image

Figure 5 Used tracer-particle trajectory and it’s corresponing velocity. 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 Full size image

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 semi-analytic scheme (for a more detailed description of RAPTOR, see Bronzwaer et al. (2018)). 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 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 Lorentz-invariant form and integrating this intensity from the camera to its point of origin within the plasma, i.e., Eq. (37) in Bronzwaer et al. (2018). This can then be expressed in integral form (upon including a term for the background radiation) as

$$ \frac{I_{

u,\mathrm{obs}}}{

u_{\mathrm {obs}}^{3}} = \biggl(\frac{I_{

u,\infty}}{

u_{\infty}^{3}} \biggr) {\mathrm{ e}}^{-\tau_{

u,\mathrm {obs}} ( \lambda_{\infty} )} + \int^{\lambda_{\infty}}_{\lambda_{\mathrm {obs}}} \biggl(\frac{j_{

u}}{

u^{2}} \biggr) { \mathrm{e}}^{-\tau_{

u,\mathrm{obs}}} \,d{\lambda'}, $$ (20)

where the optical depth along the ray is calculated as

$$ \tau_{

u,\mathrm{obs}} (\lambda ) = \int^{\lambda} _{\lambda_{\mathrm {obs}}}

u\alpha_{

u}\,d \lambda'. $$ (21)

Here, \(I_{

u}\) describes a ray’s specific intensity, ν its frequency, and \(j_{

u}\) and \(\alpha_{

u}\) 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_{

u,\infty}/

u _{\infty}^{3}\). The first term on the right-hand-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-hand-side of Eq. (20) is evaluated at a given observer position, \(\lambda_{\mathrm {obs}}\), and specifies the accumulated intensity of emitted radiation after taking into account the local emissivity and absorptivity of the accreting plasma. See Fuerst and Wu (2004); Younsi et al. (2012); Bronzwaer et al. (2018) for further details.

A physical description of the radiation is needed for \(I_{

u,\infty} /

u_{\infty}^{3}\). Since this quantity is projected onto the celestial sphere, it is a function of two coordinates \((\hat{\theta},\hat{\phi})\). Note that for the ray coordinates, in the limit \(r\to\infty\), both \(\theta\to\hat{\theta}\) and \(\phi\to\hat{\phi}\), i.e., space-time 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 non-zero background intensity after integration. In order to evaluate \(I_{

u}\) for a given ray, we therefore take the ray’s \((\theta,\phi )\) coordinates after the ray leaves the simulation volume, and use them as the coordinates \((\hat{\theta},\hat{\phi})\) 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

$$ x = \biggl\lfloor \frac{\hat{\phi}}{2 \pi} W \biggr\rfloor \quad\text{and}\quad y = \biggl\lfloor \frac{\hat{\theta}}{\pi} H \biggr\rfloor , $$ (22)

where \(\lfloor z \rfloor\equiv\mathrm{ 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 Scott (2008) 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_{\mathrm{ BH}} = 4.0\times10^{6}~\mathrm{M}_{\odot}\) (Gillessen et al. 2009), and a dimensionless spin parameter of \(a=0.9375\), consistent with the particle simulation. The plasma flow was simulated with the GRMHD code BHAC (Porth et al. 2017). The simulation domain had an outer radius of \(r_{\mathrm{outer}} = 1000~r_{\mathrm { g}}\). The simulation is initialised with a Fishbone–Moncrief torus (Fishbone and Moncrief 1976) with an inner radius of \(r_{\mathrm{inner}}=6~r_{\mathrm { g}}\), and with a pressure maximum at \(r_{\mathrm{ max}} = 12~r_{\mathrm { g}}\). Magnetic fields were inserted as poloidal loops that follow iso-contours of density, and the initial magnetisation was low, i.e., \(\beta= P_{\mathrm{ gas}}/B^{2} = 100\), where \(P_{\mathrm{ 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_{\mathrm { g}}/c\).

The GRMHD simulation only simulates the dynamically-important ions (protons). We, therefore, require a prescription for the radiatively-important electrons in order to compute the observed emission. Most radiative models for Sgr A* or M87 either assume that the coupling between the temperatures of the electrons and protons is constant or parameterised based on plasma variabels, see e.g. Goldston et al. (2005); Noble et al. (2007); Mościbrodzka et al. (2009); Dexter et al. (2010); Shcherbakov et al. (2012); Mościbrodzka and Falcke (2013); Mościbrodzka et al. (2014); Chan et al. (2015a, 2015b); Gold et al. (2017). In this work we use, an electron model by Mościbrodzka et al. (2014) 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 Davelaar et al. (2018b) showed that this model accurately describes the quiescent state of Sgr A*. The used model (Mościbrodzka et al. 2014) is capable of recovering the observational parameters of Sgr A*, such as radio fluxes and intrinsic source sizes (Falcke et al. 2000; Bower et al. 2004, 2014; Doeleman et al. 2008).

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 post-processing, 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.