{'fraction_non_alphanumeric': 0.030752959362160843, 'fraction_numerical': 0.01948500996623683, 'mean_word_length': 4.650080441277868, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
No evidence for protoplanetary disk destruction by OB stars in the MYStIX sample

HST images of proplyds in the Orion Nebula, as well as submillimeter/radio measurements, show that the dominant O7 star Theta1 Ori C photoevaporates nearby disks around pre-main sequence stars. Theory predicts that massive stars photoevaporate disks within distances of order 0.1 pc. These findings suggest that young, OB-dominated massive H II regions are inhospitable to the survival of protoplanetary disks, and subsequently to the formation and evolution of planets. In the current work, we test this hypothesis using large samples of pre-main sequence stars in 20 massive star-forming regions selected with X-ray and infrared photometry in the MYStIX survey. Complete disk destruction would lead to a deficit of cluster members with excess in JHKs and Spitzer/IRAC bands in the vicinity of O stars. In four MYStIX regions containing O stars and a sufficient surface density of disk-bearing sources to reliably test for spatial avoidance, we find no evidence for the depletion of inner disks around pre-main sequence stars in the vicinity of O-type stars, even very luminous O2-O5 stars. These results suggest that massive star-forming regions are not very hostile to the survival of protoplanetary disks and, presumably, to the formation of planets.


INTRODUCTION
Proplyds are protoplanetary disks that are heated and thereby illuminated by ultraviolet light from nearby massive stars. The first proplyds were discovered in the Orion Nebula's Trapezium Cluster around θ 1 Ori C (spectral type O7; Simón-Díaz et al. 2006) using radio observations (Churchwell et al. 1987). The most well-known observations of the Orion proplyds come from optical imaging by the Hubble Space Telescope (O'Dell et al. 1993;O'Dell & Wen 1994;Bally et al. 1998;Smith et al. 2005;Ricci et al. 2008). Many of the images show cometary tails directed away from θ 1 Ori C, which result from the combined action of far ultraviolet (FUV) radiation, which heats and therefore expands the disk, and the stellar wind, which creates a bow shock at the winddisk interface (García-Arredondo et al. 2001). Spectroscopic observations confirm that disks lose mass to outflows of heated gas (Henney & O'Dell 1999), which presumably also carry along small dust particles. Indeed, sub-mm measurements of dust emission from disks in the Trapezium Cluster reveal disk mass truncation within tenths of a parsec of θ 1 Ori C (Mann & Williams 2009b, 2010Mann et al. 2014), though in the Flame Nebula (NGC 2024), which is somewhat younger than the Trapezium Cluster , Mann et al. (2015) find no evidence of such an effect around the less massive IRS 2b (spectral type O8-B2; Bik et al. 2003).
These observations have motivated models of mass loss rates for externally photoevaporated disks as a function of incident ultraviolet flux and the properties of the disk, which can in turn be compared with observed disk life-times and initial masses Störzer & Hollenbach 1999;Adams et al. 2004Adams et al. , 2006Clarke 2007;Fatuzzo & Adams 2008). Predicted mass loss rates vary significantly among these studies; the shortest estimates of disk lifetimes are of order 10 5 yr (e.g., Störzer & Hollenbach 1999), which agree with the observational findings of Henney & O'Dell (1999), and are much shorter than ordinary viscous accretion timescales of 10 6 -10 7 yr (Haisch et al. 2001;Hernández et al. 2007, and references therein). By contrast, the combined external photoevaporation and viscous accretion model of Anderson et al. (2013) predicts that external photoevaporation shortens disk lifetimes by a factor of only a few, implying a minimal effect on planet formation and evolution. Observationally constraining planet formation and photoevaporative destruction timescales is therefore vital to determining whether OB-dominated star-forming regions are hostile to planet formation.
While the aforementioned theoretical works disagree somewhat as to how far external photoevaporation by massive stars proceeds into the disk (estimates are of order tens of AU from the disk host star), the destruction of the outer disk may lead indirectly to the depletion of the inner disk Anderson et al. 2013), as material accreted onto the host star is no longer replenished by the now-missing outer disk, and accretion in the inner disk occurs on timescales of 10 4 -10 5 yr (Hartmann et al. 1998). Hence, the total depletion of the outer disk may imply the quick disappearance of the inner disk. By determining whether the inner regions of disks are present in the vicinity of O stars in regions of different ages, the rate at which the outer regions are destroyed can be observationally constrained. This is the goal of the present work.
A valuable new dataset that examines disk properties in a range of UV-dominated environments is the Massive Young Star-Forming Complex Study in Infrared and X-ray (MYStIX), a survey of 20 young, OB-dominated regions that combines X-ray and infrared photometry to study clustered star formation, early cluster evolution, and protoplanetary disk evolution . The advantage of the MYStIX approach is that it captures both X-ray selected disk-free and infrared-selected disk-bearing stars. The catalog of MYStIX Probable Complex Members (MPCMs) consists of 31,784 young stellar objects (YSOs; Broos et al. 2013). Disk-bearing YSOs are identified by infrared excess in the 1-8µm spectral energy distribution (SED), while disk-free YSOs are identified through X-ray emission and an SED in the 1-8µm range that is consistent with a bare stellar photosphere. Membership criteria and YSO classification are discussed in Sections 2.1 and 2.2, respectively. The catalog also includes OB stars already identified in the literature (discussed further in § 2.3).
The MPCM catalog, with thousands of lower-mass YSOs and dozens of OB stars, therefore in principle contains the information needed to determine whether external photoevaporation by massive stars leads to the complete depletion of the inner disk. Celestial coordinates from the MPCM catalog, along with known region distances (ranging 0.4-3.6 kpc), can be used to calculate projected distances between disk-bearing YSOs and O stars in each MYStIX region. With this information, we are able to test for any tendency of disk-bearing YSOs to spatially avoid O stars, which would imply complete photoevaporative destruction of the disks. It should be emphasized that a failure to find evidence that O stars are spatially avoided by disk-bearing YSOs would not imply that disks are not externally photoevaporated; rather, it would mean that the process of external photoevaporation is not as rapid or complete as predicted by some theory.
After examining all 20 MYStIX regions, we find no evidence of disk destruction 5 , with the results for two regions (Carina Nebula and M 17) being inconclusive due to observational difficulties. It seems that the ultraviolet radiation fields found in OB-dominated regions do not present such a hostile environment for protoplanetary disks as suggested in some previous works. In Section 2, we discuss membership selection, YSO classification, and the exclusion of high-nebulosity regions in each MYStIX field. In Section 3, we present the results of our tests for spatial avoidance of O stars by disk-bearing YSOs, and compare these results with several previous works for two regions. Interpretation relating to disk destruction astrophysics and implications for planet formation are discussed in Section 4.

SFR membership
The MYStIX catalog is derived by cross-matching photometric point sources from X-ray, near-infrared, and mid-infrared bands. X-ray archival data are taken from the ACIS instrument on the Chandra X-ray Observatory; near-IR (JHK S ) data for each region come mostly from either 2MASS or United Kingdom Infrared Tele-scope (UKIRT) archival data; and mid-infrared photometry is taken from the Spitzer Space Telescope's IRAC instrument (3.6 µm, 4.5 µm, 5.8 µm, and 8.0 µm). Further details for these data and their analysis are given by Feigelson et al. (2013) and the associated MYStIX observational papers (Kuhn et al. 2013a,b;Naylor et al. 2013;Povich et al. 2013;Broos et al. 2013;Townsley et al. 2014). Broos et al. (2013) assign a probability of region membership to each X-ray source using a naive Bayes classifier, wherein each attribute is assumed to contribute independently to the probability of membership. Attributes used for this classification are J -band magnitude; X-ray median energy, variability, and local source density; 4.5 µm magnitude; infrared excess (with respect to bare stellar photosphere); and mid-IR local source density. Each X-ray source is subsequently classified as a foreground star, an MPCM, a background (Galactic) star, or a galaxy/active galactic nucleus, using training sets constructed according to the methods presented by Getman et al. (2011) and Broos et al. (2013).
Infrared and X-ray sensitivities vary between MYStIX regions due to different telescope exposures, region distances, and levels of absorption. X-ray sensitivity also varies within each MYStIX field. It is therefore difficult to reliably compare the surface densities of YSOs of each class among different MYStIX regions. However, this does not present a significant problem for studying YSO surface densities in the vicinity of O stars across the relatively small distances (tenths of a parsec) involved in disk destruction.
In each region, completeness with J -band magnitude varies, and certainly many low-mass members are missing from the MYStIX sample. This limits the scope of our conclusions, as disks around lower-mass stars are predicted to experience higher mass-loss rates from external photoevaporation due to lower escape speed for outflows of heated gas (see, e.g., Anderson et al. 2013). A failure to identify disk destruction for the stars in the MYStIX sample therefore does not preclude the possibility that the disks around stars with masses below our completeness limits are being completely destroyed. Povich et al. (2013) classify MPCMs by YSO class using both color-color cuts (to exclude contaminants and bad photometry from diffuse nebular emission) and fitting of model spectral energy distributions (SEDs) based on near-and mid-infrared photometry for each source. The stellar atmosphere models of Castelli & Kurucz (2004) and the YSO SED models of Robitaille et al. (2006) are used. Observed near-/mid-IR SEDs most consistent with bare stellar photospheres are classified as disk-free (Class III YSO), while those consistent with the model SEDs of Class I or II YSOs are classified as diskbearing. Disk-free YSOs are additionally required to be detected as X-ray sources in order to avoid contamination by the often very populous field star population.

YSO classification
For 18 of 20 MYStIX regions, in order for a given source to be classified, it must be detected at 3.6 µm and 4.5 µm, in addition to at least two of the remaining five bands: JHK S , 5.8 µm, and 8.0 µm. The Carina Nebula was analyzed in a separate work (Povich et al. 2011), where photometry for any four of the seven infrared bands were required. For the Orion Nebula, we use the catalog of Megeath et al. (2012), who classified YSOs based on midinfrared colors and spatial distribution, and, as in the case of Carina, require that disk-free YSOs be detected in 4 of 7 infrared bands.
Although the MYStIX sample includes both diskbearing and disk-free stars, it should be noted that the science goals of this paper can in principle be met using only data for disk-bearing YSOs (along with the O stars). We test for disk destruction by looking for any apparent dropoff in the surface density of disks found at small distances from O stars. We do not calculate disk fraction as a function of distance from O stars, because the estimation of ratios of Poisson-distributed variables can be very uncertain (Brown et al. 2001;Park et al. 2006). We also include a parallel analysis of disk-free YSOs as a test of reliability for any apparent avoidance effect for the disks. The disk destruction interpretation of an apparent tendency for disk-bearing YSOs to avoid O stars would be undermined by any tendency for disk-free YSOs to do the same, which would imply observational rather than astrophysical origins for the avoidance effect.

Selection of OB stars
The MYStIX catalog incorporates published OB stars based on optical spectroscopy. Spectral types are taken from Skiff (2009) and SIMBAD, and range from B3 to O2. Positions from the literature are matched to JHK S photometry as described by Broos et al. (2013).

Excluding PAH-contaminated regions
Many MYStIX regions contain extensive nebulosity where ultraviolet radiation and winds from OB stars heat dust on the periphery of the surrounding molecular cloud, leading to thermal emission from polycyclic aromatic hydrocarbons (PAHs) at mid-infrared wavelengths (particularly in the IRAC 5.8 µm and 8.0 µm bands). In principle, mid-IR source contamination at these wavelengths due to bright PAH emission could increase the apparent disk fraction in a given area on the sky, since disk-free sources will be more difficult to detect in the IR bands. A disk-free source in a PAH-dense region is also more likely to be misclassified as disk-bearing than a disk-bearing one is to be misclassified as disk-free, which again will lead to an overestimate of disk fraction. Depending on the position of such a region with respect to O stars, a disk destruction effect could potentially be falsely created or falsely obscured.
To mitigate the problem of PAH contamination, we exclude the most strongly contaminated areas of each of the 20 MYStIX regions based on 8.0 µm Spitzer /IRAC images, which also excludes those regions where infrared images are strongly contaminated by the PSFs of very bright stars, such as η Carinae (Carina Nebula). In each region, we calculate the median 8.0 µm flux around each YSO within a ∼7"×7" square. We construct histograms of these median local background fluxes (MLBFs) for disk-bearing and disk-free YSOs in each MYStIX region. Figure 1 shows this analysis for a typical case, DR 21, where the disk-bearing and disk-free distribution follow each other up to an MLBF of ∼ 500 MJy sr −1 . At higher MLBFs, more disks are found due to their higher mid-infrared fluxes, as well as possible misclassifications due to unreliable photometry in PAH-dense regions. An MLBF of 500 MJy sr −1 threshold for the DR 21 field is also validated by visual examination of the 5.6 and 8.0 µm images; fainter point sources are more difficult to identify at this and brighter nebulosity levels. In each region, at the MLBF where the two distributions seem to diverge in this way, we select a cutoff value. We then calculate contours based on this cutoff value in the 8.0 µm Spitzer /IRAC image for each region, and exclude all sources (including O stars) from our analysis. Possible complications relating to PAH regions and their exclusion are discussed throughout Sections 3 and 4.
2.5. Summary of membership For each MYStIX region, Table 1 shows distances, selected MLBF cutoff levels, and the numbers of O stars (N O ), disk-bearing YSOs (N D ), and disk-free YSOs (N ND ) in the MPCM sample after MLBF cutoffs are applied. For disk-bearing YSOs in each region, we also compute the approximate 50% mass completeness limits, M 50 , in other words, 50% of stars of mass M 50 are detected in our sample, based on a fit of the stellar initial mass function of the form given by Maschberger (2013). Regions are divided into three categories, described in Section 3. Information on subcluster ages and structures can be found in Getman et al. (2014) and Kuhn et al. (2014), respectively.

RESULTS
In order to test for disk destruction in a given region, we calculate projected distances between O stars and their disk-bearing and disk-free neighbors. For each region under study, we combine these data separately for early-and late-type O stars and construct empirical cumulative distribution functions (ECDFs) for distances within 0.5 pc, a limiting distance which is well beyond the radii of disk destruction predicted by theory 6 . Figure 2 illustrates the construction of these ECDFs using mock data. Disk destruction would be revealed by a sudden dearth of disk-bearing YSOs at short distances from O stars. The use of ECDFs avoids arbitrary radial binning, and the separate analysis of disk-bearing and disk-free YSOs avoids uncertain small-sample disk fractions.
The observed geometries of the regions under study are complex due to limited X-ray fields of view, PAHcontaminated regions and the exclusion thereof, and the intrinsic concentration of YSOs into distinct subclusters with different ages Getman et al. 2014). Therefore, it cannot be assumed a priori that the number of YSOs around a given O star will increase with distance in a predictable manner. In practice however, this distribution is close to parabolic within distances of ∼0.5 pc, as seen in the figures in Sections 3.2 and 3.3. As such, for each ECDF result, we also plot a uniform distribution based on two-dimensional geometry, i.e., a parabola, while also taking into account the effect of PAH exclusion regions and the edge of the X-ray field of view. These parabolas provide a test of statistical significance for any apparent avoidance effect. The presence of fewer disk-bearing YSOs than disk-free within some distance of O stars is not regarded as evidence for disk destruction so long as the ECDF for the disk-bearing YSOs follows its predicted parabolic increase. This is because the distance between an O star and the nearest disk-bearing YSO is expected to be large in a region with an intrinsically sparse population of disks.
We organize the results for each of the 20 MYStIX regions into one of three categories: those without O stars ( § 3.1), those with O stars but with apparent surface densities of disk-bearing YSOs that are too low to test for avoidance of O stars ( § 3.2), and those with O stars and sufficiently dense disk-bearing YSOs to test for disk destruction ( § 3.3). The term "sufficiently dense" will be explained in Section 3.2. Appendix A contains figures showing, for each MYStIX region containing O stars, the spatial distribution of O stars, disk-bearing and disk-free pre-main sequence stars, and excluded areas with nebular contamination above the MLBF cutoff superposed on the Spitzer /IRAC 8.0 µm image.

Regions without O stars
Seven of 20 MYStIX regions lack O stars in the sample used in this work. DR 21 has no O stars to begin with, with the only OB-type star in the MYStIX catalog being a B0.5V star. The other six regions each have one or two O stars, but all were excluded from our analysis due to the PAH exclusion scheme described in Section 2.4. These are the Flame Nebula, NGC 2362, RCW 36, RCW 38, the Trifid Nebula, and W 40. Except for RCW 38, these regions contain only one or two late O stars in the non-PAH excluded sample, and are therefore unlikely to produce a disk destruction effect. In some cases, the surrounding density of disk-bearing YSOs (to the extent that they can be observed on account of PAH contamination) are so low that disk destruction would not occur anyway, as in the regions described below in Section 3.2. RCW 38 contains an O5.5V star, however the apparent surface density of YSOs around it is far too low for significant disk destruction to be observed.

Regions with low apparent densities of disk-bearing
YSOs In order to determine whether disks spatially avoid O stars within some radius due to complete photoevaporative destruction, the surface density of disk-bearing YSOs must be great enough to detect a deficit of them at short distances to O stars. Astrophysically, disk destruction cannot occur if there are no disks present in the vicinity of a given O star, even if the UV radiation field within some radius of the star is in principle high enough to photoevaporate a disk. Observationally, disk destruction cannot be detected if too few disk-bearing YSOs are present to reliably measure the distribution of their distances from O stars and detect a dropoff at short distances.
In 7 out of the 13 MYStIX regions containing O stars, the observed surface densities of disk-bearing YSOs around O stars are too small to detect destruction of nearby (d 0.2 pc) disks: the Lagoon Nebula, NGC 1893, NGC 2264, NGC 3576, NGC 6334, W 3, and W 4. These regions therefore cannot be used to constrain the physics of disk destruction. ECDF results for these seven regions are shown in Appendix B. seven regions discussed in this section. We determine the local intrinsic surface densities (LISDs) for the eight O stars in those four regions, and compare them with the LISDs of the O stars in regions with higher apparent surface densities ( § 3.3). We confirm that LISDs around O stars in the four regions discussed in this section are systematically lower, suggesting that intrinsically low surface densities, convolved with relatively high stellar mass sensitivity limits of ≥ 0.6 M , are responsible for the sparsely populated ECDFs seen in Figures B1-B7.
It is interesting to note that in the case of NGC 2264, there is a relatively large number of disk-bearing YSOs located within 0.5 pc of the late O stars, however a strong disk fraction gradient in the region leads to significantly non-parabolic distributions for both the diskbearing and disk-free YSOs due to the global age gradient in NGC 2264 identified by ). This prevents the identification of disk destruction at distances less than 0.2 pc.
3.3. Regions tested for disk destruction Six of 20 MYStIX regions have O stars and a sufficient density of disk-bearing YSOs to test for disk destruction: M 17, NGC 6357, and the Carina, Eagle, Orion, and Rosette Nebulae. In all six cases, our ECDF analysis reveals no evidence that disk-bearing YSOs spatially avoid O stars in a manner consistent with photoevaporative dissipation. In the following, we discuss each of these regions. For the Eagle and Rosette Nebulae, we also compare our results with previous similar works examining disk fraction around massive stars based on infrared excesses (Guarcello et al. 2007(Guarcello et al. , 2009Balog et al. 2007).

Carina Nebula
ECDF results for the Carina Nebula are shown in Figure 3, with the caveat that we have excluded O stars in Trumpler 14, a subregion of Carina which suffers from significant crowding, resulting in unreliable infrared photometry. Figure A1 shows Trumpler 14 to the northwest next to an illuminated molecular cloud above the MLBF cutoff.
Though the disk-bearing YSOs of Carina appear to avoid early O stars at distances of ∼0.1 pc, the diskfree YSOs also appear underabundant at distances up to ∼0.15 pc. Performing this analysis without the requirement of detection in four infrared bands for the disk-free objects (i.e., designating all non-infrared excess sources as disk-free), we find that this apparent avoidance effect for the disk-free objects disappears. This strongly suggests that the origins of this effect are observational rather than astrophysical.
The Carina Nebula is 2.3 kpc away; at this distance, a linear separation of 0.1 pc between two objects corresponds with an angular separation of just under 10 arcsec, which is comparable to the observed point spread functions (PSFs) of early O stars in several infrared bands. This problem of crowding may be further exacerbated by local illumination of PAHs below the MLBF cutoff, extending the effective PSF of each O star in the mid-infrared.
Our results for the Carina Nebula are therefore inconclusive, and would best be studied through multiwavelength observations at higher resolutions to overcome the problem of crowding. For now, the hypothesis of disk destruction within a few tenths of a parsec of early O stars in the Carina Nebula cannot be conclusively supported or refuted. Figure 4 shows results for M 17. As in the Carina Nebula, disk-bearing YSOs seem to spatially avoid early O stars within ∼0.15 pc, however the disk-free YSOs appear to be underabundant at these distances as well. As with the aforementioned case of Carina, an astrophysical explanation is unlikely. Instead, a combination of crowding, large PSFs of early O stars, and PAH contamination result in a decreased detection of all infrared stars in the Spitzer images in the dense core of M 17. Less than 25% of MYStIX X-ray sources in M 17 have reliable infrared photometry in enough bands for SED classification, suggesting that crowding plays an important role. As with Carina, disk destruction is neither supported nor refuted by the available data.   Figure 5 shows results for the Eagle Nebula. For both early and late O stars, disk-bearing and disk-free YSOs follow the parabolic distributions expected for twodimensional geometry, and the disk-bearing YSOs show no sign of O star avoidance. This result conflicts with the findings of Guarcello et al. (2007) and Guarcello et al. (2009), who report a decrease in disk fraction for the YSOs most strongly illuminated by O and early B stars in NGC 6611, the central dense region of the Eagle Nebula. In the following paragraphs, we identify the sources of this discrepancy. Guarcello et al. (2007) calculate total incident bolometric flux from massive stars (spectral type earlier than B5) on each disk-bearing and disk-free YSO in NGC 6611, based on the OB spectral types of Hillenbrand et al. (1993) 7 . The YSO membership catalog was based ini-tially on BVI, 2MASS, and Chandra/ACIS data, with disks identified from color-color diagrams (as opposed to the model SED fitting used in the current work). The catalog and YSO classifications were updated by Guarcello et al. (2009, hereafter G09) using Spitzer /IRAC data. In both papers, the authors find a significant decrease in disk fraction of the most illuminated YSOs, suggesting that some disks have been photoevaporated by massive stars.

Eagle Nebula
In order to investigate the apparent disagreement between the results of G09 and the current work, we apply our ECDF-based analysis to the data of G09, and apply the bolometric flux-based analysis of G09 to MYS-tIX data within NGC 6611. Incident bolometric fluxes from massive stars in the same range of spectral type (O-B4.5) were calculated using two-dimensional projected distances. Bolometric luminosities were determined from spectral type; for O stars, we use the synthetic photometry of Martins & Plez (2006), and for B stars, we use the values given by Carroll & Ostlie (2006, Appendix G). Figure 6 shows the results of these analyses. The upper three panels show the analysis of MYStIX data using dis- tance ECDFs for early/late O stars and calculating disk fraction for bins of incident bolometric flux as in G09. The middle row shows the same set of methods applied to the data of G09. We also consider the possibility of differing mass completeness limits between disk-bearing and disk-free YSOs within the MYStIX and G09 samples. We find that the J -band luminosity functions of the disk-bearing and disk-free members are similar for the G09 sample, whereas they differ significantly in the MYStIX sample. We therefore explore the effect of truncating the MYStIX sample by J -band magnitude, imposing a cut of 10<J <16. Results are shown in the lower three panels.
The left three panels of Figure 6 show that the difference between our results and those of G09 is primarily due to membership. The MYStIX sample within NGC 6611 contains many more disk-bearing members than the G09 sample (715 in MYStIX sample versus 264 in G09 sample). This is likely due to the fact that for near-IR photometry, MYStIX uses the UKIDSS survey, while G09 uses 2MASS, whose limiting JHK S magnitudes are not as deep.
Of the six panels of Figure 6 showing ECDF-based analysis, only one shows an apparent destruction effect, which is for the two early O stars using the G09 YSO sample. Figure 7 is an RA-Dec map of NGC 6611, showing the O stars along with disk-bearing and diskfree YSOs of the G09 sample. Both of the early O stars, which dominate the high energy flux in the region, are located toward the edge of the cluster, where especially for the region to the northwest (upper right), the disk fraction appears significantly lower even in areas without O stars. In the case of the O5V star toward the southwest (lower right), MYStIX identifies at least half a dozen disk-bearing YSOs not seen in the G09 sample, which eliminates the avoidance effect.
The right three panels of Figure 6 show disk fraction as a binned function of incident bolometric flux from OB stars (hereafter simply "flux"). We have used the same range in flux used by G09 (20-900 erg s −1 cm −2 ). Interestingly, only the MYStIX data show a monotonic decrease in disk fraction with flux similar to that seen in the Figure 13 of G09. The difference between our flux analysis of the G09 data (right column, middle row) and Figure 13 of G09 likely stems from the use of different bolometric luminosities for OB stars. G09 does not discuss the calculation of these quantities. Nonetheless, the flux values fall within a similar range, and in the following paragraph we discuss the sensitivity of the binned flux distributions to arbitrary choices made in the G09 analysis, which may further explain this discrepancy.
In order to assess the reliability of the binned fluxbased analysis of G09, we test a range of values for the lower and upper bin edges (in G09 and Figure 6, these are 20 and 900 erg s −1 cm −2 , respectively). While the shape of the distribution is not sensitive to the choice of upper limit, we find that it is strongly sensitive to the choice of lower limit. Depending on the value chosen, an apparent destruction effect can be artificially produced or made to disappear completely with the same data. If a lowered disk fraction in the highest flux bin is to be interpreted as the result of disk destruction, then the overall shape of the distribution should not depend on the arbitrary choice of lower flux cutoff.
In our reanalysis of the G09 data set using the G09 flux histograms, if no lower flux cutoff is imposed, then a steep dropoff in disk fraction with incident flux is observed, as in the original G09 result. For the MYStIX data set, Figure 6 shows a steep drop-off in the incident flux distribution, however this apparent destruction effect disappears completely if a lower flux cutoff of 50 erg s −1 cm −2 is applied. Of the G09 point sources with incident fluxes less than 50 erg s −1 cm −2 , about a third of them are located at projected distances of more than a parsec from the nearest B star, and are typically even farther from the nearest O star. It would therefore seem that high disk fractions on the periphery of the G09 field of view, which are apparent in Figure 7, dominate the shape of the distribution even in the high incident flux bin. The disk fraction in these outlying regions is not relevant to disk destruction, and the size of the field of view under study, which is arbitrary in an astrophysical sense, should not dictate the outcome of a test for disk destruction among YSOs located near massive stars.
We conclude that the apparent disk destruction effect reported by G09 is not astrophysical, and instead results from sample incompleteness and the choice of binning scheme used in the G09 analysis, which is strongly sensi-   There are 264 disk-bearing YSOs (red circles), 759 disk-free YSOs (blue circles), and 10 O stars (large green circles), labeled by spectral types from Hillenbrand et al. (1993). The two earliest O stars (O5V and O5.5V) lie toward the edge of the cluster where the disk fraction is lower, suggesting that the destruction effect reported by G09 is instead the result of disk fraction gradients in the cluster, as well as incompleteness of the sample of disk-bearing YSOs.
tive to disk fractions at very large distances from massive stars. It remains unclear whether the analysis of incident bolometric fluxes by G09 is in principle capable of detecting disk destruction. We believe that the non-detection of disk destruction in NGC 6611 in the current work is more reliable, as the MYStIX sample of disk-bearing YSOs is significantly more complete (including in the vicinity of the earliest O star in the region), and our ECDF-based analysis avoids arbitrary binning schemes and is less sensitive to the surface densities of YSOs at large distances.

Rosette Nebula
Results of our ECDF analysis for the Rosette Nebula are shown in Figure 8. For both early and late O stars, the distribution of distances to disk-bearing YSOs very closely follows the distribution predicted from simple two-dimensional geometry. Given the height of the steps of the disk-bearing ECDF in both the left and right panels (early and late O stars, respectively), our data cannot probe the physics of disk destruction over distances less than ∼0.15 pc. Hence, constraints on disk destruction models derived from these data are weaker than for the other regions discussed in this subsection, which all show higher densities of disk-bearing YSOs around O stars. In the remainder of this section, we compare these results with those reported previously in the literature. Balog et al. (2007, hereafter B07) test for disk destruction in NGC 2244 by calculating the distance of each disk-bearing and disk-free YSO to the nearest O star (of which there are 5 in the region), and compute histograms of these nearest distances. They report a decrease in disk fraction for YSOs located within 0.5 pc of an O star.
However, the bin widths used in the B07 histograms are much too large to be sensitive to disk destruction on the scales predicted by theory (the smallest bin extends to 0.5 pc). Given that the distances under consideration by B07 are comparable to the size of the cluster, other possible sources of spatially varying disk fraction must be considered, such as the large cluster age gradients reported by Getman et al. (2014). Figure 9 shows a map of NGC 2244 using the point source data of B07. On visual inspection it is clear that disk fraction decreases significantly with distance from the cluster center, even in regions far away from O stars.
To further explore this problem, in Figure 10, we apply the ECDF analysis used in the current work to the data of B07. This time we extend our analysis to larger distances due to the large (∼0.5 pc) radius of effect for disk destruction claimed by B07. We find that the distributions of disk-bearing YSOs in NGC 2244 increase smoothly, showing no avoidance effect. The ECDF of disk-bearing YSOs for late O stars (right panel) is somewhat irregularly shaped, including at ∼1 pc and ∼1.6 pc. This is most likely explained by intrinsic disk fraction gradients and missing sources due to PAH contamination rather than disk destruction.
In the current work, we avoid the potential problems of binning by using ECDFs, and we restrict our analysis to distances within 0.5 pc of O stars in order to ignore large-scale gradients in the distribution of disks across the region unrelated to disk destruction. If disk destruction acts gradually over ∼0.5 pc distances as claimed by B07 (though this is very unlikely), then the effect would not be distinguishable from other astrophysical and observational sources of apparent disk fraction gradients. We therefore conclude that the disk destruction result reported by B07 is not reliable, and our ECDF analysis of their data appears to refute it. Figure 11 shows results for NGC 6357. Both the diskbearing and disk-free distributions closely follow the predicted parabolas. Disk-bearing YSOs even show a slight overabundance at short distances to early O stars, possibly reflecting an age gradient in the region, or simply observational biases.

NGC 6357
3.3.6. Orion Nebula Figure 12 shows results for the Orion Nebula. Unlike the other five regions previously discussed in this subsection, the Orion Nebula contains no early-type O stars. θ 1 Ori C is omitted from our analysis due to PAH contamination. Both the disk-bearing and disk-free YSOs show a tendency to cluster around the one non-PAH excluded late O star (spectral type O9.5) in the region, as both ECDFs are significantly steeper than the predicted parabolic CDF at short distances. This is the clearest non-detection of inner disk destruction (within MYStIX completeness limits) due to the large density of disk-bearing YSOs around the late O stars (i.e., the very small steps in the disk-bearing ECDF).
It is interesting to note that the majority of the 178 Orion proplyds identified by Ricci et al. (2008) with Hubble appear in the MYStIX catalog as infrared excess (i.e., disk-bearing) sources in the MYStIX sample, with the notable exception of about a dozen proplyds close to θ 1 Ori C not seen by MYStIX due to PAH contamination. Hence, many of the disk-bearing YSOs near O stars in the MYStIX sample are likely in the process of being photoevaporated. The test for complete disk destruction in the current work can at most only place upper limits on photoevaporative mass-loss rates. Direct estimation of mass-loss rates requires either spectroscopic analysis of individual disks (e.g., Henney & O'Dell 1999)   O8.5V( (f)) O8V((f)) O5V O9.5V Figure 9. RA-Dec map of NGC 2244 (Rosette Nebula), using data from Balog et al. (2007) and 8.0 µm Spitzer /IRAC image. There are 336 disk-bearing YSOs (red circles), 705 disk-free YSOs (blue circles), and 5 O stars (large green circles), labeled by spectral types from Balog et al. (2007). Disk fraction in the cluster is highly non-uniform, suggesting that the destruction effect reported by Balog et al. (2007) is instead the result of disk fraction gradients.
Our study of protoplanetary disk destruction by external photoevaporation has several advantages over past studies, but also several limitations. In its favor, our sample covers 20 massive star-forming regions out to ∼3 kpc distance, allowing us to examine disks around a large sample of O stars, including the hottest and most luminous O2-O5 stars where disk destruction should be most effective. The MYStIX survey provides an unusually large sample of both high and low mass stars so that even a slight but consistent disk destruction effect in a given region could be discerned when averaged over many O stars (our data, which provide only a binary de-termination of disk-bearing versus disk-free, cannot reliably detect disk destruction around individual O stars due to the relative sparseness of surrounding disk-bearing YSOs).
The size of the MYStIX dataset notwithstanding, our search for spatial avoidance of O stars by disks faces numerous difficulties: non-detection of disks due to PAH contamination; incompleteness of the disk sample; nondetection of disks due to crowding; and astrophysical gradients in the spatial distribution of disks due to clustering behavior, age gradients, etc. The problem of crowding in dense, distant regions may be the most significant difficulty in detecting disk destruction with MYStIX, as it is precisely these denselypopulated regions where disk destruction is most likely to operate (due to small interstellar distances). In the case of the MYStIX sample, Trumpler 14 in the Carina Nebula and NGC 6618 in M 17 are two of the most densely-populated regions with some of the earliest O stars, and they are both located far away (∼2 kpc), making their members even more difficult to resolve. It is for the aforementioned reasons that, using MYStIX data as well as the data of Guarcello et al. (2009) andBalog et al. (2007), we find that the apparent distribution of disk-bearing YSOs can vary significantly across a given region on spatial scales similar to those associated with external photoevaporation by massive stars. All of these effects are in principle capable of giving rise both to false positive and false negative detections of disk destruction.
There are several ways in which the use of point source data in identifying spatial avoidance of massive stars by disks could be improved in future works. While nebular contamination will always make disk detection difficult, the problem of stellar crowding, due in part to large region distances, could be alleviated with high-resolution (PSF 1") infrared photometry of these dense cluster cores using 8 m-class telescopes with adaptive optics, or with the space-based James Webb Space Telescope. The dense OB association in Trumpler 14 (Carina Nebula) is highly crowded but has relatively little nebulosity, and would therefore make an excellent target for highresolution observations. Deeper observations would also be useful for detect disks around lower mass stars, particularly at high resolutions given the large numbers of low mass stars in cores of dense, massive star-forming regions. Such observations would be particularly useful for examining disk destruction in some of the regions discussed in Section 3.2, with low apparent surface densities of disk-bearing YSOs which may in many cases be the result of observational limits rather than intrinsically low number densities.
Another way of improving spatial avoidance-based detection of disk destruction would be the incorporation of far-infrared observations. Photometry from Spitzer /MIPS and Herschel, for instance, could probe disks at greater radii (tens of AU), though the relatively low spatial resolution of these instruments would make disk detection difficult in dense, distant regions. Nonetheless, far-infrared observations could help to resolve questions of how far photoevaporation proceeds into a given disk, which potentially has important implications for the evolution of planetary systems (discussed further in the following section).

Astrophysical considerations
Our non-detection of disk destruction around dozens of O stars in several densely-populated massive starforming regions is consistent with Orion Nebula studies that suggest that external photoevaporation by massive stars does not significantly suppress Solar Systemscale planet formation (Eisner & Carpenter 2006;Mann et al. 2014). , for instance, find that the fraction of disks with masses comparable to the Minimum Mass Solar Nebula in the Trapezium Cluster is similar to the fraction found in regions not dominated by O stars. Moreover, several theoretical works predict that external photoevaporation, as it disperses the disk starting at the outer edge and working its way inward, ceases to significantly operate at radii of a few tens of AU from the disk's host star where accretion timescales may still be relatively large (Clarke 2007;Adams et al. 2004Adams et al. , 2006Fatuzzo & Adams 2008;). This particular prediction is supported by findings based on near-/midinfrared excesses, namely, the current work and Oliveira et al. (2005), who find that the disk fraction in NGC 6611 (Eagle Nebula) is consistent with that of regions of similar age but without O stars. Anderson et al. (2013) explore disk models which combine external FUV photoevaporation and viscous accretion. For a range of Shakura-Sunyaev α viscosities (Shakura & Sunyaev 1973) and local FUV field strengths comparable to or greater than those of the Orion proplyds, they predict disk lifetimes between 5 × 10 5 yr and 10 7 yr, which are not significantly different from disk lifetimes based on viscous accretion alone. Their predictions based on nominal values of disk α and FUV field strength are roughly consistent with the results of the current work, however stronger constraints on photoevaporative mass loss rates would require the study of somewhat older regions (as our results only rule out the largest values of disk viscosity and FUV field strength). The Carina Nebula, whose members range from ∼1-4 Myr in age and which contains a significant number of both early and late O stars (as well as dense clusters of O stars such as in Trumpler 14), is therefore a prime target for further studies of disk destruction.
While the effect of external photoevaporation on the inner disk appears limited, it is possible that external photoevaporation suppresses the formation of planetesimals, including comets and the seeds of massive planets, in the outer disk (beyond 10-30 AU from the disk host star; Adams et al. 2004Adams et al. , 2006Fatuzzo & Adams 2008), or halts the inward migration of Jovian planets formed in the outer disk (or even contributes to outward migration due to positive torques from gas within the planet's orbit). Adams et al. (2004) suggest that external photoevaporation of gas in the outer disk in the early Solar System may explain the relative gas poorness of Uranus and Neptune, although other explanations have been offered Thommes et al. (e.g., 2002).
The observed ubiquity of exoplanets (Borucki et al. 2011;Petigura et al. 2013;Tuomi et al. 2014) clearly suggests that the overall effect of external photoevaporation on planet formation must be limited, either because complete disk destruction by massive stars is rare, or because planets are generally able to form before disk destruction occurs (the findings of Najita & Kenyon 2014, for instance, suggest that planet formation during protostellar collapse is required to explain the observed population of exoplanets). Nonetheless, understanding the effect of external photoevaporation on planet formation and migration in detail will require stronger constraints on the timescales of external photoevaporation as well as of planet formation and migration. In terms of the former, while the methods used in the current work provide loose constraints on photoevaporation timescales, other more direct methods for constraining photoevaporative mass loss rates are available. Henney & O'Dell (1999), for instance, use high-resolution spectroscopy to measure proplyd sizes and outflow velocities. Direct measurements of disk masses using sub-mm observations (Mann & Williams 2009b, 2010Mann et al. 2015), can be in principle be coupled with age estimates to calculate mass loss rates as a function of projected distance from O stars. Mass loss rates are relatively constant over a disk's lifetime (Anderson et al. 2013), therefore disk destruction timescales can be calculated as the ratio of assumed initial disk masses to observed mass loss rates. Many infrared-excess objects in the MYStIX sample, particularly near early O stars in the regions discussed in Section 3.3, are good candidates for such observations.  O6-O9.5 stars (N=4) Figure B2. NGC 1893. ECDFs of distances from massive stars to disk-bearing (red line) and disk-free (blue line) YSOs. Left panel is for early O stars, right panel is for late O stars. Dashed black line is a parabola, indicating the predicted CDF for two-dimensional geometry. The density of disk-bearing YSOs in this region is too small to detect disk avoidance within the destruction radii predicted by theory. O6-O9.5 stars (N=1) Figure B3. NGC 2264. ECDFs of distances from massive stars to disk-bearing (red line) and disk-free (blue line) YSOs. Left panel is for early O stars, right panel is for late O stars. Dashed black line is a parabola, indicating the predicted CDF for two-dimensional geometry. The density of disk-bearing YSOs in this region is too small to detect disk avoidance within the destruction radii predicted by theory. O6-O9.5 stars (N=2) Figure B4. NGC 3576. ECDFs of distances from massive stars to disk-bearing (red line) and disk-free (blue line) YSOs. Left panel is for early O stars, right panel is for late O stars. Dashed black line is a parabola, indicating the predicted CDF for two-dimensional geometry. The density of disk-bearing YSOs in this region is too small to detect disk avoidance within the destruction radii predicted by theory.  Figure B5. NGC 6334. ECDFs of distances from massive stars to disk-bearing (red line) and disk-free (blue line) YSOs. Left panel is for early O stars, right panel is for late O stars. Dashed black line is a parabola, indicating the predicted CDF for two-dimensional geometry. The density of disk-bearing YSOs in this region is too small to detect disk avoidance within the destruction radii predicted by theory. O6-O9.5 stars (N=2) Figure B6. W 3. ECDFs of distances from massive stars to disk-bearing (red line) and disk-free (blue line) YSOs. Left panel is for early O stars, right panel is for late O stars. Dashed black line is a parabola, indicating the predicted CDF for two-dimensional geometry. The density of disk-bearing YSOs in this region is too small to detect disk avoidance within the destruction radii predicted by theory.  Figure B7. W 4. ECDFs of distances from massive stars to disk-bearing (red line) and disk-free (blue line) YSOs. Left panel is for early O stars, right panel is for late O stars. Dashed black line is a parabola, indicating the predicted CDF for two-dimensional geometry. The density of disk-bearing YSOs in this region is too small to detect disk avoidance within the destruction radii predicted by theory.