Magnetic field structure due to the global velocity field in spiral galaxies

We present a set of global, self-consistent N-body/SPH simulations of the dynamic evolution of galactic discs with gas and including magnetic fields. We have implemented a description to follow the evolution of magnetic fields with the ideal induction equation in the SPH part of the Vine code. Results from a direct implementation of the field equations are compared to a representation by Euler potentials, which pose a div(B)-free description, an constraint not fulfilled for the direct implementation. All simulations are compared to an implementation of magnetic fields in the Gadget code which includes also cleaning methods for div(B). Starting with a homogeneous seed field we find that by differential rotation and spiral structure formation of the disc the field is amplified by one order of magnitude within five rotation periods of the disc. The amplification is stronger for higher numerical resolution. Moreover, we find a tight connection of the magnetic field structure to the density pattern of the galaxy in our simulations, with the magnetic field lines being aligned with the developing spiral pattern of the gas. Our simulations clearly show the importance of non-axisymmetry for the evolution of the magnetic field.


INTRODUCTION
Radio observations have revealed that disc galaxies are permeated by large scale magnetic fields ordered on kpc scales and beyond (Beck & Hoernes 1996, Hummel & Beck 1995, Beck et al. 1985. The typical field strength, determined from polarization, Faraday rotation and energy equipartition is of the order of 10 µG (e.g. Beck 2004). The spatial structure of the B-field reflects the spiral and/or barred structure of the gas distribution within the galactic discs (Beck 2008). For example, Fig. 1 shows optical observations of the spiral galaxy M51 overlayed with contours of total synchrotron intensity (tracing the total magnetic field) and magnetic field vectors. It reveals the tight connection of magnetic field with the gas distribution in the galactic disc.
The motion of the gas within the gravitational potential of a galaxy strongly influences the strength and direction of the magnetic field in the interstellar medium. This can be seen by inspecting the well known induction equation of magnetohydrodynamics (MHD), i.e. the temporal evolution ⋆ E-mail: kotarba@usm.lmu.de equation for the magnetic field, where v denotes the gas velocity and η represents the magnetic diffusivity which is inversely proportional to the electrical conductivity. Apparently, within the frame of MHD, the role of the galaxy as a whole is simply to provide for the gas velocity field. Since the conductivity of the interstellar medium is very high, the magnetic field is closely coupled to the gas motion. It is this 'frozen-in'-property of both, the magnetic field and the gas, which determines the spatial structure of the magnetic field. In other words, a detailed investigation of the velocity field of the interstellar gas in disc galaxies is necessary for a deeper physical understanding of the evolution of galactic magnetic fields.
The gas in the disc rotates differentially within the global gravitational potential. Angular momentum transport via spiral arms, bars and gravitational interaction forces the gas to move towards the central regions, and eventually, star formation activity in the disc (superbubbles, winds etc.) drives gas perpendicular to the plane of the disc towards the galactic halo. In general, the axisymmetric rotation velocity Figure 1. Optical image of M51 (Hubble) overlayed with contours of total synchrotron intensity as measure for the total magnetic field (combined observations at Effelsberg and VLA at 6 cm) and vectors of magnetic field. From A. Fletcher & R. Beck (MPIfR) and Hubble Heritage Team (STScI), published by 'Sterne und Weltraum', September 2006. is the dominant component, followed by non-axisymmetric and radial components. The velocity components perpendicular to the disc are typically the smallest. Altogether, v in Eq. 1 represents a complex three-dimensional nonaxisymmetric velocity field strongly coupled to the global properties of the galaxy, including the dark matter halo, stellar component and internal disc activity.
Beside the large scale components of the gas velocity field there are also small scale velocity fluctuations of interstellar gas driven by all kinds of local disc activity, i.e. stellar winds, supernova explosions, cloud-cloud collisions, galactic winds, etc (see e.g. Ferriere 1992, Efstathiou 2000, Johansson & Efstathiou 2006, Kulsrud & Zweibel 2008, Gressel et al. 2008. These unordered velocity components generate two effects which are known as helicity (in terms of a convective turbulent motion perpendicular to the disc) and turbulent diffusion (magnetic field lines with antiparallel direction reconnect and annihilate partially). Helicity supports the amplification of the magnetic field, whereas turbulent diffusion reduces the field strength (see, e.g. Brandenburg & Subramanian 2005 for a review of nonlinear dynamo theory). Therefore, an incorporation of these small scale velocity components into the analysis requires some manipulation of the induction equation (Eq. 1) in terms of a mean-field theory (Steenbeck & Krause 1969, Wielebinski & Krause 1993, Sur et al. 2007. Within the frame of the mean-field description the velocity and magnetic fields are considered as superpositions of the mean and fluctuating parts (v = v + v ′ and B = B + B ′ ). The fluctuating velocity components are coupled to smallscale fluctuations of the magnetic field. The coupling terms are then given by ∇ × α B , where α = 1 3 τ v ′ · (∇ × v ′ ) (Zeldovich et al. 1983), and by ηT ∆ B , where ηT now describes the turbulent diffusion coefficient ηT ∝ v turb · l turb , where v turb and l turb are the typical velocity and length scale of the turbulent motion, respectively.
This leads to the dynamo equation where we have neglected the diffusivity and dropped the mean-brackets for convenience (here, and in the following, B and v refer to their mean values). Eq. 2 is the central equation of cosmic mean field dynamos. It describes the circle of amplification of the different components. The classical dynamo model describes the amplification of the magnetic field through the following chain of α (convective turbulence) and Ω (differential rotation) actions: The radial component Br is amplified through α-action from turbulence; then Bϕ is generated from Br through Ω-action from the shear of the galactic differential rotation. Such an αΩ mean field dynamo amplifies the magnetic field by repeating the chain of α and Ω actions (see Widrow 2002 andStefani et al. 2008 for a review of dynamo theory). However, the origin of the αeffect is still under discussion (Cattaneo & Vainshtein 1991, Vainshtein & Cattaneo 1992, Kulsrud & Anderson 1992.
We emphasize that the described classical dynamo models use only one velocity component, the differential rotation.
To be more precise the role of any deviation from axisymmetry is considered to be unimportant for the evolution of the large-scale magnetic field, which is not necessarily true in real galaxies.
On this account, there have been three-dimensional numerical simulations using an analytical turbulent velocity field, where deviations from axisymmetry were incorporated in the gas-and turbulence-profiles (Rohde et al. 1997, Rohde & Elstner 1998. These studies showed, that even accounting for the α-effect calculated out of the analytical velocity field an initial magnetic field cannot survive for more than 500 Myr. Moreover, Elstner et al. 2000 performed N -body simulations of two component (collisionless stars and gaseous clouds moving in the gravitational potential of the stellar population), self-gravitating discs embedded in an analytical bulge-and halo-potential. These simulated clouds provided an already very good approximation of the gas velocity field. However, full hydrodynamics was not incorporated. The obtained velocity field was used in an αΩ-dynamo description. Without including the α-effect, the non-azimuthal 3D gas flow alone did not provide an amplification of the magnetic field. The field got amplified by several orders of magnitude within 0.7 Myr only when the α-effect was included. In addition, they found an alignment of the magnetic field with the developed spiral pattern of the disc.
Recently, Dobbs & Price (2008) performed three dimensional, full MHD, single and two-component (cold and hot gas) simulations using smooth particle hydrodynamic (SPH) methods to treat MHD. They applied a spiral potential to the gas, thus, the self-induced formation of spiral structure was not included. Their work concentrated on structure formation in the disc, like molecular clouds and inter-arm spurs. They found that the main effect of adding a magnetic field to these calculations was to inhibit the formation of structure in the disc. They did not consider global enhancement and structure formation of the magnetic field, but nevertheless, they found that the global magnetic field was following the large scale velocity field.
It is the aim of this paper to present further steps towards a more complete dynamo model. We perform for the first time a set of self-consistent N -body calculations of a spiral galaxy including hydrodynamics as well as the induction equation via the SPH method to obtain the complex three dimensional velocity field. Compared to all previous work, we use no analytical potential for any component of the galaxy. All components (disc, gas, bulge and halo) are represented by particles which are treated as self-gravitating N -body-particles, while hydrodynamics is applied to the gas component only. We use more than one order of magnitude more particles than Elstner et al. (2000). We follow the evolution of the magnetic field according to the induction equation (eq. 1). Thus, we have implemented the SPH variant of the induction equation as well as the representation of the magnetic fields by Euler potentials in the SPH code Vine and compare the results with simulations performed using the SPH code Gadget. N -body/SPH methods are well adapted for simulating whole galactic discs as the simulated discs stay stable for at least 15 dynamical times (where we define the dynamical time for a disc galaxy as its half mass rotation period). As we will show in section 4, our discs are forming spiral structure without applying a spiral potential or any other mechanism to provide extraordinary flows.
In summary, we investigate here the kinematic reaction of a large-scale magnetic field on the complete threedimensional, large-scale velocity field of a disc galaxy obtained from the N -body SPH simulations, using two different numerical codes.
The paper is organized as follows: Section 2 gives shortly the theory of magnetic field evolution in differentially rotating systems. A summary of the SPH method and the treatment of magnetic fields including the method based on Euler potentials is given in section 3. The simulations together with a comparison of the performance of the Vine and Gadget codes are presented in section 4. The results are discussed in section 5, where we also analyse the terms of the induction equation in detail. Finally, we summarise and conclude in section 6.

THEORETICAL EXPECTATIONS
When only studying the effect of the gas velocity on the evolution of the magnetic field, we can neglect the diffusive term in eq. 1. Keeping this term, one would physically except the magnetic field to dissipate depending on the value of η and reconnect when converse magnetic field lines come together. Technically, η is not always assumed to be spatially dependent, so that the diffusive term reads −η∇ 2 B. However, this formulation leads only to an effective smoothing, and not a real diffusion of the magnetic field. Neglecting the diffusive term thus corresponds to considering an upper limit of field amplification. Additionally, η is assumed to be small except within strong shocks.
The induction equation 1 then yields and applying cylindrical coordinates eq. 3 reads: These equations can be simplified to get a first idea of how magnetic fields will evolve in a galactic disc. For a differentially rotating disc (∂vϕ/∂r ≈ 0) with a perfectly axisymmetric velocity field, v does not depend on ϕ. The same holds for an axisymmetric magnetic field. If we also assume that changes of all quantities in the z direction are small compared with those in the radial direction and Bz ≃ 0, the equations for the regular field, i.e. the field in the plane of the disc, read: where Ω is the angular velocity.
In the absence of radial flows, the last term of eq. 8 describes the generation of a toroidal magnetic field from the radial component of the already present magnetic field by differential rotation. This effect is the so called Ω-effect already mentioned above. Since vϕ ≫ (vr, vz) this term is dominant and one would expect any initial magnetic field to be first wound up by differential rotation. However, this effect alone cannot be responsible for a significant amplification of the magnetic field, as the amplification stops when all of the radial field is wound up. However, if a gas flow in the negative radial direction (i.e. towards the centre of the disc) is present the radial field can be amplified and then be converted into a toroidal field. These radial gas flows occur when angular momentum is transported in the gas out of the disc, e.g. by spiral arms or bars. The toroidal magnetic field can be amplified further if the radial gas flow velocity decreases with increasing galactocentric radius.
Therefore, a good understanding of the evolution of magnetic fields in galactic discs requires full information of the three-dimensional velocity field of the gas which is naturally provided by self-consistent numerical simulations. We will discuss the velocity field in our simulations and the re-sulting values of the different terms of the induction equation in section 5.

Vine
The equations presented below are implemented within the OpenMP parallel N -body/SPH evolution code Vine. For all details we refer the reader to Wetzstein et al. (2008) and Nelson et al. (2008).

SPH basics
Within the SPH formulation a hydrodynamic quantity A is interpolated by a kernel function W (r− r ′ , h) with R W dr = 1 and lim h→0 W = δ(r − r ′ ), where the so called smoothing length h defines the spatial extent of the function W . This interpolation is then discretised, so that where i (j) is the index of the particle at position ri (rj ) and Ai (Aj) the value of the quantity A at the position of particle i (j). ρj and mj denote the density and mass at position of particle j, respectively. The Vine code uses the common W4 kernel defined by Monaghan & Lattanzio (1985) as where values with index ij denote differences (e.g. rij = ri − rj ) and arithmetic means (e.g.hij = 0.5 · (hi + hj)), respectively, ̺ = |r − r ′ |/h, ν is the number of spatial dimensions of the system and σ is a constant of order unity. See Monaghan (1992), Monaghan (2001) or Price (2005) for more details.

Continuity equation
As long as the kernel itself is differentiable, every function A can be interpolated to a differentiable function by the procedure described above. The most common formulation of derivatives in SPH is (see e.g. Monaghan 1992, Price 2005) Using the continuity equation, the total time derivative of the density thus reads

Momentum and Energy equation
A natural ansatz to derive a conservative form of the momentum equation comprising the force due to pressure gradients (in addition to the force due to the gravitational potential) is to use the Lagrange formalism together with the first law of thermodynamics. This leads to the following SPH variant of its ideal form ( dv dt = − ∇P ρ ): In this formulation momentum is conserved exactly, since the contribution of particle j to the momentum of particle i is equal and negative to the contribution of particle i to the momentum of particle j.
The change in the thermodynamical state of the gas requires an evolution equation for a state variable corresponding to the internal energy or entropy of the gas. Vine employs an equation for the specific internal energy (u) of the gas. Without external heating or cooling terms, only compressional heating and cooling are important and the SPH variant of the ideal form ( du dt = − P ρ ∇ · v) reads: To close the set of equations, an isothermal equation of state is used throughout this paper.

Artificial viscosity
Artificial viscosity is required to model shocks and angular momentum transport properly, where the latter is important to be able to simulate spiral arm formation. The Vine code uses the most common form of the artificial viscosity. It is described by the tensor Πij as in Monaghan (1992).
Since the value of Πij depends on the difference in velocity between the considered particles (i.e. the velocity gradient) the viscosity increases with increasing velocity gradient. Moreover, the viscosity is only applied if particles are approaching each other. Balsara (1995) suggested a viscosity limiter to avoid spurious angular momentum and vorticity transport in gas disks. However, a lower viscosity leads to a higher velocity dispersion of the gas and therefore to higher divergence of the velocity and magnetic field. As will be shown and discussed in section 5, this higher divergence causes a more violent magnetic field amplification.
The viscous terms within the momentum and energy equations read: This treatment of viscous forces allows for a sensible description of the behaviour of gas in a spiral galaxy. However, eqs. 14 and 16 are not applied when using isothermal equation of state.

Induction equation
In order to follow the evolution of the magnetic field we have additionally implemented equation 3 discretised as with µ, ν denoting the spatial directions.

Euler potentials
A well known problem related to magnetic fields within SPH is the maintenance of ∇ · B = 0 throughout the simulation. Different attempts to solve this problem have been made (see Price & Monaghan 2005), examples include source term approaches (Powell et al. 1999) and projection methods. Theoretically, the problem can be avoided if the magnetic field is represented by Euler potentials (Stern 1970, Price & Bate 2007, Rosswog & Price 2007, an approach we have also implemented into the Vine code.
The magnetic field is expressed as a function of two scalar potentials αE and βE as: Taking the divergence of B we get: Thus, the divergence constraint is fulfilled by construction. Moreover, for the ideal case (η = 0) the Euler potentials for each particle (i.e. the convective values of the potentials) are direct tracers of the magnetic field and stay constant with time, thus one does not need to perform an additional integration when following magnetic fields, leading to a higher accuracy of the calculation. The variation of magnetic field is only due to the motion of the particles, which corresponds to the advection of magnetic field lines by Lagrangian particles (frozen flow) (Stern 1970). Moreover, the formulation by Euler potentials guarantees conservation of magnetic helicity, which is again reasonable for ideal treatment. Actually, the magnetic helicity is zero, since for A = αE∇βE and equivalently A = −βE∇αE, respectively, B = ∇αE × ∇βE = ∇ × A and therefore A · B = 0. The gradients of αE and βE can be expressed as where which is exact for linear functions (Price & Bate 2007), i.e. for initial conditions with an uniform field. However, the difference in using this exact-linear interpolation compared to the usual gradient operator is marginal (D. Price, private communication). Unfortunately, not all magnetic field configurations can be expressed in terms of Euler potentials easily as they enter eq. 18 in a nonlinear way, and they are not unique for certain field configurations (Stern 1970, Yahalom & Lynden-Bell 2006. The former problem restricts only the choice of the initial magnetic field, whereas the latter can be crucial. If there are field configurations which can be expressed by different sets of Euler potentials, then this implies, that some other field configurations cannot be expressed at all using Euler potentials. However, since the fields considered in this work are topologically 'simple', we do not expect to encounter these problems. Furthermore, Euler potentials do not allow to follow the winding of magnetic fields beyond a certain point. This constraint is due to the fact that using the Euler potentials, the magnetic field is essentially mapped on the initial particle arrangement. If the initial arrangement evolves too much during the simulation, particles carrying conflicting values of Euler potentials (i.e. values, which do no longer allow for a finite and unambiguous calculation of their gradients) can come close. Then, the ability of the Euler potentials to represent the magnetic field correctly is lost. This conflict is expected to occur when the magnetic field is wound up more than once, which poses a problem especially towards the central region of a simulated galaxy.

Timestepping
In Vine, there are basically three different time step criteria, based on changes in the acceleration of a particle, its velocity, or both in combination, where ǫ, a and v are the gravitational softening length, the acceleration and velocity of a particle in the previous time step (n), respectively, and τacc is an accuracy parameter. Two additional time step criteria are applied in SPH simulations: First, the Courant-Friedrichs-Lewy criterion as suggested by Monaghan (1989), where αi and βi are artificial viscosity parameters, cs is the sound speed, hi the SPH softening length for gas particle i, and µij corresponds to the velocity divergence between particles i and j with the maximum taken over all neighboring particles j of particle i (see Wetzstein et al. 2008 for more details). Secondly, there is a limit on how much the SPH softening lengths are allowed to change during one timestep: where τ h is again an accuracy parameter. Usually, we apply τacc = 1, τCFL = 0.5 and τ h = 0.15. The timestep actually employed in the simulation is the minimum of the timesteps in eqs. 26-30.

Gadget
A somewhat different treatment of hydrodynamics and magnetic fields is realised within the MPI parallel Nbody/SPH code Gadget (Springel et al. 2001, Springel 2005, Dolag & Stasyszyn 2008. There are two significant differences in the implementation relevant even for nonradiative simulations: First, Vine follows a classical implementation which is integrating the internal energy, whereas Gadget utilises what is generally called the entropy conserving formulation. The important difference thereby is not the fact that Gadget integrates the entropy instead of the internal energy. The crucial differences are rather the way in which the smoothing length hi is defined (in Gadget, hi is defined based on the mass within the kernel instead of the number of particles) and the inclusion of correction terms arising from the varying smoothing length. Also, the entropy conserving formulation uses a way of symmetrizing the kernel given by the derivation of the SPH equations, which in sum leads to conservation of energy and entropy at the same time (Springel & Hernquist 2002).
The second difference originates in an alternative formulation of the artificial viscosity. In Gadget, artificial viscosity is based on the signal velocity instead of sound speed (Monaghan 1997) and apt to incorporate magnetic waves in a natural way (Price & Monaghan 2004).
This different implementation was shown to bring measurable improvements specially for MHD applications (Dolag & Stasyszyn 2008), but should not make too much of a difference for passive magnetic fields. The implementation of the induction equation and the Euler potentials formalism is the same in both codes.
The integration in Gadget is also performed using the leapfrog integration scheme, but Gadget utilises a kick-drift-kick-scheme whereas Vine uses a drift-kick-driftscheme.
The timestep is given by where η translates to the accuracy parameter τacc in eq. 26 via τacc = √ 2η. For SPH particles, also a Courant-like condition in the form is applied, where hi is the SPH softening length for gas particle i and v sig ij the signal velocity between particles i and j as defined in Price & Monaghan (2004) with the maximum taken over all neighboring particles j of particle i. Ccour total mass Mtot = 1.34 · 10 12 M ⊙ disc mass M disc = 0.041 · Mtot bulge mass M bulge = 0.01367 · Mtot mass of the extended gas disc Mgas = 0.2 · M disc exponential disc scale length l D = 3.5 kpc scale height of the disc is an accuracy parameter which does not translate one-toone to τCFL in eq. 29 due to the different definition of the Courant criterion. We commonly use values of η = 0.02 and Ccour = 0.15 to ensure that the timestep ∆t in Gadget does not get too large compared to Vine. However, changing the accuracy parameters by a factor of two does not affect the overall evolution and amplification of the magnetic field in the simulated systems (not shown). Beside that, the codes differ in details of the tree construction for calculating gravitational forces. For more details we refer the reader to the code papers for Vine  and Gadget (Springel 2005, Dolag & Stasyszyn 2008.

Setup
The initial conditions for our Milky Way like galaxy are realised using the method described by Springel et al. (2005) which is based on Hernquist (1993) (see also Johansson et al. 2009). The galaxy consists of an exponential stellar disc and a flat extended gas disc, a stellar bulge and a dark matter halo of collisionless particles. The gas is represented by SPH particles adopting an isothermal equation of state with a fixed sound speed of cs ≈ 15 km s −1 , which corresponds to a temperature of T ≈ 2 · 10 4 K for a molecular weight of 1.4/1.1·mproton . We briefly note that by using an isothermal equation of state only one component of the ISM is modeled, typically this is a reasonably good approximation for the warm gas phase in disc galaxies (e.g. Barnes 2002, Li et al. 2005, Naab et al. 2006. Assuming an isothermal equation of state implies that additional heat created in shocks by adiabatic compression and feedback processes (e.g. by SNII) is radiated away immediately. In addition, substantial heating processes prevent the gas from cooling below its effective temperature predefined by its sound speed.
The parameters describing the initial conditions can be found in Table 1. The particle numbers and the gravitational and SPH softening lengths used in the different runs can be found in Table 2.

Figure 2.
Surface densities Σgas of the extended gas discs as a function of radius before the inclusion of the magnetic fields after 0.5 Gyr (red line) and after 2 Gyr (black lines) for simulations with Gadget (solid line) and Vine (dotted line). The gas discs are stable for more than ten half mass rotation periods.
To set up the corresponding Euler potentials, we choose We have checked the stability of our discs in independent simulations without magnetic fields. Figs. 2 and 3 show the surface densities Σgas of the extended gaseous discs and Σstars of the exponential stellar discs, respectively, as a function of radius for t = 0.5 Gyr (red), i.e. the time at which the magnetic field is switched on, and t = 2.0 Gyr (black). Fig.  4 shows the circular velocity curves of the simulated galaxies at the same times. The discs simulated with Vine (dotted line) and Gadget (solid line) show similar results and stay stable over more than ten half mass rotation periods.

Direct magnetic field simulations
Figs. 5 and 6 show the face on view of the magnetic field energy and gas density of the simulated galaxy at different  Fig. 2 but for the stellar surface densities Σstars. Both the stellar and the gas discs are stable for more than ten half mass rotation periods. output times. The magnetic field was switched on at t = 510 Myr. The viscosity limiter was not applied. Fig. 5 shows the simulation performed with Vine and Fig. 6 the same initial conditions simulated with Gadget. The magnetic field energy B 2 /8π is colour coded and normalised to the initial value of 1 8π · 10 −18 erg cm −3 on a logarithmic scale from 1 (blue) to 1.5 · 10 8 (red). The contours overplotted indicate physical densities of 23, 37 and 52 M⊙ pc −3 , respectively. We use a grid with a cell size of 0.3 kpc for the calculation of the mean values of the densities and the magnetic field energies, averaging in the vertical direction from -h to h, where h is the local height of the gas disc.
In both simulations we see that the magnetic field energy pattern is tightly connected to the density pattern of the gas. Moreover, both simulated galaxies show a very similar morphology in the gas and magnetic field distributions. The magnetic field energy in the spiral arms is amplified by up to five orders of magnitude in both codes and even more in the central region (see also Fig. 13). Furthermore, the SPH smoothing lengths hgas are similar for both codes (Fig. 9), indicating that the performance of the hydrodynamic calculations is concerted. The smoothing lengths in Figure 5. Face-on magnetic field energy and gas density as a function of time for the simulation performed with Vine using direct magnetic field description and without applying the viscosity limiter. The colours correspond to the magnetic field energy B 2 /8π on a logarithmic scale, normalised to the initial value of 1 8π · 10 −18 erg cm −3 . The contour lines indicate physical densities of 23, 37 and 52 M ⊙ pc −3 , respectively.   Fig. 5, this time the magnetic field is followed using Euler potentials implemented in Vine. In contrast to the direct simulation the magnetic field is more strongly amplified in the spiral arms than at the centre. The maximum amplification of the magnetic field energy is only three orders of magnitude. Note that the colour scaling is different to Figs. 5 and 6.  . Circular velocity curves of the simulated galaxies at two different times. The colour coding is the same as in Figs. 2 and 3. Again, the circular velocity curves are stable over more than ten half mass rotation periods. Figure 9. SPH softening lengths hgas as a function of radius shortly after the inclusion of the magnetic fields at 0.55 Gyr (red line), at 1.25 Gyr (green lines) and at 2 Gyr (black lines) for simulations with Gadget (solid lines) and Vine (dotted lines). The SPH smoothing lengths are very similar for both codes.
Vine are initially set to a constant value of hgas ≈ 0.3 kpc at the time of the magnetic field inclusion.

SPH with Euler Potentials
Figs. 7 and 8 show simulations starting from the same initial conditions as before. However, this time the evolution of the magnetic field was followed using the Euler potentials. Again, we show magnetic field energies and gas densities. This time the amplification of the magnetic field energy in the spiral arms is only three orders of magnitude for both simulations with Vine and Gadget, with both showing a remarkably similar evolution. The most notable difference to the simulations with direct magnetic field treatment shown in Fig. 5 and 6 is at the centre of the galaxies, where in the direct simulations the field amplification was strongest. With Euler potentials the magnetic field grows mostly in the spiral arms of the galaxy (see also Fig. 13).
Since the magnetic fields in our simulations are passive, the density profiles (Figs. 2 and 3) of the disc are the same for all runs. Thus, the different profiles of the magnetic Figure 10. Numerical h · |∇ · B|/|B| at t ≈ 1.5 Gyr as a function of radius for simulations without applying the viscosity limiter using direct magnetic field treatment (blue) and using Euler potentials (black) in Gadget (solid lines) and Vine (dotted lines). Direct magnetic field simulations for which the viscosity limiter was applied are also shown (orange). Using direct magnetic field description, the numerical h · ∇ · B is highest at small radii, and much larger than in the Euler potentials formalism. field energy cannot be traced back to the density profiles. In fact, it is the numerical ∇ · B which presumably causes the high amplification of the magnetic field at the centre in simulations with the direct magnetic field treatment. Fig.  10 shows the radial profile of the numerical h · |∇ · B|/|B| at time t ≈ 1.5 Gyr for simulations using direct magnetic field treatment (blue for simulations without applying the viscosity limiter and orange where the limiter was applied) and Euler potentials (black) performed using Gadget (solid lines) and Vine (dotted line). Utilising the direct magnetic field description, the numerical ∇ · B is highest at small radii, and much larger than for the Euler potential formalism. As will be discussed in the following section, high ∇ · B corresponds to high amplification of the magnetic field. Fig. 11 shows the magnetic field vectors for the normal resolution Vine simulation utilising Euler potentials at the time t ≈ 0.9 Gyr. This time the colours correspond to the gas density on a logarithmic scale from 0.3 · 10 −3 to 2.3 · 10 3 M⊙ pc −3 , overplotted with the field vectors. The length l of the vectors is normalised to the initial value and displayed logarithmically as l = 3 · log(B/B0), i.e. l = 0 corresponds to B ≈ B0 or smaller, l = 1 to B ≈ 2 · B0, l = 2 to B ≈ 5 · B0 and l = 3 to B = 10 · B0. The magnetic field lines follow the spiral structure of the gas. They have been amplified by contraction in regions of higher density and restructured by differential rotation of the galaxy. Their orientation is caused by the motion of the gas. These characteristics are very similar to typical observations of magnetic fields in galactic discs (e.g. Fig. 1).
Qualitatively, this behaviour is the same for all simulations using both codes. Only the central region in simulations using direct magnetic field treatment shows chaotic orientation of the magnetic field lines, indicating artificial amplification of the magnetic field due to high numerical ∇ · B.  Fig. 6, this time the magnetic field is smoothed every 30 timesteps. The field morphology is similar to the morphology in Fig. 6 and 5, respectively, but the magnetic field values in the spiral arms are now more similar to the values in the Euler implementation. Figure 11. Gas density (colour coded) and magnetic field vectors for the normal resolution simulation with Vine using the Euler potentials formalism at t ≈ 1 Gyr, i.e. ≈ 500 Myr after the inclusion of the magnetic field. The length of the vectors is normalised to the initial value and displayed logarithmically. l = 0 corresponds to B ≈ B 0 or smaller and l = 3 to B = 10 · B 0 .

Magnetic field growth
Figs. 5 (6) and 7 (8), respectively, reveal the differences in the magnetic field amplification for the direct magnetic field treatment and the Euler potentials formalism: Using the di-rect description, the amplification of the magnetic field energy in the spiral arms is higher by at least two orders of magnitude, and at the centre even more than six orders of magnitude compared to the Euler potentials method. This difference is probably caused by the numerical ∇ · B in these simulations (Fig. 10), but possibly also by the fact that field winding is not traced beyond a certain evolutionary state in the Euler potentials formulation (see section 3.1.6). Since the Euler potentials are free from physical divergence by construction (i.e. the divergence is zero to measurements errors), the numerical divergence in simulations using the Euler potentials is due to the SPH derivative approximation when calculating the magnetic field from the potentials (Eq. 18). In this sense, the numerical divergence found in simulations using Euler potentials reflects the ability of SPH operators to measure the gradient of a curl to zero. Thus, the fact that ∇ · B is higher by approximately one order of magnitude in the disc (i.e. within ≈ 5 to 15 kpc) and by several orders of magnitude at the centre (Fig. 10), presumably causes the different magnetic field amplification in these simulations. This is the case at least in the disc region, where the winding of the field is not strong enough to constrain the Euler potentials formulation.
To get a better idea of the influence of numerical ∇·B on the amplification of the magnetic filed, we have performed simulations applying magnetic field smoothing, a technique allowing for reduction of small scale fluctuations and therefore also the numerical divergence (Dolag & Stasyszyn 2008). Within this method, the magnetic field is smoothed periodically as suggested by Børve et al. (2001). Fig. 12 shows again the magnetic field energies and gas densities for a Gadget simulation starting from the same initial con- Figure 13. Total magnetic field (Btot = q B 2 x + B 2 y + B 2 z ) at t ≈ 1.5 Gyr as a function of radius for the normal resolution Gadget (solid lines) and Vine (dotted lines) simulations without applying the viscosity limiter using Euler potentials (black) and the direct magnetic field description (blue). Direct magnetic field simulations for which the viscosity limiter was applied are also shown (orange). Two direct magnetic field implementations including field smoothing to reduce the numerical ∇ · B contribution are shown in red and green. The simulations including smoothing have been run with a smoothing interval of 30 (red) and 5 timesteps (green), respectively. Increasing the frequency of smoothing tends to decrease the amplitude of the magnetic field between 3 and 10 kpc, but has relatively little effect for larger radii.
ditions as before and without applying the viscosity limiter. This time, the magnetic field was smoothed every 30 timesteps. Applying the smoothing scheme, the amplification of the magnetic field energy is reduced to approximately three orders of magnitude within the spiral arms, which is the same as the amplification seen in simulations using the Euler potentials, and it is also lowered towards the centre of the galaxy. The structure of the magnetic field is despite the smoothing still very similar to the other runs and again correlates well with the structure of the gas density, however, the magnetic field energy is more concentrated within the spiral arms. Fig. 13 shows the total magnetic field at t ≈ 1.5 Gyr as a function of radius for the normal resolution Gadget (solid line) and Vine (dotted line) simulations using Euler potentials (black) and the direct magnetic field description without applying the viscosity limiter (blue) and with the limiter turned on (orange), respectively. The direct magnetic field simulations including field smoothing are shown in red and green. They have been performed with a smoothing interval of 30 (red) and 5 timesteps (green), respectively, and without applying the viscosity limiter. As discussed before, the most notable difference between simulations with direct magnetic field treatment and the Euler implementation is at the centre of the galaxies. There, the amplification in the direct simulations is much stronger than in the Euler simulations. This behaviour could be at least partly physical, as there are high radial velocities and strong in-and outflows of gas in the central region ( fig. 17), resulting according to Eqs. 7 and 8 in an amplification of the magnetic field. In addition, also the azimuthal derivatives of the radial and toroidal velocity components are large at the very centre, which also could account for the violent amplification (see section 5.3). On the other hand, the second term of eq. 8 does not play an important role, since dvϕ/dr is large and therefore dΩ/dr small in the central region (by reason of solid body rotation). However, the high ∇ · B values at the centre make it difficult to distinguish between physical growth and numerical errors. Since the Euler potentials are also unreliable in this region (see section 3.1.6), it is not easy to decide which formalism is the most capable in describing the physics in the centre of the galaxy correctly. This is also true for the simulations including smoothing. Increasing the frequency of smoothing tends to decrease the amplitude of the magnetic field between 3 and 10 kpc, but has relatively little effect for larger radii. Interestingly, the large increase of B in the centre is never smoothed away, which could indicate, that this behaviour is actually partly physical. For the simulation which applies smoothing every 30 timesteps (red), the amplification of the field at r > 3 kpc is similar to the simulations with Euler potentials. Applying smoothing every 5 timesteps (green), the amplification is considerably weaker than in the Euler potentials simulations, indicating that by such strong smoothing essential physics is lost, in agreement with earlier findings by Dolag & Stasyszyn (2008).
In the following, we only consider the disc region (from 5 to 15 kpc), since the high numerical divergence in the centre makes it difficult to lower it to the value of the divergence seen in simulations with Euler potentials (i.e. h · |∇ · B|/|B| ≈ 1), without smoothing the magnetic field structure too much. Fig. 14 shows the evolution of the total magnetic field (Btot = p B 2 x + B 2 y + B 2 z ) within the disc with time for the different implementations. The colour coding is the same as in Fig. 13. As before, for the simulation which applies smoothing every 30 timesteps (red), the amplification of the field is similar to the simulation with Euler potentials. However, the performance of these simulation is not very convincing due to the "jumps" in the evolution caused by the artificial periodic smoothing. Applying smoothing every 5 timesteps (green), the amplification is as discussed before lower than in the Euler potentials simulations.
This behaviour can be understood by considering the corresponding numerical divergence of the magnetic field. Fig. 15 shows h · |∇ · B|/|B| as a function of time for all simulations. In all cases, the growth of h · |∇ · B|/|B| behaves similar to the amplification of the total magnetic field, i.e. the higher the divergence, the stronger the amplification of the field. Though the numerical divergence in the simulation using Euler potentials (black) is higher than in the simulation with a smoothing interval of 5 timesteps (green), its value does not directly correlate with the field growth. That is because the (defective) magnetic field itself is not used for calculating the magnetic field evolution within the Euler potential formalism as is the case for the direct magnetic field description (compare eqs. 17 and 18). Using the smoothing scheme lowers the divergence (in case of smoothing every 5 timesteps even below the numerical divergence of the Euler potential formalism) and lowers also the field amplification, leading (if applied not too often) to an amplification of the total field much more similar to that using the Euler potentials, which are free from physical divergence by construction. Figure 14. Total magnetic field (Btot = q B 2 x + B 2 y + B 2 z ) within the disc (between 5 and 15 kpc) as a function of time for different implementations. The colour coding is the same as in Fig. 13. Applying the smoothing scheme reduces the amplification of the magnetic field.
Interestingly, for simulations applying the viscosity limiter suggested by Balsara (1995), the magnetic field amplification using the direct magnetic field description is in both codes much higher than without applying this limiter (orange lines in Figs. 13 and 14). The reason for this higher amplification is the higher velocity dispersion in these simulations. The viscosity limiter lowers the viscosity in regions of strong shear flows, thus suppressing velocity diffusion and leading to higher velocity gradients. Consistently, also the numerical divergence of the magnetic field is higher (and considerably higher than the "unavoidable" value of approximately one) in these simulations (orange lines in Figs. 10 and 15). Applying the viscosity limiter in simulations using Euler potentials, however, does not change the evolution of the magnetic field significantly (not shown). Therefore, again, it is probable that the higher numerical ∇ · B terms lead via the induction equation (eq. 3) to an enhanced magnetic field growth.
In summary, the field amplification in case of direct magnetic field description correlates with the non-vanishing, numerical ∇ · B. The Euler potential formalism also has its shortcomings (like the non-uniqueness and the dependence of the magnetic field on two derivatives (Eq. 18) leading to lower numerical accuracy). Thus there is a strong need for simulations with different ∇ · B cleaning techniques and even higher resolution in order to be able to distinguish the best description for simulations of magnetic fields in galactic discs.
However, since the physical divergence is zero in the case of the Euler potentials, we believe this method (for the time being) to be the best one for our studies of magnetic convection in disc galaxies. The following discussion therefore concentrates on simulations using Euler potentials. Fig. 16 shows the total magnetic field as a function of time for different resolutions (see Table 2) in simulations with Figure 15. Total divergence of the magnetic field within the disc (between 5 and 15 kpc) as a function of time for different implementations. The colour coding is the same as in Fig. 13. The higher the divergence, the stronger the amplification of the magnetic field. B 2 x + B 2 y + B 2 z ) as a function of time for the Gadget (solid line) and Vine (dotted line) simulations without applying the viscosity limiter using Euler potentials. The total numbers of particles are 1.3 · 10 5 (blue), 1.3 · 10 6 (black) and 1.3 · 10 7 (red).

Numerical resolution
Gadget (solid lines) and Vine (dashed lines) without applying the viscosity limiter.
One Gyr after its initialization the magnetic field has been amplified from 10 −9 to approximately 9 · 10 −9 G in the low resolution simulation (blue), whereas the final magnetic field strength in the normal resolution simulation is slightly more than 1.5 times higher (1.5 · 10 −8 G). The final magnetic field strength in the high resolution run is again approximately 1.5 times higher than in the normal resolution run (i.e. ≈ 2.5 · 10 −8 ). The numerical h · |∇ · B|/|B| values are of the same order for all resolutions (not shown). Thus, we have not yet reached numerical convergence in the magnetic field evolution.

Inspection of the induction equation
By analyzing the velocity and magnetic field in our simulation we can identify the single terms of the induction equation responsible for the behaviour of the magnetic field. Dropping all dependencies on z, the equations for the evolution of the radial and toroidal magnetic fields read where we have labelled the single terms with numbers for easier reference. The radial and toroidal components of the velocity and the magnetic field and their corresponding derivatives are shown in Fig. 17 after approximately three half mass rotation periods after the onset of the magnetic field. The radial velocity (top left) is typically negative, leading to an effective gas inflow towards the centre of the galaxy. This negative radial velocity mirrors the angular momentum transport to large radii of the galaxy by spiral arm formation. The mean circular velocity (second row, left) is 210 km s −1 at large radii, and drops to zero towards the centre (see also Fig. 4). The toroidal magnetic field (bottom left) is wound up by differential rotation, leading to a structure of altering positive and negative magnetic field values from centre to the edge of the galaxy. Consequently the derivatives with respect to ϕ (right panel) are smaller than the radial derivatives (middle panel), mirroring the approximate axial symmetry. However, since the terms of the induction equation depend always on a product between a derivative and a velocity or magnetic field component, one cannot a priori neglect the terms depending on azimuthal derivatives.
In order to quantify the influence of the different terms 1-10 during the simulation we calculated their values in cylindrical bins within the disc (5 to 15 kpc) and their mean value at different times. We have taken the negative values of each term in case of negative magnetic field to distinguish between amplifying and attenuating terms. The result of this calculation is shown in Fig. 18. The upper plot shows the temporal evolution of the terms responsible for amplification/attenuation of the radial magnetic field (terms 1 to 5) and the lower of the toroidal magnetic field (terms 6 to 10). Positive values imply amplification, and negative attenuation of the corresponding B-component. The nonaxisymmetric terms are shown in red.
Looking at Fig. 18, the most important term for the evolution of the radial magnetic field is term 5, i.e. − vϕ r ∂Br ∂ϕ . Since the toroidal velocity dominates the velocity field, this term is most important although ∂Br ∂ϕ is comparatively small. This can be seen following the evolution of the circular velocity and the radial magnetic field more closely: The radial magnetic field is strongest where the circular velocity has its highest value, with a delay of roughly 40 Myr. All other terms lie in the same range and therefore compete with each other. Since their values are positive as well as negative, one should not expect a significant amplification on their account. This analysis shows, that even small deviations from axial symmetry are very important for the evolution of the magnetic field in spiral galaxies.
One reaches the same conclusion looking at the terms of the evolution equation for Bϕ. Except for the beginning of the simulation, the leading term here is clearly term 9, − vϕ r ∂Bϕ ∂ϕ , i.e. the only non-axisymmetric term in this equation. Term 10, which was our candidate for the most important term for axial symmetry, is only the second most important. Both terms depend on the toroidal velocity component, thus demonstrating the importance of the differential rotation for the evolution of the toroidal component of the magnetic field.
Neglecting all non-axisymmetric terms (plotted in red) one finds term 1 (−Br vr r ) to be largely dominant over term term 4 (−vr ∂Br ∂r ), in agreement with the theory for the evolution of Br. Also the term responsible for the evolution of Bϕ is as expected: Term 10 (− vϕ r Br) is the leading term and followed by term 6 (−Bϕ ∂vr ∂r ). However, term 1 and 6 are both of order 10 −13 G Myr −1 , thus not being able to account for any significant amplification of our initial magnetic field, and term 10 can only amplify Bϕ effectively if Br is amplified.
This behaviour is qualitatively the same also for runs with the direct implementation of the induction equation. We conclude that the non-axisymmetry of the system is the driving force for the observed field amplification in our simulations.

CONCLUSION AND OUTLOOK
We have presented a set of self-consistent simulations of the evolution of magnetic fields in galactic discs performed with the N -body codes Vine and Gadget. Hydrodynamics was treated using the SPH method. The evolution of magnetic fields within the framework of ideal MHD was followed by both a direct implementation of the induction equation and a formalism using Euler potentials.
The presented set of simulations shows the importance of a sensible treatment of ∇ · B when simulating magnetic fields in spiral galaxies. Since artificial magnetic monopoles can be responsible for unphysical amplification of the field, more studies of possibilities to avoid or inhibit numerical ∇ · B terms are still needed. Although the description using Euler potentials avoids (physical) magnetic monopoles by construction, the drawback in using them is that they lead to constraints on magnetic helicity. Since helicity fluxes can affect the dynamo process within a mean field theory (Brandenburg & Subramanian 2005), Euler potentials would probably not be suitable for simulations including the α-effect. Furthermore, Euler potentials do not allow for all initial field configurations, since they are not necessarily single valued and in addition, their derivation can become quite complex. Nevertheless, using topologically simple initial conditions for the magnetic field, the Euler potential formalism  seems to be the best tool to follow the ideal evolution of magnetic fields in simulations of spiral galaxies with SPH.
A possible alternative to Euler potentials is the vector potential A. The disadvantages are the need for a time integration of A when evolving magnetic fields and the occurrence of second derivatives in the force equation when calculating magnetic forces (Fmag ∝ j × B ∝ (∇ × B) × B ∝ (∇ × (∇ × A)) × (∇ × A)), both leading to lower accuracy in the calculation. On the other hand, the advantages are a somewhat easier derivation of A for a given magnetic field and that there are no constraints on magnetic helicity using a vector potential. It would be definitely interesting to study the differences between simulations utilizing a vector potential and the Euler description, although it could be hard to overcome the problems related to numerical intricacies within a SPH implementation of the vector potential. Figure 18. Values of the different terms of the induction equation with time for the normal resolution Gadget simulation using Euler potentials. Upper plot: temporal evolution of the terms responsible for amplification/attenuation of the radial magnetic field (terms 1 to 5). Lower plot: Evolution of terms 6-10, responsible for the evolution of the toroidal magnetic field. Positive values imply amplification, and negative attenuation of the corresponding B-component. The non-axisymmetric terms are shown in red, the axisymmetric terms are shown in black.
The analysis of the different terms of the induction equation applied to our simulations clearly show that the non-axisymmetry of the velocity and magnetic field cannot be ignored in any consideration of the kinematic dynamo. There are two main processes leading to angular momentum transport and hence non-axisymmetry in spiral galaxies: Internal driving due to spiral structure and bar formation (the former considered in the presented paper) and external driving due to interaction with other galaxies. Simulations of interacting systems would therefore enrich our understandings even further on how large scale magnetic fields evolve due to large scale velocity fields.
Our simulations show only a weak amplification of the initial magnetic field. Observations of spirals galaxies at high redshifts suggest that their magnetic field strengths were at least as strong as the magnetic fields at the current epoch within few Gyrs of the Big Bang (Kronberg et al. 2008). Assuming initial strengths of order BIGM = 10 −9 G an am-plifying process should therefore account for four orders of magnitude of increase within few Gyrs in order to reach the observed values of ≈ 10µG. Since our simulations of a purely kinematic dynamo account at best for one order of magnitude, there is still need for a more complete scenario with additional subgrid physics. Such subgrid physics should comprise the α-effect due to turbulent gas motions below the resolution limit, estimated from local high-resolution MHD simulations and observations of turbulent motions in nearby galaxies. Hereby, potentially the most promising ansatz is the cosmic ray driven dynamo (Lesch & Hanasz 2003, Hanasz et al. 2005. Given the fact, that the presented simulations reveal the complete three dimensional velocity field to fully account for the large-scale structure of the magnetic field, we believe that N -body SPH together with sensible subgrid physics will be apt to test our understanding of the evolution of magnetic fields in spiral galaxies numerically.