ASYMPTOTIC STEADY-STATE SOLUTION TO A BOW SHOCK WITH AN INFINITE MACH NUMBER

Bow shocks occur when a supersonic flow encounters an obstacle. Prominent examples from astrophysics include planetary bow shocks (Treumann & Jaroschek 2008), inundation of a dense molecular cloud by a supernova shock wave (McKee & Cowie 1975), and bow shocks around stars that move at a supersonic velocity relative to the interstellar medium (ISM) (Noriega-Crespo et al. 1997).

The last decades have yielded a plethora of analytic and numerical works on hydrodynamic bow shocks (as opposed to MHD bow shocks). Farris & Russell (1994) studied the different factors affecting the distance between the obstacle and the nose of the shock front. Wilkin (1996) considered a radiative (fast cooling) bow shock that forms when an object moves relative to a uniform medium, while also emitting an isotropic wind. Cantó & Raga (1998) considered a radiative bow shock that forms when a solid sphere moves relative to a uniform medium. Verigin et al. (2003) developed a semi-analytic theory for the shape of an adiabatic bow shock around a parabolic obstacle. Schulreich & Breitschwerdt (2011) developed a topological theory in order to solve the inverse problem: determining the shape of the obstacle from the shape of the shock. Brighenti & D'Ercole (1995) performed 2D simulations for the bow shock formed around a wind-emitting red supergiant (RSG) moving relative to the ISM. Mohamed et al. (2012) performed similar simulation in 3D for Betelgeuse. Miceli et al. (2006) performed simulations for the bow shock that forms when a supernova shock wave sweeps a dense cloud. Comeron & Kaper (1998) performed simulations of an OB star moving through the ISM. van Marle et al. (2011) also simulated the bow shock around a moving RSG, but included the effect of dust in the calculation.

Despite the work mentioned above, there remains an unexplored niche in parameter space. That niche is an adiabatic shock of infinite Mach number, and it is the focus of this work. We assume that the incoming matter is cold, so for any finite velocity its Mach number would be infinite. The asymptotic shape far from the obstacle is qualitatively different from the case of a finite Mach number M. In the latter, far away from the obstacle, the size of the obstacle becomes irrelevant, and the shock front coincides with the Mach cone, i.e., a cone with an opening angle $\alpha ={\sin }^{-1}\tfrac{1}{M}$ (Landau & Lifshitz 1987, p. 467). When the Mach number tends to infinity, the opening angle of the Mach cone tends to zero. Hence the bow shock in that case cannot be straight, and the curvature must still be significant at large distances from the obstacle. In this paper we show that the shape of the shock is a parabola. We develop the mathematical theory to describe the shape of the shock and the profile of the hydrodynamic variables within, and then verify our supposition using a numerical simulation.

In astrophysical bow shocks, the wind of the external medium brakes by means other than direct collision with the compact object (a star, planet etc.). The braking mechanism could be magnetic fields or winds emitted from the compact object. The effective size of the obstacle is determined by the condition that the inner and outer pressures balance. In the case of a magnetic obstacle the effective radius of the obstacle is where the magnetic pressure balances the ram pressure of the outside wind:

$$\begin{eqnarray}&&{r}_{m}=\sqrt[6]{{\mu }^{2}/{\rho }_{o}{v}_{o}}\end{eqnarray} \tag{ 1 }$$

where μ is the magnetic dipole moment, ${\rho }_{o}$ is the density of the external medium, and v_o is the relative velocity between the obstacle and the external medium. In the case of a wind obstacle, the stand-off distance is given by the balance of ram pressures of the inner and outer winds:

$$\begin{eqnarray}&&{r}_{w}=\sqrt{\dot{M}{v}_{i}/{\rho }_{o}{v}_{o}^{2}}\end{eqnarray} \tag{ 2 }$$

where $\dot{M}$ is the mass loss rate from the inner wind and v_i is its velocity (relative to the object emitting it). In the latter case the contact discontinuity between the inner and outer winds serves as an effective rigid obstacle, because energy and matter do not cross that surface. Moreover, at increasing distances from the object, the fraction of enclosed mass and energy from the inner wind disappears.

This paper is organized as follows. In Section 2 we describe the complete mathematical formulation. In Section 3 we present validation of our analytic results with a numerical simulation. In Section 4 we calculate the drag force of such a shock. In Section 5 we present a manifestation of this solution in an observation. Finally, in Section 6, we discuss the results.

In this section we develop the mathematical equation to describe the shape of the bow shock and the profiles of the hydrodynamic variables inside it. We begin with an intuitive explanation for the asymptotic shape of the shock, proceed by giving a complete description of the steady-state interior of the shock using parabolic coordinates, and end with a discussion of the growth of the bow shock.

2.1. Intuitive Explanation for the Shock's Shape

2.2. Steady State

Spherical parabolic coordinates are given by

$$\begin{eqnarray}&&x=\sigma \tau \cos \phi \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&y=\sigma \tau \sin \phi \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&z=\displaystyle \frac{1}{2}({\tau }^{2}-{\sigma }^{2}).\end{eqnarray} \tag{ 5 }$$

In these coordinates, the steady-state, azimuthally symmetric hydrodynamics equations take the following form. The conservation of mass is

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial \sigma }(\rho \sqrt{{\sigma }^{2}+{\tau }^{2}}\sigma \tau {v}_{\sigma })+\displaystyle \frac{\partial }{\partial \tau }(\rho \sqrt{{\tau }^{2}+{\sigma }^{2}}\sigma \tau {v}_{\tau })=0.\end{eqnarray} \tag{ 6 }$$

The conservation of entropy is

$$\begin{eqnarray}&&{v}_{\sigma }\displaystyle \frac{\partial s}{\partial \sigma }+{v}_{\tau }\displaystyle \frac{\partial s}{\partial \tau }=0\end{eqnarray} \tag{ 7 }$$

where $s=\mathrm{ln}p-\gamma \mathrm{ln}\rho$ is the specific entropy. Since there is no dissipation inside the shocked region, this equation is equivalent to energy conservation (Landau & Lifshitz 1987, p. 467). The conservation of momentum in each direction is

$$\begin{eqnarray}&&{v}_{\sigma }\displaystyle \frac{\partial {v}_{\sigma }}{\partial \sigma }+{v}_{\tau }\displaystyle \frac{\partial {v}_{\sigma }}{\partial \tau }+\displaystyle \frac{{v}_{\tau }({v}_{\sigma }\tau -{v}_{\tau }\sigma )}{{\tau }^{2}+{\sigma }^{2}}+\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial p}{\partial \sigma }=0\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{v}_{\sigma }\displaystyle \frac{\partial {v}_{\tau }}{\partial \sigma }+{v}_{\tau }\displaystyle \frac{\partial {v}_{\tau }}{\partial \tau }-\displaystyle \frac{{v}_{\sigma }({v}_{\sigma }\tau -{v}_{\tau }\sigma )}{{\tau }^{2}+{\sigma }^{2}}+\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial p}{\partial \tau }=0.\end{eqnarray} \tag{ 9 }$$

One of the momentum equations can be replaced by Bernoulli's equation

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}{v}_{a}^{2}=\displaystyle \frac{1}{2}{v}_{\sigma }^{2}+\displaystyle \frac{1}{2}{v}_{\tau }^{2}+\displaystyle \frac{1}{\gamma -1}{c}^{2}\end{eqnarray} \tag{ 10 }$$

where c is the speed of sound. We found it most convenient to express the equations in terms of the density ρ, the sound speed c, and the two Mach numbers ${m}_{\sigma ,\tau }=\tfrac{{v}_{\sigma ,\tau }}{c}$.

We assume that the shock front coincides with a curve on which $\sigma ={\sigma }_{s}$ is constant. The Rankine–Hugoniot boundary conditions at the shock front are (assuming strong shock/infinite Mach number)

$$\begin{eqnarray}&&{\rho }_{s}=\displaystyle \frac{\gamma +1}{\gamma -1}{\rho }_{a}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{c}_{s}=\displaystyle \frac{\sqrt{2\gamma (\gamma -1)}}{\gamma +1}{v}_{a}\displaystyle \frac{{\sigma }_{s}}{\sqrt{{\tau }^{2}+{\sigma }_{s}^{2}}}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{m}_{\sigma s}=-\sqrt{\displaystyle \frac{\gamma -1}{2\gamma }}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{m}_{\tau s}=\displaystyle \frac{\gamma +1}{\sqrt{2\gamma (\gamma -1)}}\displaystyle \frac{\tau }{{\sigma }_{s}}\end{eqnarray} \tag{ 14 }$$

where subscript a indicates a property of the ambient medium (upstream) and subscript s the value at the shock front (immediate downstream). Although an infinite Mach number was assumed, the correction for a finite Mach number scales as ${M}^{-2}$ (Landau & Lifshitz 1987, p. 467); we contend that the boundary conditionscan be applied to moderate values of Mach number without great loss of precision. At very high values of $\tau \gg {\sigma }_{s}$, the equations can be simplified. We assume that the variables vary with τ as the shocked values do, i.e., ρ is independent of τ, $c\propto \tfrac{1}{\tau }$, ${m}_{\sigma }$ is independent of τ, and ${m}_{\tau }\propto \tau$. Using this approximation, τ can be eliminated from the hydrodynamic equations, and the problem reduces to a set of ordinary differential equation in σ. These equations can be integrated numerically, and we can obtain curves for ρ, $\tau \cdot c$, m_σ, and $\tfrac{{m}_{\tau }}{\tau }$. The asymptotic equations for the interior of the shock are

$$\begin{eqnarray}&&2{m}_{\sigma }+2\left({m}_{\tau }\displaystyle \frac{\sigma }{\tau }\right)+{m}_{\sigma }\displaystyle \frac{d\mathrm{ln}\rho }{d\mathrm{ln}\sigma }+{m}_{\sigma }\displaystyle \frac{d}{d\mathrm{ln}\sigma }\mathrm{ln}\left(c\displaystyle \frac{\tau }{\sigma }\right)=0\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}2{m}_{\sigma } & - & 2\left({m}_{\tau }\displaystyle \frac{\sigma }{\tau }\right)-(\gamma -1){m}_{\sigma }\displaystyle \frac{d\mathrm{ln}\rho }{d\mathrm{ln}\sigma }\\ & + & 2{m}_{\sigma }\displaystyle \frac{d}{d\mathrm{ln}\sigma }\mathrm{ln}\left(c\displaystyle \frac{\tau }{\sigma }\right)=0\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}-{m}_{\sigma }^{2}-\displaystyle \frac{2}{\gamma } & + & {m}_{\sigma }\displaystyle \frac{{{dm}}_{\sigma }}{d\mathrm{ln}\sigma }+\displaystyle \frac{1}{\gamma }\displaystyle \frac{d\mathrm{ln}\rho }{d\mathrm{ln}\sigma }\\ & + & \left({m}_{\sigma }^{2}+\displaystyle \frac{2}{\gamma }\right)\displaystyle \frac{d}{d\mathrm{ln}\sigma }\mathrm{ln}\left(c\displaystyle \frac{\tau }{\sigma }\right)=0\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}-{m}_{\sigma }^{2}+\displaystyle \frac{2}{\gamma } & + & {m}_{\sigma }\displaystyle \frac{d}{d\mathrm{ln}\sigma }\left({m}_{\tau }\displaystyle \frac{\sigma }{\tau }\right)\\ & + & {m}_{\sigma }\left({m}_{\tau }\displaystyle \frac{\sigma }{\tau }\right)\displaystyle \frac{d}{d\mathrm{ln}\sigma }\left(c\displaystyle \frac{\tau }{\sigma }\right)=0.\end{eqnarray} \tag{ 18 }$$

These equations can be integrated numerically from the shock front to the axis of symmetry to produce the complete profile. An example of such a profile, for $\gamma =\tfrac{5}{3}$, can be seen in Figure 1, along with a comparison to a numerical simulation, which will be described in the next section.

Zoom In Zoom Out Reset image size

2.3. Growth of the Shocked Region

We describe the shock in a cylindrical coordinate system, where the $\hat{z}$ axis coincides with the axis of symmetry. In this section we consider a shocked region of finite extent, so there is some point where the r coordinate attains a maximum, and we denote that maximum by ϖ. The rate at which ambient matter accretes onto the shocked region is $\pi {\rho }_{a}{v}_{a}{\varpi }^{2}$. Matter can only be compressed by a factor of a few (for reasonable values of γ), so the volume of the bow shock $V\propto {\varpi }^{2}z\propto {\varpi }^{4}$ has to grow at the same rate:

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}({\varpi }^{4})\propto {\varpi }^{2}\Rightarrow \varpi \propto {t}^{1/2}.\end{eqnarray} \tag{ 19 }$$

Therefore, the farthest point from the obstacle on the axis of symmetry (on the lee side) moves at a constant velocity $z(\varpi )\propto {\varpi }^{2}\propto t$. In the absence of other scales, that velocity has to be of the same order of magnitude as the incident velocity. For shocks with a finite Mach number, $z(\varpi )\approx \min ({v}_{a}t,R({M}^{2}-1))$. Thus the age of a bow shock can be inferred from its size.

These results were verified with a numerical simulation (see next section). We examined the history of the mass enclosed within the bow shock. From the discussion above, it is clear that the mass scales with time as $M\propto V\propto {\varpi }^{4}\propto {t}^{2}$. Figure 2 shows that indeed the mass is parabolic in time.

Zoom In Zoom Out Reset image size

3.1. Solid Obstacle

To verify our results, we ran a simulation in the RICH code (Yalinewich et al. 2015). The simulations are in 2D cylindrical coordinates, assuming azimuthal symmetry. The boundaries of the computational domain are 180 > z > −20 and 25 > r > 0. We used a Cartesian grid with 100 cells in the r direction and 800 cells in the z direction. The velocity of the ambient medium is 1 in the positive z direction. The pressure was 10⁻⁹ and the density 1, which translates to a Mach number of more than 10⁴. The obstacle is a circle of radius 1, located at the origin $(0,0)$. The adiabatic index used was γ = 5/3. We ran the simulation to time t = 1000 and compared the numerical results with our analytic predictions. A density map from the last time step of the simulation can be seen in Figure 4. First, we fit the shock front to an even parabola (Figure 3), i.e., $y={{ax}}^{2}+b$. As can be seen from the figure, the shock front does converge to a parabola when the distance from the obstacle is much larger than the radius. The coefficients obtained from the simulation are a = 0.474 and b = −9.67. The ratio between the curvature radius of the shock front and the radius of the obstacle in this simulation is $\xi =\tfrac{R}{{R}_{o}}=1.06$. Next, we compared the numerical profiles inside the bow shock to the analytic prediction (Figure 1). In order to do that, we interpolated the variables on the curve τ = 1.5. The agreement between the two is reasonable. The deviations are due to two main reasons. The first is that since we are working in cylindrical coordinates, the relative error at each radius is equal to the ratio between the cell's width and the radius. The second is that the slope of the density is very steep next to the shock front, so a very fine resolution is required to resolve it. In the previous section it was shown analytically that very far from the obstacle the shock is parabolic, but the analytic theory did not constrain the shape of the shock front close to the obstacle (henceforth referred to as the "nose"). Since the only length scale in this problem is the size of the obstacle, the shape of the bow shock should scale as $\tfrac{z}{R}=f\left(\tfrac{r}{R}\right)$. The dimensionless function f can be obtained from the simulation. The ansatz chosen to represent f is a six-parameter Padé approximant in even powers of its argument:

$$\begin{eqnarray}&&f(\psi )=\displaystyle \frac{{c}_{0}+{c}_{1}{\psi }^{2}+{c}_{2}{\psi }^{4}+{c}_{3}{\psi }^{6}}{1+{c}_{4}{\psi }^{2}+{c}_{5}{\psi }^{4}}.\end{eqnarray} \tag{ 20 }$$

Odd powers were excluded to avoid a cusp at the z axis. The choice of powers in the numerator and denominator reproduces the analytic asymptotic behavior $z\propto {r}^{2}$ for $r\to \infty$. The calibrated coefficients are given in Table 1. A comparison between the numerical and analytic results is presented in Figure 3, and the deviations between the two are shown to be quite small.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Coefficients for the Ansatz (20), Calibrated from the Simulation

c0	−1.30
c1	−0.29
c2	−0.0548
c3	0.0215
c4	0.392
c5	0.0404

Download table as:  ASCIITypeset image

3.2. Wind Obstacle

The discussion so far has concentrated on a flow around a rigid obstacle. In astrophysics, however, it is extremely rare for a solid object to collide directly with diffuse matter. A more common scenario is that diffuse matter is halted at a distance much larger than the object itself, by magnetic field or by the wind emanating from the object. In this section we deal with the second case, which we refer to as wind-termination bow shock. In this case the effective radius of the obstacle (up to a dimensionless prefactor) is given by Equation (2). We claim that the formalism described above can also be used to describe bow shocks around wind obstacles. In order to put this idea to the test we ran another simulation, quite similar to the previous one, but instead of a solid obstacle we kept injecting matter into a small area around $(0,0)$. The velocity of the inner wind was ten times higher than the outer wind, and the mass loss rate was $\dot{M}=\tfrac{1}{10}$. At the other extreme, where the inner wind is much slower than the outer wind, the flow becomes unstable (Comeron & Kaper 1998; Decin et al. 2012). The shape of the forward shock front is given in Figure 5, and a density map from the last time step is given in Figure 6. The results of the simulation support our notion that the nature of the obstacle has very little effect on the shape of the bow shock around it.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this section we calculate the drag force on a rigid obstacle. In realistic cases, even dense clumps cannot withstand bow shocks indefinitely, and they eventually break up due to hydrodynamic instabilities. We can still justify this approximation by restricting it to the early times before the obstacle is disrupted.

From dimensional analysis one can infer an expression for the drag force,

$$\begin{eqnarray}&&{F}_{d}={c}_{d}\pi {R}_{0}^{2}{\rho }_{a}{v}_{a}^{2}\end{eqnarray} \tag{ 21 }$$

where c_d is a dimensionless drag coefficient. It is possible to obtain an analytic upper bound on c_d. The highest pressure on the obstacle is at the stagnation point. This is the point on the surface of the obstacle that is closest to the nose of the bow shock. The component of the velocity normal to the obstacle vanishes on the surface of the obstacle, and due to symmetry considerations, the tangential velocity also vanishes at that point. Therefore, using Bernoulli's equation and the adiabatic relations, the pressure can be evaluated there (McKee & Cowie 1975):

$$\begin{eqnarray}&&{P}_{o}=\displaystyle \frac{2}{\gamma +1}{\left(\displaystyle \frac{5{\gamma }^{2}-2\gamma +1}{4{\gamma }^{2}}\right)}^{\tfrac{\gamma }{\gamma -1}}{\rho }_{a}{v}_{a}^{2}.\end{eqnarray} \tag{ 22 }$$

So the upper limit on the drag force is $\pi {R}^{2}{P}_{o}$, and the upper limit on the drag coefficient is c_d = 0.827 (for $\gamma =\tfrac{5}{3}$). The drag coefficient can be evaluated numerically, using the pressure on the obstacle $P(\theta )$, where θ is the angle measured relative to the direction of the velocity of the ambient medium. If the pressure is known, then the drag force is given by

$$\begin{eqnarray}&&{F}_{d}=\pi {R}^{2}{\int }_{0}^{\pi }P(\theta )\sin (2\theta )d\theta .\end{eqnarray} \tag{ 23 }$$

For the simulation described above, the numerical drag coefficient is c_d = 0.426, which is approximately half of the upper limit discussed above. The assumption of a steady state can be justified only if the obstacle travels a distance much larger than its own radius before it experiences a significant change in velocity. This is equivalent to the condition that the density of the object be much larger than that of the ambient medium.

RX J1856.5-3754 is a neutron star that plows through the ISM at a speed of about 100 km s⁻¹, and is adorned by a bow shock (van Kerkwijk & Kulkarni 2001). The bow shock sports Balmer lines, indicating that the ISM is not ionized, hence the temperature is below 10⁴ K, so the Mach number is in excess of 10 (the magnetic field is bounded by equipartition (de Gouveia Dal Pino 2006) so the magnetic Mach number cannot be lower). Even if the Mach number is just 10, our solution for the shock should be valid within a distance of 100 times the size of the object. The distance between the neutron star and the nose of the bow shock is about 1'', and its tail extends to as far as 25''. Its distance is about 400 lt-yr and its inclination (angle between the line of sight and its velocity) is about 60°. The latter was inferred from measurements of proper motion and radial velocity. A comparison between our model and Hα observations (taken with the Very Large Telescope) can be seen in Figure 7. In order to perform the fit, we first manually digitized several points along the shock front.

van Kerkwijk & Kulkarni (2001) fit the observations to Wilkin's model (Wilkin 1996). The shock was shown to be adiabatic, whereas Wilkin's model assumes fast cooling. The figure below shows that our adiabatic model fits the observation. The analytic form used for our fit is the ansatz from (20). We note that in this fit we have only one free parameter, which is the size of the obstacle. We find that the best fit is obtained when this size is 09, which agrees with a previous estimate 1 ± 02 (van Kerkwijk & Kulkarni 2001). We note, however, that due to the uncertainties in the measurement, the fit to Wilkin's model is no worse than ours. In principle, it is also possible to use the fit between our model and the observations to estimate the inclination. Unfortunately, the data do not provide very stringent constraints, and only allow us to set a lower limit on the angle at 40°. From Figure 7 it is clear that inclination affects the wings of a bow shock more than it does the nose.

Zoom In Zoom Out Reset image size

We discussed the problem of a solid body traveling through a cold medium. We derived the steady-state, azimuthally symmetric hydrodynamics equations in cylindrical parabolic coordinates. We were able to use the approximation of very large τ to reduce these equations to a set of ordinary differential equations in σ. We numerically integrated these equations, and compared them to a full hydrodynamic simulation.

In the case of a large, but finite, Mach number, the bow shock would start out parabolic, but eventually straighten out and converge with the Mach cone. The transition will occur when the parabolic curve $z=\tfrac{1}{2}\tfrac{{r}^{2}}{R}$ intersects the Mach cone $z=\sqrt{{M}^{2}-1}r$:

$$\begin{eqnarray}&&{r}_{t}\approx {R}_{o}\sqrt{{M}^{2}-1}\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{z}_{t}\approx {R}_{o}({M}^{2}-1).\end{eqnarray} \tag{ 25 }$$

Therefore, the parabolic solution will only be valid in the region ${{RM}}^{2}\gg z\gg R$.

Using these analytic results, we modeled the bow shock around RX J1856.5-3754. The fit is as good as that of Wilkin's model. However, Wilkin's model assumes fast cooling, whereas the presence of Balmer lines around that bow shock indicates slow cooling (Heng 2010).

Further work is required in order to apply this model to most astrophysical bow shocks. For example, planetary bow shocks require magnetohydrodynamics, supernova–cloud interactions involve a time-varying flow, and the complete description of the bow shock around Betelgeuse requires radiative cooling (Mohamed et al. 2012).

However, this model may apply to "knots" in supernova remnants (Fesen et al. 2006). These are bright, filamentary structures observed inside supernova remnants, but also sometimes outside the shock front. It is theorized that a fraction of the ejecta is in the form of small, dense gas clumps, and the filaments are the shocks they leave in their wake. Since they are much more compact than the rest of the ejecta, they are less affected by deceleration, and can therefore overtake the shock front. While inside the supernova remnant, their velocity is of the same order of magnitude as the bulk and thermal velocity of the ejecta. However, once they emerge from the shock front, the bow shock that will develop around them will be the epitome of the solution presented here.

Another utility of our solution is that it provides a lower bound on the width of an adiabatic bow shock. Bow shocks with lower Mach numbers are bound to be wider than the solutions described here.

This research was partially supported by ISF and iCore grants. We would like to thank Marten van Kerkwijk for allowing us to use the image from his paper. We thank Elad Steinberg, Nir Mandelker, and Prof. Chris McKee for enlightening discussions.