  In this paper we develop a new well-balanced discontinuous Galerkin (DG)
finite element scheme with subcell finite volume (FV) limiter for the numerical
solution of the Einstein--Euler equations of general relativity based on a
first order hyperbolic reformulation of the Z4 formalism. The first order Z4
system, which is composed of 59 equations, is analyzed and proven to be
strongly hyperbolic for a general metric. The well-balancing is achieved for
arbitrary but a priori known equilibria by subtracting a discrete version of
the equilibrium solution from the discretized time-dependent PDE system.
Special care has also been taken in the design of the numerical viscosity so
that the well-balancing property is achieved. As for the treatment of low
density matter, e.g. when simulating massive compact objects like neutron stars
surrounded by vacuum, we have introduced a new filter in the conversion from
the conserved to the primitive variables, preventing superluminal velocities
when the density drops below a certain threshold, and being potentially also
very useful for the numerical investigation of highly rarefied relativistic
astrophysical flows. Thanks to these improvements, all standard tests of
numerical relativity are successfully reproduced, reaching three achievements:
(i) we are able to obtain stable long term simulations of stationary black
holes, including Kerr black holes with extreme spin, which after an initial
perturbation return perfectly back to the equilibrium solution up to machine
precision; (ii) a (standard) TOV star under perturbation is evolved in pure
vacuum ($\rho$=$p$=0) up to t=1000 with no need to introduce any artificial
atmosphere around the star; and, (iii) we solve the head on collision of two
punctures black holes, that was previously considered un--tractable within the
Z4 formalism.
