HELIOSPHERIC STRUCTURE: THE BOW WAVE AND THE HYDROGEN WALL

The interaction of the supersonically expanding solar wind with the surrounding interstellar medium (ISM) is mediated significantly by interstellar neutral gas. The charge exchange of neutral hydrogen (H) with supersonically flowing solar wind protons yields a supra-thermal pickup ion (PUI) population that removes momentum and energy from the solar wind, leading to its gradual deceleration (Holzer 1972; Isenberg 1986; Richardson et al. 1995; Zank 1999). In the outer heliosphere, the PUIs form a thermally dominant proton component co-moving with the solar wind (e.g., Williams & Zank 1994, see Zank 1999 and Zank et al. 2009 for extensive reviews). PUIs further mediate the heliospheric termination shock (HTS) and experience preferential heating compared to solar wind protons (Zank et al. 1996a, 2010; Richardson et al. 2008; Richardson 2008; Burlaga et al. 2008; Burrows et al. 2010). This yields a complicated proton distribution downstream of the HTS (i.e., in the inner heliosheath; Zank et al. 2010) that in turn creates complex energetic neutral atom (ENA) spectra that can be observed by the NASA Interstellar Boundary Explorer (IBEX) spacecraft (McComas et al. 2009; Funsten et al. 2009; Schwadron et al. 2009; Desai et al. 2012; Livadiotis et al. 2012). The ability to utilize observations of ENAs created in the boundary regions of the interacting solar wind and local interstellar medium (LISM) has yielded considerable insight into the global structure of the heliosphere. Particularly notable was the discovery of the IBEX Ribbon (McComas et al. 2009) and the subsequent development of models (McComas et al. 2009; Schwadron et al. 2009; Heerikhuisen et al. 2010) that, depending on the ribbon generation mechanism, may allow for estimates of the LISM magnetic field strength and direction to be inferred (Heerikhuisen & Pogorelov 2011).

Most recently, observations derived from the IBEX-LO instrument have provided a detailed parameterization of the local interstellar medium (LISM) properties (Lee et al. 2012; Möbius et al. 2012; Saul et al. 2012; Bochsler et al. 2012; Bzowski et al. 2012; Hłond et al. 2012). In particular, the IBEX measurements revised earlier values of the interstellar flow vector derived from Ulysses observations (Witte 2004), suggesting a value of 23.2 km s⁻¹ rather than the previous value of 26 km s⁻¹. McComas et al. (2012) resolved a small discrepancy between the inferred LISM parameters (i.e., the parameterized plasma velocity and temperature) and combined them with estimates of the density and magnetic field strength and orientation to conclude that there is a good possibility that the interaction of the solar wind with the LISM is sub-fast magnetosonic and that a bow shock is likely not present. This is a provocative and challenging set of observations for theoretical models of the heliospheric interaction with the LISM to address.

Although some of the earlier original models had considered the possibility of a one-shock model (i.e., an HTS that bounds the supersonic solar wind and no bow shock), including Parker's original gas dynamic model (Parker 1961) and the Zank et al. (1996b) one-shock models that incorporated neutral H self-consistently, the general earlier consensus (circa 2000) was that the interaction of the solar wind with the LISM was of a two-shock character, i.e., the HTS and a bow shock existed due to the supersonic relative motion of the Sun with respect to the LISM (Baranov et al. 1971; Baranov & Malama 1993; Pauls et al. 1995; Zank et al. 1996b; Williams et al. 1997; Pauls & Zank 1996, 1997), even for an exceptionally strong LISM magnetic field (Florinski et al. 2004).

The revised plasma parameters provide a challenge in understanding the nature of the LISM transition upwind of the heliopause due to the very non-local character of the charge exchange process in the immediate interstellar neighborhood of the heliopause. The creation of either fast or hot neutral H in the supersonic solar wind or inner heliosheath respectively leads to their propagation into the very local ISM, where they experience secondary charge exchange (Zank et al. 1996b). This process may well be responsible for the formation of the IBEX Ribbon (McComas et al. 2009; Heerikhuisen et al. 2010; Chalov et al. 2010). More importantly, secondary charge exchange in the LISM is dynamically important because it acts to heat the interstellar plasma (Zank et al. 1996b) to the extent that even a nominally super-fast magnetosonic (as measured at a large distance) relative heliosphere–LISM flow can become sub-fast-magnetosonic close to the heliopause and the bow shock may appear to be smoothed or even absent (Pogorelov et al. 2006, 2008, 2009b). Furthermore, even if a weak bow shock exists, it is possible that charge exchange may be the primary dissipation mechanism and the overall shock transition becomes smoothed over the charge exchange length scale, giving an adiabatic-like interaction region that is essentially "bow-shock-free."

Given that we can expect the nature of the bow-shock transition to be different in light of the revised LISM parameters, it is entirely possible that the nature of the hydrogen wall (H-wall) may be different from our conventional expectations. Many models, beginning with that of Gayley et al. (1997) based on the observations of Linsky & Wood (1996), have tried to use Lyα observations toward nearby stars to constrain models of heliospheric structure. These models rely on the additional absorption of photons by the H-wall. Models that account for the observed Lyα observations have almost inevitably been of the two-shock kind, so revisions to the basic structure of the heliosphere to account for the IBEX observations require us to reconsider the constraints imposed by the Lyα observations along multiple lines of sight.

In this work, we consider three models using the current best values for the LISM flow speed and plasma temperature and vary the LISM magnetic field strength and plasma and neutral H density. Our choice of parameters (Table 1) yields an interstellar plasma flow that is either (1) easily super-Alfvénic and super-fast magnetosonic, (2) barely super-Alfvénic and super-fast magnetosonic, or (3) easily sub-Alfvénic and sub-fast magnetosonic. The precise values are listed in Table 1. In the absence of shock mediation, each of these cases should correspond respectively to (1) the existence of a relatively weak but well-defined fast mode bow shock, (2) the existence of a very weak fast mode bow shock, and (3) the absence of a bow shock ahead of the heliospheric obstacle. These inferences are based of course on a purely magnetohydrodynamic (MHD) interpretation of the interaction, and, as discussed above, the non-local character of the charge exchange process introduces non-classical elements into the MHD description (both in terms of plasma heating and mediation of shock waves). The important point about the values listed in Table 1 is that all the models have an H number density n_H ∼ 0.1 cm⁻³ at the HTS and a heliocentric distance to the HTS of about 89 AU in the Voyager 1 and Voyager 2 directions—see J. Heerikhuisen et al. (2013, in preparation) for a more detailed discussion. In the following section, we discuss the models and underlying physics before using the modeled H data to determine synthetic Lyα absorption profiles along various lines of sight. The model absorption profiles are compared to observed profiles along these directions. We conclude by discussing the model implications for the structure of the heliosphere and the nature of a bow shock or wave.

The LISM parameters assumed at a heliocentric distance of 1000 AU for each of the three models are listed in Table 1. We consider a steady-state solar wind model and assign the following standard parameters at 1 AU: n_p(1 AU) = 7.4 cm⁻³; T_p(1  AU) = 51, 100 K; U_SW(1  AU) = 450 km s⁻¹, and |B|(1  AU) = 37.5 μG. In all three cases, the HTS is located at approximately the same distance, ∼89 AU, along the Voyager 1 and 2 trajectories. For these simulations, we use the Huntsville three-dimensional (3D) MHD plasma–kinetic neutral H code MSFLUKSS (Pogorelov et al. 2004, 2006, 2008, 2011; Heerikhuisen et al. 2006, 2007). This particular model uses a kappa distribution (with κ = 1.63 everywhere) for the inner heliosheath plasma (Heerikhuisen et al. 2008; see Livadiotis & McComas 2009), which allows us, at a relatively simplified level, to approximate the complex proton distribution in the inner heliosheath (Zank et al. 2010; Desai et al. 2012). The neutral H, both interstellar and that generated by the interaction of LISM neutral H with the solar wind, is determined by solving the corresponding charge exchange form of the Boltzmann equation with appropriate source terms from the various plasma distributions (Pauls et al. 1995; Zank et al. 1996b; Zank 1999; Müller et al. 2000; Heerikhuisen et al. 2006). We compute the appropriate source terms for the plasma equations from the kinetic neutral H distribution, thereby ensuring a self-consistent coupling of plasma and neutral H throughout the solar-wind–LISM interaction region. The source terms are constructed based on the statistics of individual events (see Heerikhuisen & Pogorelov 2010 for a description of the method).

Although we utilize a steady-state model for the solar wind, we should be cognizant of the importance of the variable solar wind. Factors such as solar cycle and even longer term variability in the ram pressure (McComas et al. 2000, 2008) will affect the location and asymmetry of the HTS and heliopause on corresponding timescales (Zank & Müller 2003; Scherer & Fahr 2003; Izmodenov et al. 2005b; Pogorelov et al. 2009a, 2011), and solar wind disturbances, especially during solar maximum, lead to significant movement of the HTS and heliopause and the generation of considerable variability in the inner heliosheath and even beyond (Zank & Müller 2003; Washimi et al. 2011, 2012). Indeed, the solar wind induced disturbances propagating into the LISM will heat the very localized plasma weakly and may even create surrogate bow waves or shocks under LISM plasma conditions corresponding to the one-shock model. Solar wind variability is likely to be more important for the case of a marginally super- or sub-Alfvénic LISM flow.

A much more detailed discussion of the three models used here is presented in J. Heerikhuisen et al. (2013, in preparation). As illustrated in Figure 1, Model 1 (|B| = 2 μG, n_p = 0.13 cm⁻³) very clearly possesses a bow-shock-like transition, although comparatively weak. The top row of Figure 1 plots the plasma temperature on a logarithmic scale and, although a little difficult to discern in these 2D plots, a temperature gradient is present, extending from the heliopause into LISM beyond the bow shock. The gradient, as we show below, is more apparent in the 1D cuts below (see Figure 4 where we show values along the α-Cen sightline). The presence of the temperature gradient was discussed by Zank et al. (1996b) who compared a gas-dynamic-only (i.e., no neutral H present) 2D model of the solar-wind–LISM interaction to an otherwise identical self-consistently coupled gas dynamic–neutral H model. In the case of the gas-dynamic-only model, Zank et al. found that no temperature gradient was present. However, the self-consistent model showed that hot neutral H created in the hot heliosheath propagated into the LISM where it eventually experienced charge exchange to produce hot protons in the LISM. This process leads to the formation of an extended temperature gradient. Such non-classical transport of heat across the heliopause (a tangential discontinuity), mediated by neutral H and the charge exchange process, is quite unlike standard MHD. As a result of the temperature gradient, there is a gradual change in the sound (and hence magnetosonic) speed, which changes as T^1/2, increasing as it approaches the bow shock and heliopause region. This effect can be sufficiently strong, as noted by Pogorelov et al. (2008), that the sound speed increase can eliminate the formation of a bow shock, even for super-fast magnetosonic far upstream conditions. The parameters of Model 1 indicate, however, that even though the bow shock is somewhat weakened by the increasing sound speed, a bow shock is nonetheless necessary. The shock transition illustrated in Figure 1 is uncharacteristically not abrupt. Also shown are cuts in the ecliptic and polar planes of the neutral H number density, illustrating the presence and extent of the H-wall, which begins behind the bow shock.

Zoom In Zoom Out Reset image size

We now consider Model 2 (|B| = 3 μG, n_p = 0.095 cm⁻³), plotted in Figure 2 in the same format as Figure 1. Recall that the Alfvén Mach number at 1000 AU was barely above 1. Instead of a sharp transition, as in Model 1, there is now a gradual increase in plasma temperature (and density) beginning from about 650 AU for ∼200 AU before reaching a plateau. The gradient is readily apparent in the 1D cuts shown in Figure 4. This gradual change is very reminiscent of the self-consistent one-shock models discussed by Zank et al. (1996b). The 2D hydrodynamic models (i.e., a magnetic field was not included) for a subsonic LISM exhibited a gradual adiabatic compression of the LISM gas impinging on the heliosphere. Model 2 is rather different, however, in that the relative speed of the LISM far from the heliosphere is in fact marginally super-fast magnetosonic, and the transition to a sub-fast magnetosonic, and hence bow-shock-free configuration, is the result of heliospheric-created neutral H mediating the LISM. The precise physics underlying this behavior is examined in more detail below. The H-wall, illustrated in the bottom row of Figure 2, is again present, but now obviously broader and of lower amplitude than that of Model 1. The differences in the H-wall density, temperature, and velocity yield distinguishable differences in the Lyα absorption profiles along a given line of sight for Models 1 and 2. We exploit this possibility in the following section to constrain Models 1–3 against a suite of Lyα observations toward multiple nearby stars.

Zoom In Zoom Out Reset image size

For Model 3 (|B| = 4 μG, n_p = 0.048 cm⁻³), the far upstream Alfvén Mach number is ∼0.57. For these parameters, there is no possibility that a bow shock is present. As illustrated in Figure 3 (see also Figure 4), there is a very weak adiabatic compression of the LISM plasma as it piles up against the heliopause. The weak compression of the plasma in the presence of an oblique magnetic field contrasts with the gas dynamic one-shock model of Zank et al. (1996b). In the one-shock gas dynamic models, the plasma upstream of the heliopause experiences quite significant although adiabatic compression. The reason for the difference here is that the additional compression provided by the oblique magnetic field reduces the pressure contributed by the LISM plasma. Consequently, the MHD models do not exhibit the same degree of plasma compression and heating as the gas dynamic one-shock models (Zank et al. 1996b). This in turn reduces the effectiveness of secondary charge exchange of heliospheric neutral H and consequently results in the less effective formation of an H-wall abutting the heliopause. This can be seen in the bottom row of plots in Figure 3 which show a small amplitude, very broad H-wall. An interesting point is that a less oblique magnetic field would require greater compression of the plasma for a one-shock model and therefore a correspondingly larger H-wall even if |B| were identical. This may be an interesting parameterization to exploit if trying to constrain one-shock heliospheric models using Lyα absorption measurements.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As can be discerned from the 3D plots of Figure 1 and clearly illustrated in Figure 4 for the α-Cen sightline, Model 1 exhibits a plasma transition at ∼450 AU that, although reminiscent of a typical fast mode shock, is at least ∼50–75 AU wide. The thickness and smoothness of the transition is not an artifact of numerical inaccuracies, and it varies in thickness along different sightlines. Instead, the transition resembles a shock mediated by neutral H with a structure and thickness determined partially by the charge exchange length scale. The charge exchange mean free path in the LISM is approximately 70 AU (assuming a charge exchange cross section σ_c = 5 × 10⁻¹⁵ cm² and a total LISM number density of 0.2 cm⁻³). The overall compression ratio across the structure is ∼0.21/0.135 = 1.5 meaning that the transition is very weak. If the transition of Model 1 is interpreted as a weak bow shock, it raises interesting questions about the nature of the shock dissipation mechanism, discussed further in the conclusions section. By contrast, the plasma transition from the distant LISM to the outer heliosheath for Model 2 is quite different from Model 1. Model 2 exhibits a very broad transition, beginning at ∼650 AU with a slowly increasing density that reaches an approximate plateau or a small peak at ∼400 AU. Thereafter, the density decreases slowly until the heliopause, at which there is an abrupt jump down to the inner heliosheath density. The transition structure is certainly not due to numerical inaccuracies. Finally, the Model 3 density profile shows that there is very little density increase as the LISM plasma impinges on the heliosphere. Similar conclusions hold along all upwind directions, as well as the nose direction, which is discussed in the following section.

The plasma temperature cuts along the α-Cen sightline illustrated in Figure 4 for the three models show that there is an extended gradient from the heliopause into the LISM. As discussed above, this temperature gradient is quite unlike models of the global heliosphere that do not include neutral H self-consistently (see, e.g., Figure 5.4 in Zank 1999 for a typical temperature profile for a gas-dynamic-only simulation). The extended temperature gradient is a direct result of the secondary charge exchange of heliospheric neutral H. Note also that there is a temperature decrease from the HTS to the heliopause—this is due to charge exchange of cold primary LISM neutral H with hot heliosheath plasma, so reducing the temperature (and creating in turn hot heliospheric neutrals). The temperature gradient from the heliopause to ∼400 AU for all three cases is almost identical, indicating that the primary heating mechanism for the (sub-fast magnetosonic) LISM plasma is charge exchange with hot neutral H. The temperature profile for Model 1 shows that there is a modest temperature increase in the plasma as it crosses the bow shock. There is a similar gradual heating of the plasma for Model 2, starting at the transition region ∼650 AU. The temperature increases slowly but almost uniformly from the outermost boundary for Model 3. Evidently, secondary charge exchange dominates the heating of the LISM and, since the sound speed is proportional to T^1/2_p, must therefore play a role in determining whether a bow shock is needed locally for the solar-wind–LISM interaction.

Plotted in Figure 5 are the fast magnetosonic (solid lines) and Alfvén (dashed lines) Mach numbers M_f and M_A along the α-Cen sightline. The fast (and slow) magnetosonic speed, being a combination of the sound speed, Alfvén speed, and magnetic field obliquity, is sensitive to the increasing plasma temperature as the LISM flow approaches the heliosphere. In Model 1, the LISM plasma is clearly super-fast magnetosonic (and super-Alfvénic) until the bow wave or mediated shock, at which it transitions to a sub-fast magnetosonic flow. Immediately downstream of the transition, the plasma remains super-Alfvénic (as it should) for perhaps ∼20–30 AU. Note that the plasma density peaks immediately downstream of the broad bow wave (Figure 4) and decreases rapidly until M_A = 1, after which it decreases more slowly toward the heliopause. The far upstream flow for Model 2 is just barely super-fast magnetosonic and, like Model 1, the fast Mach number M_f decreases slowly as it approaches the heliosphere. As seen from the density profile, the transition extends from ∼650 AU to ∼400 AU and M = 1 at ∼700 AU. The increase in the LISM density exhibited by Model 2 occurs only after M_f = 1. In Model 1, the LISM density increase begins before M = 1 and continues beyond. The charge exchange mean free path for interstellar plasma is ∼40–50 AU in the vicinity of M_f = 1. As we discuss below, such a flow is analogous to flow through a nozzle where a supersonic flow can be decelerated smoothly to a subsonic flow through a critical point where M = 1. In this case, the "nozzle" is due to the momentum and energy source terms in the plasma equations. Thus, the 3 μG Model 2 does indeed represent a shock-free, smooth transition or a one-shock model despite the nominally super-fast magnetosonic state of the LISM. The final example, Model 3, is sub-fast magnetosonic throughout the LISM, and, like the other models, M_f decreases as the LISM plasma flows toward the heliosphere. In this case, M_f ≃ M_A in the LISM upwind of the heliopause.

Zoom In Zoom Out Reset image size

2.1. Structure of the Bow Wave Transition

Let us consider Model 2 more closely to illuminate why the transition is shock-free yet decelerates a super-fast magnetosonic flow to a sub-fast state. To illustrate, consider a somewhat idealized problem. Since the LISM flow is super-fast magnetosonic, we approximate the upwind flow as stationary and 1D. For this subsection, we focus our attention on the nose region for simplicity. A completely general analysis is presented in the Appendix. Note that the 1D figures of this section apply to the nose region of the solar-wind–LISM interaction only, unlike the other 1D cuts depicted elsewhere which are along the α-Cen sightline.

We assume an idealized configuration in which the LISM flow velocity defines the x-coordinate and the magnetic field is taken to be orthogonal to the LISM flow, i.e.,

$$\begin{equation} {\bf U} = (U_x,0,0), \quad {\bf B} = (0, B_y, B_z). \end{equation} \tag{ 1 }$$

The 1D MHD equations with momentum and energy source terms Q_m and Q_e, respectively, are then

$$\begin{eqnarray} &&\frac{d}{dx} \left(\rho U_x \right) = 0 \; \Longrightarrow \; \rho U_x = \alpha = \mbox{const.} ; \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} &&\rho U_x \frac{dU_x}{dx} + \frac{\it d P}{dx} + \frac{1}{4\pi } \left(B_y \frac{\it d B_y}{dx} + B_z \frac{\it d B_z}{dx} \right) = Q_{mx} ; \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} &&\frac{d}{dx} \left(\frac{1}{2} \rho U_x^3 + \frac{\gamma }{\gamma - 1} U_x P + \frac{1}{4\pi } U_x B^2 \right) = Q_e; \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} &&U_x B_y = \mbox{const.} ; \; U_x B_z = \mbox{const.} ; \end{eqnarray} \tag{ 5 }$$

where ρ, U_x, P, and γ denote the LISM plasma density, velocity, pressure, and adiabatic index respectively and B² = B²_y + B_z². Provided the LISM flow remains smooth, we may rewrite energy conservation in terms of the adiabatic energy equation,

$$\begin{equation} U_x \frac{\it d P}{dx} + \gamma P \frac{d U_x}{dx} = \left(\gamma - 1 \right) \left(Q_e - U_x Q_{mx} \right). \end{equation} \tag{ 6 }$$

By introducing the sound speed C_s, the Alfvén speed V_A, and the fast magnetosonic speed for this configuration through

$$\begin{equation} C_s^2 = \frac{\gamma P}{\rho } ; \; V_A^2 = \frac{B^2}{4\pi \rho } ; \; V_f^2 = C_s^2 + V_A^2, \end{equation} \tag{ 7 }$$

and the corresponding fast magnetosonic M_f and Alfvén M_A Mach numbers,

$$\begin{equation} M_f^2 = \frac{U^2}{V_f^2} ; \; M_A^2 = \frac{U^2}{V_A^2}, \end{equation} \tag{ 8 }$$

some straightforward manipulation of Equations (2)–(6) yields

$$\begin{eqnarray} &&\frac{\alpha V_f^2 }{\gamma - 1} \left(M_f^2 - 1 \right) \frac{1}{U_x} \frac{d U_x}{dx} = Q_e - \frac{\gamma }{\gamma - 1} U_x Q_{mx} ; \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} &&\frac{1}{\gamma - 1} \frac{M_f^2 - 1}{M_f^2} \frac{d P}{dx} = \frac{ U_x^2 - V_A^2 }{U_x^3} Q_e \nonumber \\ &&\mbox{} - \frac{U_x^2 - V_A^2 - \frac{1}{\gamma - 1} C_s^2 }{U_x^3} U_x Q_{mx}. \end{eqnarray} \tag{ 10 }$$

An important point to note about Equations (9) and (10) is that the source terms Q_m and Q_e are non-zero only because of the secondary charge exchange of fast and hot heliospheric neutral H. In the distant LISM⁶ where the LISM plasma and neutral H are equilibrated, the source terms are identically zero since T_p(> 1000 AU) = T_H(> 1000 AU) and U(> 1000 AU) = U_H(> 1000 AU) (where U_H is the neutral H flow speed). In principle, by using

$$\frac{1}{M_f^2} \frac{d}{dx} \left(M_f^2 \right) = \left(1 + 2 \frac{V_f^2}{V_A^2} \right) \frac{1}{U_x} \frac{d U_x}{dx} - \frac{V_f^2}{C_s^2} \frac{1}{P} \frac{d P}{dx},$$

one can combine Equations (9) and (10) as a single ordinary differential equation in M²_f, but this is unnecessarily laborious and we use Equations (9) and (10) directly. For a critical point to exist, both the left-hand side and right-hand side of Equations (9) and (10) must be zero simultaneously. Obviously, the left-hand side of both Equations (9) and (10) vanishes for M²_f = 1. Thus, for a critical point to exist, Equation (9) implies

$$\begin{equation} Q_e = \frac{\gamma }{\gamma - 1} U Q_m \end{equation} \tag{ 11 }$$

and Equation (10) requires

$$\begin{equation} \left(U_x^2 - V_A^2 \right) Q_e = \left(U_x^2 - V_A^2 - \frac{1}{\gamma - 1} C_s^2 \right) U_x Q_{mx}. \end{equation} \tag{ 12 }$$

On using Equation (11), Equation (12) reduces to

$$\begin{equation} U_x^2 = C_s^2 + V_A^2, \end{equation} \tag{ 13 }$$

which is of course nothing more than M²_f = U_x²/V²_f = 1. Consequently, a critical point exists for M²_f = 1 simultaneously with Equation (11). Unfortunately, given the complexity of the source terms Q_m and Q_e (see Pauls et al. 1995; Heerikhuisen et al. 2008), it is not trivial to determine the nature of the critical point analytically. However, given the smooth solutions illustrated in Figure 5, the critical point would appear to be a saddle point, ensuring that the heliospheric–LISM flow transition can possess a smooth decelerating structure. Figure 5 is appropriate to an oblique magnetic field configuration rather than the perpendicular case discussed here. The Appendix provides a fully general analysis that would be appropriate to the α-Cen sightline.

Plotted in Figure 6 is the fast magnetosonic Mach number for each of the three models, analogous to the plots of M_f and M_A of Figure 5 but now along the nose direction. The differences in the location of the flow transition between the nose direction and the α-Cen sightline are apparent. For Model 1, the bow shock is located at ∼360 AU with a width of ∼40 AU and M_f = 1 at ∼330 AU. Model 2 begins its transition from a super-fast magnetosonic state to one that is sub-fast at ∼600 AU and has a width of ∼200 AU, and M_f = 1 at ∼550 AU.

Zoom In Zoom Out Reset image size

We plot in Figure 7 the radial evolution of γ/(γ − 1)U_xQ_mx and Q_e. The top plot of Figure 7 shows the source terms Q_e (blue curve) and γ/(γ − 1)U_xQ_mx (red curve) as a function of heliocentric distance along the LISM flow line for Model 1, and the bottom plot is for Model 2. Also plotted is the spatial location of the M_f = 1 point. Although the evaluated source terms Q_e and Q_mx are comparatively noisy because of their kinetic origin, we see that the relation (11) holds where M_f = 1 for both Models 1 and 2, i.e., the super-fast magnetosonic LISM flow is decelerated smoothly through the point M_f = 1, Q_e = γ/(γ − 1)U_xQ_mx for both the broad bow wave or shock-like transition (Model 1) and the smooth transition (Model 2). The important, and remarkable, conclusion to draw from Figure 7 is that hot and fast neutral H created self-consistently from heliospheric plasma leads to a smoothing of the heliospheric–LISM flow transition, whether it is expected to be a shock (Model 1) or a super-fast magnetosonic–sub-fast magnetosonic "nozzle-like" (Model 2) transition. A sufficiently super-fast LISM flow will still require a shock transition but the structure will in part be determined by fast and hot neutral H. Charge exchange of fast and hot neutral H therefore provides part of the dissipation mechanism for the weak bow shock in the 2 μG case since expression (11) holds at M_f = 1 in the middle of the transition. The structure may therefore be viewed as a bow wave rather than a shock wave. By contrast, for Model 2, a barely super-fast magnetosonic flow can be decelerated completely by fast and hot neutral H experiencing charge exchange in the LISM. The transition occurs over ∼200 AU, which is ∼2 times the charge exchange mean free path in the LISM, which is consistent with the transition length scale if neutral H were the primary dissipation mechanism. We regard the transition of Model 2 as a smoothly decelerated super-fast–sub-fast magnetosonic transition mediated by fast and hot neutral H, and we will describe the transition as "bow-shock-free" or as a bow wave. The nonlinear feedback of LISM and heliosphere can therefore have a profound influence on the structure of the heliosphere.

Zoom In Zoom Out Reset image size

In closing, we note there is some correspondence between the mediation of the LISM bow shock discussed here and cometary shocks. Biermann et al. (1967) and Wallis (1973) considered the interaction of the solar wind plasma as it interacted with the neutral gas of an active cometary coma. By introducing mass-loading source terms into the plasma equations, they found that there were two possible flow transitions, (1) a gradual deceleration of the upstream solar wind plasma flow followed by a sharp MHD shock transition or (2) the entire flow was mediated and a smooth super-fast magnetosonic–subfast magnetosonic transition through a critical point was possible. The existence of bow shocks or bow waves has since been confirmed by all spacecraft observations (Giotto, Vega-1 and 2, Suisei, Sakigake at P/Halley and P/Grigg-Skjellerup, and International Cometary Explorer at P/Giacobini-Zinner), and the observations have been reviewed in detail by, e.g., Szegö et al. (2000). Although somewhat different in terms of the fundamental role played here by hot and fast heliospheric neutral H, an analogy can nonetheless be drawn between cometary shocks and the broadened bow shock and the smooth super-fast magnetosonic–sub-fast magnetosonic transition exhibited by the solar-wind–LISM interaction models of Section 2.

2.2. Neutral H Distributions

Consider now the neutral H distributions for the three models. Figure 8 shows the neutral H number density (top), temperature (middle), and radial velocity (bottom) profiles along the α-Cen sightline. Already evident from the global density distributions illustrated in Figures 1–3, the structure of the H-wall is quite different for each model. The largest amplitude and narrowest wall is that of Model 1, and the broadest and lowest amplitude is that of Model 3. Since the column density is an important factor in determining the Lyα profile along a line of sight, it is not a priori obvious which of the three models is likely to provide a better fit to observations. In Model 1, the H number density increases just downstream of the bow wave due to charge exchange with the slowed, diverted, and heated LISM protons, and the density of 0.22 cm⁻³ at 1000 AU is compressed to ∼0.36 cm⁻³. The H-wall has a width of ∼200 AU. The bow-shock-free Model 2 also produces a fairly strongly compressed H-wall (the H density increasing from ∼0.195 to ∼0.27 cm⁻³) but with a width of ∼400 AU. By contrast, Model 3 produces a very broad, low-amplitude structure, which, although it experiences considerable filtration as it enters the heliosphere, certainly cannot be described as in an H-wall. We note that the weak H-wall for Model 3 contrasts quite markedly with the self-consistent gas dynamic–neutral H one-shock models (Zank et al. 1996b), and as discussed above, this distinction is due in part to the orientation of the LISM magnetic field (e.g., a configuration in which the LISM magnetic field was more parallel to the upwind flow vector might produce a larger H-wall).

Zoom In Zoom Out Reset image size

The middle panel of Figure 8 shows that the three temperature profiles are distinguished only in the transition regions. Between ∼200 and ∼400 AU, the H temperature in all three models is identical. Between ∼400 and ∼650 AU, the H temperature of Model 2 exceeds that of Model 1, and from ∼500 AU both Models 1 and 2 have lower H temperatures than that of Model 3. This behavior is due to the shorter charge exchange mean free paths in the higher density regions of Models 1 and 2, which allows fast and hot heliospheric neutral H to penetrate to greater distances into the LISM for Model 3.

Finally, the bottom panel of Figure 8 shows the velocity moment of the H velocity distributions; for this case, we plot the radial bulk H velocity along the α-Cen line of sight. The velocity profiles reflect the differences in the density profiles exhibited by all three models.

The neutral H number density, temperature, and velocity all play an important role in determining the absorption of Lyα by the H-wall and LISM neutral H. This can be seen from the expression for the optical depth of the H-wall (e.g., Gayley et al. 1997),

$$\begin{equation} \tau (\lambda) = \frac{7.5 \times 10^{-13} }{b_H } N_{H} e^{-(247 \lambda - U_H)^2 / b_H^2 }, \end{equation} \tag{ 14 }$$

where U_H is the neutral H velocity km s⁻¹, λ is the wavelength in Angstroms from line center in the heliocentric rest frame, and N_H is the H-wall column density. The term b_H is the "Doppler parameter" and includes the line broadening due to all mass motions, both thermal (mass-dependent) and non-thermal (mass-independent), the latter of which might be called turbulent motion, and is measured in km s⁻¹. Note that both H i and D I are dominated by thermal broadening in the ISM, and presumably even more so in the hotter heliosphere. Although the neutral H column density varies across the three models, the velocity and temperature differences can have the largest effect on optical depth.

In closing this section, we present in Figure 9 the reduced neutral H distribution function. Figure 9 shows the radial velocity distribution function at 300 AU along the α-Cen sightline for the three models considered. The 1D distribution functions very nicely illustrate the basic physics that we have described above. Overplotted on the figure is the assumed neutral H Maxwellian distribution at 1000 AU (in all cases almost identical since the temperature is equal and the number densities are only slightly different). Because the LISM plasma has been heated by secondary charge exchange, the neutral H core distribution has broadened considerably in all cases. The hot neutrals created in the inner heliosheath form an extended tail to the Maxwellian-like core from about ∼70 km s⁻¹ to ∼300 km s⁻¹. The fast neutrals created by charge exchange with the supersonic solar wind manifest themselves as the broad bump at the fast end of the distribution, from ∼300 km s⁻¹ to ∼470 km s⁻¹. It is the neutral H particles in the extended (i.e., outflow) part of the distribution function that leads to the physical effects associated with the bow shock-transition region discussed here.

Zoom In Zoom Out Reset image size

It was recognized that a very powerful diagnostic tool for constraining the structure of the global heliosphere was the absorption of stellar Lyα light by interstellar neutral H observed by the Hubble Space Telescope (Linsky & Wood 1996; Gayley et al. 1997). The absorption signatures in stellar Lyα will be mediated by the neutral H created in the solar-wind–LISM boundaries and the extent of this mediation can provide constraints on theoretical models of the global heliosphere. The first effort to constrain global heliospheric structure based on the self-consistent inclusion of neutral H was that of Gayley et al. (1997). Although using a simple 2D hydrodynamic multi-fluid neutral H model (Zank et al. 1996b), the comparison of the synthetic Lyα absorption profiles with those observed along the α-Cen sightline demonstrated convincingly that (1) the H-wall did indeed provide an observable signature, confirming the claim of Linsky & Wood (1996) that they had observed the H-wall, and (2) in principle, Lyα absorption profiles could potentially distinguish between a two-shock and a one-shock heliospheric model. Indeed, Gayley et al. suggested that a barely subsonic model of the LISM could explain the additional Lyα absorption reported by Linsky & Wood. Gayley et al. (1997) did caution, however, that a more extensive parameter study was needed, and this has been the subject of numerous investigations, based largely on hydrodynamic models, along many different sightlines (Wood et al. 2000, 2005b), and even for different stars (Müller et al. 2001; Wood et al. 2001, 2002, 2005a). Synthetic Lyα models derived from MHD-kinetic models were considered by Wood et al. (2007a, 2007b), based on the models discussed by Izmodenov et al. (2005a) and Izmodenov & Alexashov (2006). The heliospheric models from which Wood et al. (2007a, 2007b) derived the Lyα profiles are limited in including only the LISM magnetic field and not the interplanetary magnetic field. As discussed by Pogorelov et al. (2008, 2011), Heerikhuisen et al. (2006), and Washimi et al. (2007, 2011), the interplanetary magnetic field plays an important role in determining global heliospheric structure. Based on the 3D MHD-kinetic models discussed in Section 2, we synthesize Lyα absorption profiles along four sightlines, and compare these against those observed.

The precise methodology for synthesizing Lyα profiles from heliospheric models is described in detail elsewhere (Gayley et al. 1997; Wood et al. 2000, 2005b, 2007a). However, we emphasize that we compute here absorption profiles directly from the kinetically computed velocity distributions, such as those illustrated in Figure 9, and not from quantities computed by taking moments of them (see Wood et al. 2000). We consider four representative lines of sight and defer to another study the 20 or so sightlines that are available for comparison. Details regarding the stars, their spectral type, ecliptic coordinates, and distance are listed in Table 2. With the Sun at the origin, the angle θ that defines a sightline is measured with respect to the far upwind interstellar flow direction. The most directly upwind sightline that we model and measure is that toward 36 Oph, with θ = 9°. This and α Cen, with θ = 51°, are the most reliable observations that we possess since (1) the H-wall absorption toward these stars (relative to the LISM absorption) is the strongest and (2) the spectral resolution of the 36 Oph and α Cen data is superior to other lines of sight (Wood et al. 2007a). We also consider a side sightline because this direction in principle should identify the distinct spatial distribution and extent of the H-wall for the three models considered, especially Model 3. The observation we choose to use is the DK UMa line of sight with θ = 116°. Note that no heliospheric Lyα absorption is actually detected toward DK UMa, or any other sidewind line of sight for that matter, but the DK UMa data are still useful for providing an upper limit for the amount of absorption in that direction. Finally, we choose the most reliable observation in the downwind direction (i.e., along the heliotail direction), which is the χ¹ Ori line of sight, θ = 170° (Wood et al. 2007b).

Shown in Figure 10 is the Lyα absorption observed along the four lines of sight (the solid black curves) together with the absorption curves predicted by the three models (the red curves correspond to the Model 1 two-shock case, the blue curves to the Model 2 bow-shock-free super-fast magnetosonic case, and the green curves to the one-shock sub-fast magnetosonic case). Figure 10 shows the red side of the absorption profile only since this is where heliospheric absorption will be present. The dotted lines of Figure 10 identify the absorption that is expected from the LISM only. The solid colored lines result from the combined effect of absorption from the heliospheric models (up to 1000 AU) and the interstellar absorption beyond (Gayley et al. 1997). A reconstructed stellar Lyα profile for the particular star being observed is prescribed (e.g., Wood et al. 2005b). In the upwind direction, Model 3 predicts too little absorption along both the 36 Oph and α Cen sightlines, and yet a little too much absorption along the DK UMa sightline. By contrast, Models 1 and 2 fit the observed data much better, with the α Cen data being noticeably better fitted by the Model 2 profile. There is little to distinguish between Models 1 and 2 along the 36 Oph and DK UMa sightlines, and Model 2 may be marginally better along the χ¹ Ori sightline. We need to be cautious in interpreting the Lyα absorption toward χ¹ Ori, i.e., in the heliotail direction, for two reasons. First, the models themselves have a heliocentric extent of 1000 AU, meaning that the heliotail is truncated in the simulations, and consequently heliospheric absorption may be underestimated along a sightline that intersects the heliotail. This requires that ideally we should use a very long (and very time consuming) heliotail simulation. Second, an effect that may partially offset the first point is that Schwadron et al. (2011) have interpreted IBEX observations as suggesting that the heliotail is deflected from the LISM flow vector direction (see Figure 4 of McComas et al. 2012). Such a deflected heliotail is present in most MHD models and is seen clearly in Figures 1–3, although the deflection is not as strong as suggested in Schwadron et al. (2011).⁷ Thus, depending on the orientation of the heliotail, a θ = 170° sightline may not have an extremely long path through the heliotail region. Furthermore, as illustrated in Figures 1–3, especially for Models 2 and 3, the cooling of the heliotail by neutral H leads to quite a rapid pinching of the tail. Consequently, there may be reasons to accept the χ¹ Ori modeled absorption profiles as reasonable.

Zoom In Zoom Out Reset image size

The strong conclusions that can be drawn from Figure 10 are (1) the highly sub-fast magnetosonic LISM case, Model 3 with |B| = 4 μG, cannot account consistently for the observed Lyα absorption along different sightlines and (2) a marginally super-fast magnetosonic LISM can account consistently for the observed absorption along all the sightlines considered here. A weaker conclusion is that the Model 2 3 μG model may be marginally better than the two-shock 2 μG Model 1, particularly if weighted by the α Cen observations. This case would certainly be consistent with the interpretation advanced by McComas et al. (2012) of a bow-shock-free configuration of the heliosphere-LISM interaction. The interpretation is, however, more subtle since the absence of a bow shock does not necessarily imply a sub-fast magnetosonic flow but instead a super-fast magnetosonic–sub-fast magnetosonic transition through a critical point. To ensure that the H-wall is sufficiently localized spatially with a relatively large amplitude, the LISM plasma flow must be reasonably efficiently decelerated via charge exchange. It is clear from our analysis that the Lyα absorption data are consistent with a very slightly super-fast magnetosonic LISM wind, regardless of the complex physics involved in such a marginal M ∼ 1 shock transition.

Based on a set of observations from IBEX, McComas et al. (2012) concluded that the current best value for the LISM flow speed is 23.2 km s⁻¹. It was also suggested that the LISM flow might be sub-fast magnetosonic and that the heliosphere might not therefore possess a bow shock. The IBEX observations raise several challenges to theory and models.

• 1.  
  Is a LISM model that is barely super- or sub-fast magnetosonic consistent with Lyα absorption measurements along multiple sightlines, since the interpretation of the Lyα observations relies critically on the additional absorption provided by the H-wall.
• 2.  
  If the LISM flow is weakly super-fast magnetosonic and a shock transition of some kind is necessary, what then is the basic dissipation mechanism, and hence structure, of the shock? Weak collisionless shocks in the solar wind are thought to be laminar (e.g., Formisano 1977) but in a partially ionized plasma such as the LISM, does charge exchange play a role in the shock dissipation process?

To address these questions, we considered three self-consistently coupled MHD plasma–kinetic H models using current LISM parameters and varied the magnetic field strength (using values of 2, 3, and 4 μG) and the plasma and neutral number density (Table 1). The 2 μG Model 1 yielded the canonical two-shock structure for the heliosphere, with a very broad and weak bow shock (compression ratio ∼1.5). By contrast, the 3 μG Model 2 was weakly super-fast magnetosonic far upwind but the transition from a super-fast magnetosonic to a sub-fast magnetosonic flow was accomplished smoothly over a distance of ∼200 AU. The 4 μG Model 3 was sub-fast magnetosonic upwind of the heliopause and no bow shock or transition of any sort was present.

An analysis of the MHD equations with charge exchange source terms showed that a super-fast magnetosonic flow admitted a critical point in the flow when M = 1 and Q_e = γ/(γ − 1)UQ_m simultaneously. We found that for both Model 1 and Model 2, the LISM flow passed through this point in transitioning from a super-fast magnetosonic to a sub-fast magnetosonic state. Thus, we may interpret this result as indicating that fast and hot neutral H created in the heliosphere mediates the bow shock via charge exchange. In the two-shock case, the mediation is only partial since the flow is sufficiently super-fast magnetosonic that an additional dissipation mechanism is needed. By contrast, for the weakly super-fast magnetosonic Model 2 case, fast and hot neutral H completely mediates the shock transition, and imposes the charge exchange length scale on the transition that takes the super-fast magnetosonic upstream state to a sub-fast magnetosonic state. We regard this ∼200 AU thick structure as a shock-free transition or bow wave. The possibility that charge exchange can provide dissipation at a shock and mediate the structure has begun to be explored (Dastgeer et al. 2004; Borovikov et al. 2008; Blasi et al. 2012).

Both of the super-fast magnetosonic LISM two-shock and shock-free Models 1 and 2, respectively, produce an H-wall of sufficient column depth to account for the Lyα observations along the α Cen, 36 Oph, DK UMa, and χ¹ Ori sightlines. The sub-fast magnetosonic Model 3 possesses a small H-wall that is unable to account for the Lyα observations. The observations may marginally favor the 3 μG shock-free Model 2. However, further modeling and more Lyα comparisons with other sightlines is needed to definitively establish this result. We are nonetheless left with a tantalizing question: Has IBEX discovered a new class of shock wave mediated by interstellar neutral H?

We acknowledge the partial support of NASA grants NNX08AJ33G, Subaward 37102-2, NNX09AG70G, NNX09AG63G, NNX09AJ79G, NNG05EC85C, Subcontract A991132BT, NNX09AP74A, NNX10AE46G, NNX09AW45G, and NNH09AM47I, and NSF grant ATM-0904007. E.Z. acknowledges the support of an NESSF grant NNX11AP91H. This work was carried out as a part of the IBEX mission, which is part of NASA's Explorer Program.

The analysis of Section 2.1 was restricted to the simplified case of a LISM magnetic field orthogonal to the LISM flow vector. Here we consider a more general LISM flow and magnetic field orientation although still assuming a steady-state 1D model. Thus, we assume

$$\begin{equation} {\bf U} = (U_x, U_y, U_z ), \quad {\bf B} = (B_x, B_y, B_z ), \end{equation} \tag{ A1 }$$

and d/dy = 0 = d/dz. Conservation of mass yields

$$\begin{equation} \frac{d}{dx} \left(\rho U_x \right) = 0 \Longrightarrow \rho U_x = \alpha = \mbox{const.} \end{equation} \tag{ A2 }$$

Conservation of momentum,

$$\begin{equation} \rho {\bf U} \cdot \nabla {\bf U} + \nabla P - \frac{1}{4\pi } \left(\nabla \times {\bf B} \right) \times {\bf B} = {\bf Q}_m, \end{equation} \tag{ A3 }$$

yields

$$\begin{eqnarray} \alpha \frac{d U_x}{dx} + \frac{\it d P }{dx} + \frac{1}{4 \pi } \left(B_y \frac{{\it d B}_y}{dx} + B_z \frac{{\it d B}_z}{dx} \right) &=& Q_{mx} ; \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} \alpha \frac{d U_y}{dx} - \frac{1}{4 \pi } B_x \frac{{\it d B}_y}{dx} &=& Q_{my} ; \end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray} \alpha \frac{d U_z}{dx} - \frac{1}{4 \pi } B_x \frac{{\it d B}_z}{dx} &=& Q_{mz}, \end{eqnarray} \tag{ A6 }$$

while

$$\begin{equation} \nabla \cdot {\bf B} = 0 \Longrightarrow B_x = \mbox{const.} \end{equation} \tag{ A7 }$$

It is again convenient to work with the adiabatic form of the energy equation since we are not considering shocks:

$$\begin{equation} U_x \frac{\it d P}{dx} + \gamma P \frac{d U_x}{dx} = (\gamma - 1) \left(Q_e - {\bf U} \cdot {\bf Q}_m \right). \end{equation} \tag{ A8 }$$

Finally, ∇ × (U × B) = 0 yields

$$\begin{equation} \frac{U_x B_y}{B_x} = U_y ; \; \frac{U_x B_z}{B_x} = U_z. \end{equation} \tag{ A9 }$$

By eliminating P and B in (A4), we can rewrite Equation (A4) as

$$\begin{eqnarray} && \frac{\alpha }{U_x^2} \big(U_x^4 - \big(C_s^2 + V_A^2 \big) U_x^2 + C_s^2 V_x^2 \big) \frac{1}{U_x} \frac{d U_x}{dx} = -\left(1 - \frac{V_x^2 }{U_x^2} \right) \nonumber \\ && \quad\quad\times \left((\gamma - 1) Q_e - \gamma {\bf U} \cdot {\bf Q}_m \right) - {\bf U} \cdot {\bf Q}_m + U_x Q_{mx}. \end{eqnarray} \tag{ A10 }$$

The left-hand side of Equation (A10) introduces the fast and slow magnetosonic speeds

$$\begin{equation} V_{f,s} = V_{\pm } = \frac{1}{2} \left(C_s^2 + V_A^2 \pm \sqrt{ \left(C_s^2 + V_A^2 \right)^2 - 4C_s^2 V_x^2 } \right), \end{equation} \tag{ A11 }$$

where V²_A = B²/(4πρ) and B² = B²_x + B_y² + B²_z. This allows us to factor (A10) as

$$\begin{eqnarray} && \alpha \frac{V_f^2 V_s^2}{U_x^2} \big(U_x^2 - V_f^2 \big) \big(U_x^2 - V_s^2 \big) \frac{1}{U_x} \frac{d U_x}{dx} = -\left(1 - \frac{V_x^2 }{U_x^2} \right) \nonumber \\ && \quad\quad \times \left((\gamma - 1) Q_e - \gamma {\bf U} \cdot {\bf Q}_m \right) - {\bf U} \cdot {\bf Q}_m + U_x Q_{mx}. \quad\quad\quad \end{eqnarray} \tag{ A12 }$$

Note that V_x = 0 corresponds to the case discussed in Section 2.1, i.e., V²_f = C_s² + V²_A and V²_s = 0. For B_y = 0 = B_z, V²_f = C_s² and V²_s = V_x².

Similarly, Equation (A8) can be expressed as

$$\begin{eqnarray} && \frac{1}{ \gamma - 1} \frac{V_f^2 V_s^2}{U_x} \big(M_f^2 - 1 \big) \big(M_s^2 - 1 \big) \frac{\it d P}{dx} = \big(U_x^2 - V_A^2 \big) Q_e \nonumber \\ && \quad\quad - \left(U_x^2 - V_A^2 - \frac{1}{\gamma - 1} C_s^2 \frac{V_x^2}{U_x^2} \right) {\bf U} \cdot {\bf Q}_m \nonumber\\ && \quad\quad - \frac{1}{\gamma - 1} C_s^2 U_x Q_{mx}. \end{eqnarray} \tag{ A13 }$$

For a critical point to exist in the flow, we require that the left- and right-hand sides of Equations (A12) and (A13) be zero simultaneously. Thus,

• 1.  
  the left-hand side of both Equations (A12) and (A13) is zero if M²_f = 1 (or M²_s = 1),
• 2.  
  the right-hand side of Equation (A12) is zero whenever

  $$\begin{equation} Q_e = \frac{\gamma }{\gamma - 1} {\bf U} \cdot {\bf Q}_m - \frac{1}{\gamma - 1} \frac{ {\bf U} \cdot {\bf Q}_m - U_x Q_{mx} }{1 - V_x^2 /U_x^2 } , \end{equation} \tag{ A14 }$$

  and
• 3.  
  the right-hand side of Equation (A13) is zero whenever

  $$\begin{eqnarray} && \left(U_x^2 - V_A^2 \right) Q_e - \left(U_x^2 - V_A^2 - \frac{1}{\gamma - 1} C_s^2 \frac{V_x^2}{U_x^2} \right) {\bf U} \cdot {\bf Q}_m \nonumber\\ && \quad\quad\quad\quad\quad\quad - \frac{1}{\gamma - 1} C_s^2 U_x Q_{mx} = 0. \end{eqnarray} \tag{ A15 }$$

On using Equation (A14) to evaluate Q_e in Equation (A15), we find that the left-hand side of Equation (A15) can be expressed as

$$\begin{eqnarray} && \left(U_x^4 - \left(C_s^2 + V_A^2 \right) U_x^2 + C_s^2 V_x^2 \right) \left(U_x Q_{mx} - \frac{V_x^2}{U_x^2} {\bf U} \cdot {\bf Q}_m \right) \nonumber \\ && \quad\quad = \left(U_x^2 - V_f^2 \right) \left(U_x^2 - V_s^2 \right) \left(U_x Q_{mx} - \frac{V_x^2}{U_x^2} {\bf U} \cdot {\bf Q}_m \right), \nonumber\\ \end{eqnarray} \tag{ A16 }$$

which is identically 0 when U²_x = V_f². Consequently, the general condition for a critical point in the flow is M_f = 1 and Equation (A14) simultaneously.