Cosmic magnetic fields with MASCLET: an application to galaxy clusters
Vicent Quilis José Ma. Martí Susana Planelles
astro-ph.CO astro-ph.IM




 A hybrid quantum-classical approach to mitigating measurement errors 
    Hyeokjea Kwon^* 

        and Joonwoo Bae^†
 H. Kwon and J. Bae are with the School of Electrical Engineering, Korea Advanced Institute of Science and Technology, Daejeon, Republic of (South) Korea.


E-mail: ^*hyukjaekwon@kaist.ac.kr, ^†joonwoo.bae@kaist.ac.kr



    December 30, 2023
============================================================================================================================================================================================================================================================================





We describe and test a new version of the adaptive mesh refinement
(AMR) cosmological code MASCLET. The new version of the code includes
all the ingredients of its previous version plus a description of the
evolution of the magnetic field under the approximation of the ideal
magneto-hydrodynamics (MHD). To preserve the divergence-free condition
of MHD, the original divergence cleaning algorithm of <cit.>
is implemented. We present a set of well-known 1D and 2D tests, such
as several shock-tube problems, the fast rotor and the Orszag-Tang vortex. 
 The performance of the code in all the tests is excellent with 
  estimated median relative errors of ∇· B in the 2D
  tests smaller than 5 × 10^-5 for the fast rotor test, and
  5 × 10^-3 for the Orszag-Tang vortex. As an astrophysical application of the
code, we present a simulation of a cosmological box of 40  comoving
Mpc side length in which a primordial uniform comoving magnetic field
of strength 0.1 nG is seeded. The simulation shows how the magnetic
field is channelled along the filaments of gas and is concentrated and
amplified within galaxy clusters.  Comparison with the values
  expected from pure compression reveals an additional amplification
  of the magnetic field caused by turbulence in the central region of
  the cluster.  Values of the order of ∼ 1 μG are obtained in clusters at z∼ 0  with median relative errors 
  of ∇· B below 0.4%. The
implications of a proper description of the dynamics of the magnetic
field and their possible observational counterparts in future
facilities are discussed.
 


magneto-hydrodynamics - methods: numerical - galaxy formation
- large-scale structure of Universe - Cosmology






§ INTRODUCTION



A new era in cosmological observations favoured by forthcoming new
observational facilities, such as the Square-Kilometre
  Array[https://www.skatelescope.org/] <cit.> or
[http://www.the-athena-x-ray-observatory.eu]
<cit.>, and the already existing 
ones like [https://www.almaobservatory.org/es/inicio/]
<cit.>, will produce a huge amount of extra
high-quality data, undoubtedly leading to a deeper knowledge on the
processes carving the cosmic structures in the Universe.   

The formation of galaxies and galaxy clusters would imprint features
on the baryonic component of cosmic structures that eventually  would
be   observable in different wavelengths of the electromagnetic
spectrum. Related to these processes, shock waves, turbulence,
strong gradients, and all sort of emission processes associated to the
baryonic nature of the gas are expected to be observed and quantified
by the new generation of telescopes <cit.>.  

In this context of high-quality data, the models describing the 
formation and evolution of such structures must also include all the 
relevant physical ingredients. Cosmological codes
<cit.>
have vastly improved in the last years by including all sort of
processes associated to the microphysics of the gas <cit.>. 
 Although for many
years it was believed that the role of B was irrelevant in
the formation of galaxies and galaxy clusters, nowadays, there are
many observational and theoretical results pointing towards the
opposite direction <cit.>. In this line,
the description of shock waves, strong gradients and X-ray emission,
among others, would require a proper modeling of the magnetic field
within current cosmological codes  to compare with the observations
delivered by the new generation of telescopes. This
  modelling should allow for the amplification of primordial weak
  magnetic fields via small-scale dynamo during the structure
  formation, and includes other astrophysical sources of magnetic
  seeding at lower redshifts (see <cit.> 
  for a critical discussion on the origin of magnetic fields in
  clusters and a review on small-scale turbulent dynamo).

With this motivation, in recent years the most popular cosmological
codes have developed  updated versions incorporating the evolution of
seed magnetic fields <cit.>. Current
cosmological smoothed particle MHD codes include those developed by 
<cit.> and
<cit.>. <cit.> have developed the MHD extension
of AREPO (<cit.>, a finite volume code using unstructured
moving mesh). These codes have adapted the divergence cleaning
(DC) algorithm of <cit.> to remove the numerical errors which
otherwise would accumulate in non-soleinodal magnetic fields. Dedner's
DC method was also the algorithm implemented in the MHD version  of the
AMR, finite volume code ENZO <cit.>. An alternative
to the DC algorithm is the constrained transport (CT) method
<cit.>, which satisfies the magnetic divergence constraint
by construction. The CT method has been implemented in the AMR,
finite-volume, MHD upgrades of RAMSES <cit.>, ENZO <cit.>
and CHARM <cit.>. AREPO-DG <cit.> implements a
discontinuous Galerkin method for ideal MHD on an Eulerian AMR grid
using the eight-wave method of <cit.>, an alternative DC
method to control the magnetic divergence errors.

This paper presents a new version of an already well-tested and used
cosmological code, MASCLET <cit.>. This new version is
improved by adding the evolution of a magnetic field B in the
context of ideal magneto-hydrodynamics (MHD) in a full cosmological
frame. Contrary to CT techniques, which evolve face-centered
  magnetic fields, the chosen DC method of <cit.> keeps the
  cell-centered discretization of the full scheme (at the cost of
  introduzing an additional variable and the corresponding evolution
  equation), which makes it technically simpler to implement in the
  structure of the AMR code.

The paper is organized as follows. In  Section <ref> the
basic equations to be solved are described. The implementation of  
the numerical techniques to solve previous equations and the details
of the numerical code are presented in Section <ref>. In
Section <ref>, the code is tested by carrying out a battery of
numerical tests. A cosmological application is shown in Section
<ref>, where a cosmological box including a magnetic seed is
evolved. Finally, we discuss our results and offer our conclusions in
Section <ref>. 




§ EQUATIONS






 §.§ The MHD equations in an expanding background


The equations of MHD <cit.> describing the evolution
of a magnetized fluid in a gravitational field in Eulerian coordinates
(t,  r) are the following: 



  * Continuity equation:


    ∂ρ/∂ t + ∇· (ρ u) = 0,



  * Momentum conservation:

    ∂ρ u/∂ t + ∇· (ρ
      u u + p^*  I -  B B) = -ρ∇Φ,



  * Energy conservation:

    ∂ e/∂ t + ∇· ((e + p^*) 
      u - ( B· u) B) = -ρ v·∇Φ,



  * Induction equation:

    ∂ B/∂ t - ∇× ( u× B) = 0,



  * Divergence-free condition:

    ∇· B = 0.



In these equations, ρ is the mass density, u is
the fluid velocity, B is the magnetic field, and quantities e
and p^* stand, respectively, for the total energy density (internal
+ kinetic + magnetic) and total  pressure (thermal + magnetic):


    e = ρϵ + 1/2ρ u^2 +  B^2,



    p^* = p +  B^2/2.


The thermal pressure p is related to the mass density ρ  and
the specific internal energy ϵ by means of an equation of
state of the form p = p(ρ, ϵ). Quantity Φ appearing
in he source terms of the momentum and energy equations stands for the 
gravitational potential.

In the presence of an expanding background, the system of
Eqs. <ref>-<ref> is rewritten in terms of
the comoving coordinates, x≡ r/a(t), where a(t)
is the scale factor and H(t) ≡ȧ/a is the Hubble
constant <cit.>: 



    ∂ρ̃/∂ t + 1/a∇·
    (ρ̃ 
      v) = 0,



    ∂ρ̃ v/∂ t + 1/a∇· (ρ̃
      v v + p̃^*  I - B̃B̃) = -
    Hρ̃ v - ρ̃/a∇ϕ



    ∂Ẽ/∂ t   +   1/a∇· ((Ẽ + p̃^*)  
      v - (B̃· v)B̃) = 
                = - 3H p̃^*
    + H B̃^2 - H ρ̃ v^2 -
    ρ̃/a v·∇ϕ,



    ∂B̃/∂ t - 1/a∇×
    ( v×B̃) = - H/2B̃,



    ∇·B̃ = 0.


In this set of equations, all the spatial differential operators
refer to the comoving coordinates x. 
The overdensity with respect to the background density, ρ__B,
is  ρ̃= ρ/ρ_B. The peculiar velocity of the
fluid  is  v =  u -ȧ  x, being ȧ
  x  the Hubble flow velocity. The quantity Ẽ=E/ρ_B is related to the total energy density, E,
of the magneto-fluid, which includes only the kinetic energy
corresponding to the fluid's peculiar velocity, i.e., E = ρϵ + 1/2ρ v^2 +  B^2. The total pressure
and the magnetic field are  redefined as p̃^*= p^*/ρ_B 
 and B̃ =  B/√(ρ_ B), respectively. 
 
The peculiar potential, ϕ = Φ + a ä x^2/2, satisfies
the Poisson equation 

    Δϕ = 3/2H^2a^2 Ωρ̃

where Ω is the density parameter. 

Equations (<ref>-<ref>) coincide with
those of <cit.> in the case with B̃ = 0.



 §.§ Divergence cleaning


MASCLET uses the divergence cleaning algorithm of <cit.> to
control the magnetic field divergence errors. In the original mixed
hyperbolic/parabolic correction approach of <cit.> (the one
chosen in MASCLET), the divergence constraint of the magnetic field
and the induction equation (in an Eulerian frame,
Eqs. <ref>-<ref>) are coupled by introducing
an additional scalar function ψ in such a way that: i) ∇· B and ψ satisfy the same evolution equation, and ii)
this evolution equation leads to the propagation to the domain
boundaries and the decay of ∇· B (and ψ).

In the case of an expanding background, the original Dedner et al.'s
approach can be applied once the linear B̃-term in
Eq. <ref> is removed by shifting to the comoving
magnetic field B' = a^2 √(ρ_ B) B̃,
and the ψ-terms are introduced in
Eqs. (<ref>-<ref>)  

    ∂ B'/∂ t - 1/a∇× (
      v× B') + ∇ψ = 0,



    ∂ψ/∂ t  + c_h^2 ∇· B' = -
    c_h^2/c_p^2ψ.


Now, direct manipulation of these two equations leads to the telegraph equation for ψ (and ∇· B')


    ∂^2 ψ/∂ t^2 + c_h^2/c_p^2∂ψ/∂ t = c_h^2 Δψ


establishing that divergence errors propagate away from the point
where they are produced at a speed c_h, and damp at a rate given by
c_h^2/c_p^2. Quantities c_h and c_p are dimensional parameters
to be tuned (see next section for details).

Equations (<ref>) and (<ref>)
substitute the original ones (Eqs. <ref> and
<ref>) in the system of equations. Since they are
designed to ensure that ∇· B' (as ψ) verify the
telegraph equation, the errors in ∇· B' propagate
out of the numerical domain and damp as desired.  Since ∇·B̃ = 1/a^2 √(ρ_ B)∇· B' ∝ (∇· B')/√(a),
errors in ∇·B̃ dilute in the expanding
background as they propagate and damp.




§ NUMERICAL IMPLEMENTATION






 §.§ The reference version of  MASCLET code


MASCLET <cit.>, is a cosmological
  multidimensional hydrodynamic and N-body code based on an adaptive 
  mesh refinement scheme (AMR). We address interested readers to
    the original paper for a complete technical
    description of the code and offer here a short summary of its
    basic ingredients. 

The hydro solver is a high-order, finite volume Godunov scheme based on
monotonicity preserving cell-reconstruction routines and approximate
Riemann solvers. Two different cell reconstructions are implemented:
a piecewise linear reconstruction with the MINMOD slope limiter, and
the piecewise parabolic method PPM. Two Riemann solvers are also
implemented: a Roe-type linearized Riemann solver and the HLLE Riemann 
solver. Advance in time is done by TVD preserving second and third
order Runke-Kutta methods following a method of lines.

Dark matter is evolved using a Particle-Mesh scheme with a
Lax-Wendroff temporal integrator. 

Both components, gas and dark matter, evolve coupled by the total
gravitational field which is computed at each time step by solving the
Poisson equation. 

The gain in numerical resolution, both spatial and temporal, is
obtained by means of an AMR scheme. This method
refines the original coarse grid into patches whose cells are half size
of their parent cells. In these new grids, all the relevant physical
quantities are be obtained, either by evolution of the same quantities
from the previous time step, or by tri-linear interpolation from the parent
grid (lower level). This process can be repeated iteratively between
two consecutive levels, producing a whole hierarchy of nested
grids. The criteria to decide which regions of a 
given grid must be refined,  can be configured for every specific
application <cit.>.  

Once the AMR hierarchy is defined, at each level, the hydro and the
dark matter solvers, previously described, can be applied. The
different patches at the same level and the child patches and their
parents are connected through the boundary conditions and the values
used to initialize them.  

In the same manner, a multigrid SOR method is used to solve the
Poisson equation in each patch of the hierarchy.    

The current version of the code also  includes  inverse Compton and
free-free cooling, UV heating, atomic and molecular cooling for the
gas depending on the metallicity,  and a phenomenological description of
star formation. The description of feedback phenomena, both thermal
and kinematic, from stars and active galactic nuclei is also included.  



 §.§ The MHD version of  MASCLET code


The MHD version of MASCLET solves the system formed by
Eqs. (<ref>-<ref>),
(<ref>) and (<ref>)
incorporating all the numerical ingredients of the reference code with
the only evident modifications of the fluxes in
Eqs. (<ref>)-(<ref>), and the addition of
two new equations
(Eqs. <ref>-<ref>).   

These equations are evolved in conservation form with
numerical fluxes which, in general, depend on the eigenstructure of the
system. Following the discussion in <cit.>, in the
one-dimensional case (the one used to obtain the eigenstructure),
the system decouples into a subsystem formed by the equations for
the parallel component of the magnetic field, in our case
B'_∥, and ψ, and the usual 1D-MHD system (with a
factor 1/a in front of the vector of fluxes) for the variables
(ρ̃,  v, B̃_⊥, p̃) with the
magnetic field divergence-free constraint, B̃_∥
constant.  

  Along the x-direction, for sufficiently large c_h, the
corresponding eigenvalues of the full system, ordered in a
nondecreasing sequence, are:  


    λ_1 = -c_h,  λ_2 = (v^x - v_f)/a,  λ_3 = (v^x
      - v_a)/a,


    λ_4 = (v^x - v_s)/a,  λ_5 =
          v^x/a,  λ_6 = (v^x + v_s)/a,


    λ_7 = (v^x
      + v_a)/a,  λ_8 = (v^x + v_f)/a,  λ_9 = c_h  ,


where


    v_a = B̃^x/√(ρ̃)

and

    v_f,s
    = {1/2( c_s^2 + B̃^2/ρ̃±√(( c_s^2 + B̃^2/ρ̃)^2 - 4 c_s^2 v_a^2))}^1/2


are, respectively, the Alfvén speed and the fast and slow
magnetosonic speeds, and c_s is the sound speed.

On the other hand, for any right eigenvector r̃_i =
(r̃^1_ i,..., r̃^7_ i)^T of the original  1D-MHD system  
<cit.>, there is
a right eigenvector of the full MHD system with the divergence
correction (Eqs. <ref>-<ref>,
<ref> and <ref>)
R̃_i = (r̃^1_ i,..., r̃^5_ i, 0,
r̃^6_ i, r̃^7_ i, 0 )^T with respect to the same
eigenvalue λ_i, i = 2,3,...,8. Moreover, two new
eigenvectors appear associated to  eigenvalues λ _1 and
λ_9, respectively, R̃_1 = (0,0,0,0,0,1,0,0,-c_h)
and R̃_9 = (0,0,0,0,0,1,0,0,c_h).


In its present version, MASCLET incorporates a new slope limiter 
<cit.> for the piecewise linear reconstruction, 
and an HLLE Riemann solver adapted to the MHD equations. All the
results shown in this paper have been obtained with the MC
slope limiter. The 1D subsystem for the variables
(ρ̃,  v, B̃_⊥, p̃) is evolved
in time with numerical fluxes calculated with the HLLE Riemann solver 
based on upper and lower bounds, respectively, of the largest and
smallest local propagation speed of fast magnetosonic waves,
λ_8, λ_2. The decoupled subsystem for
B'_∥ and ψ (in charge of the diverence cleaning)
is then evolved in time with numerical fluxes based on a (constant)
characteristic speed c_h chosen to be the maximum speed compatible
with the time step restriction, i.e., fixed to the largest (in
absolute value) propagation speed of fast magnetosonic waves,
λ_2, λ_8, across the grid. Following <cit.>, and
attending to the dimensional nature of c_h and c_p, in several 
preliminary tests we defined c_p = √(Δ x c_h/α), with
α = 0.2, 0.5 and Δ x being the cell size of the finest
grid. The results in the AMR applications were not satisfactory and
numerical instabilities between different levels of resolution arose.
The final choice was to recover Dedner's original prescription
<cit.> as suggested recently by <cit.> in the context on
non-uniform grids, and take c_p = √(0.18 c_h)  (in code units). This choice of
c_h and c_p  is the one used in all the numerical applications
presented in this paper.




§ TESTS




In this Section, we present several classical MHD tests in order to
quantify the performance of the code. In all of them, there is no
expanding background as in the general form of
Eqs. <ref>-<ref>. Consequently,
these tests have been carried out by integrating the aforementioned
equations with a=1, H=0 and ρ_B = 1 (in the absence of any 
gravitational field).  



 §.§ Shock-tube tests



Shock tubes have become standard tests where to prove the ability of 
hydrodynamical codes to describe shocks and contact
discontinuities. In the case of MHD codes, the importance of shock
tube tests is even greater as the solution involves a larger variety
of discontinuity types <cit.>. In this section we
discuss the performance of  MASCLET in reproducing  two
representative MHD shock tube tests from <cit.>, in particular, 
test 2a and test 4a. In both tests, we set the adiabatic index of the 
equation of state to γ = 5/3, and use a one-dimensional grid
with x ∈ [0,1]. The initial discontinuity is placed at x=0.5. 
  
  The tests are performed on a fixed grid with a numerical resolution of
512 zones and a piecewise linear reconstruction with the MC slope limiter. 
Finally, in one-dimensional tests, the divergence free
constraint is fulfilled automatically and the divergence cleaning
algorithm does not operate. The numerical results are plotted as red
circles whereas the analytical solutions (kindly provided by D. Ryu)
appear as filled lines.

The first test is test 2a in <cit.> with left state {ρ_L,
v^x_L, v^y_L, v^z_L, B^y_L, B^z_L, p_L} = { 1.08, 1.2, 0.01, 0.5, 3.6/√(4
  π), 2/√(4 π), 0.95} and right state {ρ_R,
v^x_R, v^y_R, v^z_R, B^y_R, B^z_R, p_R} = { 1, 0, 0, 0, 4/√(4
  π), 2/√(4 π), 1}, and B^x = 2/√(4 π), involving a
three-dimensional magnetic field and flow velocity. The solution at t
= 0.2 is shown in Fig. <ref>. Fast shocks, rotational
discontinuities (where the transverse magnetic field changes its
direction) and slow shocks propagate from each side of the contact
discontinuity. Despite the thinness of some structures, all the
constant states and discontinuities are properly captured   and only
tiny numerical oscillations and overshootings are seen in some
quantities in the post-shock states. 







The second test is test 4a in <cit.> with left state {ρ_L,
v^x_L, v^y_L, v^z_L, B^y_L, B^z_L, p_L} = { 1, 0, 0, 0, 1, 0,
1} and right state {ρ_R, v^x_R, v^y_R, v^z_R, B^y_R, B^z_R,
p_R} = { 0.2, 0, 0, 0, 0, 0, 0.1}, and B^x = 1. The solution at
t = 0.15 is shown in Fig. <ref>. This test is a planar Riemann
problem in which the initial transverse magnetic fields and velocities
are confined in a plane. In this case, the Alfvén waves, which
rotate the fields, do not emerge. From left to right, the solution
produces a fast and a slow rarefaction, a contact discontinuity, a
slow shock and a fast, switch-on shock (where the tranverse
magnetic field turns on behind the shock). As in the previous case,
the code captures the analytical solution correctly with only small
overshootings at discontinuities and at the end points of the
rarefaction waves.










 §.§ The fast rotor



The rotor problem <cit.> aims to test the propagation of strong
torsional Alfvén waves. A dense, rapidly spinning disk rotates
rigidly in a light static ambient medium. Both disk and ambient fluids are
threaded by an initially uniform magnetic field. The rapid rotation of
the rotor causes torsional Alfvén waves to propagate into the
ambient fluid reducing the angular momentum of the rotor and hence its 
rotation velocity. At the same time, the magnetic field lines wrapped
along the surface of the rotor increase the magnetic pressure
compressing the fluid in the rotor and giving it an oblong shape.

The numerical domain for this test consists of a two-dimensional
squared grid of unit length with zero-gradient boundary conditions. A
rotating disk of radius r_0 is embedded in a homogeneous ambient
fluid with a transition layer between r_0 and r_1 (> r_0). In
Cartesian (x,y) coordinates, the initial conditions are given by
{ρ, v^x, v^y} = {10, -ω  y, ω  
x} for r = √(x^2+y^2)≤ r_0, and {ρ, v^x, v^y} =
{1,0,0} for r ≥ r_1. For r ∈ ]r_0, r_1[ a transition layer
is set with {ρ, v^x, v^y} = {1+9f, -f  ω y r_0/r, f
ω x r_0/r} where ω = 20, f = (r_1 - r)/(r_1 - r_0),
and r_0 = 0.1, r_1 = 0.115. The thermal pressure is constant with
p = 1 throughout the grid and the magnetic field is uniform and
aligned with the x axis (B^x = 5/√(4 π)). An ideal gas
equation of state with γ = 1.4 is used. As in the previous
tests, a, H and ρ_B are set to 1, 0, 1, respectively.

Figure <ref> shows the density (top row) and  magnetic pressure
(bottom row) distributions at time t = 0.15 (when the torsional
Alfvén waves are about to reach the grid boundaries) with two
numerical resolutions: a fixed grid with 512 × 512 cells and an
AMR grid with a coarse grid of 128 × 128 cells and two levels of
refinement (refinement factor 2; effective resolution 512 × 512 
cells). The piecewise linear reconstruction with the MC limiter is used.

Both simulations resolve the structure of the flow (high-density
shell at the rotor's oblong surface just behind the high-magnetic
pressure barrier; outgoing torsional Alfvén waves) with the same
detail <cit.>. 

Following, for instance <cit.> and <cit.>, we define a
  dimensionless quantity that could be interpreted as a  
relative error of the divergence of  the magnetic field, 
 |∇· B_i| Δ x_i/| B_i|, where
B_i is the magnetic field in the i cell, ∇·
  B_i is the value of the divergence computed numerically in that
cell, and Δ x_i is the cell width. In order to show the
performance of our numerical scheme, we present Figure <ref> where    
a map of the relative error of the ∇· B for our
best numerical resolution run is shown at 0.15 time units. Largest
errors appear associated to and move with magnetic field
discontinuities. 

Complementary to Fig. <ref> which offers a description of the
distribution and magnitude of the numerical error of the magnetic
field divergence, we present a statistical analysis in
Figure <ref>, where yellow, blue, orange and green lines indicate
the upper bound errors for 25, 50, 75 and 90% of the cells as a
function of time, in the AMR run. Most of the cells (more than 90%)
have divergence errors smaller than 10^-3 along the
simulation. The maximum error of the divergence (black line) remains
stable at about 10%. These errors values are in fair agreement with
those shown by <cit.> for an AMR run using divergence cleaning.
















 §.§ Orszag-Tang vortex


The initial state of the Orszag-Tang vortex <cit.> consists of
two magnetic field loops embedded into a large-scale rotating flow
structure on a two-dimensional box of unit length with periodic
boundary conditions. In the subsequent evolution, complex, small-scale
structures are formed. It has become a standard test to probe the
accuracy of MHD codes when simulating the formation and
interaction of MHD shock waves and the transition to two-dimensional MHD
turbulence. 

In Cartesian (x,y) coordinates, the initial conditions are  {ρ,
v^x, v^y, B^x, B^y, p} = {ρ_0, -v_0 sin (2 π y), v_0
sin (2 π x), -B_0 sin (2 π y), B_0 sin (4 π x), p_0}
with ρ_0 = 25/(36 π), v_0 = 1, B_0= 1/√(4 π), p_0 =
5/(12 π), and an ideal gas equation of state with adiabatic index
γ = 5/3. As in the previous tests, we set {a, H, ρ_B} =
{1,0,1} in the code equations. Figure <ref> shows
the density (top row) and  magnetic pressure (bottom row)
distributions at time t = 0.5 with three numerical resolutions: a
fixed grid of 128 × 128 cells, a fixed grid with 512 ×
512 cells and an AMR grid with a coarse grid of 128 × 128 cells
and two additional levels of refinement with a refinement factor of 2
(for an effective resolution of 512 × 512 cells). The piecewise
linear reconstruction with the MC limiter is used.

As it is clearly seen, the AMR run captures the small structures with
the same degree of detail than the 512 × 512 cells fixed grid
run. Additionally, the agreement between our result and previous
published works  is excellent <cit.>. 

The zero-gradient boundary conditions used, as well as
  the fact that a non-negligible fraction of the 
  cells are only slightly perturbed along the simulation, favour that
  the magnetic field divergence errors remain small in the rotor test. In the
  Orszag-Tang vortex test, on the contrary, all the cells are largely
  perturbed in the course of the evolution. This, together with the
  periodic boundary conditions used in this test (similar to those used
  in standard cosmological simulations), makes the Orszag-Tang vortex
  harder from the point of view of the magnetic divergence
  errors[Whereas the zero-gradient boundary 
    conditions used in the rotor test let the outgoing divergence
    errors driven by the hyperbolic term leave the numerical grid, the
    periodic boundary conditions imposed in the Orszag-Tang vortex
    reintroduce them continuously into the computational domain
    (although damped).] but
  also more meaningful.
  

Figure <ref> displays a map of the relative error of the divergence
of B for the best numerical resolution run. As in the fast rotor
test, the largest numerical errors appear associated with flow
discontinuities and do not accumulate in special locations of the
domain. In order to know the time behaviour of the numerical errors,
Fig. <ref> shows the normalized divergence of the magnetic field
for the AMR run as a function of time. As in the previous test,
yellow, blue, orange and green lines indicate the upper 
bound error for 25, 50, 75 and 90% of the cells along the
evolution. These errors, within 10^-4 and 10^-2, tend to
stabilize towards the end of the simulation. The black line represents the
maximum error of the divergence (of order 1) along time. Compared with
the divergence error in the rotor test, the median error in this test
is about one order of magnitude larger but still reasonably small
(∼ 0.4%). The maximum error compares well with the one obtained by
<cit.> for the same test with an AMR run using divergence
cleaning and is about a factor of ten the one reported by <cit.>
with an SPMHD code using divergence cleaning.


















§ A COSMOLOGICAL APPLICATION: MAGNETIC FIELDS IN GALAXY CLUSTERS




A straightforward application of a cosmological code would be to
simulate the evolution of a cosmic volume where a large number of  
cosmic structures, spanning a huge range in masses and sizes, form and
evolve. Thus, as an example of application of the new version of the code
MASCLET, we present the results for a simulation of a
computational box representing a moderate volume of the Universe. 







The simulation assumes a spatially  flat ΛCDM cosmology, with 
these cosmological  parameters:  Ω_ m = 0.31,
Ω_Λ = Λ/3H_0^2 = 0.69, Ω_ b = 0.048,
h = H_0/(100 km s^-1 Mpc^-1) =0.678,  n_ s = 0.96
and σ_8=0.82. The simulated region is discretised  with
128^3 cubical cells within a cube of comoving  side length 40
Mpc. Employing a  CDM transfer function from <cit.>,
we  set up the initial conditions at  z = 100. We applied a
constrained realization in order to generate a rich galaxy cluster in
the center of the box <cit.>. 

From the initial conditions, evolved until present  time using a low
resolution domain, we select regions satisfying some refining criteria
in order to arrange three  refinement levels (l=1, 2, and 3) for
the AMR scheme. In  these initially refined  levels, the  dark matter
(DM)  component is sampled with  DM  particles 8, 64, and 128 times,
respectively, lighter than  those used  to sample regions in the
coarse grid (l=0). As the evolution proceeds, the total density,
that is baryonic plus dark matter densities, is used to refine regions
on  the different grids by means of a pseudo-Lagrangian approach in
which a cell is flagged as refinable if its density increases in a
factor of eight.  

In the present simulation we use a maximum of seven refinement levels 
(l = 7), allowing for a peak physical spatial resolution of ∼ 3
kpc at z=0. Four different particles species are considered for the DM,
corresponding to those particles on the coarse grid and those
within the three first  levels of refinement.  The best  mass
resolution  is  ∼ 2×  10^6   M_⊙, equivalent to use
1024^3 particles in the whole computational domain. 

The simulation is adiabatic, therefore, no cooling or heating
processes are considered. The gas is described as an ideal fluid with
an adiabatic exponent equal to 5/3. 

We mimic the primordial cosmological magnetic field by seeding a
comoving homogenous magnetic field, B_0, filling the whole
computational domain and orientated along one of the axis of the
computational box. The chosen value of this seed is 0.1 nG, which 
is just below observational constraints imposed by the analysis of  CMB
data <cit.>. 
This choice of the initial magnetic field is somehow arbitrary
 and neglects many 
other possible sources of this initial field. Nevertheless, for the purpose of 
this paper and for the sake of comparison with previous works, we consider
this value reasonable enough. Although a different choice could have an impact in the final value of B if no evidences of dynamo amplification were found.
Some recent
works have also proved that the initial topology of the magnetic field
has no relevant effects on the formation of cosmic structures such as
galaxies and galaxy clusters
<cit.>. 

In terms of the simplicity of the physics involved (only cosmic
expansion, gravity and magnetohydrodynamics), the mass of the formed
cluster and the effective numerical resolution, our simulation can be 
compared with the one with the highest numerical resolution discussed
in <cit.>.  
Figure <ref> shows two thin slices (∼ 8 kpc width) of  the 
gas density (left panel) and the magnetic field strength (right panel)
at z∼ 0. The slices cut the cluster centre and have a size of
four virial radius approximately. The density plot exihibits a
prominent core, as expected in an adiabatic 
simulation, corresponding to a galaxy cluster with a mass of 4×
10^14  M_⊙ and a virial radius R_ Vir = 1.96 Mpc. The
comoving magnetic field intensity, B, shows a more complex structure 
as a consequence of the entanglement of the magnetic field lines,
specially within the virial radius. The similarity of this plot
  (characterized by entangled tongues of alternating high and low
  magnetic field) with the corresponding one in <cit.> (see
  their Figure 2) is remarkable. 

In order to visualize the topology and structure of the magnetic field
and its time evolution, we display the Fig. <ref>. The four
panels present the magnetic field lines in a cubic region centred at
the cluster position at z∼ 0 with a side length of 8 comoving Mpc
(∼ 4 R_ Vir). The lines are colour-coded according to the
value of magnetic field intensity in μG. Four epochs,
corresponding to redshifts z = 2, 1, 0.5, and 0 are displayed.
It is clearly visible how, as the cluster builds up, the magnetic
field is amplified in the region within the virial radius. Based on
the properties of the magnetic field, the cluster can be separated
into two broad regions. Within the virial radius, although the
distribution of magnetic field lines is complex and intricate, the
intensity of B is rather homogenous with values of the order of
μG. By contrast,
at the cluster's outskirts, field lines basically follow the
filaments of gas feeding the cluster, and the value of B in these
regions is remarkably lower compared with values at the inner
parts.   







The radial profile of the magnetic field is presented in
Fig. <ref>, where its mass-weighted comoving intensity is
plotted against cluster-centric distance (continuous red line) in units normalized to the
cluster virial radius, R_ Vir. For the sake of comparison, we
also plot  the magnetic field intensity that
would be expected from a pure compression of the magnetic field
lines (discontinuous red line). The values of B in this case can be computed as
B=B_0(ρ/ρ̅)^2/3, where B_0 = 0.1 nG is the comoving
magnetic field, and ρ and ρ̅ are, respectively, the gas
density and the average gas density in the whole computational box. 

Our results show a B profile which is above the
expected profile for magnetic fields from pure compression of the
gas. The computed magnetic field is in fact 2 to 10 times higher than
the one estimated from pure compression for most of the cluster volume,
r > 0.1 R_ Vir (at the very inner centre of the cluster, r<0.1
R_ Vir, the profile of B presents a more complex
structure with a peak and a central dip which we associate with the 
particular dynamical state of the simulated cluster). This magnetic
field amplification is in contrast with the results for the
adiabatic simulation displayed by <cit.> (see their
Figure 3). In this case, the magnetic field at the inner parts of the
cluster is below the expected values from B ∝ρ^2/3. The
authors attribute these results to numerical dissipation and invoke
cooling processes – which resumes gravitational contraction and
shearing motions – to amplify the magnetic field.

As we already advanced, our results are much more in
agreement with those of <cit.>. In this work, considering
also an adiabatic simulation with a comparable effective numerical
resolution, the authors manage to amplify the magnetic field to values
of the order of μG (starting from a seed field of 0.1 nG) and
typically one order of magnitude above the value expected if the field
purely follows the gas (B ∝ρ^2/3). According to these
authors, at the highest spatial resolution (≈ 4 kpc) the
numerical dissipation of their code is small enough to resolve the
small-scale turbulent dynamo process. In this process, the magnetic
field lines are stretched following the fluid particles random motions
until diffusion stops the amplification. It is then critical for the
success of the process that i) numerical viscosity remains small
enough to let the code maintain turbulent motions at the required
spatial scales, and ii) numerical magnetic resistivity remains small
enough to let the turbulent flow stretch the magnetic field lines
efficiently before non-ideal effects kill the process. The
  similarity of the magnetic intensity distribution shown in 
  Fig. <ref> with the corresponding one in <cit.>
  (Fig. 2), might point to the same process as the responsible of the
  additional magnetic field amplification seen in
  Fig. <ref>. This should be verified in future works. As a first
  step, we have confirmed that the modulus of the vorticity (not
  shown) and the magnetic field intensity are closely correlated,
  pointing to turbulence as the origin of the added field
  amplification. 

To assess the performance of our code in keeping magnetic
field divergence errors under control, and as in previous sections, we
study the normalized divergence of the comoving magnetic
field, B. In this particular test, we modify the definition
of this quantity to reduce the weight of spurious values associated to
cells where the magnetic field intensity is very low. Following
<cit.>, we change | B_i| in the denominator of the
normalized divergence by | B_i| + 0.01 | B_ max|,
where | B_ max| is the maximum of | B_i| across the
numerical grid.

Figure <ref> displays the median of the relative error (blue 
thick line), and the relative error of the 25% (yellow line), 75%
(orange line), and 90% (green line) percentiles, respectively. The
black line shows the maximum relative error of ∇· B 
according to the previous definition. Finally, let us note that with
90% of the numerical cells with a divergence relative error well
below 3× 10^-2 in ∇· B, the spurious
magnetic energy change related to the divergence error must be of the
order of ∼ 10^-3 or less, i.e., much smaller than the increase
of magnetic energy associated to the added magnetic field
amplification.







The most relevant structure formed in this computational box is
the galaxy cluster located at the box centre. The radial profile of the average value  of the error of the divergence of
B as a function of the cluster-centric distance is displayed as the blue line in Fig. <ref>. The
average values are computed as volume weighted means in spherical
shells. 
This average radial error can be compared with the average radial value of | B| from the simulation together with the expected value of the magnetic field intensity from a pure compression (red lines in Fig. <ref>).
The
relative error exhibits an acceptable behaviour with values below 4%.







To spot possible large errors that could be concealed in the radial
average values of the normalized divergence, we present in 
Figure <ref> a map of the error of the divergence of B in the
same thin slice through the cluster centre displayed in
Fig. <ref>. As seen in the figure, the point values of the
normalized divergence in the central region of the cluster are smaller
than 20% with a mean value 0.5%, in agreement with the median
error shown in Fig. <ref>. 








§ CONCLUSIONS




The role of magnetic field in cosmological scenarios has traditionally
been a matter of debate. Most of the times, they have been pushed into
the background as minor actors in the formation and evolution of
cosmic structures, whose part was completely subdominant. However,
this picture is rapidly changing in recent years. In the epoch of the
so called precision cosmology, where theory and data are widening our
knowledge on the formation of galaxy and galaxy clusters, the
magnetic fields can play a relevant role, both in the theoretical and
observational planes.   

In this line, we presented a new version of an already well-tested and
used cosmological code, that includes a proper description of the MHD
processes. Several commonly used tests, whose solutions are perfectly
known, are shown. In all of them, the performance of our code is more
than satisfactory.  

As a final application with cosmological interest, we followed the
evolution of gas and dark matter components in a cosmic volume, where
a uniform magnetic field was introduced at a very early stage of the
evolution (z∼100). In this simulation, the magnetic field was
amplified from its initial uniform value of B_0∼ 0.1 nG
until values of the order of B∼ 1  μG in the most
massive galaxy clusters. Our simulation shows how the magnetic field
is channeled along the filament of gas and it is amplified at the
structures formed at the intersections of such filaments. The topology
of the magnetic field lines clearly reveals how the field lines tangle
up in the galaxy clusters. The comparison between the values of
B obtained in the simulation and those expected from pure
compression (B ∝ρ^2/3) reveals that there is an
additional process operating at the core of the cluster responsible
for the extra amplification of the magnetic field. To elucidate the
nature of this process, essential to understand the origin of magnetic
fields in clusters of galaxies, will be the subject of future research.

It is nowadays clear that cosmic magnetic fields, although dynamically
negligible, play a crucial role in shaping the physical properties of
the intergalactic medium. From an observational point of view, current
radio observations, such as those with JVLA or LOFAR, provide a unique
tool to reveal the current distribution of magnetic fields
<cit.>. In this sense,
current and next generation of radio facilities, like the SKA and its
precursors, will achieve, mainly by means of  the Faraday tomography
technique, an unprecedented detail in the description of the evolution
and distribution of cosmic magnetic fields. From a theoretical point
of view, while in the last years cosmological simulations have been
significantly improved in terms of the development of different
feedback models to explain galaxy formation, a proper numerical
understanding of cosmic magnetic fields and gas turbulence phenomena
remains as a challenge. Therefore,   in order to  interpret future
radio observations and to shorten the existing distance between the
observational and the theoretical planes, an accurate and `realistic'
description of cosmic magnetic fields in full cosmological simulations
seems to be imperative.  
 



§ ACKNOWLEDGEMENTS

This work has been supported by the Spanish Ministerio de Ciencia,
Innovación y Universidades (MICINN, grant AYA2016-77237-C3-3-P) 
and by the Generalitat Valenciana (grant PROMETEO/2019/071). JMM
also acknowledges partial support from MICINN (grant
PGC2018-095984-B-I00). 
The authors are indebted to D. Ryu and T.W. Jones for providing them
with the analytical solutions for the one-dimensional tests shown in
the paper and ackownledge fruitful discussions with
P. Cerdá-Durán. We thank the anonymous referee for his/her
constructive criticism. Simulations have been carried out using the
supercomputer Lluís Vives at the Servei d'Informàtica of the
Universitat de València. 

mnbst
 





