Interplay between Young Stars and Molecular Clouds in the Ophiuchus Star-forming Complex

We present spatial and kinematic correlation between the young stellar population and the cloud clumps in the Ophiuchus star-forming region. The stellar sample consists of known young objects at various evolutionary stages, taken from the literature, some of which are diagnosed with Gaia EDR3 parallax and proper-motion measurements. The molecular gas is traced by the 850 $\mu$m Submillimetre Common-User Bolometer Array-2 image, reaching $\sim$2.3 mJy beam$^{-1}$, the deepest so far for the region, stacked from the James Clerk Maxwell Telescope/Transient program aiming to detect submillimeter outburst events. Our analysis indicates that the more evolved sources, namely the class II and III young stars, are located further away from clouds than class I and flat-spectrum sources that have ample circumstellar matter and are closely associated with natal clouds. Particularly the class II and III population is found to exhibit a structured spatial distribution indicative of passage of shock fronts from the nearby Sco-Cen OB association thereby compressing clouds to trigger star formation, with the latest starbirth episode occurring now in the densest cloud filaments. The young stars at all evolutionary stages share similar kinematics. This suggests that the stellar patterns trace the relics of parental cloud filaments that now have been dispersed.


Introduction
It is likely that all stars are formed in groups out of giant molecular clouds. A star group may soon be dispersed upon emerging from the natal cloud (Lada & Lada 2003). Subsequent gas dispersal, two-body relaxation between member stars, and stellar ejection exacerbate the disassembly process, leaving behind a star cluster we are witnessing now. With time, external perturbations, such as tidal forces plus differential rotation, continue to dissolve the star cluster to supply the disk field stars.
While the formation of individual stars in observation or in theory is reasonably understood, particularly for low-mass stars, certain details of stellar birth remains unclear. For instance, are stars of different masses formed coevally? Do luminous members tend to be centrally located in a star cluster as their birthplace, or are they dynamically settled with time? How do these short-lived massive stars influence the subsequent star formation on nearby clouds? While dense molecular cores, within which individual stars form, appear to have already mass segregated , tidal disruption or dynamical ejection could lead to certain spatial distribution of protostars and pre-main-sequence (PMS) stars versus cloud filaments (Stutz & Gould 2016). To address questions concerning, e.g., the intricate interplay of massive stars, molecular clouds, and young stars, a sample of the youngest stellar population recently exposed out of, but still closely associated with, molecular clouds is needed.
The Ophiuchus molecular cloud complex offers such a target. With the proximity of ∼138 pc (Ortiz-León et al. 2018), it is one of the nearest known star-forming regions (Wilking et al. 2008) to render optimal sensitivity and angular resolution to diagnose the cloud morphology and, with a wealth of optical and infrared observations, the stellar kinematics. Particularly, the densest cloud core in the region, known as L1688, is considered the nearest site of cluster formation (Wilking & Lada 1983). The region is known to be influenced, or even with the starbirth triggered at least partly, by the nearby Sco-Cen OB (Sco OB2) association (see Wilking et al. 2008), with the cloud "streamers" as the leftover material of star formation (Vrba 1977). With an age spread of 0.3 Myr in the core (Chen et al. 1995;Wilking et al. 2008) to 2-4 Myr for the rest of L1688 (Esplin & Luhman 2020), the young stellar and substellar populations in L1688 represents a comparative sample at a variety of evolutionary stages (Dunham et al. 2015) from embedded protostars, PMS objects with dusty disks/envelopes, to relatively evolved PMS objects void of inner disks, and from aging massive stars to brown dwarfs in infancy, so as to shed light on the interaction and feedback between the cloud complex and the newborn stars.
This study uses the deepest 850 μm James Clerk Maxwell Telescope (JCMT)/Submillimetre Common-User Bolometer Array (SCUBA-2) image of L1688 to trace the cold dust, as a proxy of molecular cloud structure, so as to correlate spatially and kinematically with the sample of young stars in the region compiled from the literature (Wilking et al. 2008;Dunham et al. 2015;Allers & Liu 2020;Esplin & Luhman 2020), plus a few probable member candidates identified using the Gaia EDR3 parallax and proper-motion data (Gaia Collaboration et al. 2021). The PMS evolutionary status of these young stars, either already reported in the literature or quantified by us with archival photometric measurements across a range of wavelengths, then allows us to infer the star formation history in the region.

Data and Analysis
Data used in our study include submillimeter emission tracing cold dust in molecular clouds, and a sample of young stellar objects in L1688 collected from the literature. We are particularly interested in the spatial distribution and kinematics of young stars at different evolutionary stages relative to the cloud substructures.

JCMT 850 μm Clumps
To trace the molecular substructure, we used the 850 μm emission acquired with the SCUBA-2 (Holland et al. 2013) on the JCMT. SCUBA-2 is a 10,000 pixel bolometer taking continuum data simultaneously at 450 μm and 850 μm, with effective beam sizes of 9 8 and 14 6, respectively. The data adopted for this study were taken as a part of the JCMT Transient Survey (Herczeg et al. 2017;Johnstone et al. 2018), designed to detect continuum variability of young embedded sources utilizing multiepoch observations of star-forming regions, including L1688 in the Ophiuchus cloud. Raw data were reduced using MAKEMAP , an iterative map making software. Maps of individual epochs were then stacked together, thereby creating a deeper coadded image  for each region. For our study, we employed a coadd of 29 epochs of L1688 observations, reaching a sensitivity of ∼2.3 mJy beam −1 , about a factor of 5 lower than a typical single epoch observation, and roughly three times more sensitive than the JCMT Gould Belt Survey data (Pattle et al. 2015).
We identified enhanced submillimeter substructure, or "clumps" using the STARLINK implementation of the Fell-Walker algorithm (Berry 2015), which traces boundaries of clumps by following paths from low-valued to high-valued pixels in an iterative manner. Polygonal clump boundaries are identified and their peak fluxes and total fluxes are estimated. For our analysis, the key parameter for clump detection, the "noise," was set to be 10 mJy beam −1 . A total of 108 clumps were identified with the minimal peak flux being ∼30 mJy beam −1 .
Our study focuses on the spatial and kinematic distributions of YSOs relative to cloud clumps; no attempt was made to derive the temperature and mass of individual clumps. Figure 1 exhibits the latest coadded SCUBA-2 850 μm image of L1688, marked with clump boundaries. Our JCMT data are relevant to those reported by Johnstone et al. (2000) and by Pattle et al. (2015). Johnstone et al. (2000) used SCUBA, the predecessor of SCUBA-2, to cover a 700 arcmin 2 region, and with the clfind clump-finding algorithm (Williams et al. 1994), identified 55 clumps. These authors derived the clump mass by two models, one adopting a single temperature of 20 K for every clump, and the other by assuming the clump to be an isothermal, selfgravitating, pressure-confined Bonnor-Ebert sphere. The resultant clump mass functions are similar. Pattle et al. (2015) using SCUBA-2 to analyze L1688 along with parts of L1689, L1709, and L1712, identified 70 clumps in L1688 with the CuTEx algorithm (Molinari et al. 2011). Combining the 160 and 250 μm data from the Herschel mission (Ladjelate et al. 2020), plus the 450 μm data from SCUBA-2, Pattle et al. (2015) constructed the spectral energy distribution (SED) of each clump to estimate its temperature, and hence to derive the mass. The mass function of the clumps from a few solar masses down to ∼0.1 M e is found to have a slope consistent with that of the stellar initial mass function.

Young Stellar Objects
Our young stellar sample consists of a compilation of literature lists. The primary reference is the review article by Wilking et al. (2008), in which 238 objects are within the sky area of our study on the basis of stellar youth characteristics, such as infrared excess, lithium absorption, or X-ray emission. In terms of decreasing level of infrared excess as a proxy of PMS evolutionary sequence, from the youngest, embedded protostellar stage, through a transition phase, to the T Tauri stage, and finally to the inner-disk free stage, there are 23 class I, 28 flat-spectrum, 92 class II, and 36 class III sources. There are 59 sources without classification.
In addition, Dunham et al. (2015) identified young stars in the star-forming regions in the Gould Belt, including Ophiuchus, using data from the Spitzer Telescope (Werner et al. 2004), and found 14 additional sources, eight of which are class I or flat-spectrum sources. Toward lower masses, Allers & Liu (2020) reported the young substellar objects in Ophiuchus and Perseus, seven of which (four newly found) are within our field of view. Esplin & Luhman (2020) compiled a list of YSOs in the Ophiuchus cloud complex on the basis of spectroscopic evidence of youth; among these, 13 additional sources are included in our list.
We combined all these 269 sources as our young stellar sample. Of these, a total of 99 have Gaia EDR3 counterparts, each within a matching radius of 5″, with astrometric and photometric measurements down to the Gaia brightness limit of g ≈ 21 mag. Figure 2 presents the Gaia EDR3 parallax and proper-motion measurements. In each case, the young sample stands out distinctly. However, there are 10 sources (GSS 39, Elia 2−31, V * V2675 Oph, EM * SR 24A, ISO-Oph 9, DROXO 106, EM * SR 23, VSSG 22, BBRCG 54, BKLT J162748−242204) found not to be in the parallax range; i.e., they are likely foreground or background objects. Incidentally, we recognized four sources coexisting and comoving with our YSO sample. Two of them (BKLT J162749−242522, BKLT J162818−242836) are known PMS stars, whereas the nature of other two (BKLT J162553 −244129, BKLT J162702−241639) needs to be further vindicated; but, we included them in the member sample for now. At the end, we have a total of 263 YSO sources for analysis. Cánovas et al. (2019) and Grasser et al. (2021) conducted a similar search for young stars in the whole Ophiuchus cloud complex using Gaia DR2 and EDR3 data, respectively. They also identified BKLT J162702−241639 and BKLT J162818−242836 as cluster members. All the new sources they identified as candidate young stars within our field of L1688 are already in our final list. All the YSOs have Gaia EDR3 parallax consistent with membership. Because the distance estimate as the direct reciprocal of the parallax introduces statistical bias (Bailer-Jones et al. 2021), this work focuses only on the projected sky location of stars with no consideration of the distance distribution.

Kinematic Data
In order to investigate the kinematic properties of the cluster, we gathered the tangential velocities (proper motions) and radial velocities of the YSOs, when available. The Gaia EDR3 provides proper motions of 84 YSOs in our sample, incomplete for embedded sources. A literature search using Simbad led to 32 additional sources having proper-motion data in Ducourant et al. (2017) and one in Monet et al. (2003). The list in Ducourant et al. (2017), though greatly increasing the propermotion data in our sample, is limited mainly to the Ophiuchus F core and partially to the E core. For radial velocities, the Gaia EDR3 data are available for only six objects, whereas Sullivan et al. (2019) provide data for 31 sources, and Torres et al. (2006) for additional one source. Dunham et al. (2015) reported the evolutionary status of the YSOs in the region. We conducted an independent analysis so as to compare with the previous results and also to expand the sample to include sources with no classification in the literature. We compiled for each source its SED as widely covered in wavelength as possible from archival photometric data. These In our analysis, 30 SEDs have inferior fits (σ α > 0.5) mainly arising from imperfect cross-matching among database catalogs or a lack of data in certain wavelength ranges, so were not included in further analysis. Even though our α index is not extinction corrected, of the 132 sources common with Dunham et al. (2015), 80% have consistent classifications, which is considered satisfactory given the systematic uncertainties due to inclination or aspherical geometries of the disks/envelopes, etc. For the few sources we did not have our own classifications, we adopted the classes from Dunham et al. (2015) or Wilking et al. (2008). At the end, our YSO sample consists of 30 Class 0/I sources, 56 flat-spectrum sources, 110 Class II sources, and 58 Class III sources, of which ∼21% of the sources had no classification prior to our analysis. Some ∼3% sources still remain unclassified.

Results and Discussion
In addition to distinct cloud clumps, the deep JCMT submillimeter data presented as Figure 1 reveal clearly the known filamentary morphology of the clouds, with cloud Oph-B manifest the "streamers" shaped by the passing shocks from the OB stars to the west and southwest, including but not limited to HD 147889 (B2 III), σ Sco (O9.5 V + B7 V), and α Sco (M0.5 Iab + B3 V). The extension to the northeast immediately from Oph-A may be a part of the ring/shell created by ρ Oph (B2 IV + B2 V). All these luminous stars are outside the field of our JCMT data and the readers may refer to Wilking et al. (2008) and Nozawa et al. (1991), both with a 13 CO map, or Ladjelate et al. (2020) with a Herschel map for a larger-scale illustration of the complex. Motte et al. (1998) noticed the alignment of protostars along clouds Oph-C, Oph-E, and Oph-F. This direction of clump filaments (see also Johnstone et al. 2000), orthogonal to the steamers, arches to A, which probably acts as the "hub," i,.e., the conjunction of cloud filaments, along which mass flows to seed the formation of massive stars or star clusters (Kumar et al. 2021). The scenario is consistent with the magnetic field direction, inferred by optical and infrared polarization (Vrba et al. 1976;Vrba 1977) as a consequence of dichroic extinction due to grain alignment, that by and large parallels the dark cloud filaments.
The deep JCMT/SCUBA-2 image uncovers additional small clumps, for the first time to our knowledge, to the southwest of the main filaments A-C-E-F, forming a pattern that seems to loop between Oph-A and Oph-E/F, stretching to R.A. ∼ 16 h 26 m , decl. ∼ −24°32′ where it is recognized previously as a part of cloud Oph-J (Simpson et al. 2008). This string of discrete cloud clumps is marked as the purple dashed line in Figures 1 and 3 and will be discussed later in the context of the YSO distribution. In each case, the gray symbol represents every Gaia star in the region, whereas the pink symbol marks known YSOs. In (a), the thick and thin dashed lines mark, respectively, the median and interval of five times median-absolute deviation. In (b), the ellipse is the region of five times the median-absolute deviation. Figure 3 displays the spatial distribution of the YSOs with respect to the 850 μm emission. It is seen that the class I and flat-spectrum sources are preferentially nearer the molecular clumped structure than more evolved sources are. This is consistent with the previous investigation in the same region by observations (Jørgensen et al. 2008), and by numerical simulation (Frimann et al. 2016).

Stellar Spatial Distribution
The result is quantified in Figure 4 for which the projected linear separations of YSOs to the center of the nearest 850 μm clump, adopting a distance of 138 pc to the complex, are presented for different evolutionary classes. While more evolved sources are scattered equally likely up to 0.2 pc from a cloud, objects in infancy are an order closer in the vicinity of clouds, mostly within 0.02 pc. Also exhibited is the same separation distribution but for the sample of Dunham et al. (2015), which is smaller in size than ours but with reddening corrections. The distributions are similar for class O/I and flatspectrum sources, but contrast mostly for the class III samples. Conceivably, a highly reddened class III object could be mistaken with an α index as one with significant IR excess. But the reverse is not true; that is, a genuine class I/O source would not have been recognized as an evolved class III source even without dereddening. The conclusion that less evolved YSOs, either with our sample or with that of Dunham et al. (2015), are physically more associated with dense gas clumps therefore remains valid. Mairs et al. (2016) and Stutz & Gould (2016) found similar trends of younger stars being more closely associated with molecular clouds in the Orion A star-forming region, roughly 10 times in size scales of our field. Under the assumption that the YSOs were formed at the location of the presently observed cloud clumps, these authors attributed this trend to greater relative velocities of more evolved sources so as to being decoupled from parental molecular gas.
While clouds in L1688 have a hub-filamentary morphology, the young stars in the region also exhibit a structured spatial distribution. As discussed earlier, the class I and flat-spectrum sources are clustered around the cloud emission and thus follow the distribution of the clouds. The class II/III sources, though correlated with generally weaker gas emissions, appear to form arc-like patterns. Most prominent is the western partial loop between cloud Oph-A and clouds Oph-E/F connecting the 850 μm clumps (marked as a purple dashed line in Figure 3) that are associated with mostly flat-spectrum and class II sources. Further to the southwest as the outermost arc (marked as red dashed line in Figure 3) near the edge of our JCMT field are class III young stars with no corresponding submilimeter emission, thereby to our knowledge not recognized before.
To the other side of the main cloud filaments (A-C-E-F), between Oph-B and Oph-F additional stellar patterns appear to stand out, approximately in the NS of NE-SW direction. It is not clear whether these are parts of the cloud streamers or of the southwest YSO loops which are physically closer in 3d to the OB stars.
Such organized patterns could not have been the outcome of random ejection from common sites of formation. In order to test the randomness in the spatial distribution of the evolved sources, we used pointpats implementation of Ripley's K function (Ripley 1976), which is a proxy of percentage of points (in this case class II/III sources) that are within a certain distance of each other. It is commonly used to quantify spatial randomness for discrete location data, i.e., applicable to our case to address whether YSOs are distributed randomly or in a clustered pattern (see Retter et al. 2019). Figure 5 compares the K function of our YSO sample against that of an otherwise complete spatial randomness via Monte Carlo simulation. The observed K function lies outside the envelope at all angular scales considered here (less than ¢ 3.5), leading to to a p value of 0.0001, thereby rejecting the null hypothesis of being a random distribution. We now present the motion of the YSOs relative to clouds based on available kinematic data, as discussed in Section 2.2.1. The peculiar motion of YSOs, i.e., the proper motion of individual sources subtracted by the median proper motion, is shown in Figure 6. Using these relative velocities, we inspected the kinematic distribution for different groups of YSOs. Figure 7 presents the tangential velocity distributions of YSOs at different evolutionary stages. Lacking proper-motion data for highly embedded sources, we combined Class 0/I sources and flat-spectrum sources together. One sees that these least evolved population has a spread of velocities, partly perhaps because of the uncertainties in the data, whereas the more evolved sources, class II and class III alike, consistently move slowly. This supports the notion that the different YSO classes have similar velocities, and more evolved sources are further away from the clouds not because they are moving faster. Motte et al. (1998) presented evidence, with the cloud morphology, kinematics, and magnetic field configuration, of multiepoch triggered star formation in ρ Oph by a combination of plausible mechanisms, such as the shocks from Sco-Cen OB association by supernova explosion, massive stellar winds, or ionization fronts.

Stellar Kinematic Distribution
The two southwestern arcs (or partial layers in 3d) of YSOs reported in this work provide morphological evidence of the direction of the star-forming sequence propagating from the southwest now to the dense segment extending from A to F. Kinematic analysis indicates that all the YSOs remain near their formation groups and hence retain the spatial structures of their natal clouds. The earlier episode of star formation to the southwest is imprinted in the YSOs as relic of the filamentary or sheet-like clouds. The impact of massive stars occurred at an    Figure 3. Each arrow marks the direction of the proper motion, with the length proportional to the speed. The blue arrows represent Gaia measurements, whereas the maroon arrows are for values from other data sets (see Section 2.2.1). The scale reference on the top right marks 10 mas yr −1 , which, in accord with the coordinate axes, corresponds to a motion of ¢ 6 in 3.5 × 10 4 yr.
earlier epoch on the chain of class III sources that are more evolved or depleted much of the circumstellar disks and interstellar cloud. White et al. (2015), using the C 18 O J = 3-2 emission, estimated a total molecular gas mass to be ∼515 M e , consistent with those derived from other lines (e.g., 13 CO J = 1-0; Loren 1989; or C 18 O J = 1-0; Tachihara et al. 2000). White et al. (2015) measured a line width of 1.5 km s −1 , corresponding to an isotropic one-dimensional velocity dispersion of ∼0.64 km s −1 . This value is comparable to the line-of-sight velocity dispersion between cores in Perseus (Kirk et al. 2007). Under the same assumption of virial equilibrium as in White et al. (2015), i.e., s g = M G R 5 vir 2 , where σ vir is the velocity dispersion, M is the mass of molecular gas, R is the radius of the cloud, γ is the density profile factor (taken here to be 5/3 for a uniform sphere), and G is the gravitational constant, the threedimensional velocity dispersion with the adopted distance of 138 pc, assuming spherical symmetry, would be 0.96 km s −1 .
In comparison, our young stellar group has a threedimensional velocity dispersion, each represented by fitting with a Gaussian, to be 2.42 ± 0.33 km s −1 along the line of sight, 1.28 ± 0.09 km s −1 in the R.A. direction, and 1.11 ± 0.07 km s −1 in the decl. direction. The latter two, both in the sky plane, agree with each other, whereas the higher value and dispersion of the radial velocity can be attributed to larger observational uncertainties and possible stellar binarity. In any case, the values suggest a velocity dispersion of the young stars more than that of the natal gas, and also marginally greater than that expected for a system in virial equilibrium. For a sheet-like configuration stretching along our line of sight, the actual surface mass density, while uncertain, is largely overestimated by the measured column density, hence leading to a reduced velocity dispersion.
This increase of stellar velocity can be the consequence of dynamical evolution of the young cluster. Simulations of dynamics of young embedded clusters (e.g., Proszkow et al. 2009) demonstrate that as a system evolves from a subvirial state, the individual members tend to fall toward the center of the potential well and gain kinetic energy leading to a boosted velocity dispersion. Alternatively, the increased dispersion can also be due to the dispersal of gas by stellar feedback, leaving behind a shallower potential well and thus, a loosely bound supervirial cluster.
Moreover, as a star cluster ages, two-body relaxation among stellar members leads to mass segregation, i.e., lower-mass members would occupy progressively larger volume in space, vulnerable particularly for the lowest ones to subsequent ejection from the system, whereas the massive stars "sink" to the center. Because mass information is not readily available for our YSO objects, we investigated the kinematics on the spectral type of the YSOs listed in Allers & Liu (2020) and Esplin & Luhman (2020), as presented in Figure 8. YSOs earlier than the M type, while low in numbers, appear to have consistently tangential speeds up to 3 km s −1 including the two Herbig Ae stars, SR 3 and IRAS 16245−2423. We note that the star IRAS 16245−2423 (YLS 46, or 2MASS J16273718−2430350) has uncertain spectral typing; Erickson et al. (2011), measuring a range of B5 to F2, adopted an A0 type. Those of M types or later seem to have noticeably a wider dispersion (up to ∼7 km s −1 ), notwithstanding the larger errors associated with non-Gaia proper-motion measurements Ducourant et al. 2017) for the faint, late-type sources. All the Gaia data, however, indicate a consistent tangential velocity regardless of the spectral type. Given that only ∼17% early-type stars have tangential speeds greater than ∼3 km s −1 , the overall kinematics do not seem to vary greatly across different spectral types. Better astrometric data for the faint early-type sources is needed for further analysis.

Conclusions
Using high-precision Gaia EDR3 proper-motion measurements of the young stars in L1688 of the Ophiuchus starforming region, incorporating deep 850 μm JCMT/SCUBA2 image, in a field of ∼1300 square arcmin, we diagnose the interplay between starbirth and natal clouds to trace the star formation history in the complex. Our main results are: 1. With the JCMT/SCUBA-2 850 μm image, we derive the position and boundary of molecular clumps. The image, stacked from the JCMT Transient Survey aiming to detect possible brightness outbursts related to young stellar accretion, or magnetic activity, reaches ∼2.3 mJy per beam, deeper than previously reported in the literature. 2. The class I and flat-spectrum YSOs in the region are found to be closely associated with molecular clouds. The class II and class III are located preferentially further away from clouds, and are not randomly distributed in position; instead these young stars with aged disks show a structured spatial pattern that is devoid of molecular gas. 3. Stellar kinematics suggests Class II sources to have a comparable velocity dispersion as the Class III sources do. While data are of less precision for Class I and flat-spectrum sources, they do not seem to move slower in space. This suggests that the YSO populations are comoving in space and different evolutionary classes serve to diagnose the star formation history in the region. 4. The young stars exhibit segregated spatial distribution, with a loop of the most evolved young stars (class III), to those younger and associated with a part of cloud Oph-J, to the youngest near ongoing star-forming sites in the main cloud complex. This is consistent with a compression front from the OB stars located to the west and southwest.