Abstract
High-mass multiples might form via fragmentation of self-gravitational disks or alternative scenarios such as disk-assisted capture. However, only a few observational constraints exist on the architecture and disk structure of high-mass protobinaries and their accretion properties. Here, we report the discovery of a close (57.9 ± 0.2 mas = 170 au) high-mass protobinary, IRAS17216-3801, where our VLTI/GRAVITY+AMBER near-infrared interferometry allows us to image the circumstellar disks around the individual components with ∼3 mas resolution. We estimate the component masses to ∼20 and ∼18 M⊙ and find that the radial intensity profiles can be reproduced with an irradiated disk model, where the inner regions are excavated of dust, likely tracing the dust sublimation region in these disks. The circumstellar disks are strongly misaligned with respect to the binary separation vector, which indicates that the tidal forces did not have time to realign the disks, pointing toward a young dynamical age of the system. We constrain the distribution of the Brγ and CO-emitting gas using VLTI/GRAVITY spectro-interferometry and VLT/CRIRES spectro-astrometry and find that the secondary is accreting at a higher rate than the primary. VLT/NACO imaging shows L'-band emission on (3–4)× larger scales than the binary separation, matching the expected dynamical truncation radius for the circumbinary disk. The IRAS17216-3801 system is ∼3× more massive and ∼5× more compact than other high-mass multiplies imaged at infrared wavelength and the first high-mass protobinary system where circumstellar and circumbinary dust disks could be spatially resolved. This opens exciting new opportunities for studying star–disk interactions and the role of multiplicity in high-mass star formation.
Export citation and abstract BibTeX RIS
1. Introduction
There is now solid evidence that the formation of high-mass stars (>10 M⊙) involves accretion through circumstellar disks (Cesaroni et al. 2005; Kraus et al. 2010; Boley et al. 2013; Johnston et al. 2015). This observational evidence supports recent theoretical work that suggests that the radiation pressure barrier problem can be overcome when considering more complex than spherically symmetric infall geometries (Yorke & Sonnhalter 2002; Krumholz et al. 2009). However, the classical accretion disk scenario, which assumes the monolithic collapse of a core to a single star and which is commonly applied to low-mass star formation, fails to explain the high multiplicity fraction that has been observed for high-mass stars. Multiplicity studies on main-sequence stars found that ≳80% of all O-type stars (M ≳ 16 M⊙) are close multiple systems, while this fraction rapidly drops to 20% for stellar masses of ∼3 M⊙ (Chini et al. 2012). At the same time, the number of companions per system also increases with stellar mass. For instance, the five O-/early B-type stars in the Orion Trapezium cluster have on average 2.5 known companions, which is about 5× higher than what has been found for low-mass stars (Grellmann et al. 2013). Different mechanisms have been proposed in order to explain these extraordinary characteristics. For instance, these high-mass multiples might form via fragmentation of self-gravitating massive disks (Kratter & Matzner 2006), disk-assisted capture (Clarke & Pringle 1991; Bally & Zinnecker 2005), or through failed mergers in stellar collisions (Dale & Davies 2006). Testing these scenarios requires the detection of high-mass "protobinaries" that are still in their formation phase.
Here, we report high angular resolution observations that resolve the high-mass protostar IRAS17216-3801 into a close (∼170 au) binary, where both components are associated with actively accreting circumstellar disks. IRAS17216-3801 (=G350.011–01.341) was identified as a high-mass young stellar object (YSO) by Persson & Campbell (1987) and its distance estimated to 3.08 ± 0.6 kpc (Boley et al. 2013). Various studies reported the detection of OH masers (Cohen et al. 1988; Caswell 1998; Argon et al. 2000). Boley et al. (2013) resolved the mid-infrared (8–13 μm) emitting region and found an elongated structure with a Gaussian FWHM size of 
milliarcseconds (mas) and speculated that the emitting structure might be associated with a circumstellar disk or an outflow.
We observed the source as part of a small survey on high-mass YSOs that we conducted with the Very Large Telescope Interferometer (VLTI) and its AMBER beam combiner. The AMBER measurements with low spectral resolution (R = λ/Δλ= 35) were complemented with higher spectral resolution data (R = 500) from the newly commissioned VLTI/GRAVITY combiner.
2. Observations
2.1. VLTI/AMBER Spectro-interferometry
We observed IRAS17216-3801 on 2012 May 6, 2012 June 1, and 2012 June 3 using the AMBER instrument (Petrov et al. 2007). The observations used the UT1-UT2-UT4 and UT2-UT3-UT4 triplets of the VLTI array of four 8.2 m unit telescopes, with projected baseline lengths of 30...130 m (Figure 1, top left).
Figure 1. Spectro-interferometry on IRAS17216-3801: obtained uv-coverage (top left) and the measured visibility amplitudes, plotted as a function of spatial frequency projected perpendicular (top right, upper panel) and parallel to the binary separation vector (top right, lower panel). Bottom: aperture synthesis images obtained with GRAVITY (bottom left) and GRAVITY+AMBER (bottom right).
Download figure:
Standard image High-resolution imageOur setup covered the K-band with R = 35, where we employed a 50 ms integration time. The science star observations were bracketed with observations on calibrator stars with known diameters, which allowed us to correct for atmospheric and instrumental effects. As calibrators we selected HD155259 (2012 May 6, 2012 June 1, UDD = 0.226 ± 0.016 mas) and HD165787 (2012 May 6, 2012 June 3, UDD = 0.49 ± 0.04 mas). The uniform disk diameters (UDDs) were estimated using the JMMC SearchCal tool (Bonneau et al. 2011).
The visibilities and closure phases were extracted using the amdlib software (Release 3.0.4; Tatulli et al. 2007; Chelli et al. 2009). The resulting amplitudes/phases show strong wavelength-differential modulations (Figure 2, bottom), which is characteristic for wide-separation multiple systems. Furthermore, we noticed that the observables, in particular, the closure phases, changed slightly between individual object exposures. Therefore, we decided to analyze the individual exposures (∼50 s) separately in order to avoid temporal smearing. We also checked whether our data might be affected by bandwidth smearing effects using the procedure outlined in Kraus et al. (2005). We estimate the loss of coherence to ≲2.5% for our longest baseline, which is within our measurement uncertainty.
Figure 2. GRAVITY (top panel; red lines, 300 s bins) and AMBER data (bottom panel; red lines, 50 s bins), overplotted with the best-fit model (blue lines). In the top right corner, we illustrate the model components (sketch not to scale).
Download figure:
Standard image High-resolution image2.2. VLTI/GRAVITY Spectro-interferometry
On 2016 June 23 we acquired observations with the VLTI/GRAVITY instrument (Eisenhauer et al. 2011, 2016, in preparation) as part of Science Verification. The observations combined the light from the four VLTI 1.8 m auxiliary telescopes and covered the K-band with R = 500. The telescopes were in the compact configuration (A0-B2-D0-C1) with projected baseline lengths 11...34 m. Seven object exposures were recorded on IRAS17216-3801, where each exposure consists of 30 interferograms taken with 10 s integration times.
The science star observations were bracketed with observations on calibrator HD 159868 (UDD = 0.358 ± 0.025 mas). Wavelength-differential visibilities and phases were extracted using the GRAVITY pipeline (Release 0.8.4; Lapeyrere et al. 2014). Over the 66 minutes covered by our data recording sequence, we see strong visibility and phase changes. Therefore, we splitted the data into bins and reduced them separately, where we use 50 s binning to minimize potential temporal smearing effects in our model fits.
2.3. VLT/NACO Imaging
We imaged the system on 2015 July 25 with the VLT/NACO adaptive optics system (Lenzen et al. 2003; Rousset et al. 2003) using a J-band filter (8 × 30 s integration) and a selection of narrowband filters centered on the [Fe ii] shock-tracing line (NB1.64; plus adjacent continuum with NB1.75 filter; 20 × 1 s integration) and on the H2 line (NB2.12; plus continuum with NB2.17 filter; 250 × 1 s integration). For the narrowband filters, we interlayed the individual pointings with calibrator observations on HD 155974 to determine the point-spread function. The system was also observed with an L'-band filter during our 2016 March 28 run (1500 × 0.2 s integration). To eliminate the thermal background, we dithered the object between the different detector quadrants. The images were reduced using the NACO pipeline (Release 4.4.1). The continuum-subtracted [Fe ii] and H2 images do not exhibit extended line emission and are therefore not further discussed.
2.4. VLT/CRIRES Spectro-astrometry
We obtained VLT/CRIRES high-resolution spectroscopy (Kaeufl et al. 2004) near the Brγ 2.166 μm line on 2012 July 5. The long-slit spectra offer spectral resolution R = 100,000. Our setup was optimized for the spectro-astrometry technique, where we measure the centroid position along the slit direction with high precision in order to derive the photocenter offset in the line with respect to the continuum. We recorded spectra for six slit orientations (PA = 0°, 60°, 120°, plus anti-parallel orientations; each with 4 × 45 s integrations), which allows us to remove spurious instrumental effects. Spectra were recorded on the calibrator Hip 88154 to correct for telluric lines. Details of our data processing procedure and on our spectro-astrometric observable (X) can be found in Kraus et al. (2012).
3. Aperture Synthesis Imaging
To derive the brightness distribution from the VLTI data in a model-independent fashion we used the IRBis image reconstruction algorithm (Hofmann et al. 2014). This algorithm minimizes a global cost function that includes a likelihood term (which determines the likelihood that the image is a representation of the data) plus a regularization term that prevents the algorithm from overfitting the data in the presence of an incomplete uv-coverage.
We first apply the algorithm to the GRAVITY data alone. The resulting image (Figure 1, bottom left) shows a binary with separation ∼58 mas (170 au at 3.1 kpc) and a K-band contrast of ∼0.6.
We were also able to reconstruct an image from the combined 2016/GRAVITY and 2012/AMBER data set, after accounting for the orbital motion between the two epochs. For this purpose, we use the method outlined in Kraus et al. (2005), where the uv-plane is rotated and scaled synchronously to the system motion. In order to quantify the separation change Δρ and position angle (PA) change Δθ between 2012 and 2016, we varied these parameters systematically on a grid and searched for the image reconstructions with the lowest χ2. We find a χ2-minimum with a counterclockwise rotation of Δθ = 8
5 between the two epochs, while the separation stayed constant within the uncertainties. This procedure results in the image shown in Figure 1 (bottom right) that represents the IRAS17216-3801 system at the 2016 epoch with a 3 × 5 mas beam.
The binary is also detected in the NACO adaptive optics images obtained with the J-band, NB1.64, NB1.75, NB2.12, and NB2.17 filters, where we see some indication that the contrast increases toward shorter wavelengths, indicating that the northern component has a redder color than the southern component (Figure 3).
Figure 3. NACO imaging of IRAS17216-3801 (top) and standard star HD155974 (bottom left). For the J-band/1.64 μm/2.12 μm image, we plot intensity contours at 30%, 60%, and 90% of peak intensity (white contours). The bottom right panel shows radial intensity cuts along the binary separation vector from IRAS17216-3801 (solid lines) and the standard star images (dashed lines). The arrow indicates the position of IRAS17216-3801-B, as derived from interferometry.
Download figure:
Standard image High-resolution imageIn all images, the northern component appears brighter than the southern component, which leads us to denote the northern component with IRAS17216-3801-"A" and the southern component with "B." Comparing the intensity profile of the binary components in the image with the interferometric beam clearly indicates that the emission is spatially extended, both in the GRAVITY and GRAVITY+AMBER image (Figure 1). The strong elongation seen in both components indicates that the circumstellar disks are seen under inclination. The extended flux around A is elongated along PA ∼ 70°, while B is elongated along PA ∼ 0°. Some artifacts are visible in the AMBER+GRAVITY images that result from the remaining imperfections in the uv-coverage. These artifacts reach up to 10% level of peak intensity, but are easy to discern, as they mainly cause a repetition of the source binary pattern along the southeast/northwest direction.
The emission from the circumstellar disks might also have been traced by VLTI/MIDI interferometry (Boley et al. 2013). This data could be fitted with an over-resolved flux component (52% of total flux) plus an elongated component with a Gaussian FWHM of 
mas oriented along 
(inclination 
; Boley et al. 2013). Within the large uncertainties, the size and orientation are consistent with the derived binary separation vector, which leads us to suggest that the compact MIDI component traced IRAS17216-3801-A + B, but did not fully resolve the components.
4. Modeling
4.1. Spectral Energy Distribution (SED)
To assist our further interpretation, we estimated the fundamental parameters of IRAS17216-3801 by compiling its SED using archival data. Integrating over the SED yields a bolometric luminosity Lbol = 6.1 × 104 L⊙ (assuming d = 3.08 kpc; Boley et al. 2013).
We then fitted the SED with the radiative transfer model grid by Robitaille et al. (2006) and find a best-fit model (ID #3009730) with an age of ∼105 years, a 24 M⊙ stellar mass, and a massive disk (gas+dust mass 1.2 × 10−2 M⊙) seen at intermediate inclination (50°). This SED model does not incorporate the companion, but we assume that it still provides a reasonable first-order estimate, as the primary will dominate the bolometric luminosity. In order to estimate the mass of the individual components, we assume that the measured K-band contrast is representative of the bolometric luminosity split between the two components and estimate from the mass–luminosity relation (L ∝ M3) the primary mass to 20 M⊙ and the secondary mass to 18 M⊙. Spatially resolved multi-wavelength observations will be necessary to improve on this very rough estimate.
4.2. Structure of Circumstellar+Circumbinary Disks (Continuum)
After identifying the basic morphology of the source with aperture synthesis imaging, we fit geometric models to the visibilities/phases, with the goal to determine the precise binary astrometry and to characterize the circumstellar emission.
Our model includes emission from two point sources (separation ρ; PA Θ), where we allow the algorithm to associate each component with spatially extended flux. We denote the contributions of the spatially extended circumstellar emission to the integrated flux of the northern and southern component with 
and 
, respectively. The flux ratio between the integrated fluxes is given by FA/FB. In addition, we allow the algorithm to attribute a fraction of the total flux (Fext/Ftot) to an over-resolved emission component.
We then fit the GRAVITY and AMBER data simultaneously (using 50 s temporal binning), where we allow for orbital motion (ΔΘ) between the two epochs (2012.378 and 2016.476). We explore the parameter space using a Levenberg–Marquardt algorithm and vary the initial guess parameters on a grid to avoid local minima.
For the brightness distribution of the circumstellar components, we assumed either elongated Gaussians ("GAUSS" model) or a temperature–power-law disk model ("DISK" model). The GAUSS model introduces three free parameters to define the geometry of the circumstellar emission for each component, namely, the inclination (iA;B), the PA (θA;B), and the FWHM (ΣA;B). The DISK model introduces four free parameters for each component, namely, disk inclination, disk PA, and the inner and outer disk truncation radius (
, 
). The radial temperature gradient of the disk is parameterized as 
, where 1500 K corresponds to the commonly assumed silicate dust sublimation temperature and r−0.43 to the theoretically derived temperature gradient for irradiated flared dust disks (Chiang & Goldreich 1997). Our DISK model implementation has already been successfully applied to fit interferometry and SED data of intermediate-mass and high-mass YSOs (Kraus et al. 2009, 2010).
The best-fit model (illustrated in Figure 2, top right) indicates that the disk around the southern component (θB = 159 ± 15°) is roughly aligned with the binary separation vector (Θ = 16676), while the northern component is oriented roughly perpendicular to it (θA = 67 ± 7°). For the radial intensity profile, we find that the DISK model provides a better representation of the data (
) than the GAUSS model (
), indicating that the disks are not smooth, but exhibit a central brightness depression, potentially tracing the opacity drop associated with dust sublimation near the star. The measured inner disk truncation radii of 2.78 and 2.49 mas (8.6 and 7.7 au) are consistent with the theoretical dust sublimation radii, if one assumes gray dust and dust sublimation temperature of 1300 K.
Comparing the properties of components A and B, we find that component A contributes a larger fraction to the total K-band flux, has a larger inner disk cavity, and a higher point-source flux contribution, all pointing to the conclusion that this is the more luminous and more massive component in the system.
Our best-fit model attributes 66% of the total flux to the circumstellar disks, 17 ± 4% and 1 ± 6% to photospheric/unresolved emission close to the stars, while 16 ± 6% extended emission is over-resolved on the shortest baselines, indicating a Gaussian HWHM size ≳50 mas. We speculate that this extended K-band flux might trace the same physical structure that is also seen in our L'-band image. The radial intensity cut through the L'-band image reveals that IRAS17216-3801-A + B contribute ≲20%, while the majority of the L' flux traces an extended component with an HWHM size ∼200 mas (Figure 3, bottom right). We note that this measured size is in good agreement with the expected inner truncation radius of the circumbinary disk around IRAS17216-3801-A + B. Hydrodynamic simulation suggests that binary systems with low eccentricity truncate their circumbinary disk at ∼(2–3)× the semimajor axis a. Using our lower limit on the semimajor axis (a ≳ 58 mas), we therefore expect the circumbinary disk to extend beyond ≳120–180 mas, in good agreement with the size of the L'-band structure.
4.3. Accretion-tracing and Disk-tracing Emission Lines
The GRAVITY wavelength-different visibilities/phases allow us to determine the distribution of the Brγ and CO line-emitting gas. It is interesting to note that the visibility/phase signals in Brγ and CO appear to be uncorrelated and are stronger in CO (despite their weaker line strength), which already indicates that these lines emerge from substantially different regions in the circumstellar environment.
The Brγ-line traces hot (≳104 K) ionized gas and is a well-established mass accretion/ejection tracer (Natta et al. 2004; Kraus et al. 2012; Caratti o Garatti et al. 2016). Therefore, we model Brγ as unresolved emission at the position of the stars. Our modeling allows us to exclude scenarios that associate the Brγ only with the primary or secondary component, as this results in much stronger visibility/phase signals than measured. We achieve a satisfactory fit by attributing 40 ± 10% of the Brγ-emission to the primary and 60 ± 10% to the secondary (Figure 4, right).
Figure 4. CRIRES spectro-astrometric (left) and GRAVITY spectro-interferometric data (right) on IRAS17216-3801. Top left: the Brγ line is single-peaked and shows strong astrometric signals (X) toward different PAs. Bottom left: continuum-corrected two-dimensional photocenter offsets derived from the astrometric signals, with the color indicating the wavelength of the corresponding channel. For clarity, we also indicated the position of IRAS17216-3801-A + B, as derived with GRAVITY+AMBER (Section 4.2). Right: spectra, differential phases, and differential visibilities derived from our GRAVITY data. The black lines show the visibilities/phases that correspond to the Brγ/CO model outlined in Section 4.3.
Download figure:
Standard image High-resolution image
The CRIRES spectro-astrometry provides an independent method for constraining the origin of the Brγ emission. The Brγ photocenter is displaced along PA = 161 ± 3°, which is consistent with the binary separation vector derived from the AMBER data for the same epoch (1594). The offset is ∼14 ± 1 mas, which places the photocenter between the two stars (Figure 4, left), again indicating that both components contribute to the line emission. Modeling the photocenter displacement quantitatively we find that the primary (secondary) contributes 38% (62%) of the Brγ-flux. Therefore, the GRAVITY and CRIRES data indicate independently that the (likely lower-mass) B component accretes at a higher rate.
The CO bandhead emission is believed to trace warm (≳103 K) neutral gas in the more extended disk regions (few au to a few hundred au in massive YSOs; Ilee et al. 2013). To model the GRAVITY CO line data, we approximate the CO-emitting region with a Gaussian of FWHM ΣCO. We vary the position of the Gaussian along the binary separation vector, where sCO denotes the separation of the Gaussian to the primary. We find that the CO emission is extended (ΣCO = 5.8 ± 0.8 mas) and that it is centered on a position between the two stars (sCO = 31.6 ± 3.3 mas), possibly including contributions from both circumstellar disks and the gas streams between the disks (Figure 4, right).
5. Conclusions
Our VLTI interferometry resolves the high-mass protostar IRAS17216-3801 into a close (58 mas = 170 au) binary. The system is ∼3× as massive (∼20 + 18 M⊙) and ∼5× more compact than the benchmark protobinary IRAS20126 + 4104 (separation 850 au; Sridharan et al. 2005). IRAS17216-3801 is also ∼3× more compact than the NGC7538 IRS1a/b system (separation 500 au), where Keplerian-rotating methanol maser disks (Pestalozzi et al. 2004; Moscadelli & Goddi 2014) and ammonia line absorption consistent with a circumbinary envelope were found (Goddi et al. 2015).
Our IRAS17216-3801 imaging traces the thermal dust emission in the circumprimary and circumsecondary disk and reveals that both disks are misaligned with respect to the binary separation vector. Various formation mechanisms have been proposed that might have resulted in the observed disk misalignments, including turbulent disk fragmentation, perturbation by a third component, star–disk capture, or infall of material whose angular momentum vector was misaligned to that of the gas that formed the binary initially (Bate et al. 2010). Once a highly misaligned system has formed, tidal interactions will work toward re-aligning the disks, which has been proposed to happen on the viscous timescale (Papaloizou & Terquem 1995) or on the much shorter precession timescale (Bate et al. 2000). From the deduced stellar masses and outer disk radii (Table 1) we estimate the precession timescale of the IRAS17216-3801 circumprimary disk to ∼200,000 years and ∼900,000 years for the circumsecondary disk, although these values should be considered upper limits as we derive only lower limits on the disk radii from our infrared observations. Therefore, our observations indicate that the tidal realignment is still ongoing, highlighting the young dynamical age of the system.
Table 1. Model-fitting Results (Section 4)
| GAUSS model | DISK model | |||
|---|---|---|---|---|
| Binary separation | ρ | (mas) | 57.94 ± 0.24 | 57.93 ± 0.16 |
| Binary PA | Θ | (°) | 166.8 ± 0.2 | 166.76 ± 0.2 |
| PA change 2012–2016 | ΔΘ | (°) | 7 ± 2 | 7 ± 2 |
| Flux ratio | FA/FB |
|
|
|
| Extended flux contribution | Fext/Ftot | 0.16 ± 0.09 | 0.16 ± 0.06 | |
| Circumprimary Disk, Continuum (Northern Component, A) | ||||
| Contribution extended flux |
|
0.60 ± 0.08 | 0.64 ± 0.04 | |
| Disk inclination | iA | (°) | 89 ± 10 | 60 ± 10 |
| Disk PA | θA | (°) | 64 ± 9 | 67 ± 7 |
| Gaussian FWHM size | ΣA | (mas) | 7.63 ± 0.8 | ... |
| Inner disk radius |
|
(mas) | ... | 2.77 ± 0.39 |
| Outer disk radius |
|
(mas) | ... | >12 |
| Circumsecondary Disk, Continuum (Southern Component, B) | ||||
| Contribution extended flux |
|
0.90 ± 0.05 | 0.96 ± 0.09 | |
| Disk inclination | iB | (°) | 42 ± 7 | 38 ± 10 |
| Disk PA | θB | (°) | 185 ± 32 | 159 ± 15 |
| Gaussian FWHM size | ΣB | (mas) | 4.60 ± 0.4 | ... |
| Inner disk radius |
|
(mas) | ... | 2.49 ± 0.42 |
| Outer disk radius |
|
(mas) | ... | >10 |
|
1.71 | 1.11 | ||
Note. We define inclination 0° as face-on orientation. PAs are measured east of north.
Download table as: ASCIITypeset image
The stronger orbital misalignment of the circumprimary disk might also explain its larger spatial extend, as the misalignment results in a weaker Lindblad torque that acts to truncate the circumstellar disks (Lubow et al. 2015).
The circumstellar disks are fed by a circumbinary disk/envelope that has been resolved by our VLT/NACO L'-band imaging. We detect an extended structure, whose size matches roughly the expected dynamical truncation radius of the IRAS17216-3801 binary system. Both circumstellar disks show signatures of ongoing accretion, where we measure a higher Brγ-line luminosity at the position of the secondary than at the primary. This suggests that the secondary disrupts the accretion stream on the primary and channels most of the infalling material onto the circumsecondary disk, confirming the prediction of hydrodynamic simulations (Bate & Bonnell 1997). We resolved for the first time the CO line-emitting region in a massive YSO and find that it traces primarily warm neutral gas located between the two components.
With its unique properties, IRAS17216-3801 provides an ideal laboratory for studying the formation of young high-mass multiple systems and to unravel how such systems accrete from their disks.
We thank the GRAVITY consortium and the Science Verification team, which is composed of ESO employees and GRAVITY consortium members (https://www.eso.org/sci/activities/vltsv/gravitysv.html). We acknowledge support from an STFC Rutherford fellowship/grant (ST/J004030/1, ST/K003445/1), Marie Sklodowska-Curie CIG grant (#618910), Philip Leverhulme prize (PLP-2013-110), and ERC Starting grant (grant agreement #639889).
Facilities: VLTI, VLT.
Footnotes
- *
Based on observations made with ESO telescopes at Paranal Observatory under program IDs 60.A-9174(A), 089.C-0819(A,C), 089.C-0959(D,E), 094.C-0153(A), 096.C-0652(A).