© ESO, 2012

1. Introduction

Massive stars play a crucial role in shaping their environment by ionising large regions around them, affecting the temperature and structure of the ISM and through chemical enrichment. They are the flag-posts of star formation, with large H ii regions often dominating the optical structure of their host galaxies, e.g., 30 Doradus in the Large Magellanic Cloud (LMC). Despite their prominence, we know little about how massive stars form due to multiple challenges: their large distances to us, heavy extinction, rapid formation process, and small population statistics compared to their lower mass counterparts.

The spatial distribution of massive stars can provide us with clues to their formation. Do they only form in massive clusters or can they also form in isolation? There are multiple scenarios on how massive stars form and here we consider two of them, as these represent opposite ends of the spectrum. The first is competitive accretion (e.g., Bonnell et al. 2001, 2004; Smith et al. 2009) where stars, including massive ones, only form in clustered environments. In this model, massive stars are built up from the seeds of lower mass stars through the accretion of gas from their environment. Since a large amount of gas is needed for accretion, a cluster of low mass stars is expected to be present around the higher mass stars (e.g.,  Maschberger et al. 2010). The second model of massive star formation is that of monolithic collapse (Yorke & Sonnhalter 2002; McKee & Tan 2003; Krumholz et al. 2009) where a star’s mass is set by the initial dense core from which it will form. In this scenario, massive stars can form, albeit rarely, without a surrounding cluster. Hence, determining the spatial distributions of the massive stars allows us to test these theories and will possibly lead to a better understanding of the star formation process.

In addition to testing star-formation theories, determining whether or not massive stars can form in isolation would have important ramifications on how the stellar initial mass function (IMF) is sampled (e.g. Bastian et al. 2010). Depending on how stars form they will sample the IMF in different ways. It has been suggested that in the competitive accretion scenario massive stars only form in clusters, so the IMF will be sampled in a “sorted” way, such that enough low mass stars need to be present before higher mass stars can form (Weidner et al. 2010). However, in the monolithic collapse scenario, massive stars are able to form in relative isolation, meaning that the IMF will be sampled stochastically (e.g., Oey et al. 2004; Elmegreen 2006; Parker & Goodwin 2007; Selman & Melnick 2008), likely reflecting the mass distribution of the star-forming cores (e.g., Alves et al. 2007). These different ways of sampling the IMF have important implications to the resulting mass distributions, and in particular massive star numbers, in different galaxies (e.g. Weidner et al. 2010; Bastian et al. 2010; also see Sect. 6).

These two IMF sampling scenarios lead to similar relationships between the mass of the most massive star in a group and the mass of the group in total. However, there are important differences. In the stochastic scenario, such a relation is statistically expected (with a large scatter), whereas in the sorted sampling scenario such a relationship is causally expected with little scatter. Observations have been presented which seem to favour one scenario or another (e.g., Maschberger & Clarke 2008; Weidner et al. 2010), however large homogeneous datasets are required to provide a solid answer. Alternatively, in the present work, we will search for extremes in the distribution, O-type stars that formed in relative isolation, in order to address this issue.

The term “isolated massive star-formation” does leave room for different interpretations and one could argue against a star being “isolated” in several different ways, which we list below, from most to least restrictive.

• 1.

  The massive stars formed from the same molecular cloud or filament (i.e. the stars would be associated with the same OB associations).

• 2.

  The massive stars formed in the same or a related environment in terms of radiation field or kinetic energy input (e.g., a second massive star formed due to feedback from a first star).

• 3.

  The massive stars formed within the same gravitational well of size  ≲ 3 pc that led them to be (at least initially) bound and constitute a physical cluster (see Efremov & Elmegreen 1998; Maíz-Apellániz 2001; Scheepmaker et al. 2007) and not simply part of an unbound association.

In general terms, case 1 is responsible for large-scale similarities in age, location, and composition but is otherwise irrelevant for the detailed physics of star formation. Case 2 influences the overall efficiency and the age distribution within a large region and possibly the IMF, since triggering can, in principle, lead to variations. In case 3, only within the dense (sub)parsec scale regions, can gravitational interactions between different clumps and/or stars happen within the massive-star formation timescale of  ~10⁵ years (McKee & Tan 2003). In other words, from the gravitational point of view, it is mostly irrelevant whether a large mass (molecular cloud or other clusters) is located at a distance of 10 pc because for a period of 10⁵ years this would constitute a nearly constant gravitational field. It is only nearby (~1 pc) parts of the cloud that have orbital time scales short enough to lead to multiple interactions, such as accretion, close encounters, and collisions. Therefore, from the point of view of star-formation theories and the origin of the IMF, case 3 is the relevant scale, which is what we will adopt in the current paper.

The possible effect of triggering (case 2) could work as an external agent to influencing massive star-formation in 30 Doradus. The central cluster, R136, other clusters and massive stars have been suggested to influence the massive star-formation process in the region (see Brandner et al. 2001; Walborn et al. 1999, 2002), but De Marchi et al. (2011) questions whether triggered star-formation is relevant as their investigation with the Hubble Space Telescope’s (HST) Wide Field Camera 3 (WFC3) observations do not show clear evidence of such causal effects. In either scenario, whether triggering (case 2) is influencing massive star-formation or not in 30 Doradus, it would not change a massive star’s candidacy for forming in isolation following the case 3 definition. Once triggering activates star-formation it is up to the parts of the cloud within the cluster length scale to interact with each other.

When we talk about massive stars and them forming in “isolation”, we adopt a similar¹ definition given in Parker & Goodwin (2007), where massive stars refer to those those with spectral types of O7 or earlier ( ≳ 30 M_⊙ depending on mass estimate method used – discussed later on in the paper) and isolation means that they are not found in clusters of  ≥ 100   M_⊙ and r ≲ 3 pc. Once massive stars are located, we will estimate the expected underlying cluster mass according to the (Weidner et al. 2010) theory by using the m_max − M_cl relation, which is expected to hold in the sorted sampling scenario. This allows us to assess if the underlying cluster should have been detected, and if the O-type stars are truly isolated, will enable constraints to be put on star-formation theories as well as IMF sampling scenarios.

We are not the first to look for massive stars forming in isolation. A comprehensive analysis of isolated field O-type stars (hereafter referred to as O-stars) in our Galaxy was conducted by de Wit et al. (2004, 2005) who found that 4  ±  2 percent of the O-stars in their sample (model derived value from observations) could not be traced back to clusters and hence, likely formed in isolation. Similarly, Lamb et al. (2010) looked at the Small Magellanic Cloud (SMC) and found three O-stars that are in sparse-clusters. The total mass of the sparse-clusters in relation to the mass of the O-stars is not compatible with the sorted sampling scenario. The sparse-clusters fall within the definition as we defined above for “isolated”. The SMC was also investigated for isolated star formation by Selier et al. (2011), where they reported an interesting compact H   II region, N33. The compact H   II region cannot be traced back to any nearby clusters, associations or molecular clouds and may be evidence for isolated massive star formation at the earliest stages. Here we consider observations from the VLT-FLAMES Tarantula Survey (hereafter VFTS; Evans et al. 2011) of the 30 Doradus nebula and its environs which contains 350 O-stars.

In addition to the de Wit et al. (2004, 2005) findings, there have been several theoretical investigations of the likelihood of isolated, massive star formation. Firstly, Oey et al. (2004) noted that the number of stars per cluster appears to follow a power-law like slope. They conclude that the power-law dependence extends, continuously, all the way to down to one OB star per group, association or cluster, i.e. N_⋆ = 1. Parker & Goodwin (2007) simulated a large number of stellar clusters assuming a cluster mass function with β = 2 and a universal IMF (stochastic sampling). They found that “isolated” massive star formation is not unexpected in this scenario, and that contrary to Weidner & Kroupa (2006) who found a strict m_max − M_cl relation, stochastic sampling leads to a m_max − M_cl relation with significantly more scatter and a slightly different slope. In the sorted sampling scenario, the scatter is significantly reduced, due to the explicit link between the cluster mass and the stars that it forms. However, without a large number of clusters, it is difficult to differentiate between the two (Maschberger & Clarke 2008). The two scenarios do differ significantly in the extreme end of the distributions, i.e. the presence/absence of isolated high mass stars. Lamb et al. (2010) conducted similar and updated simulations as presented by Parker & Goodwin (2007) and confirmed their results.

It is important to note that Weidner & Kroupa (2006) and Weidner et al. (2010) adopt similar scales from case 3 (concerning arguments against “isolated” stars) for defining the m_max − M_cl relation, i.e., they adopted the cluster length scale and not the OB association scale. For example, the authors use the Orion Nebular Cluster (ONC) and not the entire Orion complex (an OB association) for supporting the m_max − M_cl relation. Using the cluster scale makes sense if all stars are formed in clusters, i.e., in quantised units of star-formation. However, if star-formation is hierarchical (e.g., Elmegreen et al. 2006) then no distinct scale exists in the star-formation process (above the scale of individual stars/cores). This means that clusters are made up of sub-clusters which merge as they evolve dynamically, and that clusters/associations are themselves grouped into larger structures (e.g., Bastian et al. 2005). In a hierarchical scenario, the individual sub-clusters cannot fit the m_max − M_cl relation if the final cluster does (since the most massive star would be too massive for its sub-cluster size), so the stars must “know” about the cluster that they will finally be a part of. This same argument also holds for larger associations and cluster complexes. Hence, the m_max − M_cl relation must have a length scale associated to it, if star-formation is hierarchical, above and below which the relation breaks down.

The 30 Doradus region is a complex and dynamic region that contains multiple, although not necessarily spatially distinct, generations of stars (see Walborn & Blades 1997). The youngest population is dominated by the central cluster R136, with an age of 1–2 Myr (e.g. de Koter et al. 1998; Massey & Hunter 1998) and, most pertinently for the discussion here, there appears to be another young population to the north and west of R136, exemplified by the compact multiple systems in the dense nebular knots observed with the HST by Walborn et al. (1999, 2002). These comprise an apparently young, still embedded phase of star formation. Interestingly, Walborn et al. (2002) also presented imaging of two notable infrared (IR) sources, one of which was resolved into a small, embedded cluster, while the other was a point like source, seemingly single monolithic object (albeit at the distance of the LMC). This led the authors to note that the later object may have formed without an associated cluster or association.

R136 has a relatively shallow power-law density profile and it does not appear to have a strong truncation out to at least 10 pc (e.g. Campbell et al. 2010). All massive stars observed beyond this radius are likely to be either runaway stars, or to have formed in situ. The majority of the stars in our sample are young (2 Myr), hence the disruption of large clusters with lower-mass stars blending into the background due to “infant mortality”, the “cruel cradle effect”² or the dispersal of unbound associations (either of which would cause the most massive star to appear isolated), is not expected to influence our analysis (e.g. Bastian & Goodwin 2006; Gieles & Portegies Zwart 2011; Kruijssen et al. 2011, 2012; Girichidis et al. 2012).

In the present work, we attempt to overcome the limitations of the previous studies by removing the possibility that the stars are runaways by using multi-epoch medium-resolution spectroscopy and by cross-correlating the spatial distribution of the candidates with known gaseous filaments. We do this by combining observations of 30 Doradus with several instruments from the Very Large Telescope (VLT), HST, and the Spitzer Space Telescope (SST), which will be discussed in Sect. 2. The methods to finding isolated massive star candidates are discussed in Sect. 3, where the results and conclusions are discussed in Sects. 4, 5. In the final section we provide a summary and discuss implications.

2. Observations

In this work we used optical spectroscopy and optical to mid-IR imaging to study the O-stars in 30 Doradus. The core dataset is composed of medium resolution optical spectra from the VLT-FLAMES Tarantula Survey (Evans et al. 2011). bf VFTS allows us to identify the spectral type of each star, test whether the star is part of a spectroscopic binary and measure each star’s radial velocity. The survey employs three different modes of the Fibre Large Array Multi-Element Spectrograph (FLAMES; Pasquini et al. 2002) instrument on the VLT: Medusa-Giraffe, ARGUS-Giraffe, and the Ultraviolet and Visual Echelle Spectrograph. In this paper we only use spectra that were observed using the Medusa fibre-feed to the Giraffe spectrograph. There are 132 fibres available for observations, deployable within a 25′ diameter field-of-view and with a diameter of $\mathrm{1}\stackrel{\mathrm{″}}{.}\mathrm{2}$ on the sky. The European Southern Observatory (ESO) Common Pipeline Library (CPL) FLAMES reduction routines (v2.8.7) were used for all of the data processing. The standard reductions were then applied such as Heliocentric correction and sky subtraction. Additional information regarding the observations, reductions and survey strategy are given in Evans et al. (2011). The radial velocity and multiplicity analysis is discussed in Sana et al. (in prep.).

The second part of the dataset is composed of imaging surveys to 1) find filamentary structures of gas and dust, i.e., the sites of star formation and 2) determine if a given star is truly isolated. For the first goal we use the Hα/Hβ map from Lazendic et al. (2003), 70 μm maps from the SST Legacy survey Surveying the Agents of Galaxy Evolution (SAGE) (Meixner et al. 2006), and VLT High Acuity Wide field K-band Imager (HAWK-I) observations (Kissler-Patig et al. 2008) where ionised filaments can be detected due to the presence of the Brγ emission line within the bandpass. To determine whether there are significant clusters associated with the massive stars we use HST imaging with the Wide-Field Planetary Camera 2 (WFPC2) and the Advanced Camera for Surveys (ACS) instruments as well as the VLT HAWK-I K_s band images to search for embedded clusters.

In Fig. 1 we show a V-band image of 30 Doradus taken with ESO and the Max Planck Gesellschaft (MPG) 2.2 m telescope using the Wide Field Imager (WFI) instrument (Baade et al. 1999) to show the 30 Doradus region for context (Program ID: 076.C-0888; PI J. Alves).

For the candidates, except VFTS 089 and 849, we acquired the HST data from the Hubble Legacy Archive, and selected the science grade data for analysis. The proposal ID numbers for the HST data are: 05114, 08163, 09471. For VFTS 089 and 849 we used high-level science products from The Archival Pure Parallels Project (APPP) on the LMC. Wadadekar et al. (2006) describes the data processing and quality of the HST images. WFPC2 and ACS data typically have pixel resolutions of 0.1′′ and 0.05′′, respectively. See Fig. 2 for 2.5  ×  2.5 pc subplots of the HST data for each candidate. Further details on HST’s sensitivity around the O-star candidates is discussed in Sect. 4.2.

3. Method

To determine that an O-star is a not a runaway and is still located near its birth-site, we need to apply several constraints on the star in the line-of-sight and the plane of the sky. Several observational tools are needed to do this: precise radial velocity measurements (line-of-sight), association with filamentary structure (plane of sky), and relative isolation from lower mass stars.

The radial velocity measurements of all the single O-stars in VFTS (Sana et al., in prep.) can be well fit by a gaussian distribution. The mean radial velocity is 270.73 $\mathrm{km}\hspace{0.17em}\mathrm{s}{-1}^{}$ with a dispersion (1σ) of 10.50 $\mathrm{km}\hspace{0.17em}\mathrm{s}{-1}^{}$. These values are used in selecting the candidates (see criterion 3).

3.1. Criteria

We begin with the VFTS observations of the 30 Doradus region. The survey will provide multi-epoch radial velocity measurements for a sample of  ~800 OB-type stars in the 30 Doradus region with a typical precision better than 5 $\mathrm{km}\hspace{0.17em}\mathrm{s}{-1}^{}$. From the sample we select candidates based on the selection criteria provided below. Within brackets we report on the number of stars that fulfil all successive criteria in the order as listed and which constraint it supports (e.g. plane-of-sky, line-of-sight, or isolation). The criteria are:

• 1.

  O-type stars that do not show significant variations in their RV, meaning that they are unlikely part of a spectroscopic binary system. This allows us to determine accurately the systemic velocity of each star. In Sect. 4.3, we quantitatively estimate our detection biases and we will show that our sample cannot be heavily contaminated by undetected spectroscopic binaries (stars: 184, which include some B0-type stars from preliminary classifications) (constraint: line-of-sight).

• 2.

  Their sub-spectral type are O7 or earlier, which corresponds so stellar masses of  ≳ 30 M_⊙. See Sect. 3.2 for details on mass estimates (stars: 65).

• 3.

  Their radial velocities are within one σ of the mean of all the O-stars in the VFTS sample of 30 Doradus. This ensures that the stars are not runaways in the line-of-sight (stars: 39) (constraint: line-of-sight).

• That the candidate star be located on (projected) gaseous filament seen in either ionised gas (i.e., Hα or Brγ), warm dust (i.e., 70 μm), or cold gas/dust (i.e., based on extinction or molecular maps). These gaseous filaments are the likely sites of star-formation, and the chance projection of all ejected stars from nearby clusters that lie projected upon these filaments is exceedingly low (stars: 27) (constraint: plane-of-sky).

  Fig. 3 The 90% completeness magnitude limit as a function of radius for VFTS 385 (using HST/WFPC2). The closest detectable source to the candidate is  ~0.15 pc, where the upturn begins. VFTS 385 is the brightest observed candidate presented in this paper.

  Fig. 4 HST stellar surface density distributions (cumulative – see Lamb et al. 2010) around VFTS 208. The black line with triangles is the HST F555W band and the blue line with circles is the HST F814W band. Due to the saturation on the ACS images caused by the candidate (candidate reference point is at 0 pc), VFTS 208 represents a worst case. However, even here we see that there is no stellar count increase from 3 pc and inwards (to 0 pc) toward the candidate. This has been observed similarly for each candidate.

In addition to the above criteria, we use high resolution HST imaging (when available) of the regions around the candidate stars, along with deep ground-based near-infrared imaging from VLT to place constraints on the size of any potential cluster surrounding the candidate O-star. With this data we add the following criterium:

• 5.

  The surface density distribution of stars does not increase towards the candidate star within 3 pc, see Fig. 4. Note that archival HST data does not cover all candidates from criterion 4. Due to the lack of HST data, four sources are automatically rejected from consideration. This means that this step automatically reduced the number of stars from the previous step of 27 stars to 23 stars. Then we apply the surface density criterion on the remaining 23 stars (stars: 15) (constraint: isolation).

Where applicable, if two stars are located along the same gaseous filament, we apply the additional criterion below. Note, that the additional criterion does not exclude the 15 candidates from above.

• 6.

  Two stars have radial velocities within 5 $\mathrm{km}\hspace{0.17em}\mathrm{s}{-1}^{}$ of each other. While filaments can have flows along them, we expect the gas/stars within a filament to have very similar radial velocities (stars: 2) (constraints: plane-of-sky, line-of-sight).

With these criteria we can build up a collection of candidates that may have formed in isolation. Once a homogeneous map of 30 Doradus showing the gas/dust filamentary structures become available, e.g. Herschel Space Observatory maps (Meixner et al. 2010) or extinction maps derived from the Visible and Infrared Survey Telescope for Astronomy (VISTA) Magellenic Cloud Survey (Cioni et al. 2011), we will be able to apply simple probability tests to strongly limit the possibility of massive stars being runaways (see Appendix B).

3.2. Spectral types and ages

The ages of our isolated O-star candidates likely range from less than 1 Myr to more than 4 Myr based on their spectral types (Weidner & Vink 2010). We have therefore assigned grades to each star on the basis of their spectral types, employing estimated ages from Weidner & Vink (2010) for LMC stars. Grade 1 candidates are most likely 2 Myr old or younger, grade 2 are between 2 and 4 Myr, and Grade 3 are those older than 4 Myr. The grade scheme is provided to highlight possible issues of candidate associations, depending on their age, with filamentary gas/dust structure (see Sect. 5.4).

4. Results

The properties of the 15 isolated O-star candidates are given in Table 1, including their positions, mean radial velocities (of the gas and stars), spectral types (Walborn et al., in prep.), masses (estimated from the spectral types), V-band magnitudes, their expected cluster masses based on the m_max − M_cl relation, and a grade to indicate their likely ages. Eleven of the candidates are classified as Vz stars, dened by stronger He ii λ4686 absorption than either He i λ4471 or He ii λ4542 in their spectra (Walborn & Parker 1992; Walborn & Blades 1997). The Vz phenomenon has been proposed to be related to stars located close to the zero-age main sequence (ZAMS). Alternatively, rotational broadening may induce an apparent modification of the line strength and of He line ratios. In the latter case the Vz signature may be common for normal dwarf stars and not be tracing stars close to the ZAMS (Sabín-Sanjulián, priv. comm.). Quantitative spectroscopic analyses are needed to elucidate the true nature of the Vz signature and will be undertaken for the Vz stars in the VFTS sample.

Eleven of the candidates were classified as Vz stars, defined by strong He ii λ4686 absorption in their spectra (Walborn & Parker 1992; Walborn & Blades 1997) and thought to be close to the zero-age main sequence (ZAMS). The ZAMS classification is not absolute since the morphological hypothesis has not yet been confirmed by systematic quantitative analyses; determinations of their gravities and luminosities is underway as part of the VFTS.

The candidates are located at projected distances of between 14 and 130 pc from R136 (see Fig. 1). Twelve are associated with the active star forming environment around R136 (see Fig. 5 for example). Four are located to the southwest of R136 and are closer to the older association NGC 2060, which has an associated supernova remnant (Danziger et al. 1981; Chu et al. 1992). Note that the distance between the centre of NGC 2060 and the closest candidate in this paper is  >20 pc.

We have taken the spectral types of the candidates and estimated their respective masses from Weidner & Vink (2010), for rotating O-star models in the LMC to approximate the expected O-star parent cluster’s mass based on m_max − M_cl. The stellar masses are also used as an input in the computation of the sensitivity of our data to spectroscopic binaries (see Sect. 4.3). Estimating stellar masses from spectral types is known to suffer from a number of caveats (e.g. Weidner & Vink 2010), however the masses are precise enough for our objectives (typical errors are expected to be less than 30%). In the following subsections we discuss the candidates in more detail.

4.1. Local environment examples

4.1.1. VFTS 385, 398, 470, and 581

From east-to-west in Fig. 5 we have VFTS 581, 470, 398 and 385. We see that they are associated to filamentary structures from the VLT HAWK-I K_s band, the Hα/Hβ extinction map (Lazendic et al. 2003), and the 70 μm map. The observation of VFTS 581 will be specifically modelled in Sect. 5.1 in order to see how visible its underlying cluster would be. VFTS 398 and 470 are both classified as Vz stars.

4.1.2. VFTS 208 and 216

VFTS 208 and 216 are found to be associated with a large filamentary structure (Fig. 6) in the 70   μm Spitzer band, which is associated with NGC 2060, as seen in Fig. 1. No significant stellar populations are seen nearby, implying that they are relatively isolated. There may be a possible bow shock or bubble related to VFTS 216, which is shown in Fig. 7. If there is a bow shock associated with the star, this would strongly suggest that it is a runaway, and did not form in isolation. The stellar surface distribution (based on HST V and I-band imaging) of VFTS 208 is shown in Fig. 4. We note that VFTS 208 is classified as an nfp star (Walborn 1973). While the origins for the spectral peculiarities seen in such stars are still unknown (see Walborn et al. 2010), we retain the star in our discussion as it appears both relatively massive and young.

4.1.3. VFTS 488 and 706

VFTS 488 lies south of R136 at the edge of the dense 70 μm emission (see in Fig. 8). It appears that VFTS 488 is associated with a filament, and both the VLT K_s and HST F814W bands show its relative isolation from star clusters.

Northeast of R136 is a prominent ionised filament structure that is seen both in the V and K_s band. VFTS 706 is at the centre of the filament (see Fig. 9). It can be seen in 70 μm that the region is associated cold dust. Both candidates are Vz stars.

4.1.4. VFTS 682

Northeast of R136 we have VFTS 682, which does not fit our criteria as our best estimate of its radial velocity (from models of its emission lines) is outside the standard deviation from the mean. However, its young age and exceptionally high mass without evidence of clusters merits attention. Evans et al. (2011) discovered VFTS 682 to be a, previously unknown, Wolf–Rayet star (see Fig. 10). They also noted it was within 4 pc of VFTS 702, a candidate young stellar object (YSO) from Gruendl & Chu (2009). The star displays excess emission in the mid-IR although the cause is still unclear, however the correlation with molecular gas suggests a region of ongoing star-formation. VFTS 682 is the most massive star in our sample, with Bestenlehner et al. (2011) estimating an age ranging between 1 to 1.4 Myr (i.e. in our youngest age range) from comparisons with evolutionary models. Bestenlehner et al. (2011) estimate that the A_V > 4 for the candidate which could make it difficult to detect lower mass stars around the candidate. However, according to the m_max − M_cl relation, there should be a cluster of  ~3900 M_⊙ around it. This candidate is specifically modelled in Sect. 5.1, and such a cluster would have easily been detected. Banerjee et al. (2012) proposes that VFTS 682 could be a slow runaway, but is only possible under their model assuming complete mass-segregation in the parenting cluster (R136) and that all massive stars are in binaries. Additionally, we find no evidence of bow shocks around the candidate at 24 μm.

4.2. Photometric completeness around the candidates

In order to test our photometric completeness as a function of magnitude and distance from each candidate source, we have estimated the 90% completeness limit relative to the radius from the candidates. For this we used PSFEx² to generate artificial stars from modelled point spread functions (for each of the cameras used) and added them to the science frames (Fig. 2). The artificial sources were given magnitudes between 16.5 and 24 mag, and 49 sources were added to the image in each iteration (in order to avoid excess crowding). For each magnitude step of 0.5, 50 iterations were carried out (resulting in 2450 sources per magnitude added). SExtractor (Bertin & Arnouts 1996) was used to detect the sources, and the resulting catalogues were used to test the success rate of detections near the candidates. An example of the procedure is given in Fig. 3, where we see (for VFTS 385, which is a typical source) we are 90% complete at V = 23 mag at 0.26 pc from the candidate O-star. This corresponds to a  ~2 M_⊙ star assuming an extinction estimate of A_V = 1.

Additionally, the cumulative surface densities of sources as a function of radius away from each candidate was investigated using HST V and I-band imaging. An example is shown in Fig. 4, which shows no evidence of stellar density increase around VFTS 208. For each source we adopted a conservative completeness limit (generally in the range of V = 21 to 23 mag) and focussed on radii greater than 1 pc. If a cluster was present, we would expect a rise in the number of stars per pc² from  ~3 pc and inwards (towards 0 pc, the reference frame of the candidate). Similar results are seen for all the other candidates where no stellar density increase is observed. The information from this section will be used further in Sect. 5.1 to explore the constraints on any underlying cluster that may be present.

Upon visual inspection of higher resolution HST data from the WFC3 (De Marchi et al. 2011) of VFTS 385, 392, 577 and 706 we see no evidence of clustering in the F555W and F814W filters. We do not present the data itself in this paper.

4.3. Binary detection probability

To estimate the fraction of undetected spectroscopic binaries in our sample, we have quantified the observational biases of each candidate using the Monte Carlo methods from Sana et al. (2009). We estimated the probability (P_detect) to detect binarity in each candidate, adopting a primary mass on the basis of their spectral types and with binary properties randomly drawn from uniform cumulative distribution functions (CDFs) of mass-ratios and eccentricities, and a bi-modal CDF in log    P (as described by Sana & Evans 2011). Assuming a random orientation of the system in space and a random time of periastron passage, we then apply the specific observational sampling of each object, adding noise to the RV signal that corresponds to the error measurements at each epoch. Finally, we consider whether the simulated system would have been detected if the amplitude of the RV signal is larger than an adopted detection threshold.

For all the O-star candidates, except VFTS 706, we adopt a threshold of 20 $\mathrm{km}\hspace{0.17em}\mathrm{s}{-1}^{}$ from the mean. These levels are chosen to ensure no false detections at the 99.99% confidence level ( > 3σ). VFTS 706 shows evidence of faster rotation than the others, so we have fixed the threshold to 35 $\mathrm{km}\hspace{0.17em}\mathrm{s}{-1}^{}$, yielding no false detections at a 99.7% confidence. Table 2 provides the computed detection probabilities across three ranges of periods as well as across the full period range. This shows that it is unlikely that our sample contains undetected short and intermediate period spectroscopic binaries. For period larger than one year, our detection probabilities become relatively small. Under the adopted hypothesis, the average detection probability is 0.83    ±    0.05 up to a period of  ~ 8.6 yr. While 10 to 20% of the objects in our sample might still be undetected, likely long-period spectroscopic binaries, these simulations allow us to conclude that most of our targets are either single stars or very wide/large mass ratio pairs so that the measured systemic velocity is left unaffected by the companion. Additionally, the spectra were inspected for SB2 systems.

5. Discussion

5.1. Modelling potential underlying clusters

5.1.1. Visual examples of cluster presence

In order to check if we are sensitive to clusters that may be associated with the candidates according to the Weidner et al. (2010)m_max − M_cl relation, we use the MASSCLEAN package (Popescu & Hanson 2009) to generate illustrative examples. The package generates images (at a given distance and resolution), of simulated clusters with stars drawn from a Kroupa IMF (Kroupa 2001) and spatially distributed according to a King profile (King 1962). For the analysis we adopt King profile radii similar to those reported by Werchan & Zaritsky (2011) for the modelled clusters. Stellar masses are drawn such that the Weidner et al. (2010)m_max − M_cl relationship is recovered from these simulations. All simulations were carried out adopting the resolution of the ACS WFC.

We generated three clusters, whose most massive stars were 25, 45, and 100 M_⊙. The 25 M_⊙ star represents the lower mass estimates of the candidates given in Table 4, while the 45 and 100 M_⊙ stars represent VFTS 581 and 682 (Table 1), respectively. The total cluster masses are 420, 1000 and 4000 M_⊙, respectively, and are shown in Fig. 11. In order to account for the effects of extinction (which we assume affect all cluster members equally), we then scaled the images such that the most massive star had the same V-band apparent magnitude as the candidate which we were modelling.

In Fig. 11 the most massive star is indicated by a green diamond. Note that in the 4000 M_⊙ image, the most massive star is not in the centre of the cluster, this is simply due to the fact that the positions of each of the stars is distributed stochastically, regardless of mass. We then applied the completeness curve from Fig. 3, assuming that all stars above the 90% completeness limit (at a given radius) would have been detected. All stars that pass this criteria are circled (red) in the images. In the middle panel of Fig. 11 fewer stars would have been detected than in the left (420 M_⊙ case) panel. This is due to the difference in reddening between the two model stars (for similarity to the candidates).

In all cases we would likely have observed an underlying cluster, which is evident when comparing Figs. 2 and 11, suggesting that these stars are isolated compared to what one would expect for star formation drawn from sorted sampling. Each of the O-star candidates should have B-stars present and several should have other O-stars in their associated cluster.

5.1.2. Monte Carlo simulations

Alongside the visual examples of the typical clusters that should be present around the massive O-star candidates, we conduct three comprehensive sets of Monte Carlo simulations.The exact same parameters of the most massive stars and clusters given in the previous section are adopted, but with uniform dispersion of cluster masses for each simulation set (SS#). The cluster mass ranges are extrapolated from the 1/5 and 5/6 quantile dispersions reported from Fig. 5 in Weidner et al. (2010). Details of the SS# are listed below.

• 1.

  SS1: 10 000 runs with median cluster mass of 420 M_⊙ with its most massive star at 25 M_⊙. The cluster mass dispersion range is 10²−10³ M_⊙.

• 2.

  SS2: 10 000 runs with median cluster mass of 1000 M_⊙ with its most massive star at 40 M_⊙. The cluster mass dispersion range is 10^2.2−10^3.6 M_⊙.

• 3.

  SS3: 10 000 runs with median cluster mass of 4000 M_⊙ with its most massive star at 1000 M_⊙. The cluster mass dispersion range is 10^3.1 − 10^4.1 M_⊙.

To determine the observable number of stars in the simulations we remove stars below our estimated mass sensitivity limit of 3 M_⊙. Furthermore, we assign cluster positions to each of the stars using the same King profile parameters as mentioned in the previous section and remove any stars within 0.15 pc of the most massive star in the cluster. This reflects the completeness limit as a function of radius discussed in Sect. 4.2. In Fig. 12 the number of observable stars with a 1σ dispersion is shown. For the cases of SS1, SS2 and SS3 we observe at minimum 11, 53, and 148 numbers of stars at the 1σ dispersion level from the mean assuming A_V ~ 1.

If we place a system like the Trapezium cluster (just the inner 0.1 pc of the Orion Nebula Cluster) in the LMC we would not be able to resolve it. However, with an effective radius of 2 pc (e.g. Portegies Zwart et al. 2010) we would have readily detected the surrounding cluster with minimum of  ~40 excess stars in a 5 pc radius. This is particularly relevant for the m_max − M_cl relationship presented by Weidner et al. (2010), for which the full Orion Nebula Cluster must be used for the most massive star, θ¹Ori   C (48 M_⊙, Kraus et al. 2007) to fit the relation.

5.2. Comparing Monte Carlo simulations to observations

If the m_max − M_cl relationship is correct we should have an excess number of stars within a 5 pc radius of each of the candidates. Around each candidate we measured the number of observable stars above 3 M_⊙ in the HST data and compared them to the Monte Carlo simulation sets discussed in the previous section. Two apertures were used to estimate the number of excess stars relative to the background and foreground stars. An inner aperture of 5 pc and an field aperture of 10 pc, which are denoted as NN_inner and NN_field, respectively. We define the number of excess stars (field star corrected) in the inner aperture as N_excess = (NN_inner − NN_field) × π   r² where r = 5 pc. Upon comparing the number of excess stars to the number of expected stars if a cluster is present from the simulations (see Table 3) we see no evidence for clusters around the candidates.

5.3. Mass and age discrepancies

We use the spectral types of our candidates to derive their age and masses via Weidner & Vink (2010) as shown in Table 1, but these properties can also be obtained from other methods. To assess potential uncertainties we have also obtained physical parameters using the Galactic spectral type-T_eff calibration of Martins et al. (2005). Absolute visual magnitudes were obtained by correcting the photometry from Evans et al. (2011) for the effects of extinction, with stellar luminosities estimated from the bolometric correction calibration of Martins et al. (2005). Masses and ages follow from comparison with non-rotating, LMC metallicity stellar evolutionary models and isochrones (Meynet et al. 1994; Lejeune et al. 1997).

We acknowledge that the use of an alternative spectral type-T_eff calibration (e.g. Massey et al. 2009) may lead to increased stellar luminosities/masses. Conversely, decreased stellar masses would be inferred from contemporary evolutionary models allowing for rotation (Maeder & Meynet 2001; Brott et al. 2011). Results are presented in Table 4, although refined parameters await detailed analyses, which are currently in progress. Individual stellar masses are on average  <25% below those derived using spectral types (Weidner & Vink 2010). Nevertheless, our key conclusions are unaffected by the method used to derive masses since the parent clusters, if present, would still be observable. Except for VFTS 682, all candidates (including VFTS 208) would meet our criterion and be assigned grade 1 or 2 according their estimated ages under this method.

5.4. Filamentary structures in 30 Doradus

The 30 Doradus region is largely affected by R136. The gaseous filaments that we consider in method four for associating with the O-stars (see Sect. 3) could possibly be problematic for the slightly older stellar population. In particular, the grade 3 candidate, VFTS 208, since the filaments that such a star formed in, could have been moved/destroyed by feedback from the nearby R136 cluster (e.g. Tenorio-Tagle et al. 2006).

However, others argue that since the ISM is filamentary and heterogeneous (e.g. André et al. 2010; Bergin & Tafalla 2007), the effects of ionising sources on the gas will be less significant than what Tenorio-Tagle et al. (2006) show in their simulations (e.g. Dale & Bonnell 2011).

Note that as an additional constraint on the association between stars and gas, we measured line-of-sight gas velocities from the [NII]6583, [SII]6717, [SII]6731, and Hα nebular emission lines superimposed on the stellar spectrum (see Table 1). When multiple gas velocity components are identified, only the one that is closest to RV_star is presented. We can see that 7 of the 16 candidates have a radial velocity for the ionised gas (RV_ISM) within  ~5 $\mathrm{km}\hspace{0.17em}\mathrm{s}{-1}^{}$ of the mean RV_star, while the agreement between RV_ISM and RV_star is within 15 $\mathrm{km}\hspace{0.17em}\mathrm{s}{-1}^{}$ for most candidates.

5.5. Bow-shocks and the ISM in 30 Doradus

Gvaramadze et al. (2010) discovered bow shocks around two isolated massive stars a few hundred pc away from R136. They argue that these OB stars are consistent with being ejected from the stellar cluster. Their discovery proves that bow shocks can form in the 30 Doradus region (or at least in its surroundings) and that these structures can be detected in Spitzer 24 μm images. So, while the presence of a bow shock suggests the runaway nature of an isolated OB star (Gvaramadze et al. 2010), what can one conclude for the isolated candidates presented in this paper that do not show evidence of bow shocks? We can estimate the minimum required runaway velocity to explain the apparent isolation of the candidates considered in this paper. From the observed projected distance d from R136 (14 to 130 pc) and assuming an age t ≤ 1 Myr $\mathit{v}\mathrm{\ge }\mathrm{9.8} (\frac{\mathit{d}}{\mathrm{10} \mathrm{pc}})\hspace{0.17em}(\frac{\mathit{t}}{\mathrm{Myr}}{)}^{-1}\hspace{0.17em} \mathrm{km}\hspace{0.17em}\mathrm{s}-1$we obtain minimum runaway velocities between 14 and 130 $\mathrm{km}\hspace{0.17em}\mathrm{s}{-1}^{}$. These velocities are higher than the isothermal sound speed in a warm ISM. In the case of 30 Doradus the average gas temperature is estimated to be around 10⁴ K, with radiation pressure appearing to be negligible outside 10 pc from the central cluster R136 (Pellegrini et al. 2011; Peimbert 2003). This results in an isothermal sound speed of order 10 $\mathrm{km}\hspace{0.17em}\mathrm{s}{-1}^{}$. Inside the cavities containing hot, low density gas, the sound speed can increase up to 100 $\mathrm{km}\hspace{0.17em}\mathrm{s}{-1}^{}$. However such cavities are mostly observed in the inner region of 30 Doradus (Pellegrini et al. 2011). Therefore, if our isolated massive stars are runaways, they are likely to be moving supersonically with respect to the local ISM.

In a fraction of supersonically moving runaway OB stars, the interaction between the stellar wind and the ISM will produce observable bow shocks (e.g. Comeron & Kaper 1998). This is supported by observations of such objects, both in the Galaxy (e.g. van Buren et al. 1995; Huthoff & Kaper 2002; Gvaramadze & Bomans 2008) and in the Magellanic Clouds (Gvaramadze et al. 2010, 2011).

The absence of bow shocks in a large sample of isolated OB stars could hint to the fact that not all of them are runaways. Gvaramadze et al. (2010) found that  ~30% of the runaways in their sample in the LMC had bow shocks. This is consistent with Galactic studies of bow shocks as well where 20 − 40% of well established runaway OB stars has a detectable bow shock (van Buren et al. 1995; Huthoff & Kaper 2002). Hence, for the LMC we would expect  ~ 5 of the 16 candidates to have bow shocks if they are all runaways. While this argument is clearly not conclusive, in line with recent efforts of Gvaramadze et al. (2010, 2011), we suggest that the investigation of the presence/absence of bow shocks around large enough samples of isolated O-stars could be an important way to constrain their formation mechanism.

5.6. Number of isolated stars

What fraction of stars would we expect to form in isolation, generally? Not only does the answer depend on how the IMF is formed (e.g. sorted or stochastic sampling) but also on the intrinsic spatial distribution of stars at the time of their formation, and how this probability function is sampled. Historically, if one assumes that all stars are formed in “clusters”, then clusters are the basic unit of star-formation, from which stars are then sampled. However, if “clusters” in fact do not represent a basic unit of star-formation, but instead there exists a distribution of surface densities (where “clusters” merely represent the high-surface density tail of the distribution) then one needs to know how the spatial distribution is sampled. Thus, for each IMF sampling algorithm, one would first need to draw the mass of the star, and then the surface density distribution in order to make a stellar population. At the moment, it is unclear how the surface or volume density distribution of young stars depends on environment (see Bressert et al. 2010), hence we cannot say what fraction of isolated stars would be expected, within this scenario.

6. Summary and implications

We have presented a new method to identify massive stars that formed in relative isolation. Applying the method to the 30 Doradus region in the LMC, 15 O-star candidates are found, where 11 are Vz stars, that may have formed in isolation. The method uses precise radial velocities of the stars to rule out massive binaries and runaways along the line of sight. Additionally, high resolution imaging is used to constrain the possible presence of a stellar cluster around the stars. Using extensive Monte Carlo stellar cluster simulations we confirm that our observations are sufficiently sensitive to rule out typical stellar clusters that should be associated with the candidates if the m_max − M_cl relation is correct. Finally, we search for gas/dust filaments that are associated with the massive stars (at least in projection), i.e., the birth-sites of stars, to mitigate the possibility that the stars are runaways in the plane of the sky. The gaseous filamentary structures were identified using a heterogeneous set of techniques, wavelengths and instruments. This could be significantly improved through the use of high-resolution mid-IR imaging of the cold dust in the region, such as that provided by the Herschel Space Telescope. The observations presented here can constrain theories of massive star-formation as well as scenarios for sampling from the stellar IMF.

All the candidates, except VFTS 682, have met the criteria presented in Sect. 3. We have mentioned several caveats, including filamentary structures and stellar age issues and possible bow shock association that could be problematic for grade 2 and 3 candidates. Hence, the best candidates of the 16 that withstand the caveats are the 11 Vz (ZAMS) stars, where their likely young ages, high mass, robustness to different stellar model parameters, and general close agreement between RV_star and RV_ISM is within 15 km   s^-1, with the exception of 682 (40 km   s^-1) for which formal uncertainties in stellar RV measurement are likely too low.

Isolated massive star-formation is possible according to the monolithic formation scenario (e.g., Krumholz et al. 2009), whereas the competitive accretion scenario implies that massive stars can only form in sufficiently massive clusters (e.g., Bonnell et al. 2004). The observations presented here suggest that competitive accretion might not be the only mechanism responsible for the formation of massive stars. Whether a massive-star forms in a cluster or not can be affected by an external agent, such as triggering, but once star-formation starts in the cluster radius of the massive star, external triggering is no longer an important factor (see Dale & Bonnell 2011).

We have found 16 candidates in the LMC that most likely formed in isolation or in sparse clusters that do not follow the m_max − M_cl as reported by Weidner et al. (2010). This supports the conclusions of Lamb et al. (2010) who studied O-stars in the SMC which had sparse clusters associated with them. Both studies attempt to remove the possibility that the O-stars investigated are runaways from a nearby cluster. Chu & Gruendl (2008) found that 4% of massive YSO candidates (with some contaminant, Gruendl & Chu 2009; Evans et al. 2011), in the LMC were located far enough away from young clusters and/or associations that they most likely formed in isolation. Additionally, lower in the stellar mass range, a collection of Ae/Be stars have recently been shown to not originate in clusters (Wheelwright et al. 2010). Finally, Eldridge (2012) notes that γ-Velorum is a  ~65 M_⊙ WR+O binary system in a cluster with a total mass of 250–350 M_⊙, which is well below that expected from the m_max − M_cl relation.

An additional constraint that not all of the surrounding O-stars in 30 Doradus are runaways, can be obtained by comparing R136 and the Galactic young massive cluster NGC 3603. The clusters have similar ages, masses, and densities, yet NGC 3603 does not appear to have any massive stars outside the central  ~1 pc (Moffat et al. 1994). Since the number of runaways scales with the central cluster density, NGC 3603 and R136 are expected to have similar numbers. Hence, the observed differences are likely a reflection of the formation of the cluster and surrounding environment.

Hence, the observations presented here (and those cited above) appear to contradict the m_max − M_cl relation shown in Weidner et al. (2010), in that large clusters are not necessary to form massive stars. This favours a scenario that the stellar IMF is sampled stochastically. Amongst the candidates, we have shown that short-period binaries are very rare. There may be unresolved lower-mass stars and a handful of long-term binaries that accompany the candidates, however, if present these stars wouldn’t be enough to raise the m_max − M_cl to fit the sorted sampling scenario. The fact that the m_max − M_cl relation does not appear to be causal, but rather statistical, raises questions regarding the Integrated Galactic IMF (IGIMF) scenario (Weidner & Kroupa 2006), where galaxies that form only low mass stellar clusters will never form high-mass stars. This has many potentially large implications for extragalactic observations (cf., Bastian et al. 2010), especially for low-mass galaxies.

The observation of the formation of high-mass stars in isolation, along with the age spread within the 30 Doradus region (e.g., Walborn & Blades 1997), is consistent with a hierarchical distribution of star-formation, in time and space (e.g. Efremov & Elmegreen 1998; Bastian et al. 2009). This means that star-formation does not happen in a quantised way on the group level, i.e., not all stars form in compact coherent structures like clusters. The spatial distribution of star-formation appears to be a continuous distribution in surface density (Bressert et al. 2010) of which the stars discussed here would fall into the extreme low surface density end of the distribution. Additionally, clusters are sensitive to how one exactly defines them, as the distribution of clusters and associations seem so overlap at young ages (Gieles & Portegies Zwart 2011).

Online material

Appendix A: Hubble archival data

We make use of HST images that were retrieved through the Hubble Legacy Archive web service. The images came from a collection proposals and programs. A table is provided below to summarise the different data and their properties which are important for this paper.

Appendix B: Probability test for chance alignment of runaway stars with filaments

Currently, there are no high-resolution, continuous, large field-of-view maps of the 30 Doradus region that show where the star forming filamentary gas and dust are. In the coming year this deficiency will be resolved as the Herschel Space Observatory survey maps of the LMC, HERschel Inventory of The Agents of Galaxy Evolution (HERITAGE) (Meixner et al. 2010), will be publicly released. By combining the HERITAGE maps with the data presented in this paper one could reject the likelihood of O-stars being runaways based on probability. We propose that the definition of a filament from Herschel observations is best defined by the 2nd differential maps Molinari et al. (2010) introduced, where a physically relevant radius from the peaks of the filaments can bear association to the stars.

To do this, we take O-stars similar to the ones discussed in this paper and assume that they are runaways. Using the radial velocity measurements from the Tarantula Survey, we

approximate that there are roughly N number of runaways in total. Then we place a circular boundary around R136 at distance R, where R = R_boundary − R_cluster, and has area A_boundary. Within R distance, let’s assume that gas/dust filaments covers an area A_filament, such that A_filament ≤ A_boundary. Refer to Fig. B.1 for a visual context of the scenario we have now discussed.

With the above conditions, what is the probability that n of N runaway O-stars line up with filamentary structures in the line-of-sight? For simplicity, assume that the N O-stars are randomly distributed within the boundary such that the probability of an O-star lining up with a filament can be expressed as ρ. $\mathit{\rho }\mathrm{=}\frac{{\mathit{A}}_{\mathrm{filament}}}{{\mathit{A}}_{\mathrm{boundary}}}\mathrm{·}$(B.1)To see what probability n of N (n ≤ N) O-stars is in the line-of-sight with filaments we invoke the binomial probability equation as shown below, where we assume that each line-of-sight event is independent of one another. $\mathit{P}\mathrm{\left(}\mathit{n}\mathrm{\right)}\mathrm{=}\left(\begin{array}{}\mathit{N}\\ \mathit{n}\end{array}\right){\mathit{\rho }}^{\mathit{n}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\rho }{\mathrm{)}}^{\mathit{N}\mathrm{-}\mathit{n}}\mathit{.}$(B.2)

---

Acknowledgments

We would like to thank Nolan Walborn, Sally Oey, John Bally, Leonardo Testi, and Adam Ginsburg for their comments and suggestions. Additionally, we thank Jasmina Lazendic-Galloway for sharing Hα/Hβ extinction map with us and Carsten Weidner for interesting discussions. Based on observations at the European Southern Observatory Very Large Telescope in programme 182.D-0222. Based on observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESA) and the Canadian Astronomy Data Centre (CADC/NRC/CSA).

References

All Tables

All Figures