© The Authors 2022

1. Introduction

The IC1396 H II region (also known as S131; Sharpless 1959) is part of the large, star-forming Cepheus bubble (Patel et al. 1998), ionized by the multiple system (O5+O9) HD206267 (Peter et al. 2012; Maíz Apellániz & Barbá 2020). IC1396 contains the young cluster Tr 37 (Platais et al. 1998; Patel et al. 1995, 1998), which has been used to estimate the distance to the region, ranging from 870 pc from main sequence (MS) fitting (Contreras et al. 2002) to 945${}_{-73}^{+90}$ pc based on Gaia DR2 (Sicilia-Aguilar et al. 2019). Hundreds of known IC1396 members, mainly young stellar objects (YSOs), have been identified in the region using various techniques. These detections have been based on spectral features that are characteristic of youth, such as Li I absorption or Hα emission (Kun 1986; Contreras et al. 2002; Sicilia-Aguilar et al. 2004, 2005, 2013; Barentsen et al. 2011; Nakano et al. 2012), X-ray emission (Mercer et al. 2009; Getman et al. 2012), near-infrared variability (Meng et al. 2019), and infrared (Rebull et al. 2013), as well as detections of infrared (IR) excess signaling the presence of a disk (Reach et al. 2004; Sicilia-Aguilar et al. 2006a; Morales-Calderón et al. 2009).

The IC1396 region also contains some bright-rimmed clouds (BRCs) shaped by the UV radiation from HD206267 (Patel et al. 1995; Barentsen et al. 2011). These BRCs display evidence of ongoing star formation and are excellent laboratories for studying young stars in different evolutionary stages. The BRCs in the region display a range of velocities (V_LSR = +5 to −9 km s⁻¹) derived from CO molecular maps (Patel et al. 1995). The largest BRCs are IC1396A, IC1396B, and IC1396N (labeled as “E” by Pottasch 1956; Patel et al. 1995), all of which contain several tens of solar masses in gas and have different velocities (e.g., V_LSR = −7.9 for IC1396A, −5.4 km s⁻¹ for IC1396B, and 0.6 km s⁻¹ for IC1396N; Patel et al. 1995). There is also a velocity difference between IC1396A and Tr37 (Sicilia-Aguilar et al. 2006b) but, despite of this, the initial Gaia DR2 results revealed no significant discrepancy between the proper motion of the sources physically related to IC1396A and to Tr37 (Sicilia-Aguilar et al. 2019).

The largest globule, IC1396A, contains a moderate population of low-and intermediate-mass stars ranging from a Class 0 source (Sicilia-Aguilar et al. 2014, 2019) to many Class I, II, and III members (Reach et al. 2004; Sicilia-Aguilar et al. 2006a; Getman et al. 2007; Barentsen et al. 2011). This has been interpreted as a signature of sequential or triggered star formation (Sugitani et al. 1991; Sicilia-Aguilar et al. 2005, 2006a, 2019; Getman et al. 2007, 2009) related to the expansion of the H II region and its interaction with the surrounding cloud. IC1396N also displays recent star formation activity, featuring many Herbig Haro (HH) objects (Nisini et al. 2001), a large number of H₂ jets and knots (Caratti o Garatti et al. 2006; Beltrán et al. 2009), and embedded YSOs inside the globule, such as the sources BIMA 1, 2, and 3 – two of which are driving molecular outflows (Beltrán et al. 2009).

Gaia data (Gaia Collaboration 2016, 2021) now allow us to explore the astrometry of the known cluster members, which can then be used to identify new members (Lindegren et al. 2000; Franciosini et al. 2018; Roccatagliata et al. 2018, 2020). Gaia maps the structure of the region in a multidimensional space, including proper motion and parallax. Using Gaia has the advantage that it can detect objects for which the youth indicators may be less clear (e.g., stars without disks and no accretion and intermediate-mass stars) and can be used to distinguish populations by their astrometry rather than by their location, which is always uncertain due to projection effects.

In this paper, we use Gaia EDR3 data to provide a new, independent view of IC1396, including the IC1396A, IC1396B, and IC1396N globules. We extract the astrometric properties of the known members, using them to identify new ones. We also explore the properties of the YSOs, including disk fraction and spatial distribution, and combined 2MASS and WISE data to characterize their disks. Narrow-band [S II] and H₂ photometry are used to trace the presence of shocks, jets, and outflows produced by young embedded members. Finally, we discuss the advantages and limitations of each method and the improvements of Gaia EDR3 concerning DR2. In Sect. 2, we describe the observations and data reduction. New spectroscopic members are identified in Sect. 3. In Sect. 4, we study the new members revealed by Gaia and the kinematics of the region. Section 5 shows a final analysis of the region properties. The results are summarized in Sect. 6.

2. Observations and data reduction

We used ground-based spectroscopy, photometry, and Gaia EDR3 data to study star formation in IC1396. We obtained near-infrared broad- and narrow-band imaging to study the star-cloud interaction and optical spectra to confirm and classify objects that were young star candidates. The details of the observations are given in the following sections.

2.1. Optical and NIR imaging

Since the largest BRCs, namely, IC1396A, IC1396B, and IC1396N, have signs of ongoing star formation, we obtained deep, broad, and narrow-band IR photometry and narrow-band optical observations at the Calar Alto Observatory (Centro Astronómico Hispano en Andaluía (CAHA), Spain). They included data from proposals F-17-2.2-0.19, F18-3.5-007, and F17-3.5-018 (PI M. Pelayo-Baldárrago) and H14-3.5-002 (PI A. Sicilia-Aguilar). The observations are summarized in Table 1 and the fields covered are shown in Fig. 1. We also include data for IC1396A from Sicilia-Aguilar et al. (2013) to round out the discussion.

Fig. 1. WISE W4 (22.19 μm) image of the IC1396 H II region. The BRCs are labeled with their names. The red boxes and yellow circles enclose the fields observed by OMEGA2000 and CAFOS, respectively. Note: CAFOS FOV has been reduced to account for the fact that the Fabry-Perot etalon has only a small area where the wavelength is accurate. The HD206267 star that is ionizing the region is indicated with a yellow star.

We performed JHKs broad-band near-IR photometry using the OMEGA2000 near-infrared wide-field camera (Röser et al. 2008)¹ on the 3.5-m telescope. Its field-of-view (FOV) is large enough to fit each globule within one pointing, achieving a high spatial resolution, sensitivity, and dynamical range (∼9–20 mag). We followed the OMEGA2000 pattern of 3 second integrations, combining 20 coadds, within a 25-position dither pattern, for a total exposure time of 25 min in the center of the image. The data reduction, including dark and flat subtraction, was done using standard Python Astropy routines (Astropy Collaboration 2013, 2018). The sky emission was subtracted with Python using sky frames created from the science images. After reduction, images were aligned and combined using AstroImageJ (Collins et al. 2017) and the coordinates were derived using Astrometry.net² (Lang et al. 2010).

We also obtained narrow-band imaging at 2.122 μm (H₂) and 2.160 μm (Brγ) and continuum at 2.144 μm for all three globules to explore the presence of shocks or jets from the embedded population. We obtained 25 dithered positions with 30s exposures and 5 coadds, with a total exposure of 3750s in the center of the image. The sky emission subtraction was taken from the sky frames created from science images. The data reduction was carried out with Python, following the same procedure as for JHKs described above.

To complete the study of feedback in the BRCs, we also obtained narrow-band [S II] images for IC1396N and IC1396B using the Fabry-Pérot interferometer with the Calar Alto Faint Object Spectrograph (CAFOS) camera³ at the 2.2 m telescope and following the same setup and procedure used by Sicilia-Aguilar et al. (2013) for IC1396A. The etalon was tuned for a narrow bandwidth of 10 Å around 6717 and 6730 Å (for the [S II] lines) plus an adjacent, line-free continuum at 6750 Å.

We can only use the central part of the large CAFOS field because the wavelength changes across the Fabry-Pérot FOV. We thus obtained three pointings, one for the IC1396N globule and two for the more extended IC1396B globule. In each case, we obtained three dithered 600s exposures centered on the two [S II] lines and the narrow continuum. Data reduction, including bias and flat fielding, was done with Python (Astropy). We combined all science images with AstroImageJ and assigned the coordinates through Astrometry.net. Finally, we combined both [S II] images at 6717 Å and 6730 Å to get a better signal-to-noise ratio (S/N) and we subtracted the continuum at 6750 Å to remove the stellar contribution.

2.2. MMT/Hectospec Spectroscopy

We obtained the first set of spectra with the multifiber spectrograph Hectospec (Fabricant et al. 2005) mounted on the 6.5m Multiple Mirror Telescope (MMT). Hectospec has 300 fibers to assign to individual positions, distributed over a one-degree field of view. It covers a wavelength range from 4000 Å to 9000 Å with resolution R ∼ 3000. Our observations include data from proposals 2014b-UAO-S17S and 2015b-UAO-S300 (PI: S. Kim; see Table 2). We obtained eight fiber configurations, each of which had between ∼84 and ∼161 fibers for stellar objects and between ∼40 and ∼80 fibers for sky positions.

The data were reduced using IRAF routines following standard procedures (Fang et al. 2013). These include flat-fielding and extracting the spectra using dome flats and the IRAF specred package (dofibers task). The wavelength solution was extracted using the HeNeAr comparison spectra with the IRAF identify and reidentify tasks and calibrated using IRAF dispcor task. For each one of the eight configurations, there were several scientific exposures with exposure times in the range of 450–1500 s repeated between three and four times, depending on the magnitude of the stars in the field. Combined those belonging to the same field for the final analysis and sky subtraction (see Table 2).

The sky was subtracted using IRAF sarith, and following Sicilia-Aguilar et al. (2005, 2013). We selected the sky spectra in each setting and combined them using IRAF scombine. Because the H II region background is highly variable, we created an average sky spectrum, plus three further templates combining sky spectra visually classified as having bright, intermediate, or faint nebular emission for each configuration. Each sky template resulted from combining about 40–80 spectra for average emission or close to 10 for the bright, faint, and medium sky. These low-noise sky templates did a good job for most sources, although in some cases, the sky spectra had to be further re-scaled to remove the H II lines completely. The best result was selected by visual inspection, checking that the forbidden [N II] lines (λλ 6548 Å, 6583 Å), which normally arise from the H II region, and the telluric 5577 Å line, disappear. The final outcome of this process was a total of ∼900 source spectra with significant S/N, which were later examined (Sect. 3).

2.3. GTC / MOS and Long-Slit Spectroscopy

The second set of spectra was obtained with the Gran Telescopio CANARIAS (GTC, Rodríguez Espinosa et al. 1998) using the Optical System for Imaging and low Resolution Integrated Spectroscopy (OSIRIS, Cepa et al. 2003) instrument for doing both long-slit and Multiple Object Spectroscopy (MOS). Our GTC/MOS observations include data from proposals GTC35-14A (PI: M. Fang) and GTC/long-slit observations from proposals GTC55-12A and GTC30-12B (PI: D. García Álvarez) as part of a campaign to study unknown variable stars.

The GTC long-slit observations were taken in May 2013, and the GTC/MOS were obtained in July 2014. The wavelength range is ∼5000–10 000 Å, and the resolution is R ∼ 2000 for MOS and R ∼ 1900 for the long-slit. For the MOS spectra, we first cut the data into individual slits. After being extracted, all data reduction was done in the same way for MOS and long-slit spectra, using standard IRAF routines. In total, we extracted 63 spectra from MOS and 18 spectra from the long-slit. Their analysis is presented together with the rest of the spectral data in Sect. 3.

2.4. Gaia EDR3 data

We studied the structure of IC1396 and young star membership using EDR3 data (Gaia Collaboration 2021). We selected data within a 2° in radius, centered on HD206267, which contains nearly 1.2 million sources⁴. Gaia data is particularly important for the intermediate-mass stars since usual surveys (targeting youth spectral lines features, X-ray, and infrared excesses) have a certain bias against them, which led to them being poorly represented when exploring the initial mass function (IMF) of the region (e.g., Sicilia-Aguilar et al. 2005). We analyzed and compared the Gaia EDR3 and DR2 releases to estimate the differences, focusing the final conclusions on EDR3.

The membership identification is based on three astrometric parameters: parallaxes (ϖ) and proper motions (μ_α, μ_δ). We are limited by the reliability of the astrometric solution for faint targets, selecting only those with good-quality astrometry. Gaia EDR3 is complete between G = 12 and G = 17 mag and has a limit of G = 20.7 mag (Gaia Collaboration 2021), having improved slightly for stars fainter than G = 18 mag compared to Gaia DR2 (Fabricius et al. 2021). The magnitude limit for reliable data in dense areas is ∼18 mag (Luri et al. 2018). We set a limit on the fractional parallax error of < 0.25 (Bailer-Jones 2015), so the effective completeness limit drops to G = 17 mag. The quality and reliability of the Gaia astrometric data also depend on the Re-normalized Unit Weight Error (RUWE). The RUWE is a statistical indicator obtained from the Unit Weight Error (UWE) or astrometric chi-square that requires a renormalization depending on the magnitude and color of the source (Lindegren 2018, see their Eq. (2)). We only use sources with RUWE < 1.4, which have a reliable and consistent astrometric solution (Roccatagliata et al. 2020). Only about 38% of the known members have Gaia data that are good enough for our analysis. This introduces a bias against detecting new low-mass members but covers the intermediate-mass members that have been less studied.

We also used Gaia photometry (Riello et al. 2021) and the color-magnitude diagrams (Gaia Collaboration 2018) and isochrones to constrain the astrometric results and remove further polluting sources. The isochrones used throughout this work are those from the Padova and Trieste Stellar Evolution Code (PARSEC) release v1.2S (Bressan et al. 2012). They were obtained from the CMD 3.6 interactive service from the Osservatorio Astronomico di Padova⁵; using the EDR3 photometric system (filters G, G_BP, G_RP; Riello et al. 2021).

2.5. Ancillary data

To study the parallax and proper motion structure of the region, we compiled a collection of objects labeled as cluster members in the literature (Marschall & van Altena 1987; Kun & Pasztor 1990; Contreras et al. 2002; Mercer et al. 2009; Morales-Calderón et al. 2009; Sicilia-Aguilar et al. 2004, 2005, 2006a; Sicilia-Aguilar et al. 2013, 2019; Barentsen et al. 2011; Nakano et al. 2012; Getman et al. 2012; Rebull et al. 2013; Meng et al. 2019). From the original 1544 cluster members, we rejected 8 sources (see Table 3) that were marked as low-probability members in the literature and/or were in anomalous locations in the parallax or proper motion space. We added to this list our two newly confirmed spectroscopic candidates (see Sect. 3), summing 1538 known members in total (Table A.1). We summarize the literature member search results in Table 4. We also used 2MASS JHKs data (Cutri et al. 2003) to select potential new variable members or characterize the properties of members and mid-IR observations from Wide-field Infrared Survey Explorer (WISE, Cutri et al. 2013), to search for excess emission that could evidence protoplanetary disks. Both datasets were matched to the known member list using 0.6″ and 0.7″ radii, respectively.

3. Identification of new spectroscopic members

We use the MMT and GTC spectra to identify new members. From an initial collection of nearly 900 spectra, we filtered 121 targets with high S/N that could be potential members. Those are analyzed in more detail below, using four independent membership criteria combined to obtain the final membership. The results derived from each criterion and the final membership for each object are given in Table B.1.

3.1. Characteristic spectral lines of young stars

The most robust youth indicator is the Li I absorption line at 6708 Å (White & Basri 2003), but this criterion applies only to fully convective stars and the line may be hard to confirm in low S/N spectra. Other emission lines indicative of youth are the CaII IR triplet, He I, and strong Hα (Hamann & Persson 1992). We measured the equivalent widths (EW) of the lines and noise levels using IRAF splot.

We detected the Li I absorption line and/or Ca II triplet emission lines at > 3σ level in 52 sources. For fiber spectra, subtracting the Hα emission line produced by the H II region is problematic and may not always be accurate (more details in Sicilia-Aguilar et al. 2006b), so that weak Hα emission may not be of stellar origin and requires careful interpretation. Therefore, we classify as probable members those objects with large equivalent widths, EW(Hα) > 10 Å, plus at least one more membership criterion. This adds up to 49 objects. The 20 objects with a weak or narrow Hα emission, without Li I absorption or Ca II triplet emission lines and no other evidence of membership in the literature (no X-ray emission, no IR excess) are classified as probable non-members. Some spectra examples are shown in Fig. 2. We give a higher priority to this criterion for the final membership classification, so the 52 stars with clear indicators of youth are labeled as confirmed members.

Fig. 2. MMT (Hectospec) spectra of two confirmed members. The Hα emission line is cut off at the top for better visualization. The fluxes are given in arbitrary units.

3.2. Spectral types and extinction

We obtained the spectral type of all 121 candidates following Fang et al. (2020). We fitted the normalized spectra with two different sets of templates, including observations of pre-main-sequence stars from X-Shooter/VLT (covering a spectral type range from G6 to M9.5, see Fang et al. 2021), supplemented by a set of combination of stellar atmosphere models from Mészáros et al. (2012) and Husser et al. (2013). The spectral templates were artificially veiled and reddened, using the visual extinction law by Cardelli et al. (1989) and assuming R_V = 3.1. The spectral fit was done considering two cases, with and without veiling. The reduced χ² is used to estimate the best-fit combination of spectral type and veiling parameter.

We did not use the extinction from the spectral fitting (as pointed out in Fang et al. 2020) and instead recalculated the extinction (A₀) for the 121 sources using their spectral types and the Gaia colors, which gives more accurate results. The extinction of the new members is expected to be over a lower limit based on the cluster distance and to follow the extinction distribution of confirmed members (Contreras et al. 2002). Objects outwith the extinction distribution of spectroscopically confirmed members were considered as probable non-members (see below).

We calculated the interstellar extinction, A₀ (at λ = 550 nm, Danielski et al. 2018), using Gaia photometry and A₀ = 3.1 * E(B − V), (Gaia Collaboration 2018), so that A₀ = A_V for the mean extinction law (Cardelli et al. 1989). Among 121 candidates, the analysis was done only for 119 candidates because two sources lacked complete Gaia EDR3 counterparts. These two sources were labeled as having uncertain extinction (Table B.1). We use the effective temperature versus spectral type relation from Kenyon & Hartmann (1995) and the 4 Myr PARSEC isochrone (since the average age is ∼4 Myr, Sicilia-Aguilar et al. 2005; Getman et al. 2012, which we confirm later in Sect. 5.1) to obtain the theoretical magnitudes and colors for each spectral type. These results are used to calculate the color excesses.

E(G_BP − G_RP) and E(G − G_RP). Finally, the interstellar extinction A₀ is estimated using the relations from Gaia Collaboration (2018) and Danielski et al. (2018). The final value is obtained as the average of the A₀ values derived from the two color excesses for each source and the uncertainty results from their standard deviation. The values of A₀ are consistent for G <  16 mag, although they show a typical systematic deviation of 0.2 mag. For sources with G >  16 mag, the A₀ values are uncertain due to systematic effects on the Gaia EDR3 band photometry (Riello et al. 2021). Moreover, there is a trend for the flux in the BP band to be overestimated in faint red sources (Gaia Collaboration 2021; Riello et al. 2021; Fabricius et al. 2021), which appear much bluer in (G_BP − G_RP) than they should be. Therefore, their A₀ values need to be regarded with care. The correction of these systematic effects in the Gaia bands photometry is beyond the scope of this paper.

Fig. 3. Distribution of the extinction for the spectroscopically confirmed members (top). The mean error (0.3 mag) is shown with a red bar. Color-magnitude diagram for candidates with spectroscopic signatures of youth (bottom). Members (M) are marked by blue circles. Probable members (PM) are marked by red triangles. Probable non-members (PNM) are marked by gray circles. The number of objects in each class is also given. The magnitudes of each object are corrected by their individual extinction (full color symbols), except for objects with non-physical or uncertain extinctions, which have been corrected by cluster average extinction of 1.4 mag (open symbols). The PARSEC isochrones are corrected by the average distance of 925 pc.

We analyzed the extinction distribution for 52 objects with clear Li I and/or Ca II detection (spectroscopically confirmed members, Fig. 3) to establish the properties of the known members. We obtained an average extinction for the region of A₀ = 1.40 mag with a standard deviation of 0.52 mag (excluding two spectroscopically confirmed members with non-physical extinctions). This is in good agreement with previous results (A_V ∼1.5±0.5; A_V = 1.56±0.55 for Tr 37; A_V ∼ 1.74 mag, A_V ∼ 1.94 ±0.69 mag; Contreras et al. 2002; Sicilia-Aguilar et al. 2005, 2013; Nakano et al. 2012, respectively). We note that the distribution of extinction values is not Gaussian, extending from 0.5 mag towards larger extinctions up to < 3.1 mag, as expected due to the presence of dark clouds and to the typical lower limit galactic extinction 1 mag kpc⁻¹ distance. With this in mind, we consider as members those with extinctions in the range 0.9 mag < A₀ < 3.0 mag, which comprises our second membership criterion.

Targets in this extinction range with one additional membership criterion are considered sure members. Without an additional membership criterion, sources with a consistent extinction are considered probable members. Those outside the 0.9–3 mag range without other membership criteria are considered probable non-members. We note that an Anderson-Darling test (Stephens 1974) reveals that the extinction distribution of members and probable non-members do not exclude each other. This is likely caused by cluster members and background stars having often similar extinctions in the 1–3 mag range, which makes it hard to rule out objects following the extinction criterion alone.

We used the G versus G − G_RP diagram as a further consistency test (Fig. 3), correcting the objects by their individual extinctions. As explained before, the G − G_RP color offers better results, especially for faint red sources. We obtained three sources with negative (nonphysical) color excesses, which are labeled as uncertain in Table B.1. Two of them are confirmed members (with clear Li I and/or Ca II detection), and another is a probable member. These are faint M-type sources and with Gaia magnitudes G >  17 mag (before correcting for extinction), beyond the completeness limit of the survey (Gaia Collaboration 2021). Anomalous extinction can be caused by circumstellar material, scattering in the disk, spectral veiling, and/or uncertain spectral types. This is also a further sign that extreme care has to be taken with faint objects when determining ages or masses using the Gaia colors. For these three sources with unreliable extinction, we used the cluster average extinction (A₀ = 1.4 mag) in the color-magnitude diagram (Fig. 3).

3.3. Presence of protoplanetary disks

We used, as a third membership criterion, the presence of protoplanetary disks and/or accretion or X-ray emission (indicator of coronal activity, see Feigelson & Gaffney 2003; Mercer et al. 2009; Getman et al. 2012). YSOs with a full circumstellar disk have excess emission at λ >  2 μm. However, IR excesses alone need to be treated carefully as it could also correspond to post-main-sequence stars, background quasars, or surrounding nebulosity. We used mid-infrared photometry from the WISE catalog, considering only W1 [3.4 μm] and W2 [4.6 μm] bands because W3 [11.6 μm] and W4 [22.1 μm] are affected by the extended variable background emission for clusters at large distances (Ribas et al. 2014). The WISE data were used in combination with the JHKs 2MASS data. We also checked whether any disk signatures have been previously reported in the literature.

We consider as an IR excess > 3σ any deviations from the zero-color of stars without disks in both W1 − W2 and H − Ks. Only 67 spectroscopic sources have complete IR data with low uncertainties (< 0.05 mag). The W1 − W2 vs. H − K diagram (Fig. 4) reveals 27 targets consistent with having a disk (see Table B.1). Five have not been previously classified as disked stars, for another four have confirmed the presence of a disk. The rest had been previously classified as disked stars by other methods (Morales-Calderón et al. 2009; Barentsen et al. 2011; Getman et al. 2012; Nakano et al. 2012; Sicilia-Aguilar et al. 2013) and, in addition, our work now adds further information such as spectral types and extinction.

Fig. 4. W1 − W2 vs. H − K diagram for the 67 spectroscopic candidates with errors < 0.05 mag. The red box shows the 27 candidates with excess emission in the IR consistent with a disk.

Since most disks are expected to accrete, we contrasted this criterion with the veiling value calculated and accretion emission lines found (CaII and a strong Hα). Among the 27 sources with disk, 52% (14) show veiling and accretion lines, and 33% (9) show only accretion lines (strong Hα emission, EW(Hα) > 10 Å). Two sources clearly display the Li I absorption line, and another two sources have weak Hα emission. We complete this list with 23 more sources listed as disk-bearing in the literature. From these 23 sources, 11 show veiling plus accretion lines (Hα/Ca II), and 9 show accretion-characteristic lines only. We note that the main difference between the disks we found and the literature is that we do not use any wavelength longer than W2, which leaves sources with anemic and/or transitional disks out of our list.

The veiling measurements allow us to estimate the accretion rate for 25 sources. Assuming that the veiling is roughly constant as a function of wavelength (Hartigan et al. 1989; Dodin & Lamzin 2013), we can derive a relation between veiling and accretion luminosity (as in Calvet et al. 2000) using the data from Gullbring et al. (1998). As with other accretion estimators, the dispersion is large (∼0.5 dex). We then estimate the accretion rates (Ṁ_acc) assuming that the magnetospheric infall radius is 5R_* (Gullbring et al. 1998) and with stellar masses and radii derived from the pre-main sequence (PMS) evolutionary tracks (see Sect. 5.1). The results (Table 5) agree with typical accretion rates of Classical T Tauri stars (CTTS; 10⁻⁹–10⁻⁸ M_⊙ yr⁻¹, Fang et al. 2009; Sicilia-Aguilar et al. 2010) and with the presence of disks. Among the 25 sources with veiling measurements, 20 are confirmed members with accretion-related lines and disks and 5 have a low S/N spectrum and/or uncertain spectral types, so their veiling measurements are uncertain, but the overall veiling-disk trend is coherent.

3.4. Membership constraints from Gaia

We also used Gaia EDR3 data, with the maximum likelihood and Mahalanobis distance analysis (Mahalanobis 1936, 1969) as an independent membership test. The method is explained in detail in Sect. 4, so here we simply refer to that information to include an additional criterion for the spectroscopic candidates. Among the 120 candidates with Gaia EDR3 counterparts, we selected those with good quality astrometric data (see Sect. 2.4) and imposed the same conditions as for our probability analysis (see Sect. 4.2), and the same restrictions in the uncertainties used in Sect. 4.3, which leaves us with 83 sources.

A total of 38 of them satisfy the Mahalanobis criterion, having a probability of 95% of belonging to the region. Of the remaining 45 objects with good quality astrometric data, 28 are rejected because of their large errors, 3 are out of the age range to be considered as new members, and 14 do not belong to any subcluster of IC1396 (see Sect. 4.2). We note that the Mahalanobis criterion will lose bona-fide members due to the very stringent constraints imposed to limit contamination (see Sect. 4.3). Therefore, this criterion is used only to look for extra information for uncertain cases and not to cast doubts on the membership of objects confirmed by spectral lines, disks, or X-ray emission. The criterion can obviously not be applied to the sources with Gaia data below the quality threshold.

3.5. Combination of all membership criteria

Combining all four membership criteria for the spectroscopic candidates, we obtained a list with the final membership listed in Table B.1. As mentioned before, we give a higher priority to sources with spectroscopic lines indicative of youth. The final membership for the rest of the sources is defined depending on whether two or more criteria confirm the membership. We classified as probable members those sources that satisfy two or more criteria at the probable level.

From the initial list of 121 candidates, 66 are classified as members, 42 as probable members, and 13 as probable non-members. Two of the confirmed members and three of the probable members had not been previously classified as IC1396 members. The other 116 spectroscopic sources correspond to objects previously listed in the literature. In addition, our spectroscopic analysis supplies information on the previously unknown spectral types for 111 targets, accretion status for 25 targets, and interstellar extinction for 119 sources.

Among the 13 probable non-members, four have spectral types in the A-G type range. They also have uncertain Gaia photometry (G ≥ 18 mag). Nine sources are M-type stars without clear indications of youth (such as weak Hα emission or Li I line absorption too close to the noise level), which could be dMe stars. In addition, four of these have extinction outside the membership range we consider. From these 13 sources rejected as probable non-members, 6 were previously classified as members, 3 as uncertain, and the remaining one as a non-member.

4. Characterizing IC1396 with Gaia

In this section, we use Gaia data to study the distance and kinematics of the known members of IC1396, using this information to find new members. We also estimate the age, mass, and the presence of disks among the entire population. The analysis is done using the EDR3 data, although we compare and comment on the differences with DR2.

4.1. Selecting the best-known members in IC1396

The first step is characterizing the parallax and kinematics of the well-known members. A simple proper motion plot revealed deviations among the known cluster members suggestive of several sub-structures (see also Sicilia-Aguilar et al. 2019). We also expected differences in proper motion based on the radial velocity variations between different BRCs and the main Tr37 cluster (Patel et al. 1995).

We then studied the presence of distinct astrometric groups following Franciosini et al. (2018) and Roccatagliata et al. (2020). Significant clustering is identified by maximizing the likelihood function for the region known members, assuming that they are distributed in parallax and proper motion space following a number of multivariate Gaussians (Lindegren et al. 2000) that represent the position in parallax-proper motion space of each group. Iterating the procedure for various numbers of clusters allows us to determine which structures are significant and their corresponding parallaxes and proper motions. Further details are given in Appendix C.

The method requires two inputs: a collection of well-known members (rejecting all those that may be uncertain) and an initial subset of potential subclusters with their astrometric parameters (parallaxes and proper motions), their standard deviations, and the number of stars in each subcluster. We started with the 1538 known members (Table 6). We cross-matched them (using a 0.9″ radius and removing sources with multiple matches) to the Gaia data and imposed the same restrictions as in Sect. 2.4. We consider as reliable members those with parallaxes between 0.6–1.6 mas, obtaining 578 from EDR3 (see Table D.1), an increase with respect to the 536 objects with good DR2 data. These represent a very stringent selection of reliable members, which is needed to avoid introducing uncertainties via non-members or sources with unreliable astrometry. They were used to define the properties of the subclusters, but we must keep in mind that they represent only 53% of the known members with Gaia EDR3 counterparts. Thus, we expect that our search for new members based on them may miss a similar number of cluster members with poorer Gaia data. From the 578 known members, the mean parallax is 1.081±0.006 mas, and the average distance to the region falls to 925±73 pc. This value is consistent with the Gaia DR2 result of 943±82 pc, although the difference could reflect the known biases in the DR2 dataset (Stassun & Torres 2018).

The main difference between DR2 and EDR3 is that the uncertainties in the EDR3 data are significantly smaller (Gaia Collaboration 2021), especially on the proper motion. Proper motion uncertainties are typically reduced by half compared to DR2. We also gain 6% (85) members, and 5% of the DR2 sources are rejected because they fail to meet the goodness criteria (e.g., RUWE > 1.4, parallax values outside the range of stated values, or fractional error higher to 0.25). This number reflects the potential contamination of the members listed in the literature, although the fraction is small.

The initial subcluster parameters were obtained by examining the 3D astrometric space positions for the stringent collection of known members. The visually-identified structures allowed us to create an initial conditions file with 10 potential subclusters, the significance of which is then refined by maximizing the likelihood function. The likelihood function is computed by shifting the positions of these initial subclusters and revising them according to the data. Thus, if the initial subcluster parameters are reasonable and the sampled parameter space is large enough, this method will correct any imprecision in the initial subcluster positions. Subclusters that are not significant will also be automatically removed, as explained below.

4.2. The maximum likelihood algorithm for clustering

The significance of the subcluster structure is determined via maximizing the likelihood function using a customized Python routine. The algorithm uses the initial subclusters file, where each jth subcluster is defined, in a the multidimensional space, by seven astrometric parameters: the mean parallax (ϖ_j), the mean proper motion in right ascension and declination (μ_α, j, μ_δ, j), their intrinsic dispersion for the stars in the subclusters (σ_ϖ, j, σ_μα, j, σ_μδ, j), and the fraction of stars that belong to each subcluster (fs_j). We follow the same formulation used by Roccatagliata et al. (2018, 2020) where the probability of each ith star belongs to the jth subcluster is given by:

$$\begin{array}{c}\hfill P_{i,j}=f{s}_j\frac{L_{i,j}}{L_i},\end{array}$$(1)

where L_i, j is the likelihood function for the ith star belonging to the jth subcluster (see below), L_i is the total likelihood, and fs_j is the fraction stars of belonging to the jth subcluster. The individual likelihood L_i, j is calculated as:

$$\begin{array}{c}\hfill L_{i,j}={(2\pi )}^{2/3}{|C_{i,j}|}^{1/2}exp[-\frac{1}{2}{(a_i-a_j)}^T{C}_{i,j}^{-1}(a_i-a_j)],\end{array}$$(2)

and the total likelihood is obtained by summing all the individual ones multiplied by their corresponding fraction of stars for the total of n subclusters, namely,

$$\begin{array}{c}\hfill L_i={\displaystyle \sum _{j=1}^n}f{s}_j{L}_{i,j}.\end{array}$$(3)

In Eq. (2), C_i, j is the covariance matrix, |C_i, j| is its determinant, a_i corresponds to the vector for the ith star in the multiparameter parallax/proper motion space, a_i = [ϖ_iμ_α, iμ_δ, i], and a_j is the equivalent vector for the subcluster parameters, a_j = [ϖ_jμ_α, jμ_δ, j]. Transposition is denoted by ^T. Further details are given in Appendix C. Once the multidimensional location of the clusters is known, it can be used to identify new members using the Mahalanobis distance technique, which we describe in Sect. 4.3.

We calculate the likelihood for the well-known cluster members by varying the initial subcluster positions over a grid of initial parallax and proper motions distributed over ±2σ around the initial subcluster values (Lindegren et al. 2000). We tested various grid samplings, ranging from 14 to 34 points. In each loop, we first use the values of the initial parameter file as a guess solution, modifying them around the grid mentioned above to estimate the likelihood. The guess value is then refined once to allow for differences in the number of stars in each cluster (Franciosini et al. 2018) and to revise the central position of each subcluster parameter for objects that belong to the cluster. All the combinations of parallax and proper motion over the ±2σ space are tested, and each position of the subcluster is revised in turn. Therefore, as long as the initial parameter file is sensible and the number of grid points is reasonable, the final values are not significantly dependent on the initial cluster positions nor on the number of points in the grid.

With this procedure, subclusters that are not significant disappear since the probability of a star belonging to that subcluster vanishes as the probability becomes negligible when all stars are too far from the center parameters. Non-significant clusters also tend to merge with others with more members. For instance, we started with ten initial subclusters that merged into only seven statistically significant subclusters. We note that the sizes of these seven subclusters do not exceed the limit imposed by Kounkel et al. (2018) to avoid linking physically unrelated structures or contaminating field stars.

The combination of position and number of subclusters with the maximum likelihood is kept as the best solution for the structure of the region, with the standard deviation of the parameters of all stars assigned to each subcluster representing the distribution of subcluster members. For the final result, we consider those stars having a probability higher than 80% of belonging to one subcluster as subcluster members. We performed several additional tests by changing the subclusters initial conditions and the number of subclusters, always converging to the same results.

Table 7 shows the resulting best-fit parameters for the seven surviving subclusters, which we labeled with letters A through G. A total of 552 stars (among the 578 known members with good astrometric data, see Table 6) have a probability of > 80% of belonging to one of these subclusters. This does not mean that objects with a lower probability are less likely to be members of IC1396, but they may not be so clearly associated with one of the particular subclusters. For instance, we identify 26 sources with intermediate probabilities, 4 of them with probabilities > 40% of belonging to two subclusters (which are considered probable members of both) and the rest with probabilities less than 80% (considered probable members of one of the particular subclusters). Figure 5 shows the distribution of parallax and proper motion of the stars associated with each subcluster and their positions in the color-magnitude diagram (see also Table 7).

Fig. 5. Parallax histogram (top, including the full histogram on the left and a zoom of the subclusters on the right), proper motion and color-magnitude diagram (bottom) of the known members of the subclusters found by the maximum likelihood analysis with the Gaia EDR3 data. The symbols and colors mark the various subclusters using circles, triangles, and capital letters (see legend). The number of subcluster members is given in parentheses. Full circles represent stars with a membership probability > 80%. Triangles represent objects with probabilities of < 80%. The theoretical PARSEC isochrones on the color-magnitude diagram have been corrected for a distance of 925 pc and a lower limit extinction AV = 1 mag to reject foreground stars. The extinction vector for 1 mag is indicated with a red arrow.

The distances of subclusters A, B, E, and F are all consistent (although subcluster F has only seven members). This is also true for the distances of groups C, D, and G. Population G (64 sources) is poorly defined in parallax and proper motion (e.g., it is extended from μ_α = −20 mas yr⁻¹ to +15 mas yr⁻¹). It also has a consistently larger age than the rest, with 58% of its members below the 10 Myr isochrone (Fig. 5). Therefore, the extended population G is not considered further to search for new members, as it could include too much contamination. However, many of its members are bona-fide cluster members, with clear indicators of youth, such as disks or accretion. This shows that IC1396 contains a significant number of members that are very extended in the proper motion plane.

Subcluster D has nine members, of which four are intermediate-mass (B, A, F spectral types) stars considered as probable members (Contreras et al. 2002), and one additional source identified by Nakano et al. (2012) with low EW(Hα) < 10 Å, which makes them less reliable. Since its proper motions overlap with the densest part of the galactic plane, it is likely that many, if not most, of the objects associated with subcluster D, are contaminants previously thought to be members. Subcluster D also has a different parallax (Table 7), being more distant than the rest of the subclusters, followed by the extended population G and subcluster C. The fact that subclusters C, D, and G distances are consistently higher suggests that these groups could also contain a significant number of contaminants. Therefore, we did not use subclusters C and D to search for new members either.

Following on the previous findings, we checked if the distance distribution of the sources belonging to each subcluster is significantly different using the Anderson-Darling test (Stephens 1974). This non-parametric test examines the total difference in the cumulative distribution. We used the scipy.stats module of the Anderson-Darling test for k-samples (Scholz & Stephens 1987; Virtanen et al. 2020). The test returns the test value, the critical values, and the significance levels (25%, 10%, 5%, 2.5%, 1%), which can reject the null hypothesis. We consider the samples are statistically different if the significance level is less than 1%. Table 8 shows the test results. Although subclusters C, D, and F have few members, limiting the significance of the test, we find a significant difference in distance between subclusters A and D; A and C; B and C; B and D; C and E, and D and E. The distances of subclusters C and D are not significantly different, but both C and D are more distant than the rest, while subclusters A, B, E, and F all share a similar distance. For comparison, the proper motions of the subclusters are different from expected since, otherwise, they would have been merged by the algorithm.

We also conclude that the EDR3 data has made the identifications of the DR2 subclusters more reliable. The maximum likelihood analysis with Gaia EDR3 data confirms and recovers previously identified subclusters (and the extended population) with the Gaia DR2 data. The reliable subclusters A, B, E, and F, are essentially identical, but the results with Gaia EDR3 data make them more compact and thus better defined in the parallax-proper motion space. In contrast, extended populations D and G do not improve. The objects associated with subcluster C have a lower spread in proper motion in EDR3 than in DR2. However, the number of members of subcluster C remains so low that a further independent study (including more confirmed members around this region) is required before its status as a subcluster can be confirmed.

4.3. New candidates identified via the Mahalanobis distance

The Mahalanobis distance technique can be used to identify objects that belong to the same subclusters as the known members (Mahalanobis 1936, 1969). The Mahalanobis distance is equivalent to the Gaussian probability, generalized to standard Euclidean distance in a transformed, multidimensional space (see Appendix C). Sources with a high probability of belonging to the same distribution than the known members are considered new members. For a collection of i stars and j subclusters, we estimate the probability of each ith source to be a member of each jth subcluster, calculating the multidimensional distance between the star position and the subclusters centers in the parallax-proper motion space (Mahalanobis distance), using the Python routine scipy.spatial.distance.mahalanobis. As it happens in a 1D Gaussian space, a Mahalanobis distance lower than 2 corresponds to a ∼95% confidence interval regarding the membership of the object to any subcluster with the multivariate position considered.

We used the Mahalanobis distance to find new members of IC1396. To this end, we drew upon the reliable subclusters (A, B, E, F) from the maximum likelihood analysis to estimate the membership probability of other stars with good quality Gaia data. Our first tests showed that the parallax positions derived from the maximum likelihood method were poorly defined due to the typical large uncertainties for objects at nearly 1 kpc, which resulted in excessive contamination. We thus readjusted the parallax distribution with a Gaussian function to the best-defined subcluster (A). The fitting produced a mean parallax of ϖ = 1.100 mas, with standard deviation σ = 0.054. Since there are no significant differences in the distances of the rest of the well-defined subclusters, we used this value for all.

We explored the membership of the sources in a 2° radius field around HD206267. Imposing the same restrictions, we obtained about 68k sources with good quality astrometric data and in the parallax range (0.6–1.6 mas) of known members. We calculated the Mahalanobis distance and selected those within the 95% confidence interval as members. We further reduced contamination by restricting the errors of the magnitudes (G, and G_RP) to less than 0.05 mag, and the errors on proper motion to less than 0.1 mas yr⁻¹, to reduce contamination from those sources with large errors⁶. In addition, we imposed an age cut between 0.1 Myr and 10 Myr, assuming a fixed extinction of A_V = 1.0 mag, avoiding excessive contamination by those objects at ≥1 kpc with relatively large uncertainties. This is particularly important for IC1396 since the proper motions of the members are not so different from those of the galactic background. According to our maximum likelihood results and considering only the four reliable subclusters (Table 9), ∼91% of known members are in this age range, so there is no substantial evidence of an older population. This also implies an approximate loss of ∼9% of new members, to add to these what we excluded from the extended population (G), representing 11% of total known members. Although our final new members list is incomplete, the limits are required to avoid too much contamination by sources that are unlikely members.

The procedure described above resulted in a total of 334 new members (see Table 6) belonging to the subclusters (see Table E.1). Among these, five sources are shown to belong to both subclusters B and E (with probabilities > 30–60%). Figure 6 shows the new members distribution in the color-magnitude diagram. Table 9 shows the final result of the total number of members per subcluster from the Mahalanobis distance analysis and from the maximum likelihood analysis.

Fig. 6. Color-magnitude diagram for the new members obtained from Mahalanobis distance analysis with the Gaia EDR3 data. The 334 new members between 0.1 and 10 Myr are in full color, candidates with ages > 10 Myr are rejected and shown in semi-transparent color. The theoretical PMS tracks have been corrected for a distance of 925 pc and a lower limit of extinction of AV = 1 mag to reject foreground stars. Five members are plotted as triangles and belong to both subclusters B and E. We note there is an empty space between the 10 Myr and 15 Myr isochrones, which reinforces the criterion of considering as members only those up to 10 Myr. A vector extinction is indicated by a red arrow.

We analyzed the spatial distribution of all known and new members of IC1396 (Fig. 7). Most of the known members are clustered around HD206267, while most new members are near the BRCs and on the filaments around the region. Although the coordinates have not been used to identify the subclusters, most of the members of the subclusters A, B, and E are not uniformly distributed in space. This spatial structure and clustering, together with the differences in proper motion between subclusters, are additional confirmations that they have different physical origins. Subcluster A is the best defined, located in the center of the region, and it has few new members because the IC1396A globule and the Tr37 cluster have been extensively studied. The spatial distribution of subclusters B and E members are clearly different from subcluster A, suggesting various star-forming episodes.

Fig. 7. Spatial distribution of the known (full-color circles) and new (pentagon) members. The colors represent the four subclusters (A, B, E, F, see legend). Orange circles mark the rest of the known members not belonging to any subcluster. The age of members (derived in Sect. 5.1) is represented by the symbol sizes, with the youngest members having large symbols ranging from 0.1 to 20 Myr. Black numbers indicate the main BRCs, 1=IC1396 A, 2=IC1396 B, 3=IC1396 D, 4=IC1396 N, 5=IC1396 G. The background image is a ∼4° × ∼4° mosaic built from WISE 22.19 μm.

Subcluster B mostly runs along the northern part of the region, besides IC1396N, with some objects extending towards the northwest behind BRC IC1396D. The members of subcluster E are east of HD206267, spreading to the northwest past the ionized edge of the bubble. Subcluster F has too few members to discuss its spatial location meaningfully despite being well-defined in the parallax-proper motion space. In contrast, the members of subcluster D and the extended population G are not spatially clustered, which is expected from groups that are not real, are heavily contaminated, or are not well defined in astrometric space. The lack of clear astrometric and spatial signatures further confirms that, even if a fraction of confirmed IC1396 members belong to these groups, we cannot use them to identify further members (Sect. 4.2). However, subcluster C, also rejected, is spatially clustered in the south of the region, but as we mention above, it has too few members to be confirmed.

We used the 2MASS J − H vs. H − K diagrams to explore the properties of the new members. Figure 8 shows that most of the new members are consistent with diskless intermediate-mass and massive T Tauri stars, in contrast with the known members, which are mostly low mass T Tauri stars (see Sect. 5.1 for a detailed mass analysis). This is likely a consequence of the requirements imposed on the astrometry quality. It also reflects that the most efficient methods used to identify members (X-ray emission, Hα excess from photometry or spectroscopy, Li I absorption, and the presence of near- and mid-IR excesses) are biased towards revealing solar-and low-mass stars. For instance, we found 43 new members (23 in subcluster A, 6 in subcluster B, and 14 in subcluster E) at (J − H) < 0.2 mag. From these 43 sources, 33 are known stars with spectral types A0-B1 (Alknis 1958; Marschall & van Altena 1987; Kun 1986; Patel et al. 1998; Getman et al. 2007), which are now found to be cluster members.

Fig. 8. J − H vs. H − K diagrams for known members associated with the various subclusters obtained from the Maximum likelihood analysis (colored circles) and new members obtained from the Mahalanobis distance analysis (magenta pentagons). The number of sources with complete JHK photometry is indicated in the legend. Subclusters B and E include the new members with an intermediate probability of belonging to both subclusters. The MS and the giant branch are indicated by thin solid lines (Bessell & Brett 1988). The CTTS locus (Meyer et al. 1997) is indicated by a thick solid line. The theoretical tracks are corrected for extinction using the threshold AV = 1 mag. A reddening vector for AV = 3 mag is shown as a red arrow.

Figure 9 shows the velocity dispersion in the plane of the sky for all members relative to the mean proper motion of subcluster A. The relative motions also indicate that the subclusters members are dispersing. Further analysis of the velocity structure of the region is discussed in Sect. 4.4.

Fig. 9. Map of the velocity dispersion in proper motion of known and new members. The arrows mark the direction and size of the proper motion vectors relative to the average proper motions of subcluster A. They are color-coded by subcluster (see legend). Known members not associated with any subcluster are shown in orange. The long orange arrows are mostly from the extended population G. New members are marked by open circles. The magenta arrow in the bottom left corner indicates the physical velocity scale (at 925 pc).

Comparing with Gaia DR2, the Gaia EDR3 errors are significantly lower, especially for the proper motion (see Table 7). This produces a tangible improvement in membership selection. Gaia EDR3 allows us to reject 62 objects which had been previously identified as new members with the Gaia DR2 data. Eight of them now have RUWE > 1.4, 30 have photometric errors > 0.05 mag, and the rest were border-line cases for which the probability of belonging to the subclusters has now dropped below the 95% confidence levels due to the reduction in the uncertainties. The fact that the Gaia EDR3 data have smaller uncertainties and, still, the number of members increases is an additional confirmation that the subgroups are real. Non-existing groups affected by significant uncertainties would tend to disappear when the uncertainties are reduced.

Recent studies of the west of Tr37 (Silverberg et al. 2021) found over 400 members via X-rays, some of which were previously known members. We did a very brief test finding that only 45 of them have Gaia EDR3 counterparts, and 10 were already in our collection of known members. Finally, from these, only 22 sources have good quality astrometric Gaia data. We added these 22 sources to our 578 reliable known members (see Table 6) to check the robustness of our method, we calculated the maximum likelihood function again and obtained the same distribution of subclusters in the entire region. The only difference was a slightly greater overlap in proper motion between cluster B and cluster E with respect to previous results since some more of the new members have an intermediate probability of belonging to both subclusters B and E. This happened to a lesser degree before, and it may indicate that these two subclusters are, in fact, part of the same structure that spans a continuous range of proper motions. To test the physical relation between both subclusters, we would need more objects in each group, detected by means independent of their proper motions, which is left for future research. Looking for new members, we did not find any changes with respect to our previous analysis since the subcluster locations do not change significantly, proving that our study is robust against adding small numbers of new sources.

4.4. The velocity structure of IC1396

In this section, we present a complementary analysis of the kinematic structure of IC1396 using the Gaia EDR3 data for all members. This analysis also helps us to give a more quantitative explanation of the apparent expansion trend observed in Fig. 9. We examine the velocity structure of IC1396 using the 578 known members (Table D.1) and 334 new members (Table E.1). Since only very few objects have reliable radial velocities (and they belong to the central part of subcluster A; Sicilia-Aguilar et al. 2006b), it is not possible to analyze the kinematics in 3D, so we concentrate on the velocity structure on the plane of the sky using the proper motions.

The velocity structure is explored using the Velocity Structure Analysis Tool (VSAT, Arnold & Goodwin 2019), which is a Python routine that estimates the projected velocity differences (Δv_ij) between any pair of stars (i,j) in the region, together with its corresponding projected source-to-source distance (Δr_ij). We examine the relation between Δv_ij and Δr_ij, which gives us information about velocity gradients to reveal bulk motions within the cluster. The object pairs are binned in the projected source-to-source distance (Δr), calculating the average projected velocity difference (Δv) per bin. Provided that the number of pairs in each bin is large enough, the width of the bin should not affect the results. Observational errors (σ_Δvij) are propagated when calculating each Δv_ij, and each Δv_ij is weighted accordingly ($\omega =1/{\sigma }_{\mathrm{\Delta }{\mathit{v}}_{ij}^2}$). The errors in the spatial coordinates are negligible compared with the uncertainties in the projected velocities. The code returns Δv(Δr) showing the average Δv values per bin, against their corresponding Δr. The method tracks two aspects of the projected velocity structure of the region, the magnitude (Δv_M) and the directional (Δv_D) projected velocity difference. Here we include a brief explanation, since the details can be found in Arnold & Goodwin (2019). If each pair of stars (i,j) has vector velocities (v_i, v_j), then the magnitude of the projected velocity difference (Δv_M) is |v_i − v_j|, which is always positive. This definition is a general measurement of how similar or different the velocity vectors are. In a two-dimensional space (x, y), Δv_M is calculated as:

$$\begin{array}{c}\hfill \mathrm{\Delta }{v}_{\mathit{ijM}}=\sqrt{{(v_{\mathit{xi}}-v_{\mathit{xj}})}^2+{(v_{\mathit{yi}}-v_{\mathit{yj}})}^2}.\end{array}$$(4)

On the other hand, the directional projected velocity difference, Δv_D is the rate at which the projected source-to-source distance (Δr) between a pair of stars changes, which helps to understand if the region is expanding or collapsing. If the projected distance between stars (Δr) increases and the bulk of stars are moving away from each other, Δv_D is positive; otherwise, Δv_D will be negative. Therefore, Δv_D is calculated as

$$\begin{array}{c}\hfill \mathrm{\Delta }{v}_{\mathit{ijD}}=\frac{(x_i-x_j)(v_{\mathit{xi}}-v_{\mathit{xj}})+(y_i-y_j)(v_{\mathit{yi}}-v_{\mathit{yj}})}{\mathrm{\Delta }{r}_{\mathit{ij}}}.\end{array}$$(5)

VSAT does not require the center of the region or its radius to be defined, and the frame of reference is irrelevant since all measurements are relative. The results are not strongly affected by the observational uncertainties and biases of individual objects, which are washed out by binning the data. For better visualization, we convert RA and DEC coordinates to parsecs and the proper motions to km s⁻¹ (see Fig. 10). The velocity structure of the region is quite clean (bottom panels) if we do not consider the members of the extended population G, which has very large and significantly different proper motions compared to the clustered populations. Population G could include ejected members, although it is also likely affected by contamination at a higher level than the clustered members. If we consider the extended population (top panels), although the tendency is the same as before because the number of objects in the extended population is small, the noise increases due to the larger velocity differences. In any case, there are signs of expansion up to distances of Δr <  12 pc.

Fig. 10. Velocity structure in the IC1396 region, shown in terms of the mean difference between the velocity magnitude ΔvM (left panels) and the directional velocity ΔvD (right panels) versus the difference distance in parsec (Δr; see text). The two top panels show the region velocity structure considering all known and new members, including the extended population (G), which has considerably larger velocities and velocity dispersion than the rest of clustered members. The two bottom panels include only the clustered part of the population, excluding subcluster G, and thus offer a better view of the evolution and dispersal of the clustered members.

In analyzing the velocity structure (Fig. 10, bottom panels) without the extended population, the magnitude of projected velocity difference Δv_M(Δr) shows an increase in the velocity differences between nearby stars, rising from 2 km s⁻¹ to 6 km s⁻¹ for stars with Δr up to ∼12–14 pc (Fig. 10, left-side panels). This shows that nearby stars have similar velocity vectors, so Δv_M is small. Stars that are more distant from each other are less connected, as evidenced by their very different velocity vectors and large Δv_M. The region shows a clear velocity structure until Δr ∼ 12–14 pc. At relative projected distances Δr >  12 pc and up to Δr ∼ 30–32 pc, Δv_M remains almost constant at ∼5.5 km s⁻¹, not revealing any significant structure at these scales. When comparing stars separated by even larger distances (Δr >  30 − 32 pc), Δv_M appears to increase again, but this result is uncertain since large distance bins have fewer stars and noise increases.

For the directional projected velocity difference Δv_D(Δr), (Fig. 10, right panels), we find a positive correlation between Δv_D and Δr, which means that the stars are moving radially away from the center of their velocities distribution. We observe expansion between stars up to separations Δr ∼12–14 pc, and the expansion velocity increases from ∼0.6 km s⁻¹ to ∼4.8 km s⁻¹ respectively. Those projected velocities exceed what we could expect from gravitational attraction due to the mass in the region, a sign that the region is not gravitationally bound. This is in agreement with the values obtained by CO maps (Patel et al. 1995, ring radius = 12 pc, and V_exp = 5±1 km s⁻¹), also suggested by Getman et al. (2012) and Sicilia-Aguilar et al. (2019). Since we are using all the stars, this gives us a better perspective of the velocity of the entire region. Stars separated more than ∼12 pc belong to physically unrelated populations. This is in agreement with the stellar crossing times for sources with typical dispersion velocities 2–3 km s⁻¹ (see also Sicilia-Aguilar et al. 2006b), which would travel between 8 and 12 pc in a cluster lifetime of 4 Myr. At separations between ∼12 pc and ∼30–32 pc, Δv_D does not show any significant trend that could hint at either expansion or contraction, probably due to noise. Beyond separations of ∼30–32 pc, the noise level increases even further (see Fig. 10, right panels).

In summary, considering Δv_M, we observed a velocity structure on distance scales up to 12 pc. Considering Δv_D, the region is also expanding until up Δr ≈ 12 pc. This 12 pc scale is similar to the distance of the compressed molecular ring around the periphery of the IC1396 region (Patel et al. 1995). There is no significant velocity structure on scales beyond that.

5. Properties of the IC1396 members

5.1. Ages and masses for new and known members

We derived the age and mass for new and known members (see Fig. 11). The magnitudes were corrected using the average extinction, A_V = 1.40±0.52 mag, and its standard deviation provides the uncertainties. As in previous sections, we used the PARSEC isochrones, corrected by the 925 pc mean distance. We derived individual ages and masses by interpolating the isochrones from 0.1 to 20 Myr, using a spacing of 0.1 Myr, to fit the color-magnitude data (G vs. G − G_RP). This offers a better solution than G_BP − G_RP since the G_BP magnitudes are often underestimated, leading to an offset in age, although the relative age values are maintained in both cases.

Fig. 11. Age (left) and mass (right) distribution of known members (top panels) and new members (bottom panels) from different subclusters. We have included those known members that do not belong to any subcluster and those that belong to one more subcluster (labeled as others “O” in the legend). Some new members belong to two subclusters, marked as “B/E”. Objects with ages 20 Myr and over are all included in the last bin.

For the age calculation, we excluded 52 objects with G ≲ 13 mag since they are already in the main sequence and the separation between isochrones in this age range is significantly less than the uncertainties in magnitude and extinction. Among the known members, we considered only those with magnitude errors < 0.05 mag (537). The mean age is ∼4.0 Myr, with a standard deviation of ∼5.0 Myr. The mean age for the 296 new members is also ∼4.0 Myr (standard deviation ∼4.0 Myr), although we note that we imposed an age cut for selecting new members⁷.

We also obtained the mean ages for each subcluster: A (∼2 Myr), B (∼4 Myr), E (∼4 Myr), and the extended population G (∼13 Myr). The older age for G, together with its very large age spread, suggests once again contamination. Subcluster F has only six known members, so its mean age is uncertain. An age histogram is shown in Fig. 11. Since we had imposed a limit on the age for the new member selection, there is a bias against older ages. Considering only the known members, we found that subcluster A has 2% members with ages > 10 Myr, subcluster B has 12% members with ages > 10 Myr, and subcluster E has 7%. Therefore, it is to be expected that the new member determination is worse for the older regions.

The mean masses are 0.7 M_⊙ and 1.0 M_⊙ for the known- and the new members, respectively, so the region is dominated by low mass stars. There are only 4% intermediate-mass stars (> 5 M_⊙) among the known members, but 14% of the new members fall in this range. This confirms that we are preferentially identifying more massive sources, an expected bias considering that most surveys were more sensitive to lower-mass stars and the Gaia quality data requirements. We find six new members with M >  5 M_⊙ previously identified as MS stars with spectral types B1–B8. The most massive stars are HD 239689 (B2, Patel et al. 1998) for which we find ∼11 M_⊙, and LS III +57 12 (B1, Marschall & van Altena 1987), with ∼18 M_⊙. Other massive members confirmed by Gaia are BD+56 2622, HD 239731, BD+57 2363, and BD+58 2294. Among the previously-known members, the most massive stars are HD 239729 (B0, Patel et al. 1998; Mercer et al. 2009, ∼10 M_⊙) and LS III +57 14 (B1, Mercer et al. 2009, ∼16 M_⊙). Figure 11 also shows the final distribution of masses for known and new members.

We use the Anderson-Darling test to check the statistical differences of the age and mass distributions of known and new members. Here, we must note that although we can compare the known members of different subclusters between themselves since we have imposed a cut in the isochrone for the final selection of new members, comparisons of the age differences between known and new members can only be performed after making the same cut for both collections.

First, we tested the age and mass distribution of known members in subclusters A, B, E, and G. We do not consider subcluster F because it has too few members. Subcluster E has only 11 members, which limits the reliability of the test. The results are shown in Table 10, where the age distributions are statistically different between subcluster A and subclusters B, G, and probably E. There are no statistically significant differences between the mass distributions, except for the subclusters B and G, but this is because population G has many members with larger masses than subcluster B, which again evidences contamination. The age differences between subclusters A and B, and E (despite their low number of members), further confirms that the IC1396 region is very complex and contains various star formation episodes.

We also ran an Anderson-Darling test for the total population in subclusters A, B, and E. For this, the known members are cut the same way as new members, adopting the fixed extinction A_V = 1 mag and those between isochrones 0.1–10 Myr. Table 11 shows that the age distributions of subclusters A, B, and E are statistically different. There are no significant differences in the mass distribution other than those induced by poor statistics in the clusters with few members. Our results confirm age differences between subclusters A, B, and E, suggesting various episodes of star formation giving rise to the subclusters. This adds to the previous evidence of age differences and triggered and sequential star formation (Patel et al. 1995; Reach et al. 2004; Sicilia-Aguilar et al. 2006a, 2019; Getman et al. 2012).

Analyzing subclusters B and E together with the spatial distribution of their members, we find that they are not smoothly distributed and cluster in front of the BRCs. Subcluster B members tend to appear in front of IC1396N, and subcluster E members group in two regions, in front of IC1396G in the east and behind IC1396D towards the north-west (see Fig. 7). We tested the age distribution of the subclusters B and E members in spatial proximity to the BRCs versus the rest of the non-clustered B and E members for the entire population (known and new members). The Anderson-Darling test shows that for subcluster B, the clustered and non-clustered members are statistically different at a significance level < 1%. For subcluster E, the age distributions between the two clustered populations and the spread members are not significantly different. Therefore, there is some evidence that clustered objects are younger, on average, than the spatially extended population, despite both belonging to the same proper motion group. This is, to some extent, what is also observed between the IC1396A and the Tr37 population, which have similar proper motions but different ages (Sicilia-Aguilar et al. 2019). Nevertheless, the picture is not complete because although there are clearly younger and less evolved sources in the bright-rimmed clouds (Beltrán et al. 2009; Reach et al. 2004; Sicilia-Aguilar et al. 2006a, 2014) and we confirm this for subclusters A and B (populations related to the BRC IC1396A and IC1396N, see also Sicilia-Aguilar et al. 2019), Gaia cannot detect embedded objects, so we are missing part of the population.

The global ages observed in the region, comparing the center and the outskirts, do not point to star formation triggered by HD 206267 since subcluster A is, on average, younger. Age distributions where the central part of the cluster is younger than the outskirts have been observed (Kuhn et al. 2014, 2015; Getman et al. 2014, in NGC 2024, Orion nebula cluster, respectively). More complex, non-triggered structures may result from the coalescence or different star formation processes Smilgys & Bonnell (2017).

5.2. Disk fraction among known and new members

Since the known members had been identified with different techniques (X-rays emission, Hα excess from photometry or spectroscopy, and near- and mid-IR excess), it is necessary to characterize the presence of disks uniformly so that we can compare those properties with the newly identified members. We built a J − H vs. H − K diagram using the JHKs photometry from the 2MASS survey to characterize as many members as possible. After collecting 1538 known members, we obtained 591 JHKs counterparts (with a matching radius of 0.6″) and photometric errors < 5%. The color-color diagram is consistent with intermediate- and low-mass members (see Fig. 12). A significant part of the objects lies within the locus of the CTTS (Meyer et al. 1997) and likely harbor primordial disks, but longer wavelength data is required to verify the disk fraction independently of extinction effects.

Fig. 12. J − H vs. H − K diagram for all known members (KM) with magnitude uncertainties < 0.05 mag. The MS and the giant branch (Bessell & Brett 1988), as well as CTTS locus (Meyer et al. 1997), are indicated by a solid line. The theoretical tracks are corrected for extinction in JHKs (Cardelli et al. 1989), using the threshold AV = 1 mag. The extinction vector is indicated with a red arrow. KM in full color belong to the populations (A, B, C, D, E, F, and G) found from the maximum-likelihood analysis.

We further revised this classification using the WISE data combined with 2MASS, obtaining 518 known members with complete, good-quality data in the W1, W2, H, and Ks bands and errors up to 0.05 mag. The W1 − W2 vs. H − Ks diagram (Fig. 13) reveals 200 sources with excess emission > 3σ, consistent with protoplanetary disks, or a disk fraction of ∼39%. This fraction is slightly lower than what has been derived from other studies (e.g., 45% from Sicilia-Aguilar et al. 2006a), which likely arises from the fact that we are not able to detect disks with low excess or inner gaps. We tested this classification using the 148 objects in common with Sicilia-Aguilar et al. (2006a, 2013), which were classified using the multi-wavelength spectral energy distribution (SED), including optical and Spitzer data. With our method, we recover all the diskless objects. All those we detect as harbouring disks correspond to disky sources, and we only fail to detect as disks 10% of the weakest disks classified via the SEDs, consistent with transitional or settled/depleted disks. This means that our methods agree very well, so the disk populations, excluding transition and depleted disks, are comparable. We also note that we are studying a more extended region, where ages and disk fractions could differ.

Fig. 13. W1 − W2 vs. H − K diagram for all members with magnitude errors < 0.05 mag. Known members (KM) are shown in green and new members (NM) in magenta. The green and magenta boxes enclose the known and new members, respectively, with color excesses consistent with disk emission.

We also investigated the presence of disks around the 334 new members (see Table 6) identified with Gaia, of which 240 have W1 − W2 and H − K data. Figure 13 reveals 11 sources consistent with excess emission larger than 3σ. Therefore, the disk fraction for new members is 5%, significantly lower. This could be due to the fact that with Gaia, we are preferentially identifying more massive stars, which are thought to lose their disks faster (Contreras et al. 2002; Hernández et al. 2005). In addition, the known members were preferentially identified using methods that favor lower-mass and disked stars (e.g., photometry, spectroscopy, infrared excesses), while Gaia identifies sources independently from their disk status.

The final disk fraction for stars with all masses, putting together the objects, is 28%. We also revised the disk fraction for subcluster A, the most populous one, which is 30%. This result agrees with the disk fraction calculated by methods unbiased towards diskless sources, such as X-ray; for instance, the ∼39% fraction derived by Mercer et al. (2009), the 29% derived by Getman et al. (2012), and the 33% fraction from Silverberg et al. (2021).

5.3. Feedback evidence in the BRCs

We searched for evidence of interaction and feedback between the stars and the cloud in IC1396A, IC1396B, and IC1396N. For this purpose, we combined the [S II], H₂, and the Ks band images to unveil the places of interaction between the stars and the clouds at different depths identifying shocks, knots, and jets. We also compared [S II] (which results from combining the 6716 and 6730 Å images) with 6750 Å and H₂ with Brγ to distinguish them from continuum emission.

We detected substantial emission in H₂ and [S II] at the edge of all the BRCs (Figs. F.1–F.3). The combination of these two emission lines and the fact that the extinction in the globules is larger than A_V ∼ 1 mag, the limit for significant FUV irradiation (Burton 1992), rule out the possibility of fluorescent emission caused by far-UV field from HD206267. Some fluorescent emission may nevertheless be present at the edges of the globules. For IC1396N, we clearly detect several knots and chains of knots in the deep continuum-subtracted H₂ images, which match up with previous observations (Nisini et al. 2001; Caratti o Garatti et al. 2006; Reipurth et al. 2003; Beltrán et al. 2002, 2009; Rebull et al. 2013). A strong chain of knots (labeled B1, B2, and B3, Beltrán et al. 2002) is also clearly visible in the JHKs images, and the weakest H₂ features are also observed in Ks band due to the width of the filter. The H₂ knots lie in the same direction than the outflows found by Beltrán et al. (2002, 2012). For the N and S outflows (driven by the C and I sources, respectively; Codella et al. 2001), the trajectories followed by the knots suggest the collision of both outflows, as recently pointed out by López et al. (2022).

We also detected shocks, possibly the result of strong outflows coming from YSOs inside and breaking the wall of the cloud, impacting the material surrounding IC1396N. Bow shock 3 (HH 777, Reipurth et al. 2003) is more intense than the others, and possibly driven by the near infrared (NIR) source 331A (21:40:41.4 +58:16:38.4; Reipurth et al. 2003). The same source could be responsible for the other bow shock just below bow shock 3. Shock 2 could be driven by the 1.3 mm continuum source IC1396N-C (source C from Codella et al. 2001), and it would be associated with the chain of knots around the source. We detected 2 weaker shocks, numbered 1 and 4. All shocks are shown in Appendix F.

In IC1396A, we detected several known members inside and outside the globule, and few new members, mainly scattered around the globule, since most members in the globule were already known. We observe a small jet close the “eye” of IC1396A (Fig. F.1). In IC1396B, we observe mainly the shocked material on the small irregular edges of the cloud by the combination of H₂ and [S II] emission. We note that most of the known members are around the globule, and we do not find any new members, suggesting a lower rate of embedded star formation, as mentioned before. In IC1396B, most of the known members are found around the globule and we did not find any new embedded sources with our study, which suggests a lower rate of embedded star formation. However, our study has a bias against very low mass and highly embedded stars. For this globule, most of the excited and shocked material we observe in the figure is probably caused by external sources and appears along the irregular edges of the cloud (Fig. F.3).

6. Summary and conclusions

In this paper, we study the star formation history and the global structure of IC1396, including the Tr37 cluster, and the properties and kinematics of their young stars. For this purpose, we first identified, classified, and confirmed new members using optical spectroscopy. We then compiled the complete collection of known members, including those from the literature, and used their parallaxes and proper motions to characterize the region. This has allowed us to identify further stars consistent with cluster members. The revised collection of members permits the study of the disk fraction, ages, and masses in a more uniform manner.

The Gaia data reveal a complex star formation history. In particular, IC1396 is composed of several substructures that are clearly distinct in proper motion and, some of them, in age. This, together with the spatial distribution of members, suggests a multi-episodic formation process. We also find that the region is not gravitationally bound and is expanding as a whole. Expansion is also evident in subcluster A. Gaia EDR3 data confirms a distance of 925±73 pc, with the bulk of members being distributed around 113 pc (standard deviation), slightly closer than the value derived from DR2 (943±82 pc). A detailed summary of the results is provided below.

We obtained optical spectra for 121 sources, deriving their spectral types, interstellar extinction, and accretion rates (from veiling), and reviewed their membership. Among those, we confirm 66 as members (2 of which are new members) based on their spectral features, extinction, presence of disks, and/or kinematics. An additional 42 are classified as probable members. We also rejected 13 sources as probable non-members, of which six were previously classified as members. For 111 of them, there were no previous estimates of the spectral type. The veiling measurements allow us to estimate accretion rates for 25 sources, which are in the range of 10⁻⁹–10⁻⁸ M_⊙ yr⁻¹. These accreting sources consistently display strong emission lines and have IR excesses from disks. We derived individual extinction values and obtain an average extinction of A_V = 1.4 mag for the region.

Using Gaia EDR3, we confirm that IC1396 contains four independent subclusters (A, B, E, F) that are significantly different in proper motion but not in parallax. The various populations are located in clearly distinct areas around the bubble, with subcluster A filling most of the center, around HD206267. Subclusters B and E run along the ridges and BRCs associated with IC1396N, IC1396D, and IC1396G. Subcluster E is a new subcluster entirely identified in our analysis. In addition, we also observe a population very extended in proper motion. Although part of the extended population could correspond to contaminants, it also clearly contains bona fide members, showing that IC1396 has some high proper motion members.

We identified 334 new members using the multidimensional (parallax and proper motion) information provided by Gaia in a region that is ∼2° in radius. We confirmed the mean age of the region to be 4 Myr. We used 2MASS and WISE data to re-calculate the disk fraction in a self-consistent way for all sources. We obtain a disk fraction for all known members of 39%, consistent with previous values excluding anemic and transitional disks. Including the newly-identified members, the disk fraction is 28%, lower than previous studies, since our method detects mainly more massive and less embedded sources.

The color-magnitude diagrams show that most new members are intermediate- to solar-mass stars. This is reasonable since previous studies have identified members using techniques biased towards the lowest-mass members and confirm that Gaia data can help complete the stellar content of young regions for objects that are otherwise more difficult to identify.

Our results reveal a clear velocity structure in IC1396 on scales up to Δr <  ∼12–14 pc (which is approximately the edge of the ionized cloud). The expansion is a sign that the region is not gravitationally bound. We do not find any velocity structure on scales ∼12 <  Δr <  30 pc, as expected from the velocity dispersion in the region, ∼2–3 km s⁻¹, for a cluster age of ∼4 Myr. The velocity structure in the region, along with the spatial distribution of all members, also suggests different episodes of star formation.

We find differences between the ages of the subclusters, with the populations in the outskirts of subclusters B and E being older (on average) than subcluster A, except for the subcluster B objects spatially related to IC1396N. As found by (Sicilia-Aguilar et al. 2019), we confirm that the members of IC1396A do not show any significant difference in proper motion compared to the Tr37 cluster, which conforms most of the population of subcluster A. The same lack of differences in proper motion between the older objects and the globule population is observed in subclusters B and E. This is counterintuitive considering some triggered star formation models (e.g., Mintz et al. 2021) and indicates that even physically related groups may experience several star-forming events. Nevertheless, one limitation to further delving into this issue is that most objects physically associated with the globules are too embedded and not detected by Gaia.

The kinematics, age, and evolutionary state differences in groups of sources within the same cloud suggest that IC1396 has suffered more than one-star formation episode over time from various mechanisms, although radial velocity differences mean that triggering is not obvious (Sicilia-Aguilar et al. 2019). This, together with the spatial distribution of their members, indicates multiple episodes of star formation, which may be indicative of other processes such as coalescence (Smilgys & Bonnell 2017). The importance of Gaia data to complete cluster surveys, including intermediate-mass stars, is also demonstrated.

Acknowledgments

We thank the anonymous referee for providing valuable comments that helped to clarify this paper. M.P.B. is supported by the funding of the Peruvian National Council of Science, Technology and Technological Innovation (“Consejo Nacional de Ciencia, Tecnología e Innovación Tecnológica-CONCYTEC)”) via a predoctoral grant, contract “033-2016-FONDECYT (Peru)”. We thank Konstantin Getman for his valuable comments and for providing his unpublished list of X-ray sources for complete the collection of known members. Observations reported here include data obtained at the MMT Observatory, a joint facility of the University of Arizona and the Smithsonian Institution, and observations made with the Gran Telescopio Canarias (GTC), installed at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, on the island of La Palma. This work is (partly) based on data obtained with the instrument OSIRIS, built by a Consortium led by the Instituto de Astrofísica de Canarias in collaboration with the Instituto de Astronomía of the Universidad Autónoma de México. OSIRIS was funded by GRANTECAN and the National Plan of Astronomy and Astrophysics of the Spanish Government. This work is based on observations collected at the Centro Astronómico Hispano-Alemán (CAHA) at Calar Alto, operated jointly by Junta de Andalucía and Consejo Superior de Investigaciones Científicas (IAA-CSIC). 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 the DPAC has been provided by national institutions, in particular, the institutions participating in the Gaia Multilateral Agreement; 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, and 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. We also have made use of data products of the SIMBAD database, of the VizieR catalog access tool, operated at CDS, Strasbourg, France, and of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018).

References

Appendix A: Collection of known members

Table A.1 shows the information on known members of IC1396 obtained from the literature.

Appendix B: Spectroscopic members

Table B.1 shows the 121 spectroscopic candidates from MMT and GTC telescope observations, the results of the analysis of the different membership criteria, and the final membership derived from them all.

Appendix C: Maximum likelihood function

The presence of subclusters in the region (Sect. 4.2) was studied in a way similar to Roccatagliata et al. (2018, 2020), using the complete astrometric data for each object, including their uncertainties and correlated errors, and estimating which subcluster distribution and properties maximize the likelihood function. In this appendix, we include the details of the terms involved in the calculation of the individual likelihood function, L_i, j.

Each i-th star is defined by its astrometric data, which includes the trigonometric parallax (ϖ_i) and the proper motion in right ascension and declination (μ_α, i, μ_δ, i). Each quantity has its associated uncertainty, but because the uncertainties are not independent, we also have the cross-correlated uncertainties, ρ_i. In the multiparameter space, the position of star i will be represented by the vector a_i=[ϖ_i, μ_α, i, μ_δ, i] and each j-th subcluster by the vector a_j=[ϖ_j, μ_α, j, μ_δ, j], where ϖ_j, μ_α, j and μ_δ, j are the astrometric positions for subcluster j.

In Eq. (2), a_ij=(a_i-a_j) is the vector that represents the difference between the astrometric parameter of the star vs each subcluster, and (a_i-a_j)^T is the transposed vector. Therefore,

$$\begin{array}{c}\hfill a_i-a_j=\left[\begin{array}{c}{\varpi }_i-{\varpi }_j\\ {\mu }_{\alpha ,i}-{\mu }_{\alpha ,j}\\ {\mu }_{\delta ,i}-{\mu }_{\delta ,j}\\ \end{array}\right]\end{array}$$(C.1)

Following the formulation from Lindegren et al. (2000), the analog of the Gaussian standard deviation in the multivariate Gaussian space is given by the covariance matrix C_i, j (see Eq. (2). The covariance matrix contains the standard uncertainties derived from the observational data, which will be σ_ϖ, i, σ_μα, i and σ_μδ, i for the errors associated to the parallax and proper motion measurements for the i-th star. For each j-th subcluster, σ_ϖ, j, σ_μα, j and σ_μδ, j are the intrinsic dispersion of the parallax and proper motions. These, together with the correlated uncertainties between each pair of quantities, ρ_i(ϖ, μ_α), ρ_i(ϖ, μ_δ) and ρ_i(μ_α, μ_δ), define the covariance matrix C_i, j given by

$$\begin{array}{c}\hfill C_{i,j}=\left[\begin{array}{ccc}{C}_{ij,11}& C_{ij,12}& C_{ij,13}\\ C_{ij,21}& C_{ij,22}& C_{ij,23}\\ C_{ij,31}& C_{ij,32}& C_{ij,33}\\ \end{array}\right]\end{array}$$(C.2)

where the various terms correspond to

$$\begin{array}{cc}\hfill [C_{ij,11}]& ={\sigma }_{\varpi ,i}^2+{\sigma }_{\varpi ,j}^2\hfill \\ \hfill [C_{ij,22}]& ={\sigma }_{{\mu }_{\alpha },i}^2+{\sigma }_{{\mu }_{\alpha },j}^2\hfill \\ \hfill [C_{ij,33}]& ={\sigma }_{{\mu }_{\delta },i}^2+{\sigma }_{{\mu }_{\delta },j}^2\hfill \\ \hfill [C_{ij,12}]& =[C_{ij,21}]={\sigma }_{\varpi ,i}\times {\sigma }_{{\mu }_{\alpha },i}\times {\rho }_i(\varpi ,{\mu }_{\alpha })\hfill \\ \hfill [C_{ij,13}]& =[C_{ij,31}]={\sigma }_{\varpi ,i}\times {\sigma }_{{\mu }_{\delta },i}\times {\rho }_i(\varpi ,{\mu }_{\delta })\hfill \\ \hfill [C_{ij,23}]& =[C_{ij,32}]={\sigma }_{{\mu }_{\alpha },i}\times {\sigma }_{{\mu }_{\delta },i}\times {\rho }_i({\mu }_{\alpha },{\mu }_{\delta }).\hfill \end{array}$$(C.3)

The astrometric data, uncertainties, and correlated uncertainties for each individual object are provided by Gaia Archive, while the values for the subclusters are derived from the standard deviation of subcluster members during the iterations.

Starting with the first-guess subcluster parameters, we derive the likelihood using as initial positions those initial ones shifted along a grid to cover a large parameter space. The grid points are defined to sample the parameter space within 2σ of the initial values. The subcluster values are refined on each loop by re-deriving it again using the data from the stars that are found to belong to a cluster with a probability > 80% to calculate a final likelihood for each refined grid point. The final subcluster positions will be those that achieve the maximum likelihood.

Appendix D: Table of maximum likelihood results

Table D.1 shows the probabilities of the known members to belong to each subcluster, resulting from the maximum likelihood analysis from the Gaia EDR3 data. Objects with probabilities > 80% are considered members of one of the subclusters. Sources with a lower probability are probable members of one of the subclusters. Sources with probability > 40% for two subclusters are considered probable members of both subclusters. Objects with very small probabilities may still be members of the region but may belong to the disperse population, and they are labeled as members that do not belong to any subcluster.

Appendix E: New members table

Table E.1 shows the new members found using the Mahalanobis distance method. These new members have a probability of 95% or higher to belong to one of the four reliable subclusters, A, B, E, and F, and members with lower probabilities of belonging to two subclusters.

Appendix F: Three color images of the BRCs

In this appendix, we include some detailed images of the star-cloud interactions discussed in Sect. 5.3. Figures F.1, F.2, and F.3 show evidence of feedback in the BRCs IC1396A, IC1396N, and IC1396B, respectively.

Fig. F.1. Three-color images (Ks, H2, and S[II] as red, green, and blue, respectively) of the BRCs. This figure shows IC1396A (top, zoomed on the bottom). The white dotted-line box is zoomed at the bottom. The shocks and knots are visible in H2 and S[II] as green and blue structures and marked with white arrows. Dotted white arrows show the trajectory of jets following the shocks or knots (see text Sect. 5.3). The remaining globules are shown in Figs. F.2 and F.3

Fig. F.2. continued. IC1396N (top, zoomed on the bottom). The strong shocks are numbered from 1-4. Sources B1, B2, and B3 from Beltrán et al. (2009) are located in the middle of three chains of knots. Red circles enclose sources C, I, and 331A (21:40:41.4 +58:16:38).

Fig. F.3. continued. Shocked material on the small irregular edges of IC1396B. A source is labeled as a coordinate reference.

All Tables

All Figures

Fig. 1. WISE W4 (22.19 μm) image of the IC1396 H II region. The BRCs are labeled with their names. The red boxes and yellow circles enclose the fields observed by OMEGA2000 and CAFOS, respectively. Note: CAFOS FOV has been reduced to account for the fact that the Fabry-Perot etalon has only a small area where the wavelength is accurate. The HD206267 star that is ionizing the region is indicated with a yellow star.

In the text

	Fig. 2. MMT (Hectospec) spectra of two confirmed members. The Hα emission line is cut off at the top for better visualization. The fluxes are given in arbitrary units.
In the text

Fig. 3. Distribution of the extinction for the spectroscopically confirmed members (top). The mean error (0.3 mag) is shown with a red bar. Color-magnitude diagram for candidates with spectroscopic signatures of youth (bottom). Members (M) are marked by blue circles. Probable members (PM) are marked by red triangles. Probable non-members (PNM) are marked by gray circles. The number of objects in each class is also given. The magnitudes of each object are corrected by their individual extinction (full color symbols), except for objects with non-physical or uncertain extinctions, which have been corrected by cluster average extinction of 1.4 mag (open symbols). The PARSEC isochrones are corrected by the average distance of 925 pc.

In the text

	Fig. 4. W1 − W2 vs. H − K diagram for the 67 spectroscopic candidates with errors < 0.05 mag. The red box shows the 27 candidates with excess emission in the IR consistent with a disk.
In the text

Fig. 5. Parallax histogram (top, including the full histogram on the left and a zoom of the subclusters on the right), proper motion and color-magnitude diagram (bottom) of the known members of the subclusters found by the maximum likelihood analysis with the Gaia EDR3 data. The symbols and colors mark the various subclusters using circles, triangles, and capital letters (see legend). The number of subcluster members is given in parentheses. Full circles represent stars with a membership probability > 80%. Triangles represent objects with probabilities of < 80%. The theoretical PARSEC isochrones on the color-magnitude diagram have been corrected for a distance of 925 pc and a lower limit extinction AV = 1 mag to reject foreground stars. The extinction vector for 1 mag is indicated with a red arrow.

In the text

Fig. 6. Color-magnitude diagram for the new members obtained from Mahalanobis distance analysis with the Gaia EDR3 data. The 334 new members between 0.1 and 10 Myr are in full color, candidates with ages > 10 Myr are rejected and shown in semi-transparent color. The theoretical PMS tracks have been corrected for a distance of 925 pc and a lower limit of extinction of AV = 1 mag to reject foreground stars. Five members are plotted as triangles and belong to both subclusters B and E. We note there is an empty space between the 10 Myr and 15 Myr isochrones, which reinforces the criterion of considering as members only those up to 10 Myr. A vector extinction is indicated by a red arrow.

In the text

Fig. 7. Spatial distribution of the known (full-color circles) and new (pentagon) members. The colors represent the four subclusters (A, B, E, F, see legend). Orange circles mark the rest of the known members not belonging to any subcluster. The age of members (derived in Sect. 5.1) is represented by the symbol sizes, with the youngest members having large symbols ranging from 0.1 to 20 Myr. Black numbers indicate the main BRCs, 1=IC1396 A, 2=IC1396 B, 3=IC1396 D, 4=IC1396 N, 5=IC1396 G. The background image is a ∼4° × ∼4° mosaic built from WISE 22.19 μm.

In the text

Fig. 8. J − H vs. H − K diagrams for known members associated with the various subclusters obtained from the Maximum likelihood analysis (colored circles) and new members obtained from the Mahalanobis distance analysis (magenta pentagons). The number of sources with complete JHK photometry is indicated in the legend. Subclusters B and E include the new members with an intermediate probability of belonging to both subclusters. The MS and the giant branch are indicated by thin solid lines (Bessell & Brett 1988). The CTTS locus (Meyer et al. 1997) is indicated by a thick solid line. The theoretical tracks are corrected for extinction using the threshold AV = 1 mag. A reddening vector for AV = 3 mag is shown as a red arrow.

In the text

Fig. 9. Map of the velocity dispersion in proper motion of known and new members. The arrows mark the direction and size of the proper motion vectors relative to the average proper motions of subcluster A. They are color-coded by subcluster (see legend). Known members not associated with any subcluster are shown in orange. The long orange arrows are mostly from the extended population G. New members are marked by open circles. The magenta arrow in the bottom left corner indicates the physical velocity scale (at 925 pc).

In the text

Fig. 10. Velocity structure in the IC1396 region, shown in terms of the mean difference between the velocity magnitude ΔvM (left panels) and the directional velocity ΔvD (right panels) versus the difference distance in parsec (Δr; see text). The two top panels show the region velocity structure considering all known and new members, including the extended population (G), which has considerably larger velocities and velocity dispersion than the rest of clustered members. The two bottom panels include only the clustered part of the population, excluding subcluster G, and thus offer a better view of the evolution and dispersal of the clustered members.

In the text

Fig. 11. Age (left) and mass (right) distribution of known members (top panels) and new members (bottom panels) from different subclusters. We have included those known members that do not belong to any subcluster and those that belong to one more subcluster (labeled as others “O” in the legend). Some new members belong to two subclusters, marked as “B/E”. Objects with ages 20 Myr and over are all included in the last bin.

In the text

Fig. 12. J − H vs. H − K diagram for all known members (KM) with magnitude uncertainties < 0.05 mag. The MS and the giant branch (Bessell & Brett 1988), as well as CTTS locus (Meyer et al. 1997), are indicated by a solid line. The theoretical tracks are corrected for extinction in JHKs (Cardelli et al. 1989), using the threshold AV = 1 mag. The extinction vector is indicated with a red arrow. KM in full color belong to the populations (A, B, C, D, E, F, and G) found from the maximum-likelihood analysis.

In the text

	Fig. 13. W1 − W2 vs. H − K diagram for all members with magnitude errors < 0.05 mag. Known members (KM) are shown in green and new members (NM) in magenta. The green and magenta boxes enclose the known and new members, respectively, with color excesses consistent with disk emission.
In the text

Fig. F.1. Three-color images (Ks, H2, and S[II] as red, green, and blue, respectively) of the BRCs. This figure shows IC1396A (top, zoomed on the bottom). The white dotted-line box is zoomed at the bottom. The shocks and knots are visible in H2 and S[II] as green and blue structures and marked with white arrows. Dotted white arrows show the trajectory of jets following the shocks or knots (see text Sect. 5.3). The remaining globules are shown in Figs. F.2 and F.3

In the text

	Fig. F.2. continued. IC1396N (top, zoomed on the bottom). The strong shocks are numbered from 1-4. Sources B1, B2, and B3 from Beltrán et al. (2009) are located in the middle of three chains of knots. Red circles enclose sources C, I, and 331A (21:40:41.4 +58:16:38).
In the text

	Fig. F.3. continued. Shocked material on the small irregular edges of IC1396B. A source is labeled as a coordinate reference.
In the text