MAGI: many-component galaxy initialiser

Providing initial conditions is an essential procedure for numerical simulations of galaxies. The initial conditions for idealised individual galaxies in N -body simulations should resemble observed galaxies and be dynamically stable for time scales much longer than their characteristic dynamical times. However, generating a galaxy model ab initio as a system in dynamical equilibrium is a diﬃcult task, since a galaxy contains several components, including a bulge, disc, and halo. Moreover, it is desirable that the initial-condition generator be fast and easy to use. We have now developed an initial-condition generator for galactic N -body simulations that satisﬁes these requirements. The developed generator adopts a distribution-function-based method, and it supports various kinds of density models, including custom-tabulated inputs and the presence of more than one disc. We tested the dynamical stability of systems generated by our code, representing early-and late-type galaxies, with N = 2,097,152 and 8,388,608 particles, respectively, and we found that the model galaxies maintain their initial distributions for at least 1 Gyr. The execution times required to generate the two models were 8 . 5 and 221 . 7 seconds, respectively, which is negligible compared to typical execution times for N -body simulations. The code is provided as open-source software and is publicly and freely available at https://bitbucket.org/ymiki/magi .


INTRODUCTION
N-body simulations are powerful tools for investigating the dynamical evolution of galaxies including, e.g., galactic mergers or the development of spiral arms.In order to run Nbody simulations, it is important to set up appropriate initial conditions.The initial system of N-body particles should be in dynamical equilibrium for at least the longest dynamical timescale in the problem.For example, the longest dynamical timescale may be the crossing time of a satellite galaxy in a galactic minor merger, or it may be the rotation time of a galactic disc in a simulation to follow the development of spiral arms.However, producing a system in dynamical equilibrium that represents a galaxy is difficult, since galaxies, in general, consist of several components, e.g., bulge, disc, and halo.The construction of initial conditions for N-body simulations is an on-going research topic, and many earlier studies have tackled this tough problem (Hernquist 1993;Kuijken & Dubinski 1995;Boily et al. 2001;Widrow et al. 2003;Widrow & Dubinski 2005;McMillan & Dehnen 2007; E-mail: ymiki@cc.u-tokyo.ac.jp (YM) Widrow et al. 2008;Vasiliev 2013;Perret et al. 2014;Yurin & Springel 2014;Vasiliev & Athanassoula 2015).
A model system of particles that represents an observed galaxy is useful not only in providing initial conditions for N-body simulations, but also in creating mock observations to fit observed datasets.Recent progress with observatories like Gaia (Lindegren et al. 2016) demands ever-better models to study galactic dynamics in much greater detail.Since Gaia provides observational data in a six-dimensional phase space, direct comparisons between particle models and the Gaia data can provide information about the phasespace distribution function (hereafter, DF) of the Milky Way galaxy.
Various density-profile models have been proposed and employed to cover the diversity of simulated and observed galaxies.For example, the Navarro-Frenk-White (hereafter, NFW) profile (Navarro et al. 1995(Navarro et al. , 1996)), the Moore profile (Fukushige & Makino 1997;Moore et al. 1998), and the Einasto profile (Einasto 1965;Navarro et al. 2004Navarro et al. , 2010) ) are often used to represent the density profiles for dark matter halos in cosmological N-body simulations with Λ cold dark matter.For the bulge component, the well-known de Vaucouleurs's law (de Vaucouleurs 1948) and the Sérsic profile (Sérsic 1963;Ciotti & Bertin 1999)-which is a generalisation of de Vaucouleurs's law-are frequently employed to fit the observed surface-density profiles.The Hernquist profile (Hernquist 1990), which resembles de Vaucouleurs's law, has an analytic DF.For the disc component, an exponential disc or a Sérsic profile are often used to fit the surface-density profiles of observed disc galaxies and to investigate their dynamical evolution in numerical studies.The presence of more than one disc component in a galaxy may be a general property of late-type galaxies.The Milky Way has both a thin and a thick disc (Jurić et al. 2008), and observations by Dalcanton & Bernstein (2002); Yoachim & Dalcanton (2006) have revealed multiple disc components (thin and thick discs) in most edge-on, late-type galaxies.
Many earlier studies have been dedicated to developing initial condition generators for N-body simulations; Hernquist (1993) pioneered this research field.The approach by Hernquist (1993) approximates the velocity distribution by a local Maxwellian and calculates the velocity dispersion using the Jeans equation.Boily et al. (2001) generalised the approach for axisymmetric bulges and halos.Springel et al. (2005) also extended the work by Hernquist (1993) to include a gaseous disc, and Perret et al. (2014) provided implementation as the open-source software DICE.However, Kazantzidis et al. (2004) showed that particle systems based on the local Maxwellian approximation are not always in dynamic equilibrium, which may lead to inadequate interpretations of results from N-body simulations that use such initial conditions.Currently, the standard approach is to employ DF-based implementations.A widely used example in the community is the software GALACTICS (Kuijken & Dubinski 1995;Widrow et al. 2003).This code generates a dynamically stable system that employs a lowered Evans model (Kuijken & Dubinski 1994) for the halo, a King model (Michie 1963;Michie & Bodenheimer 1963;King 1966) for the bulge, and an exponential disc.Alternatively, one can adopt orbit-based methods like the Schwarzschild method or an iterative method by Rodionov et al. (2009) to generate self-consistent dynamical-equilibrium systems of particles.Removing assumptions for DFs-e.g., spherical symmetry of the system and isotropy of velocity distributionextends the scope of applications to non-axisymmetric systems and makes the codes more flexible than DF-based implementations.However, the execution time is much longer than for DF-based methods.Vasiliev & Athanassoula (2015) reported that the execution time required to produce an axisymmetric disc galaxy model with 10 6 particles reached an hour or much longer with both the GALIC code (Yurin & Springel 2014) and the SMILE code (Vasiliev 2013;Vasiliev & Athanassoula 2015), while GALACTICS took only a few minutes.In this study, we have accordingly adopted a DFbased method rather than an orbit-based method to reduce the time-to-solution.
To generate a model galaxy that resembles an observed galaxy, the initial-condition generator needs to be flexible.The physical properties of the components represented by the N-body particles must be easily adjustable in order to facilitate investigations of the dependence of physical processes on the mass, size, or profile of a galaxy.For example, the particle system generated by GALACTICS is dynamically stable and is useful for numerical investigations of galactic dynamics.However, changing the particle system to repre-sent accurately an observed galaxy is cumbersome.Some input parameters in GALACTICS-e.g., the potential and the velocity dispersion at the centre-depend on the density profile of the system, making it difficult for the user to control basic quantities such as the total mass of the system.Furthermore, GALACTICS internally solves Poisson's equation and modifies the mass distribution toward a dynamically stable system.Therefore, a trial-and-error procedure is necessary to generate the desired particle system.Changing the input parameters from kinematic ones (such as the value of potential and the velocity dispersion) to structural ones (such as the total mass and scale length), and fixing the mass distributions of the spherical components, would alleviate the burden on the user performing these iterations.
The following are the key requirements we set for our initial-condition generator: (1) dynamical stability of the particle system generated, (2) ability to represent various kinds of density profiles, (3) ability to produce more than one disc component, and (4) simplicity and convenience of use.We have realised these requirements in our initial-condition generator for N-body simulations, called MAGI (MAny-component Galaxy Initialiser).The generator adopts a DF-based method, and it supports various kinds of density-model inputs, including a machine-readable tabular format and the presence of multiple disc components.The code is provided as open-source software.The organisation of this manuscript is as follows: Section 2 introduces methods and assumptions adopted in MAGI.The stability of the particle systems generated is tested in Section 3. In Section 4, we examine the validity of some numerical treatments and measure the execution time of MAGI.Section 5 summarises our work.

METHODS
In this section, we describe the DF-based implementation of our initial-condition generator, MAGI.Sections 2.1 and 2.2 describe how to calculate the DF and discuss the further techniques required to generate a spherical component for a given surface-density profile.One of the strong points of MAGI is a high degree of freedom in generating disc components.We describe the theories, assumptions and numerical techniques employed for the disc components in Section 2.3.Finally, Section 2.4 contains other implementation details.

Eddington's Formula
In this section, we derive the DF for a case in which the system is spherically symmetric.Consider a system consisting of N components with a volume-density profile ρ i (r), where the subscript i specifies the i-th component.If we assume an isotropic velocity distribution, Eddington's formula (Eddington 1916;Binney & Tremaine 2008) gives the DF f i of the i-th component as where and where E and Ψ are the relative energy per unit mass and the relative potential of the system, respectively.The last term in equation ( 1) corresponds to the density gradient at infinity, which vanishes.Now the second derivative of the mass density with respect to the relative potential can be rewritten as which is much easier to implement (Kazantzidis et al. 2006) compared to equation (1).We have calculated the DF using equations ( 1) and (3) for cases in which the density profile and its first and second derivatives are given for the individual components.Since M(r) is the total mass enclosed within radius r, the mass distributions of all components, including the disc components and the mass of the central massive black hole, must be given before calculating the DF.The mass profiles of the disc components are treated as spherically averaged profiles of the enclosed mass (i.e., the zeroth order multipole approximation).
MAGI provides various density profiles for generating models commonly adopted in astrophysics.The functional forms of the volume-density profiles supported by MAGI, together with their first and second derivatives with respect to radius, are summarised in Appendix A. The supported density profiles that have a central core are the Plummer profile (Plummer 1911), the Burkert profile (Burkert 1995), the King profile (the profile and derivatives are determined numerically using equations given in Appendix B), and the King profile which is given in empirical form (King 1962).Density profiles that have a central cusp are also implemented: the Hernquist profile, the NFW profile, the Moore profile, and the Einasto profile.MAGI also supports two broken-power-law density profiles, the double-powerlaw model (Hernquist 1990;Merritt et al. 2006) given by equation (A26) and the triple-power-law model described by equation (A29).If the system contains a central massive black hole, the black hole particle is placed at the centre-ofmass of the system with zero velocity.
MAGI also accepts density profiles given in machinereadable tabular form.Since evaluating the first and second derivatives of the profile by differencing the tabulated values can result in significant loss of accuracy, especially when the profile has a central core, MAGI interpolates the profile using a cubic-spline curve, which is then used to compute tables of the derivatives.Once the density profile and its derivatives are tabulated from the input density profile, the subsequent steps are identical to the case for a predefined density profile.
To generate systems of particles in dynamical equilibrium, in some cases we need to specify an explicit cutoff radius r c .If a cutoff radius is specified, MAGI multiplies the profiles with a complementary-error-function-based smoother of the form where ∆ c is the smoothing scale.

Abel Transformation
In most cases, observations provide surface-density profiles rather than the volume-density profiles discussed in Section 2.1.MAGI also supports the input of a surface-density profile.De-projection from the input surface-density profile Σ(R) to a volume-density profile ρ(r) is performed by an Abel transformation (cf., Binney & Tremaine 2008): if ρ(r) drops faster than r −1 .MAGI constructs the volumedensity profile as a numerical table by applying the Abel transformation numerically.The derivatives of the density profile are given by cubic-spline interpolation, as in the case for which the input volume-density profile is a machinereadable table .Once the density profile and its derivatives are obtained in tabular form from the input surface-density profile, the subsequent steps are the same as in the case for which the input profile is a volume-density profile.MAGI also supports the Sérsic profile described in Appendix C, which includes de Vaucouleurs's law, and the surface-density profiles for this case are provided as machine-readable tables.

Disc Components
The disc is one of the most important characteristics of a late-type galaxy.Observations by Dalcanton & Bernstein (2002) and by Yoachim & Dalcanton (2006) established that most edge-on, late-type galaxies have multiple disc components (both thin and thick discs).Our initial-condition generator must therefore be capable of generating several disc components in dynamical equilibrium.MAGI is the first initial-condition generator to support more than one disc components.
We have extended the approach adopted in GALACTICS (Kuijken & Dubinski 1995;Widrow et al. 2003) to provide particle distributions that include multiple disc components.In GALACTICS, the disc component is assumed to have an exponential surface-density profile, with an isothermal profile in the vertical direction.The potential-density pair is derived by solving Poisson's equation using a spherical multipole expansion with a simple approximation for highorder terms.MAGI supports Sérsic profiles (the exponential profile corresponds to a Sérsic profile with Sérsic index unity), and surface-density profiles are specified in machinereadable tabular form for the disc components.Supporting Sérsic profiles having various Sérsic indices is essential for representing a large variety of observed disc components.For example, Kelvin et al. (2012) reported that a typical Sérsic index is around 0.5 − 2 for disc components obtained from the GAMA (Galaxy And Mass Assembly) database.We have also assumed that every disc component has an isothermal profile in the vertical direction, given by where z d is the scale height of the disc component.The vertical profile of a disc component thus represents the density field as a function of the potential field, while Poisson's equation gives the potential field as a function of the density field.
Therefore, we derive the potential-density pair for the superposition of all components numerically by iterating the following two procedures until convergence is obtained: (1) Determine the vertical profile of the individual disc components from the potential field, and (2) solve Poisson's equation to determine the potential field from the superposition of the density fields of all components.We solve Poisson's equation using the BiCGSTAB method preconditioned with ILU(0) (van der Vorst 1992; Itoh et al. 2012).We use a nested grid to discretize the density and potential fields, and the number of levels-typically ten-is automatically determined.The outer boundary condition at R = R max or z = z max for level L = 0 (the coarsest grid) is derived from Binney & Tremaine (2008) by assuming a constant disc-scale-height.The outer boundary condition for level L > 0 is set by the potential field in the coarser level L − 1.For the inner boundary conditions (R = 0 or z = 0), Φ(±R, ±z) = Φ(R, z) is set from the symmetry about R = 0 and z = 0.The number of grid points needed to determine the potential-density pair is 256 in the R direction and 64 in the z direction, respectively.We adopted the full multigrid method proposed by Press & Teukolsky (1991).In order to determine the velocity structure, we adopted the Schwarzschild DF (Schwarzschild & Villiger 1907;Binney & Tremaine 2008): where L z is the z-component of the specific angular momentum, and v R , v φ , and v z are the velocities in the horizontal, azimuthal, and vertical directions, respectively.Also, R g is the guiding-centre radius, v c (R) = RΩ is the circular velocity, where Ω is the circular frequency, and σ R and σ z are, respectively, the velocity dispersions in the horizontal and vertical directions.The remaining quantity, γ, is defined as 2Ω g /κ, where Ω g is the circular frequency at the guidingcentre radius, and κ is the epicycle frequency.The underlying assumption behind the Schwarzschild DF is the epicycle approximation, assuming small velocity dispersion.The radial profile of the velocity dispersion in the horizontal direction, σ R (R), is another problem.GALACTICS assumes an exponential profile where R s is the disc scale length, even though GALACTICS adopts the epicycle approximation, assuming that the orbit is nearly circular.In MAGI, the velocity dispersion σ p in the azimuthal direction is equal to min under the epicycle approximation.The velocity dispersion is thus smaller than that in GALACTICS, especially in the central region, ensuring consistency with the epicycle approximation.
In GALACTICS, the circular velocity at the particle position is determined by the derivative of the potential in the equatorial plane of the disc component; in other words, the rotation velocities of the particles do not depend on z, while the potential does.This results in too fast a rotation to maintain the initial distribution, especially in a thick disc # set cutoff (1) or not (0) 50.0 5.0 # cut off radius and width at large z, and a relaxation phase is required before dynamical stability is reached.MAGI improves on this by calculating the circular velocity from the gradient of the potential at the particle location.
Finally, in order to generate a realistic distribution, we need a recipe for suppressing a needle-like structure that occasionally forms near the rotation axis of the disc.When the scale height z d of the disc is not sufficiently small, compared to the scale length r sph of the spherical components, the density decreases gradually in the vertical direction for R r sph .This produces a needle-like structure near the rotation axis of the disc component (see Fig. 5b); it is an artefact resulting from the isothermal profile assumed in the vertical direction.

We have introduced a new quantity
and we replace the scale height This ensures that the vertical density profile vanishes for |z| ≥ min (16z d , r c ) and damps the needle-like structure.

Configuration and Implementation
In this subsection we provide some details about the structure and use of MAGI.As pointed out in Section 1, input parameters such as the mass and scale length are preferable to kinematic ones such as the central potential or velocity dispersion.In MAGI, the input parameter sets for each component are limited to the mass M, the scale length r s , and other dimensionless parameters that are specific to the adopted models, such as the Sérsic index n or the dimensionless King parameter W 0 .Listings 1-3 are samples of the input files that show the expected format and parameters.The text following each "#" symbol are brief comments describing the parameters; they do not appear in the actual input files.In Listing 1, the first two integers specify the system of units to be used in the simulation and the number of components.The configuration of each component follows: the index of the model, the corresponding parameter file, whether to specify the number of particles (1) or not (0), and the specified number of particles.For the central massive black hole, the number of particles must be unity.Also, the configurations of the disc components must come after those of the spherical components.
In each parameter file, the mass M and scale radius r s of the components appear first.If the unit ID = 1 (as specified in Listing 1), the units of mass and length are M and kpc, respectively.Since MAGI internally converts the system of units from astrophysical units (like M , kpc, and km s −1 ) to code units, users need not perform unit conversions.For the spherical components, model-specific parameters (if they exist) follow; for example, the non-dimensional parameter W 0 at the centre of the King model (Listing 2).In the case of the disc components (e.g., Listing 3), the scale height z d is an additional key parameter.The following two parameters pertain to the velocity-dispersion profile.The first sets the radial velocity dispersion at the centre, σ R,0 , and the second sets the scaling parameter f .A negative value of σ R,0 indicates that the value is identical to the vertical velocitydispersion at the centre.The first parameter is used only when the radial velocity-dispersion model given by equation (8), which is that used in GALACTICS, is adopted.By default, MAGI adopts the radial velocity-dispersion model given by equation ( 9) in the previous subsection.The next parameter is optional and allows the introduction of a retrograde component in the disc.The remaining parameters are related to the explicit cutoff of the density distribution.An integer indicates whether an explicit cutoff of the density profile is required (1) or not (0).When a cutoff is set (specified as 1), the cutoff radius r c and its smoothing scale ∆ c are passed to the software by specifying them in the next line.
MAGI employs the inverse-function method and the rejection method, respectively, to determine the particle positions and velocities.Pseudo-random numbers are generated by a SIMD (single instruction, multiple data)-oriented Fast Mersenne Twister of period 2 19937 − 1 (SFMT 1 1.5.1 by Saito & Matsumoto 2008) with a jump function 2 .The jump function guarantees the independence of multiple series of random numbers generated in parallel, while just changing the initial values of a random number generator cannot strictly ensure independence.In the current version of MAGI, the program is parallelized using OpenMP and shifts the initial state for SFMT in each OpenMP thread by 10 100 steps.The particles' positions and velocities are shifted component-bycomponent to set the centre-of-mass of the system at the origin of the coordinate system and remove the bulk motion.

RESULTS
This section examines the dynamical stability of particle distributions generated by MAGI.Section 3.1 considers the stability of a spherically symmetric system like an early-type galaxy or the bulge of disc galaxy.Section 3.2 studies the stability of an axisymmetric system having two disc components with different scale heights.The initial condition is generated on a workstation (Intel Xeon E5-2640v3, 16 cores, 2.60 GHz, DDR4-2133, 64GB, gcc 4.8.5).The code is compiled with the following options for performance optimisation: -O3 -ffast-math -funroll-loops -march=nativefopenmp.The code is run via the numactl command with the --localalloc option.Time evolution is calculated using the octree code GOTHIC (Miki & Umemura 2017), which runs on NVIDIA GeForce GTX TITAN X with nvcc 8.0.44.

Stability of a Spherically Symmetric System
We first benchmark MAGI with a spherically symmetric model in order to confirm the effectiveness of the Eddington formula.We chose the mass ratios of the spherical components in our model from well-established observed correlations.The Magorrian relation (Magorrian et al. 1998;Marconi & Hunt 2003) gives the correspondence between the galactic bulge and central massive black hole (MBH), where the BH mass is around 0.2% of that of the spheroidal component.A similar observational relation between the darkmatter halo and the central MBH (Ferrarese 2002) implies that the mass of the central black hole is around 10 7 M , if the dark-matter halo mass is 10 12 M .In addition to the typical components of early-type galaxies, we added a stellar halo resembling that of the Milky Way in the model galaxy.Recent observations of the stellar halo of the Milky Way suggest that its power-law index is around −3 at r 50 kpc and −5 in the outer halo (Keller et al. 2008;Akhter et al. 2012).The masses of the stellar halos of nearby disc galaxies with masses comparable to that of the Milky Way are 1−6×10 9 M (Harmsen et al. 2017).In summary, the model galaxy is a superposition of an Einasto halo (M = 10 12 M , r s = 10 kpc, α = 0.2, r c = 200 kpc, ∆ c = 10 kpc) with 2,084,644 particles; a stellar halo given as a triple-power-law model (M = 10 9 M , r in = 3 kpc, r out = 50 kpc, α = 0, β = 1, γ = 3, δ = 1, = 5, r c = 150 kpc, ∆ c = 20 kpc) with 2,084 particles; a stellar component obeying the de Vaucouleurs's law (M = 5 × 10 9 M , r s = 2 kpc, r c = 10 kpc, ∆ c = 1 kpc) with 10,423 particles; and an MBH particle with a mass of 10 7 M .The elapsed time for MAGI to generate this model galaxy was 8.5 seconds.
We have calculated the time evolution for this model galaxy over 1 Gyr using GOTHIC with an accuracy-controlling parameter ∆ acc = 2 −7 = 7.8125 × 10 −3 and a Plummer softening length of 15.625 pc.The resulting radial profile is shown in Fig. 1.The figure compares the volume-density profile and the enclosed-mass profile for each component at t = 1 Gyr (symbols) with the input radial profiles (curves).The figure clearly shows that the given density profile is stable over the entire domain after an integration time of 1 Gyr, which is approximately 70 times the free-fall time at r = 2 kpc.In other words, MAGI successfully generates the model galaxy as a system in dynamical equilibrium.

Stability of Multiple Disc Components
The next test involves models of late-type galaxies.In contrast to early-type galaxies, late-type galaxies often harbor   The bulge-to-total ratio of the model galaxy is 0.25, which is consistent with observed late-type galaxies (Oohama et al. 2009).Other physical properties of the model galaxy are shown in Appendix D. The elapsed time required for MAGI to generate this model galaxy was 221.7 seconds.We calculated the time evolution of this model galaxy over 1 Gyr using GOTHIC, with an accuracy-controlling parameter ∆ acc = 2 −8 = 3.90625×10 −3 and a Plummer softening length of 15.625 pc. Figure 2 compares the surface-density profile of each component at t = 1 Gyr (symbols) with the corresponding input profile (curves).The figure shows very satisfactory agreement between the imposed surface-density profile and the surface-density profile at t = 1 Gyr, which is around 220 times the free-fall time for the bulge component and around 10 times the rotation periods of the discs at their scale radii.Maintaining the respective disc thicknesses is an essential requirement for the disc components.Figure 3 compares the disc heights at t = 1 Gyr (the bold lines) with those set in the initial conditions (the regular lines), for the thin and thick discs separately (red and black lines, respectively).The disc heights are evaluated as the root-mean-square of the disc particle heights with respect to the z = 0 plane.Figure 3 shows that the thicknesses of both discs are almost unchanged after 1 Gyr (around ten-fold the rotation time at R = R s ).The drop in the thickness of the thick-disc component toward the centre is due to the scale-height reduction prescribed by equation ( 10) to remove the artificial needlelike structure (see also Fig. 5).In a series of convergence tests, it has turned out that simulation runs with a larger number of disc particles, N d , result in almost identical distributions.On the other hand, runs with a smaller number of halo particles, N h , significantly thicken the disc thicknesses.We found that N h 3 × 10 6 and N d 10 5 is a sufficient condition to stabilize the disc components in a viewpoint of reducing disc heating by halo particles.It is noteworthy that we were able to use the output from MAGI directly as the initial condition for the N-body simulation.In other words, for these tests, MAGI did not require any relaxation procedure before starting the simulations.

Accuracy of MAGI
To validate the accuracy of the numerically generated DFs, we compared the DFs generated by MAGI with those given by analytic formulas.Figure 4 shows the DFs for Hernquist and Plummer spheres, together with the differences originating from the input data format for MAGI.Here r out (E) are the radii satisfying the relation which corresponds to the apocentre for a given specific energy E. Table 1 lists the physical properties of the profiles.The tables of density profiles have 128 bins equally spaced in the logarithm from 10 −5 r s to 50r s .Figure 4 indicates that the DFs derived from the density profiles in machine-readable tabular format match well with those given by the functional form.The relative error is below ∼2% in most regions, irrespective of whether the central density profile is a cusp or a core.That is to say, MAGI properly generates DFs for the density profile given in machine-readable tabular form.
Figs. 4(a) and 4(b) also compare the DFs from MAGI with analytic expressions (equations A15 and A4).Since the DFs in MAGI have a density cutoff at a finite radius, they drop sharply around r c , while the analytic counterparts decrease continuously.The agreement between the DFs, except near the cutoff radius, confirms that the DF generator in MAGI works properly and with sufficient accuracy.Figure 5 confirms that the additional treatment using equation ( 10) removes the artificial needle-like structure protruding from the disc (Fig. 5a).As shown in Fig. 5(b), the needle-like structure appears only in the thick disc and not in the thin disc.The presence of the needle-like structure is a natural result of the assumption of an isothermal profile in the vertical direction.Because the thickness of the thick-disc component is comparable to the scale length of the bulge component, the vertical density profile near the rotation axis drops very slowly.
The upper limit to the scale heights of disc components that work properly in MAGI is another concern.Since the extension to thick discs implemented in MAGI relies partly on the workaround described in Section 2.3, there are limitations in the present method.To examine these limits, we have tested the dynamical stability of the thick-disc component by varying its scale height in the late-type galaxy model described in Section 3.2 from R s /10 = 0.3 kpc to 2R s = 6 kpc. Figure 6 displays surface-density distribution maps of thick discs with various scale heights at t = 0 Gyr and 1 Gyr.The figure shows that the surface-density distribution is stable for 1 Gyr, even when the scale height exceeds the scale length R s = 3 kpc.This indicates that the thick disc is stable, at least for z d 2R s .We note that the morphology of the thick disc is similar at z d R s = 3 kpc because the scale height of the disc in the central region was reduced to r c /16 ∼ 3 kpc to remove the needle-like structure.If the removal of this artefact is incomplete, then infalling components from the needle-like structure generate a shell structure after having been scattered by the bulge.The masses of the needle-like structure and the resulting shell structure are negligibly small; however, numerical artefacts should be explicitly removed by capping the disc scale height according to equation (10).or in machine-readable tabular form (the black solid curves) are shown, together with analytic expressions for the DFs (the blue dashed curves).Lower panel: The relative errors in the DFs based on tabular forms, as compared with the functional forms, are plotted as functions of the normalised apoapsis at the given energy.The left panels present the results for the Hernquist profile, while the right panels correspond to the Plummer profile.

Execution Time
The execution time of the software is an important performance metric.Figure 7 plots the medians of 10 measurements of the execution time required for MAGI to generate the late-type galaxy model tested in Section 3.2.The problem size-that is, the number of N-body particles N tot -is varied from 2 17 = 131,072 to 2 29 = 536,870,912.Above this limit, the problem size does not fit into the DRAM in the measurement environment.The total execution time is around 220-230 seconds for N tot 10 8 .It increases for N tot 10 8 and reaches 420 seconds for N tot = 2 29 .For all values of N tot , the dominant process is the preparation of the disc components (filled squares), which includes the time required to generate the potential-density pairs and calculate the spherically averaged density profiles, circular velocity, and velocity dispersion in the vertical direction.The execution time for this procedure is independent of N tot and is 212.0 seconds on average.If there are no disc components, the total execution time is reduced to 10 seconds for N tot 10 7 .The time required to prepare the spherical components (filled circles)-which includes setting the radial density profiles, integrating them and calculating the DF using Eddington's formula)-and the time required to write miscellaneous files for further analysis (the open triangles) are also independent of N tot : they are about 3.4 and 4.6 seconds on average, respectively.On the other hand, the execution time required to distribute the N-body particles (the filled diamonds and solid line) and to write the particle data (the open diamonds and solid line) increase with N tot , as expected.Above N tot ∼ 2×10 8 , the writing speed for particle data becomes slow.Protocol switching in HDF5 or the bandwidth of the file system (4 WD40EFRXs managed as RAID 5 by mdadm v3.4 and XFS format) might be the reason, but investigating this is beyond the scope of the present study.In summary, the measured performance results show that: (1) generating galaxies with disc component(s) takes 4 minutes for N tot 10 8 , (2) generating galax- ies without disc components take 10 seconds for N tot 10 7 , and (3) the elapsed time increases for N tot 10 7 .Since the most time-consuming procedure in MAGI is the iterative calculation of the potential-density pairs for the disc components, a more sophisticated preconditioner than ILU(0), or performance optimisations of sparse matrix-vector multiplication would accelerate this part.
In the default configuration of MAGI, pseudo randomnumbers are generated by SFMT using the jump function to obtain the benefits of OpenMP thread parallelization.Without the jump function, we observe approximately a ten-fold deceleration in the particle-distribution process (the filled diamonds and dotted line in Fig. 7).The speed-up rate delivered by the jump function reaches a factor of 10 for N tot 10 8 , which corresponds to a parallel efficiency of 63%.GSL also provides a pseudo-random number generator based on the Mersenne Twister of period 2 19937 − 1, and it can also be used in MAGI.Since GSL does not provide an equivalent to the jump function, the independence of multiple series of the random numbers it generates is not strictly guaranteed, and OpenMP parallelization must be switched off during the particle-distribution process.The performance of the random numbers generator provided in GSL is 4.6 times slower than SFMT, and the particle-distribution process requires a ∼ 30% longer execution time compared to SFMT without the jump function.
The HDF5 container supports on-the-fly file compression with Szip6 2.1, and the resulting file size of the particle data is 46% of the original.However, we have not employed Szip compression in the default operation of MAGI, since it is 10-30 times slower than without Szip compression in most cases.

Strengths and Versatility of MAGI
Initial conditions for N-body simulations are also useful for fitting observed datasets.Typically, when fitting an observed dataset, one must generate many mass distribution models and pick an optimal one.Three features of MAGI that makes this procedure very efficient are: (1) MAGI supports various kinds of density profiles including machine-readable tables; (2) the input parameters adopted in MAGI are structural ones (total mass and scale length), which helps users to adjust the generated mass-distribution models to reproduce the observed datasets; and (3) the short execution time for MAGI (10 seconds for early-type galaxy models and 4 minutes for late-type galaxy models) enables users to produce many trial models to fit the observed galaxies.Gaia provides sixdimensional phase-space information about the Milky Way's stellar component (Lindegren et al. 2016), and direct comparisons between the Gaia data and particle systems generated by MAGI will be useful for modelling the Milky Way's DF.
The flexibility of MAGI also enables one to use idealised simulations to study any kind of physics in detail.Most initial-condition generators are limited to astrophysical profiles such as the NFW profile or de Vaucouleurs's law.On the other hand, MAGI accommodates arbitrary spherically symmetric density profiles, so long as the profile is twice differentiable.Among other things, this is sometimes useful for testing predictions of analytic models.For example, Michikoshi & Kokubo (2014, 2016) investigated the development of stellar spirals in disc galaxies using linear analysis and local N-body simulations, and they estimated the number of spiral arms, assuming a constant ratio of disc mass to the total mass of the galaxy.In order to verify their predictions, one would need to carry out global N-body simulations with carefully controlled initial conditions.The flexibility of MAGI facilitates the generation of initial conditions for such highly idealised problems.
Hydrodynamic simulations based on a particle method (e.g., smoothed particle hydrodynamics) are widely employed to investigate in detail the formation and evolution processes of galaxies.Extensions to include gaseous components would be meaningful and not difficult, if we assume hydrostatic equilibrium for the gaseous component(s).

CONCLUSIONS
We have developed an initial-condition generator for Nbody simulations, called MAGI.The implementation relies on DFs to determine the velocity distribution of the target system.In order to calculate the DFs, we have exploited Eddington's formula for spherically symmetric components assuming an isotropic velocity distribution, and the Schwarzschild DF with an isothermal profile in the vertical direction for disc components.Since MAGI supports various kinds of volume-density profiles and surface-density profiles including a machine-readable tabular format, MAGI is very flexible in producing model galaxies.One of the strong points of MAGI is that the software can generate galaxy models including multiple disc components, appropriate for latetype galaxies like the Milky Way.Input parameters for each component are limited to the scale length, mass, and other model-specific dimensionless parameters to control the output intuitively.
We tested the dynamical stability of the generated particle systems using N-body simulations.The results show that the systems retain the original distribution for at least 1 Gyr both early-type and late-type galaxy models.The ex- where ρ 1 is the scale density, and σ is a parameter having the dimension of velocity (but it is not the velocity dispersion).Integrating equation (B1) over velocity space gives the density profile: where W ≡ Ψ/σ 2 is the dimensionless potential.It is determined by numerically solving Poisson's equation where x denotes the radius normalized by the King radius r 0 .When solving equation (B3), two boundary conditions are necessary to specify the solution: the non-dimensional King parameter W 0 at the centre and the requirement of a central core ( dW/dx = 0 at the centre).Using the central density ρ 0 , the King radius r 0 is given by The first and second derivatives of the density profile can then be expressed as (B6) Here, we already have dW/dx from the solution of equation (B3), and d 2 W dx 2 is obtained by solving the differential equation

APPENDIX C: SURFACE-DENSITY PROFILES
In this section, we list the surface-density profiles provided by MAGI.Here, X is defined as X ≡ R/R s , where R s is the scale length.The surface-density profile for the Sérsic model (Sérsic 1963) is given by where n is the Sérsic index and b n is a dimensionless scale factor.The scale factor b n is internally calculated using an asymptotic formula by Ciotti & Bertin (1999) that works well, at least in the range 1 ≤ n ≤ 10.Sérsic profiles with n = 1 and 4 correspond to an exponential profile and de Vaucouleurs's law (de Vaucouleurs 1948), respectively.

APPENDIX D: RADIAL PROFILES OF DISC COMPONENTS
In this section, we show radial profiles of the late-type galaxy model investigated in §3.2.curve of the model galaxy.Figure D2 exhibits the radial profiles of Toomre's Q-value for the disc components.Since the velocity dispersion of the disc components is capped in the central region to ensure that the epicycle approximation remains valid, the Q-value decreases toward the centre and becomes smaller than unity.Figure D3  profiles of the velocity dispersion for the disc components at t = 1 Gyr (symbols) with the input radial profiles (curves).The figure shows that the given profiles of σ z are stable over the entire domain after an integration time of 1 Gyr.On the other hand, the radial profiles of σ R evolve, especially in the central region.This indicates that the epicycle approximation is not satisfied in the central regions of late-type galaxies.The σ R -profile of the thin-disc component has better stability in other domains, which suggests that the thin disc is more consistent with the epicycle approximation than the thick-disc component.
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Listing 1 .
Format of the input file.1 # select the system of units 4 Format of the parameter file (King model).5.0 e +10 # mass of the component 1.0 # scale radius of the component 5.0 # model specific parameter ( s ) 0 # set cutoff (1) or not (0) Listing 3. Format of the parameter file (Sérsic disc).2.5 e +10 # mass of the component 5.0 # scale radius of the component 2.0 # Sersic index ( only for Sersic disc ) 1.0 # scale height of the component -1.0 0.1 # parameters for velocity dispersion 0.0 # retrograde fraction [0.0 , 1.0] 1

Figure 1 .
Figure1.Volume-density profile (upper panel) and enclosedmass profile (lower panel) of a model galaxy that represents an early-type galaxy with 2 21 = 2,097,152 particles.Symbols and curves show the profile at t = 1 Gyr and the input profile, respectively.Every pair of symbols and curves represents different components of the model galaxy: the dark-matter halo with an Einasto profile (the black filled circles and the solid curve), the stellar halo with a triple-power-law profile (the blue filled squares and the dashed curve), and the bulge with a profile obeying the de Vaucouleurs's law (the red open circles and the dotted curve).

Figure 2 .
Figure 2. Surface-density profile of a model galaxy that uses 2 23 = 8,388,608 particles to represent a late-type galaxy.Symbols and curves show the profile at t = 1 Gyr and the input profile, respectively.Each pair of symbols and curves represents a different component of the model galaxy: the dark-matter halo represented with an NFW profile (the black filled circles and the solid curve), the bulge as a King model (the red open circles and the dotted curve), the thick disc with a Sérsic profile (the blue filled squares and the dashed curve), and the thin disc as an exponential profile (the magenta open squares and the dot-dashed curve).

Figure 3 .
Figure3.Radial profiles of disc thicknesses evaluated as the root-mean-square of disc particle heights about the z = 0 plane.The black lines show the thick-disc component, while the red lines correspond to the thin-disc component.The bold and regular lines represent t = 1 Gyr and the initial condition, respectively.

Figure 4 .
Figure 4. Comparison of DFs.Upper panel: DFs generated by MAGI using density profiles in functional form (the red dotted curves) or in machine-readable tabular form (the black solid curves) are shown, together with analytic expressions for the DFs (the blue dashed curves).Lower panel: The relative errors in the DFs based on tabular forms, as compared with the functional forms, are plotted as functions of the normalised apoapsis at the given energy.The left panels present the results for the Hernquist profile, while the right panels correspond to the Plummer profile.

Figure 5 .
Figure 5. Removal of the artificial needle-like structure.The black, red, and blue dots show the distributions of the N -body particles representing the bulge, thick-disc, and thin-disc components, respectively.The left panel shows the particle distributions with the modification given by equation (10) and the right panel the distributions without any modification.

Figure 7 .
Figure7.Execution time for the code as a function of the total number of N -body particles.The crosses and solid line show the total running time for the code (with the use of the jump function), while the other markers and lines show the breakdown by different stages of the calculation.The filled and open diamonds and solid lines are the times required to distribute the N -body particles according to the prescribed DF and the time required to dump the output file, respectively.The filled diamonds and dotted line also show the execution time required to distribute the N -body particles, but for the case without the jump function.The circles and dashed line, and the squares and dot-dashed line, respectively, represent the time necessary to calculate the DF for the spherical components and the disc components.The triangles and solid line correspond to the time required to dump the files (the DF, density distribution, and various profiles such as the enclosed mass and potential).

FigureFigure D1 .Figure D2 .
Figure D1.Rotation curve for the late-type galaxy model of §3.2.The black solid, blue dot-dashed, magenta dotted, and green dashed curves show the contributions from the dark-matter halo, bulge, thick disc, and thin disc, respectively.The solid red curve is the sum of all components.

Figure D3 .
Figure D3.Radial profiles of the velocity dispersion for the disc components.Symbols and curves show the profile at t = 1 Gyr and the input profile, respectively.The black symbols and lines show the thick-disc component, while the red ones correspond to the thin-disc component.The circles and solid curves represent the velocity dispersion in the R direction, while the squares and dotted curves correspond to the velocity dispersion in the z direction.