© ESO 2018

1 Introduction and motivation

According to a generally accepted paradigm, the structures known as superbubbles (Chu 2008; McClure-Griffiths 2012; Ambrocio-Cruz et al. 2016) are the result of the interaction of multiple stellar winds within stellar clusters. Individual examples are the central cluster at the Galactic center (Ozernoy et al. 1997), the Arches cluster (Cantó et al. 2000; Raga et al. 2001), the Carina complex (Gull 2011), M 17 (Reyes-Iturbide et al. 2009; Mernier & Rauw 2013), and N 70 (Rodríguez-González et al. 2011).

The interest in superbubbles is triggered by the desire to understand their X-ray luminosities (e.g., Raga et al. 2001), by their potential contributions to the Galactic cosmic-ray flux (e.g., Binns et al. 2007; Murphy et al. 2016) via acceleration at multiple shocks (Bykov 2001; Ferrand & Marcowith 2010), and because they are sources of gamma rays (Cesarsky & Montmerle 1983; Domingo-Santamaría & Torres 2006).

In many cases the resulting large-scale wind outside a stellar cluster is assumed to be a spherically symmetric outflow, similar to that described by the original model of interstellar bubbles by Castor et al. (1975) and Weaver et al. (1977). This assumption is invalid inside a given cluster where the winds of several (probably dominating) O and B stars interact with each other. The mutual interaction of stellar winds has been studied over many years in great detail for the case of binary stars (e.g., Kallrath 1991; Pittard 2011; Reitberger et al. 2014, 2017, and references therein), including their additional interaction with the interstellar medium (Wareing et al. 2007; Banda-Barragán et al. 2016, 2018). Additionally, Parkin et al. (2014) and Reitberger et al. (2014), among others, have studied the effects of orbital motion in binary systems. These studies are neglected here, however, because we consider widely separated (~1 pc) stars for which the orbital timescale is much longer (~10⁷ yr) than the time taken for the interaction region of the winds to come into equilibrium (~10⁵ yr).

Compared to the existing body of research on binary winds, the investigations of the simultaneous interaction of multiple stars are still in their infancy. The first quantitative model of a large-scale cluster wind resulting from a more realistic internal flow structure was presented by Cantó et al. (2000). For numerical purposes, these authors had to assume a symmetric distribution of wind-driving stars inside the cluster. This model was subsequently generalized to a statistically homogeneous (Raga et al. 2001) and to a nonuniform (Rodríguez-González et al. 2007) stellar distribution. These hydrodynamical (HD) simulations were later applied to explain the X-ray emission of various clusters (as mentioned above) and further generalized to include stellar winds and a supernova event as well (Rodríguez-González et al. 2011).

All of the models employing a nontrivial internal structure of stellar clusters are of numerical nature and do not provide details of the colliding plasma flows inside a given cluster. The present paper intends to fill these modeling gaps. First, we present a simplified two-dimensional (2D) analytical model of the interaction of multiple stellar winds and the resulting large-scale cluster wind. Second, we analyze the principal structure of such multisource wind flow geometries. Third, we discuss the corresponding generalization to explore asymmetric 3D configurations, sometimes resorting to numerical methods. In doing so, we restrict ourselves to stellar clusters where the mutual distances between the dominating wind-driving stars are large enough to justify neglect of any orbital motion. Since the termination shock surfaces of the stars are disjoint, the entire wind interaction region is subsonic and thus, to a good approximation, can be taken as incompressible. As shown by several authors (Arthur 2007; Mackey et al. 2015; Scherer et al. 2016), at least for isolated astrospheres, the density inside an astropause is so low that it is not affected by cooling. The region between astropause and bow shock is affected, but this is not relevant to our study because there are no bow shocks, but rather flows of the intracluster medium around the astropauses. For these, we assume that this medium is in thermal equilibrium, i.e., that cooling is not effective, and thus it can be neglected as well.

2 Planar two-dimensional analytical model

The mutual interaction of stellar winds is best and often modeled within the framework of either hydrodynamics (HD) or magnetohydrodynamics (MHD; see, e.g., Kissmann et al. 2016). The purely hydrodynamical description of the subsonic interaction region allows an analytical treatment based on the theory of velocity or momentum density potentials. This concept has been introduced for the interaction of a stellar wind with the so-called interstellar wind (blowing in the rest frame of a given wind-driving star as a consequence of the relative motion between star and surrounding interstellar medium) by Parker (1961; see also Nerney & Suess 1995), and has been used recently to compute the distortion ofthe local interstellar magnetic field due to the solar wind plasma bubble (Röken et al. 2015; Kleimann et al. 2016, 2017).

In the following, we extend the modeling based on velocity potentials to the case of multiple interacting stellar winds. In order to explore the principal structure of such multisource winds, we first analyze the planar 2D case because it allows the definition of a single stream function, and thus an easy analytical representation of the stream lines of the interacting winds. While a stream function is, in general, not available in the 3D case (see, e.g., Elshabka & Chung 1999; Lee 2009), velocity potentials can be defined without problem and permit a semi-analytical study, at least.

2.1 Velocity potential and stream function of a multisource flow

The use of velocity potentials is based on the assumption (1)

according to which flow velocity u =u(r) at is irrotational. It has been demonstrated for the example of the interaction of the solar with the interstellar wind that this is a very reasonable approximation for steady-state flow configurations (e.g., Röken et al. 2015, and references therein). For a generalization to a momentum density potential see Kleimann et al. (2017). The above Eq. (1) implies the existence of a velocity potential Φ = Φ(r) with (2)

in which we adhere to the convention of omitting the minus sign. For the case of a single point source of strength m located at r_⋆, i.e., an isolated, noninteracting stellar wind, the velocity potential and the associated stream function Ψ = Ψ(r) are given by

(Landau & Lifshitz 1966). The neglect of the third Cartesian coordinate indicates that we limit our consideration in this section to the 2D planar case. The flow velocity can be computed from either Φ or Ψ:

This dependence on the stream function implies (7)

and thus that the flow is incompressible and obeys the homogeneous mass continuity equation. The significance of the functions Ψ and Φ lies in the fact that lines of constant Ψ represent stream lines to which lines of constant Φ are perpendicular, such that these lines span the grid of a curvilinear, orthogonal, flow-aligned coordinate system.

Due to the linearity of Eqs. (5) and (6), the flow from N different (point) sources of strength m_k located at r_k = (x_k, y_k) is obtained from a superposition

of the respective individual potentials and stream functions as

2.2 Complex potential

Upon noting from Eqs. (5) and (6) that Φ and Ψ fulfill the Cauchy–Riemann differential equations, one can define a complex potential at via (12)

and canthus combine Eqs. (8) and (9) into one, (13)

with ℜ(⋅) and ℑ(⋅) denoting the real and imaginary part of their complex argument. The analogous definition of a complex velocity (14)

(note the minus sign!) enables us to further combine Eqs. (10) and (11) into (15)

with a bar indicating the complex conjugate.

A crucial consequence of choosing Ψ according to Eq. (4) is that this choice implies (16)

for any , and thus prevents us from distinguishing pairs of stream lines on opposite sides of the source, since these get mapped to the same Ψ value. A viable method for circumventing this shortcoming is to replace the arctan in Eq. (4) by its two-argument form, giving (17)

where φ = arg(x + iy) ∈ [−π, π] is the unique angle satisfying (18)

for any (x, y) with . In particular, we see that for the special cas e r_⋆ = 0, the stream function Φ is simply given by the polar azimuthal coordinate φ, whose range is now [−π, π] instead of the more conventional [0, 2π]. Furthermore, since we then recover the standard definition (19)

of the complex-valued logarithm, the total complex potential (13) simplifies considerably to (20)

in full consistency with Eq. (15).

The use of the above potential theory for an exploration of the principal structure of the flow resulting from the interaction of multiple stellar winds is justified by the notion that this interaction occurs between the subsonic regions of the winds, which in a very good approximation can be considered as incompressible and irrotational (e.g., Siewert et al. 2014; Röken et al. 2015). But even at larger distances from the cluster, where we expect the flow velocity to approach a constant value, and thus mass density ρ to decrease with distance r as ρ(r) ∝ r⁻², the derived flow structure is still valid. This can formally be seen by exchanging u →U: = ρu (see also Kleimann et al. 2017) and noting that U and u share the same structure of stream lines. The possibly questionable constraint of incompressibility (7) now becomes ∇ ⋅ (ρu) = 0 and merely implies mass conservation, which is certainly not under dispute, but in fact constitutes a very desirable property of any physical flow model.

2.3 Stream lines and boundary surfaces

2.3.1 Discontinuity of the stream function

Another technical problem in using the stream function as defined in Eq. (9) is its discontinuity, which becomes evident when considering the limits

with Ψ_[k](x, y) denoting the total stream function dueto all sources except m_k, assuming that none of them is located at x = x_k as well, i.e., Ψ is discontinuous in (x_k, y) where it jumps by ± m_kπ. For the calculation of stream lines via the equationΨ = const., it is therefore necessary to add a suitable term ± m_kπ to the constant, possibly multiple times (see below). This inconvenience can be moderated (though not avoided entirely for principal reasons) by using Eq. (17) in favor of Eq. (4). Since the double-argument arctan covers the full range [−π, π] rather than just [−π∕2, π∕2] as its single-argument form does, Eqs. (21) and (22) are then replaced by (23)

such that there is just one single discontinuity to the left of r_k in the − x direction, across which Ψ jumps by a full ± 2π m_k.

For situations where the stream function is to be evaluated at some position y = y_⋆ that has multiple sources located to its right (also at y_k = y_⋆), Eq. (23) may be generalized to (24)

where the position-dependent index set (25)

serves to enumerate all sources to the right of the position in question (and sharing the same y coordinate).

The need to carefully take these discontinuities into account is illustrated in Fig. 1, in the right panel showing sample stream lines for three equally strong stellar wind sources. The stream lines are computed from Ψ (x, y) using the stream function (17) with jump corrections and labeled by their respective values of Ψ. These corrections are necessary across the negative horizontal lines y = y_k when using Eq. (17) (or, alternatively, the vertical line x = x_k when using Eq. (4)) through the sources as illustrated in the left diagram of the figure.

Fig. 1 Left panel: principal jump corrections of the stream function Ψ induced by the presence of a source of strength m (gray circle) located at rk. The use of the single-argument arctan (4) for the definition of Ψ entrails jumps by ± πm (shown in red) when crossing the vertical line x = xk either “above” (y > yk) or “below” (y < yk) the source, while in the case of definition (17), only a single jump by ± 2πm (shown in blue) occurs across the horizontal line y = yk, x < xk to the “left” of the source. Right panel: sample stream lines for three stellar wind sources with equal strengths m1 = m2 = m3 = 1 (black circle) computed from curves of constant stream function Ψ(x, y) as given in Eq. (17) with jump corrections according to Eq. (24). Colors indicate the various values of Ψ. The two stream lines emanating from the right source at r3 = (1, 0) are both singled valued because there is no other source to its right, while the stream line pointing into the − y direction from the central source at r2 = (0, 1) has a jump of − π at y = y3 = 0. Finally, the source at r1 = (−1, 0) has jumps of ± π in both the upward- and downward-pointing stream line at y = y2 = 1 and y = y3 = 0, respectively.

2.3.2 Boundary surfaces

Assuming laminar flow conditions so that the plasmas of different stellar winds cannot mix, separate regions form in a steady state; these regions contain the outflowing material from the different sources. Consequently, there are boundary surfaces (or boundary lines in a plane) that separate the different stellar wind flows. Those points on the boundaries where the flow velocity vanishes and that are stream line-connected with more than one source are called stagnation points. According to Eqs. (5) and (6) these are the critical points of the stream function, i.e., those with ∂Ψ ∕∂x = ∂Ψ∕∂y = 0. With Eqs. (14) and (15) this condition can be written as (26)

which translates into (27)

after some algebraic manipulations.

The three-star configuration shown in Fig. 1 exhibits two stagnation points, namely, where the associated stream lines seem to cross (see Fig. 2). Of course, they do not actually cross: the flows from the stars decelerate to zero velocity towards a stagnation point, so that a line connecting two sources represents two stream lines with opposite flow directions corresponding to the outflow from each star. The other two lines, not connected to the stars, indicate whereto the material is flowing away from the stagnation points: these lines represent the boundaries separating the nonmixing stellar wind plasmas.

From this symmetric configuration of sources in the corners of an isosceles triangle we can already draw some conclusions. First, it becomes clear that the resulting large-scale wind has individual regions where the outflow remains determined by the star whose material is filling it. Second, it is likely that comparatively narrow outflow channels are forming, one of which is visible in Fig. 2; some of the wind material of the “upper” star escapes from the cluster into the region towards positive y-values, and also within a narrow “channel” forming between the wind regions of the other two stars.

Fig. 2 Stream lines through the two stagnation points of the configuration shown in Fig. 1. The stagnation points are located where the stream lines seem to cross, although they actually do not.

2.3.3 Asymptotic flow behavior

The above findings motivate an analysis of the asymptotic flow behavior. To this end, it is advantageous to transform from Cartesian topolar coordinates (18), such that (28)

We note that for the last step, the use of the double-argument arctan (17) is mandatory; otherwise, arctan[tan(φ)] = φ would only hold for the left (x ≥ 0) half-plane, in which φ ∈ [−π∕2, π∕2]. Thus, for a given jump-corrected value c of a boundary stream line, we obtain from the stream line equation Ψ = c the angle of its asymptotic direction: (29)

Choosing a new origin of the polar coordinate system in the location of the cluster’s barycenter (30)

gives the asymptotic lines for the overall flow configuration shown in Fig. 3.

Fig. 3 Resulting flow of the configuration of stellar wind sources shown in Fig. 1 visualized by many stream lines(blue). In addition, the two boundary lines (green) are indicated, as are their asymptotic directions (red). The corresponding straight lines cross in the origin (Eq. (30)) of the shifted polar coordinate system.

2.4 Principal structure of the large-scale cluster wind

To conclude the treatment of planar configurations, we consider a more realistic, asymmetric configuration of five stars with one wind having only 1∕10 of the source strength of the other four (see Fig. 4).

While the actual flow field is more complicated, its principal structure is analogous to that discussed for the symmetric configuration above. There are five outflow regions, each associated with one of the five sources. The total angular extent covered by any of these regions can either form a continuous interval or be partitioned into at most N − 1 smaller subintervals, some of which may take the form of narrow outflow channels. In the example shown in Fig. 4, the source in the bottom left quadrant has two outflow channels and one large outflow region, while the small central source also has two channels but no large outflow region, and the other three sources each have a continuous outflow region and no channels at all. A simple example featuring the maximum number of outflow channels would be a set of N − 1 sources of uniform strength placed equidistantly on a circle, together with an additional source of arbitrary strength located at the circle’s center. Specifically, if the four equally strong sources depicted in Fig. 4 were located in the corners of a square , the setting’s symmetry would force the four outflow channels to be located around the Cartesian x- and y-axes, carrying material from the central source into the ± x and ± y directions.)

It is interesting to ask what the velocity field far from a cluster looks like. As in the case for the asymptotic directions discussed in the previous section, this information can be obtained from the asymptotics of the combined velocity potential (8). Using again polar coordinates, we obtain (31)

Evidently, the outflow regions are adjusting to each other in such a way that the asymptotic large-scale wind is given by that of a point source with strength , i.e., the combined source strengths of all wind-driving stars in the cluster. Consequently, at large distances the cluster outflow is approximately spherically symmetric.

Fig. 4 Flow configuration resulting from five stellar wind sources, with the one star in the origin having 1/10 of the strength of the other four.

3 Extension to three dimensions

Despite the valuable insights obtained from the analytical modeling discussed so far, the major drawback of this model lies in its limitation to a planar 2D configuration. While it could be argued that the topological structure of a 2D flow model will be very similar, if not identical, to that of a 3D model where both the distribution of sources and the region of analysis is limited to a plane (provided that only locations within this plane are considered), there will of course be quantitative differences. Most notably, the flow speed |u| of a 2D point source (which actually corresponds to a line source in 3D) decreases with radius as r⁻¹, rather than r⁻². For these reasons, we now turn to the case of asymmetric arrangements of point sources in full 3D.

While itis of course still possible to compute the overall flow field from the superposition of the individual velocity potentials of the sources distributed in 3D, the easy access to stream lines via a stream function is no longer available in that case, and a numerical procedure becomes necessary. Rather than improving the model further by actually computing the interaction of multiple stellar winds on the basis of numerical solutions of the (M)HD equations (which we intend to study in future work), we now describe this procedure, the underlying mathematical theory, and the results thus obtained.

3.1 Mathematical background

The statements from the previous section generalize to 3D as follows. For a set of N point sources with strengths m_k and locations r_k = (x_k, y_k, z_k) in 3D space, the total velocity field found from the superposition (32)

can have at most N − 1 stagnation points (also called nulls). In the vicinity of a null at r_null, u(r) can be approximated by (33)

where is the Jacobian matrix of u with elements (34)

evaluated at the null. Since u derives froma potential, we see from (35)

that J is symmetric, implying that all of its eigenvalues are real, and that the corresponding eigenvectors are orthogonal (or at least can be chosen to have that property). Furthermore, the trace (36)

of J, i.e., the sum of its eigenvalues, vanishes for incompressible flow in any dimension d. These constraints limit the topological variety that would otherwise be possible for nulls in the following cases (see, e.g., Parnell et al. 1996):

• In 2D, only X-type nulls are possible, i.e., stream lines in the direction (anti)parallel to the negative eigenvalue’s eigenvector converge towards the null (thus connecting it to either a source or another null), and an additional pair of stream lines emanate from the null into the perpendicular direction (as shown in Fig. 2). The additional stream lines form separatrices, separating regions filled by material from different sources. The second possibility of an O-type null (with circular stream lines closing in on the null in their common center) is ruled out by ∇×u =0.

• In 3D, the situation is similar. We also get individual “spine” stream lines that originate at a source and converge on the null (in the direction of the negative eigenvalue’s eigenvector), but now the remaining two eigenvectors (whose eigenvalues are both positive) span an entire “fan” plane that is orientated perpendicular to the spine. The set of stream lines starting from the null within this plane form separatrix surfaces, which again separate the regions that are filled by material from different sources.

It is also reasonable to generalize the reasoning presented in Sect. 2.3.3 for the 3D case. At large distances from the cluster of point sources, the total velocity field approaches that of a monopole (i.e., a single point source whose strength is the sum of the individual strengths of the contributing sources), and when placing a sufficiently large sphere around the cluster’s barycenter, the wind from any source j fills a solid angle μ_j on that sphere given by (37)

irrespective of how complex the spatial arrangement of both the sources in the cluster and their corresponding regions on the sphere may be. This finding, which seems surprising at first sight, is illustrated in the right panel of Fig. 6 and quantitatively supported by the data listed in Table 1.

Fig. 5 Left panel: curves on which the x component (red) and the y component (blue) of the combined velocity field of four sources (filled gray circles) vanish. Null points (filled yellow circles) are identified as those locations where both lines meet outside of any source. Straight dashed lines connect those initial centers c (crosses) of a rectangular cell that did eventually converge onto a null to this null (which apparently need not be the null closest to c’s initial position). We also note that these dashed lines merely serve to illustrate which starting point converges upon which null, while convergence itself usually does not proceed along straight lines. Right panel: same situation, indicating identified nulls at (0.5316, 1.4514), (1.4755, 0.8766), and (2.0562, 1.8008), separatrices (gray), spine stream lines (green), and unit vectors of the actual flow field u (brown arrows).

Fig. 6 Left panel: 3D rendering showing the five sources (gray spheres, volumes proportional to strength of the source), four null points (yellow spheres), spine stream lines (black lines), and four separatrix surfaces (semi-transparent blue, red, green, cyan)emanating from the nulls, together with some of their constituting fan stream lines. Only parts of objects within a sphere of radius 5 around the origin are shown. Right panel: perspective projection of a spherical surface of radius R = 10, color-coded by the source from which stream lines passing through a given point on the sphere originate, thus indicating the extent and geometrical shape of the solid angle that is filled by the wind from the corresponding stellar source. The colors of sources 1 (cyan), 2 (blue), 3 (red), and 4 (green) correspond to those of the null (as used in the left image) to which the respective source is connected by a streamline. The patch for the central source (filling the small lanes between the other four patches) is omitted.

3.2 Finding null points and separatrices in three dimensions

Since the structure of a given flow configuration tends to be more difficult to visualize in 3D, it is important for the analysis to identify the “skeleton” consisting ofcritical field lines, which are those of the spine and the fan. The first step towards this goal is to locate all the stagnation points. Since the structure of u resulting from Eq. (32) is not conducive to an analytical approach (not even in 2D, where solving Eq. (27) amounts to finding all complex roots of an Nth-order polynomial), we resort to numerical methods. The following procedure, which is a generalization of the Newton–Raphson method, hasshown to yield good results for the purpose at hand (for the sake of clarity, Fig. 5illustrates the situation in a 2D plane, but the procedure in 3D is completely analogous):

• 1.

  A rectangular region containing all sources is partitioned into an (n × n × n) array of brick-shaped cells of equal size, where typically n ≲ 10.

• 2.

  At the center c of each cell, each velocity component u_α, α ∈{x, y, z} is approximated by its first-order Taylor polynomial (38)

  and the position c₀ is determined at which p_α(c₀) = 0 for all α. If c₀ is ill-defined, which can happen if the subspaces defined by the zeroes of any two polynomials are parallel, the cell is discarded, else it is centered on c₁ and reduced in size by a factor of 1/2 in all directions.

• 3.

  Step 2 is repeated with c₀ as the new c until the cell size has shrunk below a pre-defined threshold (typically ≲ 10⁻⁸).

• 4.

  The final c₀ is written into a list of prospective nulls (if it does not coincide with a source), and the process is repeated for the next cell.

• 5.

  Finally, the list needs to be cleared of duplicates that occur whenever two cells converge upon the same null. The list must be checked that it contains exactly N − 1 distinct points at which u is indeed sufficiently close to zero.

This procedure was preferred over the topological degree method (Greene 1992) because it is easier to implement. It is similar to the trilinear method proposed by Haynes & Parnell (2007) except that the field components within a cell are assumed to be planar, rather than trilinear, and that a given null may find multiple times from different starting cells (and in fact usually does). We also skip the initial preconditioning of cells since for small cell arrays the gain in execution speed does not justify the extra coding effort.

3.3 Numerical example

From the list of nulls thus identified, we can continue along the lines of Sect. 3.1 by numerically evaluating the Jacobian and its eigenvectors at each null, and then trace stream lines along the spine and into selected directions within the fan plane. As an instructive example, we consider the configuration specified in Table 1, which consists of four unequally strong sources forming an irregular tetrahedron, together with a fifth source placed in their common barycenter at r₅ = 0. As expected, four nulls are found, situated approximately between the central source and each of the four vertices of the tetrahedron.

Having located the four stagnation points, we can generate a pair of stream lines in the spine direction of each null, as well as a range of stream lines (80 for each null in this case) starting in the fan plane by numerical second-order Runge–Kutta integration. The resulting structural skeleton is shown in the left plot of Fig. 6. As can be seen, the separatrices form umbrella-like shapes around the four outer sources, approaching the shape of wide cones at larger distances.

We use this example to confirm the prediction of Eq. (37). For this purpose, we construct the contour of points at which each of our fan stream lines pierces a sphere of radius R (see the right plot of Fig. 6, which provides a visual impression of these polygonal patches), and determine the area which these polygons encompass on the sphere. As can be seen from the data presented in the right part of Table 1, the solid angles thus computed are indeed consistent with the predicted asymptotes to an excellent degree.

4 Summary and conclusions

With the present study we addressed the problem of the mutual interaction of multiple winds from stars located within a stellar cluster. While several studies have investigated large-scale cluster winds resulting from an intra-cluster interaction of multiple stellar winds, they have as yet not provided details of the colliding flows inside a given cluster. The present work represents a further step towards a quantification of the principal structure of the combined flow within and near such cluster. To this end we have applied the theory of complex velocity potentials to determine the stagnation points and the separatrices between the individual stellar wind regions originating from a planar arrangement of point sources, as well as the structure of the asymptotic large-scale wind outside the cluster. While in two spatial dimensions the use of complex-valued potentials allows a concise analytical treatment of stream lines using stream functions, the extension to the more realistic 3D case requires recourse to numerical methods (although the general theory of vector-valued nulls may be explored analytically, as has in fact been done by various authors).

The main results can be summarized as follows. First, one can distinguish distinct regions in the large-scale cluster wind that are filled with plasma from the individual stars. Second, as a consequence of the strengths and locations of the wind sources, the outflow of a given star when projected onto a sphere does not necessarily result in a singly connected region. The wind of a given star may be forced, partly or fully, into comparatively narrow outflow channels. Third, as expected, one can demonstrateanalytically that the large-scale cluster wind asymptotically approaches spherical symmetry at large distances.

In conclusion, we state that the combined flow inside a stellar cluster resulting from the interaction of multiple stellar winds is highly structured. This confirms the expectation, and is the motivation to carry out further studies on the basis of (magneto) hydrodynamical simulations.

Acknowledgements

We are grateful to the Deutsche Forschungsgemeinschaft (DFG) for the funding project SCHE334/9-2. Furthermore, J.K. acknowledges financial support through the Ruhr Astroparticle and Plasma Physics (RAPP) Center, funded as MERCUR project St-2014-040.

References

All Tables

All Figures

Fig. 1 Left panel: principal jump corrections of the stream function Ψ induced by the presence of a source of strength m (gray circle) located at rk. The use of the single-argument arctan (4) for the definition of Ψ entrails jumps by ± πm (shown in red) when crossing the vertical line x = xk either “above” (y > yk) or “below” (y < yk) the source, while in the case of definition (17), only a single jump by ± 2πm (shown in blue) occurs across the horizontal line y = yk, x < xk to the “left” of the source. Right panel: sample stream lines for three stellar wind sources with equal strengths m1 = m2 = m3 = 1 (black circle) computed from curves of constant stream function Ψ(x, y) as given in Eq. (17) with jump corrections according to Eq. (24). Colors indicate the various values of Ψ. The two stream lines emanating from the right source at r3 = (1, 0) are both singled valued because there is no other source to its right, while the stream line pointing into the − y direction from the central source at r2 = (0, 1) has a jump of − π at y = y3 = 0. Finally, the source at r1 = (−1, 0) has jumps of ± π in both the upward- and downward-pointing stream line at y = y2 = 1 and y = y3 = 0, respectively.

In the text

	Fig. 2 Stream lines through the two stagnation points of the configuration shown in Fig. 1. The stagnation points are located where the stream lines seem to cross, although they actually do not.
In the text

	Fig. 3 Resulting flow of the configuration of stellar wind sources shown in Fig. 1 visualized by many stream lines(blue). In addition, the two boundary lines (green) are indicated, as are their asymptotic directions (red). The corresponding straight lines cross in the origin (Eq. (30)) of the shifted polar coordinate system.
In the text

	Fig. 4 Flow configuration resulting from five stellar wind sources, with the one star in the origin having 1/10 of the strength of the other four.
In the text

Fig. 5 Left panel: curves on which the x component (red) and the y component (blue) of the combined velocity field of four sources (filled gray circles) vanish. Null points (filled yellow circles) are identified as those locations where both lines meet outside of any source. Straight dashed lines connect those initial centers c (crosses) of a rectangular cell that did eventually converge onto a null to this null (which apparently need not be the null closest to c’s initial position). We also note that these dashed lines merely serve to illustrate which starting point converges upon which null, while convergence itself usually does not proceed along straight lines. Right panel: same situation, indicating identified nulls at (0.5316, 1.4514), (1.4755, 0.8766), and (2.0562, 1.8008), separatrices (gray), spine stream lines (green), and unit vectors of the actual flow field u (brown arrows).

In the text

Fig. 6 Left panel: 3D rendering showing the five sources (gray spheres, volumes proportional to strength of the source), four null points (yellow spheres), spine stream lines (black lines), and four separatrix surfaces (semi-transparent blue, red, green, cyan)emanating from the nulls, together with some of their constituting fan stream lines. Only parts of objects within a sphere of radius 5 around the origin are shown. Right panel: perspective projection of a spherical surface of radius R = 10, color-coded by the source from which stream lines passing through a given point on the sphere originate, thus indicating the extent and geometrical shape of the solid angle that is filled by the wind from the corresponding stellar source. The colors of sources 1 (cyan), 2 (blue), 3 (red), and 4 (green) correspond to those of the null (as used in the left image) to which the respective source is connected by a streamline. The patch for the central source (filling the small lanes between the other four patches) is omitted.

In the text