# A rich population of free-floating planets in the Upper Scorpius young stellar association

Núria Miret-Roig, Hervé Bouy, Sean N. Raymond, Motohide Tamura, Emmanuel Bertin, David Barrado, Javier Olivares, Phillip A. B. Galli, Jean-Charles Cuillandre, Luis Manuel Sarro, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Núria Miret-Roig, Hervé Bouy, Sean N. Raymond, Motohide Tamura, Emmanuel Bertin, et al. A rich population of free-floating planets in the Upper Scorpius young stellar association. Nature Astronomy, 2022, 6, pp.89-97. 10.1038/s41550-021-01513-x. insu-03678906

HAL Id: insu-03678906

https://insu.hal.science/insu-03678906v1

Submitted on 6 Oct 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# A rich population of free-floating planets in the Upper Scorpius young stellar association

Núria Miret-Roig<sup>1,2✉</sup>, Hervé Bouy<sup>1</sup>, Sean N. Raymond<sup>1</sup>, Motohide Tamura<sup>1,3,4,5</sup>, Emmanuel Bertin<sup>6,7,8</sup>, David Barrado<sup>9</sup>, Javier Olivares<sup>1</sup>, Phillip A. B. Galli<sup>1</sup>, Jean-Charles Cuillandre<sup>10</sup>, Luis Manuel Sarro<sup>11</sup>, Angel Berihuete<sup>12</sup> and Nuria Huélamo<sup>9</sup>

The nature and origin of free-floating planets (FFPs) are still largely unconstrained because of a lack of large homogeneous samples to enable a statistical analysis of their properties. So far, most FFPs have been discovered using indirect methods; microlensing surveys have proved particularly successful to detect these objects down to a few Earth masses<sup>1,2</sup>. However, the ephemeral nature of microlensing events prevents any follow-up observations and individual characterization. Several studies have identified FFPs in young stellar clusters<sup>3,4</sup> and the Galactic field<sup>5</sup> but their samples are small or heterogeneous in age and origin. Here we report the discovery of between 70 and 170 FFPs (depending on the assumed age) in the region encompassing Upper Scorpius and Ophiuchus, the closest young OB association to the Sun. We found an excess of FFPs by a factor of up to seven compared with core-collapse model predictions<sup>6-8</sup>, demonstrating that other formation mechanisms may be at work. We estimate that ejection from planetary systems might have a contribution comparable to that of core collapse in the formation of FFPs. Therefore, ejections due to dynamical instabilities in giant exoplanet systems must be frequent within the first 10 Myr of a system's life.

So far, most exoplanets have been detected through radial velocity or photometric modulations induced in their host stars<sup>9–13</sup>. As such, the vast majority of known exoplanets are gravitationally bound to stars. However, several free-floating planets (FFPs) have been discovered over the past two decades in astrophotometric surveys of nearby star-forming regions<sup>3,4,14–18</sup>, young associations<sup>19–22</sup>, the solar neighbourhood<sup>23,24</sup> and in gravitational microlensing surveys of the Galactic field<sup>5</sup>. These ultra-faint objects are incapable of sustaining nuclear fusion and steadily fade through time, making them easier to observe when they are very young. FFPs are compact objects of less than about 13 Jupiter masses ($M_{Jup}$) that are not bound to a star or brown dwarf, but instead wander among them. At present, four scenarios have been proposed for the formation of these extreme objects: (1) a scaled-down version of star formation via core collapse<sup>25,26</sup>; (2) within a protoplanetary disk, either like gas-giant planets through core accretion<sup>27</sup> or like companions through gravitational fragmentation of massive extended disks<sup>28,29</sup>, followed by ejection by dynamical scattering between planets in both cases<sup>30</sup>; (3) as aborted stellar embryos ejected from a stellar

nursery before the hydrostatic cores could build up enough mass to become a star<sup>31</sup>; and (4) through the photo-erosion of a prestellar core by stellar winds from a nearby OB star<sup>32</sup>. Although direct observational evidence confirms that these different processes are all at work<sup>33,34</sup>, we still do not understand their relative contributions to the overall FFP population.

Here we present a search for FFPs in the 171 deg² region occupied mainly by Upper Scorpius (USC) and Ophiuchus (Oph). We selected an elliptical area centred on RA = 243.5°, dec. = -23.1° with a semi-major axis of 8.5° in RA and a semi-minor axis of 6.4° in dec. where the spatial and temporal coverage of the observations is the best. This large complex represents the perfect hunting ground to search for young and nearly coeval FFPs owing to its proximity (120–145 pc) and youth (1–10 Myr; refs.<sup>35–38</sup>). We combined our ground-based observations in the optical and infrared with wide-field images available in various public archives (Methods). We processed and analysed a total of 80,818 individual wide-field images acquired with 18 different cameras over the past 20 years to obtain a final catalogue, the Dynamical Analysis of Nearby Clusters (DANCe)<sup>39</sup> catalogue, containing positions, proper motions and multiwavelength photometry (grizyJHK) for more than 26 million objects. We complemented the DANCe catalogue with the astrometry and photometry of the Gaia Data Release 2 (Gaia DR2)<sup>40</sup> and Hipparcos<sup>41</sup> catalogues. We used this dataset to compute probabilities of membership to USC and Oph using a probabilistic model of the distribution of the observable quantities in both the cluster and background populations (Methods). To identify the FFPs, we used the parameter space that contained the greatest number of parameters and sources with complete information, namely proper motions and iJHK photometry.

We identified 3,455 high-probability candidate members in the area covered by our study, including between 70 and 170 FFPs depending on the age assumed for the region. Approximately 20% of the members were not identified in previous studies<sup>17,42–46</sup>, and this proportion increases to 75% in the planetary-mass domain. It represents one of the largest and most homogeneous sample of FFPs identified using direct images, and constitutes an excellent benchmark for testing star and planet formation theories. Figure 1 shows the area covered by our survey in the optical and radio wavelengths. The members are distributed all over the region. The majority of the members (and FFPs) lie in the area of USC, where the extinction by

<footnote n="1">¹Laboratoire d'Astrophysique de Bordeaux, Université de Bordeaux, CNRS, Pessac, France. ²Department of Astrophysics, University of Vienna, Wien, Austria. ³The University of Tokyo, Tokyo, Japan. ⁴Astrobiology Center, National Institutes of Natural Sciences, Tokyo, Japan. ⁵NAOJ, Tokyo, Japan. ⁶CNRS, UMR 7095, Institut d'Astrophysique de Paris, Paris, France. ⁷Sorbonne Université, Institut d'Astrophysique de Paris, Paris, France. ⁸Canada-France-Hawaii Telescope - CNRS UAR2208, Kamuela, HI, USA. ⁹Centro de Astrobiología (CSIC-INTA), Departamento de Astrofísica, ESAC, Madrid, Spain. ¹⁰AIM, CEA, CNRS, Université Paris-Saclay, Université de Paris, Gif-sur-Yvette, France. ¹¹Departamento de Inteligencia Artificial, UNED, Madrid, Spain. ¹²Departamento Estadística e Investigación Operativa, Universidad de Cádiz, Puerto Real, Spain. ¹³e-mail: nuria.miret.roig@univie.ac.at</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1** | Sky distribution of stars, brown dwarfs and FFPs discovered in this study and classified assuming an age of 5 Myr. Stars are shown as yellow triangles, brown dwarfs as blue squares and FFPs as red dots. The dashed ellipses indicate the area analysed with the DANCe catalogue (Methods). The distributions are overlaid on background images taken in optical (left) and radio (right, at 857 GHz) wavelengths. Background images adapted with permission from: left, Mario Cogo (http://galaxlux.com); right, ref. ⁹⁶, EDP Sciences.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 2** | Colour-magnitude diagrams of the members of USC and Oph identified in this work. The left panel indicates the (i, i-Ks) colour-magnitude diagram and the right panel the (J, J-Ks) colour-magnitude diagram. Previously known members are shown as grey symbols and new members as black symbols. Error bars represent the uncertainty in the photometry reported in the Gaia DR2 and DANCe catalogues. The isochrones from the BHAC15⁴⁷ (solid lines) and the PARSEC-COLIBRI⁴⁸ models (dashed lines) at 3 Myr (red) and 10 Myr (blue) (as well as the extinction vector) are superimposed.</figcaption></figure>

the interstellar medium is lower (total extinction in the photometric V (550 nm) band, $A_V < 1$ mag).

We used our comprehensive membership analysis to study the origin of the FFP population in USC and Oph. We computed the number density of members as a function of their masses. This ‘mass function’ of the region constitutes a fundamental constraint for formation theories because different mechanisms predict different relative abundances of stellar, substellar and planetary-mass objects. The transformation from observed luminosities to masses requires the age to be known, and evolutionary models must be used. The age spread of the region³⁵⁻³⁸, as well as the complex and overlapping spatial and kinematic distributions of the different coexisting populations, prevented us from disentangling the various groups

and assigning ages to individual objects. Instead, we chose to assign ages of 3, 5 and 10 Myr for the entire sample and assumed that the underlying real mass distribution must be included between these borderline cases. We used the entire spectral energy distribution available for each object and the theoretical evolutionary models⁴⁷,⁴⁸ to infer the individual mass of each member (Methods). Figure 2 shows colour-magnitude diagrams of the members in two different spaces. We overlaid the theoretical isochrones at the extreme ages mentioned above to illustrate the uncertainties in the mass inherited from the lack of precise individual ages.

The *J* apparent magnitude distribution (top panel of Fig. 3) is a direct product of the observations, and thus it is not affected by the uncertainties and errors in the transformation from luminosities
<figure><img src="image_6.png" /><figcaption>**Fig. 3** | J apparent magnitude distribution and mass functions of the members of USC and Oph. Top: apparent magnitude distribution. Middle: observational mass function. Bottom: observational mass function (the mass functions from simulations$^{7,8}$ are overlaid). The shaded regions indicate the 1 and 3σ uncertainties from a bootstrap (top) and the dispersion due to the age (3–10 Myr, middle and bottom). All of the mass functions are normalized in the mass range 0.004–10 $M_⊙$. The hydrogen (75 $M_{\text{Jup}}$) and deuterium (13 $M_{\text{Jup}}$) burning limits, which indicate the boundary between stars and brown dwarfs (BD) and BD and FFPs respectively, are indicated by the vertical dotted lines according to the BHAC15 evolutionary models$^{47}$ and assuming an age of 5 Myr.</figcaption></figure>

to masses. Owing to improved statistics achieved with our larger sample, we were able to unambiguously identify a dip at planetary masses (apparent magnitude $J \approx 17.8$ mag, corresponding to absolute magnitude $M_J \approx 12$ mag and masses of 7–13 $M_{\text{Jup}}$ depending on the age assumed) that was also reported in volume-limited samples in the solar neighbourhood$^{49,50}$. Are we seeing changes relative to the formation mechanisms? Could this be the ‘real’ frontier between brown dwarfs and FFPs? The origin of this dip remains unknown and more studies are needed to answer these questions. However, the presence of this feature at the young ages of USC and Oph (1–10 Myr) suggests that it must be the result of the formation and/or early evolution of these objects.

Figure 3 (middle panel) shows the mass function at 5 Myr (which can be considered an intermediate age of the entire complex); the mass functions at 3 and 10 Myr define the upper and lower confidence intervals. A number of details and features that were clearly statistically significant in the magnitude distribution (in particular the dip mentioned above) are blurred and lost in the mass distribution, probably due to the uncertainties related to the transformation from luminosity to mass. Our magnitude and mass distributions show a rich population of FFPs. We evaluated the contamination rate in this sample, mostly due to background reddened

giant stars and background Galaxies, using two different methods. First, using synthetic data$^{51}$, we estimated a contamination rate in the planetary-mass domain smaller than 4%. Second, using sources for which both DANCe and Gaia DR2 measurements are available (in the range $9 < J < 14$ mag, hence above the planetary-mass regime) and assuming that the Gaia DR2 sample represents the ground truth, we estimated a contamination rate of approximately 8%. The real contamination rate of FFPs is probably between these two values and should be confirmed by follow-up spectroscopic observations. The number of FFPs reported in our analysis must constitute a lower limit of the true total number of FFPs as our analysis is expected to miss the objects most affected by extinction ($A_V \gtrsim 3$ mag), as well as objects displaying a large near-infrared excess related to the presence of circumstellar material.

The fraction of FFPs in our sample, meaning the relative proportion of FFPs to stars and brown dwarfs, is $0.045^{+0.023}_{-0.029}$ where the uncertainties come from the uncertainty in individual ages (between 3 and 10 Myr). We estimated this fraction by first integrating the observational mass function (middle panel of Fig. 3) in the FFP (4–13 $M_{\text{Jup}}$), brown dwarf (13–75 $M_{\text{Jup}}$) and stellar (0.075–10 $M_⊙$, where $M_⊙$ is the mass of the Sun) mass regimes and then computing the ratio between FFPs and brown dwarfs plus stars. Whereas previous studies reported a similar fraction of FFPs$^{3,4}$, our sample doubles the absolute number of FFPs in a single association, which significantly reduces the statistical uncertainties in the mass function. A recent photometric study in the central 6 deg$^2$ region of USC$^{46}$ found a ratio of 1.0–1.5 planetary-mass members per square degree; this result is broadly consistent within the statistical uncertainties with the ratio we find in the 171 deg$^2$ area covered by our survey (0.5–1.0 planetary-mass members per square degree).

In the high-mass regime ($>1 M_⊙$), our mass function has a slope of $\Gamma = -1.2 \pm 0.2$ (measured in units of logarithmic mass), compatible with the Salpeter slope$^{52}$ and with all the models and simulations considered here$^{6–8,53,54}$. In the substellar mass regime ($<75 M_{\text{Jup}}$), our observational mass function (bottom panel of Fig. 3) has a slope of $\Gamma = 0.62^{+0.13}_{-0.01}$ (measured in units of logarithmic mass), similar to values reported in the field population for L, T and Y dwarfs$^{23,24}$ and to models that include several channels of substellar object formation$^{53,54}$. The mass function over the low-mass stars and high-mass brown dwarfs regime ($0.03–1 M_⊙$) is compatible with a log-normal distribution$^{6,7}$. When integrating the analytical mass function of models that include mostly core-collapse formation$^{6,7}$ over the planetary-mass range (4–13 $M_{\text{Jup}}$), we find that they predict a fraction of only 0.009–0.019 FFPs, underestimating our measurement by up to seven times (depending on the age assumed). This excess of FFPs with respect to a log-normal mass distribution is in good agreement with the results reported in $\sigma$ Orionis$^4$. Interestingly, our observational mass function also has an excess of low-mass brown dwarfs and FFPs with respect to simulations that include both core collapse and disk fragmentation$^8$. This suggests that some of the FFPs in our sample could have formed via fast core accretion in disks, rather than by disk fragmentation. We also note that the continuity of the shape of the mass function at the brown dwarf/ planetary mass transition suggests a continuity in the formation mechanisms at work for these two classes of objects.

Hereafter, we use the current knowledge of planet and star formation to discuss the origin of FFPs. The fraction of observed FFPs ($f_{\text{FFP observed}}$) is the sum of FFPs formed by ejection from a disk ($f_{\text{FFP ejected}}$), the FFPs formed by core collapse ($f_{\text{FFP core collapse}}$) and the FFPs formed by other mechanisms ($f_{\text{FFP other}}$, including photo-evaporation and ejection from a prestellar cluster). $f_{\text{FFP ejected}}$ depends on the fraction of stars and brown dwarfs that form giant planets ($f_{\text{giant}}$), on the fraction of such planetary systems that become unstable ($f_{\text{unstable}}$) and on the number of ejected planets per unstable system ($n_{\text{ejected}}$). This formulation is similar to that of previous studies$^{30}$ and has the following expression in our study, for objects more massive than 4 $M_{\text{Jup}}$,
<table><caption>**Table 1 | Planetary occurrence rate for different studies**</caption><thead><tr><th>Planetary detection method</th><th>Planet<br>occurrence rate<br>(%)</th><th>Mass range of the<br>planets (M<sub>jup</sub>)</th><th>Separation or<br>semi-major axis</th><th>Spectral type of<br>the primary body</th><th>f<sub>g</sub> {FFP<br>ejected}<sup>c</sup></th><th>f<sub>g</sub> {FFP ejected<br>FFP observed<br>(%)</th></tr></thead><tbody><tr><td>Radial velocity<sup>12,55</sup></td><td>2.1 ± 0.5</td><td>4-13</td><td>0.1-100 au</td><td>Solar</td><td>0.016-0.021</td><td>23-130</td></tr><tr><td>Direct imaging<sup>57</sup></td><td>0.6<sup>+0.7</sup><sub>-0.5</sub></td><td>5-13</td><td>30-300 au</td><td>BAFGKM</td><td>0.005-0.006</td><td>7-38</td></tr><tr><td>Direct imaging<sup>94</sup></td><td>1.83<sup>+5.76</sup><sub>-0.62</sub></td><td>5-13</td><td>30-300 au</td><td>BAFGKM</td><td>0.014-0.018</td><td>20-115</td></tr><tr><td>Microlensing<sup>58</sup></td><td>1-2</td><td>0.007-0.02<sup>a</sup></td><td>0.2-5 θ<sub>E</sub></td><td>all</td><td>0.008-0.02</td><td>11-125</td></tr><tr><td>Combination (radial<br>velocity + direct<br>imaging + microlensing)<sup>95</sup></td><td>-4</td><td>1-13<sup>b</sup></td><td>2-1,000 au</td><td>all</td><td>0.03-0.04</td><td>44-250</td></tr></tbody></table>

The conservative ranges we provide include both the uncertainties from our observations and current knowledge of planetary systems. ΘE = angular Einstein radius. <sup>a</sup>Mass ratio between the planet and star mass (q). <sup>b</sup>Our observations are not sensitive to masses between 1 and 4 M<sub>jup</sub> and thus, the planet occurrence rate is overestimated with respect to our observations. <sup>c</sup>Assuming f<sub>unstable</sub> = 0.75-1 and n<sub>ejected</sub> = 1.

$$ 
\begin{aligned}
f_{\text{FFP observed}} &= f_{\text{FFP ejected}} + f_{\text{FFP core collapse}} + f_{\text{FFP other}} \\
&= f_{\text{giant}} \times f_{\text{unstable}} \times n_{\text{ejected}} \\
&\quad + \int_{4 M_{\text{jup}}}^{13 M_{\text{jup}}} \xi_{\text{core collapse}}(m) \, dm + f_{\text{FFP other}}
\end{aligned} 
\quad (1) 
$$

The value of $f_{\text{giant}}$ is constrained by the observed demographics of giant exoplanets measured by radial velocity, direct imaging, transit and microlensing surveys$^{11,12,55-62}$. The planetary occurrence rate depends on many astrophysical parameters (the host star and planet masses, orbital separation, stellar metallicity and others), and each of the techniques above mentioned is sensitive to a specific region of the parameter space$^{63}$. Therefore, we combined the occurrence rates obtained by different techniques to minimize the possible observational biases. In Table 1, we summarize the occurrence rates obtained by different authors and describe the properties of each study. $f_{\text{unstable}}$ must be at least 75% to match the observed distribution of giant planet eccentricities$^{64-68}$. A minimum of two planets per system is needed for instability to occur and many giant exoplanets are in fact found in multiple planetary systems or contain hints (such as radial velocity trends) of additional companions$^{61}$. $n_{\text{ejected}}$ scales with the number of planets involved in the instability$^{60}$. We consider a simplified scenario in which planetary systems contain two to four giant planets, the three cases being considered equally likely, and every time a system becomes unstable it ejects one planet. With these assumptions, we find $n_{\text{ejected}} \approx 1/3 \times 2 \times 1/2 + 1/3 \times 3 \times 1/3 + 1/3 \times 4 \times 1/4 \approx 1$.

Combining the upper and lower limits on the three factors defining the fraction of FFPs ejected from planetary systems, we obtain $f_{\text{FFP ejected}} \approx 0.005-0.021$, which represents between 10 and 130% of the FFP population we found. The large uncertainty in this percentage is due to the lack of precise masses in our observations on the one hand (which itself is related to the lack of precise individual ages), and to the uncertainty in the occurrence rate of planets and their ejection process on the other. A percentage of planets formed by ejection above 100% is obviously of no physical meaning and reflects the limitations of our assumptions and simplifications, which are based on the best current knowledge of planetary systems and the latest evolutionary models. This result nevertheless suggests that ejection from planetary systems is a significant mechanism of FFP formation as at least 10% must have formed by ejection from a disk.

The discovery of a large population of FFPs in USC and Oph also has important implications for the formation and early evolution of planetary systems and, specifically, for the timescale of the processes involved. N-body simulations indeed predict that most dynamical instabilities happen within 0.1–1 Myr of the planets’ formation, although there are configurations that produce later

instabilities$^{64-66,69}$. Our results suggest that giant planet systems must form and become dynamically unstable within an observed lifetime of the region of 3–10 Myr to contribute to the population of FFPs. Although the instability among the Solar System’s giant planets$^{70}$ was much less violent than those for the mass range of FFPs in our sample$^{71}$, current studies suggest that it may have also happened early$^{72,73}$.

Instabilities can also result from close stellar encounters: numerical simulations have shown that dynamical interactions with other stars in clustered environments may induce instability in planetary systems or even liberate planets, enriching the FFP population$^{69,74}$. Recent studies showing that the demographics of exoplanets depend on the stellar environment$^{75}$ confirm that such interactions must indeed play a role and contribute to the observed population of FFPs. Our observations suggest that these encounters might take place within the first 10 Myr of a system’s life.

The combined contributions of FFPs from core collapse (13–118%) and ejection from planetary systems (10–130%) derived from our analysis can explain the formation of the majority of FFPs. But other mechanisms are known to be at work: photo-erosion of prestellar cores$^{32,76}$ has been observed around massive B stars$^{77,78}$. Since USC is an OB association, some of the FFPs might have formed by photo-erosion. At the same time, given that photo-erosion can only occur in the direct vicinity of relatively scarce OB stars, we can reasonably assume that the contribution of this mechanism to the overall FFP population must be fairly small. The contribution of dynamical ejections of substellar embryos in the parent prestellar cluster$^{29,31}$ could, in principle, be a significant source of FFPs. However, recent hydrodynamical simulations show that while such dynamical ejections can produce a realistic population of brown dwarfs, they under-produce FFPs$^{8}$ (Fig. 3). This new sample of FFPs is larger and more comprehensive than previous samples and offers exciting opportunities to better understand their origin by performing statistically robust studies of their properties. The multiplicity, kinematics and properties of disks among FFPs are expected to depend on the processes at work and should hold important clues to their formation and early evolution. Finally, this new sample also provides an opportunity to study the atmospheres of planetary-mass objects in the absence of a blinding host star.

## Methods

**DANCe catalogue.** The COSMIC-DANCe project$^{39}$ aims to perform deep and complete censuses of young nearby associations down to the planetary-mass regime. Although the identification of young stars and even brown dwarfs has become trivial with the advent of Gaia, finding the extremely faint FFPs well beyond Gaia’s sensitivity limit remains a difficult challenge that requires deep ground-based observations. The strategy chosen for COSMIC-DANCe relies on the measurement of proper motions, multiwavelength photometry and the use of modern data-mining techniques to identify the faintest members. We therefore
combine optical and infrared wide-field images found in public archives with our own observations to obtain multi-epoch and multiwavelength catalogues of young clusters and star-forming regions in the solar neighbourhood. The relatively long time baseline encompassed by this dataset (~20 years) allowed us to derive proper motions with a precision of the order of 1 mas yr⁻¹.

We searched in the European Southern Observatory (ESO), National Optical-Infrared Astronomy Research Laboratory (NOIRLab), Palomar Transient Factory, Canadian Astronomy Data Centre, Isaac Newton Group, WFCAM Science Archive and Subaru-Mitaka-Okayama-Kiso-Archive public archives for wide-field images inside the area

$$235^{\circ} < RA < 252^{\circ}, \\ -29.5^{\circ} < \text{dec.} < -16.7^{\circ}. \qquad (2)$$

We complemented the data found in these public archives with our observations with the Dark Energy Camera (DECam) mounted on the Blanco telescope at the Cerro Tololo Inter-American Observatory, the VISTA and VST telescopes at ESO, the MegaCam camera at Canada-France-Hawaii Telescope, the NEWFIRM camera mounted on the 4m telescope at the Kitt Peak National Observatory and Cerro Tololo Inter-American Observatory, the Hyper Suprime-Cam mounted on the Subaru telescope and the Wide Field Camera mounted on the Isaac Newton Telescope. Supplementary Table 1 gives an overview of the various instruments used for this study. We used the MaxiMask and MaxiTrack software packages⁷⁹ to detect problematic pixels (cosmic rays, dead/hot pixels, satellite trails, saturated stars and associated blooming artefacts, diffraction spikes) and problematic astronomical images (for example, where tracking was lost). We also discarded a few images after a visual inspection because of their poor quality, limited sensitivity or acquisition problems. We collected a total of 80,818 individual images from 18 different instruments obtained over the past 20 years. The raw and processed data added up to almost 120 TB and were processed and analysed on a dedicated High Performance Computer (HPC) server. We extracted over 1.3 billion individual detections from these images. The image reduction and the photometric and astrometric analysis were described in a previous article³⁹. Briefly, all images were either processed using their official pipeline (for example DECam, Hyper Suprime-Cam) or using Alcamb³⁶ following standard procedures. Sources were detected, extracted and their astrometry and photometry measured using SExtractor⁸¹ and PSFEx⁸². The astrometric and photometric calibrations were then obtained using SCAMP³⁵, and nightly sets of individual images were average-combined (weighted by their exposure times) using SWarp⁸³. Sources were then extracted from these deeper nightly stacks, and proper motions and photometry measured using SCAMP again. Deep stacks were also produced, combining all available images obtained on a given camera plus filter and used for the photometry only. Details of the whole procedure can be found in the original COSMIC-DANCE article³⁹. In Supplementary Fig. 1 we show the precision of the proper motions as a function of magnitude. Because of the superiority of Gaia compared with our ground-based observations, we always used the astrometry from Gaia DR2 when available and the astrometry from DANCe elsewhere. This explains the shift in precision observed at *i* ≈ 21 mag.

Our catalogue contains proper motions and photometry (grizyJHK) for 40,882,164 unique sources. To optimize the number of sources with complete photometry (which was essential for the membership analysis, see the following section) and accelerate the computational time of the membership algorithm, we selected the area where the coverage of most instruments was best. We defined an elliptical area centred on RA = 243.5°, dec. = -23.1°) with a semi-major axis of 8.5° in RA and a semi-minor axis of 6.4° in dec. This selection roughly follows the coverage of the UKIDSS near-infrared survey⁸⁵, which we also used to define our DECam and Hyper Suprime-Cam surveys.

The globular cluster NGC 6121 (with proper motions in RA and dec. $\mu_{\text{r}} = -12.48 \text{ mas yr}^{-1}$, $\mu_{\text{s}} = -18.9 \text{ mas yr}^{-1}$; ref. ⁸⁶) is inside the area covered by our survey. To avoid contamination from its members in our sample, we discarded the sources encompassed in a circular region of 12' around the globular cluster centre (RA = 245.896°, dec. = -26.527°). The final catalogue contains 28,062,542 sources and has a median precision of <1 mas yr⁻¹ in proper motions for sources brighter than *i* < 20 mag. Supplementary Fig. 2 shows the density distribution of sources as a function of magnitude for different filters. We used the maximum of this distribution as the approximate completeness limit in each band. However, it depends on dust extinction and varies with position in the DANCe catalogue.

**Membership analysis.** We used a maximum likelihood approach to infer the parameters of the models describing both the cluster and field populations²¹,³⁷. The algorithm models first the distribution in the space of observables (parallaxes, proper motions and photometry) of the sources that belong to the field and then iteratively searches for a maximum likelihood solution for the parameters of the model that describes the distribution of observables for the cluster sources. In each iteration, the algorithm calculates the membership probabilities using Bayes' theorem and the fractions of cluster and field members as priors. The sources with missing data cannot be used to infer the cluster and field models and it is therefore extremely important to define an adequate representation space (that is, the largest possible set of parameters that contains a large fraction of sources with

complete observations). We note, however, that our algorithm uses the final model (computed with complete sources) to obtain a membership probability for the sources with partial information by marginalizing over the missing information. We searched for members in three different catalogues: the DANCe catalogue (produced in this work), the Gaia DR2 catalogue and the Hipparcos catalogue. These catalogues include very different photometric bands, and we decided to run an independent analysis for each catalogue. The parameters (photometry and astrometry) used in each case are described in the following sections.

**Initial members.** We compiled a list of 2,865 published candidate members in the literature¹⁷,⁴²,⁴³,⁴⁵,⁸⁸ in the area covered by this study. We cross-matched this list with each of our three catalogues (Hipparcos, Gaia DR2 and DANCe) to obtain the initial list to start each analysis. In the case of Hipparcos, we excluded Antares (α Sco) as it is a giant star, and therefore falls out of the empirical pre-main-sequence isochrone. For the analysis with Gaia and DANCe, we excluded the most extinguished members as they confuse our empiric isochrone.

**Representation space.** For the analysis with the Hipparcos catalogue, we searched for members in the space of pmra, pmdec, parallax, V and B-V for which all the sources in the catalogue have complete observations. For the analysis with the Gaia DR2 catalogue, we used the representation space pmra, pmdec, parallax, $G_{RP}$ and $G - G_{RP}$ excluding the $G_{RP}$ band, which is less accurate for cool dwarfs⁸⁹. In this space, 7,768,856 sources (97%) have complete observations. For the analysis with the DANCe catalogue, we used the representation space pmra, pmdec, *i*, *J*, *H*, *i* - *Ks*. We combined the *i* band in the optical, which has the largest coverage, with the infrared bands *J*, *H* and *Ks* where the ultracool dwarfs are best detected. With this representation space, 10,483,667 sources have observations in all the photometric bands, which represents 37% of the catalogue. We decided not to include the *g*, *r*, *z* and *Y* bands in the representation space because they reduced the number of sources with complete photometry, especially for the coolest objects.

**Field model.** The model of the field population is a Gaussian mixture model in the whole representation space. We explored models with different number of Gaussians and used the Bayesian information criterion (BIC) criterion to choose the final model. As the Hipparcos catalogue has a reduced number of sources, we explored models with several components between 1 and 20 and selected 6 as the optimum choice according to the BIC. For the Gaia DR2 catalogue, we explored models with between 20 and 180 Gaussians and chose 60 as the optimum choice according to the BIC. Finally, for the DANCe catalogue, we explored models with between 60 and 300 and chose 100 components.

**Cluster model.** The proper motion distribution of the region of USC and Oph is much more complex than that of open clusters previously analysed with this methodology²¹,⁸⁷,⁹⁰. Whereas open clusters have symmetric nearly Gaussian distributions in astrometry, this young region shows a rich substructure far from Gaussian (Supplementary Fig. 3) that is indicative of multiple kinematic populations. To model this complex distribution, we used a Gaussian mixture model where the Gaussians are not necessarily concentric and explored models with between 1 and 10 Gaussians. As the Hipparcos catalogue contains a very reduced number of sources, we found that a single multivariate Gaussian function sufficed to model the cluster proper motions and parallaxes. The Gaia DR2 and DANCe catalogues are much larger and the number of Gaussians selected according to the BIC criterion is 5-7 for Gaia DR2 and 5-6 for DANCe (depending on the internal probability ($p_{\text{in}}$) parameter, see below).

We ran the model with different internal probability thresholds⁷⁸, that is different degrees of completeness and contamination ($p_{\text{in}} = 0.5, 0.6, 0.7, 0.8$ and 0.9), and for each we computed the optimum threshold, $p_{\text{opt}}$ using synthetic data²¹,⁸⁷. In Supplementary Table 2, we show $p_{\text{in}}$, $p_{\text{opt}}$ and the number of members for each independent analysis (Hipparcos, Gaia DR2 and DANCe).

**Final list of members.** Membership probabilities obtained from the analysis with different $p_{\text{in}}$ values, as well as the astrometry and photometry used in the Hipparcos, Gaia DR2 and DANCe catalogues, are available at Centre de Données astronomiques de Strasbourg. Choosing the best solution (the best $p_{\text{in}}$) is a non-trivial decision that depends on the aim of the study. To study the magnitude distribution and mass function, we needed a list of members as complete as possible. For this reason, we preferred solutions with low $p_{\text{in}}$ values, which have greater completeness, although they can also be slightly more contaminated.

First, we compared the Gaia DR2 solutions obtained with different $p_{\text{in}}$ and found that 2,603 sources (94%) are the same in all the lists. The contamination rate and true positive rate computed with synthetic data are also very similar in the five studies (Supplementary Table 2), so we had no reason to prefer one list over another. Therefore, we chose the list for $p_{\text{in}} = 0.5$ as the final list from Gaia DR2 as it was the one with the largest number of members. Following an analogous procedure with the DANCe solution, we also chose the list for $p_{\text{in}} = 0.5$. With the Hipparcos study, we selected the solution for $p_{\text{in}} = 0.7$, which represents a good compromise between low contamination and high completeness. To this final list, we added the giant star Antares manually. Our final list of members contained 3,455 sources from the Hipparcos, Gaia DR2 and DANCe catalogues.
**Membership completeness.** The completeness of our membership analysis depended on the completeness of the astrophotometric catalogue and the membership algorithm. Our optical DECam and near-infrared VISTA images and archival UKIRT images covered the entire area and ensured that the instrumental i, Y, J, H and Ks sensitivities were fairly homogeneous spatially. To get an estimate of the completeness in the stellar mass regime, we propagated the apparent magnitude completeness of the DANCe catalogue to masses. The limiting magnitude to search for ultracool dwarfs was set by the i band, which we estimated to be sensitive up to $i \approx 26$ mag and complete up to $i \approx 23$ mag (Supplementary Fig. 2). This approximate magnitude limit completeness corresponded to masses between 7 $M_{\text{hup}}$ (assuming an age of 3 Myr) and 13 $M_{\text{hup}}$ (assuming an age of 10 Myr). However, our membership algorithm is expected to miss objects highly affected by extinction ($A_V \ge 3$ mag) and sources with near-infrared excess related to the presence of circumstellar material.

**Membership validation.** The membership classification mostly coincided in the Hipparcos–Gaia DR2 and Gaia DR2–DANCe studies, in the magnitude range where both catalogues are complete. The small differences between catalogues can be attributed to the different information provided by each of them. We used the Gaia DR2 membership analysis, with the additional information on the parallax, to estimate a contamination rate of 8% on the DANCe membership in the magnitude range $9 < J < 14$ mag assuming that the Gaia DR2 selection was perfectly clean. Similarly, the comparison between the Gaia DR2 and DANCe samples over the common luminosity domain shows that one-third of the objects identified with Gaia DR2 were not recovered with DANCe because of either missing photometry, high extinction or near-infrared excess probably related to the presence of a circumstellar disk. Therefore, the completeness of our census was expected to be better in USC than in Oph, as extinction is much lower and near-infrared excesses related to disks should be less frequent given the more advanced age and timescale for disk decay.

We note that at $J \approx 10$ mag (Fig. 2) there are some contaminants. These are sources identified with the DANCe membership analysis, therefore, using only the proper motions and photometry (no parallaxes). The proper motions are compatible with the proper motion distribution of the association and that is why they are classified as members. Moreover, as they are bright stars, their photometry saturates in many bands so we have little photometric information. In any case, these represent 1% of the members, which is within the contamination rate that we estimated.

We recovered the majority of members previously reported in the literature$^{17,42–45}$. Supplementary Fig. 4 shows a comparison between the members reported by previous studies and the members found in this study, using the photometry and proper motions we measured in this work; in this figure we only considered studies sensitive to FFPs. We missed around 80 stellar members previously reported in the literature, most of which are in Oph, are highly extincted or host circumstellar disks (as can be seen in the colour–magnitude diagrams). Some of the members reported in the literature were discarded by our membership analysis because of their inconsistent proper motions. This is particularly significant with respect to a recent study$^{46}$; only half of their members were identified in our astrophotometric analysis. The other half are either classified as non-members (having photometric and/or proper motion measurements inconsistent with the association) or undetected. In this study, we add ~800 of new members, 70–170 of which are FFPs, depending on the age assumed.

**The mass function. Distances.** We used Kalkayotl$^{91}$ (https://github.com/olivares-j/kalkayotl) to infer Bayesian distances for all the members with a parallax measurement in the Gaia DR2 catalogue. We used a Gaussian prior with a locus and scale of 145 pc and 45 pc, corresponding to the median and 3 s.d. of the distribution of distances obtained by inverting the parallax. For the sources in the DANCe catalogue without a parallax measurement, we sampled the distance from the cluster distance distribution obtained with all the Gaia DR2 members.

**Masses.** We combined the apparent photometry (grizyJHKs) and the distance estimate of each star to obtain absolute magnitudes. These were compared, in a Bayesian framework, to theoretical evolutionary models to infer the posterior distribution of the mass and extinction of each source with Sakam$^{21}$ (https://github.com/olivares-j/Sakam). This algorithm ignores any possible source of error related to the theoretical evolutionary model chosen by the user. The model does not include effects on the variability of the source due to binarity, activity or other factors. These effects are eventually included in the extinction estimate, enlarging its uncertainties. Finally, the mass and the extinction are degenerated, but at least half of the planetary-mass objects we found are in regions of low extinction, favouring their planetary nature.

We had to use different models for the high- and low-mass regimes as there is no single set of models covering the entire mass range of our members. We combined the PARSEC-COLIBRI (http://stev.oapd.inaf.it/cgi-bin/cmd) models$^{48}$ and BHAC15 (http://perso.ens-lyon.fr/isabelle.baraffe/) models$^{47}$, which cover the high- and low-mass range of our members, respectively. We found that both grids of models agree fairly well around $0.5 M_\odot$ and decided to use the masses inferred from the BHAC15 models below $0.5 M_\odot$ and the masses inferred with the PARSEC-COLIBRI models above $0.5 M_\odot$. In Supplementary Fig. 5 we show two

examples of the best-fit spectral energy distribution models obtained with Sakam for a brown dwarf and a planetary-mass object. The complete final list of members with the masses and extinctions inferred with Sakam is available at the Centre de Données astronomiques de Strasbourg.

**Magnitude and mass distributions.** To obtain the magnitude (mass) distribution, we sampled the individual magnitude (mass) of each source with a Gaussian centred at the measured magnitude (mass) and a standard deviation equal to the uncertainty. Then, we defined a grid between the least and most bright (massive) object in our sample and added the contribution of all the sources to each magnitude (mass) bin. We convoluted this distribution with a Gaussian kernel density estimation with a bandwidth chosen according to Scott's$^{92}$ and Silverman's$^{93}$ rules. We estimated the uncertainties in the magnitude (mass) function with a bootstrap of 100 repetitions and reported the $1\sigma$ and $3\sigma$ confidence levels.

## Data availability

The data that support the findings of this study are available at the Centre de Données astronomiques de Strasbourg at: via anonymous FTP to cdsarc.u-strasbg. fr (130.79.128.5) or via https://cdsarc.unistra.fr/viz-bin/cat/J/other/NatAs/ or from the corresponding author upon reasonable request.

Received: 6 May 2021; Accepted: 13 September 2021; Published online: 22 December 2021

## References

1. Mróz, P. et al. A terrestrial-mass rogue planet candidate detected in the shortest-timescale microlensing event. *Astrophys. J. Lett.* **903**, L11 (2020).

2. Ryu, Y.-H. et al. KMT-2017-BLG-2820 and the nature of the free-floating planet population. *Astron. J.* **161**, 126 (2021).

3. Scholz, A. et al. Substellar Objects in Nearby Young Clusters (SONYC). VI. The planetary-mass domain of NGC 1333. *Astrophys. J.* **756**, 24 (2012).

4. Peña Ramírez, K., Béjar, V. J. S., Zapatero Osorio, M. R., Petr-Gotzens, M. G. & Martin, E. L. New isolated planetary-mass objects and the stellar and substellar mass function of the $\sigma$ Orionis cluster. *Astrophys. J.* **754**, 30 (2012).

5. Mróz, P. et al. No large population of unbound or wide-orbit Jupiter-mass planets. *Nature* **548**, 183–186 (2017).

6. Chabrier, G. in *The Initial Mass Function 50 Years Later* Astrophysics and Space Science Library Vol. 327 (eds Corbelli, E. et al.) (Springer, 2005).

7. Haugbolle, T., Padoan, P., & Nordlund, Å. The stellar IMF from isothermal MHD turbulence. *Astrophys. J.* **854**, 35 (2018).

8. Bate, M. R. The statistical properties of stars and their dependence on metallicity. *Mon. Not. R. Astron. Soc.* **484**, 2341–2361 (2019).

9. Charbonneau, D., Brown, T. M., Latham, D. W. & Mayor, M. Detection of planetary transits across a Sun-like star. *Astrophys. J. Lett.* **529**, L45–L48 (2000).

10. Lissauer, J. J. et al. A closely packed system of low-mass, low-density planets transiting Kepler-11. *Nature* **470**, 53–58 (2011).

11. Howard, A. W. et al. Planet occurrence within 0.25 AU of solar-type stars from Kepler. *Astrophys. J. Suppl. Ser.* **201**, 15 (2012).

12. Mayor, M. et al. The HARPS search for southern extra-solar planets XXXIV. Occurrence, mass distribution and orbital properties of super-Earths and Neptune-mass planets. Preprint at https://arxiv.org/abs/1109.2497 (2011).

13. Howard, A. W. et al. The California Planet Survey. I. Four new giant exoplanets. *Astrophys. J.* **721**, 1467–1481 (2010).

14. Luhman, K. L., Esplin, T. L. & Loutrel, N. P. A census of young stars and brown dwarfs in IC 348 and NGC 1333. *Astrophys. J.* **827**, 52 (2016).

15. Esplin, T. L. & Luhman, K. L. A survey for planetary-mass brown dwarfs in the Taurus and Perseus star-forming regions. *Astron. J.* **154**, 134 (2017).

16. Zapatero Osorio, M. R., Béjar, V. J. S. & Peña Ramírez, K. Optical and near-infrared spectra of $\sigma$ Orionis isolated planetary-mass objects. *Astrophys. J.* **842**, 65 (2017).

17. Lodieu, N., Zapatero Osorio, M. R., Béjar, V. J. S. & Peña Ramírez, K. The optical + infrared L dwarf spectral sequence of young planetary-mass objects in the Upper Scorpius association. *Mon. Not. R. Astron. Soc.* **473**, 2020–2059 (2018).

18. Esplin, T. L. & Luhman, K. L. A survey for new members of Taurus from stellar to planetary masses. *Astron. J.* **158**, 54 (2019).

19. Liu, M. C. et al. The extremely red, young L dwarf PSO J318.5338-22.8603: a free-floating planetary-mass analog to directly imaged young gas-giant planets. *Astrophys. J. Lett.* **777**, L20 (2013).

20. Kellogg, K. et al. A targeted search for peculiarly red L and T dwarfs in SDSS, 2MASS, and WISE: discovery of a possible L7 member of the TW Hydrae association. *Astron. J.* **150**, 182 (2015).

21. Schneider, A. C., Windsor, J., Cushing, M. C., Kirkpatrick, J. D. & Wright, E. L. WISEA J114724.10-204021.3: a free-floating planetary mass member of the TW Hya association. *Astrophys. J. Lett.* **822**, L1 (2016).

22. Best, W. M. J. et al. A search for L/T transition dwarfs with Pan-STARRS1 and WISE. III. Young L dwarf discoveries and proper motion catalogs in Taurus and Scorpius-Centaurus. *Astrophys. J.* **837**, 95 (2017).
23. Kirkpatrick, J. D. et al. Preliminary trigonometric parallaxes of 184 late-T and Y dwarfs and an analysis of the field substellar mass function into the "planetary" mass regime. *Astrophys. J. Suppl. Ser.* **240**, 19 (2019).

24. Kirkpatrick, J. D. et al. The field substellar mass function based on the full-sky 20 pc census of 525 L, T, and Y dwarfs. *Astrophys. J. Suppl. Ser.* **253**, 7 (2021).

25. Padoan, P., & Nordlund, Å. The stellar initial mass function from turbulent fragmentation. *Astrophys. J.* **576**, 870–879 (2002).

26. Hennebelle, P. & Chabrier, G. Analytical theory for the initial mass function: CO clumps and prestellar cores. *Astrophys. J.* **684**, 395–410 (2008).

27. Pollack, J. B. et al. Formation of the giant planets by concurrent accretion of solids and gas. *Icarus* **124**, 62–85 (1996).

28. Boss, A. P. Formation of extrasolar giant planets: core accretion or disk instability? *Earth Moon Planets* **81**, 19–26 (1998).

29. Bate, M. R., Bonnell, I. A. & Bromm, V. The formation mechanism of brown dwarfs. *Mon. Not. R. Astron. Soc.* **332**, L65–L68 (2002).

30. Veras, D. & Raymond, S. N. Planet-planet scattering alone cannot explain the free-floating planet population. *Mon. Not. R. Astron. Soc.* **421**, L117–L121 (2012).

31. Reiputh, B. & Clarke, C. The formation of brown dwarfs as ejected stellar embryos. *Astron. J.* **122**, 432–439 (2001).

32. Whitworth, A. P. & Zinnecker, H. The formation of free-floating brown dwarves and planetary-mass objects by photo-erosion of prestellar cores. *Astron. Astrophys.* **427**, 299–306 (2004).

33. Testi, L. et al. Brown dwarf disks with ALMA: evidence for truncated dust disks in Ophiuchus. *Astron. Astrophys.* **593**, A111 (2016).

34. Fontanive, C. et al. A wide planetary-mass companion to a young low-mass brown dwarf in Ophiuchus. *Astrophys. J. Lett.* **905**, L14 (2020).

35. Greene, T. P. & Meyer, M. R. An infrared spectroscopic survey of the rho Ophiuchi Young stellar cluster: masses and ages from the H-R diagram. *Astrophys. J.* **450**, 233 (1995).

36. Sullivan, K. & Kraus, A. L. Undetected binary stars cause an observed mass-dependent age gradient in Upper Scorpius. *Astrophys. J.* **912**, 137 (2021).

37. David, T. J. et al. Age determination in Upper Scorpius with eclipsing binaries. *Astrophys. J.* **872**, 161 (2019).

38. Pecaut, M. J. & Mamajek, E. E. The star formation history and accretion-disc fraction among the K-type members of the Scorpius-Centaurus OB association. *Mon. Not. R. Astron. Soc.* **461**, 794–815 (2016).

39. Bouy, H. et al. Dynamical analysis of nearby clusters. Automated astrometry from the ground: precision proper motions over a wide field. *Astron. Astrophys.* **554**, A101 (2013).

40. Gaia Collaboration et al. Gaia Data Release 2. Summary of the contents and survey properties. *Astron. Astrophys.* **616**, A1 (2018).

41. Brown, A. G. A., Arenou, F., van Leeuwen, F., Lindegren, L. & Luri, X. Considerations in making full use of the HIPPARCOS catalogue. In *Hipparcos - Venice '97 Special Publication* 402 (eds Bonnet, R. M. et al.) 63–68 (ESA, 1997).

42. Mužić, K., Scholz, A., Geers, V., Jayawardhana, R. & Tamura, M. Substellar Objects in Nearby Young Clusters (SONYC). V. New brown dwarfs in ρ Ophiuchi. *Astrophys. J.* **744**, 134 (2012).

43. Ducourant, C. et al. Proper motion survey and kinematic analysis of the ρ Ophiuchi embedded cluster. *Astron. Astrophys.* **597**, A90 (2017).

44. Esplin, T. L., Luhman, K. L., Miller, E. B. & Mamajek, E. E. A WISE survey of circumstellar disks in the Upper Scorpius association. *Astron. J.* **156**, 75 (2018).

45. Damiani, F., Prisinzano, L., Pillitteri, I., Micela, G. & Sciortino, S. Stellar population of Sco OB2 revealed by Gaia DR2 data. *Astron. Astrophys.* **623**, A112 (2019).

46. Lodieu, N., Hambly, N. C. & Cross, N. J. G. Exploring the planetary-mass population in the Upper Scorpius association. *Mon. Not. R. Astron. Soc.* **503**, 2265–2279 (2021).

47. Baraffe, I., Homeier, D., Allard, F. & Chabrier, G. New evolutionary models for pre-main sequence and main sequence low-mass stars down to the hydrogen-burning limit. *Astron. Astrophys.* **577**, A42 (2015).

48. Marigo, P. et al. A New Generation of PARSEC-COLIBRI stellar isochrones including the TP-AGB phase. *Astrophys. J.* **835**, 77 (2017).

49. Barzalez Gagliuffi, D. C. et al. The Ultracool SpeXtoscopic Survey. I. Volume-limited spectroscopic sample and luminosity function of M7-L5 ultracool dwarfs. *Astrophys. J.* **883**, 205 (2019).

50. Gaia Collaboration et al. Gaia Early Data Release 3. The Gaia catalogue of nearby stars. *Astron. Astrophys.* **649**, A6 (2021).

51. Olivares, J. et al. Ruprecht 147 DANCE. I. Members, empirical isochrone, luminosity, and mass distributions. *Astron. Astrophys.* **625**, A115 (2019).

52. Salpeter, E. E. The luminosity function and stellar evolution. *Astrophys. J.* **121**, 161 (1955).

53. Thies, I. & Kroupa, P. A discontinuity in the low-mass initial mass function. *Astrophys. J.* **671**, 767–780 (2007).

54. Thies, I., Pflamm-Altenburg, J., Kroupa, P. & Marks, M. Characterizing the brown dwarf formation channels from the initial mass function and binary-star dynamics. *Astrophys. J.* **800**, 72 (2015).

55. Fernandes, R. B., Mulders, G. D., Pascucci, I., Mordasini, C. & Emsenhuber, A. Hints for a turnover at the snow line in the giant planet occurrence rate. *Astrophys. J.* **874**, 81 (2019).

56. Clanton, C. & Gaudi, B. S. Constraining the frequency of free-floating planets from a synthesis of microlensing, radial velocity, and direct imaging survey results. *Astrophys. J.* **834**, 46 (2017).

57. Bowler, B. P. Imaging extrasolar giant planets. *Publ. Astron. Soc. Pac.* **128**, 102001 (2016).

58. Suzuki, D. et al. The exoplanet mass-ratio function from the MOA-II Survey: discovery of a break and likely peak at a Neptune mass. *Astrophys. J.* **833**, 145 (2016).

59. Fressin, F. et al. The false positive rate of Kepler and the occurrence of planets. *Astrophys. J.* **766**, 81 (2013).

60. Wittenmyer, R. A. et al. Cool Jupiters greatly outnumber their toasty siblings: occurrence rates from the Anglo-Australian Planet Search. *Mon. Not. R. Astron. Soc.* **492**, 377–383 (2020).

61. Cumming, A. et al. The Keck Planet Search: detectability and the minimum mass and orbital period distribution of extrasolar planets. *Publ. Astron. Soc. Pac.* **120**, 531 (2008).

62. Butler, R. P. et al. Catalog of nearby exoplanets. *Astrophys. J.* **646**, 505–522 (2006).

63. Winn, J. N. & Fabrycky, D. C. The occurrence and architecture of exoplanetary systems. *Annu. Rev. Astron. Astrophys.* **53**, 409–447 (2015).

64. Jurić, M. & Tremaine, S. Dynamical origin of extrasolar planet eccentricity distribution. *Astrophys. J.* **686**, 603–620 (2008).

65. Chatterjee, S., Ford, E. B., Matsumura, S. & Rasio, F. A. Dynamical outcomes of planet-planet scattering. *Astrophys. J.* **686**, 580–602 (2008).

66. Raymond, S. N., Armitage, P. J. & Gorelick, N. Planet-planet scattering in planetesimal disks. II. Predictions for outer extrasolar planetary systems. *Astrophys. J.* **711**, 772–795 (2010).

67. Ford, E. B. & Rasio, F. A. Origins of eccentric extrasolar planets: testing the planet-planet scattering model. *Astrophys. J.* **686**, 621–636 (2008).

68. Ida, S., Lin, D. N. C. & Nagasawa, M. Toward a deterministic model of planetary formation. VII. Eccentricity distribution of gas giants. *Astrophys. J.* **775**, 42 (2013).

69. van Elteren, A., Portegies Zwart, S., Pelupessy, I., Cai, M. X. & McMillan, S. L. W. Survivability of planetary systems in young and dense star clusters. *Astron. Astrophys.* **624**, A120 (2019).

70. Nesvorný, D. Dynamical evolution of the early Solar System. *Annu. Rev. Astron. Astrophys.* **56**, 137–174 (2018).

71. Raymond, S. N., Izidoro, A. & Morbidelli, A. In *Planetary Astrobiology* (eds Meadows, V. et al.) 287–324 (University of Arizona Press, 2020).

72. Clement, M. S., Kaib, N. A., Raymond, S. N. & Walsh, K. J. Mars' growth stunted by an early giant planet instability. *Icarus* **311**, 340–356 (2018).

73. Morbidelli, A. et al. The timeline of the lunar bombardment: revisited. *Icarus* **305**, 262–276 (2018).

74. Parker, R. J. & Quanz, S. P. The effects of dynamical interactions on planets in young substructured star clusters. *Mon. Not. R. Astron. Soc.* **419**, 2448–2458 (2012).

75. Winter, A. J., Kruijssen, J. M. D., Longmore, S. N. & Chevance, M. Stellar clustering shapes the architecture of planetary systems. *Nature* **586**, 528–532 (2020).

76. Hester, J. J. et al. Hubble Space Telescope WFPC2 imaging of M16: photoevaporation and emerging young stellar objects. *Astron. J.* **111**, 2349 (1996).

77. Bouy, H. et al. A deep look into the cores of young clusters. I. σ-Orionis. *Astron. Astrophys.* **493**, 931–946 (2009).

78. Hodapp, K. W., Iserlohe, C., Stecklum, B. & Krabbe, A. σ Orionis IRS1 A and B: a binary containing a proplyd. *Astrophys. J. Lett.* **701**, L100–L104 (2009).

79. Paillassa, M., Bertin, E. & Bouy, H. MAXIMASK and MAXITRACK: two new tools for identifying contaminants in astronomical images using convolutional neural networks. *Astron. Astrophys.* **634**, A48 (2020).

80. Vandame, B. New algorithms and technologies for the un-supervised reduction of Optical/IR images. In *Astronomical Data Analysis II SPIE Conference Series Vol. 4787* (eds Starck, J.-L. & Murtagh, F. D.) 123–134 (SPIE, 2002).

81. Bertin, E. & Arnouts, S. SExtractor: software for source extraction. *Astron. Astrophys. Suppl. Ser.* **117**, 393–404 (1996).

82. Bertin, E. PSFEx: point spread function extractor (ASCL, 2013).

83. Bertin, E. Automatic astrometric and photometric calibration with SCAMP. In *Astronomical Data Analysis Software and Systems XV Astronomical Society of the Pacific Conference Series Vol. 351* (eds Gabriel, C. et al.) 112 (Astronomical Society of the Pacific, 2006).

84. Bertin, E. SWarp: resampling and co-adding FITS images together (ASCL, 2010).
85. Lawrence, A. et al. The UKIRT Infrared Deep Sky Survey (UKIDSS). *Mon. Not. R. Astron. Soc.* **379**, 1599–1617 (2007).

86. Baumgartd, H., Hilker, M., Sollima, A. & Bellini, A. Mean proper motions, space orbits, and velocity dispersion profiles of Galactic globular clusters derived from Gaia DR2 data. *Mon. Not. R. Astron. Soc.* **482**, 5138–5155 (2019).

87. Sarro, L. M. et al. Cluster membership probabilities from proper motions and multi-wavelength photometric catalogues. I. Method and application to the Pleiades cluster. *Astron. Astrophys.* **563**, A45 (2014).

88. Luhman, K. L., Herrmann, K. A., Mamajek, E. E., Esplin, T. L. & Pecaut, M. J. New young stars and brown dwarfs in the Upper Scorpius association. *Astron. J.* **156**, 76 (2018).

89. Maíz Apellániz, J. & Weiler, M. Reanalysis of the Gaia Data Release 2 photometric sensitivity curves using HST/STIS spectrophotometry. *Astron. Astrophys.* **619**, A180 (2018).

90. Miret-Roig, N. et al. IC 4665 DANCe. I. Members, empirical isochrones, magnitude distributions, present-day system mass function, and spatial distribution. *Astron. Astrophys.* **631**, A57 (2019).

91. Olivares, J. et al. Kalkayotl: a cluster distance inference code. *Astron. Astrophys.* **644**, A7 (2020).

92. Scott, D. W. Multivariate Density Estimation. *Theory, Practice, and Visualization* (Wiley, 1992).

93. Silverman, B. *W Density Estimation* (Chapman and Hall, 1986).

94. Baron, F. et al. Constraints on the occurrence and distribution of 1-20 M<sub>☉</sub> companions to stars at separations of 5-5000 au from a compilation of direct imaging surveys. *Astron. J.* **158**, 187 (2019).

95. Clanton, C. & Gaudi, B. S. Synthesizing exoplanet demographics: a single population of long-period planetary companions to M dwarfs consistent with microlensing, radial velocity, and direct imaging surveys. *Astrophys. J.* **819**, 125 (2016).

96. Planck Collaboration et al. Planck 2018 results. I. Overview and the cosmological legacy of Planck. *Astron. Astrophys.* **641**, A1 (2020).

## Acknowledgements

We thank P. Padoan, M. Bate and V.-M. Pelkonen for insightful comments on the comparison of our observational mass function to simulations, A. Howard, G. Mulders and C. Lovis for input on the occurrence rates and K. Peña Ramírez, A. Scholz and N. Lodieu for input on FFPS in star-forming regions. This research received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 682903, PI H.B.) and from the French State in the framework of the 'Investments for the future' Program, IdEx Bordeaux, reference no. ANR-10-IDEX-03-02. H.B. acknowledges financial support from the Canon Foundation in Europe. S.N.R acknowledges support from the CNRS's PNP programme. This research has been funded by Spanish State Research Agency (AEI) Projects PID2019-107061GB-C61 and MDM-2017-0737 Unidad de Excelencia 'María de Maeztu'-Centro de Astrobiología (CSIC/INTA). We gratefully acknowledge the support of the NVIDIA Corporation with the donation of one of the Titan Xp GPUs used for this research. This work is based on observations made with the INT operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. This paper makes use of data obtained from the Isaac Newton Group Archive which is maintained as part of the CASU Astronomical Data Centre at the Institute of Astronomy, Cambridge. This work is based on data obtained from the ESO Science Archive Facility and with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 065.1-0003, 065.1-0004, 065.1-0463, 065.0-0298, 071.A-9007(A), 071.A-9011(A), 075.C-0419(A), 075.D-0111(A), 075.D-0662(C), 079.A-9202(A), 079.A-9203(A), 079.A-9208(A), 079.D-0782(A), 079.D-0918(A), 080.A-9210(A), 081.A-9200(A), 081.A-9211(A), 081.A-9212(A), 082.A-9212(A), 082.C-0946(B), 083.A-9021(A), 083.A-9202(A), 083.A-9204(A), 083.C-0446(A), 085.A-9008(A), 085.A-9011(A), 085.C-0690(B), 085.D-0143(A), 086.C-0168(D), 088.C-0434(A), 091.A-0507(A), 091.C-0454(A), 093.A-9028(B), 094.A-9028(C), 096.A-9021(A), 097.A-9020(A), 097.A-9025(C), 164.O-0561(F), 60.A-9120(A), 67.A-0403(A), 68.D-0002(B), 68.D-0265(A), 69.A-0615(B), 69.C-0182(A), 69.C-0260(A), 69.C-0426(C), 69.D-0582(A), 71.C-0580(A), 71.C-0580(B), 71.D-0014(A), 081.A-0673(A), 083.A-0321(A), 085.C-0841(E), 085.C-1099(A), 089.D-0952(B), 089.C-0952(C), 089.C-0952(E), 089.D-0291(A), 091.A-0703(B), 091.C-0543(B), 091.C-0543(C), 091.C-0543(D), 091.C-0543(E), 092.C-0548(F), 195.B-0283(A), 60.A-9283(A), 60.A-9800(L), 60.A-9800(H), 083.C-0556(A), 279.C-5062(C), 093.B-0280(B), 095.D-0494(A), 096.C-0730(A), 097.C-0749(A), 098.C-0850(A), 099.C-0474(A), 177.D-3023(G), 60.A-9038(A), 088.D-0675(A), 089.C-0102(A), 089.C-0102(B), 089.C-0102(C), 095.D-0038(A), 097.C-0781(A), 179.A-2010(H), 179.A-2010(J), 179.A-2010(K), 179.A-2010(L), 179.A-2010(N), 198.C-2009(A), 198.C-2009(B), 198.C-2009(F), 198.C-2009(H), 198.C-2009(I), 60.A-9292(A). This research uses services or data provided by the NOIRLab's Astro Data Archive. The NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA), Inc. under a cooperative agreement with the National Science Foundation. The Community Science & Data Center is the place where the Astro Data Archive is being developed and operated the Community Science & Data Center is the place where the Astro Data Archive is being developed and operated. This project used data obtained with the Dark Energy Camera (DECam), which was constructed by the Dark Energy

Survey (DES) collaboration. Funding for the DES Projects has been provided by the DOE and NSF (USA), MISE (Spain), STFC (UK), HEFCE (UK), NCSA (UIUC), KICP (U. Chicago), CCAPP (Ohio State), MIFPA (Texas A&M), CNQP, PAPERJ, FINEP (Brazil), MINECO (Spain), DFG (Germany) and the Collaborating Institutions in the Dark Energy Survey, which are Argonne Lab, UC Santa Cruz, University of Cambridge, CIEMAT-Madrid, University of Chicago, University College London, DES-Brazil Consortium, University of Edinburgh, ETH Zürich, Fermilab, University of Illinois, ICE (IEEC-CSIC), IFAE Barcelona, Lawrence Berkeley Lab, LMU München and the associated Excellence Cluster Universe, University of Michigan, NOIRLab, University of Nottingham, Ohio State University, OzDES Membership Consortium, University of Pennsylvania, University of Portsmouth, SLAC National Lab, Stanford University, University of Sussex, and Texas A&M University. Based on observations at Cerro Tololo Inter-American Observatory at NSF's NOIRLab (NOIRLab Prop. 2012B-0569 (PI: Allen); 2013A-0214 (PI: Berger); 2013A-0327 (PI: Rest); 2013A-0351 (PI: Dey); 2013A-0723 (PI: Mamajek); 2013A-0737 (PI: Sheppard); 2013A-9999 (PI: Walker); 2013B-0325 (PI: Vivas); 2013B-0531 (PI: Mamajek); 2013B-0536 (PI: Allen); 2014A-0035 (PI: Bouy); 2014A-0239 (PI: Sullivan); 2014A-0306 (PI: Dai); 2014A-0327 (PI: Rest); 2014A-0412 (PI: Rest); 2014A-0479 (PI: Sheppard); 2014A-0480 (PI: Rich); 2014A-0634 (PI: James); 2015A-0151 (PI: Calamida); 2015A-0205 (PI: Mamajek); 2015A-0371 (PI: Rest); 2015A-0397 (PI: Rest); 2015A-0610 (PI: Fuentes); 2015B-0307 (PI: Rest); 2016A-0189 (PI: Rest); 2016A-0327 (PI: Finkbeiner); 2016B-0279 (PI: Finkbeiner); 2016B-0301 (PI: Rest); 2016B-0301 (PI: Zenteno); 2017A-0002 (PI: Bouy); 2017A-0388 (PI: Zenteno); 2017A-0389 (PI: Rest); 2017A-0389 (PI: Tucker); 2017A-0918 (PI: Yip); 2017B-0285 (PI: Rest); 2018A-0059 (PI: Bouy); 2018A-0251 (PI: Finkbeiner); 2019A-0060 (PI: Bouy); 2019A-0101 (PI: Hartigan); 2019A-0305 (PI: Drlica-Wagner); 2019A-0337 (PI: Trilling); 2019B-0323 (PI: Zenteno)), which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. Based in part on observations at Kitt Peak National Observatory at NSF's NOIRLab (NOIRLab Prop. ID 2013A-0399 (PI: Adam); 2014A-0642 (PI: Ronald)), which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. The authors are honored to be permitted to conduct astronomical research on Iolkam Du'ag (Kitt Peak), a mountain with particular significance to the Tohono O'odham. Based in part on observations made at Cerro Tololo Inter-American Observatory at NSF's NOIRLab (NOIRLab Prop. ID 2011A-0368 (PI: Stringfellow); 2010A-0475 (PI: Stringfellow); 2011A-0603 (PI: Tilvi); 2011A-0644 (PI: Catalan); 2010A-0036 (PI: Probst); 2005A-0183 (PI: Ridge); 2006A-0139 (PI: Ridge); 2006B-0021 (PI: Grindlay); 2007A-0180 (PI: Martin); 2007A-0514 (PI: Mardones); 2007A-0599 (PI: Huard); 2008B-0368 (PI: Zuckerman); 2008B-0909 (PI: Vaduvescu); 2010A-0260 (PI: Faherty); 2010A-0326 (PI: Allers); 2010A-0482 (PI: Miller), which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. This work is based in part on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan and obtained from the SMOKA, which is operated by the Astronomy Data Center, National Astronomical Observatory of Japan. The Hyper Suprime-Cam collaboration includes the astronomical communities of Japan and Taiwan and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA) and Princeton University. Funding was contributed by the FIRST program from Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA and Princeton University. Based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. Based on CFH12K observations obtained at the Canada-France-Hawaii Telescope (CFHT). Based on observations obtained with WIRCam, a joint project of CFHT, Taiwan, Korea, Canada, France, and the Canada-France-Hawaii Telescope (CFHT). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software at http://dm.lsst.org. The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg, the Max Planck Institute for Extraterrestrial Physics, Garching, Johns Hopkins University, Durham University, the
University of Edinburgh, Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under grant no. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under grant no. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE) and the Los Alamos National Laboratory. This work is based on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center at National Astronomical Observatory of Japan. This work is based on observations obtained with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA and Canada.

## Author contributions

N.M.-R. and H.B. led the observations, data analysis and scientific analysis. H.B. is the PI of the COSMIC-DANCE project. S.N.R. contributed to the scientific analysis and the discussion about planet ejection and formation. M.T. led the Subaru observations used in this study. E.B. wrote the software packages used to process and analyse the images. J.-C.C., P.G., D.B. and N.H. led the observations at various observatories.

J.O., L.M.S. and A.B. led the development of the probabilistic method and software used to identify members.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41550-021-01513-x.

**Correspondence and requests for materials** should be addressed to Núria Miret-Roig.

**Peer review information** *Nature Astronomy* thanks the anonymous reviewers for their contribution to the peer review of this work.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2022