### 4.3 The magnetic-field divergence-free constraint

The numerical advantage of using the integral form of the conservation law, Equation (76), for the state vector variables does not apply in a straightforward way in the case of the GRMHD equations for the magnetic field components. While the induction equation is part of the hyperbolic system, there remains an additional constraint equation to be fulfilled, namely the magnetic field divergence-free constraint. Indeed, there is no guarantee that such divergence is conserved numerically when updating the magnetic field using the numerical procedure of Equation (76) for the rest of the components of the state vector (continuity, momentum, and energy equations).

Among the methods designed to preserve the divergence of the magnetic field (see [401] for a review), the constrained transport method (CT hereafter) designed by [119] and first extended to HRSC methods by [343] (see also [229] for a recent discussion), stands as the preferred choice for most groups. We note, however, that for AMR codes based on unstructured grids and multiple coordinate systems (such as those of [20286]) there are other schemes, such as projection methods or hyperbolic divergence cleaning, which appear more suitable than the CT method to enforce the magnetic field divergence-free constraint. The projection method involves solving an elliptic equation for a corrected magnetic field projected onto the subspace of zero divergence solutions by a linear operator. More precisely the magnetic field is decomposed into a curl and a gradient as

whose divergence leads to an elliptic (Poisson) equation , which can be solved for the scalar quantity . The magnetic field is next corrected according to . We note that an alternative projection scheme can be implemented by taking the curl of Equation (95) and solving for the vector potential .

Correspondingly, the basic idea of the hyperbolic divergence cleaning approach is to introduce an additional scalar field , which is coupled to the magnetic field by a gradient term in the induction equation. The scalar field is calculated by adding an additional constraint hyperbolic equation given by

Therefore, divergence errors are propagated off the grid in a wave-like manner with characteristic speed . Further details on these two methods are given, e.g., in [401].

The approach followed by the CT scheme to ensure the solenoidal condition of the magnetic field is based on the use of Stokes theorem after the integration of the induction equation on surfaces of constant and , . Let us write the induction equation as

where we have defined the densitized magnetic field vector and .

Following [24] we can obtain a discretized version of Equation (97) as follows. At a given time, each numerical cell is bounded by six two-surfaces. Let us consider the two-surface , defined by  and , and the remaining two coordinates spanning the intervals from to , and from to . The magnetic flux through this two-surface is given by

Integrating Equation (97) on the two-surface , and applying Stokes theorem to the right hand side leads to the following equation

with being the boundary of the corresponding two-surface. This equation can be integrated to give
where the caret denotes that quantities are calculated at the edges of the numerical cells, where they can be discontinuous. At each edge, these quantities can be calculated using the solution of four Riemann problems between the corresponding faces, whose intersection defines the edge. It is important to emphasize, however, that irrespective of the expression used to compute , the method to advance the magnetic fluxes at the faces of the numerical cells satisfies, by construction, the divergence constraint (see, e.g., [24] for details). Indeed, the total magnetic flux through the closed surface formed by the boundary of the numerical cells and enclosing a volume can be calculated as
after applying the Gauss theorem and divergence constraint. This equation implies that no magnetic flux source exists inside the volume and, therefore, the magnetic flux is a conserved quantity as . As a result, every numerical scheme constructed on the basis of the CT scheme will conserve magnetic flux up to machine accuracy. If one is able to generate initial conditions, which satisfy the divergence constraint, i.e., with across the boundary of each numerical cell, then such a constraint will be preserved during the numerical evolution.

Any of the Riemann solvers and flux formulae discussed in Section 4.1.2 can be thus used to calculate the quantities needed to advance in time the magnetic fluxes following Equation (100). At each edge of the numerical cell, is written as an average of the numerical fluxes calculated at the interfaces between the faces whose intersection define the edge. Let us consider, for illustrative purposes, . If the indices denote the center of a numerical cell, an edge is defined by the indices . By definition, . Since and , we can express in terms of the fluxes as follows

where () refers to the numerical flux in the () direction corresponding to the equation for () and multiplied by to account for the correct definition of . Analogous expressions for and , can be derived in a similar way. We refer the interested reader to [401] for additional details on the constrained transport method.