3 Comparison with the analytic model

3.1 Shape

The geometry and internal structure of the bow shock nebula resulting from our simulations are shown in Figs. 1 and 2, while a grey-scale contour plot of the density is shown in Fig. 3. The positions of the termination shock, bow shock and contact discontinuity do not change among the above mentioned simulations, and also the position of the sonic point remains essentially the same. This is because the solution, in the hypersonic limit, depends only on the ratio between ram pressures of the stellar and ambient medium, but not separately on their densities or velocities: therefore we are confident that we have approximated well the hypersonic limit. The dashed lines represent three analytic solutions (Wilkin 1996), one scaled to match the contact discontinuity near the head, one to match the bow-shock near the head, and the long-dashed one, which is evaluated for the correct wind/environment ram pressure ratio. It can be seen that the analytic model reproduces quite well the shape of the nebula, although an uncertainty remains in the size. Figure 1 shows that the region of unperturbed stellar wind extends in the tail up to a "spherical shock'', whose position is strictly connected to the value of the thermal pressure in the ambient medium. In fact in the tail the pressure tends to the values of the ambient medium $P_{\rm o}$ so the distance of the spherical shock can be easily derived as:

$$d_{\rm back}\simeq\sqrt{{\cal L}/4\pi c P_{\rm o}}=d_{\rm o}M_{\rm o}\gamma_{1}^{1/2},$$ (4)

where $M_{\rm o}$ is the external flow Mach number.

We remark that the outflow boundaries we use may introduce spurious waves, and this plays a role in the structure of the subsonic region (region A) especially as it may change the position of the triple point (where the spherical shock meets the termination shock). However the head of the simulated nebula, and the supersonic region of the external shocked layer are not affected, so that the following comparison with the analytic model is still valid.

Looking at the head (Fig. 2) we notice a bump in the contact discontinuity, due to an instability that has not been damped completely by diffusion. The contact discontinuity in the tail (Figs. 1 and 3) tends to reach a limiting distance from the axis (i.e. the shocked pulsar wind is confined to a cylinder), with radius of the same order of the stagnation point distance; instead, the external bow-shock tends to reach the Mach cone.

An upper limit to the radius of the contact discontinuity can be evaluated using an asymptotic treatment. Asymptotically in the tail the bow shock degenerates in a Mach cone, so the pressure in the external layer is essentially the same as the thermal pressure of the ambient medium $P_{\rm o}$. Inside the contact discontinuity there are two different regions: the innermost (labeled as region A) corresponds to the material shocked in the spherical shock, which moves subsonically; the other region (region B), which surrounds the previous one, is filled with the material which crossed the termination shock at distances from the pulsar less than the spherical shock. As shown in Figs. 1 and 2, the material in the latter region moves supersonically. In a steady state model the pressure in these two regions must asymptotically match the conditions in the ambient medium. So the pressure inside the contact discontinuity is also $P_{\rm o}$ and this justifies Eq. (4) for the position of the spherical shock.

If we consider a fluid particle moving in region B, from a position near the stagnation point (labeled with subscript $_{\rm ini}$) towards the asymptotic region in the tail (labeled with subscript $_{\rm asy}$), the entropy conservation requires:

$$P_{\rm asy} /\rho_{\rm asy}^{\gamma_{2}}=P_{\rm ini}/\rho_{\rm ini}^{\gamma_{2}}.$$ (5)

On the other hand we have: $P_{\rm asy}=P_{\rm o}$, and:

$$P_{\rm ini}= A(\gamma_{1})\rho_{\rm o}V_{\rm o}^{2},\\ \rho_{\rm ini}=\rho_{\rm ups}B(\gamma_{2}),$$ (6)

where $\rho_{\rm ups}$ is the value of the density upstream of the wind termination shock along the axis, and:

$$A(\gamma) \left(\frac{\gamma+1}{2}\right)^{\frac{\gamma+1}{\gamma-1}}\left(\frac{1}{\gamma}\right)^{\frac{\gamma}{\gamma-1}} ,$$ =

$$B(\gamma) \frac{\gamma+1}{\gamma-1}\left(\frac{\gamma+1}{2}\right)^{\frac{2}{\gamma-1}}\left(\frac{1}{\gamma}\right)^{\frac{1}{\gamma-1}} \cdot$$ = (7)

If we consider an external wind with high Mach number, the position of the spherical shock is far from the star so we can assume that all the material of the stellar wind is confined to region B and that the distance Z of the contact discontinuity from the axis is much larger than the radius of internal region A. The matter flux conservation law allows us to write:

$$4\pi d_{\rm o}^{2}\rho_{\rm ups}V_{\rm wind} \simeq V_{\rm asy}\rho_{\rm asy}\pi Z^{2}.$$ (8)

From a two thin layer model we found that $V_{\rm asy}$, taken to be constant over the whole transverse section, is close to $V_{\rm wind}$, the stellar wind velocity. This can be easily understood because, if we assume that the tail extends to infinity (this constraint will be released below), the velocity scale for the stellar wind material, in the absence of mixing with the ambient medium, is $V_{\rm wind}$ (in the hypersonic limit the dependence on the external flow Mach number may be neglected). More specifically from the Bernoulli equation:

$$\frac{V_{\rm asy}^{2}}{2} + \frac{\gamma_{2}}{\gamma_{2}-1}\f... ...c{\gamma_{2}}{\gamma_{2}-1}\frac{P_{\rm ini}}{\rho_{\rm ini}},$$ (9)

together with Eq. (5), it follows:

$$V_{\rm asy}= V_{\rm wind}\sqrt{1-\left\{\gamma_{1}M_{\rm o}^{2}/A(\gamma_{1})\right\}^{\frac{1-\gamma_{2}}{\gamma_{2}}}}.$$ (10)

Substituting $\rho_{\rm asy}$ and $V_{\rm asy}$ in Eq. (8), using respectively Eqs. (5) and (10), we obtain:

 

$$\simeq 4d_{\rm o}^{2}\left\{\gamma_{1}M_{\rm o}^{2}A(\gamma_{1})\right\}^{\frac{1}{\gamma_{2}}}/$$ Z²

$$B(\gamma_{2})\sqrt{1-\left\{\gamma_{1}M_{\rm o}^{2}/A(\gamma_{1})\right\}^{\frac{1-\gamma_{2}}{\gamma_{2}}}}.$$ (11)

This can be considered as an upper limit to the radius of the contact discontinuity. For the stellar wind material to expand sideways beyond Z, some process is required that reduces the velocity in the layer or increases its pressure. This can be achieved via processes like mass accretion, or if some further shock develops.

The distance of the bow-shock from the star given by the analytic model is 0.175 Ul (long-dashed line) while the two limiting analytic solutions of Figs. 1 and 2 (dashed lines) are at 0.22 and 0.29 Ul. So the analytic model for the correct wind/environment ram pressure ratio gives a value which corresponds to the termination shock, while the thickness of the external layer is 0.025 Ul, in good agreement with what was expected (Paper I). Even though the shocked layers of the nebula are not thin at all, we see that the analytic solution reproduces quite well the morphology, at least of the external layer, which corresponds to the H$\alpha$  emitting region. In the tail simulation we see that the contact discontinuity has not reached the asymptotic regime, its distance from the axis being less than the maximum value given by Eq. (11).

The behaviour of the fluid in regions A and B may affect the synchrotron emission from the relativistic electrons and positrons coming from the pulsar and shocked in the termination shock. As we can see from Fig. 1, material flowing in region B becomes supersonic soon in the head. In a relativistic case, the Lorentz factor associated with the bulk motion of the fluid in the tail becomes greater than the Lorentz factor associated with the thermal motion of relativistic particles. If particles have relativistic motions their synchrotron emission can be detected only if their pitch angle is close to the observer direction. In the tail the relativistic bulk motion produces a beaming of the emission in the backward direction, so that the synchrotron emission from relativistic particles, in the case of a direction of observation perpendicular to the pulsar velocity direction, can be detected only in the head, while it fades in the tail. This holds true for the less energetic particles responsible for the radio emission and may even be true for those responsible for the X-ray emission. In region A the flow is subsonic, so the radio and X-ray emission can be detected. In this case we expect a cylindrical emitting region. If the spherical shock is far from the pulsar it may occur that the number of high energy particles confined in region A is so small that their emission might be under the threshold of detection, so that this region can be observed only in the radio band.

Figure 1: Tail simulation. TS is the termination shock, SS is the spherical shock, CD is the contact discontinuity (dot-dashed line), BS is the bow shock. The short-dashed lines represent two analytic solutions, the long-dashed line is the analytic solution for the correct wind/environment ram pressure ratio, the dotted lines indicates the sonic surfaces (the subsonic region is in the head). The circle corresponds to the region in which the stellar wind values are fixed every time step.

Figure 2: Head simulation. All curves and labels are defined as in Fig. 1.

Figure 3: Tail simulation. Density contours and shades in logarithmic scale are shown. The little spreading of the contact discontinuity is due to the diffusion we have introduced to avoid the growth of numerical instabilities.

3.2 Internal magnetic field

Using the numerical hydrodynamic solutions described above we have also evaluated the behaviour of a toroidal magnetic field wound around the symmetry axis and passively advected by the fluid. We considered the case of a pulsar wind with the energy associated with the ordered magnetic field much smaller than the ram pressure of the relativistic particles, and with the ratio between the two being independent of the direction. The value of such a field increases in the head, due to the compression in the termination shock, ranging from 7 times its upstream value, right after the shock, up to a factor of 30-35 with respect to its upstream value, near the contact discontinuity. In the tail it tends to decrease with respect to its value in the head due to the subsequent expansion of the fluid until it approaches a constant value:

$$B_{\rm asy} \simeq \frac{\pi}{2}B_{\rm ups}\left\{ \gamma_{1}M_{\... ...^{2}/A(\gamma_{1})\right)^{\frac{1-\gamma_{2}}{\gamma_{2}}}\right\}^{-1/4}\cdot$$ (12)

In fact the contact discontinuity tends to a limit section and all the magnetic field produced by the star is confined in such region. If the original magnetic field is not too far below equipartition, it can dominate the fluid dynamics both in the tail and in the head. So an hydrodynamical approach is asymptotically valid if:

$$\beta_{\rm asy}=\frac{8\pi P_{\rm o}}{B_{\rm asy}^{2}} \gg 1.$$ (13)

Following a procedure analogous to the one used to derive the asymptotic axial distance of the contact discontinuity we are able to relate the asymptotic condition to the initial values in the head:

$$\beta_{\rm asy} \simeq \left(\frac{2}{\pi}\right)^{2}\beta_{\rm i... ...}^{2}/A(\gamma_{1})\right)^{\frac{1-\gamma_{2}}{\gamma_{2}}}\right\}^{1/2}\cdot$$ (14)

So, even if initially the magnetic field's energy is below equipartition, it can reach equipartition if the external wind Mach number is high enough. Full MHD simulations are under development.

3.3 Surface density and tangential velocity in the external layer

To verify that the two thin layer model could be taken as a reasonable approximation of the true structure of the nebula, at least as far as the external layer is concerned, and to evaluate how far in the tail it could be used, we have extracted from our simulations the surface density and the tangential velocity in the external layer, between the bow shock and the contact discontinuity.

Figure 4 compares the surface density in the external layer, extracted from the tail simulation, with that evaluated by a two thin layer model. The numerical values have been obtained by integrating in the direction perpendicular to the bow shock. The difficulty in determining how good the analytic model approximation is follows from the fact that the analytic model depends on just one parameter, while the real problem has intrinsically two length scales (the stagnation point distance and the thickness in the head). In fact we have found two limit solutions (see Figs. 1 and 2), so that we could put only an upper and a lower limit to the length scale that should be used in the analytic model. The analytic curve that best reproduces the numerical points gives an external density of 0.26 cm^-3 (instead of the true value 0.25 cm^-3) and a stagnation point distance of 0.235 Ul, compatible with the range of the two limit solutions. The fluctuations of the numerical points are caused by round-off errors and uncertainties in the evaluation of the proper orthogonal direction, due to small oscillations in the shape of the bow shock, caused by discretisation errors in the evaluation of bow-shock position, eventually enhanced by the interpolation routine we use to find the orthogonal direction. In the simulation of the head we have noticed that, near the axis, the density values are lower than expected, and this effect corresponds to the bump of the contact discontinuity, that we have previously discussed. If corrected for the presence of the bump also in this zone we find reasonable values for the density.

Figure 4: Comparison of the surface density in the external layer, obtained from the tail simulation (crosses), with the values obtained from an analytic best-fit model (continuous line), and those from the analytic solution for the correct wind/environment ram pressure ratio (dashed line).

In the same way, we have evaluated the tangential velocity in the layer between the bow shock and the contact discontinuity. Figure 5 shows the numerical values of the tangential velocity in the external layer, normalized to the external wind value, and compares it to the two thin layer model. The analytic curve that reproduces the points has been obtained using a value of 0.25 Ul for the distance of the stagnation point, a value in agreement with the two limit solutions. The fluctuations in the numerical points are much smaller than those seen in the surface density, which implies that the velocity pattern is more uniform in the layer with respect to the density.

Figure 5: Comparison of the tangential velocity in the external layer, obtained from the tail simulation (crosses), with the velocity obtained from an analytic best-fit model (continuous line), and those from the analytic solution for the correct wind/environment ram pressure ratio (dashed line).

The fact that the analytic curves that best reproduce the numerical values correspond to different values of the stagnation point distance tells us that the density and velocity cannot be reproduced exactly by the same thin layer model. It should be stressed, though, that, in spite of the non-thin structure of the external layer, the thin layer model is not a bad approximation. Deviations from the continuous curves are present only in the tail where a steady state regime may have not been achieved yet.