Foivos I Diakogiannis

Geraint F Lewis

Rodrigo A Ibata

Magda Guglielmo

Mark I Wilkinson

Chris Power

Reliable mass calculation in spherical gravitating systems

Keywords: galaxies, dwarf -galaxies, kinematics and dynamics -techniques, radial velocities -methods, statistical -galaxies, statistics

We present an innovative approach to the methodology of dynamical modelling, allowing practical reconstruction of the underlying dark matter mass without assuming both the density and anisotropy functions. With this, the mass-anisotropy degeneracy is reduced to simple model inference, incorporating the uncertainties inherent with observational data, statistically circumventing the mass-anisotropy degeneracy in spherical collisionless systems. We also tackle the inadequacy that the Jeans method of moments has on small datasets, with the aid of Generative Adversarial Networks: we leverage the power of artificial intelligence to reconstruct non-parametrically the projected line-of-sight velocity distribution. We show with realistic numerical simulations of dwarf spheroidal galaxies that we can distinguish between competing dark matter distributions and recover the anisotropy and mass profile of the system.

INTRODUCTION

Whilst dark matter represents the dominant mass component of the universe, its true nature remains elusive. Astrophysical probes of the properties of dark matter in large galaxies and galaxy clusters are typically hampered by the complexities of baryonic physics, and the complex coupling of the properties of kinematic tracers and the underlying form of the gravitational potential.

In recent years, considerable focus has been given to dwarf spheroidal galaxies in the local universe. With a stellar mass of ∼ 10 7 M , these are seen to be both devoid of gas, limiting the impact of baryonic astrophysics, and sufficiently simple to allow the determination of the gravitational potential of the dominant dark matter component from the stellar motions. However, traditional approaches of determining the distribution of dark matter in dSphs are limited by both the influence of the observational uncertainties and the mathematical complexity of deriving the properties of the dark matter.

One such approach, the Schwarzschild (1979) method, attempts to determine the underlying dark matter distribution through the reconstruction of the observed luminosity and kinematic properties of a galaxy using a library of * E-mail: foivos.diakogiannis@data61.csiro.au precomputed orbits in trial potentials. Via the appropriate weighting of the components of the library for a particular mass model, the optimal fit to the data can be recovered and the mass determined. However, the computational aspects of the Schwarzschild method makes implementation highly impractical. Building a high-resolution orbit library to survey the likelihood of millions of mass models is currently computationally prohibitive.

Other approaches are based upon the Jeans equation [START_REF] Binney | [END_REF], which relates the properties of kinematic tracers to the form of the gravitational potential. When applying the Jeans equation, there are two key ingredients, the distribution of dark matter, and a velocity anisotropy, β, which describes the relationship between radial and tangential orbits within the structure. In established approaches, it is typical to assume a functional form for the dark matter distribution, such as a Navarro-Frenk-White [START_REF] Merritt | Dynamics and Evolution of Galactic Nuclei[END_REF] or a [START_REF] Plummer | [END_REF] profile, and a functional form for β, optimizing the parameters of both based upon the observational data. Given the mathematical form of the Jeans equation, however, the resultant determination of the mass depends upon the assumed form for β, with various combinations of the adopted mass profile and β providing equally acceptable fits to the data. Known as the "mass-anisotropy degeneracy" (hereafter MAD), this is generally accepted as a fundamental limitation of Jeans-based approaches (Merrifield & Kent 1990, see also Read & Steger 2017).

In this contribution, we present a new approach to address the MAD in the Jeans formalism, relying upon a parametrised functional form, known as a B-Spline, to account for the implicit relationship between the dark matter profile and the velocity anisotropy. In this latest version of the JEAnS (Diakogiannis et al. 2017) approach, the tight-JEAnS (hereafter t-JEAnS), we represent both the unknown radial and tangential velocity dispersions as Bsplines. Then, we allow the data to give them the correct geometric shape. In this way, we avoid having to assume the functional form of all, but one, of the unknown functions used in the modelling process. Then, even with competing dark matter models that have equal numbers of unknown coefficients, we end up with statistical fits of different quality. The key point is that by demanding that these curves be as simple as possible, i.e. that they are represented by a minimal number of variables, competing dark matter density models give different qualitative fits to the data. This eventually allows us to statistically discriminate between competing mass models and thus transform the MAD to a mere model inference problem. For the case of small datasets (of the order of 1000 tracer stars), we use Generative Adversarial Networks (hereafter GANs, Goodfellow et al. 2014) to reconstruct non-parametrically the underlying projected line-of-sight (LOS) velocity distribution. With this, we augment artificially the data to arbitrarily large numbers, and obtain reliable estimates for the moments of the LOS velocity distribution with an unprecedented density of points. The combination of t-JEAnS modelling with the GANs for artificial data augmentation is a powerful approach for reliable mass estimates.

In Section 2 we present a short review of the Jeans mass modelling method. In Section 3 we give the details of the datasets we used as well as the preprocessing method we followed. In Section 4 we give a detailed description of the t-JEAnS algorithm. In Section 5 we present our findings and in Section 6 we discuss the reasons behind the efficiency of the t-JEAnS . Finally in Section 7 we present our concluding remarks.

A REVIEW OF THE JEANS MODELLING METHODOLOGY

In this section we present an overview and analysis of the established [START_REF] Binney | Galactic Dynamics: Second Edition[END_REF]) methodology of Jeans modelling. We continue by providing a proof for the uniqueness of the anisotropy profile upon assuming a specific functional form for the mass density profiles of stars, ρ , and dark matter, (hereafter DM), ρ•.

The Jeans modelling approach subject to the assumption of spherical symmetry is fully contained in the following two equations:

- dΦ dr = 1 ρ d dr ρ σ 2 rr + 2 r β(r)σ 2 rr (1) σ 2 los (R) = 2 Σ (R) r vir R 1 -β(r) R 2 r 2 rρ σ 2 rr √ r 2 -R 2 dr. (2)
Here, Φ is the total potential of the system, ρ the stellar tracer density, σ 2 rr the radial velocity dispersion, Σ is the projected tracer surface density, σ 2 los is the observed line-ofsight velocity dispersion, β is the anisotropy profile defined by β(r) = 1 -σ 2 tt /(2σ 2 rr ), and R and r are, respectively, the projected and 3D distance radii from the centre of the system. Although the integral in Eq. ( 2) usually has infinity as its upper bound, here we define rvir as the distance in which the DM mass density, ρ•, profile falls to approximately ρ•(rvir) ≈ 200ρcrit. For all practical purposes, this is a useful numerical approximation that does not alter our findings. With the exception of the observed LOS velocity dispersion, σ 2 los , and the projected tracer density profile, Σ , all remaining functions (ρ , ρ•, σ 2 rr , β) are unknown and need to be determined from the data. Therefore, the system is underdetermined 1 . In practice we can make a very good approximation to the functional form of the tracer density profile, ρ given deep photometry of the dSph, and we are thus left with three unknown functions, {ρ•(r), σ 2 rr (r), β(r)} in a system of two equations.

Although there are variations 2 to the general methodology, the common established [START_REF] Binney | Galactic Dynamics: Second Edition[END_REF] starting point to solving this system of coupled integrodifferential equations with respect to the unknowns ρ•(r), β(r) and σ 2 rr (r), is to assume parametric functional forms for the DM mass density, ρ•, and the anisotropy profile, β(r). In an iterative approach (assuming for simplicity we have full knowledge of the tracer profile, ρ ), one proposes a set of values for the parameters that define ρ• and β, then solves the differential Equation (1) with respect to 3 σ 2 rr and substitutes the result in Eq (2). The validity of the numerical values of the parameters that define ρ• and β is tested by comparing the model σ 2 los with the observables. This iterative process is performed until some convergence criterion is met. The rationale behind this approach is that when we consider parametric forms for ρ• and β, the system becomes overdetermined (since Equations ( 2) and (1) are evaluated in various distinct locations, ri, Rj) and thus a solution exists.

It needs to be emphasized though that once we make an assumption for the parametric form of one of the three unknown functions, {ρ•, σ 2 rr , β}, the system of two equations with (the remaining) two unknowns is closed. That is, the remaining two functions can be fully determined without the need for their parametric representation. There exist published (Binney & Mamon 1982;Solanes & Salvador-Sole 1990;[START_REF] Chan | [END_REF]Mamon & Boué 2010) exact solutions to the system of these equations (termed inversion techniques) that make a parametric assumption for only one of the three unknown functions. These prove that it is an un-1 One though needs to be precise in the definition of the number of "unknowns". Usually, we make assumptions for the functional form of these unknown functions that depend on some parameters. It is the number of these parameters that define the necessary number of equations to close the system. Then, for either exact (numerical solutions) or overdetermined systems (statistical fitting), we evaluate each of the Equations ( 1) and (2) in a set of distinct locations, r i , R j that are equal or greater in numbers to the number of unknown parameters. 2 These include, e.g. using higher moments ( Lokas & Mamon 2003) of σ 2 los , or different assumptions on the distribution function of the system, i.e. different penalty functions when comparing the LOS velocity dispersion with observables.

3 Subject to the boundary condition limr→r vir σ 2 rr ≈ 0.

necessary assumption to assume two of the three unknown functions in parametric form. Usually, assuming more parametric forms than necessary, increases the uncertainty in the model parameters, thus making the distinction between competing mass models even more difficult1 . It is insightful to separate the process of solving the system of coupled integrodifferential Equations ( 1) and ( 2), in two distinct approaches: the exact numerical solution of the equations to perfect noiseless data and the statistical fitting to noisy data. Clearly, all conclusions we can draw from knowledge gained in exact solutions of the system of Jeans equations can be transferred to the case of statistical fitting, while the converse is not always true. In the following, we focus on the exact numerical solution.

Uniqueness of the anisotropy profile for a given mass model

In this section, we provide a theorem that upon making an assumption for the functional form of the tracer and DM mass densities, ρ , ρ•, and the LOS dispersion, σ2 los there exists a unique anisotropy profile, β. For our purposes, we consider we have full knowledge of the above-mentioned functions, ρ , ρ• and σ 2 los . We solve Eq (1) with respect to β(r) and substitute it under the integral sign of Eq (2). Then we end up with a single integrodifferential equation (subject to the virial boundary condition limr→r vir σ 2 rr (r) ≈ 0), namely

σ 2 los (R) = 2 Σ (R) r vir R KA d(ρ σ 2 rr ) d r + ρ d Φ d r + KBρ σ 2 rr d r (3)
where KA and KB are kernel functions defined by:

KA(r, R) = R 2 √ r 2 -R 2 , KB(r, R) = 2r √ r 2 -R 2
This equation has one unknown, the radial velocity dispersion, σ 2 rr . That is, assuming perfect knowledge of the LOS velocity dispersion profile, σ 2 los (R), if we could solve this equation for the unknown σ 2 rr we would obtain for each assumption of a mass model, {ρ , Φ(r)}, a radial velocity dispersion, σ 2 rr . The question arises: is the solution with respect to σ 2 rr unique?

Theorem 1. The solution of Eq. 3 with respect to σ 2 rr for given σ 2 los (R), ρ (r) and Φ(r) profiles, is unique. We provide the proof of Theorem 1 in Appendix A. This result, which complements the published inversion techniques, has the following implication: we can assume only one of the three unknown functions and thus reduce the uncertainty of the modelling parameters (in comparison with the uncertainty we get by assuming parametric forms for two unknown functions, as is customary). For the case of statistical fitting, we can use hierarchical models (e.g. smoothing splines) of varying complexity, that make model selection possible and this is the key for breaking statistically the Jeans degeneracy: the missing ingredient (an additional equation) is replaced by the model selection criterion. Thus, if we allow the total mass, M (r) to vary, i.e. if we assume a different mass model for the same σ 2 los profile, then the anisotropy profile will generally be different. However there is one important constraint we need to consider, namely, the projected virial theorem (discussed in detail in Section 4.5.1). The projected virial theorem does not depend on the anisotropy, which implies that not all mass profiles are consistent with the projected kinetic energy evaluated from σ 2 los . However, the projected virial theorem, on its own, is not sufficient to break the degeneracy (it is a single scalar equation, therefore the total number of equations is still less than the unknowns). It can only further reduce the feasible solution space of where the M (r) function resides. We will discuss this further in Section 6. In Section 5.1 we provide numerical examples of the uniqueness of the kinematic profile for an assumed mass density.

It should be stated that we can choose equally well to assume a functional form for the anisotropy profile, and leave the mass density to be deduced by the data (Mamon & Boué 2010, see also Read & Steger 2017): in this case, the total mass of the system follows from the assumptions of the anisotropy model, β, for a given data set of observables. This can be very easily seen from the following: again, assuming perfect knowledge of σ 2 los profile, once we use a specific functional form for the anisotropy β, the system of equations that describes a stellar dynamical system is:

1 Σ r vir R ρ K1(r, R)σ 2 rr + ρ K2(r, R)σ 2 tt d r = σ 2 los (R) (4) σ 2 tt -2(1 -β(r))σ 2 rr = 0 (5) 1 ρ d dr ρ σ 2 rr + 2σ 2 rr -σ 2 tt r = - GMtot(r) r 2 (6)
These are three equations with respect to the three unknowns σ 2 rr , σ 2 tt and the total mass Mtot. The system of equations, Eq. 4, 5 is complete, i.e. we have a unique solution for σ 2 rr and σ 2 tt . Then, from the last Eq. 6 we calculate the total mass, Mtot, whose value depends solely on the tracer stellar density, ρ , and the kinematic profile, σ 2 rr , σ 2 tt . Therefore, when we model by assuming a specific anisotropy profile, β, we effectively pre-specify the mass content of the system. In this reasoning, we did not need to adopt any assumptions for the parametric form of the DM mass density.

DATA

In this section we provide an overview of the datasets we used to validate our methodology. We describe how we preprocess the data and create validation and test data sets for the t-JEAnS solver, as well as how we use GANs to generate large artificial samples of data for the case of small datasets.

We test our algorithm with the Gaia Challenge 2 suite of mock simulations, in particular the spherically symmetric sets. The mock suites provide a snapshot of the full 6D information of the tracer profile, x, y, z, vx, vy, vz. For modelling each of the systems, we used only the projected positions, x, y and the LOS velocity, vz. Following the Gaia Challenge guidelines, we used the data that include velocity errors. For each datum we considered that this error is equal to the 2% of the true vz velocity. Our training data set, Dtrain, consists of the values, {x, y, vz}j, j = 1, . . . , N , as well as the second order moments, σ i los , of the projected LOS velocity.

As a sample of the various datasets, we chose the PlumCuspOM, PlumCuspIso, PlumCuspTan and NonPlum-CoreOM suites. For the first three we use the suites with 10k targets, and for the last one we use both 10k and 1k datasets. The Plummer-like family of tracer profiles was chosen based on the knowledge that most Stellar profiles observed in nature are cored; the latter three models were considered to be representative by the curator of the Gaia Challenge (Read & Steger 2017). In particular, the NonPlumCoreOM is a notoriously difficult set to model, and this is the reason why for this particular one we also include a set with only 1k targets. Each of these datasets was modelled with assumptions for the stellar and DM profiles only. We model each system with two competing models: one with the true parametric form (with parameters recovered from the fitting process), and one with an incorrect parametric assumption (again the parameters are fitted to the data). We report the combinations of Stellar and DM models we used in Tables 2 and6. In all model fits, the anisotropy profile is evaluated from the data.

The reference anisotropy profiles that these data sets were created from are two, namely constant and Ossipkov-Merritt (Osipkov 1979;Merritt 1985):

β(r) =    β0, constant r 2 r 2 + r α a , OM (7)
The mass density profile that these data sets follow, both for the stellar and the DM components, is given by a power law form Zhao (1996), for a variety of reference parameters:

ρ(r) = ρ0 r rs -γ 1 + r rs α (γ-β)/α (8)
In addition, the Gaia Challenge datasets make the approximation that the stellar tracer mass is negligible in comparison to the DM mass component. In order to account for this we normalized the total tracer mass to unity, i.e. M tot = 1M .

We report the reference model parameters for each of the datasets used in Table 1. In Tables 2 and6 we record the combinations of Stellar and DM mass models we used for the modelling process as well as the test error for the various competing models.

Data preprocessing

In this section we describe the process that we followed in order to create the LOS velocity dispersion, σ i los , values, that we use as training data.

Each bin contains Ntargets = √ N . For example, for 10 4 stars, we have N bin = 10 2 . This approach gives equal Poisson error (Read & Steger 2017) for each datum. We modelled the distribution of stars, within each i bin, as a Gaussian centred at zero. The likelihood of this model, for each bin i, is:

Pi(vj|s, D) = N targets j=1 exp{-v 2 j (s 2 + (δvj) 2 ) -1 /2} 2π(s 2 + (δvj) 2 ) . ( 9
)
where vj is the value of the LOS velocity of star j in bin i, δvj is the associated error, and s the standard deviation of the Gaussian distribution. For each bin i, we perform an MCMC process using the likelihood Eq (9), to estimate the marginalized distribution of the parameter s. It should be clear that this MCMC process is used only in the data preprocessing stage. It should not be confused with the MCMC we perform later for the estimation of marginalized distributions for the stellar, θ , and DM, θ•, parameters. The LOS velocity dispersion data values, σ i los , we use at each location Ri (centre of the ith radial bin), is the mode value of the histogram, ŝ2 = σ i los . The associated error, δσ i los , is the 1σ uncertainty of s 2 . Thus, Dtrain = {x, y, vz}j ∪ {σ i los , δσ i los }, j = 1, . . . , N and i = 1, . . . , N bin .

In addition to the above LOS moments, we draw 100 random samples, σ ij los (j = 1, . . . , 100), from the marginalized distribution of σ 2 i , that we keep for a validation data set, D val , and 200 random samples that we use for test sets, Dtest. During the Evolutionary Algorithm (hereafter EA) training, the validation set is used for the selection of the smoothing parameters, θ smooth (Section 4). During the MCMC training phase, instead of using the mode value σ i los of the LOS dispersion as the moments data, in each iteration of the solver, we select random realizations, σ i los = σ ij los (random j), from the validation data set, D val . In this way we incorporate the uncertainty of the moments data as prior information to the modelling process. The test set, Dtest, is used for the model selection between competing models after the EA phase. In Fig. 1 we plot for the case of the Plum-CuspTan model the binned σ i los values (black solid line), as well as the validation values, σ ij los ∈ D val , for each bin i.

Data augmentation for small datasets using GANs

In this section we briefly describe the application of GANs for the numerical reconstruction of the 3D projected LOS ve- 1,2,5,0.1] [0.064,1,1,3,1] β OM ,ra = 0.1 PlumCuspIso (10k) [0.25,2,5,0.1] [0.064,1,1,3,1] β 0 = 0.0 PlumCuspTan (10k) [0.5,2,5,0.1] [0.0239,2,1,4,10k) [0.25,2,5,1] [0.400,1,1,3,0] β OM , ra = 0.25 locity distribution, f (x, y, v los ) from the NonPlumCoreOM 1k dataset. Our goal is to give an intuitive understanding behind the reason that this method is so effective and not to detail the GAN methodology (see Goodfellow 2017 for a pedagogical introduction).

0• in M pc -3 . DataSet θ = [r , α , β , γ ] θ• = [ρ 0• , r•, α•, β•, γ•] β(r) anisotropy PlumCuspOM (10k) [0.
A fundamental limitation to the method of moments, in the Jeans framework, is that it requires a wealth of data to be successful. This is because the moments of the data, as a product of the summary information of the underlying distribution, are much fewer in number than the original unbinned dataset. This is more evident especially when the original dataset is small (from few hundred to 1k stars) as is often the case in astronomical datasets (e.g. of dSph galaxies). We overcome this difficulty by applying a preprocessing step, where we create synthetic data from a generative model, that resembles the true underlying distribution. That is, we create synthetic data to complement the original dataset and thus acquire a large number of LOS velocity moments. We do so only for the 1k NonPlumCoreOM dataset (although the method can be applied to the 10k as well for higher quality results). For this task, artificial intelligence actors (GANs) are excellent generative models, since they learn by "looking" at the real data, i.e. by example, and are not bound by assumptions of the mathematical form of the underlying distribution.

The general framework of the GANs consists of a set of two competing artificial neural networks (hereafter ANNs). The first, the Generator (hereafter G), takes as input a vector of random numbers and tries to create fake (synthetic) data whose distribution resembles the distribution of the true training dataset. The second, the Discriminator (hereafter D), takes as input, true data, drawn randomly from the training distribution, or fake data, created randomly from G, and tries to predict whether the data that it was given are genuine (real) or fake. During training, the goal of G is to make D perform a mistake, i.e. the goal of G is to generate as authentic looking synthetic data as possible. The goal of D is to discriminate the true data from the fake ones and debunk the efforts of G. This framework is a minimax two-player game. During training both players become proficient in their task. When this process reaches equilibrium, G is a faithful approximator of the true underlying distribution of the training dataset. This method is unsupervised training which in practice means there is no upper bound on the quality of the data approximation.

This method has been applied successfully, with impressive results, in artificial intelligence generative tasks, such as the creation of high quality images (Karras et al. 2017), for the creation of synthetic MRI scans for enhanced deep neural network training (Shin et al. 2018), for motion transfer in videos [START_REF] Chan | [END_REF]) and many more cases where the data distribution is anything but "easy" to express mathematically (if not impossible).

For our particular needs we construct a pytorch [START_REF] Paszke | NIPS-W[END_REF]) implementation of Wasserstein GANs with gradient penalty (hereafter WGAN-GP, Gulrajani et al. 2017). We chose WGAN-GP because it is one of the most reliable GAN frameworks for stability in training. The architectures we used for the G and D ANNs are summarized in Table 3. The input to the generator is a random 10 dimensional multinomial distribution, z ∼ N (0, 1) 10 ∈ 10 . In Table 4 we detail the hyper parameter values we used during GANs training. In addition, in order to avoid overfitting the NonPlumCoreOM 1k data set, we augmented the data with random rotations on the x, y plane and reflections with respect to x and y axis. In particular we followed the transformations (x, y, v los ) → (-x, y, -v los ) and (x, y, v los ) → (x, -y, -v los ). For zero mean v los stellar systems, these reflections are like observing the target from the opposite direction of the initial observer: clearly the physics of the system should not change. This type of information should be viewed as "prior knowledge encoding" of the modelling process with neural networks.

In Fig 3 we plot on the (R, v los ) plane the synthetic data generated from the GANs against the 1k and 10k Non-PlumCoreOM datasets. We generated ∼25k synthetic data points by training the Discriminator, D, on the NonPlum-CoreOM 1k dataset. This resulted in approximately 160 σ 2 los binned values for the LOS velocity dispersion profile. In Fig. 4 we plot the LOS velocity dispersion profile from the GAN data as well as the true 1k and 10k dispersion profiles. In all panels the reference profile (dashed curve) is overplotted. Clearly, the GAN generated profile is of high quality. In fact, the uncertainty of the data points around the reference profile is smaller than even the case of the original 10k dataset. This happens because the GAN system learns more information of the underlying distribution from the NonPlumCoreOM 1k dataset than what the moments of the 10k sample can describe. As a result, with higher number of targets (25k) we end up with a LOS velocity dispersion profile of smaller uncertainty than the 10k original dataset. A small bias is apparent in the last two σ 2 los data points, probably because the GANs overfit the outliers at the edges of the radial distance of the 1k dataset. This bias may also be due to the system of GANs not having reached the optimum equilibrium when we terminated training. Finally in Fig 5 we compare the projected density (brightness for ΥV = 1) of the tracer population. It should be noted that we did not experiment with new architectures, training schemes or hyperparameter optimization. We just used the proposed implementation scheme from Gulrajani et al. (2017) for their toy model of 25 2D Gaussian distributions. There is huge scope for improvement and adaptation for individual datasets of this technique for data augmentation in astronomy in various sub-disciplines. Here, we are merely scratching the surface of the potential of this technology.

THE t-JEAnS SOLVER

In this section we present an overview of the t-JEAnS algorithm we developed for accurate mass estimates in spherically symmetric self-gravitating systems.

The JEAnS (Diakogiannis et al. 2017) algorithm is a numerical solver that estimates the mass content and the kinematic profile of spherically symmetric gravitating systems. It models independent of anisotropy, β(r), assumptions and it requires parametric functional forms for the mass density profiles. The best mass model is selected with the use of model selection criteria (Diakogiannis et al. 2017) Table 3. Generator and Discriminator network architectures. We follow pytorch semantics to denote the dimensionality and type of the layers and non-linear activations we used. Here LDIM=10 is the dimensionality of the latent space that we sample and feed into the Generator, DIM=512 is the number of features in the linear Layers and XDIM=dim(x, y, v los ) = 3 is the dimensionality of the projected observation space.

Layer

Generator Discriminator (Akaike Information Criterion, [START_REF] Sugiura | [END_REF]Burnham & Anderson 2002, hereafter AICc). The radial velocity dispersion profile, σ 2 rr , is represented as a "free form" B-spline function, σ 2 rr (r) = i a i Bi(r). The correct kinematic profile is inferred from the data. The solver uses information of brightness and line-of-sight velocity moments, σ 2 los , to estimate marginalized distributions of the mass model parameters as well as the coefficients, a i , that describe the radial dispersion profile, σ 2 rr . The algorithm consists of three distinct phases. In the first phase it evaluates the simplest kinematic profile that gives a satisfactory representation1 to the data, as well as the most probable mass model. This is achieved with the use of evolutionary optimization and quadratic programming. In the second phase, JEAnS evaluates the optimum smoothing parameters from ideal theoretical models. Finally, in phase three the algorithm performs MCMC inference, for the determination of marginalized distributions of the model parameters.

The new version t-JEAnS is significantly modified compared to the previously-published version (Diakogiannis et al. 2017). In the first phase we again evaluate the optimum B-spline basis, as well as the statistically most favoured mass model. We introduce a new quadratic programming formalism -the Dynamic Moments Solver (hereafter DMS) -for the numerical solution of the system of coupled integrodifferential Equations (1,2). In the latest version of the JEAnS we expand both the radial, σ 2 rr , and tangential, σ 2 tt , profiles in a B-spline basis of order k = 4 (degree = 3) 2 , i.e. σ 2 rr (r) = c i r Bi(r), σ 2 tt (r) = c i t Bi(r). This allows us to treat the Jeans equation as a local, ri, constraint in the quadratic optimization problem of estimating the velocity moments, σ 2 rr , σ 2 tt . In combination with the local support of B-spline functions, this translates to more equations for the unknown coefficients c i r , c i t that further reduce the feasible solution space. In comparison with the old version of the JEAnS, by solving the Jeans equation (Eq. 16) with respect to σ 2 tt and substituting under the integral sign of the σ 2 los definition (Eq. 14), we loose the local equations that c i r and c i t coefficients participate after the last datum. By keeping the Jeans equation as a constraint we can evaluate equations for c i r and c i t in all space r ∈ [0, rvir]. This has a direct positive impact on the quality of the recovered anisotropy profile, β(r).

In a similar fashion to the first version of the JEAnS, we do not invert the dynamical equations, thus we avoid the problem of having to integrate/differentiate noisy numerical functions. We also include additional global and local constraints that guarantee that the kinematic profiles lead to physically acceptable solutions (σ2 rr , σ 2 tt 0, ∀r ∈ [0, rvir)). The fitness function is modified in order to include information from the full line-of-sight kinematics. The optimum smoothing parameters are now evaluated directly from the data according to the best bias-variance tradeoff using a validation data set, D val . The model selection is performed using a hold out test data set, Dtest. In phase two we per-form MCMC inference for the unknown stellar and DM mass model parameters, θ , θ•. In this phase, the kinematic profile is treated as a nuisance parameter. Finally in the third phase, we perform stochastic programming (SP) in order to determine confidence intervals for the velocity dispersion profiles, σ 2 rr , σ 2 tt , σ 2 los . In more detail (Fig. 2), the distinct phases of the t-JEAnS are the following:

(i) An evolutionary optimization (EA) phase. In this phase we determine: a) the simplest (best) B-spline basis1 for the representation of the unknown radial, σ 2 rr , and tangential, σ 2 tt velocity dispersions, (b) the best candidate mass models and, (c), the best smoothing 2 parameters, θ smooth = {λ1, β1, λ2, β2}. For the evaluation of the smoothing parameters we use a validation data set, D val , created from random sampling from the LOS σ i los marginalized distributions (see Fig. 1). The optimum smoothing parameters are the ones that minimize the validation error for all random samples, σ i los . We give more details of this process in the section where we describe the fitness function. For the model selection, we use a "hold-out" LOS moments test data set, Dtest (Section 3.1), and we perform model selection (Section 4.8) based on the out-of-sample prediction error (generalization test error ). This approach gives more robust model selection (in comparison with predictive information criteria (Gelman et al. 2014)), since it heavily penalizes models that do not generalize well on unseen data. It should be stressed however, that the efficiency of the model selection process depends crucially on the number of available data points. (ii) A Markov Chain Monte Carlo (MCMC) analysis, keeping the B-spline basis and the smoothing parameters fixed, for the best mass model. In this scheme, the radial and tangential coefficients, c i r , c i t are treated as nuisance parameters: they are estimated at each iteration from the DMS. This phase produces marginalized distributions of the parameters of stellar, θ , and DM, θ•, mass densities, {ρ , ρ•}. (iii) A stochastic programming (SP) phase, where the θ , θ• parameters are used iteratively in the DMS. This produces marginalized distributions for the radial and tangential coefficients, c i r , c i t , subject to local and global dynamical constraints. This last phase gives the required uncertainty of LOS and radial and tangential velocity dispersions.

Mass models

For our modelling purposes we used the following candidate mass models:

ρ(r) =                      ρ0 [1 + (r/rs) 2 ] 5/2 Plummer ρ0 (1 + r/rs)[1 + (r/rs) 2 ] Burkert r 3 s ρ0 r(r 2 + r 2 s ) 2 NFW Eq (8)
generalized Hernquist (10)

We model each dataset with two different mass model assumptions, the correct one and an incorrect one. Our goal is to demonstrate that given sufficient data it is possible, in principle, to statistically infer the most probable model using model selection criteria.

Dynamic Moments Solver (DMS)

In this section we describe the mathematical representation of the problem, i.e. the dynamic equations that enable us to recover the radial and tangential velocity moments, from knowledge of the LOS velocity dispersion, σ 2 los , the tracer, ρ , and the DM, ρ•, mass densities. The DMS solves the system of coupled integrodifferential equations (Eq. 1,2) by discretizing the solution space using B-splines. This is achieved by expanding the unknown radial, σ 2 rr , and tangential, σ 2 tt , velocity moments in a B-spline basis 3

σ 2 rr (r) = c i r Bi(r) (11)
σ 2 tt (r) = c i t Bi(r) (12)
The DMS takes as input the knots, ξi, the stellar parameters, θ , the DM parameters, θ• and the smoothing penalty variables, θ smooth = {λ1, β1, λ2, β2} and gives as output the coefficients c i r , c i t that fully describe the radial and tangential velocity moments. Using the approximation Eqs. ( 11) and ( 12), the task is transformed to a convex optimization problem (quadratic programming). The software library we use in t-JEAnS for the quadratic optimization is IBM's CPLEX 4 .

For clarity in notation, it is convenient to represent the DMS as a function:

DMS(θ|Dtrain) → (c i r , c i t ),
where θ ≡ {ξ, θ , θ•, θ smooth } are the parameters that define the B-spline basis, the tracer and DM profiles, as well as the smoothing penalty regularization. The goal of the DMS is to minimize the training error of the LOS velocity dispersion:

χ 2 train = N bins i σ 2 los (Ri) -σ i los δσ i los 2 (13)
subject to various local and global dynamic equations (constraints). We separate these constraints into local, boundary and global constraints. In addition we will impose some regularization conditions (smoothing) in the minimization process, in order to reduce the condition number of the linear system and avoid oscillatory solutions. We formally define the objective function of the DMS in Section 4.6. In Table 5 we summarize the system of equations and the objective function that fully describe the DMS. We proceed by stating exactly the mathematical equations we use in the t-JEAnS.

The LOS velocity dispersion under the B-spline approximation of the velocity moments is given by:

σ 2 los = 1 Σ (R)
K1(r, R) = 2(r 2 -R 2 ) r √ r 2 -R 2 K2(r, R) = R 2 r √ r 2 -R 2 K3(r, R) = 2r √ r 2 -R 2
are kernel functions. Applying the B-spline approximation (Eq. 11 and 12) and defining:

I r i (R) = 1 Σ r vir R ρ K1(r, R)Bi(r) d r I t i (R) = 1 Σ r vir R ρ K2(r, R)Bi(r) d r
the linearized LOS velocity dispersion takes the form:

σ 2 los (R) = c i r I r i (R) + c i t I t i (R) (15)
This is the model function that we compare with observables, subject to physical constraints. It is linear with respect to the unknown coefficients, c i r , c i t , something that simplifies the solution and allows for convex optimization.

Local constraints

These constraints are termed local, because they are valid in the whole extent of the system, r ∈ [0, rvir]. We evaluate these at the positions of the Greville abscissae of the B-spline basis.

Jeans constraints

The spherically symmetric Jeans equation (SSJE) is:

-ρ d Φ d r = d(ρ σ 2 rr ) d r + ρ (2σ 2 rr -σ 2 tt ) r (16)
The linearized form of SSJE that results from the B-spline approximation is:

-ρ d Φ d r = d(ρ Bi) d r + 2ρ Bi r c i r + ρ Bi r c i t (17)

Sign constraints

We demand the velocity moments to be positive in all solution space:

σ 2 rr (r) 0 σ 2 tt (r) 0.
In terms of the kinematic coefficients:

c j r Bj(r) 0 (18) c j t Bj(r) 0 (19)

Boundary constraints

These apply at the origin and at the virial radius of the system.

σ 2 rr (0) = σ 2 tt (0)/2 σ 2 rr (rvir) = σ 2 tt (rvir) = 0 The reasoning for the σ 2 rr (0) = σ 2 tt (0)/2 boundary condition is the following: we expect that all tangential motions at the limit r → 0 become radial. That is, if we draw the tangent line to a circle of radius r, as the radius approaches zero, the tangent line approaches the origin r = 0 of the coordinate system. In the limiting case where r → 0 the tangent line passes from the origin (it is actually a degenerate case: all directions are equivalent). In this respect it is our understanding that in this limit the tangential and radial motions are indistinguishable. This is why we expect that their dispersions will be equal at r → 0. With regards to the second boundary constraint, it is proven (Dejonghe & Merritt 1992) that for a self consistent system in virial equilibrium the radial and tangential velocity dispersions vanish in the limit of the virial radius.

Global constraints

In this category fall constraints of local functions integrated over all space.

Projected virial theorem

The virial theorem states [START_REF] Binney | Galactic Dynamics: Second Edition[END_REF][START_REF] Merritt | Dynamics and Evolution of Galactic Nuclei[END_REF]) that if K is the total kinetic energy of a system, and W its total potential energy, then for a system in dynamic equilibrium:

2K + W = 0 (20)
For a sperically symmetric system,

W = 4π r vir 0 ρ - d Φ d r r 3 d r
The total kinetic energy of a system, defined via the line-ofsight velocity dispersion is:

K los = 3 2 r vir R=0 d R 2πRΣ (R) σ 2 los (R)
Substituting σ 2 los (R) from Eq. ( 15), we have

K los = c i r K Rlos i + c i t K Tlos i ( 21
)
where

K Rlos i = 3π r vir R=0 RΣ (R)I r i (R) dR K Tlos i = 3π r vir R=0 RΣ (R)I t i (R) dR
Substituting in Eq. ( 20) yields:

2(K Rlos i c i r + K Tlos i c i t ) + W = 0. ( 22
)
This is an additional constraint on the c j r , c j t coefficients. From the perspective of linear/quadratic programming algorithmic structure, Equation ( 22) is a hyperplane equation with respect to the unknown coefficients, c j t , c j r , that further reduces feasible solution space.

The projected virial theorem is also a hard bound on the value of the total gravitational energy of the stellar and dark matter ( •) interaction. Furthermore, it is clear that since the K los value is independent of the anisotropy profile (i.e. it is an observational fact), then it is impossible to have only the stellar component with some peculiar anisotropy profile to represent the observables. In other words, the total gravitational energy of the system is fixed from the total kinetic energy as this is estimated from the LOS dispersion, σ 2 los . That is, the constraint of virial equilibrium does not allow one to vary the anisotropy profile β to fit any desired mass profile.

Objective Function

The DMS objective function that relates observables, σ 2 los , with the model function (Eq. 15) is given by:

F = i σ 2 los (Ri) -σ i los δσ i los 2 + λ1[β1 n coeff -1 i=1 (∆ 1 c i r ) 2 + (1 -β1) n coeff -2 i=1 (∆ 2 c i r ) 2 ] + λ2[β2 n coeff -1 i=1 (∆ 1 c i t ) 2 + (1 -β2) n coeff -2 i=1 (∆ 2 c i t ) 2 ] (23)
where the difference operators ∆ 1,2 are defined by:

∆ 1 c i r,t = c i+1 r,t -c i r,t ∆ 2 c i r,t = ∆ 1 (∆ 1 c i r,t ) = c i+2 r,t -2c i+1 r,t c i r,t + c i r,t .
The coefficients, λ1,2 regulate the ammount of smoothing penalty on each of the velocity dispersions. The coefficients β1,2 regulate the relative contribution of the first and second derivative penalties for each velocity dispersion. This smoothing penalty is efficient and very fast to evaluate in comparison with previous efforts (Diakogiannis et al. 2014b(Diakogiannis et al. , 2017)). It is the same penalty used in the P-splines (Eilers et al. 1996) formulation in statistical smoothing.

Fitness function

The EA phase of the t-JEAnS solver evaluates the simplest B-spline basis that best represents the observables. This is a nested optimization: the EA parameters consist of the stellar, θ , the DM, θ• and the smoothing penalty variables, θ smooth = {λ1, β1, λ2, β2}. Once these parameters, θ = {θ , θ•, θ smooth }, are proposed, then the problem is a quadratic programming optimization problem, with respect to the c i r , c i t unknown constants. The optimal variables, ĉr i , ĉt i for the proposed θ parameters are evaluated with the DMS. The evaluation of the model though, takes into account information from both the DMS and the full kinematics. For the full kinematics, we use definitions (Mamon et al. 2013) based on assumptions of a Gaussian distribution function for the velocities (in 3D space), truncated at the escape velocity of the system.

The fitness function is defined with the usage of model selection criteria (BIC, AICc) and the following penalty functions:

AICc = 2 i 1 2 σ2 los (Ri) -σ i los δσ i los 2 + 2n + 2n(n + 1) N data -n -1 χ smooth = 1 N sample N sample j=1 N bins i=1 σ2 los (Ri) -σ ij los δσ i los 2 χ virial = 2K los W - W 2K los BIC = -2 N batch i=1 log q(Ri, v i los )) + n log(N batch )
where

q(R, v los ) = 2πR M tot g(R, v los ) (24) g(R, v los ) = r vir R rρ (r) √ r 2 -R 2 h(v los |R, r)dr (25) h(v los |R, r) = exp[- v 2 los 2σ 2 z (R,r) ] 2πσ 2 z (R, r) erf{vesc(R)/ 2σ 2 z (R, r)} (26) σ 2 z (R, r) = σ 2 rr (1 -(R/r) 2 ) + σ 2 tt (R/r) 2 /2 (27)
are the full kinematics definitions from the MAMPOSSt (Mamon et al. 2013) algorithm. The projected virial theorem is satisfied by the quadratic programming solver, within some numerical tolerance. We found that we got slightly faster convergence by also penalizing this explicitly in the EA solver, with the χ virial term. The values σ2 los (Ri) are the solutions from the Dynamical Modelling Solver (DMS) for the given input parameters θ = {ξ, θ , θ•, λ1,2, β1,2}. The data values, σ i los and σ ij los are produced from a binning scheme as described in Section 3.1. We remind the reader that the values σ ij los ∈ D val are j sampled values for each bin i. They are used as a validation set for determining the smoothing parameters, θ smooth = {λ1, β1, λ2, β2}. N batch is the size of a random sample (without replacement) of full kinematics data of stars from the population. We use N batch = 1000: this provides a good approximation to the full kinematics likelihood and allows for faster convergence.

The fitness function, treated as a maximization problem, is the product of four components, namely:

fDMS = 1 1 + AICc f full kin = 1 1 + BIC f smooth = 1 1 + χ smooth fvir = 1 1 + χvir .
Then, the fitness function is:

ftot(θ) = fDMS f full kin f smooth fvir (28)

Model selection

Model selection takes place in two distinct processes inside the t-JEAnS. Once we select a set of tracer and DM mass densities, we use the EA in order to find the simplest B-Spline basis for the radial and tangential velocity dispersions. This task is a hierarchical model selection problem

σ 2 tt = c i r B i , σ 2 rr = c i t B i .

Mathematical formula

Objective function min

F = i (σ 2 los (R i ) -σ i los ) 2 + λ 1 [ n coeff -1 i=1 β 1 (∆ 1 c i r ) 2 + (1 -β 1 )(∆ 2 c i r ) 2 ] +λ 2 [ n coeff -2 i=1 β 2 (∆ 1 c i t ) 2 + (1 -β 2 )(∆ 2 c i t ) 2 ] Model function σ 2 los (R) = c i r I r i (R) + c i t I t i (R) Local constraints Jeans: -ρ d Φ d r = d(ρ B i ) dx + 2ρ B i x c i r + ρ B i x c i t sign: c j r B j (x) 0, c j t B j (x) 0 boundary: σ 2 rr (0) = σ 2 tt (0)/2 σ 2 rr (r vir ) = σ 2 tt (r vir ) = 0 Global constraints projected virial: 2(K Rlos i c i r + K Tlos i c i t ) + W = 0
(where the various competing models are the ones that have different number and locations of knots, but the same mass density parametric form). For this task, AICc or BIC, based on the training error measure (likelihood) prove to be good choices.

However, when one needs to compare competing mass models that were trained in distinct EA phases, it is best to use out-of-sample data and test how well the model generalizes on unseen (during training) data (Section 3). Once the EA phase is complete, for competing mass models, we evaluate the best model using the hold-out LOS moments test set in a cross-validation manner. The average error on unseen moments test data that we use is:

χ 2 test = 1 2Ntest N test j=1 N bin i=1 σ2 los (Ri) -σij los δσ i los 2 (29)
where σij los are out of sample test data, created for each bin i by random sampling from the marginalized distribution of the data preprocessing MCMC chains. The model with the smallest test error is selected as the best candidate. In our experiments this method has proven to be more robust than predictive training error methods (e.g. AICc), which can have bias from overfitting (Gelman et al. 2014).

Likelihood function

In Phase II of the t-JEAnS we perform an MCMC exploration using the following likelihood (Mamon et al. 2013;Ibata et al. 2013;Diakogiannis et al. 2017Diakogiannis et al. , 2014a,b),b):

L =    N bin j=1 exp{- (Σ (R j )-Σ j ) 2 2(δΣ j ) 2 } 2π(δΣ j ) 2    ×     N bin i=1 exp - (σ 2 los (R i )-σ i los ) 2 2δ(σ i los ) 2 2πδ(σ i los ) 2     × (30) × λ1e -λ 1 W 1 λ2e -λ 2 W 2 N batch j=1 q(Rj, v j los ) (31)
where Σ (Rj) is the projected tracer density at location Rj, Σ j and δΣ j the observed projected mass density and its uncertainty, σi los is a random sampled value (at each iteration, we use values from D val ) from the ith MCMC binned histograms and W1, W2 are given by:

W1 = β1 n coeff -1 i=1 (∆ 1 c i r ) 2 + (1 -β1) n coeff -2 i=1 (∆ 2 c i r ) 2 (32) W2 = n coeff -1 i=1 β1(∆ 1 c i t ) 2 + (1 -β1) n coeff -2 i=1 (∆ 2 c i t ) 2 (33)
The B-spline knots, and the coefficients λ1,2, β1,2 are kept fixed to the values of the best EA solution. The full kinematics likelihood is calculated on each iteration on a random sample (without replacement) of N batch = 1000 stars. This is sufficient for the algorithm to converge in an excellent trade-off between computational efficiency and parameter constraints. We use random samples, σi los , as data in each MCMC iteration in order to avoid overoptimistic constraints for the marginalized distributions of parameters. In this way we incorporate the uncertainty of the binned LOS dispersion 5). In order to obtain this solution we assumed perfect knowledge of the σ 2 los , ρ and ρ• profiles.

values in the marginalized distributions of the θ , θ• parameters.

Stochastic Programming

Once we have clean (after burn-in phase) MCMC chains of the θ , θ• parameters, we estimate distributions of cr, ct by applying the DMS solver iteratively to each pair of MCMC values θ j , θ j • . Here, the index j indicates the jth MCMC chain. For this computation we keep the smoothing penalty parameters, as well as the B-spline basis, fixed to the best EA values. In functional form:

PDF(c i r , c i t ) = DMS θ j , θ j • ∼ PDF(θ , θ•)|D (34)
where the symbol θ j , θ j • ∼ PDF(θ , θ•) denotes that θ j , θ j

• are sampled at random from their marginalized distribution PDF(θ , θ•) (estimated from the MCMC chains). Finally, from the marginalized distributions of cr, ct, θ , θ•, we can estimate 1σ uncertainty intervals for the velocity moments and the various mass model functions.

RESULTS

In this section we summarize our findings for both the exact solution of the Jeans system of equations and the statistical fitting of the Gaia Challenge dataset.

Exact solutions

In Figures 6 and7 we plot the exact solutions of the system of the Jeans equations (Table 5) using the DMS solver, for the case of the PlumCuspOM and PlumCuspIso reference profiles. Our aim here is to provide numerical "proof of concept" examples of Theorem 1. That is, by assuming full knowledge of the LOS velocity dispersion profile, the tracer ρ and ρ• mass densities, we recover a unique kinematic profile as this is described by the second order radial, σ 2 rr , and tangential, σ 2 tt , velocity moments. In this approach we are not using smoothing penalty coefficients (λ1,2 = 0 in the objective function Eq 23) since we have a wealth of data points. For the exact solution we use a large number of B-Spline basis, dim{Bi(r)} ∼ 150. For each of the two figures, from top to bottom panels: data σ 2 los and recovered solution, reference and recovered tangential velocity dispersion (σ 2 tt ) and reference and recovered radial velocity dispersion profile (σ 2 rr ).

10 8 M • (r) (M ) 2 4 M fit • /M true • 10 -2 10 -1 10 0 10 1 R (kpc) 0 100 σ 2 los (km 2 s -2 ) fit true 1σ σ i los 0 250 500 σ 2 rr (km 2 s -2 ) σ i rr 10 -2 10 -1 10 0 10 1 r (kpc) 0 100 σ 2 tt (km 2 s -2 ) σ i tt 10 -2 10 -1 10 0 10 1 r (kpc) -1 0 1 β = β/(2 -β)
Figure 8. Fit and true reference profiles for the GAIA Challenge dataset PlumCuspOM, for 10k targets. In the left panels we plot (top) the DM mass, (middle) the ratio of the estimated mass over the true mass and (bottom) the anisotropy profile. In the right panels we plot (top) the fit to the LOS observables, (middle and bottom) the recovered radial and tangential profiles. Overplotted are the true σ i rr , σ i tt dispersions, as estimated from the data; these were not used in the fitting process. The blue region corresponds to 1σ uncertainty for all of the quantities.

10 4 10 8 M • (r) (M ) 1 2 3 M fit • /M true • 10 -2 10 -1 10 0 10 1 R (kpc) 0 200 σ 2 los (km 2 s -2 ) fit true 1σ σ i los 0 200 400 σ 2 rr (km 2 s -2 ) σ i rr 10 -2 10 -1 10 0 10 1 r (kpc) 0 100 200 σ 2 tt (km 2 s -2 ) σ i tt 10 -2 10 -1 10 0 10 1 r (kpc) -1 0 1 β = β/(2 -β)

Statistical fitting

10k datasets

Our results are summarized in Figures 8, 9, 10 and11 andTable 2. We fully recover the mass content and the anisotropy profiles in a representative sample of synthetic data sets from the Gaia Challenge1 suite of mock simulations. In Fig. 8 we plot the best fitting model, as well as [ht] Figures 9, 10 and 11 are as Fig. 8 for the Gaia Challenge datasets with: a Plummer like tracer profile with cuspy dark matter halo and isotropic velocity anisotropy (PlumCuspIso), a Plummer like tracer profile with cuspy DM halo and tangential velocity anisotropy (PlumCuspTan) and a cuspy like (non-Plummer) tracer profile with a cored halo and Ossipkov-Merritt (Osipkov 1979;Merritt 1985) velocity anisotropy profile (NonPlumCoreOM) datasets. In all four cases, our algorithm selects the correct model and reconstructs robustly the mass content and kinematic profile of the underlying stellar distributions from LOS data only.

10 4 10 8 M • (r) (M ) 1.0 1.5 M fit • /M true • 10 -2 10 -1 10 0 10 1 R (kpc) 0 200 σ 2 los (km 2 s -2 ) fit true 1σ σ i los 0 100 200 σ 2 rr (km 2 s -2 ) σ i rr 10 -2 10 -1 10 0 10 1 r (kpc) 0 100 200 σ 2 tt (km 2 s -2 ) σ i tt 10 -2 10 -1 10 0 10 1 r (kpc) -1 0 1 β = β/(2 -β)
10 3 10 8 M • (r) (M ) 1 2 M fit • /M true • 10 -2 10 -1 10 0 10 1 R (kpc) 0 250 σ 2 los (km 2 s -2 ) fit true 1σ σ i los 0 500 1000 σ 2 rr (km 2 s -2 ) σ i rr 10 -2 10 -1 10 0 10 1 r (kpc) 0 200 σ 2 tt (km 2 s -2 ) σ i tt 10 -2 10 -1 10 0 10 1 r (kpc) -1 0 1 β = β/(2 -β)
In Table 2 we report the results of the mass model selection during Phase I of the t-JEAnS. The model selection is performed using the average test error on unseen (during training) data (Eq 29) as it has proven to be a more robust discriminator (in comparison with AICc or BIC). We perform model selection after Phase I, in order to reduce computation time. In general, better discrimination results between competing models can be achieved by performing the MCMC process (Phase II) for both competing models and then evaluating the test error, χ 2 test (Eq 29). In Table 2 we report the average error, χ 2 test , on unseen test data, Dtrain (section 3). In all four cases, the t-JEAnS finds the true underlying models from which the synthetic data were created. In Figures 8 -11 we plot the best candidate models as these were selected from the t-JEAnS. The plotted results were obtained after Phase III of the t-JEAnS. In all cases, our algorithm achieves excellent performance and reconstructs the true underlying profiles.

The 1k NonPlumCoreOM dataset

In this section we discuss our findings for the 1k NonPlum-CoreOM dataset as well as the efficiency of the GANs for synthetic data generation. The latter is judged by the quality of the fits.

The NonPlumCoreOM 1k dataset, besides being a very difficult dataset due to its strong radial anisotropy profile (Read & Steger 2017), also presents a challenge for all Jeans moments based solvers due to its small number of data. Binning 1k data, we end up with as few as 30 binned LOS velocity dispersion values. For a small model, with only 3 knots for the definition of the B-spline basis, we end up with 5 (σ 2 rr ) +5 (σ 2 tt )+ 4(smoothing penalty)+2 (DM) +2 (Stellar) = 18 unknown parameters. In addition, the uncertainty of the σ 2 los binned values is much larger, as is evident from Fig 4.

For these reasons, we fit both the true 1k profile, as well as the augmented GAN profile. For the case of the 1k dataset, the test error based on the sampled mcmc values of the σ 2 los bins fails to recover the correct model. The augmented GAN dataset (approximately 160 binned values) selects the correct model, thus underlining the importance of this data augmentation approach. We summarize the results of the model selection, during the EA phase, in Table 6.

In Fig. 12 we present the fit to the GAN generated data. In the top right panel the σ 2 los data values are the ones created from the 25k GAN generated synthetic data. In the middle and bottom right panels, the σ 2 rr and σ 2 tt data values were not used in the fit. They were estimated from the true NonPlumCoreOM 10k dataset and they are placed there for reference only. We used these because the GAN data do not have the full 3D information for us to create these data val- 11 and 12). This can be quantified, as can be seen in Fig. 13: in the left panel we plot the true σ i rr profile as this is estimated from the 10k NonPlumCoreOM dataset, as well as the highest likelihood fitted profiles, for the GAN data and the 10k NonPlumCoreOM datasets. In the right panel we do the same for the tangential dispersion, σ 2 tt . In order to quantify the quality of the fits in the unseen latent space of radial and tangential dispersions we estimate the mean square error for the radial and tangential profiles, between the best fitted profiles and the data:

χ 2 rr,10k = 1 N bin N bin i=1 [σ 2 rr10k (ri) -σ i r ]/δ(σ i r ) 2 χ 2 rr,GAN = 1 N bin N bin i=1 [σ 2 rrGAN (ri) -σ i r ]/δ(σ i r 2
and similarly for the tangential profile, σ 2 tt . We find for the ratios:

χ 2 rr,10k /χ 2 rr,GAN = 0.99, χ 2 tt,10k /χ 2 tt,GAN = 54.26
Therefore, the quality of the σ 2 rr fit is similar if we train t-JEAnS with either the 10k dataset, or the GAN generated synthetic data. However, the quality of the σ 2 tt fit is much worse when t-JEAnS is trained with moments from the true 10k dataset. This should not come as a surprise. What this means, is that from the 1k of data, the GAN system manages to recover more information than what is hidden in the moments of a 10k dataset. Then, with ∼160 binned σ 2 los data points, it passed more information to the t-JEAnS solver, than the moments of the 10k dataset can.

The principal criticism that is levelled at the Jeans approach is that one may find solutions to the Jeans equations that require a distribution function that is not positive at all phase-space locations, and is hence unphysical. However, one can always check that the results of our algorithm give a positive DF by testing the solution with a single Schwarzschild model. Since the solutions presented above recover the cor-rect input dynamical models from the Gaia Challenge , this step is not necessary here.

DISCUSSION

We suspect that the astronomy community's definition of the Jeans degeneracy would be: many choices of functional forms for M (< r) and β(r) (or equivalently σ 2 rr ) result in a σ 2 los profile that is arbitrarily close to the data. Therefore it is not possible to derive a unique mass and anisotropy profile. This is indeed the case, the system of equations is not closed (we need additional constraints that we do not have). However, when it comes to statistical model selection the situation is different. We can provide the additional necessary condition that closes the system of equations by selecting the "simplest" solution that describes well the observable data. The point of emphasis above in bold, that the profile should be arbitrarily close to the data, resembles a χ 2 "selection" criterion, which is, however, not a proper model selection method. The key point in t-JEAnS to statistically break the degeneracy is the realization that we can use hierarchical1 models that eventually result in different quantitative fits to the data (i.e. different test error). In other words, different assumptions of functional forms for mass, M (r), and anisotropy, β, are no longer quantitatively equivalent.

A special note needs to be made about the fact that the notion of the mass anisotropy degeneracy, when it comes to statistical fitting, is reinforced by the fact that for the majority of stellar systems, the observables are few in number. This makes model selection even more difficult and sustains the belief that, given the availability of data, it is not always possible to discriminate between competing mass models. This is more evident for moment-based mass estimators that rely on summary statistics of the initial dataset. The modern semi-supervised machine learning techniques that are actively being developed by the community, such as the Generative Adversarial Networks for synthetic data generation, are a remedy to this problem. Some of the main differences of the t-JEAnS that allow more efficient treatment of the degeneracy problem, assuming sufficient available data, in comparison with other approaches are:

(i) We do not assume two unknown parametric functional forms for both the DM mass density profile, ρ•(r), and the anisotropy profile, β(r). This reduces the uncertainty of the parameters and allows for more robust model selection. (ii) With our choice of hierarchical parametric models (Bsplines) for σ 2 rr and σ 2 tt we can better statistically discriminate between competing models. This is achieved because hierarchical models find a trade-off between test and train error, and are thus more resilient to over-fitting. (iii) We incorporate a set of physically plausible constraints (Section 4.2) that further reduces the feasible solution space and pushes to the limit the model selection process. (iv) For the moments solver (DMS, Section 4.2), we are not using a single σ 2 los value for each bin. In contrast, we are using the full MCMC chains to get additional information from the binning scheme. This allows for the estimation of train, validation and test error, as it is used in modern machine learning supervised training techniques.

Although we have not performed a detailed numerical comparison by switching on and off all the constraints we used, we have the following understanding of the effect of each as well as our modelling approach:

(i) The choice of the assumed DM mass model: as with all model selection processes, our effort relies on the assumption that, if we try a large set of competing mass models, then one (or some) of them will not be very far from the truth. Then our best solution should approximate reality at a satisfactory level. Our contribution is demonstrating that with the use of a hierarchical basis, satisfactory model selection is possible. Obviously, if our mass model assumptions are away from the truth, we expect that the kinematic fits will also be away from the true anisotropy profile. In our numerical experiments, even with different mass model assumptions, the kinematic profiles tend to be similar. However we cannot conclude, due to the limited number of mass models and data sets we tried, that this is a general feature.

In addition, we cannot quantify the "anisotropy similarity" in terms of similarity between competing DM mass models. This is something that requires further investigation.

We also note that we have found that the choice of tracer profile affects significantly the derived anisotropy profile. (ii) The boundary condition at the origin σ 2 rr (r = 0) = σ 2 tt (r = 0)/2 can result as a limiting case of the Jeans Eq. ( 1), as r → 0, for non divergent DM potentials. However, it helps numerically inside the solver to keep it separate. The boundary condition at the virial radius was used mainly for the domain of definition of the B-spline basis (it requires a closed finite interval). As we cannot deduce the profile further than the last datum, this constraint contributes in combination with the projected virial theorem. (iii) The MAMMPOST-style LOSVD helps to constrain more robustly the kinematic profile beyond the half light radius. It proved helpful in the case of the difficult NonPlumCoreOM dataset. In the other three datasets, even without it, the recovered fits were excellent. (iv) The projected virial theorem can alter the solution space significantly, for a given mass model assumption. For example we find that if we run an MCMC exploration with and without it the parameter chains for the same model converge at different non-overlapping regions. It should also be noted that this is a very difficult constraint to implement numerically in an MCMC scheme, because it is a hard bound and does not allow efficient mixing of the chains. It is possible that there is a connection between the projected virial theorem constraint, and the approach of the virial shape parameters taken by Read & Steger (2017, see also Richardson & Fairbairn 2014), however we have not verified this. It is also interesting to note that despite the fact that the kinematic profile is essentially "free" after the last datum, the virial theorem still helps reducing the feasible solution space.

A special note needs to be made on the particular choice of representation: in t-JEAnS we represent the kinematic profile with the variables σ 2 rr and σ 2 tt instead of σ 2 rr and β. This is because in the former representation, with the use of B-splines, we can linearize the system of equations (thereby greatly simplifying the solution). In contrast if we use σ 2 rr and β, then from Equations ( 1) and ( 2) it is apparent that due to the product term, σ 2 rr β, the system of equations is not linear. It should be emphasized that the choice of representation on its own is not adequate to statistically break the degeneracy. By linearizing the system of equations, however, we gained additional insight to the problem. The linearized equations were the key ingredient that led us to seek additional constraints (e.g. virial theorem) that further reduce the feasible solution space.

Finally, we need to emphasize again the importance of using large datasets for model selection. When these are not available, GANs can be one starting point towards the correct solution. t-JEAnS -or any other algorithm -will fail in the absence of sufficient data.

The case of multiple stellar population dynamics

The linearization of the system of equations in the Jeans formalism yields some useful insights for the case of multiple stellar populations. When it is feasible to separate the stellar population into multiple stellar sub-populations (assuming two for simplicity) that are evolving under the influence of the same DM potential, the system of equations describing the system becomes: This is a set of four equations, with five unknowns (assuming, for simplicity, that the stellar tracer densities, ρ (1,2) , are known), namely, σ 2 (1,2)rr,tt , ρ•(r). The system of equations is still not closed (in fact, irrespective of the number of sub-populations, we will always have one more unknown function than equations). However, from the insight we get from the linearized equations (say, using B-splines), we understand that if the profiles of the stellar populations are significantly different (i.e. the determinant of the linearized system is not zero), then the solution space is reduced significantly. Depending on the statistical uncertainty of the observables, this may be enough to accurately describe the underlying DM structure. In contrast, if the profiles of the sub-populations are identical, the linear systems are identical (their determinant is zero) and no additional reduction of the feasible solution space is possible. Clearly, the linearization of the equations with the use of B-splines (or other suitable complete bases, e.g. wavelets), besides being a useful numerical scheme, also allows us to gain further insight into the degeneracy problem.

σ 2 1los = 1 Σ1 (R)
For systems with multiple stellar populations where we are trying to deduce more than one kinematic profile from scarce data, the GAN synthetic data generation can be a game changer for the estimation of the different brightness and LOS velocity dispersion profiles. The reason being that it can construct robust velocity dispersion data with small uncertainties over the extent of the system under investigation.

CONCLUSIONS

In this work we describe a new method for reliable mass determination independent of the mass-velocity anisotropy degeneracy. The efficiency of our method is tested on synthetic data from the Gaia Challenge suite of mock simulations. In all cases our algorithm reconstructs accurately the underlying kinematic profile as well as the mass content of the datasets. Our method includes: a) a new way of solving numerically the Jeans equations, subject to physically plausible local and global constraints, using quadratic programming. b) a new way for performing supervised learning in the framework of Jeans mass modelling, using samples from line-of-sight velocity dispersion MCMC chains as "unseen" validation and test data sets. Based on this, we present a new approach in performing regularization and model selection. c) The application of Generative Adversarial Networks for augmenting datasets, thereby making the t-JEAnS moments solver method reliable in situations where the available samples possess a relatively small number of stars.

Figure 2 .

 2 Figure 2. The t-JEAnS algorithm.

1Figure 3 .

 3 Figure 3. Comparison of NonPlumCoreOM 1k, 10k true datasets against the GAN generated synthetic data on the (R, v los ) plane.

Figure 4 .Figure 5 .

 45 Figure 4. Comparison of NonPlumCoreOM 1k, 10k true dispersion profiles against the GAN generated dispersion profile.Overplotted is the true reference profile. The GAN generated profile was created from information from the NonPlumCoreOM 1k dataset only.

3

  We use Einstein summation convention, where double repeated indices indicate summation. E.g. σ 2 rr (r) = a i B i (r) ≡ n basis i=1a i B i . 4 Free academic license.

Figure 6 .

 6 Figure 6. Exact numerical solution of the system of Jeans equations for the PlumCuspOM model using the Dynamic Moments Solver (DMS, see Table5). In order to obtain this solution we assumed perfect knowledge of the σ 2 los , ρ and ρ• profiles.

Figure 7 .

 7 Figure 7. As Fig. 6 for the PlumCuspIso reference profile.

Figure 9 .

 9 Figure 9. Same as Fig 8 for the PlumCuspIso 10k dataset.

Figure 10 .

 10 Figure 10. Same as Fig 8 for the PlumCuspTan 10k dataset.

Figure 11 .

 11 Figure 11. Same as Fig 8 for the NonPluCoreOM 10k dataset.

Figure 12 .

 12 Figure 12. Same as Fig 8 for the NonPluCoreOM 1k dataset. In the top right panel the σ 2 los data are from the GAN synthetic data generator. In the two bottom right panels the σ 2 rr and σ 2 tt moments are from the true NonPlumCoreOM 10k dataset (we do not have 3D motions for the GAN generated data).

Figure 13 .

 13 Figure13. Difference in fitted profiles in latent profiles σ 2 rr , σ 2 tt for the models trained on the GAN generated synthetic dataset and the 10k dataset, for the NonPlumCoreOM model.

Table 1 .

 1 Synthetic data sets: parameters α, β, γ are dimensionless numbers. Distance parameters r ,• are in kpc, while ρ

Table 2 .

 2 Competing

	DataSet	Stellar model DM model	χ 2 test
	PlumCuspOM, 10k	Plummer	Burkert	223.972
	PlumCuspOM, 10k	gH	NFW	217.337
	PlumCuspIso, 10k	Plummer	Burkert	207.001
	PlumCuspIso, 10k	gH	NFW	204.812
	PlumCuspTan, 10k	Plummer	Burkert	192.115
	PlumCuspTan, 10k	gH	gH	192.017
	NonPlumCoreOM, 10k	Plummer	NFW	348.126
	NonPlumCoreOM, 10k	gH	gH	340.255

mass models for the various 10k data sets. We report the average error on unseen test data, Dtest. The true models from which the data were produced are with bold fonts. In all cases the test error, χ 2 test , selects the correct model.

Table 5 .

 5 Summary of the DMS solver (DMS(θ|D) → (cr, ct)) in the JEAnS modelling approach. Assumptions: (a) Spherical symmetry, (b) virial equilibrium, (c) parametric form for the stellar and DM mass profiles. Input: θ ≡ {θ , θ•, ξ, λ 1,2 , β 1,2 } and training data, D = {R i , σ i los , δσ i los }. Output:

Table 6 .

 6 Competing mass models for the NonPlumCoreOM 1k data set, with (GAN) and without (1k) data augmentation. We report the average error on unseen test data, Dtest. The true models from which the data were produced are with bold fonts. The lower test error is also designated with a bold font. The recovery of the dataset is much better than what we would get by using only the NonPlumCoreOM 1k dataset. The recovered profile is of lower uncertainty than the one with the NonPlumCoreOM 10k dataset, especially close to the outer regions of the data. That is, the GAN generated dataset gives a better fit than the original True NonPlumCoreOM 10k dataset (note that a different range is displayed on the vertical axis in all right panels of Figures

	DataSet	Stellar	DM	χ 2 test
	NonPlumCoreOM, 1k	Plummer NFW	67.4542
	NonPlumCoreOM, 1k	gH	gH	69.8274
	NonPlumCoreOM, GAN Plummer NFW	497.068
	NonPlumCoreOM, GAN	gH	gH	492.709
	ues for reference.

In addition, there are well known methods for solving numerically systems of coupled integrodifferential equations, such as finite differences and finite element methods (

Ŝolín 2005; Jalali & Tremaine 2011), wavelets (Bertoluzza et al. 2008) and B-splines discretization (Höllig 2003).

http://astrowiki.ph.surrey.ac.uk/dokuwiki/doku.php

That is, the best B-spline basis, B i (x), according to the biasvariance trade-off(Hastie et al

. 2001).2 The lower the degree of the B-spline basis, the smaller the condition number of the system of equations.

Equivalently, the knots ξ i that define the simplest basis.

The description of each of the four smoothing parameters, θ smooth , is given in section 4.6.

http://astrowiki.ph.surrey.ac.uk/dokuwiki/doku.php

With the term hierarchical we mean models that result from the same general equation, but with possibly different complexity. Examples of hierarchical models are a Fourier expansion of a function:f (x) = n i=1 cn, cos(nx), or a B-spline basis, f (x) = n i=1 cnBn(x).As n increases we get models of increasing complexity that are derived from the same general equation.

ACKNOWLEDGMENTS

FID, GFL and CP acknowledge support from Australian Research Council Discovery Project (DP140100198). FID also thanks the University of Western Australia for support through a Research Collaboration Awards (PG12105204). GFL thanks the European Southern Observatory (ESO) for support as a visiting astronomer and for hosting him in Garching where the final stages of the preparation of this publication were undertaken. RAI gratefully acknowledges support from a "Programme National Cosmologie et Galaxies" grant.

APPENDIX A: PROOF OF THEOREM 1

Let us assume that there exist two radial profiles, σ 2 1rr and σ 2 2rr that give the same LOS dispersion profile. Then:

Subtracting the above equations, yields:

where ∆σ 2 rr = σ 2 2rr -σ 2 1rr . In order for this integral to be identically zero for all values of the parameter R, the integrand must be zero, i.e.

KA df dr

where we have set f (r) = ρ ∆σ 2 rr . For the case of R = 0 the result is trivial f = 0, i.e. σ1rr = σ2rr. For the case r > R > 0, we manipulate Eq (A2):

where A is the constant of integration, that will be determined from the virial boundary condition: since the last equation holds for all r, R, it will also hold for r = rvir and R = rvir/2, where rvir is the virial radius of the system. However for r = rvir, it is lim

Hence, limr→r vir σ 2 1rr (r) = limr→r vir σ 2 2rr (r) = 0, i.e. f (rvir) = 0. Then A exp(-4) = 0, i.e. A = 0, then f (r) = 0 and σ1rr = σ2rr for all r.

This proof is also valid for spherically symmetric systems subject to an external gravitational field: in this case as r → rvir both the radial and tangential velocity dispersions approach the same constant value [START_REF] Chan | [END_REF], thus again at the virial radius of the system f → 0. This paper has been typeset from a T E X/ L A T E X file prepared by the author.


Évelin Assis¹
email: evelin_assis@hotmail.com

Luciana Vellinho Corso²
email: l.corso@terra.com.br

Intervention in counting principles with first-grade students

Keywords: Counting principles, Intervention, Experimental study

The objective of this article is to present an experimental study of intervention in counting principles, aiming to investigate its efficacy. The experimental study involved 136 first-grade students, from 3 public schools located in Porto Alegre, Brazil. The students were divided into control group (n=76) and experimental group (n=60). The intervention occurred 2 times per week, during 2 weeks, with the experimental group, while the control group remained with regular instruction. The results point to the efficacy of the intervention program in counting principles, with students in the experimental group showing higher statistically significant advances in the post-test (counting principles evaluation task) than the students in the control group. Some educational implications are addressed.

Number sense and counting principles have been evidenced as important predictive factors of posterior mathematical achievement [START_REF] Jordan | The importance of number sense to mathematics achievement in first and third grades[END_REF][START_REF] Martin | Number and counting skills in kindergarten as predictors of grade 1 mathematical skills[END_REF][START_REF] Passolunghi | The precursors of mathematics learning: Working memory, phonological ability and numerical competence[END_REF][START_REF] Stock | Mastery of the counting principles in toddlers: A crucial step in the development of budding arithmetic abilities? Learning and Individual Differences[END_REF]. Besides, both, when less developed, appear as characteristics of the group of children that presents mathematical learning difficulties -MD [START_REF] Bull | Children's arithmetical difficulties: contributions from processing speed, item identification and short-term memory[END_REF][START_REF] Cirino | Cognitive and Mathematical Profiles for Different Forms of Learning Difficulties[END_REF][START_REF] Dorneles | Obstáculos cognitivos na aprendizagem matemática inicial: a contagem, as operações iniciais e os diferentes sentidos de número[END_REF]Vukovic, 2012). In regard to counting, studies show that children with MD use immature counting strategies (fingers or verbal counting) for a longer period [START_REF] Bull | Children's arithmetical difficulties: contributions from processing speed, item identification and short-term memory[END_REF][START_REF] Dorneles | Obstáculos cognitivos na aprendizagem matemática inicial: a contagem, as operações iniciais e os diferentes sentidos de número[END_REF] and have difficulties to retrieve basic facts from memory [START_REF] Cirino | Cognitive and Mathematical Profiles for Different Forms of Learning Difficulties[END_REF]. Research results such as these emphasize the importance of helping students to consolidate the counting principles in the beginning of the school years, which in turn, will help them develop efficient strategies and procedures to count efficiently and, thereby, prevent the development of arithmetical difficulties.

Counting is an important knowledge that begins to be developed informally by children, through everyday experiences. [START_REF] Nunes | Crianças fazendo matemática[END_REF] point that many teachers recognize their students get to the classroom knowing a lot of things, including counting. Its relevance is highlighted by several authors and [START_REF] Dorneles | Princípios de Contagem: uma construção progressiva[END_REF] synthesizes well this idea by stating that counting is a cognitive tool important not only to comprehend posterior contents, but also to develop more elaborated mathematics skills. Knowing the relevant role of counting, it is necessary to pay attention to the development of this ability. While children learn to count, they need to put in practice some principles that guide this action, as well as comprehend what is the purpose of counting and the way of doing it correctly [START_REF] Nunes | Crianças fazendo matemática[END_REF]. [START_REF] Gelman | The Child's Understanding of Number[END_REF] studied this learning process and established five guiding principles of counting, which must be consolidated by children in order to develop other mathematical knowledge. The principles are: stable order, one-to-one correspondence, cardinality, abstraction and order irrelevance. The first three principles are classified as "how to count", the fourth is denominated as a definer of what is countable and the fifth one consists of a synthesis of the application of the other ones. Those five principles are the bases of the children's posterior numerical construction [START_REF] Dorneles | Princípios de Contagem: uma construção progressiva[END_REF].

The current study

Considering the extreme relevance the counting principles have in relation to mathematical learning and that most children who have difficulties in this area present problems, specifically, with counting strategies and the retrieval of basic facts, justifies the need of the development and implementation of intervention programs focusing on counting principles. Such an approach, put into practice in the beginning of formal schooling, can facilitate the students' learning process. Besides, most intervention studies developed in the mathematical field focused on number sense (Aragón-Mendizabal, Aguilar-Villagrán, [START_REF] Aragón-Mendizábal | Improving number sense in kindergarten children with low achievement in mathematics[END_REF][START_REF] Bryant | Early numeracy intervention program for first-grade students with mathematics difficulties[END_REF][START_REF] Dyson | A number sense intervention for low-income kindergartners at risk for mathematics difficulties[END_REF][START_REF] Fuchs | The effects of strategic counting instruction, with and without deliberate practice, on number combination skill among students with mathematics difficulties[END_REF][START_REF] Praet | Enhancing young children's arithmetic skills through non-intensive, computerised kindergarten interventions: A randomized controlled study[END_REF], which points to the lack of intervention instructions in counting principles.

Method

This is an experimental intervention study focusing on teaching the counting principles to firstgrade students. It has the objective of evaluating the effectiveness of a short intervention, comparing the achievement of the experimental and the control groups.

Participants

The study included 136 first-grade students (63 girls and 73 boys), aged between 6 and 7 years old, from 10 classes, belonging to 3 public schools located in Porto Alegre, Brazil. They were divided in control group (n=76) and experimental group (n=60), both groups were statistically equivalent and, through the analysis of students' performance case by case, were formed by children with different levels of performance in a counting principles task. The parents of all participating students gave their authorization through the signing of the Free and Informed Consent Form and the Dissent Form. Approval for this study was obtained from the Research Ethics Committee of Universidade Federal do Rio Grande do Sul. The intervention was offered to the experimental group, by the researcher, while the control group remained with regular instruction. Unfortunately, it wasn't possible to follow teachers' regular instruction during the execution of this study, so we do not have data about the control group's learning, although, we hypothesize that there was less focus on mathematics and more on literacy, because this is the main concern of first year teachers in Brazil.

Procedures

In august 2017, a pilot study was conducted to evaluate the applicability of the intervention program, that culminated in its second version, the one used in this study. From September to December 2017, the pre-test, intervention and post-test procedures were done. The interval between the pre-test and the intervention was one month and between the intervention and the post-test was 20-30 days.

Instruments

Counting principles task (as pre and post-test): this task was administered individually, by the researcher, in an interview format (see Table 1). The instrument was used to evaluate the construction of the counting principles by the students. It is the same task used by [START_REF] Dorneles | Princípios de Contagem: uma construção progressiva[END_REF][START_REF] Dorneles | Obstáculos cognitivos na aprendizagem matemática inicial: a contagem, as operações iniciais e os diferentes sentidos de número[END_REF]. For each principle, the students were classified into three groups: principle consolidated (Yyes); in construction (IC); not consolidated (N). Some criteria were established to classify the children in each group because some of the questions (like 2, 3 and 5) are formed by more than one query. In this case, the child's answer was classified as "Y" if he/she answered all the solicitations in each question correctly; "IC" if only one part of the question was correctly answered; "N" if the student did not answer any part of the question correctly.

Principle:

Solicitation:

Stable order

The experimenter asks the child: "until how much can you count?" and tell him or her to count out loud until the number said.

One-to-one correspondence

The experimenter shows 10 chips in a row and asks "how many chips are there?".

Then, shows 10 chips misaligned and asks "how many chips are there?". This procedure is repeated with 15 chips.

Cardinality

By the end of counting 15 items, the experimenter asks "how many are there in total?

Can you give me 10?"

Abstraction

The experimenter asks "if you were counting 15 candies, would you count the same way you counted the chips?".

Order irrelevance

The experimenter asks the child to count the 15 chips, arranged in a row, in a different order, that is, initiating by another chip. Then, the child is asked to say how many chips would remain if the row was undone. After, the experimenter asks the child to count 8 chips of the set, separating them aside, and then, to count the remaining 7 chips of the row. By the end, the child has to say how many chips there are in the total.

Table 1: Counting principles task Intervention in counting principles

In regard to the intervention, whose development is described in detail in the poster submitted to the TWG2, some guiding principles are on its bases, as the ones presented by [START_REF] Clarke | Effective instructional strategies for kindergarten and first-grade students at risk in mathematics[END_REF]: engaging children's previous knowledge , through "warming up" exercises in the beginning of the sessions to help them establish connections between previously learned contents and the new ones; promoting interactions, through the gradual responsibility of the students with their learning processes by considering their previous knowledge and the knowledge needed to complete the tasks by themselves successfully; and providing mathematical verbalizations, in a way that it was possible for the students to use specific mathematical language and have the opportunity to speak and think mathematically.

Besides those ideas, the principles established by [START_REF] Fuchs | Principles for the prevention and intervention of mathematics difficulties[END_REF] also contributed to the development of the intervention program in counting principles: fast pace, varied activities and engaging, challenge the patterns of achievement, self verbalization and visual and physical representations (primary level, to be done with whole classes); adaptations easy to adopt so the teacher can incorporate it in the class routine, adaptations that are adequate for the targeted students as well as for other students (secondary level to help children with more difficulties); focus on the student as unity for instruction and decision making, intensive instruction and explicit instruction of basic skills (tertiary level, as an individual and intensive intervention).

The intervention sessions consisted mainly in providing children with opportunities to work with and practice the counting principles, helping them to be aware of how counting works and how to do it correctly. They occurred twice a week, during two weeks, for 20-35 minutes, including the dislocation time of the children. The sessions were made with small groups (5 students maximum), in a room provided by the schools and were conducted by the researcher see table 2 below.

Guiding questions:

Tasks:

Session 1: stable order, one-to-one correspondence and cardinality.

-where do we find numbers? -what are their purpose? -how can we count? -is there a right way to count? Which one? -what can we do to know how many things there are in a set?

-observation of the room to search for numbers and their purposes; -analyses of journals and magazines to search for numbers and then ordering them; -game: "which one is missing?"

Session 2: stable order, one-to-one correspondence, cardinality and order irrelevance.

-how can we count without getting lost? -how can we know if we already counted one item? -how can we count a set that is not organized?

-demonstration of counting marbles in figures by the researcher; -demonstration of counting marbles in figures using little papers; -individual practice.

Session 3: stable order, one-to-one correspondence, cardinality, abstraction and order irrelevance.

-what can we count? -where do we begin to count?

-game board: "color game"; -game board: "two hands".

Session 4: stable order, one-to-one correspondence, cardinality, abstraction and order irrelevance.

-how can we count in our heads? -how can we organize numbers when we'll count? -how can we find out what number is without knowing its name? -bingo of numbers; -closing questions: (can you say one thing you learnt during the intervention? Which task was easy? Which task was difficult?).

Table 2: Intervention in counting principles

Results

The data analyses were made through the Z test of proportion comparison. This test was chosen because the objective was to compare two proportions: the proportion of progress of the intervention group and of the control group. The data related to the intervention result will be presented next -value=0,0138). In activity C (abstraction), the students in experimental group also got superior improving rate (25,37%) when compared to the control group (17,07%), but this improvement did not show statistical superiority (p-value=0,1093). In activity D (order irrelevance), the experimental group once again evidenced a bigger improving rate (71,62%) than the control group (43,18%), being statistically significant (p-value=0,0001).

Discussion

Both groups showed advances, which was expected, since even the children that did not participate in the intervention process were having regular instruction, and so, they were learning. The results demonstrated that the improving rates of the experimental group were superior to the rates of the control group in all principles analyzed. Besides that, in three counting principles (one-to-one correspondence, cardinality and order irrelevance) this improvement was statistically superior. This result points to the effectiveness of the intervention program focusing on counting principles with short duration (2 sessions per week during 2 weeks). Besides that, it is important to highlight that the program effectiveness can be attributed not only to the tasks, but also to the type of teaching that was offered and the principles underlying the development of the intervention.

This finding is in agreement with other intervention studies that also compared control and experimental groups. [START_REF] Fuchs | The effects of strategic counting instruction, with and without deliberate practice, on number combination skill among students with mathematics difficulties[END_REF] developed an intervention in counting practice and evidenced that the three groups of participants (experimental with practice, experimental without practice and control) improved their achievement, with the children that received some kind of intervention showing more advances than the ones in the control group.

In the same way, [START_REF] Dyson | A number sense intervention for low-income kindergartners at risk for mathematics difficulties[END_REF] indicated that students with low income participating in a intervention condition in number sense showed more improvements than the subjects of the control group in the post-test. [START_REF] Praet | Enhancing young children's arithmetic skills through non-intensive, computerised kindergarten interventions: A randomized controlled study[END_REF] also pointed to similar findings: both experimental groups in their study demonstrated better achievement than the control group after getting intervention. [START_REF] Mononen | Counting skills intervention for low-performing first graders[END_REF] conducted an intervention research on counting abilities to children with low achievement. Their results indicated that the students in the experimental group improved their counting abilities significantly in comparison to the students in the control group, after the intervention. Other authors evidenced that specific interventions, designed to small children, involving one or more aspects of number sense, are effective and promote improvements in mathematical achievement (Aragon-Mendizábal et al., 2017;Bryant et al., 2010;[START_REF] Dowker | Target interventions for children with arithmetical difficulties[END_REF][START_REF] Dyson | A number sense intervention for low-income kindergartners at risk for mathematics difficulties[END_REF][START_REF] Fuchs | The effects of strategic counting instruction, with and without deliberate practice, on number combination skill among students with mathematics difficulties[END_REF][START_REF] Mononen | Counting skills intervention for low-performing first graders[END_REF][START_REF] Praet | Enhancing young children's arithmetic skills through non-intensive, computerised kindergarten interventions: A randomized controlled study[END_REF].

Implications for education

Our findings have important implications for teaching and learning. It highlighted a fundamental aspect of mathematical learning (counting principles), how to teach them (intervention) and how to help students to participate actively in their learning processes. It is urgent to help teachers to promote pedagogical actions that support students' learning process in relation to mathematics, area in which a lot of students show difficulties. It is important to act positively in this scenario: helping both students, through the teaching of specific strategies to facilitate learning, and teachers, through encouraging their monitoring of teaching and learning mathematics.

Considering these aspects, we believe the study presented here could reached some teachers' needs.

The intervention was developed based on guiding principles that need to be incorporated into regular instruction. Besides that, the intervention program used is brief, being developed in a short period of time, and could be included in the teachers' planning of instruction. The materials used by the researcher to conduct the intervention sessions were given to all teachers participating in the study, making it possible their further use by them. The results demonstrated that even a short investment in counting principles caused improvements in the children's performance, showing that it is possible to promote significant learning through concise and quality interventions dedicated to certain content area in the pedagogical planning.

Limitations of the study

The research results described in this paper should be considered in the context of some limitations. First, the time available to conduct this research was limited due to the fact that this study is part of the masters dissertation work of one of the authors and there was a deadline for its development and conclusion. Second, the counting principles task used was selected because it is a short and practical instrument, however, we believe that a broader task could elucidate more aspects to be analyzed and discussed. Another limitation concerns the choice of statistical analysis that was conducted. It is possible that different findings would have emerged if, besides the Z-test, the study included the association of different statistical tests, which could have made it possible to reinforce or question some of the results found or even to elucidate other possible outcomes.

Conclusions

The study demonstrated improvement in the children's performance in the counting principles task, from pre to post-test, with the participants showing good advances rates. This improvement was attributed to the intervention program because the students in the experimental group showed a bigger proportion of progress than the ones in the control group, showing statistical superiority in three of the four principles analyzed. Also, it is important to highlight the duration of the intervention: it was a short investment of time specifically devoted to the work with the counting principles, with carefully planned opportunities for the students to deal with and practice counting using the principles established by [START_REF] Gelman | The Child's Understanding of Number[END_REF]. The results of this study point to the need of continuing to investigate interventions that target counting principles which consist of a fundamental knowledge to be learnt in the beginning of formal schooling.

Table 3 .

 3 The stable order were not considered in the analysis because all children had consolidated this principle.

	Principle	Group	Progress (%) p-value
		Control	29	35,37 0,0016
	One to one correspondence
		Experimental 39	59,09
		Control	9	11,69 0,0138
	Cardinality
		Experimental 17	26,15
		Control	14	17,07 0,1093
	Abstraction
		Experimental 17	25,37
		Control	38	43,18 0,0001
	Order irrelevance
		Experimental 53	71,62

Table 3 : Results of the intervention

 3

It's possible to observe that in activity A (one-to-one correspondence), the improving rate in the principles construction by the experimental group was of 59,09%, being statistically superior (p-value=0,0016) in relation to the control group (35,37%). In activity B (cardinality), again the experimental group demonstrated a bigger improving rate (26,15%) than the control group (11,69%), being statistically significant (p