IC 4665 DANCe I. Members, empirical isochrones, magnitude distributions, present-day system mass function, and spatial distribution

Context. The study of star formation is extremely challenging due to the lack of complete and clean samples of young, nearby clusters, and star forming regions. The recent Gaia DR2 catalogue complemented with the deep, ground based COSMIC DANCe catalogue offers a new database of unprecedented accuracy to revisit the membership of clusters and star forming regions. The 30 Myr open cluster IC 4665 is one of the few well-known clusters of this age and it is an excellent target where to test evolutionary models and study planetary formation. Aims. We aim to provide a comprehensive membership analysis of IC 4665 and to study the following properties: empirical isochrones, distance, magnitude distribution, present-day system mass function, and spatial distribution. Methods. We use the Gaia DR2 catalogue together with the DANCe catalogue to look for members using a probabilistic model of the distribution of the observable quantities in both the cluster and background populations. Results. We obtain a final list of 819 candidate members which cover a 12.4 magnitude range (7<J<19.4). We find that 50% are new candidates, and we estimate a conservative contamination rate of 20%. This unique sample of members allows us to obtain a present-day system mass function in the range of 0.02-6 Msun, which reveals a number of details not seen in previous studies. In addition, they favour a spherically symmetric spatial distribution for this young open cluster. Conclusions. Our membership analysis represents a significant increase in the quantity and quality (low-contamination) with respect to previous studies. As such, it offers an excellent opportunity to revisit other fundamental parameters such as the age.


Introduction
The Initial Mass Function (IMF), i.e. the frequency distribution of stellar masses at birth, is a fundamental parameter in the study of stellar formation and evolution. It was first introduced by Salpeter (1955) in the form of ξ(log 10 m) = dn/d log 10 m and in this study we adopt the same formalism. Since then, strong efforts have been put on trying to constrain its shape in various environments. While it is fairly well known for intermediate masses, the extremes of the IMF remain uncertain. For the high mass domain the main difficulty is the low number of stars and their fast evolution, while for the low mass domain the main diffi-culty is the high level of contamination and incompleteness even with the best photometric and astrometric surveys.
The study of the IMF requires an accurate and comprehensive census of the cluster members. In this study, we propose to derive such a census of one nearby young cluster namely, IC 4665. For this purpose, we use wide-field deep catalogues encompassing a large area around the cluster, and analyse them using a powerful classification algorithm capable of identifying the few hundreds of cluster members within the millions of interlopers. Hereafter, we refer as members to the candidate members resultant from the membership analysis. There are few well-known, nearby (< 500 pc), pre-main sequence open clusters in the age interval 10-50 Myr. One of them is IC 4665, located in the Ophiuchus constellation and first reported by Philippe Loys de Chéseaux in 1745. Its age was estimated using the lithium depletion boundary at 27.7 +4.2 −3.5 Myr by Manzi et al. (2008). From pre-main sequence isochrone fitting and upper-main-sequence turn-off fitting, Cargile & James (2010) derived an age and a distance of 36 ± 9 Myr and 360 ± 12 pc and 42 ± 12 Myr and 357 ± 12 pc, respectively.
The first study of the IMF of IC 4665 was carried out by de Wit et al. (2006). They selected their members from photometric observations in the optical obtained at the Canada-France-Hawaii Telescope (CFHT). They estimated a contamination by foreground and background stars of up to 85% using control fields, which can be explained by its low galactic latitude. They reported a mass function best described by a power law with an exponent of -0.6 for the low mass objects down to ∼0.1 M . Later, Lodieu et al. (2011) performed a similar analysis adding near-infrared photometry from the UKIRT Infrared Deep Sky Survey (UKIDSS, Lawrence et al. 2007) to the previous observations of de Wit et al. (2006). They revised the members of previous studies as well as proposed new candidate members. They reported a mass function best represented by a log-normal function with a peak at 0.25-0.16 M . The differences between the mass functions obtained with these two studies can be mainly attributed to the large contamination rate of field stars, as we shall see later.
Recently, the second Gaia data release (Gaia Collaboration et al. 2018b, hereafter Gaia DR2) was public providing the fiveparameter astrometric solution (positions on the sky, parallaxes, and proper motions) as well as G, G BP and G RP magnitudes for more than 1.3 billion sources, with a limiting magnitude of G ≈ 21 mag at the faint end and G ∼ 3 mag at the bright end. The average astrometric precision is of the order of the mas yr −1 in proper motion and below the mas in parallax, and the average photometric precision is at the milli-magnitude level. This constitutes an unprecedented accurate astrometric+photometric dataset ideally suited to study the census of nearby open clusters. A first demonstration of the power of the Gaia data is presented in Gaia Collaboration et al. (2018a). They studied the fine structures of the Hertzsprung-Russell diagram (HRD) in the field and in open and globular clusters. IC 4665 was among their targets and they provided a list of 174 high-probability members up to magnitude G < 18. Soon after, Cantat-Gaudin et al. (2018) presented another study of open clusters using Gaia DR2 data. They derived another membership list (with the same magnitude limit) made of 175 high-probability members, 146 of which are in common with Gaia Collaboration et al. (2018a). Both studies used only the Gaia data, applied a strict filtering, and discarded sources fainter than G = 18 mag, thus delivering a clean but yet highly incomplete sample.
Over the past few years, Bouy et al. (2013) started a survey program, the DANCe project, with the aim of deriving a comprehensive and homogeneous census of stellar and sub-stellar sources in the nearby (< 1 kpc), young (< 500 Myr) clusters. This survey combines deep, wide-field, multi-epoch images obtained at various observatories to build a catalogue of accurate proper motions and multi-wavelength photometry with a sensitivity up to 5 magnitudes deeper than Gaia and including near-infrared data. Here we present the DANCe catalogue for the region of IC 4665. After identifying candidate members, we study the cluster properties and in particular its empirical isochrones, distance, magnitude distribution, present-day mass function (PDMF), and spatial distribution. The assumptions in this work are based on the properties of the dataset and the cluster. We strongly recommend the reader to look at Sect. 2 of Olivares et al. (2019) for more details.
This paper is structured as follows. First, we introduce the two datasets used in this study (Sect. 2). Second, we present the algorithm we use for the membership analysis (Sect. 3) and discuss the results obtained including a comparison with previous studies (Sect. 4). Then, we provide the empirical iscohrones of this young cluster and compare them with evolutionary models (Sect. 5). In addition, we compute the apparent magnitude distribution and the PDMF (Sect. 6). Finally, we study the spatial distribution of the cluster (Sect. 7), and we draw our conclusions (Sect. 8).

The data
In this work we used two different datasets with different origins and properties to look for members in the IC 4665 open cluster. In this section, we describe how we obtained each of them.

The Gaia DR2 dataset
We queried a circular area of 3 • radius around the centre of the cluster (RA = 266.6 • , Dec = 5.7 • ), from the Gaia DR2 catalogue (see Appx. A.1). Then, we only kept the sources with a full five-parameter solution available. Some quality checks have been suggested in the literature. We did not apply any filtering techniques because we want to be as complete as possible. The filtering recommended by the Gaia team is based on the renormalised unit weighted error (RUWE) and is described in detail in a publicly available technical note 1 . The RUWE criterion is a quality indicator which might be used when the aim is to have only the most precise, reliable, and consistent astrometric solutions. However, it also leads to a higher degree of incompleteness. For instance, since the Gaia DR2 catalogue does not deal with binaries, their solution is likely to be "inconsistent", and thus the RUWE filter will remove most of the already few binaries included in Gaia DR2. We therefore have no strong scientific argument to cut our sample by this or any other kind of filtering. In addition, the sources with problematic astrometric solutions, may be rejected later on based on complementary observations and/or subsequent Gaia DRs.
This sample contains positions, proper motions, parallaxes, and G, G BP , G RP photometry for 1 217 725 sources and, hereafter, we refer to it as the GDR2 catalogue 2 . The mean errors of this catalogue are ∼ 0.5 mas for parallaxes, ∼ 1 mas yr −1 for proper motions, and < 0.1 mag for the photometry. According to Gaia Collaboration et al. (2018b), the catalogue is mostly complete down to G = 7 mag. On the faint side, Lindegren et al. (2018) reported that the five-parameter solution is 94.5% complete up to G = 19 mag (see their Table B.1). In the following, we therefore assume that the GDR2 catalogue is complete between 7 < G < 19 mag.

The DANCe dataset
We searched the European Southern Observatory (ESO) archive, the National Optical Astronomy Observatory (NOAO) archive,  (b) the chip layout has large gaps between detectors, and the coverage of the focal plane is only partial (c) one of the 12 detectors is dead  Thibault et al. (2002) the PTF archive hosted at the NASA/IPAC Infrared Science Archive (IRSA), the Canadian Astronomy Data Centre (CADC) archive, the Isaac Newton Group (ING) archive, the WFCAM Science (WSA) archive, and the SMOKA science archive for wide field images within a circular region of 3 • radius, centred on IC 4665. The data found in these public archives was complemented with our own observations using the Las Campanas Swope telescope and its Direct CCD camera, the Dark Energy Camera (DECam) mounted on the Blanco telescope at the Cerro Tololo Inter-American Observatory (CTIO), the NEW-FIRM camera mounted on the 4m telescope at the Kitt Peak National Observatory (KPNO), the Hyper Suprime-Cam (HSC) mounted on the Subaru telescope at the National Astronomical Observatory of Japan, and the Wide Field Camera (WFC) mounted on the Isaac Newton Telescope (INT). A number of observations found in the archives were discarded after a visual inspection, because of their poor quality, limited sensitivity or acquisition problems. Table 1 gives an overview of the various cameras used for this study.
The airmass and the Full-Width at Half Maximum (FWHM) measured in the images using point-like sources are two important parameters influencing the achievable astrometric accuracy. About 90% of the observations were obtained at airmass ≤ 1.5 (see Fig. 1 top). IC 4665 is located at a declination of δ ∼ 5 • and we gathered data from both hemispheres. About 82% of the images have FWHM ≤ 1 , and 90% have FWHM ≤ 1 . 2 (see Fig. 1 bottom).
In all cases except for MegaCam, WIRCam, PTF, DECam, UKIRT and HSC images, the raw data and associated calibration frames were downloaded and processed using standard procedures using an updated version of Alambic (Vandame 2002), a software suite developed and optimised for the processing of large multi-CCD images. In the case of CFHT/MegaCam, the images processed and calibrated with the Elixir pipeline were retrieved from the CADC archive (Magnier & Cuillandre 2004).
The WIRCam images processed with the official 'I'iwi pipeline were retrieved from the CADC archive. In the case of DE-Cam, the images processed with the community pipeline (Valdes et al. 2014) were retrieved from the NOAO public archive. The pipeline processed PTF images were downloaded from the IPAC archive. UKIRT images from the UKIDSS and UHS surveys (Dye et al. 2018) processed by the Cambridge Astronomical Survey Unit were retrieved from the WFCAM Science Archive. Finally, the HSC raw images were processed using the official HSC pipeline (Bosch et al. 2018).

Astrometric analysis
After a visual rejection of problematic images (mostly due to loss of guiding, tracking or electronics problems), the dataset included 6 774 individual images originating from 13 instruments. The total amount of data (scientific images, associated calibrations, and intermediate products) added to almost 20TB.
The astrometric calibration was performed as described in Bouy et al. (2013). The recently released Gaia DR2 catalogue was used as external astrometric reference instead of the 2MASS catalogue, leading to a much improved astrometric solution.
The final average internal and external 3σ residuals add up to ∼25 mas for high signal-to-noise (photon noise limited) sources. As explained in Bouy et al. (2013), the proper motions computed are relative and display an offset with respect to the ICRS. We estimate the offset by computing the median offset between our and the Gaia DR2 proper motion measurements after rejecting outliers using the modified Z-score (Iglewicz & Hoaglin 1993). We find offsets of (∆µ α cos δ, ∆µ δ ) = (1.70, 4.48) mas yr −1 . The uncertainty on this offset is estimated using bootstrapping and is found to be negligible (<0.003 mas yr −1 ).
Given the superiority and robustness of Gaia measurements compared to our ground-based measurements, the Gaia DR2 proper motion measurements are always preferred when avail- able. Therefore, we cross-matched our catalogue with the Gaia DR2 using a 1 radius and we kept all the Gaia proper motions when available. The median uncertainties in proper motions are ∼2 mas yr −1 .
We found that about 1.3% of the sources (∼60 000) were duplicated in the final catalogue. A visual inspection showed that they were almost all very low signal-to-noise and that the SExtractor deblending algorithm resolved them as two sources instead of one, in one (or a few) images. These resolved sources later fooled the cross-identification algorithm and ultimately resulted in two independent sources instead of one. There is no straightforward solution to this problem as for now, but given their very small number we treated them as regular sources in the rest of the analysis and simply looked for duplicated entries in the final members list.

Photometric analysis
The photometric calibration was performed only for the g, r, i, z, y and J, H, K s images. It was not attempted for the INT images obtained in any other filter, the ESO2.2m WFI images, the PTF images (because the camera has a significantly coarser pixel scale and the images reach a depth shallower than Pan-STARRS), and the CPAPIR I-band images.
The photometric zeropoint of all individual images was computed by direct comparison of the instrumental SExtractor MAG_AUTO magnitudes with an external catalogue: -J, H, K s images were tied to the 2MASS catalogue, g, r, i, z, y images were tied to the Pan-STARRS PS1 first release.
The procedure followed to derive the individual zeropoints is described in Olivares et al. (2019). Briefly, the zeropoints are computed as the median of the difference between the instrumental magnitude and the measurements of the closest match within 1 in the reference catalogue after rejecting outliers using the modified Z-score criterion. We find typical 1σ dispersions of the order of 0.03-0.08 mag depending on the filter. We median-combined all the images obtained with the same camera and in the same filter to build a deep stack and extracted the corresponding photometry. This allowed us to significantly improve the sensitivity in all filters and recover or improve the photometry of faint sources obtained in the individual images.
As in Bouy et al. (2013), we complemented the photometry extracted from the images with that of external catalogues: Gaia DR2 (G, G BP , G RP ), Pan-STARRS (grizy), 2MASS (JHK s) and ALLWISE (all 4 bands) to improve the spatial and wavelength coverage of the final dataset (see Appx. A.1-A.4 for the queries used). The corresponding photometric measurements were either added to our catalogue when no measurement was available in our data, or the weighted average of all measurements (our and external catalogues) was computed after rejecting outliers using the modified Z-score criterion.   The differences between the transmission of the various instruments contribute to the photometric dispersion (e.g. see Fig. 2). Moreover recent variability studies of young clusters found typical amplitudes of 0.03 mag (e.g. Rebull et al. 2016Rebull et al. , 2018. As explained in Olivares et al. (2019) we add quadratically 0.05 mag to all the photometric measurement uncertainties in our catalogue to take into account these sources of uncertainties in our membership analysis. The final mean uncertainties in photometry depend on the photometric band, and are of the order of 0.07 mag.
Given the low extinction in that area of the sky, the maximum of the magnitude distributions ( Fig. 3) gives an estimate of the completeness limit of the survey. It is nevertheless important to remember that the spatial coverage of the various instruments is not homogeneous and the depth of the survey varies spatially. The limit of sensitivity of the external catalogues merged with our data (2MASS and Pan-STARRS) are sometimes visible as secondary maxima. In Table 2 we give, for each photometric band, the number and percentage of measurements as well as the completeness limits we use in this study. This final catalogue (hereafter the DANCe catalogue) includes 2 358 937 sources.

Membership analysis
To select candidate members we used the methodology originally developed by Sarro et al. (2014) and updated by Olivares et al. (2019). Briefly, this algorithm separates all the sources within two populations, namely the cluster and the field. The field model is computed once at the beginning and fixed thereafter, based on the assumption that the few hundreds of cluster members do no affect significantly the field population model. The cluster model is built iteratively and, at each iteration, both models are used to reclassify the sources until convergence. This algorithm takes into account the covariance matrix of the astrometric parameters when available. To model the cluster and field populations, the algorithm only uses the complete sources, i.e. the sources with measurements available for all the observables. The coverage and sensitivity of the different photometric bands is therefore a key issue in our analysis. The final model allows us to compute a membership probability to incomplete sources after marginalisation over the missing values. To start the analysis and build the first model we need a catalogue of sources for the region of interest and an initial list of members. The latter can be slightly contaminated and incomplete, and serves only to define the cluster locus in the multi-dimensional space in the first iteration.
Because the field and cluster models are built from sources with complete photometric and astrometric measurements, a simultaneous analysis of the GDR2 and DANCe catalogues would not be optimum. Many faint DANCe sources would have missing parallaxes and Gaia photometry while bright GDR2 sources, saturated in our images, would have missing DANCe photometry, making it hard (if at all possible) to define a proper representation space in which a sufficient number of sources have complete measurements. As described in Olivares et al. (2019) we therefore decided to analyse both catalogues independently.

Initial members
Recently, two studies have published members of IC 4665 using the Gaia DR2 data. The work of Gaia Collaboration et al. (2018a) published a list of 174 members and the work of Cantat-Gaudin et al. (2018) published a list of 175 members. Both studies have a magnitude limit of G = 18, and most of the sources in common. We combined their results and obtained a list of 203 members which we used as initial list for our membership analysis of the GDR2 catalogue.
To start the membership algorithm for the DANCe catalogue, we used the members we obtained with the GDR2 catalogue which have a counterpart in DANCe. In this case, the initial list does not cover the full magnitude range of the catalogue (DANCe goes fainter than the initial list of GDR2 members). However, our algorithm extends the initial principal curve (see the cluster model for the photometry in Sect 3.4) by progressively and iteratively extrapolating the photometric sequence to fainter regions. These fainter regions of extrapolation are small and the new candidate members found with them are added and used in subsequent iterations to better define (or correct if necessary) the extrapolation. The extrapolation of the photometric curve is guided by the astrometry which does not change with magnitude. This extrapolation of the principal curve is further explained in Sarro et al. (2014) where we presented for the first time the membership algorithm.

Representation space
The representation space is the set of astrometric and photometric variables we use for the membership analysis. We always use proper motions, as well as parallaxes in the case of the GDR2 catalogue. The photometric variables are chosen according to their importance calculated from a random forest algorithm, as in Olivares et al. (2019). The largest the amount of features, the more information there is to classify the sources between members and non-members. At the same time, the bands with a large number of missing observations are avoided.
With the representation space established, the whole data set is split between complete and incomplete sources. A source is said to be complete when it has a measurement for all the variables of the chosen representation space. In consequence, different representation spaces lead to different complete/incomplete ratios.
For the analysis with the GDR2 catalogue, the representation space we used is pmra, pmdec, parallax, G RP , G BP − G, G −G RP . With this representation space 1 184 922 sources (97%) have complete data.
For the analysis with the DANCe catalogue, the representation space we used is pmra, pmdec, J, i − z, i − y, i − J. With this representation space, 1 627 593 sources have observations in all the photometric bands, which represents a 69% of the catalogue. We decided not to include the g, r, H, and K s bands in the representation space because of the large number of sources with missing photometry (see Table 2).

Field model
The field population is modelled with a Gaussian Mixture Model (hereafter GMM) in the whole representation space. The field model is computed once at the beginning and fixed thereafter. We explored GMMs with different number of components (60, 80, 100, 120, 140, 160, and 180). We choose the optimum model as the simplest one which minimises the Bayesian information criterion (BIC). This results in 100 components for the analysis with the GDR2 and the DANCe catalogue.

Cluster model
The cluster model is a product of two independent models: a GMM for the astrometry and a principal curve in photometry. The astrometric model is a GMM and we choose the number of components that optimises the BIC between 1 and 4.
The cluster model is computed iteratively, starting from the initial list of members. At each iteration, we compute independently a model for the astrometry and a model for the photometry. Then, we assign Bayesian probabilities of membership to all the complete sources. These probabilities together with a probability threshold, p in , are used to reclassify the complete sources between members and non-members. The p in is a free parameter of the model which defines the degrees of completeness and contamination that we desire for the training set (and as a consequence, for the final list of members). We refer the interested reader to Sarro et al. (2014) and Olivares et al. (2019) for a more detailed description of this parameter. Then, the cluster model is recomputed based on the new members list and we repeat this process until convergence.
When the model has converged, we generate a synthetic dataset from the model learnt with observed data. Therefore it has similar properties as the latter (e.g. missing values, frequency of members). As a consequence, the results derived from them are restricted to the used RS and learnt model.
We use this synthetic dataset to analyse the goodness of our classification as well as to compute the optimum probability threshold, p opt , which is used to do a final classification. The optimum threshold will of course depend on the scientific goal behind the membership analysis. In our case, to study the mass function, we are interested in reaching a compromise between the contamination and the true positive rate. For that, we choose as p opt the one which minimises the distance to the perfect classifier (DST). This distance is defined in terms of the contamination rate (CR) and the true positive rate (TPR) which in turn depend on the the confusion matrix: true positives (TP), false positives (FP), false negatives (FN), and true negatives (TN). This indices are defined as follows: As we have mentioned, the estimations which can be obtained with this synthetic data set are restricted to the same conditions as the observations and to the assumption that the model correctly represents the observed data. Thus, the measured CR and TPR can be underestimated and overestimated respectively, with respect to those obtained with better quality data and realistic models. As we shall see in Sect. 3.6, this is observed in the CR of DANCe. The value estimated using synthetic data appears to be underestimated with respect to the one obtained using the GDR2 members as reference. For more details we refer to Olivares et al. (2019). Table 3. Internal probability threshold (p in ), optimum probability threshold (p opt ) and number of complete and incomplete members (CM, IM) for each of the membership analyses, namely GDR2 and DANCe. We run the full model considering several p in thresholds (0.5, 0.6, 0.7, 0.8, and 0.9) and, for each, we compute the optimum threshold p opt , using synthetic data. In Table 3 we show the p in , p opt , the number of complete and incomplete members, for each of the two independent analysis (GDR2 and DANCe).

Classification of incomplete sources
We used the field and cluster models described in the previous subsections to compute membership probabilities for the incomplete sources, i.e. the sources which lack one or more magnitudes of the representation space. Then, we used the optimum threshold to classify all the sources between members and nonmembers.
For the GDR2 catalogue, there are very few incomplete sources and none of them is classified as member. For the DANCe catalogue, the number of incomplete sources classified as members is 4-8 depending on the p in . In general, they lack z and/or y photometry and the brightest ones are also classified as members by the analysis with GDR2.

Final members lists
Tables 4 and 5 (available at CDS) give the GDR2 and DANCe catalogues used in this work, respectively. For the GDR2 catalogue, we provide the Gaia DR2 source ID and the sky positions, and for the DANCe catalogue, we compile all the astrometry and photometry presented in Sect. 2. In addition, for each catalogue, we provide the posterior membership probabilities obtained with the different p in discussed in this section. Here, we describe our strategy to choose the members list most convenient to our necessities. However, we encourage the interested reader to chose the members list most convenient to his/her goal.
The membership probabilities obtained with different p in values have to be compared with care. The relation between different membership probabilities (obtained with different p in ) is not linear, and lower p in tend to provide higher membership probabilities. In general, the models computed with lower p in permit a larger inclusion of sources initially classified as field into the cluster class during the training of the model. This results in lists of members which can include a significant amount of contaminants. On the contrary, models computed with higher p in are more restrictive in including additional sources into the cluster model and thus, tend to have lower contamination but at the same time can become incomplete.
The membership probabilities we compute are not absolute but tightly related to the model used which at the same time depends on the representation space and on the p in (desired degree of completeness and contamination). In consequence, the com-parison between GDR2 and DANCe membership probabilities is not straightforward, specially due to the different representation spaces of each catalogue (DANCe does not have parallaxes).
To study cluster properties such as the mass function, we want a unique list of members, the cleanest and most complete possible. This means first, choosing a p in for each study (GDR2 and DANCe), and then combining both catalogues. We begin with the analysis of GDR2, which is expected to be more robust since it includes a very discriminating variable: the parallax. We used a Kolmogorov-Smirnov (KS) test in different variables of the representation space (proper motions, parallax, and photometry) to see if the distributions of these variables obtained with different p in were compatible among them. The goal is to see if we find signs of a strong contamination or a strong incompleteness in one or several of the lists with respect to the others. We started by taking the sources classified as members (i.e. those with p > p opt ) obtained with the model trained with p in = 0.9 as a reference. This is the most conservative and the least contaminated, but also probably the most incomplete list of members. Then, we compared this list of members with all the rest (p in = 0.5, 0.6, 0.7, and 0.8, one at each time). The KS test showed that for the lists of p in = 0.5 and 0.6, there is evidence that the distribution of their proper motions and parallax do not come from the same distribution as the one obtained with the list of p in = 0.9 with a p-value lower than the significance level 0.01. On the contrary, for the distribution of the astrometric variables coming from the lists of p in = 0.7 and 0.8, the KS test shows no evidence to reject they come from the same distribution as the one obtained with p in = 0.9 with p-values of 0.4-0.5. Then, we investigated which is the reason of the incompatibility of the lists of p in = 0.5 and 0.6 with respect to the rest. We found that the parallax and proper motions distributions obtained with the lists of p in = 0.5 and 0.6 have significantly more extended wings than the distributions obtained with the lists of higher p in . We interpreted this as contamination and therefore, we discarded these two lists. The remaining lists are compatible according to the KS test so we chose the list of p in = 0.7 which has the largest number of members.
To select the optimum p in for the DANCe analysis we also applied a KS test to find which distributions were compatible with the one obtained with p in = 0.9. In this case, there was no evidence to reject that the distributions of all the variables analysed for the lists from all p in come from the same distribution as the ones of p in = 0.9 since all the p-values were > 0.3. Then, to check the consistency between the GDR2 and DANCe lists, we took the members of GDR2 p in = 0.7 as a reference and compared them with the members recovered in the different DANCe lists, in the region where both studies are complete (14.5 G 19). We found that all the DANCe lists recovered roughly the same number of GDR2 members (250-266 from 285, ∼ 90%). On the contrary, the number of members in DANCe which were not in GDR2 decreased with increasing p in . These sources have parallaxes which are in general incompatible with the GDR2 members (further than 3σ), and thus we believe that most of them are contaminants (they represent a 30-35% of the DANCe members). In short, for the DANCe analysis we do not find any strong argument to discard any list. Therefore, we decided to keep the one with largest amount of members, i.e. the one with p in = 0.5 keeping in mind that it includes a contamination of the order of 30-35%, estimated from the comparison with the GDR2 members.
In the rest of this work, we will be using the GDR2 list with a p in = 0.7 and the DANCe list of p in = 0.5. Those two lists add Table 4. GDR2 catalogue of IC4665 (only the first 10 rows are displayed as example). Columns indicate: (1) Gaia DR2 source ID; (2-3) right ascension and declination; (4) Gaia DR2 G-band magnitude; (5-9) posterior probabilities obtained with p in from 0.5 to 0.9.  to a final list of 819 members which is analysed in more detail in the next section.

Members comparison
In this section we analyse and discuss our final members list. We compare the results obtained with the GDR2 and DANCe catalogues, and we also compare our final list with the members already reported in the literature.

Comparison GDR2 vs. DANCe members
Here we compare the two membership analyses obtained with the GDR2 and the DANCe catalogues. We cross-matched the two catalogues (which contain both members and field stars) and found 1 211 272 sources in common. In Figure 4 we compare the membership probabilities obtained with the two catalogues. The diagonal line represents the one-to-one relation, and the vertical and horizontal dashed lines represent the optimum thresholds. We see that most of the sources are clustered in the bottom left (field) and top right (cluster) regions of the diagram. Nonetheless, there are sources which are classified as members by one of the two studies and not by the other.
To understand the differences between the two classifiers we represented the number of members as a function of the magni-tude (Fig. 5). We distinguish between the members obtained with both classifiers (red), the members only from the GDR2 analysis (blue), and the members only from the DANCe analysis (green). Here we discuss the 4 possible cases regarding the results of the two membership analyses.

Members in GDR2 and DANCe
There are 331 sources which appear as members in both analyses (red in Fig 5). In the magnitude range where both catalogues are complete, the majority of members are classified as so by the two analyses.

Members in GDR2 only
There are 103 sources which appear as members in the GDR2 analysis but not in DANCe (blue in Fig 5). The majority of these are the brightest sources which are saturated in the DANCe catalogue. Most of them have J, H and K s photometry but lack i, z and y which are essential bands for the representation space we use in DANCe.
In the magnitude range where both analyses are complete, we see that the members obtained only with GDR2 have a low and flat distribution. This can be interpreted as the GDR2 members list having a very low contamination, which does not depend on the magnitude. The reason is that the parallax is the most discriminating variable to classify these members. Although the uncertainties on the parallax depend on the magnitude, they are at the level required to distinguish the cluster from the field along all the magnitude range. At magnitudes > 18 we see a slight increase but it is not significant.

Members in DANCe only
There are 385 sources which appear as members in DANCe but not in GDR2 (green in Fig 5). From these, 120 objects (31%) do not have the five-parameter solution in Gaia and 186 (48%) have parallax uncertainties > 10%. Aside, we discussed in the previous section that we find a ∼ 30% of contamination in the region where both studies are complete. This is significantly larger than what we found in other clusters (i.e. the Pleiades and Ruprecht 147, see Sarro et al. 2014;Olivares et al. 2019, respectively) but this is expected given the lower galactic latitude and significantly smaller proper motions of IC 4665.
In addition, we see that the amount of members only recovered by DANCe increases as a function of magnitude in the region where both analyses are complete. We interpret this as a dependence of the contamination on the magnitude. The DANCe analysis does not use the parallax and thus it is expected that the photometry plays a major role, specially in this cluster with small proper motions.

Non-members in GDR2 nor in DANCe
All the remaining sources are classified as field stars by both studies. Most of them have extremely low membership probabilities which clearly identifies them as field population. There are several sources which have rather high probabilities but fall below the threshold. This means we can not definitely discard them as members and could be considered as candidate members depending on the scientific case. The sources which are spread along the rest of the diagram may suffer from the problems already discussed or simply, the observables in the two catalogues are too different. To clarify the membership of the uncertain cases we would need either a longer temporal base-line to improve the proper motions or spectroscopy to study their properties (i.e. low gravity due to youth).
In short, we see that in general the two independent analysis agree rather well, specially in the magnitude range were both are expected to perform well. The members obtained with both catalogues occupy the same space in the vector point diagram (see Fig. 6 left). The members coming from the DANCe catalogue typically have a larger dispersion and larger uncertainties, expected by the different precision of both catalogues. In the space of parallaxes (Fig. 6 right) we see that the members from the GDR2 analysis are very well concentrated around the median value (2.84 mas with a standard deviation of 0.36 mas). Some members from the DANCe analysis, which do not use this parameter for the classification, also have parallaxes compatible with the cluster distribution. Others have parallaxes incompatible with the cluster (at 3σ level) and they are either problematic measurements (because they are very faint) or contaminants. Future releases of the Gaia catalogue will help to clarify these cases. When we introduced the GDR2 catalogue in Sect. 2.1, we mentioned that we have not filtered the data in any manner in order to be the most complete possible. Here, we discuss the RUWE goodness of fit indicator of the members found in this study. Our sample contains sources with a RUWE in the range 0.8-16.3, and only 9% of them have a RUWE larger than the recommended threshold (1.40). However, we insist that all the sources with a RUWE larger than the recommended one do not always have a wrong solution and future releases of Gaia or complementary observations will tell.

Comparison with other studies
In this section we compared our list of 819 members with other studies in the literature and found that 409 (50%) are new members. We have cross-matched each of the members lists reported in the literature with ours using a maximum separation of 1 .
In Figure 7 we compare the members we find with two of the most representative membership studies of IC 4665: the one from Lodieu et al. (2011) based mainly on photometry, and the one from Cantat-Gaudin et al. (2018) using the Gaia DR2 astrometry. As a general trend, purely photometric studies tend to have more contamination than spectroscopic ones or the ones based on Gaia DR2 astrometry.

de Wit et al. (2006)
These authors photometrically selected 691 low-mass stellar and 94 brown dwarf candidate members over an area of 3.82 square  degrees centred on the cluster. In addition, they applied a filter for bright stars based on the proper motions from Tycho-2 and UCAC2 public catalogues.
We detected some astrometric offsets between their positions and ours and consequently extended the cross-match search radius to 2 . We confirmed 195 of their members and rejected the rest of their candidates which have very low membership probabilities in our study. Therefore, we estimate a contamination up to 75% in their study compatible with their own estimate. We believe that one of the reasons of their large contamination is a problematic photometric calibration (their i and z bands photometry display an offset of ∼ 1 mag compared to Pan-STARRS).

Manzi et al. (2008)
These authors did not attempt to do a comprehensive census of the cluster. Instead, they photometrically selected candidates from the literature and then spectroscopically confirmed 37 of them. Their aim was to determine the age of IC 4665 using the lithium depletion boundary method. We confirmed 29 of their members (78%) and discarded the remaining 8 (two of them were classified as not fully secure members by the authors). These 8 members are discarded because their Gaia DR2 parallaxes and/or proper motions are far from the cluster distribution, although their photometry falls on the cluster sequence. Therefore, these sources are either interlopers or have a problematic astrometric solution in Gaia DR2.

Jeffries et al. (2009)
These authors aimed to study the pre-main-sequence lithium depletion for low-mass stars in IC 4665. For this purpose, they selected 40 members according to several spectroscopic criteria. We confirmed 30 of their members (75%) and rejected the remaining 10. This study has 12 members in common with Manzi et al. (2008) and from those, only one source is rejected by our study. Again, the 10 members are rejected by our analysis because their Gaia DR2 astrometry is incompatible with that of the cluster, but the same reasoning of Manzi et al. (2008) applies.

Cargile & James (2010)
These authors used a photometrically selected sample of members in the central region of the cluster (1 square degree) to study the age and distance of IC 4665. Their sample contained 382 candidates members, 49 of which are confirmed by our study. From this, we estimated their contamination to be 87%.

Lodieu et al. (2011)
These authors used photometry from UKIDSS and CFHT to identify members in IC 4665. They presented a sample of 1 372 members in the magnitude range 15 < i < 20.4 which they used to study the luminosity and mass functions.
Only 240 of their candidates (17%) were classified as members in our work (see Fig. 7). The majority of the rejected candidates have extremely low membership probabilities in our analysis (both in GDR2 and DANCe). We believe the reason of their large contamination (∼80%) is the same photometric offset as for de Wit et al. (2006) since they used the same data.
These works constituted the most exhaustive study, specially regarding the low-mass regime, previous to the results we present here. Given the high levels of contamination found by the present analysis, we hereafter do not attempt any comparison of their luminosity and mass function. 4.2.6. Bravi et al. (2018) These authors used the Gaia ESO Survey to study the IC 4665 open cluster. They carried out spectrosopic observations of 567 sources in the region of the cluster. They used spectroscopic criteria to exclude obvious contaminants and then, they computed membership probabilities using the radial velocity distribution of the cluster and the field. They ended up with 29 sources with membership probability > 0.5, 24 of which greater then > 0.8. From these 29 candidates, 20 were confirmed by our study (15 have probabilities > 0.8 according to their study), and the remaining 9 were definitely rejected in our study. As for the previous spectroscopic surveys, the reason we discard these 9 members is the Gaia DR2 astrometry, which is incompatible with that of the cluster.

Gaia Collaboration et al. (2018a)
With the purpose of demonstrating the power of Gaia DR2 in highlighting the fine structures of the Hertzsprung-Russell diagram, these authors selected members for a number of open clusters. Their ambitious goal requested to select only the sources with the highest precision in astrometry and photometry, and among other filters, they restricted to sources brighter than G = 18. One of the clusters of their study is IC 4665, for which they provided a list of 174 members based only on the astrometric solution of Gaia DR2. They claimed their list not to be complete but to contain potential members, i.e. to have an extremely low contamination rate.
To do a fair comparison with this study, we only considered our members which are in the same magnitude and spatial range (brighter than G = 18 and 2.4 • radii around the centre of the cluster). This results in 267 members, 215 of which are classified as members by our analysis with the GDR2 catalogue and the rest come from the analysis with the DANCe catalogue only.
The study of Gaia Collaboration et al. (2018a) and ours have 162 members in common which is a 93% of their list. From the 12 objects classified as members by these authors and not by our study, there are 4 which have probabilities > 0.5 but fall below the optimum threshold we adopted, and 8 which have lower probabilities. From these 8 members, only 1 was also classified as member by a similar study (Cantat-Gaudin et al. 2018). These small differences are part of the Poissonian noise of the membership analysis. In addition, we find 53 members not detected by these authors which are spread along all the parameter spaces (proper motions, parallax, and magnitude), following the cluster distribution. Some of these 53 members could have been discarded by the authors in their data filtering.

Cantat-Gaudin et al. (2018)
These authors provided a membership analysis for a large number of clusters making use of the recent Gaia DR2 data. In order to avoid large uncertainties, they restricted to the sources brighter than G = 18. They used an unsupervised membership algorithm to derive membership probabilities using only the astrometry of Gaia DR2, and they found 175 members of IC 4665.
To do a fair comparison with this study, we only considered our members which are brighter than G = 18 and occupy the same spatial region of the sky (∼ 2 • radius around the centre of the cluster). This results in 244 potential members, 205 of which come from the analysis of the GDR2 catalogue and the rest only from the analysis of the DANCe catalogue.
The study of Cantat-Gaudin et al. (2018) and ours have 170 members in common which is a 97% of their list. From the 5 objects classified as members by these authors and not by our study, there are 4 which have probabilities > 0.5 but fall below the optimum threshold we adopted, and 1 which has a probability of 0.2. Again, these small differences are part of the uncertainties of the membership analysis. In addition, we find 35 members not detected by Cantat-Gaudin et al. (2018), 17 of which are also classified as members by Gaia Collaboration et al. (2018a). These members are randomly distributed within the proper motions and parallax distributions. We found 3 very bright members of G magnitudes 7.5, 9.5 and 10.5 and the rest are fainter than G = 14.5. The DANCe members not classified by GDR2 in this magnitude range (39 sources with 14.5 < G < 18) are likely to be contaminants, as discussed in Sect. 3.

Empirical and theoretical isochrones
In this section we provide the empirical isochrones obtained with our members of IC 4665. Then we compare them with theoretical evolutionary models in several apparent colour magnitude diagrams. Finally, we convert the apparent colour magnitude diagrams to absolute colour magnitude diagrams and discuss them.

Empirical isochrones
The empirical isochrones provide a key information to compare and constrain the theoretical evolutionary models. In this study, we use the membership analysis of IC 4665 to report the empirical isochrone of a 30 Myr old open cluster (see Fig. 8).
To obtain the empirical isochrones, we fitted a principal curve to the members in several apparent colour magnitude diagrams. Then, we manually shifted the principal curve to reach the lower edge of the distribution which is supposed to correspond to the single star zero age main sequence (ZAMS). In addition, we applied manual offsets where needed to better fit the lower edge of the cluster sequence. The empirical isochrones we provide are thus the lower envelope sequence of the members. This does not correspond to the principal curve which indicates the mean position of the sequence.
In Tables B.1 and B.2 we give the apparent magnitudes for the IC 4665 empirical sequence in the G, G BP , G RP , i, Y, J, H, K s over the dynamic range of our dataset. We reiterate that they are just an estimation of the ZAMS locus at the age of IC 4665 which could be used to compare with other clusters. For more quantitative analysis, we provide the full list of members and let the interested user decide the most convenient way to use them.

Evolutionary models
In Figure 9 we compare the observed sequence of IC 4665 and the empirical isochrones described in Sect. 5.1 to the 30 Myr models of PARSEC-COLIBRI 3 (Marigo et al. 2017) and BT-Settl 4 (Allard 2014) in several colour magnitude diagrams. The distance modulus applied to the models uses the median parallax of the GDR2 members (2.81 mas). We did not include the effect of any systematic because here we only intend to do a qualitative comparison. We also corrected the models for the median extinction measured with Gaia (A G = 0.62 mag which corresponds to A V = 0.72 mag using the A G /A V = 0.85926 from the PARSEC website, see footnote 3).
As a general result, we see that the models show a major improvement with respect to previous versions, specially in the y, J, H, K s bands (see e.g. the comparison of the Pleiades by Bouy et al. 2015) even at such a young age. The brightest stars are only covered by the PARSEC isochrones while the faintest are only covered by the BT-Settl models. Between i = 11 − 15 mag, both models agree fairly well between them and with the observations. However, the PARSEC models start to differ from the observations at i > 15 mag, and in this magnitude range the BT-Settl models are believed to be more accurate. Despite the global improvement of the models in all the photometric bands, we still find a space for improvement in some of them, specially the ones involving the redder bands (see middle left and bottom left panels of Fig. 9). For this, low contaminated samples combined with accurate photometric measurements along a wide magnitude range are essential.
Regarding the Gaia DR2 photometry, it is noticeable that the G BP band shows a larger spread for magnitudes > 18 mag. In the near-infrared, our measurements come mostly from 2MASS that has relatively large errors beyond 14 mag explaining the larger dispersion between 14 < J < 17. Beyond, the measurements come from our own deeper images and both the uncertainties and the dispersion of the isochrone are significantly smaller.

The absolute colour magnitude diagram
To build the absolute colour magnitude diagram, we first converted individual apparent magnitudes to absolute magnitudes. This transformation requires the distance and extinction of each source. Given the different origin of our two catalogues (Gaia provides individual parallaxes and DANCe does not) we decided to follow two different approaches.
For the members obtained with the GDR2 membership analysis we used individual parallaxes to compute distances. We used the Kalkayotl 5 code as in Olivares et al. (2019) which performs a Bayesian probabilistic inference to compute posterior probability distributions for the distance of each member. We chose a Cauchy prior which is the recommended for clusters by the manual. The location of the prior was set to 350 pc (the approximate distance of the cluster), and the scale to 100 pc (in order to have a loose prior). In Figure 10 we show the individual posterior distance distribution of each of the GDR2 members. The median distance is 351 pc and the standard deviation is 55 pc.
We estimated the extinction in two independent ways. First, we used the Gaia extinction estimate (a_g_val) for the 88 members which have it available. This value is recommended not to be used individually but statistically on a set of stars (Andrae et al. 2018). Thus, we compute the median A V = 0.72 mag and the standard deviation 0.38 mag. Aside, we inferred the individual absorption of each source using a Bayesian model as in Olivares et al. 2019 (see Sect. 6.2). In this case, we also compute the median of all the individual maximum a posteriori probabilities (MAP) which is A V = 0.66 mag, compatible with the extinctions from Gaia.
To compute the absolute magnitude we sampled the apparent magnitude with a Gaussian centred at the observed magnitude and a standard deviation equal to the uncertainty. Then, each sample was converted to absolute magnitude by sampling the posterior distance distribution obtained with Kalkayotl. We added the median absorption of the cluster to each member. To compute absolute magnitudes for the DANCe members we followed a similar approach. The only difference is that instead of sampling the distance from the individual posterior distributions, we sampled the distance from the cluster distribution obtained with all the GDR2 members.
In Figure 11 we show the absolute colour magnitude diagram of IC 4665 where we have overplotted the PARSEC and the BT-Settl models. Thanks to the precision of the Gaia DR2 parallaxes, we find that the isochrone has broadened little with respect the the one observed in the apparent colour magnitude diagrams. In addition, we have included a mass scale. We have candidate members down to masses of ∼ 0.02 M , well within the sub-stellar regime. We see that the PARSEC models start to differ for masses < 0.7 M and in this low-mass regime the BT-Settl models represent better the observations. For this reason, to convert magnitudes to masses, we use the PARSEC models for the high mass stars and the BT-Settl models for low-mass stars (see Sect. 6.2).

The apparent magnitude distribution
The apparent magnitude distribution is a direct measurement of the number of sources observed at different brightnesses. The importance of this function lies in the fact that it does not depend on evolutionary models nor on distance estimates and thus, its validity does not expire (unless selection problems are present). The magnitude distribution of IC 4665 was obtained applying a Gaussian kernel density estimation (KDE) independently to the GDR2 and the DANCe members. The two samples were treated independently because of the different validity range of each catalogue. To estimate the optimal bandwidth of the KDE we considered the Scott's rule (Scott 1992) and the Silverman's rule (Silverman 1986), which gave similar results. The optimal bandwidth of the KDE was 0.3 mag (both in the GDR2 and DANCe members), and the uncertainties were estimated by means of a bootstrap with 100 repetitions. We have estimated the effect of contamination and completeness as a function of the magnitude using synthetic data as in Olivares et al. (2019). Given that the contamination rate estimated this way is less than 15%, we realised that when we correct for these two effects the magnitude distribution we obtain is compatible with the original one within the uncertainties. For this reason, we decided to work with the magnitude distribution we obtain directly from the observations. In Fig. 12 we show the magnitude distribution of IC 4665 in the G band for the GDR2 members and in the i band for the DANCe members. These functions are available in Tables B.3 and B.4, respectively.
At G ∼ 14.5 mag there is a flattening of the apparent magnitude distribution which corresponds to the Wielen dip (Wielen et al. 1983). This feature has been reported in other open clusters such as the Pleiades (Lee & Sung 1995;Belikov et al. 1998), Praesepe and Hyades (Lee et al. 1997), NGC 2516(Jeffries et al. 2001, NGC 2547 (Naylor et al. 2002), and Ruprecht 147 (Olivares et al. 2019). Kroupa et al. (1990) explained this dip as the result of a change in the opacities in the corresponding mass range.
The apparent magnitude distribution peaks at G = 18.2 mag for the GDR2 members and at i = 17.6 mag for the DANCe members which in both cases correspond to 0.2 M , according to the PARSEC and BT-Settl models and assuming an age of 30 Myr. This result is in agreement with what Bouy et al. (2015) found in the Pleiades.
A change of slope seems to happen around i ∼ 21 mag, which could indicate that different formation mechanisms are at work for ultracool objects in this mass range. This change of slope is nevertheless very close to our estimated limit of completeness and not statistically significant yet.

The present-day system mass function
We estimated the mass of each source using the Sakam 6 code (Olivares et al. 2019). This algorithm infers the posterior distribution of the mass together with the A V extinction, given the absolute photometry of each source (computed in Sect. 5.3) and a theoretical evolutionary model. The model does not include effects in the variability of the source due to binarity, activity or others. These effects are eventually included in the extinction estimate, enlarging its uncertainties. We used the PARSEC model to infer masses for the GDR2 members and the BT-Settl model for the DANCe members.
To compute the present-day system mass function, we took samples of the a posteriori distribution inferred with Sakam and applied a Gaussian KDE with a bandwith of 0.3 (in log 10 m where m is the mass, and the bandwidth of the KDE was estimated from the Scott's and Silverman's rule). We estimated the uncertainties from 100 bootstrap repetitions and reported the 1 and 3σ confidence levels. The completeness limits were propagated from the catalogue completeness in apparent magnitude (see Table 2). The mass function obtained with the DANCe analysis was renormalised so that the mass distribution functions had the same area in the region where both studies are complete, i.e. 0.15-0.8 M . Figure 13 shows the present-day system mass function of IC 4665 for the GDR2 (blue) and DANCe (green) members. We see that the two functions overlap reasonably well. There are some deviations (even inside the complete range) but they are smaller than 3σ. The robustness of our methodology, specially in the error propagation, results in a mass function with an accuracy significantly better than in the past (i.e. de Wit et al. 2006;Lodieu et al. 2011).
A number of noticeable details are present in the mass function. At 3 M we observe a feature which has not been reported in the literature before. It is not clear whether this is a real feature of the mass function or an artefact. Several sources of error could be responsible for that, in particular: 6 https://github.com/olivares-j/Sakam the uncertainties or errors of the transformation from apparent magnitudes to masses, since it is not observed in the magnitude distribution (Fig. 12), multiplicity: the Gaia DR2 catalogue excluded a number of binary stars. Since massive stars are more often in multiple systems than their lower mass counterparts (e.g. Lada 2006), we might be missing a larger fraction of massive members because of multiplicity, variability: massive stars can also display photometric variability, which is not included in our algorithm to determine individual masses. Slowly Pulsating variable stars appear at three solar masses and beyond. However, these are small amplitude variables (0.1 in V) and should not have a major impact in our selection.
Additionally, and as mentioned in Section 2.1, we assumed that the GDR2 catalogue is complete for G > 7 mag (∼ 5 M ) but a few sources between 7 < G < 12 mag (5 m 1.6 M ) could be missing (Gaia Collaboration et al. 2018b). Nevertheless, this would only increase the number of members in this range. The Wielen dip reported in the magnitude distribution (Fig. 12) would be expected around 0.75 M in the mass distribution but is not observed. If confirmed, this result would support the hypothesis of Kroupa et al. (1990) explaining this feature by a change in opacity rather than a change in the mass function. We nevertheless note that the Wielen dip may have been masked by the KDE bandwidth. Olivares et al. (2019) indeed reported a Wielen dip in their mass function between 0.6-0.8 M with a typical scale of ∆ log 10 m ∼ 0.13, smaller than our bandwidth of 0.3 (in log 10 m).
The function is rather flat between 0.1 and 1 M having a maximum at 0.28 M . For masses < 0.1 M the distribution drops. The change of slope at the very low mass end mentioned above is not visible in the mass function.
The highest mass object has a MAP estimate of 6.2 M and the lowest mass object has a MAP estimate of 13 M J according to the PARSEC and BT-Settl models respectively, and assuming an age of 30 Myr. To compute the brown dwarf-to-star ratio we sampled the posterior mass distribution of each member. Then, we used these samples to compute the ratio of brown dwarfs and stars within the completeness region of our sample (6-0.05 M ) and using a mass threshold of 0.08 M . We did a bootstrap over all the members with 100 repetitions and, we obtained a median ratio of 0.067 ± 0.005. This value is lower to what has been seen in other nearby young clusters as in IC 348 and Taurus (Scholz et al. 2012, and references therein). However these studies are complete down to lower masses (∼ 0.02M ).

Comparison to other clusters and theoretical models
In Figure 14, we compare the PDMF we obtained for the 30 Myr open cluster IC 4665 to the one from the Pleiades (120 Myr; Bouy et al. 2015). To facilitate the comparison, we normalised the mass function of IC 4665 over the whole mass range where it is complete. Then, we normalised the Pleiades mass function so that it had the same area between 0.05-0.6 M , range where both functions are complete. We see that in general both functions match fairly well within the uncertainties and the main differences are observed at the extremes of the distributions. For the high mass domain, IC 4665 has more massive stars than the Pleiades. In this range the number of members is quite small (only 12 objects have masses >3 M in IC 4665) leading to rather large statistical uncertainties. In addition, multiplicity (more frequent among high mass stars) and variability affect the luminosities and might contribute to the differences observed. Regarding the low-mass regime, we see that both functions are compatible (within 3σ uncertainties) down to the IC 4665 completeness limit (∼0.05 M ). Nonetheless, we observe that between 0.046 and 0.16 M the mass function of IC 4665 might display a slight over-density with respect to the Pleiades at the 1σ level only. For masses lower than 0.05 M , the Pleiades mass function exhibits a change of slope which the authors related to a different mechanism of star formation for this regime of masses. In the case of IC 4665, we do not detect this change of slope but this could be because it is beyond the completeness limit of the catalogue.
In Figure 14 we have overplotted the sytem IMF of two models, namely Chabrier (2005) and Thies et al. (2015), normalised in the same mass range as the mass function of IC 4665. In the high mass regime (> 1 M ), both models assume a power law IMF with Salpeter slope which is compatible within the uncertainties with the empirical mass function of IC 4665.
For intermediate and low masses, we see that the mass function of IC 4665 is compatible with the model of Chabrier (2005) between 0.1-1 M . For lower masses, the model predicts too many stars compared to our results. The model of Thies et al. (2015) is compatible with the empirical mass function between 0.2-1 M but between 0.05-0.2 M it also predicts too many stars. Below 0.05 M the model approaches the empirical mass function and beyond this limit our survey is not complete.

Projected spatial distribution
The spatial distribution of open clusters provides relevant information on the formation and early evolution of these systems. In Figure 15 we show the spatial distribution of the members of IC 4665 in galactic coordinates. At first glance we can intuit some structures which depart from a pure spherical symmetry (e.g. the cluster seems elongated towards the Galactic south). In this section, we apply a statistical treatment to quantitatively asses the probability that the structures we might see are significant.
We follow the same approach as Olivares et al. (2018) and we fit a series of parametric models to the projected spatial distribution of the cluster (i.e. in the plane of the sky). The algorithm they used, PyAspidistra 7 , computes the Bayesian evidence of each model and the posterior distribution of the parameters which characterise the model. Eventually, they compare the Bayesian evidence of each pair of models by means of the Bayes Factor. Here, we consider the same set of models as those authors: the Elson et al. 1987 (hereafter EFF), the Generalized Density Profile (hereafter GDP, also known as Nukker, Küpper et al. 2010 Fig. 12. Top: G magnitude distribution of IC 4665 with the members found with the GDR2 catalogue. Bottom: i magnitude distribution of IC 4665 with the members found with the DANCe catalogue. In both cases, the shaded regions indicate the 1 and 3σ uncertainties estimated from bootstrap (1σ dark and 3σ faint). The dashed lines indicate the region of incompleteness in each catalogue.    tion to infer the coordinates of the cluster centre, its ellipticity, and mass segregation. Using the equatorial coordinates (J2000), the median distance of the cluster (350 pc), and the J band we ran the PyAspidistra code and obtained the Bayesian evidence for each model and the Bayes Factor for each pair of models. The RGDP model is the one that shows the largest evidence in all the family models considered (spherical, elliptical and segregated). The family of models with largest evidence are the spherical models, as it is expected for such a young open cluster. However, the strength of this evidence is weak according to the criterion from Jeffreys (1961) so these results should be taken with care and we can not definitely discard the possibility of ellipticity or mass segregation. In addition, our results could be biased due to the contamination in the members and the size and shape of our initial catalogues as we discuss below.
The median parameters of each spherical model are reported in Table 6. The parameters α c and δ c correspond to the central coordinates of the cluster in RA and Dec, respectively. The core radius (r c ) is the unit scale of the density profile, therefore it differs for each model. The α, δ, and γ parameters correspond to the exponents of the different models. The tidal radius (r t ) is only defined for the family of King's models. We refer the interested reader to Olivares et al. (2018) for a detailed discussion of these parameters. We see that the centre of the cluster is well determined by all the models at RA = 266.6 • , Dec = 5.4 • . All the models also agree to a core radius of ∼ 2 pc and the small dispersion is expected since this parameter has a different interpretation in each model. Only the family of King models are defined in terms of a tidal radius. The median values of the tidal radius reported in Table 6 vary from one model to another and have extremely large uncertainties. The King model with larger evidence is the OGKing model which predicts a r t = 55 pc, 3 times larger than the radius analysed in this study (18 pc). This is probably the main reason why we fail to well constrain this parameter (see also the discussion on the main caveats at the end of the section). However, we note that up to the present, this study of IC 4665 is the one with largest radius. Our estimate of the tidal radius is remarkably larger than previous values (e.g. de Wit et al. 2006 reported a tidal radius of 1 • corresponding to ∼ 6 pc at the distance of the cluster; however, this value results from a highly contaminated sample).
Here we enumerate some of the caveats and limitations of our study of the spatial distribution.
-Our members come from a catalogue which was circularly selected and the ends of the catalogue can clearly be seen in Fig. 15. This can bias our results to favour a circular over an elliptic model. -The spatial coverage of the DANCe catalogue (see Sect. 2.2 and Bouy et al. 2013) implies that the faintest members are more likely found in the centre. This region is where we have not only the highest number of images, but also the deepest ones and the ones with the longest time baseline. As a consequence, the proper motions also have in general better precision in this area, which in turn has an influence on the membership probabilities. This can have an impact both in the study of the shape (circular or elliptical) and the study on the segregation. -The tidal radius and the contamination rate are degenerate: a large contamination rate increases the density of members and the models need a larger tidal radius to explain the observations. This is specially critical at the outskirts, where we believe we might have a larger contamination. As we have discussed in Sect. 3.6 and Sect. 4, the majority of our contaminants come from the DANCe catalogue. -Our members come from a catalogue which was truncated to a radius of ∼18 pc, similar to the expected tidal radius. It is highly difficult to estimate the tidal radius without having information beyond it. All this difficulties are reflected on the Bayesian evidence of each model: the family of King's models have lower evidences. The truncation in radius could also bias the study on the elipticity and segregation of the cluster.

Conclusions
We presented an exhaustive study of the properties of the IC 4665 open cluster. We combined the recent Gaia DR2 data with the deep, on ground observations of the COSMIC DANCe project to search for members. We used the same methodology as Olivares et al. (2019) to derive Bayesian posterior probabilities for all the sources and found 819 members, 50% of which are new members. Our members have magnitudes in the range 7 < J < 19.4, which correspond to masses between 6.2 M -13 M J , according to the PARSEC and BT-Settl evolutionary models, and assuming an age of 30 Myr. Using this sample, we provided the empirical isochrones of the cluster, an estimate for the distance, the magnitude distribution, the present-day system mass function with unprecedented accuracy for this cluster, and a study of the spatial distribution. Comparing our members with previous studies in the literature, we found that most of the previous studies were based on highly contaminated (up to 80%) or incomplete samples. The low motion of this cluster with respect to the field (< 10 mas yr −1 ) complicates the membership analysis. For this reason we found a larger contamination rate in this study compared to others which use the same methodology applied to clusters with larger proper motions (Sarro et al. 2014;Olivares et al. 2019). Using synthetic data, we estimated a CR of 10% for GDR2 and a 13% for DANCe. Comparing the two studies, we estimate that the DANCe contamination reaches up to 30% in the region of completeness. The main reason of this underestimated CR in the DANCe study is the lack of parallaxes. Anyway, this study provides the most accurate membership analysis up to date by far and thus, offers the possibility to revisit other fundamental parameters such as the age.
We found that the PDMF of IC 4665 in the intermediate mass range (0.1-1 M ) is comparable to that of the Pleiades (Bouy et al. 2015) and to models of the IMF (Chabrier 2005;Thies et al. 2015). For higher masses, the observations have a steeper slope than that of the models (Salpeter slope). In addition, for the case of IC 4665, we find a plateau at 3 M which if further confirmed would represent a new feature of the mass function. In the mass range 0.05-0.2 M the models predict too many low-mass stars. For masses lower than 0.05 M the Pleiades have more members than in IC 4665, but at this mass regime our study is not complete so we might be missing members.
Combining our comprehensive census of the cluster with the Gaia DR2 parallaxes we estimated the distance of the cluster to be of 350 pc, this value being similar to what other studies recently derived (Gaia Collaboration et al. 2018a). We found that the best surface density profile for IC 4665 is the RGDP model with a core radius of 2 pc, and sperical symmetry. However, we can not definitely discard the possibility of ellipticity nor mass segregation. In the future, we aim to include a study on the velocity distribution which would allow us to characterise the kinematic and dynamic state of the cluster in the 6D space phase.
Acknowledgements. We thank the referee for his/her thorough review and highly appreciate the comments and suggestions, which significantly contributed to improving the quality of the publication. We are grateful to Ingo Thies for sharing with us his models of the Initial Mass Function, to France Allard and Isabelle Barraffe for providing us with the latest version of the BT-Settl models, and to Paul Price for his help with the LSST/HSC pipeline. This research has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 682903, P.I. H. Bouy), and from the French State in the framework of the "Investments for the future" Program, IdEx Bordeaux, reference ANR-10-IDEX-03-02 . This Project has been funded by the Span- . This research draws upon data distributed by the NOAO Science Archive. NOAO is operated by the Association of Universities for Research in Astronomy (AURA) under cooperative agreement with the National Science Foundation. 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 work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa. int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the VizieR and Aladin images and catalogue access tools and of the SIMBAD database, operated at the CDS, Strasbourg, France. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This research has made use of GNU Parallel (Tange 2011), Astropy (Astropy Collaboration et al. 2013), Topcat (Taylor 2005), STILTS (Taylor 2006).