© GRAVITY Collaboration et al. 2020

1. Introduction

The black hole at the Galactic center provides a unique laboratory for exploring the physics of accretion and supermassive black holes. As it is the closest supermassive black hole, Sgr A* has the largest angular size of any black hole observable from Earth, with a Schwarzschild radius of R_s = 10.022 ± 0.020_stat. ± .032_sys. μas (Gravity Collaboration 2019). Additionally, the characteristic timescale of infalling matter near the black hole is on the order of minutes to hours (e.g., Baganoff et al. 2003; Meyer et al. 2009; Witzel et al. 2012, 2018; Hora et al. 2014), allowing variability and motion to be observed within a single night.

The GRAVITY instrument at the VLT Interferometer (Eisenhauer et al. 2008; Gravity Collaboration 2017) allows for high-precision astrometry of the Galactic center. This instrument has proven capable of an astrometric precision of 10 − 30 μas during bright flares (Gravity Collaboration 2018a). The observation of coherent circular motion of the centroid of emission during infrared flares lends new insight into the geometry and physical properties of matter in the accretion flow, as well as the spacetime near the black hole.

In the past, efforts to predict near-infrared flares of the accretion onto Sgr A* have fallen into a few broad categories. The first of these involves detailed general relativistic magnetohydrodynamic simulations of a small volume around the black hole (e.g., Dolence et al. 2009; Dexter & Fragile 2013; Chan et al. 2015; Ball et al. 2016; Ressler et al. 2017). These simulations allow for self-consistent modeling of the hydrodynamics of the accretion flow in the gravitational field of the black hole. However, the computational expense of these models preclude exploration of a large parameter space, making it difficult to constrain models based on observations. Simulations are usually able to produce only a limited number of orbits, making them sensitive to the initial conditions. Moreover, these models often fail to reproduce observed properties of Sgr A* such as the wavelength-dependent variability or the polarization fraction or angle during flares.

Yuan et al. (2003) developed semi-analytic models of the spectral energy distribution and multi-wavelength variability of Sgr A*. These models successfully match the global observed properties of the black hole, but they do not explain the mechanism causing flares nor do they predict the distribution of matter near the black hole.

Another class of recent models consist of Particle-in-Cell (PIC) simulations of a small region of the accretion flow (e.g., Sironi & Narayan 2015; Rowan et al. 2017, 2019). These models provide insight into the microphysical properties of the accretion flow, but they are computationally limited to very small regions. Predicting global observables, therefore, is not yet possible.

The last class of models consists of semianalytic calculations of hotspots in orbits around the central black hole (e.g., Broderick & Loeb 2005, 2006; Paumard et al. 2005; Meyer et al. 2006; Trippe et al. 2007; Hamaus et al. 2009; Zamaninasab et al. 2010). These models are computationally cheap enough to make it possible to explore a large parameter space of orbital configurations and black hole properties. However, they generally do not include any magnetohydrodynamic effects and, therefore, they cannot lend any insight into, for example, the mechanism that causes flares or the microphysical properties of the accretion flow.

In the past, semianalytic hotspot models have been used to perform numerical calculations of photon geodesics through the curved spacetime near the black hole, but they have been limited to particle orbits that can be calculated analytically. In practice, this has restricted calculations to a small family of orbits–generally, circular orbits in the equatorial plane of the black hole. However, there is no a priori reason to expect the transiently heated, relativistic electrons responsible for infrared flares to be confined to the midplane or to circular orbits.

In this paper, we describe the NERO code, which generalizes the hotspot model to allow for arbitrary particle orbits in the vicinity of a Kerr black hole. Models with non-circular, unbound, or outflowing orbits are similarly allowed in our simulations. We discuss the methods for calculating photon (Sect. 2) and particle (Sect. 3) geodesics to generate synthetic flare observables incorporating all relativistic effects. We use the resulting model to analyze the GRAVITY data from 2018 (Sect. 4). We consider both circular orbits as in past works (Gravity Collaboration 2018a) and obtain new constraints on out-of-plane motion (Sect. 5), as well as the size of the emission region due to shear (Sect. 6).

2. Ray tracing

The first step in calculating the appearance of the region near the black hole at the Galactic center is to find the paths followed by the photons as they travel towards a distant observer. We use the GEOKERR code (Dexter & Agol 2009) to calculate null geodesics between the region around the black hole and a distant image plane. We assume a Kerr metric characterized by the line element

$$\begin{array}{cc}\hfill \mathrm{d}{s}^2=& -(1-\frac{2Mr}{\mathrm{\Sigma }})\mathrm{d}{t}^2-\left(\frac{4Mar{sin}^2\theta }{\mathrm{\Sigma }}\right)\mathrm{d}t\mathrm{d}\varphi \hfill \\ \hfill & +\left(\frac{\mathrm{\Sigma }}{\mathrm{\Delta }}\right)\mathrm{d}{r}^2+\mathrm{\Sigma }\mathrm{d}{\theta }^2\hfill \\ \hfill & +(r^2+a^2+\frac{2M{a}^{2}r{sin}^2\theta }{\mathrm{\Sigma }}){sin}^2\theta \mathrm{d}{\varphi }^2\hfill \end{array}$$(1)

where Σ = r² + a²cos²θ, Δ = r² − 2Mr + a², M is the mass of the black hole, and a = J/M² is the dimensionless angular momentum. Here and throughout this paper, we use the convention G = c = 1.

There are two possible methods for finding particle trajectories in curved spacetimes. The more general approach is to directly solve the geodesic equations numerically. Alternatively, in certain astrophysically relevant metrics, the equations of motion can be reduced analytically to elliptical integrals (Cunningham & Bardeen 1973; Rauch & Blandford 1994). These integrals can then be calculated semi-analytically to find the particle trajectories.

The GEOKERR code adopts the latter approach to solve the geodesic equations in the Kerr metric (Dexter & Agol 2009). Starting from an image plane at the location of a distant observer, the code solves the geodesic equations backwards in time to find the trajectories of the photons originating near the black hole that terminate on the image plane.

3. Relativistic orbits

The second component of the NERO code is the calculation of the relativistic orbits of test particles near the event horizon of the black hole. We are particularly interested in calculating a broad range of orbits without constraining particles to be in circular orbits or in the equatorial plane of the black hole.

Since we are interested in highly relativistic orbits close to the black hole, the usual Keplerian elements are a poor choice for characterizing orbits. Instead, we follow Hughes (2001) and identify each orbit by its energy E, z angular momentum L_z, and Carter constant Q in addition to choosing the spin a of the black hole. While the procedure we use is capable of calculating general orbits with arbitrary inclination and eccentricity, we confine our initial analysis to circular orbits. In this case, we use the orbital radius R as an orbital parameter along with L_z. We sample R between the innermost stable circular orbit and an arbitrary large value (we choose R = 14 for the analysis below).

In order to sample the range of allowed angular momenta, we first calculate the specific angular momentum of a prograde and retrograde orbit, respectively,

$$\begin{array}{cc}& L_z^{\mathrm{pro}}=\sqrt{R}\frac{1-2a{R}^{-3}+a^2{R}^{-4}}{\sqrt{1-3R^{-2}+2a{R}^{-3}}},\hfill \end{array}$$(2)

$$\begin{array}{cc}& L_z^{\mathrm{ret}}=\sqrt{R}\frac{1+2a{R}^{-3}+a^2{R}^{-4}}{\sqrt{1-3R^{-2}+2a{R}^{-3}}}·\hfill \end{array}$$(3)

We can then create a family of orbits at a given radius by sampling L_z evenly between these bounds to obtain orbits with a range of inclinations to the spin vector of the black hole.

The condition of circularity on the orbits introduces the constraints,

$$\begin{array}{c}\hfill \frac{\mathrm{d}\mathcal{R}}{\mathrm{d}r}=\mathcal{R}=0,\end{array}$$(4)

where ℛ is the effective potential in the radial direction.

$$\begin{array}{c}\hfill \mathcal{R}=T^2-\mathrm{\Delta }(r^2+{(L-aE)}^2+Q).\end{array}$$(5)

These constraints allow us to find the specific energy and specific Carter constant for each value of L_z and R. We can now find the energy and Carter constant of the orbit (Hughes 2001),

$$\begin{array}{cc}\hfill E=& (a^2{L}_z^2(R-M)+R{\mathrm{\Delta }}^2)/(a{L}_{z}M(R^2-a^2)\hfill \\ \hfill & +\sqrt{R^5(R-3M)+a^{4}R(R+M)+a^2{R}^2(L_z^2-2MR+2R^2)}),\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill Q=& \frac{{[(a^2+R^2)E-a{L}_z]}^2}{\mathrm{\Delta }}-(R^2-a^2{E}^2-2aE{L}_z+L_z^2).\hfill \end{array}$$(7)

Once we have calculated the constants of motion, we check to make sure the orbit is stable. The criterion for stability is

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2\mathcal{R}}{\mathrm{d}{r}^2}<0.\end{array}$$(8)

Orbits in the region of the parameter space that does not satisfy Eq. (8) are unstable and therefore astrophysically irrelevant.

For each orbit we wish to model, we have now calculated the constants of motion E, L_z, and Q. The equations of motion for a test particle following a timelike geodesic through the Kerr metric in terms of these constants are given by (Carter 1968):

$$\begin{array}{cc}& \frac{\mathrm{d}t}{\mathrm{d}\tau }=\frac{1}{\mathrm{\Sigma }}(-a(aE{sin}^2\theta -L)+\frac{r^2+a^{2}T}{\mathrm{\Delta }})\hfill \end{array}$$(9)

$$\begin{array}{cc}& \frac{\mathrm{d}r}{\mathrm{d}\tau }=\pm \frac{\sqrt{\mathcal{R}}}{\mathrm{\Sigma }}\hfill \end{array}$$(10)

$$\frac{\text{d}\theta }{\text{d}\tau }=\pm \frac{\sqrt{\times }}{\Sigma }$$(11)

$$\begin{array}{cc}& \frac{\mathrm{d}\varphi }{\mathrm{d}\tau }=\frac{1}{\mathrm{\Sigma }}(-a(E-\frac{L}{{sin}^2\theta }+\frac{\mathit{aT}}{\mathrm{\Delta }}))·\hfill \end{array}$$(12)

Here, τ is the proper time, Θ is the effective potential in the angular direction,

$$\begin{array}{c}\hfill \mathrm{\Theta }=Q-{cos}^2\theta (a^2(1-E^2)+\frac{L^2}{{sin}^2\theta }),\end{array}$$(13)

and T is

$$\begin{array}{c}\hfill T=E(r^2+a^2)-La.\end{array}$$(14)

We calculate the solution to these equations of motion using the publicly-available YNOGKM code (Yang & Wang 2014). This code adapts the GEOKERR algorithm to semi-analytically find particle trajectories in the Kerr-Newman metric for a given set of initial conditions. Since we are concerned with orbits about Sgr A*, we set the charge of the black hole to zero, reducing the metric to the Kerr solution given in Eq. (1).

The YNOGKM code takes as input the initial position and velocity of a test particle. As an initial position, we arbitrarily choose a point in the equatorial plane (θ = π/2) with a phase ϕ = π/2, with r set to the orbital radius R. The initial velocity in the Boyer-Lindquist frame can be calculated from the constants of motion:

$$\begin{array}{cc}& v_{\mathrm{BL}}^{\mathrm{r}}=\pm \frac{\sqrt{\mathcal{R}}}{{\gamma }_{\mathrm{p}}E\sqrt{\mathrm{\Sigma }\mathrm{\Delta }}},\hfill \end{array}$$(15)

$$\begin{array}{cc}& v_{\mathrm{BL}}^{\theta }=\pm \frac{1}{{\gamma }_{\mathrm{p}}E\sqrt{\mathrm{\Sigma }}}\sqrt{Q-{cos}^2\theta (a^2(E^2-1)+\frac{{\lambda }^2}{{sin}^2\theta })},\hfill \end{array}$$(16)

$$\begin{array}{cc}& v_{\mathrm{BL}}^{\varphi }=\frac{\lambda }{{\gamma }_{\mathrm{p}}E}\sqrt{\frac{\mathrm{\Sigma }}{A{sin}^2\theta }},\hfill \end{array}$$(17)

where

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{p}}=\sqrt{\frac{A}{\mathrm{\Delta }\mathrm{\Sigma }}}\frac{1-\lambda \omega }{E}·\end{array}$$(18)

We transform these velocities from Boyer-Lindquist coordinates into the locally non-rotating frame (LNRF, see Bardeen et al. 1972) for use by the YNOGKM code:

$$\begin{array}{cc}& v_{\mathrm{LNRF}}^{\mathrm{r}}=\sqrt{\frac{A_{\mathrm{r}}}{{(r^2-2a+a^2)}^2}}{v}_{\mathrm{BL}}^{\mathrm{r}},\hfill \end{array}$$(19)

$$\begin{array}{cc}& v_{\mathrm{LNRF}}^{\theta }=\sqrt{\frac{A_{\mathrm{r}}}{(r^2-2a+a^2)}}{v}_{\mathrm{BL}}^{\theta },\hfill \end{array}$$(20)

$$\begin{array}{cc}& v_{\mathrm{LNRF}}^{\varphi }=\sqrt{\frac{(1-{\mu }^2)A_{\mathrm{r}}^2}{(r^2-2a+a^2){(r^2+a^2{\mu }^2)}^2}}(v_{\mathrm{BL}}^{\varphi }-\frac{2ar}{A_{\mathrm{r}}}),\hfill \end{array}$$(21)

where

$$\begin{array}{c}\hfill A_{\mathrm{r}}=r^4+a^4{\mu }^2+a^{2}r(2+r-2{\mu }^2+r{\mu }^2).\end{array}$$(22)

The result of the YNOGKM and GEOKERR codes is a set of position 4-vectors corresponding to the trajectories of the test particle as well as photons propagating to the observer. The intersections of these trajectories correspond to a point in spacetime at which a photon is emitted by the test particle that reaches the observer. To simplify the task of finding these intersections, we assign the hotspot a Gaussian emissivity profile with a size s (set to a fiducial value of s = 0.5M):

$$\begin{array}{c}\hfill ϵ\propto e^{\frac{-{\mathit{x}}^2}{2s^2}},\end{array}$$(23)

where x is the distance from the center of the hotspot.

Given the positions and velocities of a particle as well as the constants of motion of the photon trajectory, we can calculate the Doppler shift of the intersection as

$$\begin{array}{c}\hfill D=\mathit{u}·\mathit{v}\end{array}$$(24)

In terms of the metric elements defined above, we write this Doppler shift as

$$\begin{array}{c}\hfill D=\gamma \sqrt{\frac{A}{\mathrm{\Delta }{\rho }^2}}(1-\lambda \mathrm{\Omega }\pm v_{\mathrm{LNRF}}^{\mathrm{r}}\sqrt{\frac{\mathcal{R}}{A}}\pm v_{\mathrm{LNRF}}^{\theta }\sqrt{\frac{\mathrm{\Delta }\mathrm{\Theta }}{A}}),\end{array}$$(25)

where $\gamma =1/\sqrt{1-{\mathit{v}}^2/c^2}$ is the usual relativistic factor,

With our assumption that the emission is optically thin, we can now find the flux observed at the image plane by integrating the emissivity scaled by the Doppler term backwards along the photon trajectories calculated from the GEOKERR code:

$$\begin{array}{c}\hfill F={\displaystyle \int }{D}^2ϵ\mathrm{d}s,\end{array}$$(26)

The time of arrival of the photons at the image plane is the sum of orbital time and the light-travel time from the hotspot to the observer.

The end result of our simulation is a series of snapshots at increasing observer timesteps of the region immediately around the central black hole. For each image, we can then calculate a centroid position,

$$\begin{array}{c}\hfill \mathit{c}={\displaystyle \int \int }\mathit{x}F(x,y)\mathrm{d}x\mathrm{d}y,\end{array}$$(27)

and a flux,

$$\begin{array}{c}\hfill F={\displaystyle \int \int }F(x,y)\mathrm{d}x\mathrm{d}y,\end{array}$$(28)

Figure 1 shows an example series of images of a hotspot passing behind the black hole. As the spot moves behind the black hole, the photon trajectories are strongly lensed in the direction of the observer, resulting in the corresponding peak in the lightcurve. Simultaneously, the centroid of the emission moves towards the location of the black hole.

Fig. 1. Sample image of a hotspot at different points in the orbit. Top row: image of the hotspot before, during, and after a strong lensing event. Top-right plot: image averaged over the whole orbit. Lower left panel: lightcurve produced over several orbits, lower right panel: motion of the centroid of emission. This figure was produced with a hotspot in a circular orbit at R = 5M around a black hole with spin a = 0.9. The spin axis of the black hole is inclined i = 15° from the observer’s line of sight and the relative inclination between the spin axis and the orbital plane is θ = 50.54°.

4. Orbit fitting

The GRAVITY instrument has found evidence of orbital motion around Sgr A*. Here we describe the process used to fit the observed motion with the NERO code to obtain constraints on the properties of the emitting region (Gravity Collaboration 2018a).

We construct a grid of model centroid motions across our parameter space. For the purposes of this work, we limit ourselves to fitting only for the orbital radius and inclination to the observer. Therefore, we calculate a range of models between R = 4 and R = 14 in increments of ΔR = 0.5 across all inclinations (1° < i <  180°, Δi = 5). Inclinations between 0° and 90° correspond to counterclockwise orbits on the sky, while inclinations above 90° correspond to clockwise orbits. At each point in this grid, we further vary the position angle of the orbit on the sky, the orbital phase at the start of the observation, and the x- and y-position of the central mass. For each parameter combination, we calculate the χ² difference between the model and the astrometric data. The location of the minimum χ² value determines the best fit to the observation.

Figure 2 shows the best-fit models for each of the three flares observed during the summer of 2018. The reduced χ² for each fit is 1.6, 4.5, and 2.3 for July 22, July 28, and May 27, respectively. In each case, the best-fit model lies mostly or entirely within the astrometric data points. This is due to our assumption that the orbiting hotspot follows a geodesic of motion, meaning that its orbital speed is determined by its distance from the black hole. At face value, the data from all three nights show motion that is faster than the expected orbital speed or equivalently, a radius that is larger than expected.

Fig. 2. Data and NERO fits to the flares from May 27, July 22, and July 28, 2018. Data points correspond to the centroid position observed by GRAVITY over the course of the flare. The solid lines show the orbits that best fits each of three flares. Colors indicate the time since the start of the flare.

One effect that partially accounts for this tension is related to the scatter in the observed centroid points. Since the distance from the center of the orbit is a positive definite quantity, any scatter in the measurements will tend to increase the observed radius by a factor of $\sqrt{1+{\sigma }^2/R^2}$. In the case of an orbit at R = 7M observed with 30 μas error bars, this effect implies a 10 μas bias in the observed radius. Other possible explanations are that the motion of the hotspot is not in the plane of the black hole but rather moving outward towards the observer (see Sect. 6 below) or that other forces, such as a pressure gradient or the magnetic field, have a significant effect on the motion of the hotspot.

We marginalize over the various parameters to find the 1-, 2-, and 3-σ allowed regions, corresponding to the area containing 39%, 86%, and 98.8% confidence regions. Given the large values of the χ² of the best-fit models, we scale the uncertainties of the model parameters such that the reduced χ² of the best-fit model is 1. The first three panels of Fig. 3 show the constraints in the R − i plane for all three flares. The clockwise motion during the flares constrains the inclination of the orbit to be between ∼120° −175° for all but the May 27 flare. In the latter case, the best-fit inclination is also in this range (140°) but the allowed region extends over a larger fraction of the parameter space due to the poor astrometric precision of the data. Unless otherwise noted, we include a sin i prior on the inclination of the orbit to account for the geometric bias favoring edge-on orientations. We find that the astrometric fit favors models with moderate to high inclinations and radii between 7 and 12M. The allowed parameter spaces of all three flares overlap, indicating that a single model is consistent with all the flares.

Fig. 3. Confidence contours corresponding to 1-, 2-, and 3σ limits on the orbital radius and inclination for all three flares. The last panel shows the combined constraints, assuming the radius, inclination, and position angle are the same for all three flares.

The contours shown in Fig. 3 are slightly less constraining than those found in Gravity Collaboration (2018a). This difference is due to the additional degrees of freedom allowed in the fit corresponding to the uncertainty of the position of the black hole. In previous fits, the center of the orbit was placed at the average of all the astrometric points over the course of the flare. Figure 4 illustrates the degeneracy between the x-position of the central mass and the radius of the orbit. This degeneracy indicates that the flare astrometry allows for models in which the flare represents only a small fraction of an orbit with a large radius, where the black hole is correspondingly offset from the center of emission. However, these models can be ruled out based on our knowledge of the position of the center of mass from the measured orbits of the S-stars. In particular, the orbit of the star S2 constrains the position of the center of mass to within ∼50 μas (Gravity Collaboration 2018b).

Fig. 4. Confidence contours showing the allowed values for the orbital radius and the x-offset of the center of the black hole for the July 22, 2018 flare. An offset of zero corresponds to the black hole lying at the mean location of emission averaged over the entire flare.

The fits presented here only make use of the measurements of the centroid position of the emission throughout the flare. Two additional sources of information are available that can further constrain the parameters of the fit. The first is the absence of clear Doppler beaming during a portion of the flare. For orbits that are close to edge-on to the observer, the relativistic speed at which the hotspot moves should lead to significant brightening on the approaching side and dimming on the receding side. These effects are, of course, convolved with the intrinsic light curve of the flare, making it difficult to place precise constraints on the inclination. However, orbits with inclinations between 40° and 140° would lead to Doppler beaming comparable to the overall brightening of Sgr A* during the flare, making these scenarios unlikely (Gravity Collaboration 2018a).

Since we do not know the underlying mechanism for generating flares, it is impossible for us to construct light curves that can be compared directly to the data. However, the NERO code does calculate the Doppler boosting, lensing, and other relativistic effects of the motion of the hotspot on the observed flux. By combining these predictions with the observed lightcurves, we can find the intrinsic flare lightcurve as generated at the source. Figure 5 shows both the observed and inferred intrinsic lightcurves for each of the three flares discussed in this paper.

Fig. 5. Observed and inferred intrinsic lightcurves for the three GRAVITY flares. The flux is measured relative to the nearby star S2. Each dotted line corresponds to the lightcurve during the respective flare as observed by GRAVITY. The solid line shows the intrinsic lightcurve once all relativistic effects predicted from the best-fit orbital model have been remove.

A further constraint arises from the observed change in the polarization angle during the course of the flare. The polarization angle changes on a timescale similar to the orbital timescale at the radii favored by the astrometric data. A detailed analysis of the constraints arising from the polarization will be addressed in a forthcoming paper (Jiménez-Rosales et al., in prep.).

If we make the assumption that all of the flares share the same parameters, we can obtain a combined constraint for all of the flares. In order to calculate the combined constraint, we allow the orbital phase at the start of each observation as well as the x- and y-position of the central mass to vary between observations. We hold the radius, inclination, and position angle of the model constant across all three observations. The last panel of Fig. 3 shows the constraint on the R − i parameter space from all three flares combined. The preferred radius is between 8 and 10M, and the inclination lies between i = 120° and i = 150°.

Figure 6 shows the best-fit region from the astrometry combined with the constraints from the lightcurve and the polarization. The contrast constraint is based on a conservative upper limit on the amount of Doppler boosting that is consistent with the observed lightcurve during the flare. The two hatched regions correspond to the radius range determined from the rate of change of the polarization angle and the behavior of the polarization angle in the Q−U plane. We note that this constraint is based on the assumption that the magnetic field is vertically oriented in the emitting region of the accretion flow (see Gravity Collaboration 2018a as well as Jiménez-Rosales et al., in prep. for details).

Fig. 6. Combined constraints together with information from the polarization and the Doppler beaming. The blue astrometric constraints are as shown in Fig. 3. The green shaded region corresponds to models where the maximum contrast of the lightcurve is less than 2.5. The hatched polarization period constraint arises from the assumption that the change in polarization angle during the flare is caused by the motion of the hotspot through the magnetic field. The polarization loop constraint is based on a model with a vertical magnetic field and the requirement that the polarization angle matches the behavior observed during the flare.

We also show the constraints on the orientation of the orbiting hotspots in Fig. 7. Interestingly, the allowed region of orbital orientations is broadly consistent with those of the stars in the clockwise stellar disk, as marked by the black cross in the figure (Gillessen et al. 2017, see also Bartko et al. 2009; Yelda et al. 2014). This clockwise disk consists of O-stars with significant stellar winds (Paumard et al. 2006). These winds have been conjectured as a possible source of matter to fuel the accretion onto the central black hole. If this is the case, the accretion flow is expected to roughly align with the stars in the clockwise disk (Ressler et al. 2018; Calderón et al. 2020). While our constraints cannot definitively confirm this picture, improved constraints from future observations may be able to pinpoint whether the accreting material originates in the winds from massive stars in the clockwise disk.

Fig. 7. Combined constraints on the orientation of the hotspot orbit from all three flares. The black marker indicates the approximate orientation of the clockwise disk of stars.

5. Orbital shear

If the circular centroid motion during infrared flares is indeed due to the orbit of a hotspot, the hotspot is expected to shear over the course of its orbit due the different orbital velocity at different radii. Although the shape of the hotspot is unresolved, this shearing would manifest itself as an inward spiral of the centroid over the course of the flare. By incorporating the orbital shear into our model of the centroid motion and comparing this prediction to the observed astrometry, we can, in principle, obtain an upper limit on the size of the emitting region.

We allow for the possibility of shearing by modifying our model to account for the differential orbital velocity at different points in the emitting region. We define a parameter s corresponding to the radial extent of the hotspot. For a given orbital radius R, we calculate a series of particle trajectories with starting positions ranging from R − s to R + s. We then sum the emission from a Gaussian hotspot on each of these geodesics to find the total predicted lightcurve and centroid track.

We use the same procedure described above to fit the resulting sheared-spot models to the data. In this case, we fix the inclination and restrict ourselves to orbital radii between 5M and 12.5M. We let the diameter of the emitting region vary between 1.0M and 6.5M. Given the poor astrometric precision of the flare observations on May 27 and July 28, we have restricted this analysis to the July 22 flare. Figure 8 shows contours of constant χ² corresponding to 1-, and 2σ limits on the orbital radius and spot extent. As the size of the spot increases, the quality of the fit decreases. Spots larger than 5M in diameter can be ruled out at the 1-σ level, although larger spots are still within the 2-σ limits from the current observations.

Fig. 8. 1- and 2-σ constraints on the orbital radius and spot diameter for the July 22 flare.

6. Out-of-plane motion

Although the detection of circular motion during flares would seem to favor the interpretation that flares originate in accreting material, it is, in principle, also consistent with material launched at the base of a jet. This material would retain some angular momentum from the accretion flow and therefore spiral outward from the jet-launching region.

To investigate this possibility, we simulate the motion of material launched vertically from an accretion flow. In particular, we modify our existing models of orbital motion by adding a constant vertical component v_z to the velocity of the hotspot. We then renormalize the four-velocity of the hotspot to conserve the relation u ⋅ u = −1. The result is a hotspot on a helical track originating in the equatorial plane. Since the volume of the region we simulate is limited for computational reasons, we have restricted ourselves to scenarios where v_z <  0.3c to prevent the hotspot from leaving the bounds of our simulation during a single orbit.

We again use the procedure described above to fit the parameters of the model to the July 22 flare. If the axis of a jet is significantly inclined with respect to the observer, the hotspot should display a systematic linear motion over time. The relatively circular motion measured by GRAVITY rules out this scenario. Therefore, we restrict the inclination to be greater than 160°. Given the face-on inclination inherently required by this scenario, we do not weight the probability distribution by sin i as above. We allow the vertical velocity to vary between 0.05c and 0.3c.

Figure 9 shows the allowed region of the parameter space for the July 22 flare. As expected, the combination of high inclination and large vertical velocity has been ruled out. For a sufficiently low inclination, however, the primary effect of the vertical velocity is to change the observed orbital period due to the propagation time of the photons and the decreasing distance to the hotspot. For face-on orbits, the observed centroid motion is consistent with a significant vertical velocity. In fact, adding a vertical velocity component marginally improves the quality of the fit, with the best-fit parameters being R = 9.0M, i = 175°, v_z = 0.15c.

Fig. 9. Constant χ2 contours corresponding the 1-, and 2σ regions of the vertical velocity and inclination for the July 22 flare. We note that in this case we use a flat prior for the inclination of the orbit.

7. Prospects for improved constraints

There are two possible avenues for improving the constraints on the flare properties from astrometric measurement. The first consists of continued observations of flares at the current level of precision. As shown in the last panel of Fig. 3, co-adding several flares reduces the statistical error on the inferred parameters. Additional flare observations will further tighten the constraints. However, these improvements are contingent on the assumption that some parameters of the flares (i.e., the radius, inclination, and position angle) are consistent between multiple flares.

The second avenue for improved constraints stems from increasing the astrometric precision of measurements of individual flares. To simulate this scenario, we have generated simulated data sets at various levels of precision and used the NERO code to find the best-fit model and resulting constraints on the flare properties. Figure 10 shows the confidence contours for three such simulations. In each case, the underlying model is a hotspot orbiting at a radius of R = 8M and an inclination of i = 155°. Each simulated data set consists of ten data points chosen at time intervals corresponding to the observed times of the July 22nd flare. The x- and y-data points are drawn from a normal distribution centered on the model value with a standard deviation corresponding to an error of 10, 20, and 30 μas, respectively.

Fig. 10. Constraints from simulated data with improved error bars. In each case, the underlying model is based on a hotspot orbiting at R = 8.0M with an inclination of i = 155°, as marked by the red cross in each panel. The astrometric error is 10 μas, 20 μas, and 30 μas, respectively.

We then apply the same fitting procedure to our simulated data points as we use on the real data. The constraints from these fits are shown in Fig. 10. For small error bars, the radius and inclination are tightly constrained. We correctly recover the parameters of the underlying simulation, as marked in the figure. For a 30 μas error, which most closely corresponds to the current GRAVITY precision, the simulated data are consistent with radii between 8M and 14M. The inclination is only loosely constrained.

8. Summary and discussion

The NERO code combines the ability to numerically calculate relativistic orbits near a black hole with efficient ray-tracing between the vicinity of the black hole and a distant observer. The result is a hotspot model for fitting both astrometric and light-curve data with the ability to account for the full general relativistic treatment of the spacetime near the black hole. It also takes into account a variety of physical scenarios, including hotspot orbits that are inclined to the spin axis, eccentric orbits, out-of-plane motion, and orbital shear from hotspots with a significant radial extent.

We apply this model to the three flares observed during the summer of 2018 with the GRAVITY instrument at the VLT. We find that all three flares are well fit by hotspots on circular orbits with radii between 7 and 10M and inclinations greater than ∼90°. A single model with R = 9M is consistent with all three observed flares at a range of inclinations.

While the observed centroid motions are consistent with circular orbital motion, the size of the observed motion on the sky seems to be systematically larger than what has been predicted by our models. Equivalently, the period of the motion is smaller than predicted by the observed size of the centroid track. One possible explanation for this discrepancy is that the hotspot is moving towards the observer at a significant fraction of the speed of light, causing the observed orbit to appear faster than it really is. We find that models with an additional vertical velocity component are consistent with the observations, but they are not significantly preferred over orbits in the equatorial plane.

We also investigate the effects of orbital shear on the predicted motion of the centroid. For a hotspot with a size comparable to its orbital radius, the differential velocity at different radii should lead to a centroid track that spirals inward over the course of the orbit. We find that smaller hotspots provide a better fit to the observations, but hotspots with diameters as large as 6M are allowed by the current astrometric precision.

Additional constraints on the viability of a model of orbital motion comes from the observed light curve and polarization of the flare. A hotspot orbiting at a high inclination to the observer is expected to show significant variation in brightness due to the difference in Doppler beaming between the approaching and receding parts of the orbit. However, this constraint is significantly complicated by the unknown emission mechanism of infrared flares, which leads to a large uncertainty in the expected light curve. In addition to astrometric motion during infrared flares, the polarization angle of the emission also varies on a timescale consistent with the orbital period of a hotspot. Using these observations to constrain the orbit requires a model of the magnetic field near the black hole and will be presented in a future work.

At the current level of precision, astrometric measurements during infrared flares provide strong evidence that matter is in orbital motion around the central black hole. The NERO code allows for relativistic orbits to be fit to the astrometric data in order to constrain the parameters of the orbits. The combination of multiple flare observations can reduce the uncertainties on the parameters of the orbiting material. Further observations of greater precision will be able to improve these constraints even further.

References

All Figures

Fig. 1. Sample image of a hotspot at different points in the orbit. Top row: image of the hotspot before, during, and after a strong lensing event. Top-right plot: image averaged over the whole orbit. Lower left panel: lightcurve produced over several orbits, lower right panel: motion of the centroid of emission. This figure was produced with a hotspot in a circular orbit at R = 5M around a black hole with spin a = 0.9. The spin axis of the black hole is inclined i = 15° from the observer’s line of sight and the relative inclination between the spin axis and the orbital plane is θ = 50.54°.

In the text

	Fig. 2. Data and NERO fits to the flares from May 27, July 22, and July 28, 2018. Data points correspond to the centroid position observed by GRAVITY over the course of the flare. The solid lines show the orbits that best fits each of three flares. Colors indicate the time since the start of the flare.
In the text

	Fig. 3. Confidence contours corresponding to 1-, 2-, and 3σ limits on the orbital radius and inclination for all three flares. The last panel shows the combined constraints, assuming the radius, inclination, and position angle are the same for all three flares.
In the text

	Fig. 4. Confidence contours showing the allowed values for the orbital radius and the x-offset of the center of the black hole for the July 22, 2018 flare. An offset of zero corresponds to the black hole lying at the mean location of emission averaged over the entire flare.
In the text

	Fig. 5. Observed and inferred intrinsic lightcurves for the three GRAVITY flares. The flux is measured relative to the nearby star S2. Each dotted line corresponds to the lightcurve during the respective flare as observed by GRAVITY. The solid line shows the intrinsic lightcurve once all relativistic effects predicted from the best-fit orbital model have been remove.
In the text

Fig. 6. Combined constraints together with information from the polarization and the Doppler beaming. The blue astrometric constraints are as shown in Fig. 3. The green shaded region corresponds to models where the maximum contrast of the lightcurve is less than 2.5. The hatched polarization period constraint arises from the assumption that the change in polarization angle during the flare is caused by the motion of the hotspot through the magnetic field. The polarization loop constraint is based on a model with a vertical magnetic field and the requirement that the polarization angle matches the behavior observed during the flare.

In the text

	Fig. 7. Combined constraints on the orientation of the hotspot orbit from all three flares. The black marker indicates the approximate orientation of the clockwise disk of stars.
In the text

	Fig. 8. 1- and 2-σ constraints on the orbital radius and spot diameter for the July 22 flare.
In the text

	Fig. 9. Constant χ2 contours corresponding the 1-, and 2σ regions of the vertical velocity and inclination for the July 22 flare. We note that in this case we use a flat prior for the inclination of the orbit.
In the text

	Fig. 10. Constraints from simulated data with improved error bars. In each case, the underlying model is based on a hotspot orbiting at R = 8.0M with an inclination of i = 155°, as marked by the red cross in each panel. The astrometric error is 10 μas, 20 μas, and 30 μas, respectively.
In the text