© G. F. Gahm et al. 2022

1 Introduction

There are in our galaxy numerous H II regions surrounding young stellar clusters that contain several massive O and B stars. Far- and extreme-UV radiation and stellar wind from O and B stars heat, ionise, and shape the surrounding gas, and these expanding bubbles drive compressed molecular gas outwards. The evolution of such regions has been subject to a large number of theoretical studies, starting with the pioneering work by Strömgren (1939) (and e.g. the reviews in Henney 2007; Krumholz et al. 2014; Dale 2015; Haworth et al. 2018).

Once the first O stars have formed, the surrounding molecular cloud is compressed into a shell that is accelerated outwards. The shell breaks up into filaments, elephant trunks, and globules, and eventually it slows down and disperses. The expansion velocity of the shell obtained from simulations depends on the physical processes assumed, the initial conditions, and the source driving the expansion. Numerical studies show that the shells surrounding single stars rarely expand with velocities greater than 5 km s⁻¹ (e.g. Hosokawa & Inutsuka 2006; Arthur et al. 2011; Bisbas et al. 2015), although shell velocities of the order ~10 km s⁻¹ can be attained for very massive stars (Geen et al. 2020) ora low-density ambient interstellarmedium (ISM; Freyer et al. 2003). A compact cluster of massive stars can drive a single, common shell, and these shells can reach velocities in excess of 10 km s⁻¹ (e.g. Dale et al. 2014; Lancaster et al. 2021a), but this requires a large population of high-mass stars.

From existing estimates of the central velocities of various H II regions in comparison to those measured for molecular line emission in the surrounding shells, it appears that the expansion velocities as a rule are of the order of a few km s⁻¹,in congruence with the predictions from some of the above-mentioned theoretical model simulations. However, we have found several regions for which one could suspect much higher velocities. A closer inspection shows that published estimates of central velocities and shell velocities are in fact rather uncertain. In this paper we focus on the Sharpless 171 (S 171) complex, which has the young cluster Berkeley 59 (Be 59) at its centre.

Literature data on the kinematics of the stars, ionised gas, and shell structures in this complex suggest that the molecular shell expands at a very high velocity. However, existing maps of the molecular line emission from shell structures are rather coarse, some line profiles are complex, and only a few stars have so far been measured for radial velocity (RV). We therefore decided to more closely inspect the kinematics of stars and nebular features in the complex. We collected high resolution spectra, at the Nordic Optical Telescope (NOT), of a number of stars in an area at the central cluster to obtain a mean central velocity of the cluster and information on spectroscopic binarity and membership. Moreover, we mapped the molecular line emission from the associated shell surrounding the cluster at Onsala Space Observatory (OSO). The molecular line velocities observed from different parts of the shell in comparison to the cluster velocity should provide an estimate of the expansion velocity of the complex.

The S 171 complex is described in Sect. 2 with an account of existing information on the velocities of the stars, the ionised gas, and shell structures. The observations and the data reduction are described in Sect. 3. The general results are given in Sect. 4 and discussed in more detail in Sect. 5. We compare these results with numerical simulations in Sect. 6 and end with a summary in Sect. 7.

Fig. 1 Optical image from the Digitized Sky Survey (left) and mid-IR image from Wide-field Infrared Survey Explorer (WISE, right), featuring the H II region S 171 in relation to the cluster Berkeley 59, which is located in the central white box. The northern part of the nebula is called NGC 7822. The large circle has a radius of 1.6° and shows the approximate extent of the H II region. The green contours show the 4 mm (70 GHz) continuum map from Planck. The smaller circles in the optical image mark the positions of the hottest members of the extended Cep OB 4 association; light blue circles represent spectral types O9 and B0, and yellow circles types B1 and B2. North is up, and east to the left. The WISE image is a colour composition of red (22 µm), green (4.6 µm), and blue (3.4 µm). Figure 2 is a zoomed-in view of the area in the white box.

Fig. 2 Locations of program stars in the Berkeley 59 region according to the designations in Table 3 on an image from the Digitized Sky Survey. The cluster centre is put between the bright stars A (BD+66° 1674) and B (BD+66° 1675). The white arrow points at BD+66° 1673, V747 Cep. The red arrow marks the position of a bright IR source. North is up, and east is to the left.

2 The S 171 complex

Figure 1 shows large-scale optical and infrared (IR) images of the S 171 nebula surrounding Berkeley 59 (galactic coordinates l^II = 118.2°, b^II = 5.0°). The nebula has a rather circular shape with an approximate extent as depicted by the large circle. The emission is less intense in the south-eastern part, and foreground obscuring clouds are scattered over the nebula. The arc of enhanced emission to the north and north-west have been listed as separate nebulae, NGC 7822 and the Cepheus Loop, but are clearly part of the same H II region. In the right panel of Fig. 1 we have marked the location of the cluster NGC 7762, which lies at the north-west border of the H II region and seems to have a similar distance as Berkeley 59 (Liu & Pang 2019), but according to Patat & Carraro (1995) it is an intermediate-age open cluster of 1.8 Gyr, and has a very different RV (Casamiquela et al. 2016), for this reason we regard it as unrelated to Berkeley 59. A more widespread association of early type stars, called Cep OB4, was recognised by MacConnell (1968), and members of spectral types earlier than B3 are marked in the left image of Fig. 1 (reviewed in Kun et al. 2008). Berkeley 59 is well confined within the area outlined by the central box, which is enlarged in Fig. 2. Many fainter pre-main-sequence stars are associated with the cluster (Cohen & Kuhi 1976; Pandey et al. 2008; Rosvick & Majaess 2013; Getman et al. 2017; Panwar et al. 2018). After submitting this paper we became aware of the work by Mintz et al. (2021), a major IR study of the young stellar population and its spatial distribution, where they find the youngest sources located in the surrounding shells and pillars and the older ones in the bubble surrounding the OB association, from which they suggest that the expanding HII region may have triggered star formation in Cep OB4.

Maps of the continuous radio emission show that S 171 has developed into a more or less spherical H II region of thermal emission but with a void in the south-eastern part. Figure 1 shows the 4 mm Planck map as contours overlaid on the IR image, and these are morphologically very similar to the 9 cm radio continuum map from Rosano et al. (1980). NGC 7822 and the bright patch seen close to the south-east from the centre in the optical image in Fig. 1 show enhanced radio emission (Churchwell & Felli 1970). The central, most intense emission comes from an area offset by about 8′ to the south-west of the cluster and was mapped in more detail in Angerhofer et al. (1977) and Harten et al. (1981). Isophotes of forbidden line emission are congruent with the radio maps (Lozinskaya et al. 1987). A review of works related to the stellar content and the nebula is given in Kun et al. (2008).

The general properties of the S 171 complex are given in Table 1, where we have adopted a central position for the complex as between the two bright stars in Be 59 (A and B in Fig. 2). The radius of the H II region given in Col. 4 is based on the distance to the complex cited in Col. 3 and corresponds to the circle in Fig. 1 with a radius of 1.6°. Column 6 gives different estimates of the mean heliocentric RV of the stellar cluster, and the corresponding range in velocity is expressed in local standard of rest (lsr) in Col. 7. The last column gives published central velocities of the ionised gas, mostly from radio recombination lines but also from Ha and [N II] lines.

Published cluster velocities are not very accurate and are based on only two or three stars, and they differ considerably from the central velocities obtained for the ionised gas. The latter are also discordant, and very different velocities were obtained for different positions over the nebula (e.g. Pedlar 1980; Rosano et al. 1980; Lozinskaya et al. 1987). The lines are broad, full width at half maximum (FWHM) > 30 km s⁻¹, with complex line profiles in some directions, which makes it difficult to assign a precise value of the central velocity.

The northern molecular cloud related to NGC 7822 was mapped in CO emission in Elmegreen et al. (1978) and Wootten et al. (1983), and the shell structures related to the S 171 nebula were mapped in CO by Leisawitz et al. (1989); Yang & Fukui (1992); Yonekura et al. (1997), and in H₂CO by Rossano et al. (1983). These maps cover most of the complex with grid spacings of several arcminutes (Yang & Fukui 1992, 2′ for the central part in). The associated velocity patterns turned out to be rather complex with multiple velocity components at several positions, ranging in RVs from −1 km s⁻¹ to −19 km s⁻¹ (lsr) in some cloudlets. Considering also the viewing aspect angle of these cloudlets it appears from some of the central velocities listed in Table 1 as parts of the shell may expand at very high velocities.

We can assume that the present structure of the S 171 complex is the result mainly of the interaction between the hot members in Be 59 with the surrounding molecular shell. Since the central velocity of the complex is not well established, and existing CO maps are rather coarse, we found it motivated to obtain new RVs of several stars in Be 59 and of the molecular line emission from the shell at a higher spatial resolution than obtained before.

3 Observations, Reductions, and Measurements

3.1 Optical Spectroscopy

We aimed at including all OB stars in the Berkeley 59 cluster in our sample and furthermore reaching as low a mass as possible, but due to extinction and the signal-to-noise ratio (S/N) needed in the spectra to determine RVs, a practical limit was set around V ~ 14.5 mag. A few stars had determined spectral types from previous studies, and for the rest we made rough estimates based on IR and optical photometry (Pandey et al. 2008) to select the program stars (Sect. 4). Optical spectra were collected for 27 stars in the Be 59 region at the 2.6 m NOT during two observing runs in October 2015 and August 2016 plus spectra obtained as filler programs in service nights. We used the high-resolution FIber-fed Echelle Spectrograph (FIES; Telting et al. 2014), with a resolution of R ~ 25 000 covering the wavelength range from 3700 to 7300 Å when used with Charge Capture Device 13 (CCD13) and 3620–8580 Å when used with CCD 15, which was installed in October 2016. Our program stars are marked in Fig. 2 with numbers as defined in the subsequent section.

Most targets were observed on at least two different dates (separated by a minimum of four days) to check for spectroscopic binaries. Because our targets are mainly early type stars with few available lines, often strongly rotationally broadened, we aimed for a S/N of ~50 in the final spectra. However, our targets are heavily reddened (typically by 5–7 visual magnitudes). This lowers the S/N significantly in the blue part of the spectra, and thereby the number of useful lines. Exposure times were up to 1800 s, and for the faintest targets two consecutive exposures were taken. From the list in Fekel (1999) we selected as early-type RV standards HR 1810 (B2.5 V), HR 7512 (B8 III), HR 8404 (B9.5 V), and HR 7903 (A0 III), stars not found to vary in velocity and with accuracies in the range 0.05 to 0.4 km s⁻¹. We also observed Vega (A0 V), and for later spectral types we used the RV standard HD 31253 (F8). In addition, we collected spectra of a number of MK standards used for spectral classification of our targets.

Standard calibrations (21 halogen flats, seven biases, and one ThAr) were obtained each afternoon. Data reduction was made using the FIEStool pipeline (Stempels 2005). Over the course of our observations, FIES was upgraded with a new CCD, in October 2016, and new octagonal fibres, in June 2017. The merging of the orders in FIEStool was particularly troublesome in the spectra of blue stars with the old setup but was enhanced with the new one. Thus, the merging was improved manually by cutting the noisy edges of individual orders. We then normalised the spectra by fitting a continuum.

We measured the stellar RVs using both cross-correlation techniques (‘fxcor’ in IRAF) as well as direct measurements of line centres by fitting a Gaussian profile to each line (‘rvid-lines’ in IRAF). For most OB-type stars, we mainly used the He I and He II lines when available, while for later B and early A stars, and for targets where few lines were available, we also included the hydrogen lines. When both methods were applicable, we compared the results and typically found consistent values within the errors. The instrument drift is practically negligible compared to the typical measurement errors in the RVs of our targets. Multi-epoch RV results for the program stars are listed in Tables A.1.

3.2 Molecular Line Observations

We observed ¹³CO J = 1–0 at 110.201 GHz from the region surrounding Be 59 during several runs in 2016 and 2017 using a 3 mm SIS receiver at the 20 m telescope at OSO. Frequency-switching mode was used with a 1600-channel hybrid digital autocorrelation spectrometer with a channel spacing of 25 kHz (0.07 km s⁻¹), a bandwidth of 40 MHz, and the velocity resolution is about 0.2 km s⁻¹ after smoothing. The single-sideband system temperature was typically around 150–400 K. The pointing was checked regularly, and we estimate the pointing error to be less than a few arcseconds. At 110 GHz, the FWHM beam size of the 20 m antenna is 34″, and the main beam efficiency ~0.30 (for an average elevation of approximately 30°). The chopper wheel method was used for the intensity calibration. The data reduction was performed with the spectral line software package xs¹.

Since the complex extends over a very large area in the sky, we had to select primarily areas where distinct shell structures are seen as dark foreground features in optical images. Some of these features were listed in the catalogue of dark nebulae in Lynds (1962). A total of 489 positions were observed including both observations over grids and at single positions, and some areas have a denser spatial coverage than others. In addition, observations of positions located outside obscuring clouds provided reference spectra.

4 Results

4.1 Program Stars and Parallaxes Plus Proper Motions

Most of our program stars in Be 59 are common to the stars observed by Pandey et al. (2008), and we chose their sequential number as prime designation (P). The stars are marked accordingly on the image in Fig. 2. The program stars were selected with the purpose of determining the stellar RVs and are therefore a brightness limited sample. However, we searched the mid-IR WISE point source catalogue (Cutri et al. 2013) for early type stars extinguished behind dark clouds within a radius of 8 arc minutes of the cluster centre. The star we name IR star (identical with WISE J000126.63+672346.1) was located in the direction of very opaque shell structures obscuring the bright nebulosity and although faint in the visual, its brightness at mid-IR wavelengths was at the same level as for the well known O stars in the cluster, for this reason we included it in our program, and our multi-epoch spectroscopy will reveal it to be a binary O9 V star.

The program stars are listed according to increasing right ascension (RA) in Table 2 and include designations from MacConnell (1968, MaxC), Sharma et al. (Sharma et al. 2007, SOP), and the Tycho mission (Høg et al. 2000, TYC) with the parallaxes and proper motions, when available, from the Gaia Data Release 2 (proper motion) and 3 (parallaxes; Gaia Collaboration 2016, 2018, 2021). Numbers are in milliarcseconds (mas) and quoted errors in parentheses. The parallax of the IR star is consistent with membership in Be 59. The parallaxes are well confined around a mean of 0.9 ± 0.1 mas, corresponding to a distance of 1.1 kpc excluding six stars, which appear to be foreground stars. We adopt this distance for the present study. The corresponding radius of the H II region is 31 pc as listed in Table 1.

From the proper motions, excluding the six stars as above, we obtain a mean value for the cluster of Δα = −1.7 and Aδ = −2.0 mas yr^–1. The spread around this mean is small. Our cluster mean values, also regarding the parallax, are in good agreement with the values obtained by Kuhn et al. (2019) based on 225 stars in Be 59.

4.2 Stellar Properties and Radial Velocities

We found a number of program stars to be spectroscopic binaries, including the IR star, and for most of these we extended the observations to get a hand on orbital data. Examples of normalised spectra collected at NOT of some program stars for two spectral regions are shown in Fig. 3 with some spectral features discussed in the text marked. The third and fourth spectra from top show the double-lined spectroscopic binary BD+66° 1674 at two phases when the stellar lines from the two components merge and split. For stars considered as members the diffuse interstellar bands (DIBs) are of about the same strength but weaker in foreground stars as evident for the star P316 at bottom of the panels. As seen in the bottom panel narrow Ha line emission, flanked by [N II] emission in some cases, are superimposed on the broad Hα absorption lines. These emission lines originate in the surrounding emission nebula.

Most program stars had not been classified on the MK system before, which we did by comparisons to standard stars based on traditional line ratios in the blue spectral region (Gray & Corbally 2009), also for the secondary components in binaries observed at large line splittings. For OB stars (Walborn 2009) we used He II to He I equivalent width ratios in the blue whenever possible, and, despite multiple epoch spectra, there are uncertainties in the line ratios due to blending of the binary or triple components. The blue part of the spectrum of the reddened IR star is noisy, but using ratios based on equivalent widths of lines in the yellow and red parts we conclude that it is of spectral type O9 V (lines of e.g. He I, He II, O III 5592 Å, C IV 5801, and 5812 Å). These line ratios were measured also for all early type stars for comparison. The most luminous star is the eclipsing and spectroscopic binary V747 Cep of spectral type O5.5 V.

In Table 3, we present the results from our spectroscopic survey, where the stars are ordered according to spectral type. Photometric data from the literature are in Column 2. MK spectral classes as listed in the Set of Identifications, Measurements, and Bibliography for Astronomical Data (SIMBAD) Astronomical Database are in Column 3 and our estimates in Column 4. Column 5 gives information on the strength of the DIBs. Column 6 gives the visual extinction derived from the Two Micron All Sky Survey (2MASS) photometric data combined with intrinsic J – H colours inferred from the spectral types obtained from Koornneef (1983). We assume that the excess in the J – H colour is a measure of the extinction and use the relation A_V = R_V × E(J – H)/R(J – H) =11.9 × E(J – H), for R_V= 3.1 and the empirical 2MASS extinction coefficient R(J – H) =0.26 (Yuan et al. 2013). Column 7 lists heliocentric RVs from the literature and in Column 8 are our measurements with corresponding rms errors. For stars observed in only one night, the errors are from the determination of the RV from that single spectrum. Spectroscopic binaries are denoted as single-lined (SB 1), double-lined (SB2), and triple-lined (SB3), and these will be discussed in more detail in the subsequent section. In Column 9 are the number of nights when FIES spectra were taken if more than one night. Individual measurements are listed in Tables A.1.

Based on the distance derived from the Gaia parallaxes and, if not available, from the photometry combined with spectral type and extinction, the strength of the diffuse interstellar absorption bands, and RVs, we find that all stars of spectral type A1 and earlier are likely members of Be 59. Those that we regard as non-members, are marked in Column 10, which also contains notes on stars with broad He or metal lines (br), for which we measured v sin i > 150 km s⁻¹.

Figure 2 might give the impression that some members in Be 59 are not located behind foreground shell structures. Yet, they suffer from relatively large extinctions. The optical image indicates fragments of foreground clouds seen in absorption against a cavity illuminated by the Hα nebulosity, which is even more clearly displayed in the IR in Figs. 1 and 10, where the contour overlay shows radio continuum emission enhanced at the location of the central cluster. Towards the west there is a pronounced ridge of dust behind which the IR star is located, and as shown in Table 3 the measured extinction towards individual stars is highly variable over the cluster area, but A_V > 3 mag for all cluster members.

4.3 Completeness of OB stars

Next we wanted to determine how complete our sample of OB stars was. B3 V stars have absolute magnitudes M_V = −1.6 ± 1.3 mag (Wegner 2006) and apparent magnitudes m_V = M_V + dm + A_V, where dm = 10.21 mag for Be 59. Our magnitude limit for spectroscopy was set at V ~ 14.5 mag, which means we can reach completeness for B3 V stars only for regions with A_V ≤ 6 mag, and completeness for O7 V stars for A_V ≤ 9 mag. Among the most reddened sources in our sample is P16 and the IR star with visual extinctions slightly >9 mag while the median A_V value is 5.1 mag. From the WISE catalogue search within a radius of 8′ (2.6 pc from the cluster centre) we find 38 sources brighter than 9.7 mag in the W1 (4.3 µm) band. This is the apparent magnitude of a B3 V star if extinguished by A_V = 5.1 mag, equivalent to A_W1 = 0.30 mag (Yuan et al. 2013). Twelve of these are among our optically selected program stars (all of which are later found to be cluster members), and the brightest of the remaining 26 sources is the IR star (W1 = 7.15 ± 0.04 mag). The rest constitute a mix of intermediate-mass stars and background giants, but a few could be reddened B stars belonging to the cluster, although all are fainter than the IR star by 1 mag or more at 3.4 µm. Based on this we conclude that our sample may not be complete all the way to B3 stars, but most likely includes all O stars in the area we studied.

Fig. 3 Examples of normalised spectra obtained for some program stars over two selected wavelength regions. Certain spectral features discussed in the text are marked. The spectra are ordered in steps of 0.3 (intensity scale) in the left panel and 1.0 in the right panel and with a smaller spacing between the third and fourth spectra; the double-lined spectroscopic binary BD+66° 1674 is shown at two different phases. To improve the S/N for the faint IR star, its left spectrum is a combination of many observations, Doppler-corrected for the binary motion, which smears out the DIBs.

4.4 Spectroscopic Binaries

The frequency of spectroscopic binaries can be expected to be high for the most massive stars. According to Chini et al. (2012) stars with masses >16 M_⊙ as many as 82% are close binary systems, while this fraction drops to about 20% for 3 M_⊙ stars. In order to check for binarity we obtained repeated observations for the earliest type stars (Table 3).

BD+66°1675 and BD+66°1674 were regarded as spectroscopic binaries in Liu et al. (1989) and Crampton & Fisher (1974). We found that BD+66° 1675 is a triple-lined spectroscopic binary (SB3) and our data points suggest a period of 74 days and a systemic velocity of −10 km s⁻¹, as well as semi-amplitudes of 6, 150, and 280 km s⁻¹ for components 1, 2, and 3, respectively. While components 1 and 2 have He II lines, only the He I lines are detected in the faintest component, indicating spectral type B (Fig. A.2).

BD+66°1674 is a double-lined spectroscopic binary (SB2). However, individual lines are often heavily blended, and we could not extract associated periods and amplitudes. There are perhaps three components.

V747 Cep (BD+66°1673) is an eclipsing binary with a photometric period of P = 5.33146 days, an orbital eccentricity of e = 0.3, and an inclination i ~ 75° according to Majaess et al. (2008). Maíz Apellániz et al. (2019) assigned it a spectral type of O5.5 V(n)((f)) and found that it is also a spectroscopic binary. We can confirm that V747 Cep is a spectroscopic binary. In the data from two nights we note an indication of line splitting, but more observations are needed to confirm whether it is of type SB2 (Tables A.1 and Fig. A.1). The RV data fit well with the photometric period and suggest a semi-amplitude of K₁ = 90 km s⁻¹ around a systemic velocity of −26 km s⁻¹ (Fig. 4, top panel), giving a close orbit with a semi-major axis of only 0.04 AU or ~10 R_⊙. At the time of submission, we became aware of a more detailed study of this target in Trigueros Páez et al. (2021).

Multiple spectra of the hidden O star we call the IR star (2MASS J00012663+6723461) demonstrates that it is a spectroscopic binary, for which we estimate a period of 9.51 days and K₁ = 45 km s⁻¹, assuming a circular orbit, giving a projected semi-major axis a sin i around 0.04 AU.

For P16, with only seven observations, only preliminary values can be assigned: a period of 13 days with a semi-amplitude of 12 km s⁻¹. Examples of phase diagrams are shown in Fig 4. Finally, the late type foreground star, P201 (G5 V), show two velocity components and is likely an SB2.

4.5 Mean Cluster Velocity

Excluding spectroscopic binaries and non-members from Table 3, we obtain a mean RV of the cluster of −18.75 ± 3.3 km s⁻¹ (−9.5 km s⁻¹ in lsr), which implies a standard error on the mean of 0.9 km s⁻¹ if we assume a Gaussian distribution. Hence, our derived mean cluster velocity is close to that derived for the central velocity of the ionised gas by Pedlar (1980) as quoted in Table 1.

Fig. 4 Phase diagrams for (from the top): V747 Cep for a period of 5.33146 days; BD+66° 1675, period of 74 days; the IR star, period of 9.51 days; and P16, period of 13 days.

4.6 The Molecular Line Emission

We obtained ¹³CO spectra from a total of 489 positions in the areas labelled in Fig. 5. The left panel shows a large area with the complex at the centre. The circles mark single positions observed along two dusty striations of dust that extend from the central parts to the east, and which are seen as distinct features, especially on the photographic version of the Palomar Sky-Atlas. Field E4 is located about 3° east of the cluster and was originally included as a reference position outside the complex. Since we discovered a rather strong signal at high negative velocity here, we extended the observations over what appears to be isolated dust clouds. The right panel shows the areas observed in the central region of S 171. The circle at ‘DQ’ marks the position of an elephant trunk, called the Dancing Queen, mapped in detail in ¹²CO and ¹³CO in Gahm et al. (2006).

Examples of spectra are shown in Figs. 6 and 7. Signals from the molecular gas related to the foreground shell is seen in most panels for W8 in Fig. 6 at v_lsr~ −18 km s⁻¹ and at v_lsr~ −13 km s⁻¹ for E2 in Fig. 7. These components are blue-shifted relative to the mean cluster velocity as derived in Sect. 4.5 and are present in the bulk of our spectra. In the following, we refer to components at velocities ≤−10 km s⁻¹ as V_b. Especially in the central fields these lines can be split into two or three components well separated in velocity or blended as in Fig. 7 (shell front). There is another component present at certain positions at velocities around −4 km s⁻¹ (as in Fig. 7), slightly red-shifted relative to the cluster, which we refer to as V_r. As will be discussed in Sect. 5 these could be related to the remote side of the complex. Representative spectra for each region are presented in Figs. C.1 and C.2.

Our spectra also contain additional components, which are not related to the S 171 complex. In most directions we got signals at velocities around −2 to +4 km s⁻¹, which we identify with the local gas, called feature A, related to the Gould’s Belt (e.g. Lindblad et al. 1973). Components at v_lsr ~ − 6 km s⁻¹ are related to gas in the local spiral arm, the Orion arm, sometimes referred to as ‘the other local feature’. As expected this component declines statistically in intensity with increasing galactic latitude. In our optical spectra of cluster members, interstellar absorption lines of Na I D are present in practically all stars. Also, these lines share the velocity of this local gas. Both local components are present in the majority of our ¹³CO spectra (Fig. 7). Emission from the background Perseus arm, at expected velocities of around ≤−30 km s⁻¹, does not enter our spectra since the S 171 complex is at a galactic latitude of ~+5º.

We made Gaussian fits to the line profiles, from which we extracted peak antenna temperature $T_{\text{A}}^{*}$ (outside the earth atmosphere), central velocity, FWHM line widths, and for more complex profiles we decomposed the signals into separate components. The 1σ noise level at 0.2 km s⁻¹ resolution is normally <40 mK, but reaches above 200 mK at several offset positions in W1, 6, 7, and E2. The V_b and V_r components are as rule rather narrow, and in most cases the line widths are <1.5 km s⁻¹ (FWHM).

For most areas in Fig. 5 we defined several central positions from which offset observations were made. The details of observed central positions and offsets can be found in Table B.3, where Col. 3 gives the total number of observed positions for each sub-area. In Col. 4 all offsets in RA are given in arcminutes, and for each of these the observed offsets in declination (Dec) are given in Col. 5. In Cols. 6 and 7 we list those offsets in Dec (for each offset in RA) where we have detected a V_b component (≤ −10 km s⁻¹) and a V_r component (velocity interval −2 to −6 km s⁻¹). Non-detections are marked with a dots.

The amount of data so obtained is extensive since V_b components are present in the majority of the positions, and for an overview we have chosen to present the ranges of different parameters for each area in Table 4. Here, the total number of observed positions is given in Col. 2. The range in velocity for the V_b component is given in Col. 3 with the number of positions where it was detected in parentheses. Column 4 gives the maximum peak antenna temperature for the V_b component observed in the area. Within each area the range in $T_{\text{A}}^{*}$ can be large. Cols. 5 and 6 gives the corresponding values for the V_r component.

Finally, another circumstance to consider when interpreting the data is that there is a system of molecular clouds that move at high negative velocities extending several degrees along the galactic plane. This system was mapped in ¹³CO by Yonekura et al. (1997), and in Fig. 8 we have plotted their positions for clouds moving at ≤−12 km s⁻¹ together with the corresponding positions in our survey for two areas located outside the defined extent of the H II region. Hence, this system of clouds have velocities in the same range as the V_b components in the central area. This therefore opens the question of which clouds are associated with the S 171 nebula and which are not.

None of these clouds at significant negative velocities appear as distinct dark clouds in optical images, and none is listed in the catalogue of dark clouds in Lynds (1962). It appears that this system of clouds are background objects, possibly at the far side of the Orion arm, where gas can be expected to stream into the spiral arm. We suppose that the most eastern cloudlets in our survey, areas E4 and Esouth:E marked in Fig. 8, may belong to the same system.

For the central clouds in Fig. 5, residing in areas from W10 to E3 (west to east) and S4 to N1 (south to north), there is no doubt that all belong to the dusty shell surrounding Be 59. Some of these clouds are seen in silhouette against the bright nebulosity and have bright rims and associated elephant trunks, and most individual cloudlets are connected to each other with dark nebulosity.

Fig. 5 Map of the fields we obtained CO spectra in. Left: large region, 4.4 × 4.4 degrees, surrounding the emission nebula S 171. The arrow points at the position of the central cluster. The northern arc of bright nebulosity is NGC 7822. In this panel the more peripheral areas observed for 13CO emission are marked and labelled. The circles (Enorth and Esouth) mark observations at selected positions along two dusty striations that extend to the east from the centre. The dashed box marks the central region shown in the right panel. Right: image spanning 1.9 × 1.2 degrees. The core of the stellar cluster is inside the dashed box. The position of an elephant trunk, called the ‘Dancing Queen’ is marked with a black circle. North is up and east to the left in both images, which are from the red Deep Sky Survey.

Fig. 6 Part of the map of 13CO spectra obtained for area W8. The offsets in RA and Dec from the central position are in steps of 60″. The velocity interval, in km s−1, and the antenna temperature, TΑ in Kelvin, are indicated in the upper-right panel. The component at around −18 km s−1, present in most panels, is related to the foreground shell, whereas the component at around −7 km s−1 has a local origin.

Fig. 7 13CO spectra from one position in area E2. The strongest signal at around −13 km s−1 comes from the foreground shell, while the components at −6 and +2 km s−1 are related to gas connected to the Orion arm and Gould’s Belt. The component at around −2 km s−1 may trace the remote part of the molecular shell.

5 Discussion

5.1 Kinematics of the Molecular Shell

The RV pattern over the molecular gas in the S 171 region is complex, as noted also in previous investigations referred to in Sect. 2. In several areas listed in Table 4 we find two distinct V_b components (marked with an asterisk) at certain positions indicating that there are at least two systems of expanding clouds in the foreground of the complex. The lines are narrow but their velocities can change over small areas over our grids. When observed with a larger beam these components merge, and this appears to be the case for the broad plateau emission observed with a beam of 2′.7 in Yang & Fukui (1992) in Wl and W2 (their Fig. 3), which is resolved in our maps into narrow components drifting in velocity in RA and Dec.

As described in Sect. 4.6 individual spectra usually show several velocity components. Besides the two groups of components that we considered as not related to the S 171 complex we distinguish four kinematically separated groups related to the expanding shell. In Fig. 9, these components are plotted as RV relative to the mean cluster velocity of −9.5 km s⁻¹ against angular distance from the centre for different sub-areas (listed in Table B), where a given component was detected. Each point represents an average from all detections within a given sub-area. The spread in velocity within a given sub-area is small, ≤1 km s⁻¹.

The blue dots represent structures expanding at the highest speed in our direction when assuming a spherical shell expanding at v_exp = 12 km s⁻¹ and a radius of 31 pc (dotted curve). The yellow dots are consistent with a more moderate expansion velocity of v_exp = 4 km s⁻¹. The red squares in Fig. 9, represent components that are red-shifted relative to the cluster by a few km s⁻¹. Finally, we identify yet another component present in the central region (W1, W2, and W3) with velocities between the high and moderate velocity groups (marked with green circles). In these spectra the components associated with moderate velocities (yellow dots) are also present. The positions of clouds with implied high expansion velocities are shown as blue circles in both panels in Fig. 10 and those at intermediate velocities are all in the right panel as green circles.

While the blue dots in (Fig. 9) are reasonable well described by a spherically system of clouds expanding at high velocity and located at the outskirts of the H II region, the case is different for clouds moving at moderate velocities (yellow points). Not only is this system less extensive than the high-velocity system, also the outermost points fall below the curve assuming spherical expansion and a similar radius. This cut off is indicated by the dashed curve, where RVs approach 0 km s⁻¹ at a radius of ~25 pc. Hence, it appears that this system does not extend as far out from the centre as the high-velocity system. Moreover, it appears to have a counterpart in clouds moving at similar velocities but red-shifted relative to the cluster (red squares in Fig. 9). These red-shifted clouds are also confined within a projected distances of 25 pc from the centre. A natural interpretation is that they are part of the same expanding shell and that we detect molecular emission both from its near and far sides. It is hard to judge if there is also a remote counterpart to the high-velocity component because it would fall in the spectral region where the local components contribute. However, we note that at some positions such components are split in two.

Clouds moving at moderate negative velocity are found in all directions over the nebula including also the dark dust cloud associated with NGC 7822. In contrast, the high-velocity clouds are found only along a band extending from west to east and crossing the centre. Sub-area Enorth:C is located just outside the adopted radius of the complex (beyond 31 pc in Fig. 9), and there is a small cloud in the foreground of the central cluster (sub-area C:C). The area over which the intermediate velocity component is detected coincides in part with the peak intensity of the continuous radio emission as mapped in Harten et al. (1981), outlined as yellow ovals in the Fig. 10.

To summarise, the bulk of the molecular gas in front of the bright nebulosity expands at moderate velocity. However, we also localised a system of smaller cloudlets, which appear to expand with a higher velocity than expected from simulations, although slower than could be inferred at start from some estimates of the central velocity of the complex (Table 1). These are the cloudlets marked with blue dots in Fig. 9 and blue circles in Fig. 10.

Fig. 8 Positions in galactic coordinates of molecular cloudlets moving at high negative speed (≤−12 km s−1) in an extended area outside the central part of the S 171 complex (left image from Fig. 5 inserted for reference). Open squares are cloudlets from Yonekura et al. (1997) and filled squares are E4 and Esouth:E from our survey.

Fig. 9 Radial velocities relative to the cluster velocity for different sub-areas against angular distance from the centre (in all directions). Blue dots mark areas that are expanding at high velocity, vexp ≈12 km s−1, and yellow circles (smaller) areas expanding at moderate speed, vexp ≈ 4 km s−1, when assuming spherically expanding systems (dotted curves). Green circles mark areas with intermediate velocities (‘DQ’ is the Dancing Queen). Red squares show areas with components whose velocities are redshifted relative to the cluster.

5.2 Masses of shell Structures

Our ¹³CO observations do not cover the entire area of the S 171 complex. Nevertheless, we made an attempt to obtain rough estimates of the mass of individual cloudlets in the shell by the same method as used for deriving the masses of elephant trunks in Gahm et al. (2006). In short, the total mass of a given cloud is $$M\left[M_{\odot }\right]=1.35\times {10}^{-11}F\left(T_{\text{ex}}\right)D^{2}A\langle I_{13}\rangle ,$$(1)

where D (pc) is the distance, A (arcsec² ) is the cloud area, and 〈I₁₃〉 is the average integrated ¹³CO intensity over the area. We have assumed an excitation temperature T_exp = 23 K, which gives a dependence F(T_exp) = 30 (for instance Nikolić et al. 2001).

We then treated clouds with high and moderate expansion velocities separately. For some clouds we included also areas outside our grids with similar prominent extinction and assuming the same value of the average integrated intensity. Most of the mass in the foreground shell is concentrated in the very dark clouds obscuring the areas of maximum continuous radio emission just to the west of the cluster, including areas W1, S1, W2, and W3 (Fig. 10). These areas are well sampled and covered also by Yang & Fukui (1992), who derived masses for two ‘clumps’ C1 (our positions W2 and W3) and C2 (our positions W1). In comparison with their results our estimates are similar for C2 but larger by a factor of 3 for C1, when the different adopted distance to the complex is taken into account. We regard this as a reasonable agreement, and also for area N1 in NGC 7822 our derived mass of ~300 M_⊙ is similar (a factor of 1.5 larger) to the value quoted in Elmegreen et al. (1978). Our higher values are more consistent with column densities inferred from the large visual extinction of these clouds. For the clouds emitting the intermediate velocity components we estimate a mass of ~200 M_⊙.

Based on our estimates the total mass of the shell in front of the H II region is of the order of 2200 M_⊙, of which ~1500 M_⊙ is contained in the central areas described above. There is a considerable uncertainty in our mass estimates and more complete surveys of other molecular transition useful for mass determinations are warranted. Areas W9, S3, S4, and N3 do not show any components at high or moderate velocity. Either these are yet undisturbed clouds in the foreground of the nebula, or the lines are too weak for detection, which can be the case also in the relatively large areas E1 and E3, where lines at moderate expansion velocity are detected just over smaller parts of the maps.

All cloudlets that we identify as moving at the highest expansion velocities have low masses. Among these, W10A and N2 are the most massive, ~70 M_⊙ and ~50 M_⊙, while the others have very low masses, ≤ 10 M_⊙. As will be discussed in the subsequent sections, this may be the reason why these cloudlets obtained a larger speed and reached farther out from the centre than the heavier clouds moving at moderate velocities.

Fig. 10 Map of cloudlets. Left: overview of the region shown in Fig. 1 in 12 µm from WISE (grey scale) with the Planck 857 GHz map (black contours) overlaid to emphasise the dust emission. A black circle spanning a radius of 1.6° is centred on the Berkeley 59 cluster. The blue dots from Fig. 9 are shown as blue circles and mark the locations where the high negative expansion velocities are found. Likewise, the green circles refer to the intermediate and the yellow the low negative velocity components, while the red circles mark positions where the positive component is found, possibly related to the expansion on the far side. Right: close-up 55′× 40′ view of the central region where the red, green, blue coding is WISE 12 µm (red), 4.6 µm (green), and Digitized Sky Survey 2 red plates (blue). The positions of V747 Cep and the IR star are marked with white diamonds.

5.3 Dynamical Evolution of the Complex

The expansion of the S 171 region is driven by the interaction from the hot members in the Be 59 cluster. Some quantities related to the physical conditions in the complex are listed in Table 6, where the first entries contain estimates of the total luminosity, L_t, and the total flux of Lyman continuum photons, N_Lyc, for all cluster members earlier than B1. The five O stars dominate the fluxes, and the contributions from members later than B0.5 are negligible. According to Pandey et al. (2008); Lata et al. (2011); Panwar et al. (2018) the ages of pre-main-sequence stars in the general area around Be 59 indicate a spread from 0.5 to 5 Myr. However, it is also noted that the massive stars in the cluster are not significantly evolved and with inferred ages consistent to what is assumed in the present study.

The middle panel contains estimates of the electron temperature, T_e and electron number density, n_e, related to the physical state of the plasma in the central part of the H II region followed by properties of the molecular shell in front of the nebula as estimated above. Here, the total mass of the expanding molecular shell on the remote side of the complex could not be estimated.

As shown in both panels of Fig. 10, the areas of maximum continuous radio emission are coinciding with some highly obscured foreground clouds, as W1 and W2, which are offset to the west relative to the central cluster. This enhanced emission is most likely related to photodissociation regions (PDRs) hidden behind the dark clouds and excited by UV light from massive stars in Be 59. Possibly, a main contribution comes from the most luminous O star (V747 Cep) and the heavily obscured O star (the IR star), which are also offset to the west (by ~4 pc in projected distance). It is also over this general area our OSO spectra contain lines from the intermediate velocity component (green circles in Fig. 10).

A similar type of interaction takes place at the periphery of the H II region north of the cluster, where the bright nebula NGC 7822 borders a massive and opaques cloud, N1, only that in this case the PDR is visible as a bright nebula. The star marked with a circle just south of NGC 7822 in Fig. 1 was noted as a member of the Cep OB4 association by MacConnell (1968). This star (MacC 9 = V399 Cep) of spectral type B2 is at the right distance to the complex according to Gaia DR 2, but we conclude that it plays a minor role for the excitation of NGC 7822 in comparison to the central O stars.

The ionised gas is thermal (Angerhofer et al. 1977; Pedlar 1980), and it is unlikely that the large expansion velocities observed has been caused by an early supernova explosion, which besides would have disrupted the entire molecular shells.

In Sect. 5.1, we distinguished three systems of expanding shell structures, one containing rather massive clouds moving at 4 km s⁻¹ relative to the cluster and presumably with a counterpart at the remote side of the complex, another system with predominantly smaller cloudlets expanding at 12 km s⁻¹, and one central area where some clouds move with velocities between these two.

The 4 km s⁻¹ component is in line with expectations from some numerical simulations of evolving H II regions with comparable sizes at an age of 2 Myrs (e.g. Freyer et al. 2003; Hosokawa & Inutsuka 2006). On the other hand, the 12 kms⁻¹ component traces a system of clouds moving outwards at a considerably higher speed compared to theoretical predictions assuming similar total luminosities and ionising fluxes as in Table 6. As discussed above, this system of clouds extends farther out from the centre and contains cloudlets of smaller mass than those expanding at moderate velocities. For the high-velocity clouds the crossing time to reach from the centre to the periphery of the H II region at constant speed is 2.5 × 10⁶ yr, which is in line with the cluster age.

Below, we consider as a working hypothesis that the highvelocity clouds were confined and detached from the more massive structures in the expanding shell early in the expansion of the shell. In our picture, the shell is very clumpy and contains small overdense clumps. The entire shell is decelerated by accumulation of material from the ISM, but dense clumps have more inertia and are more difficult to slow down. The tiny cloudlets retained a higher speed and are now located farther out than the shell, even outside the border of the H II region (Fig. 9). The clouds are confined within a relatively narrow belt extending from the western edge of the H II region over the central part to the eastern edge. The cause for this asymmetry is unclear.

6 Numerical Simulations of the Complex

In this section, we explore to what extent we can match the complex velocity pattern observed for the molecular gas in S 171 with a theoretical model simulation. As discussed in Sect. 5.3, we have identified a system of high-velocity cloudlets of relatively small masses that have reached farther out from the centre than the more massive cloud structures that form the shell (Fig. 9). We attempt to include both these cloud systems in our simulation.

Our simulation describes how a cloud of gas is affected by the stellar winds of an embedded star cluster. This requires coupling between hydrodynamics (which describes the dynamics of the gas), gravity (which describes the dynamics of the stellar cluster), and stellar evolution (which describes the mass evolution and stellar wind output of the individual stars), plus a prescription for stellar winds.

6.1 Multiphysics Stellar Wind Bubble Model

6.1.1 Model

We simulate hydrodynamics using the fi code (Pelupessy et al. 2004, based on Hernquist & Katz 1989; Gerritsen & Icke 1997), which is based on the smoothed particle hydrodynamics (SPH) formalism. Stellar winds are implemented using the prescription by van der Helm et al. (2019), which inserts new SPH particles corresponding to the stellar wind into the hydrodynamics code. The stellar parameters required to describe the winds’ mass loss rate and velocity are prescribed by the seba code (Portegies Zwart & Verbunt 1996, with the adjustments by Toonen et al. 2012), a parameterised stellar evolution code. The stars’ dynamics are simulated using the huayno code, a second-order symplectic, direct N-body code (Pelupessy et al. 2012). All of these codes were brought together in the Astrophysical MUltipurpose Software Environment (AMUSE; Portegies Zwart & McMillan 2018). Finally, the gas and star cluster were able to affect each other through Bridge interactions (Fujii et al. 2007), which resolved gravitational interactions between the two components every 1 kyr, while leaving the internal evolution of each component unaffected.

Our initial conditions consist of a number of gas and stellar components. We include a main gas cloud of 3000 M_⊙ (to take into account unobserved material behind the H II region) and a virial radius of 6 pc distributed following a Plummer profile (Plummer 1911), which represents the centrally concentrated gas left over from star formation. We place two cloudlets of 10 and 50 M_⊙ (Plummer spheres with virial radii of 0.5 pc) at rest at x = ±3 pc, representing two cloudlets that have already detached from the main gas distribution. The velocity fields of these Plummer spheres are virialised and non-turbulent.

The stellar component consists of 10³ stars between 0.5 and 16 M_⊙ following a Salpeter initial mass function, and distributed according to a 2 pc virial radius Plummer sphere. Additionally, we include a collection of massive stars representative of the O-type (and giant B-type) stars observed. We used Sternberg et al. (2003) to associate a mass with the stellar type of each massive star. The mass assumed for every star is provided in Table 7. We do not initialise these stars as binaries, but we do concentrate them within the inner 0.1 pc of the cluster.

6.1.3 Results

We ran this simulation for 2 Myr at a resolution of about 3 × 10⁻³ M_⊙ per initial SPH particle (and 10⁻⁵ M_⊙ for the stellar wind particles). Figure 11 presents a number of snapshots of the simulation. The stellar winds from the centrally concentrated massive stars open a cavity in the gas cloud that expands with time and sweeps up the detached cloudlets. In the bottom panel, we can see how the core of the more massive cloudlet (50 M_⊙) lags behind the shell and deforms the shell immediately around it. This implies that detached cloudlets cannot be accelerated to velocities beyond that of the shell purely due to stellar winds.

Deviations from a spherical shell have also arisen where stellar winds break through due to slight asymmetries, but these do not strongly affect the estimate of the shell radius.

Our simulation does not include an ambient ISM, which typically dampens the expansion of a stellar wind bubble as more ISM is swept up in the bubble’s outer shell (Weaver et al. 1977). Without this damping, the shell reached a terminal velocity of ~60 km s⁻¹ around 0.3 Myr.

Our multiphysics model converts gas particles into sink particles when they reach a density of 3.3 × 10⁶ amu cm⁻³ to prevent numerical errors. These sink particles interact with the gas component only through gravity, not through any hydrodynamical interaction. In the bottom panel of Fig. 11, we show all sink particles as white points. A number of sink particles form radially directed lines, with the most prominent trailing the 50 M_⊙ globule. This is reminiscent of the pillars observed at the interface between some H II regions and molecular shells, which Tremblin et al. (2012b) propose form when the curvature of the shell becomes high and the shell collapses in on itself. The globules we inserted into the initial conditions are natural places for such a curvature to develop. As the pillar collapses our model converts it into sink particles. In contrast with a sink particle, a real pillar is able to interact hydrodynamically with its surroundings, but the fact that they lag behind the shell implies that they decouple from the processes that drive the expansion and interact with their surroundings in a similar way to sink particles.

Tremblin et al. (2012a, 2013) also discuss another structure observed near such shells, namely dense globules just within the H II region. These structures form from interactions with turbulent material outside the shell. Because our simulations do not include turbulent gas these structures do not form.

Because stellar wind bubbles expand into the ISM supersonically, the ISM beyond a certain shell has not influenced its evolution up to that point. As a result, we could use the state of the simulation at any moment as the initial conditions of a new simulation with an arbitrary ISM distribution outside the shell. Re-simulating a large number of simulation snapshots with additional ISM outside the shell would be prohibitively expensive in terms of computer resources, because the entire ISM must be represented by SPH particles. The volume required to represent S 171 and to reduce edge effects (due to the absence of SPH particles beyond the simulation box) would be so great that the number of SPH particles required would be much greater than that of the shell.

Fig. 11 Gas column density maps of the numerical simulation at 0.01 Myr (top), 0.12 Myr (centre), and 0.30 Myr (bottom). The top panel also shows the projected positions of the stars more massive than 16 M⊙ (black points) and those less massive than 16 M⊙ (white points). The bottom panel shows the projected positions of the sink particles as white points. The arrows point out the 50 M⊙ (left) and 10 M⊙ (right) cloudlets. The centre and bottom panels show the position of the shell as a black dotted circle.

6.2 Simplified Numerical Wind Bubble Model

6.2.1 Model

Instead, we turn to a simplified (and computationally cheaper) numerical approach to a similar problem, which we can apply to the state of our simulation at different times. In our simulation, stellar winds quickly form a mostly spherical bubble in the gas cloud that is much larger than the collection of massive stars. We can approximate this as a wind bubble blown by a star with the stellar wind output of all massive stars combined (although losses in wind shocks between stars can decrease the effective input of energy and momentum into the shell).

The classic solution to this problem is the similarity solution of Weaver et al. (1977), which describes the spherical expansion of a stellar wind bubble in a uniform ISM. More recently, Lancaster et al. (2021a) presented a modified solution taking into account the fractal structure of the bubble’s boundary and the subsequent efficient radiative losses. We cannot use their solutions directly because our shell consists of a swept-up Plummer sphere, not a swept-up uniform ISM, but we can use the differential equations to which they are a solution.

Both models are based on the conservation of momentum, equating the momentum of the shell to that imparted by stellar winds. Although they differ in the prescription of how momentum is transferred to the shell, they are both of the following form: $$\frac{\text{d}}{\text{d}t}\left(M_{\text{S}}\left(R_{\text{S}}\right){\dot{R}}_{\text{S}}\right)=\eta {\dot{p}}_{\text{W}},$$(2)

where M_s (R_s) is the mass of the shell once it has reached radius R_s, ${\dot{p}}_{\text{w}}$ is the shell’s momentum input rate due to the stellar winds, and η is an efficiency factor. Lancaster et al. (2021a) introduce two different factors, α_p and α_R, that account for enhanced momentum transfer and non-spherical expansion, respectively. Both of these factors are variable but of order unity (as demonstrated by Lancaster et al. 2021b). For simplicity, we combine these into one value, which can additionally account for losses in momentum transfer due to stellar wind shocks between stars.

In a uniform ISM of density ρ₀, the mass of the shell is simply $M_{\text{s}}=\frac{4\pi }{3}{R}_{\text{s}}^3{\rho }_0$, resulting in a neat power law solution. In the case of our model, the density distribution follows a Plummer sphere interior to the initial shell radius, R_s,0. The density distribution outside this radius is undefined as of yet. For the sake of simplicity we assume a uniform ambient ISM density distribution. The shell mass is then $$M_{\text{s}}\left(R_{\text{s}}\right)=\{\begin{array}{l}{M}_{\text{P}}\left(<R_{\text{s}}\right)\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}if\hspace{0.17em}}{R}_{\text{s}}>R_{\text{s},0}\hfill \\ M_{\text{P}}\left(<R_{\text{s},0}\right)+\frac{4\pi }{3}\left(R_{\text{s}}^3-R_{\text{s},0}^3\right){\rho }_0\text{\hspace{0.17em}if\hspace{0.17em}}{R}_{\text{s}}>R_{\text{s},0}\hfill \end{array}$$(3)

where M_P (< R) is the mass of a Plummer sphere within a radius R and the ambient ISM density is ρ₀.

The resulting differential equations lack neat power law solutions, and require numerical integration. We use fourth-order Runge–Kutta integration for this, with time steps equal to 0.01 · min $\left(x_i/{\dot{x}}_i\right)$, where x_i denotes a vector of all integrated variables (e.g. shell radius and velocity).

Given a simulation snapshot, we can now continue the evolution of the wind bubble using this method. As initial conditions, we need the shell’s radius and velocity, and in the case of the Weaver model, the internal energy contained within the shell radius. We define the shell’s radius as the location of the peak of the histogram of all gas particle’s distance from the origin (and the shell’s velocity in a similar way). The mechanical luminosity of the central wind source is fixed by the stellar population. Finally, we set the ISM density exterior to the shell to be equal to the density of the main Plummer sphere at the shell radius. This ensures continuity in the density distribution, and uniquely relates each shell radius to an ISM density.

If the evolution of the bubble using the simplified numerical model reaches the observed shell radius and velocity at the same moment in time, the combined model is consistent with the observations. This moment in time is then our estimate of the system age, and the density of the main Plummer sphere at the initial shell radius is our estimate of the ambient ISM density. Because we have a limited set of initial conditions, the shell radius and velocity will never match the observed values at the exact same moments. However, for every set of initial conditions (each corresponding to an ISM density), we can record the times at which the observed shell radius and velocity are reached, and find the intercept of these two curves.

6.2.2 Results

Figure 12 shows the evolution through time of the shell radius (orange) and velocity (blue) following the Weaver et al. (1977, solid) and Lancaster et al. (2021a) (dashed) models. We note that all times quoted in this section are with respect to the start of the multiphysics simulation. For comparison, we have also plotted the analytic solutions of the Weaver (dotted) and Lancaster (dash-dotted) models using the same stellar wind output and ISM density. The initial conditions of the numerical solutions are offset from these solutions, but the numerical solutions do converge to the analytic solutions. This is more obvious for the Weaver solution than the Lancaster solution; in the latter case, the shell radius is initially too large, and the velocity must be below the analytic shell velocity so the radius can eventually converge to the analytic shell radius. Further evolution of the system does show that the numerical solution converges to the analytic solution.

The shaded regions denote the radius and velocity of the massive, moderate velocity shell component (25 pc and 4 km s⁻¹, respectively) with 10% errors. For the initial conditions used for this plot, the combined multiphysics and Lancaster models agree with the observations because both shell radius and velocity cross the observed values at the same moment, at about 2 Myr. Notably, the velocity of the Weaver solution is much greater than observed, and does not reach the observed value even at 10 Myr. For this instance of the simplified model, we used the state of the multiphysics model at 0.23 Myr, when the shell radius was 6.5 pc, the shell velocity 41 km s⁻¹, and the ISM density at the shell edge was 19.4 amu cm⁻³.

Figure 13 shows, for a number of simulation snapshots, the density of a Plummer sphere at the shell radius and the moments that the numerically integrated (using the Lancaster model for the top panel, and the Weaver model for the bottom panel) shell radius and velocity reach values of 25 pc and 4 km s⁻¹ (corresponding to the massive, moderate velocity shell). The shaded regions denote the range of values allowed by the efficiencies shown and by letting the observed radius and velocity vary by 10%. We point out that the two panels cover different ranges of efficiency η. This is expected to be of order unity, but in the case of the Weaver model there were no solutions with η = 1 or η = 3 across the available range of densities. The radius and velocity curves did approach each other at high densities, but did not cross at the greatest density possible, which corresponds to the density of the Plummer sphere’s core (661 atomic mass units (amu) cm⁻³).

From the Lancaster model, the ISM density inferred is in the range of 10–30 amu cm⁻³, slightly below the typical density of molecular clouds. The age is in the range of 2–2.5 Myr, similar the age assumed for the system.

From the Weaver model, the ISM density inferred is in the range 30–300 amu cm⁻³, of the order of a typical molecular cloud. The age is in the range 2.5–4 Myr, slightly older than estimated earlier. It should be noted that the lower ISM density and age correspond to a very low stellar wind efficiency, one that we have not explored with the Lancaster model.

We do not expect such a low efficiency to be plausible; Lancaster et al. (2021b) only finds an efficiency of similarly low order of magnitude for a very dense, very low star formation efficiency cloud. In their work, they use the cloud mass and star formation efficiency to compute an effective stellar wind mechanical luminosity. Inverting this, with our actual cloud mass and mechanical luminosity, we get an effective star formation efficiency of ~50%, which implies a stellar wind efficiency close to unity.

For comparison, we also present the results of matching the observed radii and velocities to the analytic solutions of Weaver and Lancaster. These results are in Table 8. We match both the massive, moderate velocity shell and the high-velocity cloudlets. This disregards the presence of leftover cluster material, but is independent of potential numerical errors in our simulation.

Fitting the Weaver model to the moderate velocity shell results in a relatively old age and a very dense ISM (about twice the density of the main Plummer sphere in our simulation), unless η is very low. Fitting the Lancaster model, on the other hand, leads to a somewhat younger system age (though still much older than estimated for the system), and an ISM more diffuse than a GMC. For completeness, we also included the ages and densities resulting from fitting to the high-velocity cloudlets, although it is difficult to explain them as a shell when the moderate velocity shell is much more massive and closer in.

We did not observe a clear velocity bimodality in the simulation. The amount of material with velocities greater than 2.5 times that of the shell itself (which would be an analogue to the observed high-velocity structures) is only ~0.01% of the total shell material, or ~0.3 M_⊙, less massive than the individual high-velocity cloudlets. Additionally, as can be seen in Fig. 11, the cloudlets that were detached from the main cloud at initialisation typically lagged behind the shell and have similar velocities to the shell itself, and so did not end up as a high-velocity component.

Fig. 12 Evolution through time of the stellar wind bubble radius (orange) and velocity (blue) using the Weaver et al. (1977) (solid) and Lancaster et al. (2021a) (dashed) models, starting from the bubble’s radius and velocity at 0.23 Myr (with a corresponding density at the shell of l6.0 amu cm−3). The efficiency α parameters in the Lancaster model have both been set to 1. The dotted and dash-dotted lines denote the analytical Weaver and Lancaster solutions, respectively. The shaded regions show 25 pc and 4 km s−1, with 10% errors.

Fig. 13 For a number of snapshots of the simulation, the density of a Plummer sphere at the shell position (assumed to be the ISM density) versus the moment at which in the numerical integration a shell velocity of 4 km s−1 and a radius of 25 pc is reached. Different symbols denote different stellar wind efficiencies η. Shaded regions denote the envelope allowed for all efficiencies plotted and a 10% error in the observed velocity and radius. The top panel shows the results obtained using the Lancaster model, and the bottom panel shows the results from the Weaver model.

6.3 Simplified Numerical Cloudlet Model

6.3.1 Model

Lancaster et al. (2021b) show that stellar wind bubbles blown in a uniform but turbulent medium are fractal shaped rather than uniform spherical shells. Even in our smooth² and virialised Plummer sphere asymmetries arise, from secondary bubbles where the winds break through the shell into low-density material (shown in the top left and bottom right of the bottom frame of Fig. 11) to substructures in the shell itself.

In the models above, the fundamental mechanism slowing down the expanding bubble is the accumulation of ISM material in the shell. The magnitude of this deceleration is related to the ratio of the masses of the shell and of the swept-up ISM. In a non-uniform shell, different sections will have different ratios, and as a result, different rates of deceleration. Viewed differently, a denser section of the shell will have more inertia and be able to maintain its velocity for longer than a less dense section. In this section we investigate whether this varying deceleration can result in the multiple components observed in S 171.

We derive the equation of motion of a section of a shell, or cloudlet (which we use from now on), through a uniform medium, using a similar method to Lancaster et al. (2021a), that is, exploiting conservation of momentum.

This momentum equation is of the following form: $$\begin{array}{l}\frac{\text{d}}{\text{d}t}\left(\left(M_{\text{c,0}}+{\rho }_0{\displaystyle {\int }_{t_0}^t{\sigma }_{\text{c}}}\left(t^{\prime }\right){\dot{R}}_{\text{c}}\left(t^{\prime }\right)\text{d}{t}^{\prime }\right){\dot{R}}_{\text{c}}\left(t\right)\right)\hfill \\ =\eta {\dot{p}}_{\text{W}}\frac{{\sigma }_{\text{c}}\left(t\right)}{4\pi R_{\text{c}}{\left(t\right)}^2},\hfill \end{array}$$(4)

where M_c and R_c are the cloudlet mass and distance to the cloud centre, respectively, and σ_c is the cloudlet’s cross-section as seen from the cloud centre. This is both the surface being struck by stellar winds, and the surface accumulating ISM material. Also, in this expression, and all hereafter, a subscript 0 denotes the initial value of a quantity.

The evolution of σ_c is dictated by hydrodynamical processes, but for simplicity we consider two interesting edge cases: constant cross-section (corresponding to a bound, non-expanding cloudlet), and quadratically growing with distance (corresponding to a shell structure growing with the shell). These have the respective forms $${\sigma }_{\text{c}}\left(t\right)={\sigma }_{\text{c,0}},$$(5) $${\sigma }_{\text{c}}\left(t\right)={\sigma }_{\text{c,0}}{\left(\frac{R_{\text{c}}\left(t\right)}{R_{\text{c,0}}}\right)}^2.$$(6)

They reduce to the following differential equations in R_c, for the non-expanding and expanding cases, respectively $${\ddot{R}}_{\text{c}}=\left(\frac{\eta {\dot{p}}_{\text{ω}}}{4\pi R_{\text{c}}^2}-{\rho }_0{\dot{R}}_{\text{c}}^2\right)/\left(\frac{M_{\text{c,0}}}{{\sigma }_{\text{c,0}}}+{\rho }_0\left(R_{\text{c}}-R_{\text{c,0}}\right)\right)$$(7) $$\begin{array}{lll}{\ddot{R}}_{\text{c}}\hfill & =\hfill & \frac{\eta {\dot{p}}_{\omega }}{4\pi R_{\text{c,0}}^2}-{\rho }_0{\dot{R}}_{\text{c}}^2{\left(\frac{R_{\text{c}}}{R_{\text{c,0}}}\right)}^2/\hfill \\ \hfill & \hfill & \left(\frac{M_{\text{c,0}}}{{\sigma }_{\text{c,0}}}+\frac{\rho 0}{3R_{\text{c,0}}^2}\left(R_{\text{c}}^3-R_{\text{c,0}}^3\right)\right).\hfill \end{array}$$(8)

We note thatthe cloudlet’s initial mass and cross-section only appear as the ratio $\frac{M_{\text{c,0}}}{{\sigma }_{\text{c,0}}}$. This implies that a cloudlet’s evolution is determined by its initial average (mass) surface density, and its actual dimensions are not important.

Another point to note is that the expanding case is equivalent to the Lancaster model when the initial average surface density of the cloudlet is the same as that of the shell.

Fig. 14 Distance to the centre of the complex (orange) and velocity (blue) of cloudlets as a function of initial average surface density, after evolution to 1.9, 2.17, and 2.3 Myr (for η = 0.3, 1, and 3, respectively). Solid lines represent cloudlets that do not expand as they travel, i.e. bound structures. Dashed lines represent cloudlets that expand with the shell, with their cross-sections growing with the distance squared, i.e. features on the shell that grow with it. Symbols overplotted on these lines indicate the stellar wind efficiency η. The hatched shaded regions denote 31 pc and 12 km s−1, the radius and velocity of the high velocity component of S 171, and the clear shaded regions denote 25 pc and 4 km s−1, the radius and velocity of the moderate velocity component. All observed values are shown with 10% errors. The surface densities and line-of-sight velocities (with respect to the cluster centre) of the cloudlets in Table 5 are shown as upward blue arrows. This velocity is a lower limit of their expansion velocity. The initial conditions were taken from the simulation snapshot at 0.22, 0.23, and 0.25 Myr (again for η = 0.3, 1, and 3, respectively), which resulted in a bubble radius of 25 pc and a velocity of 4 km s−1 at the respective end time using the Lancaster model. The two vertical black lines denote the average surface density of the shell at 0.23 Myr (left) and the surface density at the core of a 50 M⊙ (and 0.5 pc) Plummer sphere cloudlet.

6.3.2 Results

Figure 14 shows the distance and velocity of cloudlets of varying surface density, using the differential equations described above. The initial conditions for each efficiency were taken from the snapshots best coinciding with the observations of the shell after continued integration, and evolved until the moment the shell coincided with the observations.

Non-expanding cloudlets always reach farther distances and retain higher velocities than expanding cloudlets because they sweep up less of the ISM. This effect is most pronounced near the centre of the range of densities explored, where non-expanding cloudlets can reach distances about twice as far as expanding cloudlets and velocities about four times as fast. Compared to the average shell, non-expanding cloudlets with high surface densities can reach distances almost four times the shell radius and velocities about ten times as great.

We note that there is no density for which the distance and velocity match the observed high-velocity cloudlets. The configuration closest to matching the observations is a non-expanding cloudlet of about the same initial average surface density as the shell. However, this cloudlet’s distance is more than 10% greater than the observed distance, and its velocity is more than 10% smaller than observed. Allowing a cloudlet to expand more slowly than quadratic would put it in between the expanding and non-expanding cases; this would alleviate disagreement with the distance, but increase disagreement with the velocity. Allowing the stellar wind efficiency to vary cannot alleviate this disagreement either, for similar reasons. However, assuming a slightly younger system age (at which the radius was smaller, and the velocity greater) could put the cloudlets in the observed range.

We also show the high-velocity cloudlets. We derive their average surface densities from Table 5 using a distance of 1.1 kpc, and correct the line-of-sight velocity for that of the cluster itself (−9.5 km s⁻¹). This is a lower limit on their expansion velocity, which can also have a component in the plane of the sky. In our spherical expansion scenario, the line-of-sight velocity is closer to the expansion velocity the closer a cloudlet is to the cluster centre. The cloudlets C:C and N2A are very close to the centre, while EnorthB, EnorthC, W10A, and W10B are nearly at the 30 pc outer circle of the projected shell. In terms of the surface density, the cloudlets are in the range of 20–200 M_⊙ pc⁻². This combination of velocity and surface density puts the cloudlets closer to the expanding limit than the non-expanding limit.

6.4 Summary

To summarise, through a hybrid numerical modelling approach we demonstrate that stellar winds from the massive stars of Be 59 can give rise to the complex gas structure of S 171. We used a coupled hydrodynamical, gravitational, and stellar evolutionary model to simulate the initial formation of a stellar wind bubble in the concentrated gas distribution left over from the formation of a star cluster. We then used simplified numerical models to continue the evolution of the wind bubble through uniform ISM of different densities in an effort to match the observed radius and velocity of the massive shell of S 171, and so provide an estimate of the system age and ambient ISM density. Finally, we introduced a new simplified numerical model to describe the dynamic evolution of shell substructures through a uniform ISM, and used it to show that dense and detached shell fragments are able to reach farther distances and higher velocities than the main shell.

We show that using the recent stellar wind bubble model by Lancaster et al. (2021a) to continue the bubble evolution results in more reasonable values for the system age and ambient ISM density compared to the classical model of Weaver et al. (1977), assuming reasonable values for the efficiency of momentum transfer. As derived from the Lancaster model, the system age is consistent with earlier estimates, and the ambient ISM density is reasonable for the larger surroundings of a star forming regions. The Weaver model requires very inefficient momentum transfer from stellar winds to the shell to reproduce the moderate velocity shell, and then results in a system age older than estimated earlier and a very dense ambient ISM.

We also show that variations in density within the shell can lead to components propagating through the ISM at different velocities. We are unable to reproduce the observed ratios between the distances and velocities of the observed structures, but the case of non-expanding cloudlets with initial surface density similar to the shell comes closest. These variations can arise naturally from fragmentation of the shell, as is demonstrated in our multiphysics simulation. If the shell is only able to partially fragment into small cloudlets we would naturally obtain highvelocity, low-mass cloudlets and a moderate-velocity, high-mass shell.

Our results qualitatively show that stellar winds can reproduce the gas structure of S 171, but we refrain from drawing strong quantitative conclusions. The SPH method of modelling hydrodynamics is infamously bad at modelling shocks, being unable to properly resolve thin shock layers. This makes our estimates of the shell radius and velocity that serve as the initial conditions of the simplified numerical models uncertain. Notably, the shell velocity is greater than expected from the Weaver and Lancaster models, even though the shell mass is greater due to the swept-up Plummer sphere. Although our model slightly favours the Lancaster model over the Weaver model, we also refrain from drawing conclusions about this.

A fully consistent reproduction of the system would require a grid based hydrodynamical model (such as the one used by Lancaster et al. 2021b) or a hybrid method (such as moving mesh, or meshless finite mass). Models that start from the collapse of a giant molecular cloud (such as by Wall et al. 2020; Grudić et al. 2021) can improve the authenticity of the gas distribution at the moment the bubble begins being blown. A full hydrodynamical simulation, including an ambient ISM, can also self-consistently demonstrate the origin of the dense cloudlets that end up at high velocity.

7 Conclusions

We have investigated the kinematics of an expanding H II region powered by a young stellar cluster for which the inferred expansion velocity appears to be higher than predicted in theoretical model simulations. However, published estimates of the central velocity of this complex can be highly divergent, and more detailed information on the motion of the surrounding shell is warranted. In this paper we draw attention to the nebula Sharpless 171 (with NGC 7822), which surrounds the cluster Berkeley 59. Optical spectra of 27 stars were obtained in the cluster area, providing a new estimate of the RV of the cluster. The velocity pattern over the molecular shell was mapped in ¹³CO(1–0) at OSO. From these data we deduce the expansion velocities of different cloud fragments in the shell relative to the cluster. As described below, we could define three systems of such clouds expanding at different rates. One system expands at a high rate, although not as high as one could expect from some of the published estimates of the central velocity.

We assigned a spectral class to each star, including the separate components in double-lined spectroscopic binaries. The majority of the stars had not been classified before, including one central O star hidden behind foreground dust.

From spectral classifications and RVs, plus existing photometric data and Gaia parallaxes, we selected 19 of the 27 stars as likely members of Be 59.

Repeated observations were done, especially for members, which revealed that five of the massive stars are spectroscopic binaries, of which one is double-lined and one triple-lined. We followed these stars over longer times and derived periods, semi-amplitudes, and systemic velocities. The most luminous O star is the eclipsing binary V747 Cep, and our RV curve agrees with the photometric period.

From the RVs of members we obtained a mean heliocentric velocity of − 18.75 km s⁻¹ (−9.5 km s⁻¹ in lsr), which we define as the central RV of the complex. None of the members deviate significantly from the mean value.

The S 171 complex extends over 3.2° in the sky, and for our ¹³CO mapping we primarily covered areas with distinct dark clouds. The velocity pattern over the nebula is complex, and in most directions several velocity components are present in the spectra. In addition to components that originate in gas not related to the complex, we identified two systems of shell structures expanding at different velocities. These components can be described as located at the peripheries of spherically expanding shells. One system, composed primarily of rather massive shell fragments, expands at a moderate velocity, 4 km s⁻¹, and includes two velocity components that are related to the front and back sides of the complex. This system is confined within a radius of 25 pc. The second system expands at a much higher velocity, 12 km s⁻¹, and is concentrated in a belt that crosses the entire H II region, extending to a radius ≥30 pc. In addition, we find some cloudlets in the central area with expansion velocities in between those of these systems.

The luminosity and UV flux is dominated by the central O stars. Our interpretation of the observed complex kinematic structure is based on the assumption that the expansion started about 2 million years ago, driven by radiation and winds from the central massive stars. Our model simulations show that cloudlets of smaller masses obtained very high initial velocities, while more massive shell fragments dragged behind.

The velocity and radius of the moderate velocity shell are generally consistent with theoretical predictions for the expansion of a stellar wind bubble driven by the cluster’s population of massive stars.

The cloudlets observed with higher velocity and at greater distance from the massive shell can arise through interaction with the ambient ISM, which decelerates dense and bound substructures of the shell less efficiently than low-density substructures that co-expand with the shell.

Although our models are able to qualitatively reproduce the structure of S 171 for plausible parameters, we refrain from drawing quantitative conclusions about the value of these parameters. This requires further work using more accurate methods.

Acknowledgements

We thank the anonymous referee for helpful comments that led to an improvement of the paper. This work was supported by the Magnus Bergvall Foundation and the Längmanska Kulturfonden in Sweden. We made use of the SIMBAD database operated at Centre de Donnees Astronomique de Strasbourg and data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This publication made use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. We also used a map obtained with Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. We made use of the Binary Star Solver v1.2.0 python program (Milson et al. 2020). In this work, we use the matplotlib (Hunter 2007), NumPy (Harris et al. 2020), AMUSE (Portegies Zwart & McMillan 2018), SeBa (Portegies Zwart & Verbunt 2012), Huayno (Pelupessy et al. 2012), and Fi (Pelupessy et al. 2004) packages. We thank John Telting for very useful help, comments and suggestions related to FIES, and in addition we thank the following observers at the NOT for contributing to obtaining the FIES data set: Pere Blay Serrano, Rosa Clavero, Abel de Burgos, Sune Dyrbye, Grigori Fedorets, Emanuel Gafton, Francisco Jose Galindo-Guil, Nicholas Emborg Jahnsen, Raine Karjalainen, Anni Kasikov, Marcelo Aron Keniger, Jyri Lehtinen, Teet Kuutma, Illa Losada, Julia Martikainen, Shane Moran, Christina Konstantopoulou, Viktoria Pinter, Tapio Pursimo, Rene Tronsgaard Rasmussen, Lauri Siltala, Auni Somero, and Joonas Viuho. Based on observations made with the Nordic Optical Telescope, owned in collaboration by the University of Turku and Aarhus University, and operated jointly by Aarhus University, the University of Turku and the University of Oslo, representing Denmark, Finland and Norway, the University of Iceland and Stockholm University at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofisica de Canarias.

Appendix A Individual stellar radial velocities

Fig. A.1 He I 5875 Å line at 12 epochs for V747 Cep (left). Note the inverse P Cygni profile in the top spectrum.

Fig. A.2 He I 5875 Å line at 13 epochs for the spectroscopic triple BD+66° 1675 (right).

Appendix B Positions observed for ¹³CO

Appendix C Representative ¹³CO spectra for each region

Fig. C.1 Selected 13CO spectra from each area. Components associated with the expanding shell are typically at < −10 km s−1, components related to the local gas (Gould’s Belt) are found around 0 km s−1, and those related to the Orion arm at around −6 km s−1. Offsets are W10:A (0,0), W9 (+240,0), W8 (0, +120), W7 (0,0), W6:A (+240,0), W5 (0,0), W4 (0,+60), W3:A (+120,+540), W2 (0,−60), W1 (−150,+180), C:C (−90,−60), N3:B (0,0), N2:A (0,0), N1:B (0,0), S1 (0,−60), S2:A (+180,0), S3:A (0,0), and S4 (0,0).

Fig. C.2 As in Fig. C.1. Offsets are E1:B (0,0), E2:B (+300,+60), E3:A (−540,−480), Esouth:E (0,0), Enorth:B ((0,0), and E4:B (−360,−360).

References

All Tables

All Figures

Fig. 1 Optical image from the Digitized Sky Survey (left) and mid-IR image from Wide-field Infrared Survey Explorer (WISE, right), featuring the H II region S 171 in relation to the cluster Berkeley 59, which is located in the central white box. The northern part of the nebula is called NGC 7822. The large circle has a radius of 1.6° and shows the approximate extent of the H II region. The green contours show the 4 mm (70 GHz) continuum map from Planck. The smaller circles in the optical image mark the positions of the hottest members of the extended Cep OB 4 association; light blue circles represent spectral types O9 and B0, and yellow circles types B1 and B2. North is up, and east to the left. The WISE image is a colour composition of red (22 µm), green (4.6 µm), and blue (3.4 µm). Figure 2 is a zoomed-in view of the area in the white box.

In the text

	Fig. 2 Locations of program stars in the Berkeley 59 region according to the designations in Table 3 on an image from the Digitized Sky Survey. The cluster centre is put between the bright stars A (BD+66° 1674) and B (BD+66° 1675). The white arrow points at BD+66° 1673, V747 Cep. The red arrow marks the position of a bright IR source. North is up, and east is to the left.
In the text

Fig. 3 Examples of normalised spectra obtained for some program stars over two selected wavelength regions. Certain spectral features discussed in the text are marked. The spectra are ordered in steps of 0.3 (intensity scale) in the left panel and 1.0 in the right panel and with a smaller spacing between the third and fourth spectra; the double-lined spectroscopic binary BD+66° 1674 is shown at two different phases. To improve the S/N for the faint IR star, its left spectrum is a combination of many observations, Doppler-corrected for the binary motion, which smears out the DIBs.

In the text

	Fig. 4 Phase diagrams for (from the top): V747 Cep for a period of 5.33146 days; BD+66° 1675, period of 74 days; the IR star, period of 9.51 days; and P16, period of 13 days.
In the text

Fig. 5 Map of the fields we obtained CO spectra in. Left: large region, 4.4 × 4.4 degrees, surrounding the emission nebula S 171. The arrow points at the position of the central cluster. The northern arc of bright nebulosity is NGC 7822. In this panel the more peripheral areas observed for 13CO emission are marked and labelled. The circles (Enorth and Esouth) mark observations at selected positions along two dusty striations that extend to the east from the centre. The dashed box marks the central region shown in the right panel. Right: image spanning 1.9 × 1.2 degrees. The core of the stellar cluster is inside the dashed box. The position of an elephant trunk, called the ‘Dancing Queen’ is marked with a black circle. North is up and east to the left in both images, which are from the red Deep Sky Survey.

In the text

Fig. 6 Part of the map of 13CO spectra obtained for area W8. The offsets in RA and Dec from the central position are in steps of 60″. The velocity interval, in km s−1, and the antenna temperature, TΑ in Kelvin, are indicated in the upper-right panel. The component at around −18 km s−1, present in most panels, is related to the foreground shell, whereas the component at around −7 km s−1 has a local origin.

In the text

	Fig. 7 13CO spectra from one position in area E2. The strongest signal at around −13 km s−1 comes from the foreground shell, while the components at −6 and +2 km s−1 are related to gas connected to the Orion arm and Gould’s Belt. The component at around −2 km s−1 may trace the remote part of the molecular shell.
In the text

	Fig. 8 Positions in galactic coordinates of molecular cloudlets moving at high negative speed (≤−12 km s−1) in an extended area outside the central part of the S 171 complex (left image from Fig. 5 inserted for reference). Open squares are cloudlets from Yonekura et al. (1997) and filled squares are E4 and Esouth:E from our survey.
In the text

Fig. 9 Radial velocities relative to the cluster velocity for different sub-areas against angular distance from the centre (in all directions). Blue dots mark areas that are expanding at high velocity, vexp ≈12 km s−1, and yellow circles (smaller) areas expanding at moderate speed, vexp ≈ 4 km s−1, when assuming spherically expanding systems (dotted curves). Green circles mark areas with intermediate velocities (‘DQ’ is the Dancing Queen). Red squares show areas with components whose velocities are redshifted relative to the cluster.

In the text

Fig. 10 Map of cloudlets. Left: overview of the region shown in Fig. 1 in 12 µm from WISE (grey scale) with the Planck 857 GHz map (black contours) overlaid to emphasise the dust emission. A black circle spanning a radius of 1.6° is centred on the Berkeley 59 cluster. The blue dots from Fig. 9 are shown as blue circles and mark the locations where the high negative expansion velocities are found. Likewise, the green circles refer to the intermediate and the yellow the low negative velocity components, while the red circles mark positions where the positive component is found, possibly related to the expansion on the far side. Right: close-up 55′× 40′ view of the central region where the red, green, blue coding is WISE 12 µm (red), 4.6 µm (green), and Digitized Sky Survey 2 red plates (blue). The positions of V747 Cep and the IR star are marked with white diamonds.

In the text

Fig. 11 Gas column density maps of the numerical simulation at 0.01 Myr (top), 0.12 Myr (centre), and 0.30 Myr (bottom). The top panel also shows the projected positions of the stars more massive than 16 M⊙ (black points) and those less massive than 16 M⊙ (white points). The bottom panel shows the projected positions of the sink particles as white points. The arrows point out the 50 M⊙ (left) and 10 M⊙ (right) cloudlets. The centre and bottom panels show the position of the shell as a black dotted circle.

In the text

Fig. 12 Evolution through time of the stellar wind bubble radius (orange) and velocity (blue) using the Weaver et al. (1977) (solid) and Lancaster et al. (2021a) (dashed) models, starting from the bubble’s radius and velocity at 0.23 Myr (with a corresponding density at the shell of l6.0 amu cm−3). The efficiency α parameters in the Lancaster model have both been set to 1. The dotted and dash-dotted lines denote the analytical Weaver and Lancaster solutions, respectively. The shaded regions show 25 pc and 4 km s−1, with 10% errors.

In the text

Fig. 13 For a number of snapshots of the simulation, the density of a Plummer sphere at the shell position (assumed to be the ISM density) versus the moment at which in the numerical integration a shell velocity of 4 km s−1 and a radius of 25 pc is reached. Different symbols denote different stellar wind efficiencies η. Shaded regions denote the envelope allowed for all efficiencies plotted and a 10% error in the observed velocity and radius. The top panel shows the results obtained using the Lancaster model, and the bottom panel shows the results from the Weaver model.

In the text

Fig. 14 Distance to the centre of the complex (orange) and velocity (blue) of cloudlets as a function of initial average surface density, after evolution to 1.9, 2.17, and 2.3 Myr (for η = 0.3, 1, and 3, respectively). Solid lines represent cloudlets that do not expand as they travel, i.e. bound structures. Dashed lines represent cloudlets that expand with the shell, with their cross-sections growing with the distance squared, i.e. features on the shell that grow with it. Symbols overplotted on these lines indicate the stellar wind efficiency η. The hatched shaded regions denote 31 pc and 12 km s−1, the radius and velocity of the high velocity component of S 171, and the clear shaded regions denote 25 pc and 4 km s−1, the radius and velocity of the moderate velocity component. All observed values are shown with 10% errors. The surface densities and line-of-sight velocities (with respect to the cluster centre) of the cloudlets in Table 5 are shown as upward blue arrows. This velocity is a lower limit of their expansion velocity. The initial conditions were taken from the simulation snapshot at 0.22, 0.23, and 0.25 Myr (again for η = 0.3, 1, and 3, respectively), which resulted in a bubble radius of 25 pc and a velocity of 4 km s−1 at the respective end time using the Lancaster model. The two vertical black lines denote the average surface density of the shell at 0.23 Myr (left) and the surface density at the core of a 50 M⊙ (and 0.5 pc) Plummer sphere cloudlet.

In the text

	Fig. A.1 He I 5875 Å line at 12 epochs for V747 Cep (left). Note the inverse P Cygni profile in the top spectrum.
In the text

	Fig. A.2 He I 5875 Å line at 13 epochs for the spectroscopic triple BD+66° 1675 (right).
In the text

Fig. C.1 Selected 13CO spectra from each area. Components associated with the expanding shell are typically at < −10 km s−1, components related to the local gas (Gould’s Belt) are found around 0 km s−1, and those related to the Orion arm at around −6 km s−1. Offsets are W10:A (0,0), W9 (+240,0), W8 (0, +120), W7 (0,0), W6:A (+240,0), W5 (0,0), W4 (0,+60), W3:A (+120,+540), W2 (0,−60), W1 (−150,+180), C:C (−90,−60), N3:B (0,0), N2:A (0,0), N1:B (0,0), S1 (0,−60), S2:A (+180,0), S3:A (0,0), and S4 (0,0).

In the text

	Fig. C.2 As in Fig. C.1. Offsets are E1:B (0,0), E2:B (+300,+60), E3:A (−540,−480), Esouth:E (0,0), Enorth:B ((0,0), and E4:B (−360,−360).
In the text