3 Results

We first checked our work against the solutions obtained by Dyson (1973) for WBBs with no mass loading. For this comparison it was required that $\lambda = -3$ and that $\phi$ was large. For $\lambda = -3$ the radius of minimum mass loading increases as t. Excellent agreement with Dyson (1973) was found over a large range of $x_{\rm is}/x_{\rm cd}$. For large $\phi$ (i.e. negligible mass loading) the value of $x_{\rm ml}$ has no effect on the resulting solution. Figure 1 shows a sample solution.

Figure 2: Results for $\lambda = -3$, $\phi = 100$, and various values of $x_{\rm ml}$. For these parameters the mass loading is small but non-negligible, and hence demonstrates $x_{\rm ml}$ having a minor influence on the results. Results for five values of $x_{\rm ml}$ are shown: 0.1 (solid), 0.2 (dots), 0.33 (dashes), 0.5 (dot-dash), and 0.9 (dot-dot-dot-dash). The panels are the same as in Fig. 1. In Fig. 3 we show solutions where the value of $x_{\rm ml}$ has a much greater influence.

In Fig. 2 we show results for $\lambda = -3$, $\phi = 100$, with different values of $x_{\rm ml}$. The mass loading in these results is small, but not negligible, so that the precise value of $x_{\rm ml}$ has minor consequences for the solution obtained. In Table 1 we tabulate important parameters from these solutions. In particular, we find that values for the ratio $x_{\rm is}/x_{\rm cd}$ and the energy fractions are very similar. The fractional mass loading of the bubble, $\Phi_{\rm b}$, is somewhat more sensitive to the value of $x_{\rm ml}$, and increases when the "turn-on'' radius for mass loading is decreased, as one would expect. The preshock Mach number also reflects the degree of preshock mass loading, again as expected. Interestingly, the profiles of postshock Mach number are, however, very insensitive to the preshock Mach number. With the assumption of constant $\phi$, we note that since $x_{\rm is}/x_{\rm cd}$ varies for the above results it was important to see whether this was a possible factor for some of the differences in these solutions. Therefore we also investigated the effect of different values of $x_{\rm ml}$ whilst keeping $x_{\rm is}/x_{\rm cd}$ constant and varying $\phi$ as necessary. We find that the resulting solutions are again fairly similar, and differ at about the same level as the results in the top half of Table 1. Therefore it seems that varying $x_{\rm ml}$ has a direct effect on the results, and not just an indirect influence through changes induced in $x_{\rm is}/x_{\rm cd}$ and/or $\phi$. In the following, therefore, we will vary $x_{\rm ml}$ whilst keeping $\phi$fixed.

Figure 3: Results for $\lambda = -3$, $\phi = 10$, and various values of $x_{\rm ml}$. For this value of $\phi$, the mass loading is appreciable, and hence the solutions shown demonstrate $x_{\rm ml}$ having a major influence. Results for four values of $x_{\rm ml}$ are shown: 0.02 (solid), 0.1 (dots), 0.5 (dot-dash), and 0.9 (dot-dot-dot-dash). The panels again represent density, velocity, internal energy and Mach number.

At the other extreme, Fig. 3 shows results where the value of $x_{\rm ml}$ can have large consequences for the solutions obtained. This occurs for values of $\phi$ for which the resulting mass loading of the bubble is important. Parameters from these solutions are tabulated in the lower half of Table 1. Although the ratio $x_{\rm is}/x_{\rm cd}$ remains relatively insensitive to the value of $x_{\rm ml}$, the fractional mass loading $\Phi_{\rm b}$, the ratio of the shell mass to the bubble mass $M_{\rm sh}/M_{\rm b}$, and the energy partition are all significantly affected, as Table 1 shows. This is expected given that most of the mass loading occurs near the minimum mass loading radius for $\lambda = -3$.

Figure 4: Results for $\lambda = -3$, $\phi = 50$. Here mass loading is so severe that the flow directly connects to the contact discontinuity. For the solution shown, the onset of mass loading occurs at $x = 0.359 x_{\rm cd}$. Once the flow starts mass loading, its velocity drops, its relative density increases, its kinetic energy drops, and its thermal energy increases.

An interesting consequence of mass loading the wind prior to the inner shock is that if the mass loading is large, the Mach number of the flow can be reduced so much that the flow directly connects to the contact discontinuity without the presence of an inner shock. In Fig. 4 we show one solution where this happens. In this example, the flow is continuously supersonic with respect to the clumps, and we obtain $\Phi_{\rm b} = 1.25$ and $M_{\rm sh}/M_{\rm b} = 0.07$. This perhaps indicates that this morphology is most favoured at early stages of the evolution of the bubble before a significant amount of ambient medium is swept-up. In another example of "direct connection'', we find a solution where the Mach number of the flow with respect to the clumps drops below unity (see panel f of Fig. 6). In this case, $\Phi_{\rm b} = 3.4$ and $M_{\rm sh}/M_{\rm b} = 480$, suggesting a much older bubble. We also find that low values of $\lambda$ ( $\approx -3$) are apparently needed to obtain directly-connected solutions.

Figure 5: Comparison of solutions with negligible (solid) and high (dots) mass loading rates, for a given value of $\lambda$ (-2) and $x_{\rm is}/x_{\rm cd}$ (0.606). The ratio of swept-up mass to bubble mass is 12.3 and 0.96 respectively, whilst the ratio of ablated mass to wind mass is 7.54 for the mass-loaded solution. The preshock Mach number for the mass-loaded solution is 1.006. In both cases, the similarity solution has been scaled such that the current wind parameters are appropriate for a WR star: $\hbox{${\dot M}$ } \,_{\rm c} = 10^{-5}~M_{\odot}~{\rm yr}^{-1}$ and $v_{\rm w} = 1500$ km s-1.

In Fig. 5 we compare solution profiles for the case where $\lambda = -2$ and for negligible and high mass loading rates. A $\lambda = -2$ clump distribution is the most physically plausible for many astrophysical objects. For this distribution the minimum mass loading radius increases as t^2/3, which is not far removed from a linear increase with t. The solutions have been scaled from the similarity results by assuming current values for the wind parameters appropriate for a WR star: $\hbox{${\dot M}$ } \,_{\rm c} = 10^{-5} M_{\odot}~{\rm yr}^{-1}$ and $v_{\rm w} = 1500$ km s^-1. The mass-loss rate increases with time as t^2/3 whilst the wind velocity decreases as t^-1/3, so the wind has evolved from a case more suitable to an O-star. The assumed bubble age is t = 10⁴ yr. The interclump medium varies as $\rho \propto r^{-1/2}$ (i.e. $\beta = -1/2$) which is a flatter distribution than for a constant velocity wind (r^-2). However, since the wind mass-loss rate is increasing and the velocity decreasing, an r^-1/2 interclump density profile is physically realistic for such a scenario. In comparison with the negligibly mass-loaded solution, the increase in density inside the bubble for the highly mass-loaded solution can be clearly seen, along with a decrease in preshock velocity, preshock Mach number, and postshock temperature. Interestingly, the Mach number of the postshock flow is higher for the heavily mass-loaded solution. Note also that the heavily mass-loaded solution extends to a larger bubble radius. This is a consequence of the different preshock interclump densities: for the heavily mass-loaded solution $n_{\rm e} = 0.066$, whilst $n_{\rm e} = 0.195$for the negligibly mass-loaded solution.

Figure 6: The Mach number as a function of x for a number of different solutions. With respect to the clumps, the shocked region can be either entirely subsonic, entirely supersonic, or have one (or maybe more) sonic points. The panels show: a) shocked region with a monotonic decrease in Mach number ( $\lambda = -2$, $\phi = 0.5$, $x_{\rm is}/x_{\rm cd} = 0.585$, $x_{\rm ml} = 0.1$); b) shocked region entirely supersonic with a monotonic increase in Mach number ( $\lambda = -3$, $\phi = 100$, $x_{\rm is}/x_{\rm cd} = 0.902$, $x_{\rm ml} = 0.9$); c) postshock Mach number initially subsonic and decreasing, before levelling off and subsequently increasing through a sonic point ( $\lambda = 4$, $\phi = 10$, $x_{\rm is}/x_{\rm cd} = 0.297$, $x_{\rm ml} = 0.1$); d) shocked region entirely subsonic with a decreasing Mach number followed by a modest increase close to the contact discontinuity ( $\lambda = -3$, $\phi = 10^{5}$, $x_{\rm is}/x_{\rm cd} = 0.181$, $x_{\rm ml} = 0.1$); e) postshock flow continuously supersonic though with an initial decrease in Mach number ( $\lambda = -3$, $\phi = 10$, $x_{\rm is}/x_{\rm cd} = 0.816$, $x_{\rm ml} = 0.02$); f) continuously supersonic flow (solid) directly connected to the contact discontinuity ( i.e. no inner shock) ( $\lambda = -3$, $\phi = 50$, $x_{\rm ml} = 0.359$) and a flow becoming subsonic (dots) directly connected to the contact discontinuity ( $\lambda = -3.2$, $\phi = 12.5$, $x_{\rm ml} = 0.246$).

In PDH it was shown that there is a negative feed-back mechanism caused by the evaporation of mass from embedded clumps, which set a maximum limit to the amount that a bubble could be massloaded. Additionally, for the bubble mass to be greater than the mass of the swept-up shell, a large radial dependence of the mass loading was required ( $\lambda \geq 4$). In this work we instead find that small values of $\lambda$ are required to satisfy this condition. Our simulations show that it is impossible to obtain a similarity solution with $M_{\rm sh} < M_{\rm b}$ for $\lambda \geq -1$ (e.g. for $\lambda = 1$, ). However, for smaller values of $\lambda$, solutions satisfying $M_{\rm sh} < M_{\rm b}$ can be found (e.g. for $\lambda = -2$ (-3), we can obtain $M_{\rm b} \approx 8$ (29) $M_{\rm sh}$, where $\Phi_{\rm b} \approx 7.6$ (15)). For a given value of $M_{\rm sh}/M_{\rm b}$ there appears to be a maximum value for $\Phi_{\rm b}$. It occurs when mass loading starts very close to the bubble centre (i.e. $x_{\rm ml} \rightarrow 0$) and when the flow remains supersonic relative to the clumps over the entire bubble radius.

Figure 7: Radial profiles of emissivity per unit volume and temperature, normalized to values of 1.0 at the limb. The upper panels show results where an inner shock is present, whilst the bottom panels show results where the mass loading region directly connects to the contact discontinuity. For the upper panels the results shown are for: $\lambda = -2$, $\phi = 0.5$, $x_{\rm is}/x_{\rm cd} = 0.59$, $x_{\rm ml} = 0.1$ (solid); $\lambda = -3$, $\phi = 100$, $x_{\rm is}/x_{\rm cd} = 0.90$, $x_{\rm ml} = 0.9$ (dots); $\lambda = -3$, $\phi = 10$, $x_{\rm is}/x_{\rm cd} = 0.82$, $x_{\rm ml} = 0.02$ (dashes). For the lower panels the results are for: $\lambda = -3$, $\phi = 50$, $x_{\rm ml} = 0.36$ (solid); $\lambda = -3.2$, $\phi = 12.5$, $x_{\rm ml} = 0.25$ (dots).

We have also investigated the Mach number profile of the solutions that we obtain, of which some are already shown in Figs. 2 and 3. As for the evaporative mass loading in PDH, we note that a number of different forms are possible, and summarize these in Fig. 6.

Finally, radial profiles of the X-ray emissivity per unit volume and temperature have been calculated (Fig. 7). We assume that the emissivity $\Lambda \propto n^{2} T^{-1/2}$, which is a good approximation over the temperature range (cf. Kahn 1976). The upper panels show results where an inner shock is present, and we will discuss these first. The solution with the solid line displays relatively constant values of $\Lambda$ and T, and is characterized by a high degree of mass loading ( $\Phi_{\rm b} = 9.7$). Note that the highest temperature in the bubble does not occur postshock, but rather shortly after mass loading begins. The very high emissivity at the onset of mass loading is a consequence of the low temperature at this point and the assumed T^-1/2 dependence of our emissivity, and is not physical. The solution with the dotted line has an interior dominated by free-expansion of the wind source. Mass loading only occurs over the final 20 per cent of the radius, and is relatively weak, and the maximum temperature in the bubble occurs immediately postshock. The solution with the dashed line has some characteristics of each of the other solutions. It is mass-loaded from small radius which again leads to the maximum bubble temperature occurring near the center. However, the mass loading is less severe ( $\Phi_{\rm b} = 6.5$) and its radial dependence steeper ( $\lambda = -3$) than for the solution with the solid line, so that there is a large fall in density between the switch-on radius for mass loading and $x_{\rm is}$. This leads to a steep rise in emissivity from the inner shock to the bubble center in a similar fashion to the solution with the dotted line.

The bottom panels of Fig. 7 show results where there is no inner shock and where the mass loading region directly connects to the contact discontinuity. The Mach number profiles for these solutions are displayed in the bottom-right panel of Fig. 6. As can be seen, the temperature profiles are approximately flat over their respective mass loading regions, and are very similar in shape to each other. This contrasts with the behaviour of their emissivity profiles, where it is clear that the solution represented by the solid (dotted) line has $\Lambda$ generally decreasing (increasing) with radius.

A general result from Fig. 7, which was also found in the evaporatively mass-loaded bubbles of PDH, is that higher mass loading tends to decrease the central emissivity and increase the central temperature relative to the limb.

   
3.1 Comparison with observations and other theoretical work

Although our work cannot be directly compared to the hydrodynamic simulations of Arthur et al. (1993, 1996), a number of similarities exist. These include a sharp increase in temperature at the onset of mass loading, and often $\partial \rho/\partial x > 0$ over the mass-loaded region. A noticeable achievement of the work of Arthur et al. (1996) was the reproduction of the observed correlation between the velocity and ionization potential of ultraviolet absorption features towards the central star (Smith et al. 1984). Our solutions also often reproduce this correlation (e.g. see Fig. 5 which shows the velocity and temperature decreasing together for the heavily mass-loaded solution), as do those of Hanami & Sakashita (1987). Interestingly, solutions with mass-loading from conductive evaporation of the swept-up shell do not reproduce this correlation (see Weaver et al. 1977).