3 there are two sets of first order constraints, obtained from Eqs. (4) and defined as which were introduced when performing the reduction to first order [2, 31]. On the other hand, there are the Hamiltonian and momentum constraints, that in the Generalized Harmonic formulation show up in terms of a four-vector Z$_{a}$ , which is defined as It can be shown that the Hamiltonian and momentum constraints are satisfied if Z a = ∂$_{t}$Z a = 0 [32]. In order to dynamically control the violation of the constraints, we have included certain terms proportional to these con- straints (8)–(9). These additional terms depend on free parameters σ$_{0}$ and σ$_{1}$ , allowing one to dynamically damp constraint—including the Hamiltonian, momentum, and first-order constraints (8)—violating modes on a time scale proportional to − σ$_{i}$ ([31, 33]). We evolve the gravitational field equations shown in Eqs. (5-7). These equations rely on the computation of the 4-dimensional Christoffel symbols from the metric g$_{ab}$ While we evolve the D$_{iab}$ functions, the set D$_{0}$$_{ab}$ are not evolved, but are calculated from evolved quantities as This description suffices to explain the gravitational evolution, and the following section describes the evolu- tion of the matter. However, we note here that the MHD equations are written in the standard 3+1 decomposi- tion of spacetime and thus require the spatial metric h$_{ij}$ , the lapse α , shift β $^{i}$, and ADM extrinsic curvature, K$_{ij}$ . These quantities can be written in terms of our evolved fields using Conversely, the Hamiltonian and momentum constraints are usually written in terms of spatial derivatives of the metric D$_{kij}$ and the extrinsic curvature K$_{ij}$ . In fact, we use these 3+1 quantities (and similar expressions for their derivatives) to calculate the residuals of the Hamiltonian and momentum constraints expressed in their standard form. We now briefly introduce the perfect fluid equations. Additional information can be found in our previous work [27, 28] as well as in topical review articles [34, 35]. leads to the evolution equation of the rest mass density ρ$_{o}$ . Closure of the equations is achieved by introducing an equation of state (EOS) relating the pressure with provides 4 evolution equations for the fluid variables, namely the velocity and the internal energy. Conserva- tion of the baryon number Conservation of the stress energy tensor in Eq. (13), where $^{∗}$F ab ≡ ϵ $^{abcd}$F$_{cd}$ / 2 is the dual of the Maxwell ten- sor and we have introduced a real scalar field Ψ to control the divergence constraint. This technique is known as di- vergence cleaning [36] and allows for a convenient way to control the constraint violation by inducing a damped wave equation for the scalar field Ψ. The other Maxwell equation, in the ideal MHD limit, only gives the definition for the current density, since the electric field is given in terms of the velocity of the fluid and the magnetic field, that is, The evolution of the magnetized fluid is described by different sets of conservation laws. The magnetic field, in the ideal MHD limit, follows the Maxwell equation where E a and B a are the electric and magnetic fields measured by a “normal” observer n $^{a}$. Consequently, both fields are purely spatial, i.e., E $^{a}$n$_{a}$ = B $^{a}$n$_{a}$ = 0. The Maxwell tensor F$_{ab}$ can be written as where W is the Lorentz factor between the fluid frame and the fiducial ADM observers. and construct the standard spatial coordinate velocity of the fluid v i as The fluid is described by rest mass density ρ$_{o}$ , the specific internal energy density ϵ , the isotropic pressure P and the four velocity of the fluid u $^{a}$. With these quantities we can construct the enthalpy The stress-energy tensor for the perfect fluid in the presence of a Maxwell field is given by C iab ≡ ∂$_{i}$g$_{ab}$ − D$_{iab}$ = 0 , ijab ∂$_{i}$D$_{jab}$ ∂$_{j}$D$_{iab}$ = 0 , (8) 2 Z a ≡ − Γ$_{a}$ − H $^{a}$( t, x $^{i}$) . (9) Γ$_{abc}$ = 1 2 ( D$_{bca}$ + D$_{cba}$ − D$_{abc}$ ) . (10) D$_{0}$$_{ab}$ = − αQ$_{ab}$ + β $^{k}$D$_{kab}$ . (11) h$_{ij}$ = g$_{ij}$ , α = √ − 1 /g 00 , β i = γ $^{ij}$g$_{0}$$_{j}$ , K$_{ij}$ = 1 2 Q$_{ij}$ + 1 α ( D$_{(}$$_{ij}$$_{)0}$ − β $^{k}$D$_{(}$$_{ij}$$_{)}$$_{k}$ ) . (12) ∇ $_{a}$( ρ$_{o}$u $^{a}$) = 0 (20) ∇ $_{a}$T ab = 0 , (19) E$_{i}$ = − ϵ$_{ijk}$v $^{k}$Bk . (18) ∇ $_{µ}$( $^{∗}$F µν + g $^{µν}$Ψ) = κn $^{ν}$Ψ , (17) F ab = n $^{a}$Eb − n $^{b}$Ea + ϵ abcd B$_{c}$ n$_{d}$, (16) W ≡ − n $^{a}$u$_{a}$ , v i ≡ 1 W h $^{i}$$_{j}$u$^{j}$, (15) h$_{e}$ = ρ$_{o}$ + ρ$_{o}$ϵ + P, (14) T$_{ab}$ = [ ρ$_{o}$ (1 + ϵ ) + P ] u$_{a}$u$_{b}$ + P g$_{ab}$ + F$_{a}$ $^{c}$F$_{bc}$ − 1 4 g$_{ab}$F $^{cd}$F$_{cd}$ . (13) B. MHD equations