def METHOD_NAME(self):
    """
    (`~astropy.units.Quantity`) (x, y, z) array of the Alfvén speed at
    every point in the domain.
    """
    B = self.magnetic_field
    rho = self.density
    return np.sqrt(np.sum(B * B, axis=0) / (mu0 * rho))