{'fraction_non_alphanumeric': 0.031125856725018552, 'fraction_numerical': 0.015017243637316104, 'mean_word_length': 4.912375790424571, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamical mass inference of galaxy clusters with neural flows

We present an algorithm for inferring the dynamical mass of galaxy clusters directly from their respective phase-space distributions, i.e. the observed line-of-sight velocities and projected distances of galaxies from the cluster centre. Our method employs neural flows, a deep neural network capable of learning arbitrary high-dimensional probability distributions, and inherently accounts, to an adequate extent, for the presence of interloper galaxies which are not bounded to a given cluster, the primary contaminant of dynamical mass measurements. We validate and showcase the performance of our neural flow approach to robustly infer the dynamical mass of clusters from a realistic mock cluster catalogue. A key aspect of our novel algorithm is that it yields the probability density function of the mass of a particular cluster, thereby providing a principled way of quantifying uncertainties, in contrast to conventional machine learning approaches. The neural network mass predictions, when applied to a contaminated catalogue with interlopers, have a logarithmic residual scatter which goes down to 0.043 dex for the most massive clusters. This is nearly an order of magnitude improvement over the classical cluster mass scaling relation with the velocity dispersion, and outperforms recently proposed machine learning approaches. We also apply our neural flow mass estimator to a compilation of galaxy observations of some well-studied clusters with robust dynamical mass estimates, further substantiating the efficacy of our algorithm. Such sophisticated approaches would undoubtedly be relevant for robust and efficient dynamical mass inference from upcoming surveys covering unprecedented volumes of the sky.


INTRODUCTION
The abundance of galaxy clusters is one of the most fundamental predictions of cosmological model and structure formation (Kravtsov & Borgani 2012). Counts of galaxy clusters as a function of dynamical mass are traditionally used to constrain the properties of dark matter (Rozo et al. 2010;Bocquet et al. 2019) and dark energy (Vikhlinin et al. 2009;Mantz et al. 2015), and test Einstein's theory of general relativity via the measurements of the linear growth (Rapetti et al. 2010). This is where the accuracy and precision of cluster mass estimates directly impact the robustness of cosmological inference, providing strong motivation for developing new techniques capable of handling constantly growing observational data.
Upcoming large-scale imaging and spectroscopic sur-ramanah@nbi.ku.dk veys, such as Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration et al. 2016), the Vera C. Rubin Observatory (Ivezic et al. 2008) and Euclid (Racca et al. 2016), will open up unprecedented opportunities for cluster cosmology with the potential of measuring the dark energy equation of state parameter to percent precision (Sartoris et al. 2016). The surveys will revive and challenge traditional methods of inferring dynamical cluster masses from relatively simple observables comprising the projected positions and velocities of cluster galaxies. The vast amount of future observations as well as the ongoing development of mock data from cosmological simulations will enable more effective dynamical mass inference based on neural network algorithms.
Masses of N-body systems have been traditionally measured from kinematic data (projected positions and velocities) using methods based on the virial theorem (e.g. Heisler et al. 1985). Keeping the core assumption of dynamical equilibrium and spherical symmetry, more complex approaches were successively developed. More general mass estimators derived from the Jeans equation allowed for incorporating broader assumptions about the underlying density profiles and the orbital anisotropy (An & Evans 2011). The increasing amount of observations enabled fitting complete solutions of the Jeans equation to the observed radial profiles of the projected velocity dispersion. This approach allowed constraining more detailed features of the systems such as density profiles or the orbital anisotropy (Biviano & Katgert 2004). The Jeans modelling was finally extended by including higher velocity moments ( Lokas & Mamon 2003) or by complete models describing the distribution of galaxies in the projected phase space (Wojtak et al. 2009;Mamon et al. 2013). This effectively exhausted theoretical possibilities permitted under the assumption of dynamical equilibrium and spherical symmetry.
Dynamical equilibrium and spherical models are highly idealized assumptions. They are not fully applicable to galaxy clusters which are not in perfect equilibrium (Natarajan et al. 1997), surrounded by infalling structures (Gunn & Gott 1972) and quite often exhibiting signatures of recent mergers. A way to circumvent this problem is to use nonequilibrium models linking the observed infall velocities to dynamical mass, e.g. the caustic method (Diaferio & Geller 1997) or the Jeans equation with non-vanishing streaming motion (Falco et al. 2014). Other approaches are based on exploiting phenomenological scaling relation between cluster dynamical mass and some easily measurable properties such as cluster richness (Rykoff et al. 2012). Further complications in dynamical mass estimations can arise from the fact that galaxy clusters are not spherical, neither in position space (de Theije et al. 1995) nor in velocity space (Wojtak 2013). This can likely degrade precision of the mass estimators based on kinematic data (Svensmark et al. 2015). Additional systematic errors are expected to be caused by substructures (Old et al. 2018) and the presence of interlopers (Wojtak et al. 2018). The latter effect appears to be the primary factor degrading the quality of cluster mass estimates. According to several tests performed on mock observations of galaxy clusters, dynamical masses or rich clusters with ∼ 10 2 member galaxies can be currently measured to a precision no better than 0.15 dex, quite often with a mass-dependent accuracy (Old et al. 2015;Armitage et al. 2019a), while the minimum attainable uncertainty for ideal measurements free of systematic effects related to imperfect interloper removal can be as low as 0.05 dex (Wojtak et al. 2018).
Machine learning (ML) approaches are becoming increasingly popular for applications involving the estimation of dynamical cluster masses. The initial implementations relied on non-parametric methods, where an ML model is optimized using a large training data set typically consisting of simulated mock observations. Once trained, the ML model may then be utilized to estimate cluster masses from unlabelled data sets. For instance, Ntampaka et al. (2015Ntampaka et al. ( , 2016 used support distribution machines (Sutherland et al. 2012) to directly infer the cluster mass from line-of-sight velocities and positions of galaxies. Armitage et al. (2019b) recently implemented a series of simple ML regression models, such as linear, ridge and kernel ridge regression, on a set of features extracted manually from cluster observations. This was followed by the work of Calderon & Berlind (2019), who implemented three more complex ML algorithms, namely XGBoost, random forests and neural network. The above ML approaches all led to similar significant improvements relative to the standard virial (M − σ) scaling relation, effectively reducing the prediction scatter by a factor of two, thereby highlighting the potential of ML-based approaches as promising alternatives to classical methods of cluster mass estimation.
Convolutional neural networks (CNNs) were recently designed by Ho et al. (2019) to further improve upon the previous ML approaches. Since CNNs are especially effective for image recognition purposes, they made use of a kernel density estimator to generate 2D phase-space maps of individual cluster dynamics, which are fed to the CNN as input images. The training rationale is a standard regression over logarithmic cluster mass with a realistic mock simulation providing the training and test sets. This CNN approach yields an overall improvement by a factor of three over the classical M − σ estimators, with the network prediction scatter reducing to ∼ 0.127 dex for the more massive clusters. Ho et al. (2019) also provides an excellent in-depth review of the above different classes of ML algorithms applied to cluster mass measurements.
In this work, we present a novel dynamical mass inference algorithm, inspired by the recently developed framework of neural flows (Jimenez Rezende & Mohamed 2015; Germain et al. 2015;Papamakarios et al. 2017), with our method being complementary to the ML approaches outlined above. Additionally, however, our neural flow model yields the conditional probability distribution of the dynamical mass of individual clusters, given their respective phase-space kinematics, thereby providing a way to quantify uncertainties, rather than predicting only single point estimates. Our training, validation and test data sets are drawn from the simulated cluster catalogue from Ho et al. (2019), which emulates the physical artefacts encountered in practice, described in Section 2.3.
The remainder of this paper is structured as follows. In Section 2, we briefly describe the general problem of mass estimation of galaxy clusters and outline the generation of the mock galaxy cluster catalogues employed in the training and validation of our neural network. We review the conceptual underpinnings of the neural density estimators in Section 3, with particular emphasis on neural flows, followed by a description of our network architecture and training procedure in Section 4. We subsequently validate and showcase the performance of our neural network in Section 5, and follow-up by illustrating a few applications on real galaxy cluster data sets in Section 6. Finally, in Section 7, we summarize the salient aspects of our work and outline potential extensions to further refine our neural flow approach.

ESTIMATING DYNAMICAL MASS OF GALAXY CLUSTERS
In this section, we provide a brief overview of the general problem of cluster mass estimation, including the virial scaling relation, and outline the generation of the mock cluster catalogues used in subsequent sections.

Classical M − σ relation
The virial scaling relation provides a simplified means to estimate the cluster mass M via only a single summary statistic, the so-called galaxy velocity dispersion σ v , and is, hence, referred to as the M − σ relation. The classical form of this relation may be derived from the equivalence of kinetic and potential energy, as encoded in the virial theorem, under the assumption of static and spherically symmetric clusters. The M − σ relation is typically expressed as where M 200c is the cluster mass definition adopted in this work, corresponding to the mass enclosed in a spherical region of density 200ρ c , i.e. 200 times the critical density of the Universe, with M denoting solar mass units, h(z) is the dimensionless Hubble rate, and the two scaling parameters: σ v,15 characterizes the velocity dispersion of a cluster with mass M 200c = 10 15 h −1 M and α is the power law scaling exponent globally describing the spatial distribution of mass in a particular cluster. These two scaling factors are generally set to their best-fit values obtained via simulations (Evrard et al. 2008). The velocity dispersion, as a summary statistic, may be estimated as the standard deviation of galaxy velocities projected along a single line of sight, but this does not capture and exploit all the information from the dynamical phase-space distribution.
Galaxy clusters are not perfectly homologous systems. The presence of various effects breaking the assumption of homology gives rise to substantial scatter around the M − σ scaling relation. Some examples include physical features, such as dynamical substructure (Old et al. 2018), cluster triaxiality (Svensmark et al. 2015), halo environment (White et al. 2010) and cluster mergers (Ribeiro et al. 2011). Moreover, this prediction scatter is exacerbated by selection effects, such as incomplete cluster observations or presence of interlopers (Wojtak et al. 2018). To mitigate these selection effects when applying the classical M − σ relation, complex and sophisticated membership modelling and interloper removal techniques are required (e.g. Mamon et al. 2013;Farahi et al. 2016;Abdullah et al. 2018).

Dynamical phase-space distribution
In terms of observations, the galaxy cluster data typically comprise the essential features of the cluster kinematics. For each cluster, the observables are the positions and velocities of its member galaxies to the cluster centre, computed by considering a line-of-sight axis. The positions of all the galaxies around the cluster centre are projected onto the plane (x, y) of the sky and are denoted by (x proj , y proj ). The net line-of-sight velocity, v los , for each galaxy corresponds to the sum of its relative peculiar velocity and the Hubble flow along the line-of-sight axis. This set of three observables are expressed as relative values to the cluster centre. The projected radial distance, R proj , defined as the Euclidean distance from the cluster centre, is derived from the planeof-sky positions as R proj = (x 2 proj + y 2 proj ) 1/2 . For the computations of x proj , y proj , R proj and, v los for a given cluster-  Figure 1. A sample of a galaxy cluster containing both true members and interloper galaxies. Left panel: The joint distribution of projected galaxy distances from cluster centre, R proj , and line-ofsight velocities, v los , in dynamical phase space. The true cluster members and interlopers are indicated by blue and red dots, respectively. Right panel: The corresponding 2D binned distribution which serve as inputs to our neural network (cf. Section 4.2).
galaxy pair, we refer the interested reader to Appendix A in Ho et al. (2019).
In this work, we will employ the set of (R proj , v los ) observables, which define the dynamical phase-space distribution of a galaxy cluster. An example is displayed in the left panel of Fig. 1, which also shows some interloper galaxies. Interlopers are non-member galaxies positioned along the line of sight, with similar observed line-of-sight velocities to the host cluster. The contamination induced by such interlopers is one of the main difficulties involved in cluster mass estimation (e.g. ). Since the clusters will have a varying number of galaxy members or interlopers, we preprocess the mock cluster catalogue, generated as outlined in Section 2.3, by first performing a 2D binning of the joint phase-space distribution of {R proj , v los }, as illustrated in the right panel of Fig. 1. This 2D binned distribution subsequently serves as inputs to our neural network (cf. details of architecture and training in Section 4). The extents of the 2D histogram are as follows: R proj ∈ [0, 1.6] h −1 Mpc, v los ∈ [−2200, 2200] km s −1 , with 50 bins along a given axis.

Mock cluster catalogues
The mock cluster catalogues are derived from a snapshot (z = 0.117) of the MultiDark Planck2 (MDPL2) N-body simulation 1 (Klypin et al. 2016). MDPL2 is a large dark matter simulation, tracing 3840 3 particles from an initial redshift z = 120 to z = 0, in a cosmological volume of (1 h −1 Gpc) 3 and mass resolution of 1.51 × 10 9 h −1 M , carried out using gadget2 (Springel 2005). The underlying cosmology is set to Planck ΛCDM best-fit values (Planck Collaboration et al. 2014): A detailed description of the mock generation procedure is provided in Section 2 in Ho et al. (2019). The main steps are summarized as follows: • Clusters and their galaxy members are modelled as host halos and subhalos, respectively, which are identified in the MDPL2 simulation using the rockstar halo finder  (Behroozi et al. 2013). Clusters are assigned to host halos above a certain mass threshold (M 200c ≥ 10 13.5 h −1 M ), thereby defining their physical properties such as mass, radius, position and velocity. Galaxies, in turn, are associated to subhalos via the galaxy assignment scheme of Uni-verseMachine (Behroozi et al. 2019), with a selection criterion on their mass at accretion (M acc ≥ 10 11 h −1 M ), inheriting their corresponding positions and velocities.
• The mock cluster observations incorporate realistic systematic effects typical of dynamical mass measurements. This includes physical effects such as cluster mergers and triaxiality in the distributions of cluster members, and selection effects, such as the presence of interloper galaxies. Modelling such effects requires selecting member galaxies by extracting a fixed cylindrical volume positioned at the centre of the cluster, aligned with the line-of-sight axis. This is usually referred to as an observational cut.
• Given the above line-of-sight axis, cluster membership may be determined by computing the position (x proj , y proj ) and velocity (v los ) observables, introduced in Section 2.2, for each cluster-galaxy pair. The cylindrical cut is specified in terms of its radial aperture size and half of its length along the line-of-sight axis, with their adopted values being R aperture = 1.6 h −1 Mpc and v cut = 2500 km s −1 , respectively.
• Two catalogues are produced from the MDPL2 simulation: pure and contaminated, with the primary difference being the presence of interlopers in the latter more realistic version. For the pure catalogue, galaxies located within the virial radius of a particular cluster are included in the mock observation. For the contaminated version, galaxies falling within the cylindrical cut are associated to the host cluster, irrespective of whether they are gravitationally bound to the system. For both catalogues, as a minimum richness criterion, clusters with fewer than ten galaxies are discarded.
The cluster mass function, describing the abundance of clusters as function of mass, of the respective training and test sets, extracted from the simulated pure and contaminated catalogues, is illustrated in Fig. 2. In order not to induce a selection bias during training, we generate a training data set with a flat mass distribution by downsampling the clusters at low masses, which are the most abundant relative to the scarce most massive ones. This may be interpreted as a maximally agnostic prior so as not to encode any cosmological information in our neural network. The test sets, however, follow the theoretical halo mass function that was used in the mock generation to provide a realistic setting to evaluate the performance of our neural network.

NEURAL DENSITY ESTIMATORS
A neural network, in essence, is a trainable and flexible approximation of a model, M(θ, γ) : d → τ, to map some input data d to an estimate of the desired label or target τ associated with the data. It is parameterized by a set of weights θ, which are tuned via stochastic gradient descent to optimize a given cost or loss function, and a set of hyperparameters γ, which comprise the choice of network architecture, initialization of the weights, type of activation and loss functions.
Density estimation, in its simplest form, entails the problem of estimating the joint probability density P(x) of a set of variables x from a set of examples {x i }. This joint density is crucial to perform a variety of tasks, such as prediction, inference and data generation, and therefore, constitutes a key aspect of probabilistic unsupervised learning and generative modelling. The use of neural networks for density estimation is becoming increasingly popular by virtue of their flexibility and learning capacity and this has led to the emergence of the so-called neural density estimators. These tools have been particularly successful at modelling natural images (e.g. Dinh et al. 2016;Salimans et al. 2017;van den Oord et al. 2016a) and audio data sets (van den Oord et al. 2016b).

Conditional density estimation
Neural density estimators provide a flexible parametric framework to model conditional probability densities, P(x| y; θ), 2 where θ corresponds to the weights of the neural network, trained on a set of simulated data(-parameter) pairs {x, y}. Unlike other approaches to generative modelling, such as variational autoencoders (Kingma & Welling 2013) and generative adversarial networks (Goodfellow et al. 2014), neural density estimators readily provide exact density evaluations. Ideally, such an estimator should be sufficiently flexible to represent complex distributions and straightforward to train, while having a tractable density function. There are two main classes of neural density estimators which satisfy these conditions: autoregressive models (Uria et al. 2016) and normalizing flows (Jimenez Rezende & Mohamed 2015).
The latest advances in deep learning have resulted in the emergence of a particular class of techniques, the so-called neural autoregressive flows, combining normalizing flows and autoregressive models (e.g. Kingma et al. 2016;Papamakarios et al. 2017;Trippe & Turner 2018;Huang et al. 2018), which have produced state-of-the-art performance in density estimation. We intend to explore a potential extension of our work, as outlined in Section 7, to encode such autoregressive flows in a future investigation.
This work primarily deals with a conditional density estimation problem, which has garnered significant interest in the field of deep learning. Conditional density estimation, in contrast to the above, is a supervised learning problem, encompassing both classification and regression, where the underlying aim is to predict the distribution of a set of labels x conditional on observing an associated set of features y. Here, we employ normalizing flows for conditional density estimation for dynamical mass inference, which are reviewed in Section 3.2.
Given their very recent conception and subsequent development in the deep learning community, the use of neural flows for potential astrophysical and cosmological applications remain as yet relatively unexploited. Cranmer et al. (2019) employed normalizing flows to model the probability distribution in stellar colour-magnitude space from parallax and photometry measurements. Neural autoregressive flows encoded in likelihood-free inference framework have been used to derive cosmological constraints from supernovae data (Alsing et al. 2019) and to constrain binary black hole systems using gravitational wave data (Green et al. 2020). Such neural flows have also been used within a Bayesian hierarchical model for gravitational wave population inference (Wong et al. 2020).

Normalizing flows
A normalizing flow encodes a smooth invertible mapping between probability density functions via a differentiable, monotonic bijection between the two spaces in which the functions live, i.e. R n → R n (Tabak & Turner 2013). The two fundamental requirements are that the transformation must be invertible and that the associated Jacobian must be tractable. The key property of a normalizing flow is that it is composable, such that the composition of a series of relatively simple invertible transformations, applied to a given base distribution, may sufficiently characterize more complex distributions. Jimenez Rezende & Mohamed (2015) demonstrated that a series of certain families of parametric transformations may warp a standard Gaussian base distribution into a relatively complex target density.
A normalizing flow (NoF), as introduced above, characterizes P(x) as an invertible differentiable transformation F of a base distribution Ψ(u), such that x = F (u), where u ∼ Ψ(u). The rationale behind the choice of the base distribution is that it should be easily evaluated for any input u and hence, it is typically chosen to be a standard Gaussian distribution. 3 Assuming that the transformation F is invertible, we can compute P(x) via a change of variables as which implies that F must be invertible and that it should be possible to compute the determinant of its Jacobian for equation (2) to be tractable, substantiating the two key requirements laid out above. An interesting point is that if two transformations, F 1 and F 2 , satisfy both conditions, then so does the composition F 1 • F 2 . As such, subsequently stacking multiple instances of F would still preserve the properties of NoF, while resulting in a deeper and more flexible model.
The invertibility of F implies that a NoF allows both sampling and probability density evaluation of P(x), as long as these tasks are feasible for the base distribution Ψ(u). To sample from P(x), we first draw samples u ∼ Ψ(u) to then obtain the desired samples via x = F (u). For density evaluation, we must inject the inverse mapping u = F −1 (x) in equation (2).

Conditional neural flows
Conditional density estimation, in essence, involves estimating the conditional density distribution P(x| y) from a given set of data pairs {x n , y n }. The unconditional NoF framework may be adapted to parameterize a conditional distribution simply by adding y to the set of input variables. Our choice of NoF design is inspired by the Masked Autoencoder for Distribution Estimation 4 (MADE, Germain et al. 2015), which is typically used in autoregressive models.
The NoF framework, in a nutshell, learns the means and variances of the conditional distributions. In terms of architecture, the intermediate or hidden layers of one NoF block may encode non-linear activation functions (e.g. tanh, sigmoid, ReLU). The output nodes, however, must obey certain constraints, with the nodes corresponding to the conditional means having a linear activation, whilst those associated with conditional variances must have exponential activations to ensure positivity. Our neural flow model is a stack of multiple NoF blocks of the same type to yield an overall NoF with higher flexibility than the original one, i.e. F ≡ F 1 •F 2 •. . . F k for k th component in the neural flow. The output of each NoF block is fed to the next one, along with the conditional inputs (cf. Fig. 3). Since density evaluations are possible from equation (2), the neural network parameterizing the overall NoF may be trained using stochastic gradient descent to maximize the likelihood over the set of network weights θ that the data pairs {x n , y n } emanate from the model.
Our implementation reduces to a simplified version, 4 Our neural network architecture reduces to a particular implementation of NoFs inspired by the masked autoregressive flow (MAF, Papamakarios et al. 2017) which involves a stack of several MADEs. In our special variant of MAF, we discard the binary mask in the MADE block, introduced to preserve the autoregressive (chain rule of probability) property, since we work in terms of a single autoregressive conditional distribution (due to cluster mass M being a scalar). since dim(x) = 1 in our case, since we are inferring only the dynamical cluster mass M from the phase-space distribution {R proj , v los } ≡d, such that a given NoF block must learn solely the mean and variance of a single conditional (Gaussian) distribution. The conditional NoF may be interpreted as learning the transformation of the random variate M to the latent space u where we set the unit normal distribution, u(M,d; θ) ∼ N (0, 1), via where µ * and σ * correspond to the mean and variance, respectively, of the 1D conditional distribution, with µ * ∈ R n , σ * ∈ R n + , and θ are the trainable weights of the neural network. The Jacobian of the above invertible mapping is, therefore, trivial, such that we can express the conditional density estimator of a given NoF component, from equations (2) and (3), as Driven by the above rationale, we can implement the conditional neural flow model by stacking multiple such conditional NoFs, i.e. there is now an additional inputd to every layer (cf. Fig. 3). From a straightforward extension of the NoF conditional density estimator given by equation (4), we can express the conditional density estimator characterized by the overall neural flow model as where N nof denotes the number of NoF blocks in the overall neural flow, with k labelling each NoF component, and u out corresponds to the output from the final NoF block.

NEURAL FLOWS
In this section, we describe our implementation of neural flows, the network architecture and training rationale.

Neural network architecture
The underlying objective of our neural flow framework is to model the conditional probability density function (PDF) of the dynamical mass of a galaxy cluster, given its phasespace distribution characterized by the projected radial distance from cluster centre and the line-of-sight velocity, i.e. P(M |{R proj , v los }). The network takes as input a vectord, which is a flattened 2D histogram of the phase-space distribution, i.e.d ≡ {R proj , v los }, with an example illustrated in Fig. 1. The neural flow, therefore, models P(M |d), with the training data set consisting of pairs of {M,d}.
A general schematic of our neural flow architecture is depicted in Fig. 3. We employ three NoF components, with a Gaussian base distribution and tanh activation in the three hidden layers of each component, with 512, 64 and 512 neurons, respectively. As explained in Section 3.3, the output nodes associated to the conditional mean and variance must, respectively, have linear and exponential activations to guarantee positivity.

Training methodology
In order to fit a neural density estimator to our set of data pairs {M,d}, we must optimize the weights of the neural network to minimize the KullbackâȂŞLeibler (KL) divergence between the parametric neural density estimator, P(M |d; θ), and the target distribution,P(M |d) (cf. Alsing et al. 2019, for more details), as the loss function. A Monte-Carlo estimate of this KL divergence yields the negative log loss function as a sum over the samples in our training data set, as follows: where N samples is the size of our training set. Note that this is equivalent to the negative log-likelihood of the simulated data {M,d} under the conditional density estimator P(M |d; θ). In the case of a conditional NoF, with a Gaussian base distribution as given by equation (5), we may express equation (6) as where u out is the output of the final component, u out = (M − µ out )/σ out .
We train our neural flow model by minimizing the negative log loss function given by equation (7), with respect to the network weights θ = {θ j } for the j th hidden layer. In order to obviate risks of overfitting, we adopt a standard regularization method of early stopping in our training routine. To this end, we split the mock cluster catalogue into a training and validation data set. We designate an early stopping criterion of 200 weight updates, such that training is terminated when the validation loss no longer shows any improvement for this chosen number of consecutive training iterations, and the previously saved best fit model is restored.
The neural network and training procedure are implemented in TensorFlow (Abadi et al. 2016). We use the Adam (Kingma & Ba 2014) optimizer, with a learning rate of η = 10 −4 and first and second moment exponential decay rates of β 1 = 0.5 and β 2 = 0.999, respectively. The batch size is set to 100. We train the network for ∼ 2 × 10 3 weight updates (i.e. training iterations over a given batch), requiring around five minutes on an NVIDIA V100 Tensor Core GPU. We train two different networks, but with the same architecture, on the respective pure and contaminated catalogues. The pure catalogue provides the ideal scenario without the presence of interlopers, the main contaminant in the mass determination. Comparing the performance of the model in the ideal and realistic settings would provide some insights pertaining to the effectiveness of the neural flow to account for the spurious contaminations.

Normalizing flow (NoF)
... ... The first NoF component takes as input the cluster mass M and the corresponding binned phase-space distributiond to yield the conditional mean and variance, µ * and σ * , respectively. The output of each NoF block is fed as input to the next one, along with the conditional inputsd.

VALIDATION AND PERFORMANCE
We evaluate the performance of our neural flow mass estimator (hereafter NF mass estimator) on an independent test set extracted from the mock observations from the simulated pure and contaminated catalogues described in Section 2.3. The test set follows the cluster mass function (cf. Fig. 2) employed in the mock generation, unlike the flat mass distribution of the training set, in order to assess the performance of the model under realistic conditions. Fig. 4 depicts the model predictions for three particular galaxy clusters of low, intermediate and high masses selected from the mass range of our catalogue. We find that the inferred PDFs have remarkably small 1σ uncertainties, at the level of ∼ 0.05 dex and ∼ 0.11 dex for the most massive clusters in the pure and contaminated catalogues, respectively, and closely match the corresponding ground truth masses. The relatively larger uncertainties for the contaminated clusters result primarily from selection effects such as the presence of interlopers, encountered in practice, as explained in Sections 2.1 and 2.3.
An essential feature of our approach is that it yields robust uncertainties for individual cluster mass predictions. In contrast, the standard machine learning methods for cluster mass estimation, as reviewed in Section 1, yield point estimates and typically assess the scatter in their predictions relative to the ground truth. In Fig. 5, we illustrate the 1σ and 2σ confidence regions of our model predictions, against the ground truth masses, for 5000 clusters sampled randomly from the test set across the complete mass range. For the sake of comparison, for instance, with Fig. 7(b) in Ho et al. (2019), and for a fair representation of our network performance, the model predictions are binned in logarithmic mass intervals, labelled by m, with the respective means and standard deviations, averaged per bin, displayed in the top row of Fig. 5 via the coloured stars and shaded regions, respectively, for both the pure and contaminated catalogues. The corresponding bottom panels indicate the variation of the average 1σ and 2σ uncertainties across the different mass bins. Note that if we consider only point predictions from our model, then the logarithmic residual scatter in our model predictions would be given by the dash-dotted lines, which correspond roughly to the 1σ band depicted. For the most massive clusters, this scatter is smaller by a factor of three relative to that obtained with the 2D CNN in Ho et al. (2019). As expected, our neural flow model assigns larger uncertainties to the contaminated samples, but the mean predictions are all within the 1σ band, further substantiating the reliability and accuracy of our NF mass estimator even when dealing with the presence of interloper galaxies, on top of the physical contaminations which distort the galaxy cluster shape or mass distribution.
In Fig. 6, we compare the precision of the NF mass estimator and seventeen different methods of cluster mass estimation based on galaxy data. The estimates of precision for the reference methods are from the Galaxy Cluster Mass Comparison Project, which delivered an extensive test of a wide range of cluster mass estimation techniques using two contrasting mock galaxy catalogues (Old et al. 2015). Cluster masses were estimated using different variations of methods based on measuring the number of galaxy cluster members (richness), the velocity dispersion, the distribution of galaxy positions (projected radii), the galaxy distribution in the projected phase space and abundance matching (see more details in Old et al. 2015). We show the results obtained for the mock galaxy catalogue generated with a semi-analytic model of galaxy formation, which is expected to resemble closely the mock contaminated cluster data used in our study. The precision of cluster mass estimation is quantified in terms of the total scatter about the best-fit power-law relation between the true and estimated cluster masses. Following Wojtak et al. (2018), we split the total scatter σ into statistical and systematic errors, as given by where N mem is the number of cluster members. Both σ stat and σ sys are determined by fitting the above equation to the logarithmic residuals in the cluster mass measurements from the test runs. Fig. 6 demonstrates that the NF mass estimator devised in our study outperforms the traditional methods. For the contaminated data, we find a statistical error of σ stat = 0.03 dex and a systematic error of σ sys = 0.124 dex, which are much below typical values attained by the tradi-  As a consistency test, this validates the capability of our neural network to properly recover the actual cluster masses and demonstrates the potential constraining power of our approach, with the 1σ uncertainties gradually reducing to ∼ 0.05 dex and ∼ 0.11 dex, for the most massive clusters in the pure and contaminated scenarios, respectively. As expected, the relatively larger uncertainties for the latter catalogue reflect the contamination by interloper galaxies.  Precision of cluster mass estimation methods NF (this work) − pure NF (this work) − contaminated abundance matching radius phase space velocity dispersion richness Figure 6. Comparison between the precision of our NF mass estimator and seventeen traditional methods deriving cluster masses by means of measuring richness, velocity dispersion, radial distribution of galaxies, projected phase-space distribution of galaxies and abundance matching (Old et al. 2015). The precision is quantified in terms of statistical error σ stat (statistical error for a cluster with 100 galaxies) and systematic error σ sys , as defined by equation (8). The neural network approach devised in this study, as indicated by the stars, outperforms the traditional techniques.
tional techniques. The statistical error is smaller by a factor of three relative to 3/( √ 2N mem ln 10) = 0.09 dex expected for the mass estimation based solely on the scaling relation with the velocity dispersion. This shows how much the constraining power increases when the full information on the projected phase-space distribution of galaxies is exploited instead of relying merely on the velocity dispersion. The main source of systematic error σ sys appears to be due to interlopers. When applying our NF mass estimator to the pure galaxy catalogues, the systematic error decreases to σ sys = 0.04 dex.

APPLICATIONS
We now apply the neural network trained on the contaminated catalogue to redshift data of several real galaxy clusters to infer their dynamical masses and eventually make a comparison with corresponding masses available from the literature. The main goal is to demonstrate the potential of the NF mass estimator to recover well confirmed dynamical mass measurements obtained for observed nearby galaxy clusters. For this purpose, we select well-studied low-redshift galaxy clusters with large samples of spectroscopically measured redshifts and robust dynamical mass estimates. The selected clusters include the Coma cluster, the lensing cluster A1689 and six rich galaxy clusters for which dynamical models yield consistent mass estimates for a few independent methods of removing relatively strong interloper contamination ).
Using the NASA/IPAC Extragalactic Database, 5 we compile all spectroscopic redshifts of galaxies observed in the field of the selected galaxy clusters. The primary sources of redshifts measurements include the 2dF Galaxy Redshift Survey (Colless et al. 2001), the Sloan Digital Sky Survey Abazajian et al. (2009), NOAO Fundamental Plane Survey (Smith et al. 2004), the ESO Nearby Abell Cluster Survey (Katgert et al. 1998), catalogue from Durret et al. (1998) and redshift data for A1689 from Lemze et al. (2009);Czoske (2004). The reference cluster mass measurements include the most updated dynamical estimates based on the Jeans analysis ( Lokas & Mamon 2003;Munari et al. 2014) and the virial mass estimator (Abdullah et al. 2020). They were converted to the overdensity parameter adopted in this study.
The individual galaxy clusters are preprocessed in a similar way to the preparation of the training set (cf. Section 2.2), i.e. their respective phase-space distributions are binned in a 2D histogram. The 2D binned distributions are then fed to our NF mass estimator, resulting in the inferred dynamical masses listed in Table 1. The masses are indicated as log 10 [M 200c (h −1 M )], with the literature values also shown for comparison. The inferred dynamical masses are remarkably in agreement with previous measurements, as illustrated in Fig. 7, although such network predictions rely on the accuracy and fidelity of the mock observations in the training set to replicate actual observations (e.g. Cohn & Battaglia 2020). Nevertheless, the fact that the network can reproduce well-established masses, such as for the Coma cluster, further showcases the efficacy of our NF mass estimator. Since we do not apply any interloper removal scheme to the contaminated catalogue for training nor the above real cluster observations, this implies that the neural flow is capable of building an internal representation of features to account for interlopers to an adequate extent. While a detailed investigation of the cluster masses is beyond the scope of this work, these preliminary results are very promising and demonstrate the potential of neural flows for dynamical mass inference from next-generation galaxy surveys.

CONCLUSIONS AND OUTLOOK
We have presented a first attempt at dynamical mass inference of galaxy clusters with robust uncertainty quantification using neural flows. We obtain very promising results by fitting our model to the 2D joint phase-space distribution of galaxy clusters, which consists of the projected radial distance, R proj , from cluster centre and the galaxy line-of-sight velocity, v los . Our neural network architecture is inspired by the novel advances pertaining to sophisticated neural density estimators. We have employed a normalizing flow with relatively few parameters, such that the model can be trained within a few minutes on a standard GPU, with the subsequent predictions of cluster mass PDFs being nearly instantaneous.
To demonstrate the potential of our neural network, we applied the trained model to a selection of real galaxy clusters and infer their corresponding dynamical masses. We find  Figure 7. Comparison between mass measurements obtained with our NF mass estimator and the reference mass estimates from the literature, based on the Jeans analysis or the virial mass estimator (cf. Table 1), for eight clusters of galaxies. The error bars represent their corresponding 1σ uncertainties. The golden star indicates the inferred mass of the Coma cluster.
that the neural network performs remarkably well, yielding inferred masses consistent with past measurements available from literature (cf . Table 1). This undoubtedly serves to highlight the efficacy of our neural flow model, showcasing it as an extremely promising tool for dynamical mass inference from upcoming galaxy surveys.
While we have implemented a network of relatively low complexity, there are, nevertheless, a series of interesting possibilities with a further level of sophistication that are worth exploring. Instead of viewing the problem in a condi-tional density estimation setting, we can use neural autoregressive flows to model the likelihood, P({R proj , v los }|M) ∼ P(d|M), whered corresponds to a set of informative and sufficient summary statistics extracted by a separate neural network, such as the Information Maximizing Neural Network (IMNN) developed by Charnock et al. (2018). The use of the IMNN would circumvent the use of standard summary statistics such as the velocity dispersion σ v , which, as mentioned previously, does not account for all the relevant physical effects, such that it is inadequate in practice. This would, subsequently, allow us to work in a Bayesian framework to infer the desired posterior as P(M |d) ∝ P(d|M)×P(M), where P(M) is the corresponding prior on the dynamical mass M which could be physically motivated to account for effects such as presence of interlopers.
Another plausible and exciting avenue is to make use of a conditional pixelwise probability estimator (Lanusse et al. 2019) to work directly at the level of 3D cluster dynamics, i.e. {x proj , y proj , v los }. By utilizing 3D convolutions in this uncompressed space, the neural network will, in principle, inherently account for the presence of interlopers in a given galaxy cluster more effectively than the NF mass estimator presented in this work. As we develop more sophisticated neural mass inference algorithms, it also becomes essential to model the (epistemic) uncertainty associated with the weights of the neural network to better account for the regions of the parameter space with scarce training data (such as the most massive galaxy clusters in a catalogue). To this end, it is worth exploring recently developed techniques involving variational inference to model the network weights with some probability distribution (e.g. Blundell et al. 2015;Wen et al. 2018) instead of using only their maximum likelihood estimates when performing the mass inference. By eventually marginalizing over the learned distributions of the weights, we would ensure that the uncertainties associated to the inferred masses are not underestimated.