Supersonic turbulence in Giant HII Regions: clues from 30 Doradus

The tight correlation between turbulence and luminosity in Giant HII Regions is not well understood. While the luminosity is due to the UV radiation from the massive stars in the ionizing clusters, it is not clear what powers the turbulence. Observations of the two prototypical Giant HII Regions in the local Universe, 30 Doradus and NGC604, show that part of the kinetic energy of the nebular gas comes from the combined stellar winds of the most massive stars - the cluster winds, but not all. We present a study of the kinematics of 30 Doradus based on archival VLT FLAMES/GIRAFFE data and new high resolution observations with HARPS. We find that the nebular structure and kinematics are shaped by a hot cluster wind and not by the stellar winds of individual stars. The cluster wind powers most of the turbulence of the nebular gas, with a small but significant contribution from the combined gravitational potential of stars and gas. We estimate the total mass of 30 Doradus and we argue that the region does not contain significant amounts of neutral (HI) gas, and that the giant molecular cloud 30Dor-10 that is close to the center of the nebula in projection is in fact an inflating cloud tens of parsecs away from R136, the core of the ionizing cluster. We rule out a Kolmogorov-like turbulent kinetic energy cascade as the source of supersonic turbulence in Giant HII Regions.


Introduction
Besides the obvious difference in size, the fundamental distinction between Giant HII Regions (GHRs) and their smaller (mostly) Galactic counterparts (e.g. Orion) lies in the velocity width of their integrated emission-line profiles. Whereas small HII regions have thermal velocity widths, the integrated widths of GHRs are highly supersonic.
Yet, despite decades of research and discussions, there is still no agreement as to the origin of the supersonic turbulence in GHRs. As so many other human endeavors of our days, opinions are divided into two seemingly irreconcilable camps: the "gravity" camp, championed by the authors of this paper in various combinations many with other collaborators (see e.g. Fernández Arenas et al. 2018 for a recent discussion and historical perspective); and the "stellar-winds" camp led by Y.H. Chu and R. Kennicutt (Chu & Kennicutt 1994).
The defenders of gravity argue that the surprisingly tight correlation between the global Balmer line-profile widths (σ; Hα or Hβ) and the integrated flux of these lines -the so-called L-σ relation -(shown in Figure 1) is the unequivocal telltale of gravity, although the slope of the correlation (L ∝ σ 5 ) cannot be deducted from first principles. The fact that HII region-like galaxies, or HII Galaxies also define a tight L − σ relation of similar slope, is also an argument for the gravity camp.
The stellar winds camp, on the other side, base their arguments on the morphology of GHRs, which is dominated Send offprint requests to: Jorge Melnick e-mail: jmelnick@eso.org Partly based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programmes 072.C-0348; 182.D-0222; and 60.A-9700(G), and programme ID 076.C-0888, processed and released by the ESO VOS/ADP group. by prominent arcs and filaments, and on high-resolution spectroscopy that shows that the emission lines associated with these morphological features are double or multiple as expected for expanding bubbles of ionized gas propelled by stellar winds. Chu & Kennicutt (1994) showed that when the ragged shapes of the spatially resolved profiles are added-up, the resulting integrated profiles have a smooth supersonic Gaussian core with broad unresolved wings matching surprisingly well observations of unresolved nebulae. The idea that winds from individual stars can generate the observed bubbles has been disputed by Silich et al. (2020) who argue that the close proximity of the massive stars within the ionizing clusters of GHRs implies that their winds must necessarily merge into a hot cluster wind that percolates the pre-existing interstellar medium generating thus the hierarchy of expanding shells that we see today. Although the difference between stellarwinds and a cluster-wind is fundamental, it does not change the original argument of Chu & Kennicutt (1994) to explain the multiple shell structure of GHRs.
Much of the work of the "stellar-wind" camp is based on observations of the two nearest and therefore best studied GHRs in the Local Group: 30 Doradus in the LMC and NGC604 in M33, but specially of 30 Doradus, which has become a sort of Rosetta Stone for deciphering the astrophysics of giant HII regions and their ionizing stars. In fact, the ionizing cluster of 30 Doradus (sometimes referred to as NGC2070; others as the R136 cluster), contains the largest concentration of resolved massive stars in the local Universe, including some of the few supermassive stars known (see Castro et al. 2018 for a recent review). Chu & Kennicutt (1994) argued that about 50% of the kinetic energy of the nebular gas in 30Dor is in the expanding shells and 50% in some unspecified form of general turbulence. The mass of stars+gas in the nebula, (M ∼ 5 × 10 5 M ) they argued, is not enough to generate the observed turbulence through the gravitational potential.
The aim of the present paper is to use new public FLAMES/GIRAFFE data to investigate in (renewed) detail the nature of the supersonic turbulence of the nebular gas in 30 Doradus.

The kinematics of the nebular gas in 30 Doradus
The large angular size of 30 Doradus (henceforth 30Dor) has made it difficult to fully map the kinematics of the nebular gas. Smith & Weedman (1972) used a single channel Fabry-Perot (FP) interferometer to observe a number of positions in 30Dor with spatial resolutions of 13 ; 30"; 120 . Melnick (1977) (reproduced in Melnick et al. 1987) followed-up on these observations and showed that the structure function of turbulence in 30Dor was flat, albeit with a spatial resolution of only 13 = 3.25pc . Long before that, however, Feast (1961) had already found that there was very little structure in the turbulence of the nebula. Chu & Kennicutt (1994) improved upon these previous studies by dissecting 30Dor using long-slit Echelle spectrograms that provided a significantly better spatial coverage, particularly in the outer regions of the nebula. They found that the velocity field of the nebula is dominated by expanding shells of a wide range of diameters, and concluded that the supersonic velocitywidth of the integrated Hα profile was by and large the result of these expansion motions, although the shells accounted for only about 50% of the kinetic energy of the gas. They also dismissed the contribution of gravity in stirring-up the gas kinematics. Melnick et al. (1999) refined the observations of Chu & Kennicutt (1994) using higher spatial and spectral resolutions, albeit at only one slit position. They confirmed that the supersonic Hα profiles are resolved into multiple components, but at all positions along the slit they found a very broad unresolved component of unknown origin. More recently, Torres-Flores et al. (2013) used the multiplexing capabilty of FLAMES/GIRAFFE on the VLT to address some of the issues raised by Chu & Kennicutt (1994) and Melnick et al. (1999). While they largely confirmed the findings of these previous studies, they questioned the reality of the broad unresolved component identified by Melnick et al. (1999).
For the ensuing discussion it is important to recall that several phases of the interstellar medium coexist within 30 Doradus: It seems to contain significant amounts of molecular and atomic Hydrogen (Ochsendorf et al. 2017 and references therein); there is the ionized gas at 10 4 K, which we will call the nebular gas; and there is a very hot X-ray emitting medium at 10 7 K that we will call the hot-gas. Here, when we talk about supersonic motions we are referring to the gas at 10 4 K and a sound speed of c s = 9.1 km s −1 for hydrogen. -Flores et al. (2013) raised doubts about the reality of the broad component of the integrated Hα profile of 30Dor identified by Melnick et al. (1999) , so it seems relevant to begin our analysis by revisiting this issue in order to set the ensuing discussion on self-consistent grounds. We used a subset of the data discussed by Torres-Flores et al. (2013) that we downloaded pipeline-processed from the ESO Archive (Phase3). We discarded from this dataset all points showing discernible contamination of the nebular profiles by stellar features, leaving a total of 1668 individual spectra in our full sample. Figure 2 shows the integrated profiles determined from these data in two different ways. The top panel shows the average of the 1668 individual FLAMES/GIRAFFE spectra. Torres-Flores et al. (2013) (see also our Figure 7) provide graphical illustrations of how the data points are spatially distributed in the nebula showing that the central regions are more densely sampled than the outskirts, mostly by the VLT Flames Tarantula Survey observations (Evans et al. 2011). Thus, in order to check for effects that this inhomogeneous sampling may introduce, the lower panel shows the integrated profile using only the data from the uniform grid of Torres-Flores et al. (2013) for a total of 879 points.

Torres
Our Gaussian fits for both the Full-Sample and the Uniform-Grid profiles are surprisingly similar given the differences in the two samples. The widths of the profiles also agree with those of Melnick et al. (1999) and Torres-Flores et al. (2013) for the core of the line, but the broad component of Torres-Flores et al. (2013), σ 2 = 49.1 km s −1 , is significantly broader than ours (σ 2 = 41.3 km s −1 ; lower panel) basically for the same data. The discrepancy is much larger than the formal fitting errors (< 1 km s −1 ), but it is well known that the parameters of multi-Gaussian fits are rather sensitive to the statistical errors of the data, which probably explains the discrepancy. We used the errors provided by the ESO Phase3 pipeline propagated through the various manipulations of the data, but we ignore how Torres-Flores et al. (2013) estimated the corresponding errors. Melnick et al. (1999) fitted two Gaussians to the total profile of Chu & Kennicutt (1994) and found σ 1 = 26 km s −1 and σ 2 = 44 km s −1 , without including errors. Notice that two Gaussians do not fit the center of the line for which a third component would be required.

Radial bins
The outer regions of the nebula have been characterized by Chu & Kennicutt (1994), so here we will concentrate in the inner (r < 40pc) part of 30Dor that contains the bulk of the ionizing Top. Integrated Hα profile of 30 Doradus using all the available data (1668 points). Two Gaussian components are required to fit the core and the wings of the line shown by the thin red lines. The data are shown as blue points and the sum of the two Gaussians by the red dashed lines. The legend reports the measured (left) and corrected for instrumental and thermal broadening (right) widths (σ) of each component plus the width from a single Gaussian fit (σ 0 shown by the black lines). Bottom. Same as top, but using only 879 points from the uniform grid of Torres-Flores et al. (2013). The heliocentric radial velocity of 30Dor obtained from the full-sample is V 30Dor = 264.8 km s −1 for an Hα rest wavelength of 656.2804nm. stars and has a complex structure that was not well covered in the study of Chu & Kennicutt (1994).
To begin, we present in Figure 3 profiles integrated over 7 radial rings of radii indicated in the legend. For each ring, the figure shows multiple Gaussian fits with the minimum number of components required to reproduce the observations. We will insist throughout this work that multiple-Gaussian fits are seldom unique being increasingly sensitive to data errors and to the initial guesses as the number of Gaussian components increase, so we must be conservative when interpreting the results from these fits. Nevertheless, Fig. 3 shows several robust features relevant to the present study. One is that the innermost ring (r < 10pc) shows two intense peaks that tell us that the kinematics of the very center of the nebula is dominated by a massive wind that is blowing away the gas at high speeds. Despite this raging wind, however, there is still plenty of gas at the center of 30Dor as indicated by the intensity scale.
A second feature is that there is a dramatic change in the kinematics at r = 25pc. Beyond this radius, the profiles become broader and increasingly show the emergence of a second redshifted component. Thus, beyond 25pc we see the onset of the "large expanding-shells" regime. A third feature is that the surface brightness drops dramatically as the structure of the gas becomes increasingly dominated by the large shells described by Chu & Kennicutt (1994). . Integrated profiles of the core of 30Dor (r < 40pc) and 7 radial bins identified in the figure legend that quotes the inner and outer radii of each ring in pc. The legend also gives in the right side widths of the multi-Gaussian fits of each ring and on the left, denoted as 0, the width and radial velocity of the single Gaussian fit (that is not plotted) relative to the integrated velocity of the nebular core that defines our rest-frame for 30Dor.

Profile taxonomy
The atlas of gray-scale spectrograms of Chu & Kennicutt (1994) shows that virtually everywhere in the nebula the Hα profiles tend to break-up into multiple components. However, the high contrast of the atlas prints makes it difficult to precisely discern the profile shapes. Thus, we were surprised to discover, upon visual examination of our 1668 Hα profiles in the computer screen, that the majority of the profiles can actually be divided into three well defined groups illustrated schematically by Figure 4.
The vast majority of the profiles are either single or double, although multiple profiles also contribute a fair fraction. We emphasize that we here are describing the Hα profiles. We will discuss the forbidden lines in the next section. Like the integrated profile, single profiles have Gaussian cores and broad non-Gaussian wings. Double and multiple profiles are characterized for having two or more smooth components of comparable (within factors of a few) peak intensities that are well separated in radial velocity. The radial velocity of the stronger component generally appears to be closer to zero, whereas the weaker peak is either redshifted or blue shifted relative to the stronger component.
We mentioned at the beginning that Chu & Kennicutt (1994) concluded that 50% of the kinetic energy of the gas is in the expanding shells and 50% in some form of unspecified turbulence. Assuming that double and multiple represent pure expansion, the single profiles must convey information about that "unspecified" form of turbulence. Therefore, it is important to quantify the distinction between single and multiple profiles.
We used an algorithm based on first and second derivatives to detect multiple peaks using the S /N of the continuum to weedout spurious peaks. Most of the profiles show some degree of low-level asymmetry, which does not make them multiple in the sense of having more than two peaks of comparable intensities. Therefore we introduced a finer distinction in order to separate single from multiple profiles. We defined as "double" profiles where the ratio of peak-intensities between the first (strongest) and second components P 2 /P 1 > 0.1, and "multiple" profiles having three or more components with P 3 /P 1 > 0.15.
The peak finding algorithm yields 1067 singles; 453 doubles and 148 multiples, confirming our initial impression that multiple profiles are less abundant but still numerous. Visual examination revealed that a good fraction of singles were actually doubles where the radial velocities of the peaks are too close to be detected by our peak finding algorithm.
Since the Hα profiles are blurred by thermal motions that hide components with a small difference in radial velocity, we refined our separation of single and multiple profiles using the lines of [NII] and [SII] that have much smaller thermal broadenings (2.4 and 1.6 km s −1 respectively). This "solution", however, turned out to be complicated as discussed below.

[NII] and [SII] kinematics
Using the [NII] and [SII] lines (henceforth the forbidden lines) to characterize the nebular gas kinematics has two problems. The first most obvious one is that the lines are much weaker than Hα and often their S/N is too low for a reliable characterization of the profiles. The second, less obvious problem, is that in general the peak intensities of the multiple components of Hα differ from those of the forbidden lines, as in the example shown in Figure 5.
It is not appropriate, therefore, to use the morphology of the forbidden lines to characterize the kinematics of Hα (and viceversa). So we only used the forbidden lines to flag single and multiple profiles, and we searched for multiple components in all profiles thus flagged by subtracting single Gaussian fits from the observed Hα profiles and fitting multiple Gaussians, using the positions of the residual peaks as initial guesses, to measure the peak intensities.
In this way we found that 599 of the 1067 lines classified as "single" in Hα by the peak finding algorithm were double in [NII] and/or [SII] and confirmed to have P 2 /P 1 > 0.1 at Hα , while 30 profiles previously classified as double were confirmed to be multiple (P 3 /P 1 > 0.15).
Relevant statistics from the three families of Hα profiles are presented in Table 2 that lists average radial velocities of the strongest components < V 1 >; the average peak intensity ratios; and the average velocity-widths from single-Gaussian fits (uncorrected for instrumental and thermal broadening; σ ). The ± values give the rms dispersions from the means. Using the partition of single and double profiles we can now quantify our second visual impression that the stronger component is generally at rest in the frame of the nebula while the weaker component is either redshifted or blue shifted relative to the strong component. Figure 6 shows the histograms of radial velocities of the two components of our "double" profiles. As in Table 1, P 1 corresponds to the strongest peak and P 2 to the weaker one.
It should be obvious from these distributions and the data of Table 1 that the integrated profile of 30Dor can be deconstructed into a sum of Gaussians representing the components of the double peaks plus the single peaks and their rms radial velocity dispersion.
Finally, let's look at the spatial distribution of our three families of profiles. This is presented in Figure 7 where we have color coded the positions according to morphology. Doubles tend to concentrate inside the large expanding shells; singles tend to . Distribution of radial velocities relative to the rest frame of 30Dor of the two components of the double profiles. P 1 (magenta) represents the component with highest peak intensity and P 2 (green) the weaker component.
lineate the high surface brightness borders of these shells; while multiples share the distribution of doubles.
The image illustrates that the profiles inside the large Chu-Kennicutt bubbles are double or multiple, while at the edges of these giant shells the profiles are mostly single. We have color coded the single profiles according to their peak radial velocities relative to the systemic velocity of 30Dor (265.4km s −1 ) from our Hα data. The edge of Chu-Kennicutt#2, that has been very well mapped in search for photodissociation regions, is seen to be blue-shifted indicating that the shell is expanding towards us, while the inside shell that was identified by Chu & Kennicutt (1994) but not numbered is completely delineated by red dots.
The single profiles, therefore, delineate the edges of the expanding shells as well as the dense filaments and loops that characterize the core of the nebula, and which will be described in detail below. Notice that the blue and the red dots are rarely mixed showing that the densest parts of the nebula as a whole are expanding.

Single Profile Kinematics
Assuming that the kinematics of the nebular gas represented by double and multiple profiles is dominated by expanding motions, single profiles should provide us with information about the 50% of the total kinetic energy of unknown forms of turbulence mentioned by Chu & Kennicutt (1994). We have just seen that the radial velocities of the single profiles are dominated by expansion motions, but the individual profiles also have supersonic widths. Table 2 summarizes the average intrinsic width (i.e. corrected for instrumental and thermal broadening) and heliocentric radial velocity of the single profiles for Hα and for the forbidden lines.
The table shows that the forbidden lines are on average ∼ 1.3 km s −1 narrower than Hα. Figure 8 shows the stacking of the single [NII] and Hα profiles after rebinning the data to z = 0. We have repeatedly whined about the dangers of multi-Gaussian fits, but in the case of Fig. 8 a minimum of three components is required to fit the core and the wings of the lines. In fact there are good physical reasons to fit three components; we have already encountered two: σ 1 , the main peak, and σ 3 , the broad unresolved component that pervades the nebula. The third component (σ 2 ) corresponds to the sum of all the low-intensity components that do not meet our P 2 /P 1 > 0.1 and/or P 3 /P 1 > 0.15 conditions.
A second interesting feature of Fig. 8 is that, as already shown in Table 2, the Hα profile is broader than [NII] by 1.6 km s −1 . The difference is somewhat larger than the value of 1.3 km s −1 from the table that assigns the same weight to all profiles independently of peak intensity, so the agreement is actually quite good.
We verified that the systematic difference in the line-widths is not due to a geometrical effect (see Appendix A). In fact, Hippelein (1986) found a similar systematic difference between [OIII]500.7nm and Hα for the integrated profiles of a small sample of Giant HII Regions and HII Galaxies. This difference has been confirmed by several authors notably by Melnick et al. (2017) for a substantial sample of HII Galaxies. A similar effect was also observed in the Orion Nebula by García-Díaz et al. (2008) but there is no consensus in the literature as to the origin of the discrepancy. Notice, however, that while the intrinsic width of Hα (σ 1 ) is the same as the sound speed for hydrogen at 10 4 K, the [NII] profiles are highly supersonic, suggesting that perhaps they could be resolved into even more components.
We also measured the intrinsic widths of each of the individual components of double Hα profiles -taken separatelyusing two Gaussians to fit the doublet plus a third Gaussian to fit the extended wings, using an automated procedure. We find an average width of 11.7 ± 4.7 km s −1 for these components, comparable to the values given above for the average width of single profiles (11.9 ± 3.4km s −1 ).

The [OIII] profiles with HARPS
We speculated above that a possible explanation for the reason the intrinsic width of single [NII] profiles is highly supersonic while Hα is not, could be that at higher resolution the profiles would be resolved into multiple narrower components. We therefore observed two characteristic regions of 30Dor with HARPS: one at the very center (R136), and the second at the position of Filament F from the NTT observations of Melnick et al. (1999), located 8pc due west of R136. The original purpose of this pilotproject was precisely to resolve the profiles further, and to verify if the forbidden [OIII] lines also showed the broad unresolved component seen by Melnick et al. (1999) at Hα.
HARPS provides a spectral resolution of R=100000 covering the spectral range between 370nm and 690nm with a spatial resolution, determined by the entrance fiber, of 3 . Figure 9 shows the HARPS spectrum of the center of 30Dor (R136) that nicely illustrates the extremely broad "feet" of the Balmer lines from the massive WN7h stars in R136 that provide a significant Article number, page 5 of 15 A&A proofs: manuscript no. girafR30    The Hα profiles shown in the bottom panels of the figure are double in the sense discussed above: they can be fitted reasonably well using two Gaussian components for the two main peaks plus a third for the broad extended wings. The [OIII] lines, on the other hand, have higher S/N and much smaller ther- mal broadening (2.3 vs. 9.1 km s −1 at 10 4 K), allowing to discern a significantly more complex structure at the same positions, which requires at least 7 components, plus one for the broad wings, at the center, and at least 5 components plus the broad wings for Filament F. The corrected velocity-widths of the individual components are given in the figure legends.
Although the [OIII] line resolves into several components, these are still supersonic at 10 4 K, but again we must return to the "mantra" of this paper: beware of multi-Gaussian fits. The solutions shown in Figure 10 result from leaving all parameters free. It is actually possible obtain solutions by forcing the widths of all components to be narrow, say, σ = 3 km s −1 , albeit only after carefully fine-tuning the initial guesses of all other parame-ters. So we can neither prove nor rule out that the [OIII] profiles are the sum of many subsonic components, for which we would need much higher S/N data.
Nevertheless, our results clearly show that the individual components of double profiles are resolved into a number of narrower components indicating that there is structure in the expanding shells at the spatial resolution of HARPS (0.75pc). The radial-velocity dispersion of these components is "typically" (we only have 4 examples) ∼ 11 km s −1 , quite similar to the value measured from the unresolved shells and filaments.

The kinematical core of 30 Doradus
3.1. Structure Figure 3 shows that there are two inflection points in the integrated profiles of the central part of 30Dor: the region r < 10pc, where the profile is double and the kinematics dominated by the massive cluster wind emerging from the stars in the cluster core (R136); and the region r > 25pc where expanding motions increasingly dominate the profile shapes. We used finer bins to verify that the latter change is already visible between 25pc and 27pc.
The radial Hα brightness distribution, shown in Figure 11, also shows features (i.e. changes of slope) at r ∼ 10pc and r = 25pc amid a broad peak that encompasses 33% of the total Hα flux of the nebula and the bulk of the ionizing stars, and where the cluster wind that drives the hierarchy of Chu-Kennicutt shells originates.  Figure 3 shows a distinct inflection point. Figure 12 shows that the width of the integrated profile of the region r ≤ 25pc (henceforth the kinematical core) is virtually identical to that of the nebula as a whole (the only significant difference being the width of the broad component). Therefore, the ratio of kinetic energies of the nebula as a whole and the kinematical core is just the ratio of the Hα luminosities; the core contains 23% of the total kinetic energy of the nebular gas in 30Dor. Perhaps tellingly, the shape of the profile of the kinematical halo (r > 25pc; the realm of the large expanding bubbles) is virtually identical to those of the core and of the overall nebula.
Article number, page 7 of 15 A&A proofs: manuscript no. girafR30 Fig. 12. Comparison between the integrated Hα profiles of the whole nebula (bottom) with the profile of the kinematical core (r ≤ 25pc), and the kinematical halo (r > 25pc). The blue points are the data and thickdashed red line is the fit of two Gaussian components shown in thin red. The single-Gaussian fit is shown by the black line. The parameters of these components, observed and corrected for instrumental and thermal broadening are shown in the legend. Figure 13 presents a zoom to the central 30pc of 30Dor seen in X-rays (right) and Hα (left). Overlayed on these images are colored dots coded according to profile shapes and kinematics.
Single profiles delineate bright (dense) Hα filaments, and lie at the edges of bright X-ray zones. These dense filaments define the "head" of the Tarantula, which, as shown by the distribution of blue-shifted and red-shifted profiles (blue and red dots), consists of several expanding bubbles that are the mouth and the eyes of the spider. The large elliptical bubble that surrounds R136 is asymmetrical with the eastern side being more luminous in Xrays. This "mouth of the Tarantula" opens towards the east indicating that the hot gas is escaping through the mouth to inflate the giant bubble Chu-Kennicutt#2, as one would inflate a balloon for a birthday party! There is another expanding bubble due North of R136 that is very bright in X-rays and that is open to the north toward the giant Chu-Kennicutt#5 shell. Similarly the bubble due south of R136 seems to be inflating Chu-Kennicutt#1. So we get the impression that the dense filaments that define the head of the Tarantula act as nozzles that channel the cluster wind from the source (R136) out to the giant Chu-Kennicutt shells.
The clumpy structure of the hot cluster wind seen in the Xray image is mostly due to shadowing by the dense dust+gas structures that can be clearly appreciated in the Hα image. The X-rays are faint in places where the Hα emission is bright and blue-shifted, indicating that the nebular component is in front of the hot gas, and viceversa. Receding filaments seem to lie behind that hot gas, but not always; the bright loop directly NW of R136, where the X-rays are faint, is redshifted. However, the giant shell Chu-Kennicutt#3 is adjacent so it is possible that the hot wind has already escaped from that region.

Kinematics
To compare with the results of Chu & Kennicutt (1994), Table 3 presents an estimation of the kinetic energy budget of the nebular gas. The total kinetic energy is 2K tot = M gas σ 2 0 where σ 0 is the integrated velocity dispersion of the gas estimated from the single-Gaussian fit of figure 2 and M gas = 4 × 10 5 M is the total mass of nebular gas according to Chu & Kennicutt (1994). The values for the kinematical core and halo are obtained assuming that the mass of gas is proportional to the Hα luminosity.
The kinetic energy of the "general" turbulence was calculated as follows. We saw that the single Hα profiles as well as each of the two components of double Hα profiles have average widths of ∼ 11.8 km s −1 . We do not know the source of energy that sustains these motions, but we do know that when we observe the same profiles in the forbidden lines, they are resolved into many narrower components. We also ignore the mass of these turbulent shells and filaments, but, since they comprise the brightest regions of 30Dor, it is reasonable to assume that they contain most of the nebular mass (assuming that the mass is proportional to the Hα luminosity) . We calculated the velocity dispersion in two ways: directly from the average profile widths (min) and convolving the profile width with the radial-velocity dispersion of the single profiles (max), σ max = √ 11.8 2 + 11.6 2 = 16.5 km s −1 . Since at least some, and probably much, of the radialvelocity dispersion of the single profiles is due to expansion motions, the velocity dispersion of 16.5 km s −1 given in the table is a maximum. Thus, the fraction of the total kinetic energy contained by general turbulence is closer to 25% than to 45%, compared to 50% estimated by Chu & Kennicutt (1994).

Kolmogorov turbulence
Observations of the other prototypical GHR in the local group, NGC604, by Medina-Tanco et al. (1997) appear to show a well developed Kolmogorov-like kinetic energy cascade in that giant nebula. They found that the structure function, -the square of the difference in radial velocity between two points separated by a distance r averaged over all pairs of positions separated by r -is a power-law of slope ∼ 5/3 for r ≤ 10pc and flattens at larger separations.
Since by definition σ ≡ v 2 − v 2 , if v = 0 the structure function S 2 = 2σ 2 , so it is convenient to normalize by the square of the the radial-velocity dispersion σ rv such that the normalized function peaks at S 2 ∼ 2 (Arthur et al. 2016). Figure 14  reproduces the normalized structure function of NGC604 calculated using the same radial velocity data of Medina-Tanco et al. (1997) kindly provided by Gustavo Medina-Tanco.  Medina-Tanco et al. (1997). The red bar shows the radius (4.1pc) corresponding to the seeing disk (1 ) of the TAURUS-II observations used to construct the structure function.
We used our full sample to construct the structure function for 30Dor as for NGC604, but using three different ways of measuring the radial velocities. The top panel of Fig.15 corresponds to the method adopted by Medina-Tanco et al. (1997) of fitting single Gaussians to all profiles. The second, shown in the middle panel, uses the radial velocity of the strongest peak at each position; and the third uses only the single profiles. Typically each point in the top two panels is the average of 16000 to 78000 pairs, except for the smallest separation that has only 6500 pairs. The vertical bars represent the rms dispersion around each average. Contrary to that of NGC604, the structure function of 30Dor has no structure: it is flat! Since the L − σ relation of Fig. 1 implies that the kinematics of GHRs depends only on their luminosities, we expected 30Dor and NGC604 to have similar, if not identical, structure functions, which turned out not the case. Why?
Unfortunately, there are several problems and inconsistencies with the results of Medina-Tanco et al. (1997). The first and most obvious is seeing. The red bar in Fig. 14 at 4.1pc corresponds to the seeing (1 ) reported for the Taurus-II observations of Sabalisck et al. (1995) upon which the structure function is based. The separation of the first two points in the structure function is less than the seeing, while the next two points at < 2 are also affected by seeing. These four points define the slope of the function for r < 10pc. Removing these four seeing affected points the remainder of the structure function of NGC604 is basically flat.
A second issue that may blurr the comparison between 30Dor and NGC604 is beam-smearing. At the distance of M33, 1 = 4.1pc" while the data of Medina-Tanco et al. (1997) were binned to 2 × 2 , which results in a significant smoothing of the radial velocities. There are also some problems with the TAURUS-II data for NGC604. For example, the TAURUS-II profiles of Muñoz-Tuñón et al. (1996) are substantially broader than the TAURUS-I (plus Echelle) results of Yang et al. (1996). We were unable to recover the original data to check these issues, so for the moment it seems safer to ignore the results of Medina-Tanco et al. (1997), although we cannot overemphasize the importance of re-observing NGC604 hopefully at a higher spatial resolution.

Gravity
Chu & Kennicutt (1994) dismissed gravity as an important source of kinetic energy because "...even the largest imaginable mass for 30Dor would be far too small to explain the velocity dispersion of the gas inferred from the width of the integrated profile", σ 0 = 23.5 km s −1 (Fig. 12). But we know, as they knew, that a major part of this width is furnished by expanding shells. So let us revisit the issue.
We calculated the structure function of an N-body model of a 10 3 M virialized cluster containing 2000 particles of equal mass (kindly provided by Sverre Aarseth) shown in Figure 16, Fig. 16. Structure function of a virialized cluster of 2000 particles. All particles have equal masses and the total mass of the cluster is 10 3 M . For a virialized system S 2 (r) 2σ 2 and the velocity dispersion σ is constant in the core.
The match is not perfect, but still provides a much better match than a Kolmogorov turbulent energy cascade. This does not prove that gravity is driving the turbulence in 30Dor, but at least provides a consistency check.
Another way of probing the effect of gravity is to estimate the mass required to generate the observed motions through the virial theorem. The morphology of 30Dor and the results of Chu & Kennicutt (1994) show that the outer parts of the nebula are dominated by expanding motions, so here we will restrict our analysis to the kinematical core. Figure 17 shows the radial-velocity dispersion σ rv (in red) and the mean velocity width (in blue) of the single profiles plotted as a function of radius. Each bin contains 50-70 points except for the first bin (10 < r ≤ 25pc) that contains 98 single profiles. We excluded the central 10pc that is completely dominated by expansion motions. The average width of the profiles in the kinematical core is < σ core >= 11.2 ± 3.0 km s −1 and the radial-velocity dispersion is σ rv = 15.3 km s −1 . The respective minimal and maximal turbulent velocities in the core, therefore, are σ min = 11.2 km s −1 and σ max = √ 11.2 2 + 15.3 2 = 19 km s −1 . Assuming an effective radius of 10pc (instead of the 100pc used by Chu & Kennicutt 1994), the corresponding mass limits from the virial theorem are, 2.9 × 10 5 ≤ M/M ≤ 8.4 × 10 5 (2) The observed (photometric) mass of the kinematical core is quite uncertain, the stellar mass.being probably the least uncertain. Selman et al. (1999) quote a mass of M * = 68000M that results from integrating the observed (power-law) IMF between 3M and 120M . The low-mass (m<3M ) IMF of 30Dor is not known, but the IMF of most young star forming regions turns over at masses ∼ 0.5M (Bastian et al. 2010). Extrapolating the IMF of Selman et al. (1999) down to 0.5M , yields a total stellar mass of M = 1.5 × 10 5 M for 30Dor, 50% larger than the value generally quoted in the literature (e.g. Doran et al. 2013 and references therein), The mass of ionized gas is more uncertain. Here we will adopt the value of M HII = 4.5 × 10 5 M from Chu & Kennicutt (1994), which according to these authors is uncertain by factors of 2-4. Ochsendorf et al. (2017) derive a similar mass with an uncertainty of "only" a factor of 2 (0.2dex).
The most uncertain component is the mass of neutral and molecular gas also discussed by Ochsendorf et al. (2017). A comprehensive review of this contribution is well beyond the scope of the present paper, but in Appendix B we present a summary and justify our optimistic value of M HI+H 2 = 9 × 10 4 M for the amount of neutral and molecular gas within the nebular core.
The core of 30Dor contains the bulk of the stars and 23% of the Hα luminosity, so the total mass within the kinematical core is, which is comfortably above the minimum value required by Eq.2. So, again, we may have consistency, but not proof. A simple and direct proof would be to show that the velocity dispersion of the stars is consistent with that of the gas. Unfortunately, however, this test turns out to be neither simple nor direct. Bosch et al. (2001) and Bosch et al. (2009) found a surprisingly large velocity dispersion for NGC2070 that can be explained if a large fraction of the massive stars in the cluster are binaries. Subsequently, binary orbits have been measured for many of these stars confirming that conclusion (Hénault-Brunet et al. 2012a). There are several other complications related to mass segregation and rotation (Hénault-Brunet et al. 2012b) that compromise the putative simplicity of the test, but the observations show that, at least in the cluster core (r < 5pc), the stars are in virial equilibrium (Hénault-Brunet et al. 2012a,b).
To conclude this discussion, let us revisit the mass required to explain the integrated velocity dispersion by gravity. Adding up the HII mass (4.5 × 10 5 ), the stellar mass (1.5 × 10 5 ), and the molecular mass (1.8 × 10 5 ) the total mass of the nebula is M 30Dor = 7.8 × 10 5 M . The half-mass radius, therefore, is close to the radius of the kinematical core (r = 25pc), but much smaller than the value of 100pc assumed by Chu & Kennicutt (1994). The viral mass for a velocity dispersion of 23.5 km s −1 and an effective radius of 25pc is 3.2 × 10 6 M , still significantly larger than our optimistic value for M 30Dor , but by no means "the largest imaginable"!

Small scale turbulence
We have seen that there is no large scale supersonic hydrodynamical turbulence in 30 Doradus, and that the gravitational potential is not strong enough to furnish all the kinetic energy of the gas: The large-scale structure of 30Dor is organized by the cluster wind, but the wind-driven shells are themselves turbulent.
This turbulence could be due to Kolmogorov-like kinetic energy cascades powered by the wind, such as observed in the Orion Nebula (Arthur et al. 2016). However, even within the shells and filaments of 30Dor, we are looking spatial scales orders of magnitude larger than the turbulent regions in Orion, so we have no hope of resolving the turbulent cascade. However, our data (mostly HARPS but also some FLAMES profiles) show that the forbidden profiles can be resolved into many components with radial-velocity dispersions similar to the widths of the profiles. Such mutiplicity would explain why at the same positions the Hα profiles are subsonic but highly supersonic in the forbidden lines.
An intriguing possibility to explain the small scale turbulence is that sections of the leading shocks composing the large structures loose energy by radiation as the medium they sweep piles-up in front, so necessarily secondary shocks emanate from the hot gas to transmit the interior pressure. These "shockletts" slam first onto swept-up gas and eventually onto the leading shocks giving them small jolts that here we relate to the structure we observe in the forbidden lines.
The ionization of the 30Dor nebula is by and large provided by the massive stars, with a small contribution from fast shocks that, in principle, should be more important in the denser parts of the nebula (Pellegrini et al. 2011). Figure 18 shows the shock sensitive diagnostic diagram that can be obtained using the FLAMES spectra, color coded according to the distance to R136. Fig. 18. Point-to-point nebular diagnostic diagram for the different types of profiles described in the text. The points for the single profiles are colored according to the distance to R136 shown by the scale on the right. The solid curve shows the low-velocity shock models of Dopita & Sutherland (1996) for a range of velocities between 150 km s −1 and 500 km s −1 ; two values of the magnetic field density (0 and 1 µGcm −3 ); and solar abundances. The dashed curve shows the same models scaled to the metallicity of 30 Dor from Pellegrini et al. (2011). The figure shows that emission from low velocity shocks is clearly visible in the low-density external regions of the nebula, but not in the densest regions in the core of the nebula. This does not mean that shocks are not present in this region; it only means that, if present, they are completely overshadowed by photoionization at optical wavelenghts, reflecting, perhaps, how the cooling process operates behind the leading shocks as these continue to ram into a cloudy medium of very different densities. This also determines the speed of the shockletts driven into the nebular swept-up gas as this cools and approaches the nebular equilibrium temperature.

The broad unresolved component
Our analysis thus far has ignored the mysterious broad unresolved component that pervades the nebula and that is seen both in the Balmer lines and in the forbidden lines of [NII], [SII], and [OIII], precisely because we ignore its origin. We have seen, however, that like the narrow components, the broad unresolved component of the forbidden lines is narrower than Hα, so it seems reasonable to conjecture that the physical conditions of the gas that emits the broad unresolved lines are similar to those of the nebula as a whole, and therefore, that it is probably due to a large collection of fragments of shells and filaments left behind by the cluster-wind as it successively breaks out of shells of swept-up material.
We have not explicitly included this component in our estimates of the kinetic energy budget of the nebular gas, for which we used single-Gaussian fits. Although the width of the broad component is large, its residual Hα flux relative to the single Gaussian fit is negligible, so we can safely ignore the contribution of the broad component to the kinetic energy budget.

Discussion
There is broad disagreement in the literature about the structure of 30Dor and its genesis (Pellegrini et al. 2011;Lopez et al. 2011). Below we present a summary of our own view, which is the result of a series of in-depth investigations about the formation and evolution of super-star clusters and their associated Giant HII Regions (Silich et al. 2020 and references therein), and which we will use in the interpretation of the observations presented in this paper.

The genesis of Giant HII Regions
The progenitors of GHR are Giant Molecular Clouds (GMCs) that are clumpy and highly turbulent as a consequence of multiple prior generations of star formation, the gravitational field of the parent galaxies, and their own gravitational potential, which is substantial given their large masses (e.g. Heyer & Dame 2015;Krumholz et al. 2019). Even before any stars form, the turbulence within these GMCs is supersonic (with respect to the sound speed of the molecular gas (∼ 1 km s −1 ), causing the condensations to accrete more molecular material as it cools behind their leading shocks. This sets the initial conditions at the onset of star formation that occurs within the largest and most massive condensations.
The arrival of the most massive stars on the main sequence (MS) produces a sudden ample supply of UV photons that changes the sound speed to ∼ 10 km s −1 as the gas becomes photoionised. At the same time, the strong stellar winds from these stars (that reach the MS first) collide with neighboring winds generating multiple shocks that thermalize the winds causing a large overpressure in the star forming volume that ends up driving the shocked stellar wind gas as an isotropic supersonic cluster wind (v cw ∼ 1000 km s −1 ) into the surrounding medium (Tenorio-Tagle et al. 2006). This cluster wind immediately encounters a large number of neighboring condensations, some of them fully or partially ionised, and the interaction leads to another (global) reverse shock that heats up the incoming gas to very high temperatures (T ∼ 10 7 K). The resulting hot gas rapidly expands along multiple paths of least resistance and flows away into lower pressure regions engulfing the dense ionized condensations that were originally part of the parent GMC. The expansion of this continuous wind generates secondary shocks that sweep-up the low density medium between condensations producing large bubbles and super-bubbles that eventually break-up allowing the hot gas to escape and create the hierarchy of shells, tubes, and filaments that characterize GHR. The expansion eventually reaches the edge of the GMC and the cluster wind escapes into the intergalactic medium carrying with it the products of stellar evolution (Tenorio-Tagle et al. 2006).

The genesis of 30 Doradus
The morphology of the gas (nebular and X-rays) confirms that the structure of 30 Doradus was carved by the cluster wind and not by the winds of individual stars; the coherent Hα structures that we see are seldom, if ever, centered on individual stars.
We have shown that the clumpy structure of the X-ray emission is due to shadowing by dense clouds of dust and nebular gas, and also by the giant molecular cloud located NNE of R136, which we claim is in front of the nebula on the basis of the radial velocities of all its individual clumps.
The innermost dense regions of the nebular gas have been pushed into coherent structures by the cluster wind in such way that they form nozzles through which the hot gas flows towards the outer regions of the nebula thus creating the giant expanding shells that characterize 30 Doradus. The Chandra X-ray images of Townsley et al. (2006) show regions -notably the large Northern plume -where the cluster wind is already escaping from the nebula. The X-ray gas occupies regions of low Hα surface brightness where the profiles are inevitably double, clearly showing that the expanding shells are being "inflated" by the hot cluster wind.
30 Doradus has a clear core-halo structure evidenced mainly by the brightness. However, the overall kinematics of the core and the halo are surprisingly similar; the integrated profiles are virtually identical and only differ in the broad wings albeit slightly. This may just be a coincidence given that the way the cluster winds carve the original molecular cloud must depend on the original configuration, which certainly must vary from one nebula to the next. It would be interesting, however, to verify whether other GHR in the Local Group present similar "coincidences".
The remnants of the dense clumps of the placental GMC have been piled together by the cluster wind to form the densest parts of the nebular gas, mostly in the form of long twisted filaments. At small scales, the forbidden line-profiles of these shells and filaments are supersonic and composed of many narrower components, which could be due to the detailed hydrodynamics of the expanding shells.
The source of energy for the whole process, therefore, is the ionizing cluster: the kinematics through the cluster winds and the emission-line luminosity through the ionizing radiation. Ultimately, this explains why L and σ are so tightly related, although it remains hard to explain the steep slope of the L − σ relation. The relation implies that the velocity dispersion of the gas increases very slowly with the mass of the ionizing cluster (σ ∼ M 1/5 cl ), so explaining the L − σ relation remains a challenging problem in astrophysics. Pure gravity would have been much simpler!

Summary and conclusions
Theoretical studies predict that the formation and evolution of Giant HII regions is determined by hot cluster winds that result from the merging and thermalization of the winds from the individual ionizing stars. Our observations largely support this scenario: none of the large nebular structures that characterize the nebula, with the exception of the central bubble that surrounds R136 are centered around individual massive stars. The hierarchical structure of the shells, on the other hand, is exactly what is expected from a cluster wind that originates in the core of the ionizing cluster (R136) and breaks out of successive shells built by the material swept up by the same wind.
We have confirmed previous results that the kinetic energy of the nebular gas is mostly in the form of the expansion of the multiple shells. Not only are the Hα profiles in the low surface brightness regions inside these shells double or multiple, but also the profiles of the rims of these shells, which are mostly single, clearly show organized expansion motions. Expansion accounts for ∼ 70% of the kinetic energy while the remainder is manifested by the supersonic intrinsic widths of the individual profiles.
The radial dependence of the kinematics and of the Hα surface brightness clearly show a well defined core-halo structure. The nebular core (R ≤ 25pc) contains the bulk of the ionizing stars and 23% of the total Hα luminosity. The massive cluster wind that carves the nebular structures originates in the core and is channeled to the external giant expanding bubbles through break-out structures within the core, which act as nozzles.
The clumpiness of the (hot thermalized cluster wind) X-ray emitting gas is due to shadowing by foreground nebular gas and dust, so it is reasonable to conclude that the hot gas uniformly fills the nebular core, and is expanding through the break-out nozzles to inflate the characteristic giant (Chu-Kennicutt) external bubbles of 30 Doradus.
The structure function of the turbulence of the nebular gas in 30Dor is flat. The idea that the supersonic turbulence of the gas in Giant HII Regions is due to a turbulent Kolmogorov-like kinetic energy cascade put forward by Medina-Tanco et al. (1997) can be definitively ruled out. Instead, we find that the structure function of a virialzed cluster is more compatible with our observations, albeit only within the cluster core where the velocity dispersion is constant.
The virial mass corresponding to the observed global velocity dispersion of σ 30Dor = 23.5 km s −1 and an effective radius of R e f f = 25pc is M vir = 3.2 × 10 6 M . The main source of the supersonic velocities observed in the nebular gas is the hot cluster wind, with a small but significant contribution from the gravitational potential of ionized gas, molecular gas, and stars. The total (photometric) mass of 30 Doradus is a few times 10 5 M with large uncertainties but certainly less than 10 6 M . all ionizing stars are in R136 (Castro et al. 2018), so it is safer to assume that N ν ∝ R −α where α measures the central concentration of the ionizing flux. Thus, the ratio of [NII]/Hα scales as R α−0.7 . Figure A.2 plots this ratio as a function of radius for the single profiles. The least-squares fit shows that F([NII]/F(Hα) ∼ R 0.2 in the halo and ∼ R 0.7 in the nebular core (we used only [NII] that has more points). This indicates that α = 1.4 in the core and α = 0.9 in the halo. As an exercise, keen readers can compare these predictions with the observed distribution of massive stars using the catalog of 30Dor stars (Castro et al. 2018). In principle, therefore, one could expect the difference between the widths of [NII] and Hα to depend on radius, but this is not the case. Figure A.3 plots the difference in the widths of [NII] and Hα for the single lines in our sample. We observe no obvious trend. We exclude, therefore, that the Hα lines are broader due to geometrical effects where [NII] (and [SII]) probe mostly the kinematics of nebular halo whereas Hα probes the entire nebula.   Townsley et al. (2006). The infrared emission from Spitzer is shown in inorange and the X-ray emission in blue. Right: Chandra 20-year anniversary 3-channel X-ray image of the same region. The images have been astrometrically matched using point-like X-ray sources and the software package ALADIN. The Chandra image is rotated ∼ 15 deg clockwise relative to the others due to the ignorance of the lead author about the intricacies of ALADIN. As in Fig.13 in the main text, the red and blue dots show the positions of red and blue shifted single profiles; the green open rhombs mark positions of all the CO clumps identified by Indebetouw et al. (2013) in the giant molecular cloud 30Dor-10