Abstract
We have modified our particle-by-particle adaptation of the made-to-measure (M2M) method, with the aim of modelling the Galactic disc from upcoming Galactic stellar survey data. In our new particle-by-particle M2M algorithm, primal, the observables of the target system are compared with those of the model galaxy at the position of the target stars, i.e. particles. The mass of the model particles are adjusted to reproduce the observables of the target system, and the gravitational potential is automatically adjusted by the changing mass of the particles. This paper builds upon our previous work, introducing likelihood-based velocity constraints in primal. In this paper, we apply primal to barred disc galaxies created by an N-body simulation in a known dark matter potential, with no error in the observables. This paper demonstrates that primal can recover the radial profiles of the surface density, velocity dispersion in the radial and perpendicular directions, and the rotational velocity of the target discs, along with the apparent bar structure and pattern speed of the bar, especially when the reference frame is adjusted so that the bar angle of the target galaxy is aligned to that of the model galaxy at every timestep.
INTRODUCTION
The made-to-measure (M2M) method developed by Syer & Tremaine (1996) has seen increasing use in the past few years, including multiple papers involving nmagic (e.g. de Lorenzi et al. 2007, 2008; Das et al. 2011; Morganti & Gerhard 2012), and the series of three papers (Long & Mao 2010, 2012; Long et al. 2013) based on the M2M algorithm described in Long & Mao (2010). This increasing level of interest is now highlighting the potential of the M2M method to be used in many different applications, from probing the dark matter halo of NGC 4649 (Das et al. 2011) to tailoring N-body initial conditions (Dehnen 2009). While previous work has primarily focused on theoretical models and external galaxies, Bissantz, Debattista & Gerhard (2004) and Long et al. (2013) have applied their versions of M2M to the Milky Way which is also our eventual goal.
The Milky Way is known to be a barred spiral (Spergel & Blitz 1990; Weinberg 1992; Binney, Gerhard & Spergel 1997) but there is still disagreement over whether the bar is formed of a single structure, or if it is comprised of a separate long flat bar in addition to the barred bulge (e.g. Martinez-Valpuesta & Gerhard 2011; Athanassoula 2012; González-Fernández et al. 2012). The bar angle, defined as an offset between the axis of the bar and the Sun–Galactic Centre line is also still debated. For example Dwek et al. (1995) suggested a bar angle of 20° ± 10°, Alcock et al. (2000) found 15°, Rodriguez-Fernandez & Combes (2008) found 20°–35°, Minchev, Nordhaus & Quillen (2011) found 20°–45° and Green et al. (2011) found 45°. There is also still discussion of the pattern speed of the bar, Ωp. Among others Fux (1999) calculated a value of Ωp ≈ 50 km s− 1 kpc− 1, Dehnen (1999) calculated Ωp = 53 ± 3 km s− 1 kpc− 1, Rodriguez-Fernandez & Combes (2008) calculated Ωp = 30–40 km s− 1 kpc− 1, Wang et al. (2012) concluded Ωp ≈ 60 km s− 1 kpc− 1 and Gerhard (2011) provided a review of previous work, concluding a likely pattern speed for the bar of Ωp ≈ 50–60 km s− 1 kpc− 1.
Bissantz et al. (2004) were the first to apply the original M2M algorithm from Syer & Tremaine (1996) to the Milky Way, which was an important first test of M2M on the Milky Way, matching observed bar kinematics in several fields. The M2M algorithm has undergone significant improvements since then much more observational data have become available, leaving the door open for more up to date M2M modelling of the Milky Way.
Long et al. (2013) have taken the next step, applying their M2M algorithm from Long & Mao (2010) to observed radial velocity data from the Bulge RAdial Velocity Assay (BRAVA)1 (Rich et al. 2007; Kunder et al. 2012) and the particle mass distribution of an N-body barred galaxy model (Shen et al. 2010) that has been built to reproduce the Milky Way disc and earlier BRAVA data. Long et al. (2013) used the gravitational potential calculated from the N-body model of Shen et al. (2010). Then, they rotate this potential with a fixed pattern speed and assumed bar angle from the line of the Galactic Centre and the observer, i.e. the position of the Sun. They ran many models with different pattern speeds and bar angles and explored the models which best fit their observables. They found their best model recovered the bar angle and pattern speed of the Shen et al. (2010) N-body model, and reproduced the mean radial velocity and radial velocity dispersion of the BRAVA data very well.
Although Long et al. (2013) have taken an important step forwards, there will be much room for improvement for applying M2M to the future observational data. The European Space Agency's upcoming Gaia mission along with ground-based surveys, e.g. the Panoramic Survey Telescope And Rapid Response System (PanStarrs, e.g. Kaiser et al. 2010), the Visible and Infra-red Survey Telescope for Astronomy (VISTA, e.g. Minniti et al. 2009), the Large Synoptic Survey Telescope (LSST, e.g. Ivezić et al. 2008), the Sloan Extension for Galactic Understanding and Exploration (SEGUE, e.g. Yanny et al. 2009), the Apache Point Observatory Galactic Evolution Experiment (APOGEE, e.g. Allende Prieto et al. 2008), Subaru-PFS (e.g Ellis et al. 2012), The Radial Velocity Experiment (RAVE, e.g. Steinmetz et al. 2006), the Large sky Area Multi-Object fibre Spectroscopic Telescope, LAMOST Experiment for Galactic Understanding and Explanation (LEGUE, e.g. Deng et al. 2012) and the Gaia-ESO public spectrographic survey (e.g. Gilmore et al. 2012) will provide us with an unprecedentedly large amount of data about the Milky Way, which we can use as observational constraints for our M2M algorithm. As such we have made modifications to our algorithm with this goal in mind.
In Hunt & Kawata (2013), hereafter Paper 1, we developed a particle-by-particle M2M algorithm now called primal. Because the Galactic stellar-survey data, such as the ones Gaia will produce, are in the form of the position and velocity of individual stars, primal is designed to compare the observables at the position of each star, i.e. not binned data as in previous M2M modelling. Another major difference between primal and other M2M algorithms is that the gravitational potential is calculated via self-gravity of the model particles. The potential is thus altered by the changing particle masses induced by the M2M algorithm. In Paper 1, we apply primal to the target system of a smooth axisymmetric disc created by N-body simulations, and demonstrate that primal can reproduce the density and velocity profiles of the target system well, even when starting from a disc whose scalelength is different from the target system.
We are encouraged by the success of Paper 1 and intend to build upon it. In this paper, we apply primal to barred disc galaxies again generated by N-body simulations with gcd+ (Kawata & Gibson 2003; Kawata et al. 2013). We introduce a new form of velocity observable constraints as described in de Lorenzi et al. (2008), based on the likelihood function as described in Romanowsky & Kochanek (2001). We also introduce a rotating reference frame in a similar, although not identical fashion to Long et al. (2013).
We use target systems whose information is known without any error. Ultimately, we wish to apply primal to real observational data, where the information will be provided for a limited region of the sky, with a more complicated selection and error function due to the dust extinction, crowding and stellar populations. However, we think that in the development stages it is important to test the algorithm against the ideal target. In this paper, we demonstrate the successful application of primal to the barred galaxy targets, and this is a significant step forward to modelling the Milky Way with M2M.
This paper is organized as follows. Section 2.1 describes the M2M methodology of our original M2M algorithm as shown in Paper 1. Sections 2.2 and 2.3 describe the new improvements we have made to the algorithm, and Section 2.4 describes the setup of our target systems. Section 3 shows the performance of our updated method for recreating the target disc systems. In Section 4, we provide a summary of this work.
THE M2M ALGORITHM: primal
The M2M algorithm works by calculating observable properties from the model and the target, and then adapting particle masses such that the properties of the model reproduce those of the target. The target can be in the form of a distribution function, an existing simulation or real observational data. The model can be a test particle simulation in an assumed fixed or adaptive potential, or a self-gravity N-body model.
A particle-by-particle M2M
Maximum likelihood for velocity constraints
We have again performed a parameter search to determine differences in the likelihood-based velocity constraints, as we did in Paper 1. There are four important parameters, ϵ, α, ζ and μ, which must be calibrated for M2M. ϵ provides the balance between the speed of convergence, and the smoothness of the process. Note that ϵ = ϵ′ϵ′′, where ϵ′′ is defined by equation (11). We have chosen ϵ′ = 0.1 as an appropriate balance between accuracy and simulation time. With more computing power available to us we would consider running a lower value of ϵ. However if ϵ′ ≪ 0.1, it is possible that the model might not converge as the mass change is too slow. The choice of α, which controls the strength of the temporal smoothing, should depend upon the choice of ϵ (α ≥ 2ϵ). We find that our modelling is not overly sensitive to α as long as the condition α ≥ 2ϵ is met and we set α = 2.0 in this paper. We set ζ = 1 as mentioned before. μ controls the strength of the regularization and is essential in reducing the oscillation in particle masses and ensuring smooth convergence. We discuss the importance of μ in much greater detail in Paper 1. In this paper, we adopt μ = 105. All different models presented in Section 3 use this same parameter set, and have not been individually tailored to the target or model in question. This demonstrates the robustness of the method.
Rotating reference frames
In Paper 1, it was sufficient to use a fixed reference frame as we were investigating smooth axisymmetric discs. However, if the target has some non-axisymmetric structure, such as a bar, the target bar angle is fixed, but the bar of the model rotates in the fixed reference frame. For example, if there is a bar, we expect the density and kinematics to be very different at the different azimuth angles at a fixed radius. Then, if the bar of the model is not aligned with the target bar, the observables of the model are evaluated in the different dynamical states from the target observables. Hence, if the bar of the model keeps rotating in the fixed frame, the model particles receive the different constraints from the target depending on the bar angle at each timestep, and the model never settles to the solution. This is discussed in Section 3.
Long et al. (2013) have proposed using a reference frame with a fixed bar angle, and comparing multiple simulations with different bar angles to find the best fit. This was trivial for their model, because they used a fixed shape of the bar potential and rotating with a fixed bar pattern speed, Ωp. We however have not assumed any pattern speed prior to the beginning of our simulations, nor have we placed any explicit constraints on it. Instead, we start with a smooth disc as the initial condition, allowing the pattern speed to evolve with the model galaxy due to self-gravity, and compare Ωp for the model and target galaxies at the end of the run.
We therefore calculate the angle of the bar in our target, and the angle of the bar in the model at each step. Then we rotate the model to match the bar angle of the target for the purposes of calculating the observables in the same reference frame. It is our hope that this method will allow the pattern speed to be recovered along with the density and velocity profiles. When applying this to the Milky Way we will not know the exact bar angle. But here, we assume that the bar angle is known for our first step of modelling the bar. We call this reference frame change the rotating reference frame hereafter. In Section 3, we present a comparison of our method with and without this rotating reference frame, and also present the results from cases where we have chosen an incorrect bar angle.
Target system setup
N-body target parameters. M200 is the mass of the halo, Md is the mass of the disc, c is the concentration parameter, zd is the scale height, |$\sigma _r^2/\sigma _z^2$| is the ratio between radial and tangential velocity dispersion and Ωt, p is the pattern speed of the bar, measured at 2 Gyr.
| Target | M200 (M⊙) | Md (M⊙) | c | zd (kpc) | |$\sigma _r^2/\sigma _z^2$| | Ωt, p (kms− 1 kpc− 1) | Notes |
|---|---|---|---|---|---|---|---|
| I | 1.75 × 1012 | 3.0 × 1010 | 20.0 | 0.35 | 9.0 | N/A | Smooth disc |
| II | 2.0 × 1012 | 5.0 × 1010 | 9.0 | 0.3 | 2.0 | 27.5 | |
| III | 1.5 × 1012 | 5.0 × 1010 | 7.0 | 0.3 | 2.0 | 31.7 | |
| IV | 1.75 × 1012 | 5.0 × 1010 | 9.0 | 0.3 | 2.0 | 28.9 |
| Target | M200 (M⊙) | Md (M⊙) | c | zd (kpc) | |$\sigma _r^2/\sigma _z^2$| | Ωt, p (kms− 1 kpc− 1) | Notes |
|---|---|---|---|---|---|---|---|
| I | 1.75 × 1012 | 3.0 × 1010 | 20.0 | 0.35 | 9.0 | N/A | Smooth disc |
| II | 2.0 × 1012 | 5.0 × 1010 | 9.0 | 0.3 | 2.0 | 27.5 | |
| III | 1.5 × 1012 | 5.0 × 1010 | 7.0 | 0.3 | 2.0 | 31.7 | |
| IV | 1.75 × 1012 | 5.0 × 1010 | 9.0 | 0.3 | 2.0 | 28.9 |
N-body target parameters. M200 is the mass of the halo, Md is the mass of the disc, c is the concentration parameter, zd is the scale height, |$\sigma _r^2/\sigma _z^2$| is the ratio between radial and tangential velocity dispersion and Ωt, p is the pattern speed of the bar, measured at 2 Gyr.
| Target | M200 (M⊙) | Md (M⊙) | c | zd (kpc) | |$\sigma _r^2/\sigma _z^2$| | Ωt, p (kms− 1 kpc− 1) | Notes |
|---|---|---|---|---|---|---|---|
| I | 1.75 × 1012 | 3.0 × 1010 | 20.0 | 0.35 | 9.0 | N/A | Smooth disc |
| II | 2.0 × 1012 | 5.0 × 1010 | 9.0 | 0.3 | 2.0 | 27.5 | |
| III | 1.5 × 1012 | 5.0 × 1010 | 7.0 | 0.3 | 2.0 | 31.7 | |
| IV | 1.75 × 1012 | 5.0 × 1010 | 9.0 | 0.3 | 2.0 | 28.9 |
| Target | M200 (M⊙) | Md (M⊙) | c | zd (kpc) | |$\sigma _r^2/\sigma _z^2$| | Ωt, p (kms− 1 kpc− 1) | Notes |
|---|---|---|---|---|---|---|---|
| I | 1.75 × 1012 | 3.0 × 1010 | 20.0 | 0.35 | 9.0 | N/A | Smooth disc |
| II | 2.0 × 1012 | 5.0 × 1010 | 9.0 | 0.3 | 2.0 | 27.5 | |
| III | 1.5 × 1012 | 5.0 × 1010 | 7.0 | 0.3 | 2.0 | 31.7 | |
| IV | 1.75 × 1012 | 5.0 × 1010 | 9.0 | 0.3 | 2.0 | 28.9 |
As mentioned above, in this initial stage of development, we assume that the dark matter halo potential is known and there is no other external potential such as the bulge or stellar halo. We use the same number of particles, 106, and the same initial dark matter halo and disc parameters for the model and target galaxies, except for the disc scalelength: Rd = 2 kpc for the models and Rt, d = 3 kpc for the targets.
RESULTS
M2M model results at the final timestep. Ωt, p is the target pattern speed, Ωp is the model pattern speed, |$\chi ^2_{\rho }$| is a measure of accuracy of the density, |$\mathcal {L}_{r,z,{\rm rot}}$| are the likelihood values for the radial, tangential and rotational velocity.
| Model | Target | Ωp (kms− 1 kpc− 1) | |$\chi ^2_{\rho }$| | |$-\mathcal {L}_r/10^6$| | |$-\mathcal {L}_z/10^6$| | |$-\mathcal {L}_{{\rm rot}}/10^6$| | Notes |
|---|---|---|---|---|---|---|---|
| A | I | N/A | 0.123 | 6.06 | 3.61 | 5.27 | Smooth disc |
| B | II | 27.9 | 0.254 | 4.72 | 3.62 | 4.43 | |
| C | III | 30.4 | 0.235 | 4.96 | 4.40 | 5.28 | |
| D | IV | 28.3 | 0.189 | 5.14 | 5.15 | 5.02 | |
| E | II | 27.3 | 0.230 | 4.77 | 3.65 | 4.47 | Partial data |
| F | II | 23.6 | 0.276 | 5.77 | 3.91 | 5.21 | Bar angle −30o |
| G | II | 28.0 | 0.250 | 4.87 | 3.67 | 4.55 | Bar angle −10o |
| H | II | 24.3 | 0.21–0.49 | 6.82–9.99 | 4.59–5.48 | 6.63–8.77 | No rotating frame |
| Model | Target | Ωp (kms− 1 kpc− 1) | |$\chi ^2_{\rho }$| | |$-\mathcal {L}_r/10^6$| | |$-\mathcal {L}_z/10^6$| | |$-\mathcal {L}_{{\rm rot}}/10^6$| | Notes |
|---|---|---|---|---|---|---|---|
| A | I | N/A | 0.123 | 6.06 | 3.61 | 5.27 | Smooth disc |
| B | II | 27.9 | 0.254 | 4.72 | 3.62 | 4.43 | |
| C | III | 30.4 | 0.235 | 4.96 | 4.40 | 5.28 | |
| D | IV | 28.3 | 0.189 | 5.14 | 5.15 | 5.02 | |
| E | II | 27.3 | 0.230 | 4.77 | 3.65 | 4.47 | Partial data |
| F | II | 23.6 | 0.276 | 5.77 | 3.91 | 5.21 | Bar angle −30o |
| G | II | 28.0 | 0.250 | 4.87 | 3.67 | 4.55 | Bar angle −10o |
| H | II | 24.3 | 0.21–0.49 | 6.82–9.99 | 4.59–5.48 | 6.63–8.77 | No rotating frame |
M2M model results at the final timestep. Ωt, p is the target pattern speed, Ωp is the model pattern speed, |$\chi ^2_{\rho }$| is a measure of accuracy of the density, |$\mathcal {L}_{r,z,{\rm rot}}$| are the likelihood values for the radial, tangential and rotational velocity.
| Model | Target | Ωp (kms− 1 kpc− 1) | |$\chi ^2_{\rho }$| | |$-\mathcal {L}_r/10^6$| | |$-\mathcal {L}_z/10^6$| | |$-\mathcal {L}_{{\rm rot}}/10^6$| | Notes |
|---|---|---|---|---|---|---|---|
| A | I | N/A | 0.123 | 6.06 | 3.61 | 5.27 | Smooth disc |
| B | II | 27.9 | 0.254 | 4.72 | 3.62 | 4.43 | |
| C | III | 30.4 | 0.235 | 4.96 | 4.40 | 5.28 | |
| D | IV | 28.3 | 0.189 | 5.14 | 5.15 | 5.02 | |
| E | II | 27.3 | 0.230 | 4.77 | 3.65 | 4.47 | Partial data |
| F | II | 23.6 | 0.276 | 5.77 | 3.91 | 5.21 | Bar angle −30o |
| G | II | 28.0 | 0.250 | 4.87 | 3.67 | 4.55 | Bar angle −10o |
| H | II | 24.3 | 0.21–0.49 | 6.82–9.99 | 4.59–5.48 | 6.63–8.77 | No rotating frame |
| Model | Target | Ωp (kms− 1 kpc− 1) | |$\chi ^2_{\rho }$| | |$-\mathcal {L}_r/10^6$| | |$-\mathcal {L}_z/10^6$| | |$-\mathcal {L}_{{\rm rot}}/10^6$| | Notes |
|---|---|---|---|---|---|---|---|
| A | I | N/A | 0.123 | 6.06 | 3.61 | 5.27 | Smooth disc |
| B | II | 27.9 | 0.254 | 4.72 | 3.62 | 4.43 | |
| C | III | 30.4 | 0.235 | 4.96 | 4.40 | 5.28 | |
| D | IV | 28.3 | 0.189 | 5.14 | 5.15 | 5.02 | |
| E | II | 27.3 | 0.230 | 4.77 | 3.65 | 4.47 | Partial data |
| F | II | 23.6 | 0.276 | 5.77 | 3.91 | 5.21 | Bar angle −30o |
| G | II | 28.0 | 0.250 | 4.87 | 3.67 | 4.55 | Bar angle −10o |
| H | II | 24.3 | 0.21–0.49 | 6.82–9.99 | 4.59–5.48 | 6.63–8.77 | No rotating frame |
Smooth disc
First, we demonstrate that the newly introduced likelihood velocity constraints can reproduce the smooth featureless disc target used in Paper 1. Model A applies primal to Target I, which was the target used in Paper 1, but using a larger number of particles. Note the high value of |$\sigma _r^2/\sigma _z^2$| of Target I in Table 1 was used to deliberately suppress structure formation. Fig. 1 shows that the radial profiles of the density, the radial velocity dispersion, the vertical velocity dispersion and the mean rotation velocity for the target galaxy, the initial galaxy and the final output after primal is applied. The figure demonstrates that primal with the likelihood-based velocity constraints equally or even more accurately reproduces the target galaxies compared to our old version of the particle-by-particle M2M (Paper 1). However, a quantitative comparison between the old and new version is not the main focus of this paper. As discussed above, we introduced the likelihood-based velocity constraints, because we can compare the velocity more directly and also introduce different errors for individual velocity components and individual particles. Therefore, the likelihood-based velocity constraints are a necessary update, and a comparison with the old version is not an important issue. Note that the properties shown in Fig. 1 are not explicitly constrained by primal. As discussed previously in Paper 1, it is interesting to note that although our particle-by-particle M2M uses only the first moment of the velocity components as observables, because primal tries to reproduce the velocity of individual particles, the velocity distribution becomes close to the target and therefore the velocity dispersion can be reproduced as well.
Initial (red dotted), final (green dashed) and target (black solid) density profile (upper), radial velocity dispersion (upper middle), vertical velocity dispersion (lower middle) and rotation velocity (lower) for Model A with Target I.
Barred disc galaxies
In this section, we present the results of Models B, C and D, where we apply the same parameter set for primal to model three different target barred galaxies. Target II is a barred disc galaxy showing faint spiral structure. Fig. 2 shows the face-on and edge-on views of Target II and the final state of Model B (two top panels). The final model reproduces the bar feature very well. The observables are only constrained within 10 kpc of the galactocentric radius and hence the areas outside this radius are reproduced with less accuracy.
Face-on (left) and edge-on (right) density maps of Target II, Models B, E, F, G and H, (from top to bottom) plotted for comparison. The white line indicates the angle of the bar, rotated for comparison. The density scale is the same for all panels.
Fig. 3 shows the radial profiles of the surface density, radial and tangential velocity dispersion, and the mean rotation velocity for the target and the final model compared to the initial model. As in Paper 1 and Model A, these radial profiles are not directly constrained by primal, but are reproduced remarkably well. Fig. 3 shows a substantial increase in radial velocity dispersion and a corresponding decrease in mean rotational velocity from the initial to the final model. We believe that this is due to heating from the bar which leads to an excellent agreement with the velocity dispersion of the target.
Same as Fig. 1, but for Model B with Target II.
The pattern speed of the bar, Ωp, is also reproduced very well, as shown in Tables 1 and 2. We calculate the pattern speed of the bar, Ωp, by calculating the change in angle of the bar between timesteps, divided by the difference in time between the steps. We take Ωp to be the mean value from the final 10 steps. We found that the bar pattern speed of the model is Ωp = 27.9 km s−1 which is close to that of the target, Ωt, p = 27.5 km s−1. This is probably due to our self-consistent calculation of the gravitational potential, because once the mass distribution and kinematic properties of the target disc are reproduced, a bar with a similar shape and pattern speed to those of the target is expected to develop. This is certainly helped by our use of a known, fixed dark matter halo potential. We are pleased to see a spiral arm developing in the model, which looks similar to the one seen in the target.
Model C applies primal to Target III, which is also a barred disc galaxy, but with a smaller bar than Target II, and a boxy and peanut shaped bulge (e.g. Pfenniger 1984; Athanassoula & Misiriotis 2002; Debattista et al. 2005; Bureau et al. 2006; Saito et al. 2011), as can be seen in the top panel of Fig. 4. Rather surprisingly primal reproduces the boxy structure of the target as shown in the bottom panel of Fig. 4. Fig. 5 shows the radial profiles for Model C. We see a slight inaccuracy in the inner 2 kpc of the radial velocity dispersion, and also in the rotational velocity in the inner 4 kpc, which corresponds roughly with the length of the bar. In addition, σz is systematically higher than the target at all radii. As such the bar pattern speed is not as well reproduced as with Model B, with Ωp = 30.4 km s−1 compared to Ωt, p = 31.7 km s−1. However, we still think that this is a reasonably good recovery of the target, and it is encouraging for further developments to apply primal to more complicated observational data.
Same as Fig. 2, but for Model C with Target III.
Same as Fig. 1, but for Model C with Target III.
Model D takes Target IV which is morphologically similar to Target III, with a small bar and boxy peanut feature, as can be seen in Fig. 6. We see a slightly larger bar in the model than in the target. Fig. 7 shows slight inaccuracies in the recovery of the radial and vertical velocity dispersion and mean rotational velocity in the inner 3 kpc roughly consistent with the radius of the bar. However, the pattern speed is still recovered well with Ωp = 28.3 km s−1 for the final model compared to the target of Ωt, p = 28.9 km s−1.
Same as Fig. 2, but for Model D with Target IV.
Same as Fig. 1, but for Model D with Target IV.
Working with partial data
Even with the huge amount of data returned by Gaia and related stellar surveys, due to our position within the Milky Way's disc, we will not even come close to having a complete data set of the disc stars. Therefore, it is important to make sure our method is still applicable when we do not have access to the complete picture of the disc. Our previous models have used all data within 10 kpc from the galactic centre. However, Model E was performed with a simple selection function restricting the observable volume to a 10 kpc sphere around a point in the plane 8 kpc from the galactic centre, i.e. at (x, y, z) = (8, 0, 0) in Fig. 2, roughly emulating Gaia's observable area, while ignoring effects such as extinction and errors. This is merely the first step towards using primal with realistic data.
The third panel of Fig. 2 shows the face-on and edge-on view of Model E, which has a similar bar to the target (top panel), with a hint of a spiral arm in the lower-left quadrant matching the one visible in the target. Fig. 8 shows that an excellent agreement of the final model with the target radial profiles is still obtained with the restricted data set. This is an improvement on Paper 1, where we saw loss of accuracy when the observable field was restricted. We believe that this is helped by both the likelihood form of velocity observable and the higher resolution with which the simulations have been carried out. The bar pattern speed is recovered very well with Ωp = 27.3 km s−1 compared to the target of Ωt, p = 27.5 km s−1. This shows the ability of primal to produce reasonable results when supplied with a partial data set of the disc particles. However, we are aware that this selection function is crude and the next stage of our work will deal with more realistic selection functions and expected observational errors.
Same as Fig. 1, but for Model E with Target II, performed with a partial data set.
Same as Fig. 1, but for Model F with Target II, performed with the assumed angle of the bar offset from 30° compared to the real value.
Working with an incorrect bar angle
As mentioned in Section 1, the bar angle of the Milky Way is still debated. Ultimately, we aim to recover the dynamical state of the Milky Way with primal from the future stellar survey data, and recovering the bar angle is also one of our targets. In the previous sections, we assumed that the bar angle of the target is known and we align the bar of the model galaxy to that of the target at every timestep to evaluate the observables. If we do not know the bar angle of the target, like with the Milky Way, we could try different bar angles and hope that the models with the lowest χ2 and/or the maximum likelihood values recover the bar angle of the target, which is the strategy taken by Long et al. (2013). In this section, we examine the effects of running primal with an incorrect bar angle. Models F and G are performed with the bar angles deliberately set to be incorrect by 30° and 10°, respectively. In this section, we again use all data within 10 kpc from the galactic centre.
Model F has been performed while assuming that the bar angle is 30° less than the real angle of the target. The fourth panel of Fig. 2 shows a poor reproduction of the target bar morphology in Model F, which is significantly shorter than that of the target (top panel). We also see no evidence of the spiral structure seen in other cases. Fig. 9 shows the radial profiles for Model F. There is a discrepancy in the inner 4 kpc of the model compared to the target in both the density profile and the radial velocity dispersion. This is in agreement with the weaker bar shown in Fig. 2. The average rotational velocity is also lower across the disc. This is also reflected in the final pattern speed of Ωp = 23.6 km s−1 compared to the target of Ωt, p = 27.5 km s−1. However, in the real Milky Way case, we cannot know the correct profiles or the bar pattern speed in advance. On the other hand, we can evaluate the goodness of fit by |$\chi ^2_{\rho }$| or the values of the likelihood, |$\mathcal {L}_v$|. In Table 2, Model F shows significantly worse values of |$\chi ^2_{\rho }$| and |$\mathcal {L}_v$| than those of Model B which assumes the correct bar angle. Therefore, we should be able to tell easily if the bar angle is off by 30°, at least in this simple target case.
Same as Fig. 1, but for Model G with Target II, performed with the assumed angle of the bar offset from 10° compared to the real value.
Model G has been performed while assuming that the bar angle is 10° less than the bar angle of the target. The fifth panel of Fig. 2 shows a barred disc which is morphologically similar to the target (top panel). The bar is reproduced well whereas the spiral structure is barely visible. Fig. 10 shows the radial profiles for Model G, which again reproduces very well those of the target. The bar pattern speed is still well recovered with Ωp = 28.0 km s−1 compared to the target of Ωt, p = 27.5 km s−1. In Table 2, Model G shows similar values of |$\chi ^2_{\rho }$| and |$\mathcal {L}_v$| to those of Model B, although the velocity likelihood values are slightly worse. These results may indicate that primal does not have the power to determine the bar angle within 10° accuracy, but can recover it with better than 30° accuracy. However, our ultimate target is much more complicated than this ideal target, and at this stage we do not explore further the expected accuracy of recovering the correct bar angle for this ideal target. At least we demonstrate that with this type of exercise we can examine how accurately the dynamical model, such as primal, can recover the bar angle. In our future study, we will construct more realistic mock observational data from N-body barred simulated discs and ‘train’ primal to recover the bar angle as accurately as possible, and finally evaluate the expected accuracy of our recovered bar angle using the comparison demonstrated in this section.
The importance of the rotating reference frame
In this section, we show a brief comparison between the resulting models with and without the rotating reference frame. Model H was performed under identical conditions to Model B, but without the reference frame corrections detailed in Section 2.3. The bottom panel of Fig. 2 shows that the resulting disc contains a less prominent bar, and no evidence of spiral structure in a similar fashion to Model F. Fig. 11 shows the radial profiles of Model H. In the inner 4 kpc region, the radial density and radial velocity profiles are lower for the model than for the target. The average rotation velocity is lower than that of the target across the whole disc. The pattern speed is also too low with Ωp = 24.3 km s−1 compared to the target of Ωt, p = 27.5 km s−1. Fig. 12 shows a comparison of the evolution of the |$\chi ^2_{\rho }$| of the density between Model B and Model H. The |$\chi ^2_{\rho }$| in Model H experiences periodic oscillations in time with the bar rotation which are not seen in Model B. This lack of a smooth model convergence along with the poor accuracy on the recovered radial profiles shows the importance of having a rotating reference frame.
Same as Fig. 1, but for Model H with Target II, which is performed without the rotating reference frame.
Time evolution of |$\chi ^2_{\rho }$| for density for Model B (black line) compared to Model H (red line).
SUMMARY
We have demonstrated that our updated particle-by-particle M2M algorithm, primal, can recover a target disc system with a bar, including boxy/peanut features, in a known dark matter halo potential. In primal, the observables are compared with the model at the position of the target particles. The mass of the model particles are adjusted to reproduce the target observables, and the gravitational potential is calculated self-consistently from the model particle mass distribution. We have introduced the likelihood-based velocity constraints to primal, which allows us to compare the velocity of the target particle more directly than the smoothed velocity field used in our previous algorithm. To apply this method to a barred disc, we evaluate at every timestep the density and velocity likelihood after the reference frame of the model disc has been corrected, so that the bar of the model is always aligned with the bar of the target. Our fiducial model recovers the radial profiles of the surface density, the radial and vertical velocity dispersion and the mean rotation velocity of the target system very well. In addition, because of our self-gravity implementation of M2M, we can reproduce the bar morphology and pattern speed. We have demonstrated that primal performs well even when the observables are restricted to within a sphere of radius 10 kpc around a point in the disc plane 8 kpc from the centre.
We admit that these applications are still simplified cases. The excellent recovery of the properties of the target galaxies is not surprising because we have not applied any observational errors or selection functions. Our ultimate goal is to further improve primal to be applicable to the future stellar survey data, including the Gaia data. While Gaia will return an unprecedentedly large amount of data, for approximately one billion stars, the accuracy of this data will be highly variable due to distance, extinction, location in the sky etc. In a forthcoming paper, we will apply primal to more realistic mock observational data from N-body simulations, taking into account the observational errors and selection functions. This paper also assumes a known, static, spherical dark matter halo potential for simplicity. In reality the dark matter halo remains very much unknown and yet has a significant effect on the dynamics of its inner galaxy. Thus, we need to explore different dark matter halo potentials and to consider the possibility of using a model which includes a live halo.
We remain optimistic that we can continue to improve primal, and develop a unique tool to recover the dynamical properties of the Milky Way from the future large-scale stellar survey data.
We thank the anonymous referee for their time and effort. We thank Victor Debattista for recommending their likelihood-based velocity constraints. We gratefully acknowledge the support of the UK's Science & Technology Facilities Council (STFC Grant ST/H00260X/1). The calculations for this paper were performed on Cray XT4 at Center for Computational Astrophysics, CfCA, of the National Astronomical Observatory of Japan and the DiRAC Facilities, Legion and COSMOS, jointly funded by the STFC and the Large Facilities Capital Fund of BIS. The authors also thank the support of the STFC-funded Miracle and COSMOS consortium (part of the DiRAC facility) in providing access to the UCL Legion and Cambridge COSMOS High Performance Computing Facilities. We additionally acknowledge the support of UCL's Research Computing team with the use of the Legion facility. This work was carried out, in part, through the Gaia Research for European Astronomy Training (GREAT-ITN) network. The research leading to these results has received funding from the European Union Seventh Framework Programme ([FP7/2007-2013] under grant agreement number 264895). HM is supported by the Natural Sciences and Engineering Research Council of Canada, and the Canada Research Chair program. HM is thankful to the Department of Physics and Astronomy, University of Victoria, for its hospitality.
