A MODEST review

We present an account of the state of the art in the fields explored by the research community invested in 'Modeling and Observing DEnse STellar systems'. For this purpose, we take as a basis the activities of the MODEST-17 conference, which was held at Charles University, Prague, in September 2017. Reviewed topics include recent advances in fundamental stellar dynamics, numerical methods for the solution of the gravitational N-body problem, formation and evolution of young and old star clusters and galactic nuclei, their elusive stellar populations, planetary systems, and exotic compact objects, with timely attention to black holes of different classes of mass and their role as sources of gravitational waves. Such a breadth of topics reflects the growing role played by collisional stellar dynamics in numerous areas of modern astrophysics. Indeed, in the next decade, many revolutionary instruments will enable the derivation of positions and velocities of individual stars in the Milky Way and its satellites and will detect signals from a range of astrophysical sources in different portions of the electromagnetic and gravitational spectrum, with an unprecedented sensitivity. On the one hand, this wealth of data will allow us to address a number of long-standing open questions in star cluster studies; on the other hand, many unexpected properties of these systems will come to light, stimulating further progress of our understanding of their formation and evolution.


Introduction
MODEST, which is an abbreviation for "Modeling and Observing DEnse STellar systems", was established in 2002 as a collaboration between groups working throughout the world on the dynamics of dense, collisional stellar systems. Originally, the emphasis was exclusively on theoretical and computational investigations, with the high-level goal of providing "a comprehensive software framework for largescale simulations of dense stellar systems, within which existing codes for dynamics, stellar evolution, and hydrodynamics could be easily coupled and compared to reality" [1] . Such a vision has quantitatively shaped the activities of the collaboration in the past fifteen years, with several projects pursued by a number of close-knit [1] The original MODEST mission statement and the description of many initiatives are recorded at http://www.manybody.org/modest/ working groups, during regular small-scale workshops. Many results and tools which have emerged from these initiatives have had a crucial role in defining the current interpretative paradigm for the formation and dynamical evolution of collisional stellar systems.
The MODEST community has also progressively inspired and enabled the creation of a stimulating arena for a broader discussion on the theoretical understanding and empirical characterisation of several classes of stellar systems, especially globular star clusters and galactic nuclei. Such a broadening of the scope of the collaboration was reflected in a steady growth of the number of participants to its periodic meetings, which have now become large-scale international conferences and have been organized in Europe, North and South America, and Asia.
The present article follows the trail of the long-standing, yet somehow intermittent, tradition of papers summarising the MODEST activities (Davies et al, 2006;Hut et al, 2003;Portegies Zwart et al, 2008;Sills et al, 2003). We, as a group of ten young MODEST members, wish to present an account of the state of the art in fields explored during the MODEST-17 conference [2] , which was held at Charles University in Prague, Czech Republic, in September 2017, thanks to the generous hospitality of the Prague Stellar Dynamics Group. The conference represented the primary gathering of the MODEST community in 2017, attended by more than a hundred participants, with a programme structured in nine sessions, spanning from numerical methods for the solution of the gravitational N -body problem to the next observational frontiers. Such a wide range of topics demonstrates the growing role of collisional stellar dynamics in a number of astrophysical areas.
The next decade will witness several new instruments which will deliver astrometric, photometric and spectroscopic data of unprecedented sensitivity and accuracy. The recent detection of gravitational waves opens a unique window on the realm of compact objects. Numerical simulations are expected to tackle progressively larger N -body systems and to include more realistic approaches to model a variety of multi-scale, multi-physics problems, benefiting from both software and hardware advances. This wealth of data and tools will allow us to address a number of longstanding open questions in star cluster studies, but also to formulate new ones. Many unexpected properties of these systems are already coming to light, challenging our understanding of their formation and evolution processes.
Despite the richness of the scientific programme [3] we are reporting on, this summary should not be considered as a complete assessment of the state of the art in collisional stellar dynamics. For a comprehensive view on theoretical and computational aspects we refer the reader to the books by Heggie and Hut (2003), , Merritt (2013), and Portegies Zwart and McMillan (2018). Several review articles are also available on: the internal dynamics of star clusters (Heggie, 2011;McMillan, 2015;Vesperini, 2010), young massive star clusters (Longmore et al, 2014;Portegies Zwart et al, 2010), multiple stellar populations in globular clusters (Bastian and Lardo, 2017;Gratton et al, 2012), star clusters in merging galaxies (Renaud, 2018), as well as the summary of a recent conference on the formation of globular clusters in a cosmological context (Forbes et al, 2018).
[2] https://modest17.cuni.cz/ [3] The electronic version of the book of abstracts, which includes also almost all posters, is available at: https://modest17.cuni.cz/doc/modest17-boa.pdf 2 Stellar dynamics 2.1 The new phase space complexity of old globular clusters The study of the internal dynamics of globular clusters (GCs) is traditionally pursued under a relatively stringent set of simplifying assumptions, chiefly isotropy in the velocity space, absence of internal ordered motions, and spherical symmetry. Thanks to new astrometric measurements provided by Gaia (see Section 6) and decades-long observational campaigns with the Hubble Space Telescope (HST), supplemented by targeted spectroscopic surveys (e.g., see Ferraro et al, 2018b;Gilmore et al, 2012;Kamann et al, 2018), we are currently undergoing an observational revolution, which is finally starting to reveal to us the phase space properties of many nearby Galactic globular clusters (e.g., see Gaia Collaboration et al, 2018b;Watkins et al, 2015). As a result, a new degree of kinematic richness is emerging in this class of stellar systems (e.g., see Bellini et al, 2017Bellini et al, , 2018Bianchini et al, 2018;Libralato et al, 2018). Along with this, the availability of powerful computational tools (see Section 3) calls for theoretical efforts on some forgotten aspects of collisional stellar dynamics, especially regarding the exploration of the role of angular momentum and anisotropy in the velocity space.
In this perspective, some recently addressed topics concern the implications of non-trivial initial conditions for numerical simulations of star clusters, specifically the evolution of spherical isolated systems with primordial anisotropy , and the effects of relaxing the assumptions of synchronization and coplanarity between internal and orbital angular velocity vectors (Tiongco et al, 2018). The interplay between the internal evolution of globular clusters and the interaction with the tidal field of their host galaxy also affects the evolution of anisotropy and rotation, generating interesting complexity in the kinematic properties of the clusters (Tiongco et al, 2017).
The kinematic properties of the outskirts of globular clusters have also been gaining growing attention. Recent observations show the existence of diffuse, spherical stellar "envelopes" surrounding a number of Galactic globular clusters (e.g., see Kuzma et al, 2018;Olszewski et al, 2009) and extending several times the nominal truncation radius (as estimated on the basis of simple, lowered isothermal models). The assessment of reliable star cluster members and the measurement of their kinematics, both in the plane of the sky and along the line of sight, poses great challenges due to the background confusion limit, but, in the few cases in which it has been accomplished at sufficiently large distances from the cluster centre, it has revealed an apparent lack of a truncation (e.g., see Marino et al, 2014). The nature of such stellar structures is therefore arousing great interest, and several possible formation (and disruption) scenarios have been proposed, including the presence of a population of "potential escapers" (see Claydon et al, 2017;Daniel et al, 2017;Heggie, 2001;Küpper et al, 2010), of tidal debris of accreted dwarf galaxies (e.g. Sollima et al, 2018), or stellar structures associated with a possible small dark matter halo (e.g. Ibata et al, 2013;Lee et al, 2017;Peñarrubia et al, 2017, Breen et al. in prep.). Such kinematic measurements are made particularly difficult also by the presence of binary stars. Unresolved binaries may determine an inflation of the velocity dispersion profile (e.g. Bianchini et al, 2016), which, in turn, can cause an overestimation of the dynamical mass of the system (e.g. McConnachie and Côté, 2010, Moyano Loyola, in prep.).

Recent progress on some fundamental aspects of collisional dynamics
In addition to the role played by the emerging "kinematic richness", our understanding of several fundamental aspects of collisional stellar dynamics is also far from complete. Fresh attention has been recently devoted to the study of the moment of core collapse in idealised N -body models. Approaches for a new definition of core collapse based on self-similar solutions have been developed (Pavlík andŠubr, 2018). Single component models of core collapse result in radial density profiles similar to those of previous simulations and in agreement with theoretical predictions. Multi-component models, on the other hand, result in time scales for core collapse that correlate with the times of formation of the first hard binaries (see also Fujii and Portegies Zwart, 2014). In this context, we note that Tanikawa et al (2012) demonstrated by means of direct N-body simulations that dynamically formed hard binaries usually originate from a strongly interacting group of four or more stars. This is in contrast with earlier claims that hard binaries originate from triple encounters.
New characterisations of the processes of relaxation and mixing have also been proposed (Meiron and Kocsis, 2018). While "relaxation time", as associated with the energy diffusion process, is a useful definition when considering the cluster as a whole, the concept of "mixedness" may give a better picture of motion diffusion when particles of particular orbital families are considered.

Towards the large N regime
Recent hardware and software advances (Nitadori and Aarseth, 2012) have finally allowed us to "solve" the gravitational million-body problem for selected globular clusters (see especially Heggie, 2014;Wang et al, 2016), but several challenges are none the less still present in realistic GC modelling, such as the number of particles moving towards N = 10 7 and primordial binaries contributing to the extreme time scale differences. Performing state-of-the-art numerical simulations such as those presented by Wang et al (2016) still require years of computing time on GPU (Graphics Processing Unit) supercomputers. The development of new algorithms for an efficient treatment of systems at the interface between collisional and collisionless dynamics is therefore urgently needed. Symplectic Particle tree and Algorithmic Regularization Code for Star Clusters (SPARC-SC) is a new particle-particle particle tree code designed for realistic massive star cluster simulations (Wang et al., in prep.). More details about the numerical developments can be found in Section 3.
Alternatively, "frozen" N -body models, which follow the dynamics of a tracer particle in the potential generated by N fixed particles, may serve as a tool to study the validity of the continuum collisionless limit (N → ∞). In this context, some recent results (Di Cintio et al., in prep) show that, consistently with previous work, the orbits evolved in frozen N -body potentials qualitatively resemble orbits integrated in the parent, smooth potentials. However, the dependence on N is nontrivial, and, concerning N -body chaos, the continuum limit may be questioned. Regular and chaotic orbits for large N models cannot always be safely extrapolated from systems with a lower number of degrees of freedom.
In the purely collisionless regime, a classical topic of current relevance is represented by the characterisation of a gravitational field by means of an expansion over basis functions. A new biorthogonal family of density-potential pairs has been recently proposed by Lilley et al (2018), with application to an efficient spectral decomposition of non-spherical potentials. Finally, two recent studies of barred galactic potentials have been conducted, with the identification of a class of threedimensional non-periodic orbits, which have boxy projections in both their face-on and side-on views , and the numerical exploration of the threshold for global disc instabilities (Valencia-Enríquez et al, 2017).

Globular cluster formation in a cosmological context
The current understanding of GC formation is still only partially satisfactory, both because of several outstanding issues arising from the modelling of the early phases of star formation in a clustered environment and because of the severe computational difficulties in connecting the cluster scale to a cosmological one.
Given the multi-physics and multi-scale characteristics of this problem, concerted computational efforts are needed to bridge the gap between stellar dynamical and hydrodynamical codes. The software environment AMUSE Portegies Zwart and McMillan, 2018;Portegies Zwart et al, 2013) was indeed developed to fulfil such a need (see Section 3). A novel coupling between AMUSE and FLASH [4] (Fryxell et al, 2000;Wall et al, 2017) provides new aspects to the traditional approach of star formation in GCs (McMillan et al, 2015). FLASH provides modules for magnetohydrodynamics, bringing feedback from radiation, supernovae, and stellar winds together with the dynamics and stellar evolution modules provided by AMUSE.
Several key questions that arise in theoretical studies of GC formation and evolution regard their link to the history of our own Galaxy, through the physical mechanisms that shape the population of clusters we see at z = 0. Could the products of regular cluster formation at high redshift have survived until the present day, and are these relics consistent with the properties of local globular cluster populations? In other words, are GCs old young massive clusters that have survived? To address these open issues, several groups are targeting the challenging problem of the formation of GCs in a cosmological context, with a variety of computational approaches (e.g., see the recent results by Carlberg 2018;Creasey et al 2018;Li et al 2017Li et al , 2018Renaud et al 2017;Ricotti et al 2016), or within broader galaxy formation projects (e.g., see E-MOSAIC -MOdelling Star cluster population Assembly In Cosmological Simulations within EAGLE, Crain et al 2015;Kruijssen et al 2011Kruijssen et al , 2012Schaye et al 2015;FIRE -Feedback In Realistic Environments simulations, Grudić et al 2017;Hopkins et al 2014;Kim et al 2018;CosmoGrid, Ishiyama et al 2013;Rieder et al 2013) 3 Numerical methods 3.1 NBODYX NBODY6 is a state-of-the-art direct N -body code specifically designed to study the dynamical evolution of collisional stellar systems (for details, see . The current release of the code has been upgraded to include the SSE (Streaming SIMD (Single Instruction Multiple Data) Extensions) and GPU compatibility (Nitadori and Aarseth, 2012).
The equations of motion of the individual stars are integrated by a fourth order Hermite predictor-corrector scheme (Makino, 1991) with individual time-steps coupled to a Ahmad-Cohen (Ahmad and Cohen, 1973;Makino and Aarseth, 1992) neighbour scheme. The major feature of this code is the implementation of regularisation techniques (Aarseth and Zare, 1974;Kustaanheimo and Stiefel, 1965;Mikkola andAarseth, 1990, 1993), which are essential for a proper treatment of close encounters. In principle, multiple methods may be combined to achieve an even higher performance, however, it should be noted that decision-making is the bottleneck of each calculation. For realistic modelling of stellar clusters, the code evolves single star and binary evolution from synthetic stellar evolution tracks (Hurley et al, 2000(Hurley et al, , 2002. The latest release of the code, NBODY7, contains a post-Newtonian treatment of the force calculation, which means that orbit variations and merging of compact object binaries due to general relativity can be simulated. Over the years, members of the MODEST community have also developed separate branches of the code with prescriptions to include a variety of additional physical ingredients, such as the effects of an arbitrary external tidal field (NBODY6tt [5] , Renaud et al 2011) and dynamical friction (NBODY6df [6] , Petts et al 2015).
In order to perform massive star cluster simulations, hybrid parallelization methods (SIMD + OpenMP, GPU) have recently been implemented (Nitadori and Aarseth, 2012). With a desktop equipped with state-of-the-art NVIDIA GPUs and 4 − 8 cores CPU, the simulations of N ≈ 10 5 systems are now achievable. In addition, NBODY6++GPU [7] adds the MPI support, thus million body-scale simulations become possible (see Wang et al, 2016) by using GPU computer clusters.
During the conference, one special session (which has been summarised in a dedicated document [8] ) has been devoted to an open floor discussion about maintaining and further developing NBODY6. Here we only mention that two tools have already been established to encourage the community to actively contribute to the development of the code and the support of the userbase: a wiki page [9] and a github repository [10] .
3.2 P 3 T Most of the dynamical simulations of star clusters in previous studies use the direct N -body methods with individual time steps (see previous Section). With a few hundred CPU cores and a few ten GPUs, realistic modelling of million-body GCs is now achievable (Wang et al, 2016), however, the high computational cost (i.e., approximately half a year per half-mass relaxation time) makes this kind of approach still impractical for general studies, such as specific simulations of rich clusters (especially in the case of highly concentrated systems) and broad initial parameter sampling. It is therefore crucial to devise new numerical strategies to [5] https://github.com/florentrenaud/nbody6tt [6] https://github.com/JamesAPetts/NBODY6df [7] https://github.com/nbodyx/Nbody6ppGPU [8] https://modest17.cuni.cz/doc/nbody_discussion_nup.pdf [9] https://github.com/nbodyx/Nbody6/wiki [10] https://github.com/nbodyx successfully tackle the "post-million body problem", e.g. by initially targeting the efficient investigation of stellar systems with N > 10 7 . Oshino et al (2011) invented and Iwasawa et al (2015) further explored the Particle-Particle Particle-Tree (P 3 T) method, which uses particle trees (Barnes and Hut, 1986) for long-range forces and individual time-step Hermite integrator for shortrange forces, based on an Hamiltonian splitting. They found that the P 3 T method can be much faster than a pure Hermite integrator. This result suggests a new direction for the development of the next generation of star cluster simulation codes. By adding regularization methods, the P 3 T code (e.g., SPARC-SC, which is under development by Wang et al.) may become an efficient alternative code for massive star cluster simulations, making it possible to reach not only N > 10 7 particles, but also more realistic stellar densities.
3.3 AMUSE AMUSE [11] (Astrophysical Multipurpose Software Environment; Pelupessy et al, 2013;Portegies Zwart et al, 2009) is a software framework for computational astrophysics. Existing codes from different domains, such as stellar dynamics, stellar evolution, hydrodynamics and radiative transfer can be easily coupled using a highlevel python script. It is based on a "kitchen sink" model, i.e. a realistic multi-physics and multi-scale simulation can be computed by using dedicated solvers for each physical process or component. AMUSE is not a monolithic code, but it interfaces existing mature numerical codes used by the community as modules to perform the actual calculations. The AMUSE interface handles unit conversions, provides consistent object oriented interfaces, manages the state of the underlying simulation codes and provides transparent distributed computing. Different codes in the same domain are incorporated by using the same interface specification, so that users can change a solver/integrator by changing one line of source code in the highlevel python script. This makes the architecture highly modular, and allows users to treat codes as building blocks and to choose the optimal ones according to the problem under consideration. A complete, hands-on introduction to this software environment can be found in Portegies Zwart and McMillan (2018).

Monte Carlo, Fokker-Planck and other methods
A number of approximate schemes have offered very valuable alternatives to the computationally demanding direct N -body techniques. Two main Monte Carlo approaches, which use a statistical method of solving the Fokker-Planck equation, were originally developed by Spitzer and Hart (1971) and Hénon (1971). The latter, so-called "orbit-averaged method" has been subsequently improved by Shapiro (1985) and Stodolkiewicz (1986), and, more recently, Joshi et al (2000) and Giersz (1998), respectively. Currently, two efficient numerical implementations rest on this distinguished legacy: CMC (Cluster Monte Carlo, Pattabiraman et al, 2013;Rodriguez et al, 2016c) and the MOCCA [12] code (MOnte Carlo Cluster simulAtor, Hypki and Giersz, 2013). Both implementations have been widely tested against direct N -body techniques, and have been used to explore the dynamical evolution [11] https://github.com/amusecode/amuse [12] https://moccacode.net/ of collisional systems within large parameter spaces of initial conditions (see the first results of the "MOCCA SURVEY Database", by Belloni et al 2016), an endeavour which is still unfeasible with direct methods. In recent years, these rapid, approximate schemes have been intensively exploited in particular to assess the astrophysical properties and local merger rate densities for coalescing binary black holes in GCs (see Askar et al 2017;Rodriguez et al 2015, and Sections 10, 11).
Finally, another approach based on the direct numerical integration of the orbitaveraged Fokker-Planck equation in energy and angular momentum space has been essential in the early formulation of the current evolutionary paradigm of collisional systems (see especially Cohn, 1979;Goodman, 1983). Along this line, Vasiliev (2017) has presented a new scheme for simulating the collisional evolution of spherical isotropic stellar systems based on the one-dimensional Fokker-Planck equation, with the publicly available code PhaseFlow [13] . This code implements a high-accuracy finite-element method for the Fokker-Planck equation, and can handle multiplecomponent systems (optionally with a central black hole, by taking into account loss-cone effects).

Dynamics of planetary systems in star clusters
Modelling planetary systems in star clusters poses many different challenges. Indeed, these systems have very different temporal and spatial scales, which leads to a hierarchical architecture. If one attempts to integrate the coupled system (i.e., coevolving planetary systems in star clusters), the integrator will be forced to use very small time steps (comparable to 1/10 of the orbital period of the shortest period orbits), which essentially stalls the simulation. Monte Carlo approaches can indeed effectively circumvent this difficulty, but the results will strongly depend on the quality of initial sampling. Besides, it is usually not possible to take into account various physical processes present in star clusters, such as stellar evolution, mass segregation, and tidal disruption. Also, the fact that the potential in star clusters is typically lumpy and may vary quickly in time makes it unreliable to impose a static tide to the target planetary system. Planet-planet scattering, a process which can greatly affect the dynamical evolution of a planetary system, should also be considered, since multi-planet systems are common. It is highly challenging to apply analytical treatments to multi-planet systems due to their chaotic nature.
Planetary systems can be affected by a dense stellar environment in two major ways: first, during the planet formation process (Stage 1), protoplanetary disks may be subject to truncation (e.g., Vincke et al, 2015) due to stellar encounters and/or photoevaporation due to the incident FUV photons from nearby O/B stars (e.g., Adams et al, 2006;Anderson et al, 2013). Second, as the disk dissipates (Stage 2), planets are no longer protected by the eccentricity damping mechanism of the disk, and their orbital eccentricities and inclinations can be effectively induced by stellar encounters (e.g., Cai et al, 2017;Hao et al, 2013;Li and Adams, 2015;Shara et al, 2016;Spurzem et al, 2009).
Despite the intrinsic challenges, much progress has been made in recent years. At Stage 1, simulations of disk truncations have been greatly simplified by the use of [13] PhaseFlow is provided within the Agama library by Vasiliev (2018a), available at https://github.com/GalacticDynamics-Oxford/Agama. empirical recipes generated from heuristic fits of hydrodynamical simulations (e.g., Breslau et al, 2014). First-order simulations that ignore the viscous evolution of protoplanetary disks already yield reasonable agreement with the observed disksize distribution , and a few follow-up studies that take into account the viscous evolution and photoevaporation are currently in progress.
At Stage 2, pioneer work has been carried out by Spurzem et al (2009). They use KS regularization (Kustaanheimo and Stiefel, 1965) implemented in NBODY6++ (Spurzem, 1999) to model single-planet systems in dense star clusters. Subsequently, Hao et al (2013) investigate multi-planet system in open clusters using Monte Carlo scattering experiments. They use the MERCURY package (Chambers, 1999) to integrate the long-term secular evolution of planetary systems, and the CHAIN package (Mikkola and Aarseth, 1993) to handle close stellar encounters. Li and Adams (2015) employ a similar approach to derive the cross-section for planetary systems interacting with passing stars and binaries.
A census of free-floating planets in star clusters is presented in Zheng et al (2015). By using NBODY6, Shara et al (2016) conclude that dynamical encounters excite the orbital eccentricities of planets to a sufficiently high level such that tidal circularization becomes particularly efficient at the perihelion, which in turn leads to the formation of hot Jupiters.
Finally, a recent study by Cai et al (2017) points out that external perturbations due to stellar encounters and internal planet-planet scattering are two major mechanisms responsible for the instability of planetary systems in star clusters. Furthermore, their results show that planet ejection is a cumulative process: only 3% of encounters are strong enough to induce the orbital eccentricities of the outermost planets by ∆e ≥0.5, and most encounters only result in minute changes in orbital eccentricities. A follow-up study by Cai et al (2018) argue that field planetary systems bear signatures from their parental clusters. For example, planetary systems with more planets tend to have cooler dynamical temperatures (i.e., lower orbital eccentricities and mutual inclinations), because they likely spend most of their time in the low-density outskirts of the parental clusters. As such, the orbital elements of field planetary systems can be used to constrain their birth environments.

Star forming regions
While young massive stars disperse and heat the surrounding gas, suppressing further star formation locally, their feedback might enhance or even trigger another star forming event further away. Expanding bubbles powered by HII regions, stellar winds or supernovae are able to sweep neutral or molecular gas over length scales of parsecs or tens of parsecs, collecting the gas into clouds, which cools down, gravitationally collapses, forms molecules and eventually stars. This mechanism, "collect and collapse", was first proposed by Elmegreen and Elmegreen (1978). This process also propagates star formation: under suitable conditions (mainly on the density), one single event of star formation may trigger a sequence of star forming events. Whitworth et al (1994) and Elmegreen (1994) provide analytical estimates for luminosity of young stars and density of the interstellar medium (ISM) under which star formation can propagate.
To study how widespread triggered star formation is, it is necessary to distinguish it from spontaneous star formation. A strong evidence for triggering is considered to be an observation of two or more star forming regions separated by appropriate time intervals, and a morphological structure suggesting that the regions are causally connected. Another piece of evidence comes from the comparison of the mass of the observed objects with analytical estimates predicted for observed density, luminosity and age. Observations (e.g., Deharveng et al 2010;Simpson et al 2012, Seleznev, in prep.) and numerical simulations (e.g., Dale et al, 2015) suggest that unambiguous identification of triggered star formation suffers from several limitations, mainly the uncertainty in dating young stellar objects and gas displacement due to feedback; however they indicate that approximately 1/4 to 1/3 of massive star formation is triggered.
Although the presence of filaments in molecular clouds had been known for decades, it was the Herschel Space Observatory which showed their ubiquity (André, 2015). Detailed maps provided by Herschel revealed several puzzling properties of filaments. There exists a critical line mass M lcr ≃ 2c 2 s /G for a filament of sound speed c s , above which the filament becomes gravitationally unstable. Non-star forming filaments with line mass < M lcr have been shown to be very common in molecular clouds, with some molecular clouds (e.g., Polaris Flare) being entirely composed of these filaments, and not forming stars at all. Another piece of evidence linking star formation to supercritical filaments stems from the observation of prestellar cores. The majority (≃ 70 %) of prestellar cores are found to lie within supercritical filaments. An unexpected property of filaments is their almost universal width of ∼ 0.1 pc (e.g., see Arzoumanian et al 2011, although see Panopoulou et al 2017 for a different point of view). The physical processes responsible for the observed properties of filaments are currently a subject of active research, with three suggested solutions: magnetic fields (Nakamura and Li, 2008;Seifried and Walch, 2015), large-scale turbulence (Mac Low and Klessen, 2004) and global self-gravity (Heitsch et al, 2008).
Supercritical filaments have a low Mach (c s σ NT 2c s ) non-thermal velocity component (Arzoumanian et al, 2013). Understanding of this component might shed light on dynamical state and perhaps also on the origin of filaments. Toci and Galli (2015) suggest that the non-thermal velocity component is caused by small amplitude Alfvén waves. However, more recently Di Cintio et al (2018) offer a simpler scenario, without invoking magnetic fields. They propose two models: cold collapse of an initially static filament and a perturbed filament. After initial virialisation, the models assume a state with an almost constant virial ratio throughout several hundreds of free-fall time. Both the models agree with the observational data.

Young stellar objects and the connection to their natal gaseous structures
Since prestellar cores form in supercritical filaments, the following questions are immediately prompted: for how long do prestellar cores trace the underlying gaseous structure? Which physical mechanism unbinds them?
With a multiwavelength study of the star forming region surrounding NGC 1333, Hacar et al (2017) recently found a complicated structure of fibers characterized by a non-trivial topology. At some places, the fibers are intertwined to bundles of filaments. As for the more general case of filaments, supercritical fibers are star forming, while subcritical fibers, which are the majority, are non-star forming. The age and evolutionary state of a young stellar object can be estimated from the slope of their spectral energy distribution. There are four classes of young stellar objects, with "class 0" being the youngest, and "class III" the most evolved (i.e., weak-line T Tauri stars). While class 0 sources closely trace the gaseous structures, the correlation with gas becomes less pronounced for more evolved sources; class III sources are located randomly with the respect to the gas distribution.
Even low-mass, young stellar objects can drive outflows, which may be sufficient to stop the accreting flow from setting the final mass of the young star. Alternatively, Stutz and Gould (2016) have proposed the "slingshot mechanism", which can also explain the termination of the accretion onto a protostar. In this mechanism, protostars are formed within an oscillating filament. The youngest protostars are well embedded within the filament because of their large gaseous envelopes. As the protostars evolve and become more compact, they progressively decouple from gas. When the decoupling occurs near the maximum of oscillation, where the acceleration is the highest, the protostar freely moves away from the filament, as confirmed by recent N -body simulations (Boekholt et al, 2017).

Young star clusters
Approximately 80% of young stellar objects are found in groups and clusters, with the remaining 20% being distributed in the field. The observed state of young star clusters can help improving our understanding the process of star formation by providing an empirical assessment of the boundary conditions of such process. This is particularly important in the case of massive stars, the formation of which is still poorly understood.
Two main channels for formation of young star clusters have been proposed: by means of the monolithic collapse of a massive, dense, molecular cloud, or by merging of many smaller subclusters (for a recent review, see e.g. Longmore et al 2014). The latter possibility implies that the resulting cluster may retain some memory of the original substructures, at least on a time scale which depends on the mean density of the region encompassing the subclusters themselves. In some cases, the boundary between these two scenarios might be blurred as, in sufficiently dense environments, subclusters can already start merging while star formation is still taking place.
Evidence of a hierarchical structure in the stellar component of the 30 Doradus region was reported by Sabbi et al (2012). Interesting results on the nearest young clusters containing massive stars were then obtained within MYStIX (Massive Young Star-Forming Complex Study in Infrared and X-ray Feigelson et al, 2013) survey. The youngest clusters are often elongated in shape, with typical eccentricities in the range 0.3 to 0.5, and they often contain several subclusters. Tentative evidence (Kuhn et al, 2014) suggests that young clusters then progressively evolve towards configurations that are approximately characterized by spherical symmetry and a more homogeneous density. Recent studies have also shown that hierarchical substructures, covering up to three orders of magnitude in surface density, are conspicuous in the Orion A and B molecular clouds (Gutermuth et al, 2011;Megeath et al, 2016). The empirical evidence of the hierarchical nature of Orion A and B also suggests that the precursors of massive clusters may likely be highly hierarchical as well.
The merging scenario of cluster formation has been recently investigated in a number of theoretical studies, with different numerical techniques. With direct Nbody models, Fujii et al (2012) have explored the hierarchical formation of young massive clusters via mergers of smaller clusters. They subsequently extended this study with an exploration of the products of mergers of turbulent molecular clouds, modelled with SPH simulations (Fujii, 2015;Fujii andPortegies Zwart, 2015, 2016). By means of collisional N -body models, Banerjee and Kroupa (2015) have set limits on the compactness of an initial configuration which can then form, on the time scale of Myr, smooth, spherically symmetric massive clusters similar to NGC 3603. By means of collisionless simulations, Grudić et al (2017) found that subcluster mergers tend to produce a cluster of projected density µ ∼ R −2 , which is in relatively good agreement with observational studies of massive star clusters, which have profiles of slope µ ∼ R −2.5 to µ ∼ R −3 .
Star forming regions have higher binary fraction of low-mass stars (typically by factor of 1.5 to 2), if compared to field stars or "exposed' star clusters. Thus, binary evolution is important even during the relatively short-lived embedded phase. This is mainly due to the fact that binaries are more easily destroyed in denser environments and that the stellar density abruptly decreases as the gas is expelled terminating the embedded phase. Using direct N -body simulations to model the destruction rate of binaries of different semi-major axes, Kroupa (1995) found that almost all solarmass stars are formed in binaries. A similar idea was also used to estimate the initial density in star clusters as a function of the cluster mass (Marks and Kroupa, 2012). The binary destruction rate in small clusters is faster than in more massive ones. In principle, the observed binary fraction of young star clusters can therefore offer a possible tool to distinguish between the two cluster formation scenarios mentioned above (Dorval et al, 2017).
If massive stars are formed predominantly at the cluster centre, this can leave an imprint in terms of the degree of primordial mass segregation. The knowledge of the birthplace of massive stars within a cluster is valuable for testing different theories of massive star formation (i.e., "competitive accretion', see Bonnell et al 1997), and even the origin of the initial mass function (Bonnell et al, 2007). Some young clusters (e.g. the Orion Nebula Cluster, ONC, see Hillenbrand 1997; NGC 3603 see Pang et al 2013; Westerlund 2, see Zeidler et al 2017) indeed show signs of mass segregation. Despite recent progress on observational side, this is still a matter of continuing debate because mass segregation can also arise dynamically from initially non-segregated clusters (either due to energy equipartition in the case of lower mass ONC, or due to merging of subclusters in the case of more massive star clusters).
The time scale for formation of a star cluster is another important open question. The age of individual stars in an embedded cluster can be determined by comparing observed colour-magnitude diagrams with pre-main-sequence evolutionary models. However, differential extinction, unresolved binaries and activity of premain-sequence stars introduce systematic errors leading to large uncertainties in the age determination. As a result, the age spread can range from ≃ 1 Myr up to 10 Myr, with a possible gradual increase of the star formation rate over time (Palla and Stahler, 2000). Recently, Beccari et al (2017), have detected three distinct pre-main-sequences in the ONC, which is still a very young and embedded cluster (see also Section 7).
Both increasingly more detailed observations and more sophisticated numerical tools, particularly regarding the interplay between stellar and gas dynamics (see Section 3) and the treatment of realistic feedback from massive stars, are needed in order to clarify the issues described above.
6 New observational frontiers 6.1 The (more than a) billion star surveyor: Gaia Among the datasets that are going to be delivered in the upcoming years, the most exciting for the star cluster community is arguably the one provided by the European Space Agency (ESA) astrometric mission Gaia [14] . The census of sources for which Gaia is expected to provide phase-space coordinates (positions, parallaxes, proper motions) is estimated to include ∼ 2.5 × 10 9 objects (∼ 2.5 times larger than what originally planned). In addition, other information will be provided, including stellar parameters, metallicity, etc., thus increasing the number of dimensions of the phase-space that we will be able to probe.
The Gaia-ESO survey [15] (Gilmore et al, 2012) was designed to complement Gaia, by providing accurate radial velocity measurements obtained with the GIRAFFE and UVES spectrographs at the European Southern Observatory (ESO) for more than 10 5 stars in both the Galactic field and star clusters. By considering Gaia and Gaia-ESO data simultaneously, it will be possible to improve the membership determination for stars, and thus to clean up the colour-magnitude diagrams of many nearby clusters. This will improve the determination of many global properties of star clusters and will enable a more detailed study of their internal dynamics.
The astrophysical complexity of these systems has become evident with these accurate data. For example, two kinematic components are identified in the γ 2 Velorum open cluster (Jeffries et al, 2014) with roughly equal numbers of stars having velocity offset of about 2 km s −1 : the population having a broader velocity distribution appears to be younger by 1-2 Myr and less concentrated than the other one. Another example is the detection of substructures in the open clusters Chameleon 1 (Sacco et al, 2017) and NGC 2264 (Venuti et al, 2017), which supports a formation scenario where clusters form from the evolution of multiple substructures rather than from a monolithic collapse (see also Section 5).
The first release of Gaia data (DR1-TGAS, Tycho-Gaia Astrometric Solution) provided positions, proper motions and parallaxes for ∼ 2 × 10 6 stars originally included in the Tycho mission in the solar neighbourhood. In spite of the relatively large uncertainties, many interesting results have been obtained from this release and its correlation with other surveys (SDSS, LAMOST, RAVE; Helmi et al, 2017;Myeong et al, 2017;Robin et al, 2017). For many open clusters, the mean parallaxes and proper motions have been determined using TGAS data (van Leeuwen et al, 2017), and they are generally in very good agreement with the earlier determination [14] https://www.cosmos.esa.int/web/gaia [15] https://www.gaia-eso.eu/ based on Hipparcos data (with the exception of the mean parallax for the Pleiades cluster). One of the main issue which emerged from the early analysis of DR1 data concerned the stellar parameters determination: covariances are apparent among the uncertainties of all parameter. Intense efforts have been devoted by several groups to the study of the correlation in the noise.
The second data release DR2 [16] has just provided data for a larger sample of stars, with an accuracy surpassing the expectations. DR2 contains positions (α, δ), G and integrated BP and RP photometric fluxes and magnitudes for all sources (with typical uncertainties ǫ G ∼ 0.002 mag), five-parameter astrometric solutions for all sources with acceptable formal standard errors (30 < ǫ π /µas< 700; for > 10 9 stars) and radial velocities for sources brighter than 12 mag (with ǫ v ∼ 2 km/s). Moreover, for stars brighter than G = 17 mag estimates of the effective temperature and, where possible, line-of-sight extinction are provided.
The community is currently busy mining and interpreting such a transformative dataset and several demonstration and early scientific results have already been published. As an incomplete list of possible applications, here we mention the quantification of the lumpiness of the local halo (Koppelman et al, 2018) (Malhan et al, 2018) and hypervelocity stars (Boubert et al, 2018;Hawkins and Wyse, 2018;Lennon et al, 2018;Marchetti et al, 2018). The quality of DR2 data is high enough to determine the fine structure of the Hertzsprung-Russell diagram of stars in the local neighbourhood and beyond (Gaia Collaboration et al, 2018a), an endeavour which is revealing many surprises for stellar population studies (e.g., see El-Badry et al, 2018;Jao et al, 2018;Kilic et al, 2018). As for the study of the internal kinematics of structures in the Local Group, this has been proved feasible even for individual point sources in M33, M31 and the Large Magellanic Cloud Vasiliev, 2018b), paving the way for the characterisation of many other objects, including selected young and old star clusters (Bianchini et al, 2018;Kuhn et al, 2018;Milone et al, 2018). Beside the primary interest on Galactic stellar populations and dynamics, Gaia data enable to measure exoplanet sizes (Fulton and Petigura, 2018), gravitationally lensed systems (Krone-Martins et al, 2018), and even to test general relativity using positional displacement of Sun and Jovian limb. Moreover, it can also be used to develop alternative astrometric search methods for gravitational waves sources (Moore et al, 2017).

The electromagnetic spectrum is not enough: gravitational wave detectors
The groundbreaking gravitational waves (GW) detection from Advanced LIGO [17] (Abbott et al, 2016a) provided the first evidence for massive black holes (BHs) mergers, confirming many predictions of general relativity.
[16] https://www.cosmos.esa.int/web/gaia/dr2 [17] https://www.ligo.org/ The continued improvement of sensitivity over time will result in more detections at greater distances. The advent of additional detectors around the world (such as KAGRA [18] , scheduled for late 2018) will result in better angular resolution and better chances for the identification of the electromagnetic counterparts of the detected events. Indeed, the joint detection of a GW event by LIGO and Virgo (The LIGO Scientific Collaboration et al, 2017) has provided a much better angular resolution than the detections by LIGO alone. Thanks to the joint collaboration between the different GW observatories it will be also possible to use polarization to break the degeneracy between luminosity distance and inclination.
New insights in this field are also expected from the forthcoming missions involving GW detectors in space. In this regard, it is noteworthy that LISA Pathfinder [19] , the first high-quality laboratory launched in December 2015 which conducted highprecision laser interferometric tracking of orbiting bodies in space, surpassed its performance requirements and expectations. LISA [20] will provide information about dynamical systems over a wide range of time, length, and mass scales. The advent of such mission will shift the boundary of GW astrophysics to low-amplitude/large distant events like those produced by Galactic white dwarf binaries (Finn et al. 2015), extreme mass ratio and supermassive BH inspirals (Gair et al, 2008), extra-Galactic binaries in the field and black hole binaries in GCs up to distances of 30 Mpc (Sesana, 2016).

A new eye on the invisible Universe: Lynx
A new era is also beginning for X-ray astronomy, thanks to Lynx, [21] a large area, high angular resolution, X-ray mission concept for the next decade. This mission will combine a large gain in collecting area, an angular resolution of 0.5 arcsec, and very high resolution spectroscopy over a large field of view. With its two orders of magnitude leap in sensitivity over Chandra and ATHENA [22] , it will provide insights into many different scientific problems.
One of the main science drivers of Lynx is the determination of the origin of supermassive BHs. It will be indeed possible to discriminate between the two main hypotheses regarding the seeds of these objects: Population III remnants and intermediate-mass black holes. Lynx will have enough sensitivity to directly detect 10 4 M ⊙ objects, pinning down the luminosity function to distinguish between these seeds. In this respect, Lynx would be a useful complement to JWST (James Webb Space Telescope), ATHENA, GMT (Giant Magellan Telescope) and other future missions.
Lynx is also expected to find and characterize the first generation of groups of Milky Way-sized galaxies around z ∼ 4. In this regard, it is crucial to investigate the content and ionization fraction of the circumgalactic (CGM) and intergalactic medium ( in galaxies at z = 0, but very different CGM and IGM properties. By observing feedback from star formation on scales from individual young star forming regions to the entire galaxy, Lynx will map the IGM, allowing us to uniquely pin down the physics of structure formation. Furthermore, medium-deep Lynx Surveys will expose the emergence of black hole populations in galaxies after z ∼ 6, in a wide range of galaxy types and density environments. Moreover, Lynx will provide a detailed view of every aspect of the BH feedback process, by measuring the energy and momentum flux in BH-generated outflows. This will make it possible to study where, how, and how much of the active galactic nuclei (AGN) outburst energy is dissipated in galaxies, groups and clusters.
Finally, with Lynx it will be possible to find and characterize the low L x population of GCs. With Chandra, we can probe below few ×10 30 erg/s only for a few nearby clusters (M4, 47 Tuc, NGC 6440, NGC 6752, NGC 6397), while with Lynx it will be possible to do this for dozens of Galactic GCs.
The contribution from the MODEST community to the development of the science cases within the upcoming 2020 Decadal survey has been encouraged. The Lynx observatory will offer a complementary view to dense stellar systems with respect to the facilities illustrated above and will thus be crucial to tackle the many open questions that keep emerging in this era of revolutionary observations.

Initial and present-day mass function
The original mass distribution of stars, commonly referred to as the initial mass function (IMF), represents one of the central questions in the theory of star formation and has strong relevance for many areas of astrophysics. The universality of the IMF, its shape and the parameters driving its hypothetical variation are still far from being completely understood.
Some suggestions on possible IMF variations have been recently put forward on the basis of the shape of the present-day MF in unevolved stellar systems (t age << t rh ). While Grillmair et al (1998) and Wyse et al (2002) derived MFs for Draco and Ursa Minor dwarf spheroidals which are consistent with a Salpeter (1955) IMF, Geha et al (2013) found evidence of MF variations correlated with the mean metallicity in a sample of ultrafaint dwarf galaxies. This study has been however questioned by El-Badry et al (2017), who found that evident MF differences cannot be detected with the currently available photometric data. Weisz et al (2013) analysed a large sample of young clusters and associations whose MFs are available in the literature: in spite of the large cluster-to-cluster differences, a careful revision of the associated errors indicates that the hypothesis that they are consistent with a single IMF slope cannot be ruled out.
Interesting information can be derived also by the study of dynamically evolved stellar systems. The analysis of high-resolution HST photometry for a sample of 35 Galactic GCs revealed tight correlations of the slope of the presentday MF with the half-mass relaxation time and with the fraction of remnants . The observed trends are compatible with the natural concept of dynamical evolution, with highly evolved clusters characterized by flatter MFs and large fractions of remnants, and the small observed spread leaves little room for IMF variation. However, none of the N -body simulations run for comparison are able to reproduce the flatness of the measured MFs (α > −1 for 27 out 35 analysed GCs) for their relatively low number of elapsed half-mass relaxation times (3 < t age /t rh < 10; Baumgardt and Sollima, 2017). Moreover, less evolved GCs (t age /t rh < 3) show a depletion of low-mass stars with respect to a single power-law, which should be reminiscent of their IMF shape (as predicted by Kroupa 2001 andChabrier 2003). In a similar comparison between observed MF and N -body simulations Webb and Vesperini (2016) found that while for the majority of clusters analysed in their sample the present day MF and the degree of mass segregation are consistent with the prediction of simulations starting from a universal IMF, some non-standard initial conditions should be present in NGC 5466 and NGC 6101, which are non-segregated and characterized by a flat MF. The need of IMF variation has been also put forward to interpret the evidence coming from mass-to-light ratios estimated through integrated light analyses of M31 GCs (Haghi et al, 2017) and the fraction of low-mass X-ray binaries in Virgo GCs and Ultra Compact Dwarf galaxies (UCDs, Dabringhausen et al, 2012).

Mass-to-light ratios
An accurate determination of the mass-to-light ratios (M/L) of star clusters is important because it provides direct information on their present day stellar population and important clues on their IMF. Understanding the values of M/L from a dynamical point of view is also crucial to determine the evolutionary history of these systems. However, estimates of the M/L from dynamics and from population synthesis do not fully agree, and the effects of several factors need to be considered in order to account for this discrepancy. In particular, the global mass function and the presence of binaries and dark remnants need to be properly assessed, and their determination is particularly tricky because they depend both on the unknown initial conditions and on the dynamical evolution history of the system. Not only the global M/L, but also their radial profiles within the system are crucial to uncover the stellar population content in star clusters. Indeed, mass segregation and the evolution towards a state of partial energy equipartition cause the M/L radial profile not to be constant, and to show a prominent central peak due to the presence of dark remnants . The study of the evolution of numerical simulations show that the exact shape of these profiles depend on the evolutionary stages of clusters. A decrease in M/L is seen for highly evolved systems, and it is primarily driven by the dynamical ejection of dark remnants, rather than by the escape of low-mass star.
A particularly intriguing puzzle is related to the decreasing trend of M/L K with metallicity observed in M31 clusters (Strader et al, 2011), which cannot be reproduced by means of simple single stellar population models. A proposed explanation for this is related to the role of mass segregation (Shanahan and Gieles, 2015), but this factor alone is not sufficient to completely account for the discrepancy. It has been recently showed that by taking into account dynamical evolution, together with a top heavy IMF, and a 10% retention fraction for remnants, a better agreement to the data can be obtained (Zonoozi et al, 2016), and to incorporate the age-metallicity relation further improves the fit (Haghi et al, 2017).

Multiple stellar populations
For decades GCs have been considered as a collection of stars all characterised by the same age and initial chemical composition, and treated as a prototype of a simple stellar population. However, numerous observational studies have put this simple picture to test: spectroscopic observations have unveiled different chemical composition for stars in the same cluster (Carretta et al, 2009;Gratton et al, 2004Gratton et al, , 2012, and accurate photometric data obtained with HST have shown the splitting of the evolutionary sequences in the colour-magnitude diagrams (CMDs) of GCs (Gratton et al, 2012;Milone et al, 2015;Piotto et al, 2015). One of the suggested interpretations of these observations is that clusters have been the site of multiple generations of stars with the second-generation (enriched) stars formed from material polluted by the ejecta of some first-generation (pristine) stars.
Several scenarios have been put forward to explain the formation of secondgeneration stars, each proposing a different source for the polluting material. However, none of these scenarios is able to reproduce the observed detailed chemical properties of clusters, and their origin is thus still unknown (Gratton et al, 2012). Each one of the formation scenarios proposed so far should have imprinted a typical signature in the spatial distribution and kinematic properties of different populations, even though two-body relaxation determines a progressive mixing (Vesperini et al, 2013). In particular, observations show that enriched stars (i.e., with low O and high Na abundances) are more concentrated (Lardo et al, 2011) and characterized by a lower fraction of binaries (possibly linked to their high concentration at early epochs; Lucatello et al, 2015) with respect to pristine stars.
Recently, observational efforts have been dedicated to other systems, such as clusters in nearby galaxies and young stellar clusters (see Bastian and Lardo 2017), in order to determine the importance of the environment and of the epoch of formation of these systems in shaping their populations. Theoretical studies have also explored the dynamical properties of such multiple populations, in view of the possibility of using high-quality kinematics to unravel clues about their origin (e.g., Hénault-Brunet et al, 2015;Perets, 2013, 2016;Vesperini et al, 2013). Li et al (2016) reported an age difference of a few Myr among stellar populations in three clusters of the Magellanic Clouds (NGC 1783, NGC 1806, andNGC 411). In the case of NGC 411, they have reported a reversed radial behaviour of the two populations with respect to what is observed for the majority of the Galactic GCs, with the enriched population dominating in the outskirts of the cluster. A minor merger (mass ratio less than 0.1) scenario seems to be appealing to reproduce this feature (see also Amaro-Seoane et al, 2013;Carretta et al, 2010;Gavagnin et al, 2016). A study of the dynamical properties of the merger remnants (Hong et al, 2017) shows that they are expected to be characterized by some degree of rotation and radial anisotropy; unfortunately, these kinematic properties are still beyond the current observational capabilities. It has been discussed that, by considering a different background subtraction compared to the approach by Li et al (2016), Cabrera-Ziri et al (2016) found no reversed radial behaviour for the two populations, but such a feature has been detected in other systems (see for example the case of M15, Larsen et al, 2015).
Moving on to even younger systems, three discrete sequences have been recently discovered by Beccari et al (2017) for pre-MS stars in the ONC. Among these populations, the older stars appear to be less concentrated than younger ones, and to rotate more slowly, hinting to a possible slowdown as they evolve. The origin of this discreteness is still uncertain. Several possible explanations have been proposed: the presence of unresolved binaries, rotation, or a genuine age difference. Only in the last case, the three sequences would correspond to different episodes of star formation.
The empirical characterization of the dynamics of the multiple population phenomenon is arousing progressively greater interest, and, when possible, attention is being devoted to the observational study of the phase space properties of distinct populations (Bellini et al, 2017Cordero et al, 2017;Libralato et al, 2018;Milone et al, 2018, see also Section 2.1). One additional recent example is represented by the case of M54, for which MUSE observations (Alfaro-Cuello et al., in prep.) provided an exceptional sample of kinematic and chemical data, allowing the identification of two distinct stellar populations (see also Bellazzini et al, 2008), with different metallicities, velocity dispersion and possibly rotation profiles. Such an object is indeed particularly complex, not just because it might constitute a "link" between classical globulars and nuclear star clusters, but also given its proximity to Sagittarius, which makes the membership assessment particularly challenging (especially regarding possible contamination of the younger population).

Stellar population in extra-Galactic systems
Recently, much efforts have been put also into the analysis of the dynamical and general properties of stellar populations in extra-Galactic environments. Full spectral synthesis analysis of nuclear star clusters (NSCs) in several galaxies reveals generally old ages for their stellar contents (∼ 10 Gyr), with some evidence of prolonged star formation until the present day (Kacharov et al., in prep.). However, every NSC appears to be characterised by a complex star formation history (e.g., see the case of NGC 5102 and Section 8).
Simple population synthesis has been recently used also to trace back the properties of progenitors of GCs and ultra-compact dwarfs (Jerabkova et al, 2017). From these preliminary studies it seems possible that these objects could be as bright as quasars at high redshifts and thus detectable with the upcoming JWST. In this regard, exciting evidence of young (t age < 100 Myr), low-mass (M < 10 7 M ⊙ ) and low-metallicity (Z < 0.1 Z ⊙ ) stellar systems compatible with the expected properties of the progenitors of present-day GCs have been recently revealed at redshift z ∼ 3.2 thanks to the high magnification produced by the gravitational lensing of the foreground galaxy cluster MACS J0416 (Vanzella et al, 2017).
Finally, considering GC systems as a whole, differences between the GC content of the Fornax and Virgo galaxy clusters are apparent, with GCs in the Virgo cluster being more concentrated, when observed in g-band, than those in Fornax, while the two distributions overlap when near-infrared luminosities are considered (Dabringhausen et al., in prep.).

The Milky Way Nuclear Star Cluster
The supermassive black hole (SMBH) at the centre of our Galaxy (M SMBH ∼ 4 · 10 6 M ⊙ ) is surrounded by a massive, very compact star cluster which is commonly referred to as the Milky Way nuclear stellar cluster (MWNSC), with structural properties that are reminiscent of a massive globular cluster. Due to its proximity, such a star cluster provides us with a wealth of data which are not available for extragalactic nuclei.
On the other hand, N -body simulations of the NSC formation from cluster infall give as result a system which is very similar to the MWNSC (similar mass, central BH mass, flattening ∼0.7, density profile and kinematic properties, see Tsatsi et al 2017). The observational characterization and dynamical interpretation of kinematical substructures can help shedding light on the discrepancy between these two models (Cole et al, 2017;Feldmeier-Krause et al, 2015).
Furthermore, the MWNSC shows a spread in metallicity distribution (e.g., see Feldmeier-Krause et al, 2017). The radial distribution of low metallicity stars is consistent with the one of higher metallicity . The observed stellar metallicity in the MWNSC indicates that only a small fraction of stars have metallicity close to the one typical of GCs (Do et al, 2015), which suggests that in-situ star formation must have occurred. On the basis of these chemical and kinematic constraints, it appears that, at least for the MWNSC, the two formation channels are not mutually exclusive.

Galactic centre objects: G1 & G2
Since its discovery (Gillessen et al, 2012), the gas and dust cloud G2 orbiting Sgr A* has caught the attention of the astronomical community. The recently recognized objects G1 and G2 are self-luminous objects with cold dust photospheres. They follow eccentric orbits around the SMBH at the Galactic Centre, and come close enough at periapse to suffer tidal removal of some of their extended material. Given that both G1 and G2 have survived their periapse passages, it has been strongly suggested that they are stellar objects (e.g., see Witzel et al, 2014) rather than compact dust clouds (e.g., see Pfuhl et al, 2015), but their nature is still heavily debated.
On one hand these objects could be the result of relatively recent stellar mergers of binary pairs (Stephan et al, 2016) whose internal orbits underwent rapid evolution toward extremely high eccentricities in response to the gravitational influence of a third body, the Galactic SMBH. The mergers should equilibrate to normal single stars on a Kelvin-Helmholtz time scale, but such induced mergers should be sufficiently common that other G-type sources are likely present. On the other hand a compact source model, where G2 is the outflow from a mass-losing star, has been presented (Ballone et al, 2017). By means of 3D adaptive mesh refinement hydrodynamical simulations with the grid code PLUTO [23] , this work provided a more detailed and realistic comparison to the observed position-velocity diagrams, showing that a slow (i.e., 50 km/s) outflow, with parameters roughly similar to those of a T Tauri star, can reproduce G2, and that in few years the central source should decouple from the previously ejected material.

Nuclear star clusters in external galaxies
Nuclear star clusters are a common feature in the majority of early-and late-type galaxies, with occupation fractions of 60 − 70% (e.g. Böker et al, 2002Böker et al, , 2004Carollo et al, 1998;Durrell, 1997). As for SMBHs, their masses correlate well with other properties of their host galaxies (such as velocity dispersion and total mass) possibly indicating a common evolutionary process for the host galaxy, the NSC, and the central BH; for this reason, often they are indistinctively referred to as Central Massive Objects (see e.g. Ferrarese et al, 2006). Georgiev et al (2016) explored scaling relations between NSC mass and host galaxy total stellar mass using a large sample of NSCs in late-and early type galaxies, including those hosting a central BH. They found differences in the slope of the correlation between NSCs in late-type and early-type galaxies among other trends depending on the galaxy type. This indicates possible different (ongoing) evolutionary processes in NSCs, depending on the host galaxy type.
9 Blue stragglers / stellar collisions and their products 9.1 Detection of blue stragglers Collisions involving unevolved main sequence (MS) stars or mass transfer in binaries with MS secondary components could lead to a rejuvenation of the resulting object (McCrea, 1964), which will appear brighter and warmer than other MS stars (blue straggler stars; BSS).
From the observational point of view, the need of conducting a systematic analysis of the BSS population of GCs in ultraviolet bands, where they emit the largest fraction of their light, is emerging. This is particularly true in dense environments, where crowding effects reduce the photometric completeness of faint sources. For this reason, the contribution of the HST UV Legacy Survey , thanks to the combined effect of high-resolution and UV sensitivity, appears as the most promising dataset to explore this population in the centre of Galactic GCs. However, future facilities providing similar performances will be needed in the near future when HST will end its mission.
While the main formation channels of BSS have been clearly identified (McCrea, 1964) and their efficiency as a function of the main cluster parameters has been both explored observationally  and interpreted in theoretical studies [23] http://plutocode.ph.unito.it/ , still other observational evidence requires a proper interpretation. In particular, the evidence of parallel BSS sequences in the colour-magnitude diagram of M30 and a few other clusters Ferraro et al, 2009) still awaits for a comprehensive explanation (i.e., can BSS formed through different formation channels contribute to produce the claimed bimodality? see e.g. Jiang et al 2017;Xin et al 2015).

Blue stragglers as dynamical clocks
BSS are more than just exotic objects in a star cluster: their unique mass range and time scale of formation make them an ideal tool to study the dynamical state of their host systems. As a result of mass segregation process, BSS are often more centrally concentrated than other GC stars. Computing the radial distribution of BSSs in clusters with different dynamical ages, Ferraro et al (2012) found that clusters can be categorised in three families based on their radial distribution profiles of BSSs. For the intermediate-age GCs (family II) the BSS distribution exhibits a minimum at position r min which is directly correlated with the clusters relaxation time scales, making BSS distributions ideal observables for the cluster's dynamical age. However, a tension between observation and theory has been put forward regarding this bimodal radial distribution which is not expected according to extensive surveys of Monte Carlo simulations (Hypki and Giersz, 2017).
Recently, Alessandrini et al (2016) and Lanzoni et al (2016) suggested A+ (i.e., the area between the cumulative radial distribution of BSS and that of a reference population) as a dynamical observable. This high-signal feature is easy to measure and does not require binning, making it a more stable observable than the radial distribution profile . A+ also shows the same strong correlation with the cluster's dynamical time scales (Ferraro et al, 2018a).

Simulations of stellar collisions
From the computational perspective, enormous progresses have been made in the last years in the simulation of close interactions in star clusters. Still, direct integration of three-four body interactions, which are crucial to study the formation of such exotic species, severely slow down performances of all the N -body and Monte Carlo codes. This can constitute a bottleneck in particular when rare events (like BBH mergers and interactions between neutron stars) require simulations of realsize (N > 10 6 ) clusters.
Furthermore, while simple recipes are already included in all simulation codes (Hurley et al, 2001;Spera et al, 2015), the complex interplay between stellar and dynamical evolution is still far from being properly accounted with a fully consistent SPH-stellar evolution approach, in particular when binary stars are considered. Nevertheless, the development of computational environments where different codes simultaneously interact, allows today to move the first steps in facing complex collisional processes where N -body, gas accretion and UV radiation must be taken into account, e.g. to model the formation of very massive stars in the first stellar complexes (Hosokawa et al, 2012) or the evolution of triples in a common Roche lobe overflow phase (de Vries et al, 2014).

Compact objects in star clusters
10.1 Binaries in star clusters X-ray sources are good tracers of compact binaries, especially low mass X-ray binaries (LMXBs), that are systems in which a compact object, i.e., a white dwarf (WD), a neutron star (NS) or a BH, accretes matter from a low-mass companion star. Accretion occurs through Roche-lobe overflow and disk formation around a compact object, or, in the case of red giants, wind-fed accretion partially captured by the compact object. NS binaries are mainly millisecond pulsars (MSPs), which form in dense clusters and are subsequently ejected as a consequence of tidal disruption and evaporation of the cluster (see e.g. Fragione et al, 2018). WD binaries, (WD-WD or WD-MS stars) can be found in nearby GCs.
In particular, WD-MS binaries are known as cataclysmic variables (CVs). The optimal way to identify possible CVs in GCs is by combining different techniques, to measure their optical variability, blue colour, Hα excess, and X-ray emission (Knigge et al, 2011). Belloni et al (2016) analysed the population of CVs in a sample of 12 GCs evolved with MOCCA, by considering two initial binary populations. They found that a population of CVs is mainly found in later stages of the evolution of GCs, and that selection effects can drastically limit the number of observable CVs. They also found that the probability of observing CVs during the outburst is extremely small, and they concluded that the best way of detecting such objects is by searching for variabilities during the quiescent phase. In addition, magnetic fields might be needed to explain the rare frequency of outbursts amongst bright CVs (∼ 10%).

Stellar-mass black holes
A relevant aspect of the dynamics of stellar-mass BHs in dense stellar environments is related to the retention of massive objects in star clusters. Indeed, the retention of BHs (and NSs) can have a strong influence on the global evolution of globular clusters (e.g., see Breen and Heggie, 2013a;Contenta et al, 2015;. Initial mass function and BH formation mechanisms (especially kicks) play a major role in determining the subsequent evolution of the BH population in a GC (Chatterjee, 2016;Mandel and de Mink, 2016), and a clear distinction is seen between the mass loss due to stellar evolution (connected with metallicity, Spera et al, 2016) and relaxation. BH subsystems can be formed and preserved in dense environments, provided the cluster has a sufficiently long relaxation time (Breen and Heggie, 2013a,b).
Massive binary black holes (BBHs) are preferentially formed in low-metallicity and dynamically active stellar environment (e.g., see Mapelli, 2016). Recent works carried out with Monte Carlo simulations pointed out that more massive clusters are more likely to trigger BBH mergers (Rodriguez et al, 2016a,b) although these events may also take place in stellar systems with a lower density (such as open clusters, see Section 11). The recent detections of gravitational waves from merging BBH have the potential to revolutionize our understanding of compact object astrophysics, but to fully utilize this new window into the universe, these observations must be compared to detailed theoretical models of BBH formation throughout cosmic time. Tanikawa (2013) and, more recently, Fujii et al (2017) calculated the detection rates of gravitational waves emitted from merging BBHs in star clusters modelled as direct N-body systems. Hurley et al (2016) have also reported quantitative confirmation of the merging of two stellar-mass BHs in a binary system which was dynamically formed in a moderately-sized direct N-body model. Merger rates determined on the basis of Monte Carlo modelling approaches have also been intensively explored (e.g., see Askar et al, 2017;Hong et al, 2018;Rodriguez et al, 2016b). More details are discussed in Section 11.

Intermediate-mass black holes
Intermediate-mass black holes (IMBHs) are defined as covering a mass range of 10 2 − 10 5 M ⊙ and have become a promising field of research. With their existence it could be possible to explain the rapid growth of SMBHs, which are observed at high redshifts (Fan, 2006), by assuming that IMBHs act as SMBH seeds (e.g., Ebisuzaki et al, 2001;Tanaka and Haiman, 2009). Recent discoveries of black holes in the centres of dwarf galaxies (Reines and Volonteri, 2015) have shown that the mass range between supermassive and stellar-mass black holes is by far not empty. However, whether or not IMBHs exist in ordinary GCs is still under debate.
Several formation scenarios have been proposed for these objects, but conclusive evidence on which one is preferred is still missing. Madau and Rees (2001) proposed that IMBHs could be the remnants of Population III stars, obtained after an evolution of a few Myr. However, if this is indeed the route to build up IMBHs, we should not expect them to be found in clusters of Population I stars. Portegies Zwart et al (2004) suggested that an IMBH could be the end product of a runaway collision in the centre of star clusters. This process needs specific initial conditions, and requires the time scale of the mass segregation of the most massive stars to be shorter than the evolution time-scale for those stars, to avoid them to evolve before they start to collide. Recently, another scenario has been proposed, also indicating that these objects form in star clusters. Giersz et al (2015) proposed that an IMBH is formed as a consequence of the build-up of BH mass due to mergers in dynamical interactions and mass transfers in binaries; in spite of the ones described before, this scenario does not require the onset of particular initial conditions, but the process of IMBH formation is highly stochastic. A larger formation probability is obtained for clusters with larger concentration, and for earlier and faster BH mass build-up. A great effort has been devoted to numerical simulations, not only to test these formation scenarios, but also to understand which properties of the host system are mainly determining the presence and the mass of an IMBH in their centre.
Direct detection of an IMBH is extremely challenging because GCs are almost gas free. Radio observations (Strader et al, 2012) of the cores of three Galactic GCs (M15, M19 and M22) do not unveil any compact source: this non-detection sets the upper limit of the masses of IMBHs in these systems to ∼ 3 − 9 × 10 2 M ⊙ . Such a result has been recently confirmed by a more extensive radio continuum survey conducted on Galactic GCs (Tremou et al, 2018). These limits suggest that either IMBHs more massive than 10 3 M ⊙ are rare in GCs, or that if they are present, they accrete in a very inefficient manner. On the other hand, an X-ray outburst has been detected in a star cluster located off-center of a large lenticular galaxy; such an event has been interpreted as providing strong evidence that the source contains an IMBH of 10 4 M ⊙ . Kains et al (2016) proposed a different method to look for IMBHs in GCs, by means of gravitational microlensing. From a suite of simulations, they estimate the probabilities of detecting such an event for Galactic GCs: as an example, they consider M22, and conclude that if it hosts an IMBH with mass 10 5 M ⊙ , there is a probability of 86% of detecting an astrometric microlensing event over a baseline of 20 yr.
Nevertheless, even though direct detections of IMBHs are challenging, signatures of the presence of an IMBH are imprinted in the phase-space distribution of stars in its immediate surrounding (Bahcall and Wolf, 1976). The effects of the presence of an IMBH on the structural and kinematic properties of the host star clusters have subsequently been explored in detail by means of direct N-body simulations (Baumgardt et al, 2004a(Baumgardt et al, ,b, 2005. In particular, two signatures are at the basis of the observational claims for the detection of IMBHs in Galactic GCs: the detection of a shallow cusp in the surface brightness profile (e.g., see Noyola et al, 2008) and the presence of a rise in the projected velocity dispersion profile towards the centre (e.g., see Anderson and van der Marel, 2010;Lützgendorf et al, 2013b). However, these signatures can also be due to other processes: core collapse, mass segregation, or a population of binaries in the centre also produce a shallow cusp in the brightness profile, as shown with dedicated N -body simulations (Vesperini and Trenti, 2010), and the central rise in the velocity dispersion profile is also not unique (Zocchi et al, 2017(Zocchi et al, , 2018. The controversial case of ω Cen is an example of this degeneracy: isotropic spherical models only reproduce the observed rise in the projected velocity dispersion when a central IMBH of mass ∼ 4 × 10 4 M ⊙ is included (Noyola et al, 2008), while by comparing anisotropic models (with radial anisotropy in the core and tangential anisotropy in the outer parts) to proper motion measurements an upper limit to the IMBH mass of only ∼ 7 × 10 3 M ⊙ is obtained (van der Marel and . This degeneracy in the signatures of the presence of an IMBH has also been explored by means of numerical simulations. Lützgendorf et al (2013a) presented a set of direct N-body simulations of GCs in an external tidal field, considering several values for IMBH masses, BHs retention fractions, and binary fractions. Their results show that the presence of an IMBH, or of a central population of binaries or stellarmass BHs increases the escape rate of high-mass stars; these simulations show a good agreement with observational mass functions and structural parameters of GCs. A similar result has been found by Arca-Sedda (2016), who proposed an analysis of numerical simulations showing that the excess of mass in the centre of a cluster could be due to the presence of a subsystem of heavy remnants orbitally segregated, instead of being due to an IMBH (see also Zocchi et al 2018). Finally, the coexistence of an IMBH and of a population of stellar-mass BHs has been explored by Leigh et al (2014), by means of direct N-body simulations.
In addition, some controversy has recently arisen when comparing data obtained by means of integrated light spectroscopy to measurements of line-of-sight velocities of single stars (see for example the emblematic case of NGC 6388, Lanzoni et al, 2013;Lützgendorf et al, 2011Lützgendorf et al, , 2015. Studies of mock observations of numerical simulations of star clusters with and without a central IMBH have been carried out to determine the magnitude of this effect. Bianchini et al (2015) showed that luminosity-weighted IFU observations can be strongly biased by a few bright stars introducing a scatter in the measurement of the velocity dispersion up to ∼ 40% around the expected value: this prevents any sound assessment of the central kinematics, and does not allow for an interpretation of the data in terms of the presence of a central IMBH. de Vita et al (2017) estimated that in 20% of the cases the analysis of data did not allow for a statistically significant detection of IMBHs having mass equal to 3% of the mass of their host, because of shot noise due to bright stars close to the IMBH, and that when considering a smaller fractional mass for the IMBH (∼ 0.1) the rate of non-detections corresponds to 75% of the cases. The combination of data from new-generation facilities (such as Gaia proper motions and MUSE spectroscopic data of the central regions of clusters) will provide further constraints to the GC dynamics, allowing us to determine if IMBHs are hiding in the cores of GCs, and how massive they are.
Leaving the uncertainties and difficulties of detecting IMBHs in GCs aside, other systems have proven to host IMBHs at their centre. In the last decade, detections of AGN in nearby dwarf galaxies have provided great evidence of the existence of SMBHs in the lower mass regime (e.g. Baldassare et al, 2015;Barth et al, 2004;den Brok et al, 2015;Filippenko and Ho, 2003;Reines et al, 2011). Inspired by those recent discoveries, researchers have branched out to search for BHs in smaller stellar systems such as ultra compact dwarf galaxies (Seth et al, 2014) and even our closest neighbour, the Large Magellanic Cloud (Boyce et al, 2017).
In addition, there have been recent suggestions of the presence of candidate IMBHs embedded in gas clouds in the region of our Galactic Centre, as based on emissions in the millimetre (Oka et al, 2016(Oka et al, , 2017 and infra-red bands (Tsuboi et al, 2017). More generally in this context, a powerful diagnostic tool is offered by the phenomenology of tidal disruption events. Indeed, there have been recent reports of events that may involve IMBHs (Kuin et al, 2018;Lin et al, 2018), with fresh developments also on the theoretical side (e.g., see Anninos et al, 2018;Rosswog et al, 2009;Tanikawa, 2018).
Over the next years, this spectrum of studies will greatly contribute in filling the lower mass range of the BH mass/host galaxy property relations and enhance our understanding of BH evolution and occupation fractions.

Astrophysical sources of gravitational waves
Currently, the detection by LIGO of GWs emitted by systems identified as coalescing binary black holes (BBH) suggests that these systems are primary sources of GW emission (Abbott et al, 2016b(Abbott et al, , 2017a. This long-awaited empirical evidence has immediately sparked a debate of great theoretical importance to understand the origins of these BBH systems and to predict their frequencies over cosmic time. So far, two main formation channels have been proposed: dynamical interactions in dense stellar environment (Portegies Zwart and McMillan, 2000) and evolution of isolated stellar binaries (Belczynski et al, 2002). Stellar environments may also play an important role in shaping the properties of the merger. Using three-body scattering experiments to model the scattering process between stars and BBHs, Rasskazov and Merritt (2017) suggest that galactic nuclei may harbour supermassive BBHs. Systems in that mass class are expected to be primary sources for GW emission, especially in the detection range relevant for the upcoming LISA observatory. Similarly, Antonini and Rasio (2016) and Bartos et al (2017) suggest that the densest population of stellar-mass BHs is expected to be found in galactic nuclei. A fraction (∼ 30%) of these BHs can reside in binaries.
An extensive survey of simulations, carried out with the MOCCA code (see Section 3), followed the long-term evolution of GCs, assessed the retention of their stellarmass BHs, and allowed to estimate the local merger rate densities for the BBHs originating from GCs Hong et al, 2018). Detection rates of BBH mergers have been determined also on the basis of direct N-body models (Fujii et al, 2017;Tanikawa, 2013). In other models with an IMBH, two-body collisions between an IMBH and a stellar mass BH are very frequent (Leigh et al, 2014;Morawski et al, 2018).
Although open clusters and young massive clusters, in general, produce less binary BH mergers per cluster than the GCs, they are more numerous, which makes them a competitive formation environment (Banerjee, 2017a). By means of N -body simulations carried out with the code NBODY7, Banerjee (2017b) studied the evolution of BBH populations in young massive and starburst clusters estimating the merger rate for each model. Their overall contribution to the observable BBH inspiral rate in the universe could be at least comparable to that from classical GCs, thereby potentially favouring the BBH detection rate from the dynamical channel. In addition, this study shows a tendency of the BBHs to coalesce within the clusters. In particular, the general-relativistic BBH mergers continue to be mostly mediated by triples that are bound to the clusters, rather than happen among the ejected BBHs. In fact, the number of such in-situ BBH mergers, per cluster, tends to increase with the introduction of a small population of primordial stellar binaries.
These BBH mergers are well consistent with the LIGO detection window and suggest a full-sensitivity LIGO detection rate of up to hundreds of BBH mergers from stellar clusters, per year. Rastello et al. (in prep.) also studied, by means of N -body simulations carried out with NBODY7, the evolution of massive BBHs (of the class expected to be detectable by LIGO) in open clusters (10 2 − 10 3 M ⊙ ), deriving a merger rate of the order of 2.1 yr −1 Gpc −3 . Most importantly, it has been recently shown that general relativistic corrections are crucial to appropriately assess the rates of mergers of BBH on eccentric orbits (Samsing, 2018;Samsing and Ramirez-Ruiz, 2017).
Concerning the origin of other specific LIGO events, Belczynski et al (2017) argues that GW170104 (Abbott et al, 2017a) may be formed via classical isolated binary evolution. The DRAGON simulations (Wang et al, 2016), which follow the dynamical evolution of a globular cluster (N = 10 6 ) using state-of-the-art highperformance GPUs and post-Newtonian treatment for BH binaries, found black hole mergers producing waveforms similar to GW150914 (Abbott et al, 2016b). The results of Park et al (2017), who also performed direct N -body simulations, also suggest that GW150914 is the consequence of a dynamically formed BBH. We are, therefore, just scratching the surface of the theoretical implications of the first LIGO detections.

Conclusions
With their elegant analogy to the classical gravitational N -body problem and tantalizing access to their resolved stellar populations, star clusters have been a prime target for both theoretical and observational astronomers. Resting on such a distinguished legacy, they have often been considered as a "solved problem" in astronomical research, but, in fact, a series of recent discoveries about their chemical, kinematic and structural properties have revolutionised our traditional picture of star clusters, making them baffling chemodynamical puzzles, prolific black hole cradles, precious galactic beacons, and even unexpected tools for galaxy evolution and near-field cosmology. Now, standing at the crossroads of Gaia's blooming era of "precision astrometry" and LIGO's revolutionary beginnings of "gravitational wave astronomy", the rich internal dynamics of star clusters truly brings them back to the centre stage, and the surprises are far from over.
We leave to the interested reader to compare our current understanding of these topics with the open questions (and their envisaged timeline) identified more than a decade ago by Davies et al (2006) and to assess the progress made in each frontier area. Many milestones have been achieved, especially on the computational side, thanks to the development of numerical algorithms which enabled us to finally "solve" the gravitational N -body problem Heggie, 2014;Nitadori and Aarseth, 2012;Wang et al, 2016), at least in the sense of following the entire evolution of models of individual clusters, with a realistic number of particles and fundamental internal and external effects. One of the original goals of providing a numerical framework for realistic simulations of the many aspects of the evolution of dense stellar systems [24] is about to be fulfilled (Portegies Zwart and McMillan, 2018). But many new challenges have also emerged, mostly driven by the outstanding progress on the observational side. We do not dare providing a comprehensive list, but we nonetheless wish to conclude by highlighting some pressing questions which have emerged during the 17th annual meeting of the MODEST community: • How can we leverage the emerging phase space richness of old globular clusters to decipher the dynamical signatures of their formation and evolutionary processes?
• Which is the best numerical approach to efficiently, yet accurately, attack the "post-million body problem"?
• How can we explore, from first principles, the multi-scale, multi-physics problem of the formation of globular clusters in a cosmological context?
• What is the role of star cluster environments in shaping the dynamical evolution of planetary systems?
• How can we bridge the gap between the current knowledge of star formation processes, the origin of the first stellar aggregates, and their subsequent longterm evolution?
• How can we exploit at best the synergy between Gaia, LIGO, JWST and other upcoming facilities to finally unveil when, where, and how globular clusters formed?
• What is the nature of the complex nexus between globular and nuclear star clusters and their formation channels?
• How can we observationally characterise and dynamically interpret the role of dark remnants in the evolution of collisional systems?
• Do globular clusters really harbour intermediate-mass black holes, and, in such a case, how did they form?
• What is the origin of the stellar-mass black holes binary revealed by the first LIGO events?
The MODEST-17 conference has provided evidence that the community invested in "Modelling and Observing DEnse STellar systems" is creative, active, collaborative and diverse. Since its early years, such a group of researchers has much evolved, in size and scope, but the overarching goal of continuing to explore the richness of collisional stellar dynamics with a multi-disciplinary approach is more timely than ever. As the broader astrophysical community is progressively appreciating the many roles played by star clusters throughout cosmic time (e.g., as contributors to the reionisation epoch, tracers of galaxy formation, cradles of gravitational wave sources), an additional challenge for the MODEST collaboration will be to find new synergies with appropriate theoretical and observational research communities, in order to formulate a comprehensive, modern view of these fascinating stellar systems in the evolving landscape of contemporary astrophysics.

Declarations
Availability of data and material Data sharing not applicable to this article as no datasets were generated or analysed during the current study.

Competing interests
The authors declare that they have no competing interests.
Funding ALV acknowledges support from a Marie Sklodowska-Curie Fellowship (MSCA-IF-EFRI NESSY 658088) and the Institute for Astronomy at the University of Edinburgh; MXC from the Poles Programme (initiated by the Belgian Science Policy Office, IAP P7/08 CHARM) and the European Union's Horizon 2020 research and innovation programme under grant agreement No. 671564 (COMPAT project); FD from the DFG Priority Programme 1573 "The physics of the interstellar medium"; VP from Charles University, grants SVV-260441 and GAUK-186216; SR from Sapienza, University of Rome, grant AR11715C7F89F177; AZ from Royal Society (Newton International Fellowship Follow-up Funding). These funding bodies had no direct role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Author's contributions All authors are equal contributors, hence are listed in alphabetical order; ALV coordinated the manuscript writing.