Full-3D relativistic MHD simulations of Bow Shock Pulsar Wind Nebulae: dynamics

TL;DR

This study uses 3D relativistic MHD simulations to explore how pulsar wind properties influence the diverse shapes and behaviors of bow shock pulsar wind nebulae, focusing on fluid structure, magnetic fields, and turbulence.

Contribution

It provides a comprehensive analysis of how different pulsar wind injection conditions affect nebular dynamics and morphology in BSPWNe through detailed simulations.

Findings

Variations in wind anisotropy and magnetization impact nebular shape.

Turbulence levels correlate with injection parameters.

Magnetic field structures and current sheet survival depend on wind properties.

Abstract

Bow shock pulsar wind nebulae (BSPWNe) are know to show a large variety of shapes and morphologies, both when comparing different objects, and for the same object in different energy bands. It is unclear if such a variety is related to differences in the pulsar wind properties, or to differences in the conditions of the ambient medium. We present here a set of full three-dimensional, relativistic and magneto-hydrodynamic simulations of BSPWNe, with the intention of determining how differences in the injection conditions by the pulsar wind reflect in the nebular dynamics. To achieve a good coverage of the available parameter space we have run several simulations varying those parameters that are most representative of the wind properties: the latitudinal anisotropy of the wind energy flux with respect to the pulsar spin axis, the level of magnetization, the inclination of the pulsar spin axis with respect to the pulsar direction of motion. We have followed the dynamics in these systems, not just in the very head, but also in the tail, trying to assess if and how the system retains memory of the injection at large distances from the pulsar itself. In this paper we focus our attention on the characterization of the fluid structure and magnetic field properties. We have tried to evaluate the level of turbulence in the tail, and its relation to injection, the survival of current sheets, and the degree of mixing between the shocked ambient medium and the relativistic pulsar wind material.

Tables

Table 1. Table 1: List of the runs. Subscript i 𝑖 i in the run name indicates the value of the initial magnetization (in the exponent notation specified in the Table note).

Run	Dimensions	$\alpha$	${\varphi }_M$	${\sigma }^{\left[a\right]}$	$t_f\left[t_0\right]$
2D-H0	2D	0	0	0	1000
H0	3D	0	0	0	740
I{0,i}	3D	0	0	${10}^{-i}$	1152
I{π/4,i}	3D	0	$\pi /4$	${10}^{-i}$	1152
I{π/2,i}	3D	10	$\pi /2$	${10}^{-i}$	1152
A{0,i}	3D	10	0	${10}^{-i}$	1152
A{π/4,i}	3D	10	$\pi /4$	${10}^{-i}$	1152
A{π/2,i}	3D	10	$\pi /2$	${10}^{-i}$	1152
[a] For MHD runs values of the magnetization are given by the exponent $i=[2,1,0]$, such that $\sigma =[{10}^{-2},{\mathrm{\hspace{0.17em}10}}^{-1},{\mathrm{\hspace{0.17em}10}}^0]$.

Equations

$$d_0=\sqrt{\dot{E}/(4\pi c{\rho }_{\mathrm{ISM}}{v}_{\mathrm{PSR}}^2)},$$

$$R_w(\theta )=d_0\csc \theta \sqrt{3(1-\theta \cot \theta )},$$

$$\rho =\{\begin{array}{ll}({\rho }_0/r^2)\mathcal{F}(\psi )& \mathrm{if}r\le 2d_0\mathrm{and}r\le 2R_w(\theta );\\ {\rho }_0/d_0^2& \mathrm{if}r\ge 2d_0\mathrm{and}r\le R_w(\theta );\\ {\rho }_{\mathrm{ISM}}& \mathrm{otherwise},\end{array}$$

$$\mathcal{F}(\psi )=N_0\left(1+\alpha {\sin }^2\psi \right),$$

$$p=\{\begin{array}{ll}{p}_0{\left(\rho /{\rho }_0\right)}^{\Gamma }& \mathrm{if}r\le 2d_0\mathrm{and}r\le 2R_w(\theta );\\ {\rho }_0/[6(d_0{\gamma }_0{)}^2]& \mathrm{if}r\ge 2d_0\mathrm{and}r\le R_w(\theta );\\ p_{\mathrm{ISM}}& \mathrm{otherwise},\end{array}$$

$$\mathbf{B}=\left(\frac{B_0\sin \psi }{\mathcal{R}}\right)\left(\frac{d_0}{r_0}\right)\times \{\begin{array}{l}\left(-y_0\cos {\varphi }_M+z_0\sin {\varphi }_M\right){\mathbf{e}}_{\mathbf{x}};\\ \left(x_0\cos {\varphi }_M\right){\mathbf{e}}_{\mathbf{y}};\\ -\left(x_0\sin {\varphi }_M\right){\mathbf{e}}_{\mathbf{z}};\end{array}$$

$$S_2(\mathbf{\ell })=⟨|v_{x,y}(\mathbf{r}+\mathbf{\ell })-v_{x,y}(\mathbf{\ell }){|}^2⟩,$$

Full text

1 Introduction

Pulsar Wind Nebulae (PWNe) are synchrotron emitting bubbles powered by a pulsar wind. They are commonly found inside the remnant of their parent supernova explosions, but in older systems they can also arise due to the direct interaction with the ISM. The pulsar wind is a magnetized, cold and ultra-relativistic outflow, with predicted Lorentz factors in the range $10^{4}-10^{7}$ (Goldreich & Julian, 1969; Kennel & Coroniti, 1984a, b)), and it is thought to be mainly composed by electron-positron pairs (Ruderman & Sutherland, 1975; Arons & Scharlemann, 1979; Contopoulos et al., 1999; Spitkovsky, 2006; Tchekhovskoy et al., 2016; Hibschman & Arons, 2001; Takata et al., 2010; Timokhin & Arons, 2013; Takata et al., 2016). It is launched at the pulsar magnetosphere (at the typical distances of the light cylinder) at the expenses of the stellar rotational energy. As a consequence of the interaction with the ambient medium, the supersonic wind is forced to slow down, and it does so by a strong termination shock (TS), where particles are likely accelerated to a non-thermal distribution. The observed radiation is produced via non-thermal processes (synchrotron and inverse Compton scattering) arising from the interplay of these particles with the nebular magnetic field and the background photon field.

Exemplary PWNe involving young pulsars inside their supernova remnant (SNR) are the Crab and Vela nebulae. Systems like these are referred as plerions (Gaensler & Slane, 2006; Bucciantini, 2008; Olmi et al., 2016). However, given that a relevant fraction of all the pulsars, between 10% and 50%, is born with high kick-velocity, of the order of 100-500 km s*-1*, (Cordes & Chernoff, 1998; Arzoumanian et al., 2002; Sartore et al., 2010; Verbunt et al., 2017)), while on the other hand the remnant expansion is decelerated (Truelove & McKee, 1999; Cioffi et al., 1988; Leahy et al., 2014; Sánchez-Cruces et al., 2018), they are fated to escape their progenitor SNR over timescales of a few tens of thousands of years, much shorter than the pulsars typical ages of the order of $10^{6}$ years. Once outside, given the typical sound speeds in the ISM of the order of $10-100$ km s*-1*, these pulsars turn out to move at a supersonic speed. No longer observed as plerion-like bubbles, their associated nebulae acquire a cometary like shape due to the balance of the ram pressure of the pulsar wind with the surrounding ISM (Wilkin, 1996; Bucciantini & Bandiera, 2001; Bucciantini, 2002). The pulsar is now found at the head of an elongated tail that extends in the direction opposite to the pulsar motion. These objects are then known as bow shock PWNe (BSPWNe).

This scenario was confirmed by numerical simulations in different regimes (Bucciantini, 2002; Bucciantini et al., 2005a; Vigelius et al., 2007; Barkov et al., 2019) and it is illustrated in Fig. 1.

BSPWNe have been observed in detail at different wavelengths in the last years, revealing a puzzling variety of structures at the different scales: different morphologies in the bow shock head, different shapes and elongation of the tails, with astonishing outflows misaligned with the pulsar velocity. They have been mainly revealed as non-thermal radio and X-ray emission (Arzoumanian et al., 2004; Gaensler et al., 2004; Yusef-Zadeh & Gaensler, 2005; Li et al., 2005; Gaensler, 2005; Chatterjee et al., 2005; Kargaltsev et al., 2008; Misanovic et al., 2008; Ng et al., 2009; Hales et al., 2009; Ng et al., 2010; De Luca et al., 2011; Ng et al., 2012; Marelli et al., 2013; Jakobsen et al., 2014; Auchettl et al., 2015; Klingler et al., 2016; Posselt et al., 2017), with polarimetric informations only available for a few cases, suggesting a large variety of magnetic configurations (Ng et al., 2012; Yusef-Zadeh & Gaensler, 2005; Ng et al., 2010; Klingler et al., 2016; Kargaltsev et al., 2017). Extended TeV halos have also been recently detected surrounding some BSPWNe (Abeysekara, 2017), attracting much attention because these nebulae are thought to be one of the major contributors of leptonic anti-matter in the Galaxy, in competition with possible dark matter sources (Blasi & Amato, 2011; Amato & Blasi, 2017). If the pulsar is moving through a partially ionized medium the collisional and/or charge exchange excitations of neutral hydrogen atoms in the tail can also allow one to reveal the BSPWN as Hα emission (Chevalier et al., 1980; Kulkarni & Hester, 1988; Cordes et al., 1993; Bell et al., 1995; van Kerkwijk & Kulkarni, 2001; Jones et al., 2002; Brownsberger & Romani, 2014; Romani et al., 2017), or alternatively in the UV (Rangelov et al., 2016), and IR (Wang et al., 2013).

Neutral contamination can also modify the tail dynamics and morphology, as theoretical and numerical models have recently predicted (Morlino et al., 2015; Olmi et al., 2018). The nature of this changes can be very important as an indicator of the typical speeds in the tail and of the nature of the turbulence that, together with the magnetic configuration, are determinant in understanding how particles escape from these objects (Bucciantini, 2018), fundamental to properly assess the efficiency of pair contamination in the ISM.

The earliest attempts to model BSPWNe date back to more than one decade. Apart from analytical and semi-analytical works (Wilkin, 1996; Bandiera, 1993), the first numerical models have been presented in the classical hydrodynamical regime by Bucciantini (2002) and van der Swaluw et al. (2003). Due to their intrinsic limitations, first of all the absence of the magnetic field, these models can only account for the properties of the outer layer of shocked ISM, whence the Hα emission comes and where the magnetic field is not expected to be very important, and for the global structure of the system. The first effort to include magnetic field was done by Bucciantini et al. (2005b). Here results from a series of 2D MHD relativistic and axisymmetric numerical simulations have been presented, with particular attention to the effects of different magnetization (the ratio of the magnetic to the kinetic energy fluxes) on the PWN structure, and its emission signatures. The authors found that the level of magnetization does not influence the external layer of the shocked ISM, confirming the previous findings of HD models. The growth of magnetization causes an enhanced emission near to the symmetry axis and modulations of emissivity are also observed in association with shear instabilities in the nebula.

The axisymmetric assumption of those simulations however forced the authors to consider only the case of a pulsar spin axis aligned with the pulsar kick velocity (and also to include ad hoc magnetic dissipation in the bow shock head to avoid numerical artifacts related to the enforced symmetry of the system). However there is no reason to expect such alignment (Johnston et al., 2005; Ng & Romani, 2007; Johnston et al., 2007; Noutsos et al., 2012, 2013). To account for the spin-kick inclination one needs to work in full 3D. The third dimension is particularly important for the correct description of the magnetic field configuration (Bucciantini, 2017), and it is crucial in the study of the development of turbulence, which can be strongly affected by geometric constraints. Limits of 2D models of PWNe in that sense have been largely discussed in the last years (Del Zanna et al., 2006; Olmi et al., 2013, 2014; Porth et al., 2013), since the first 3D MHD models of PWNe become available (Porth et al., 2014; Olmi et al., 2016). They confirm that the full 3D modeling allows the development of a complex structure of the magnetic field which is artificially prevented in 2D, due to geometrical limitations. Moreover, the growth of 3D instabilities leads to an higher magnetic dissipation, such that higher values of the injected magnetization can be reached (Porth et al., 2013). 2D models were proved to be still robust in reproducing the properties of the inner region of PWNe (Del Zanna et al., 2006; Volpi et al., 2008; Camus et al., 2009; Olmi et al., 2014), since deviations from these simplified models become mainly important in the outer region of the PWN. This might however not be satisfactory in the case of BSPWNe, where the tail regions extend far away from the pulsar location.

The first attempt to 3D modeling of BSPWNe has been done by Vigelius et al. (2007) but limited to the classical HD regime. They present a series of simulations of bow shocks not very extended in the tail region, considering different latitude variations in the pulsar wind and various density gradients in the ISM. In particular they found that the wind anisotropy does not directly influence the bow shock morphology in the head, which is on the contrary more deeply affected by the interaction with the anisotropic ISM. More recently 3D MHD simulations have been presented in Barkov et al. (2019). Different inclinations of the magnetic field with respect to the pulsar spin axis, and various magnetizations are considered, together with a possibly non uniform external density of the ISM. The authors conclude that the magnetic field has a fundamental influence in accounting for the bow shock morphology, while anisotropies in the ISM density only lead to marginal variations of the external surface (Barkov et al., 2019; Toropina et al., 2018). However some of the choices in these models (the initial wind Lorentz factor is less than 3, TS is not fully detached from the simulation boundary) raise a few questions of the robustness of the results.

Aim of the present work is to investigate in the full 3D relativistic MHD regime the morphology of BSPWNe, varying the pulsar parameters: inclination of the magnetic axis with respect to the pulsar spin axis, anisotropy of the wind momentum flux and level of magnetization. We are in particular interested in determining how the large scale structure of the bow shock and its tail are influenced by the morphology of the magnetic field and by the level of turbulence. In the present paper we will then focus our attention on the analysis of the dynamics. A detailed analysis of the emission, polarimetry and variability properties would be postponed to successive papers.

This paper is organized as follows: in Sec. 2 the physical model of the pulsar wind, the numerical tool and setup used for our simulations are described; in Sec. 3 we present and discuss our findings. Conclusions are finally drawn in Sec. 4.

2 Pulsar wind model and numerical setup

Our simulations have been performed with the numerical code PLUTO (Mignone & McKinney, 2007). PLUTO is a shock-capturing, finite-volume code for the solution of system of hyperbolic and parabolic partial differential equations, particularly oriented to astrophysical fluid dynamics. It supports the use of adaptive mesh refinement (AMR) thanks to the CHOMBO libraries, a fundamental requirement for the present study (Mignone et al., 2013). The numerical grid must in fact account for a sufficient resolution to capture simultaneously the pulsar wind injection region and the large scale of the cometary tail of the pulsar bow shock nebula.

A second order Runge-Kutta time integrator and an HLLD Riemann solver (the Harten-Lax-van Leer for discontinuities, see Miyoshi & Kusano 2005; Mignone & Bodo 2006) have been used, in order to better treat shear layers in the nebula. The HLLD solver is automatically relaxed to the HLL one in the regions of the domain occupied by the relativistic pulsar wind, in order to be able to use a Lorentz factor $\gamma=10$, high enough to ensure that the TS jump conditions are in the relativistic regime, even in those region of the TS that are strongly inclined with respect to the wind itself. This has been made possible due to a customization of the code. The solenoid condition of the magnetic field, $\mathbf{\nabla}\cdot\mathbf{B}=0$ is maintained using a divergence-cleaning method (Dedner et al., 2002).

The typical length scale of BSPWNe is the so called stand off distance $d_{0}$ (Wilkin, 1996; Bucciantini & Bandiera, 2001; van der Swaluw et al., 2003; Bucciantini et al., 2005a), where the wind momentum flux and the ISM ram pressure balance each other

[TABLE]

where $\dot{E}$ is the pulsar spin down luminosity, $\rho_{\mathrm{ISM}}$ is the ISM density, $v_{\mathrm{PSR}}$ the speed of the pulsar with respect to the local medium, and $c$ the speed of light. The typical values of the stand-off distance are of the order of $\sim 10^{16}-10^{17}$ cm (Chevalier et al., 1980; Kulkarni & Hester, 1988; Cordes et al., 1993; Chatterjee & Cordes, 2002; Romani et al., 2010; van Kerkwijk & Kulkarni, 2001; Brownsberger & Romani, 2014; Romani et al., 2017). Similarly a convenient choice for parametrizing time in the inner tail is the characteristic flow time $t_{0}=d_{0}/c$, while the relaxation time with respect to the interaction with the ISM is typically $t=d_{0}/v_{\mathrm{PSR}}$.

The simulation domain in cartesian coordinates is centered on the pulsar with the $z$-axis aligned with the pulsar kick velocity (the reference frame is moving with the pulsar, thus the ISM is seen as a uniform, unmagnetized flow moving along the negative $z$ direction). The domain extends in the range $[-17d_{0},\,+17d_{0}]$ along the $x$ and $y$ direction and $[-28d_{0},\,5d_{0}]$ along $z$. This choice have been made in order to allow for the development of a long tail behind the pulsar and to reduce as much as possible the domain occupied by the un-shocked ISM. The base grid have $128^{3}$ equally spaced grid points. A set of four AMR levels is then used in order to reach the required resolution around the pulsar, in order to resolve at best the bow shock head, which is a strongly dynamic region, corresponding to an effective resolution of $2048^{3}$ cells at the highest level. The configuration of AMR levels can be seen in Fig. 2.

In order to speed up the simulation towards a steady relaxed regime, at the beginning the domain is divided into two regions with separatrix given by the analytical bow shock shape (Wilkin, 1996)

[TABLE]

where $(r,\,\theta,\,\phi)$ are the spherical coordinates, centered on the pulsar, and with the polar axis aligned to the $z$ direction. In $r>R_{w}$ the domain is initialized with a cold slowly moving ISM, with a speed $v=-v_{PSR}=0.1c$. In $r\leq R_{w}$ and $r>2d_{0}$ we set a uniform density and pressure moving backward with a speed $v=-v_{PSR}=0.1c$ to allow the bow shock tail to develop without the eventual formation of back-flows towards the pulsar. The value of the PSR velocity in the ISM has been set in order to speed up the relaxation, thus imposed to be $v_{PSR}=0.1c$, in the $-z$ direction. This value was shown to be small enough to recover the proper dynamics of the non-relativistic ambient medium (Bucciantini, 2002; Bucciantini et al., 2005b). In $r\leq R_{w}$ and $r\leq 2d_{o}$ we initialize the relativistic pulsar wind, with radial velocity $v_{r}=(1-1/\gamma_{0}^{2})^{1/2}/c$. For simplicity, given that we are here mostly interested in the dynamics of the relativistic component, the ISM is assumed to be unmagnetized. The external magnetic field will be important for the escape of particles from the bow shock into the ambient medium, while being energetically subdominant in the ISM, its effect on morphology of the bow shock nebula will be negligible. This setup is shown in Fig. 1, where a 2D slice of the density profile illustrates the initial conditions. The pulsar wind is continuously injected from a radius $r_{\mathrm{inj}}\simeq 0.2d_{0}$. Such value for the injection boundary ensures that the TS radius is always detached during the entire evolution. The possibly anisotropic injection of energy from the pulsar wind is accounted for with a modulation of the wind density

[TABLE]

with $\rho_{0}=[v_{\mathrm{PSR}}d_{0}/(c\gamma_{0})]^{2}\rho_{\mathrm{ism}}$. The wind Lorentz factor is kept constant at the value $\gamma_{0}=10$, high enough to ensure that the post-shock dynamics is independent on its value (Del Zanna et al., 2004). Modulation is governed by the function $\mathcal{F}(\psi)$, connected to the wind energy flux by the expression

[TABLE]

where $\psi$ is the colatitude of the generic point $(x_{0},\,y_{0},\,z_{0})$ from the pulsar spin axis, defined from the relation $\cos{\psi}=(y_{0}\sin{\phi_{M}}+z_{0}\cos{\phi_{M}})/r_{0}$, where $r_{0}^{2}=(x_{0}^{2}+y_{0}^{2}+z_{0}^{2})$ and $\phi_{M}$ is is the spin-kick inclination (configurations with $\phi_{M}=0$ have the spin axis aligned with the bow shock symmetry axis $z$). $N_{0}$ is a normalization factor, fixed by the requirement $\int_{\Omega}\mathcal{F}(\psi)d\Omega=4\pi$ and $\alpha$ a dimensionless parameter governing the level of anisotropy of the wind, with the isotropic distribution recovered for $\alpha=0$.

The pressure in the wind is simply given by entropy conservation with adiabatic index $\Gamma$, while is set fixed the other regions

[TABLE]

with $p_{0}=0.01\rho_{0}(c/d_{0})^{2}$ and $\Gamma=4/3$, the appropriate value for the relativistic shocked pulsar plasma. In Bucciantini (2002) more sophisticated simulations with a multi-fluid treatment were realized: different equations of state are considered for describing the relativistic component of the wind and the non-relativistic material of the ISM. The authors found that finally this only leads to minor deviations to the overall geometry of the bow shock. A comprehensive treatment of the two components can be also achieved by the use of adaptive equations of state, like Taub’s one (Taub, 1948; Mignone et al., 2005). However they should not be used in the presence of strong mass contamination, since they are based on the assumption of instantaneous thermalization between the two components. This can easily happen in turbulent flow due to numerical diffusion, but physically this is expected to happen on timescales much longer than the typical flow time in the tail of a BSPWN, given the mass and energy difference of the two components.

In the wind the magnetic field in a generic point $(x_{0},y_{0},z_{0})$, is defined in cartesian coordinates as

[TABLE]

where $\mathcal{R}=(x_{0}^{2}+y_{0}^{2}\cos^{2}{\phi_{M}}+z_{0}^{2}\sin^{2}{\phi_{M}}-y_{0}z_{0}\sin{2\phi_{M}})^{1/2}$. The previous formula is applied in the region $r\leq 2r_{b}$, while the magnetic field is imposed to be 0 in the rest of the domain.

The runs are summarized in table 1. We investigated various magnetizations, with $\sigma=[0.01,1,0.1,\,1.0]$, and inclination of the pulsar spin axis $\phi_{M}=[{0,\,\pi/4,\,\pi/2}]$, both for isotropic and anisotropic distribution of energy flux in the wind. In order to speed up simulations we evolve a pure hydrodynamic isotropic case up to $t_{H,f}=740t_{0}$, sufficient for the configuration to be completely relaxed. This is then used as the starting point for all the other runs. Anisotropy and wind magnetization are gradually increased to avoid spurious jumps in the injection conditions. Once the injection conditions have reached their desired final value, the system is evolved again up to a quasi-stationary configuration until the time $t_{f}=1152t_{0}$, enough for the dense slow ISM material to stream from the head of the bow shock to the far tail. We have verified in each case that at $t_{f}$ the tail is relaxed to a final state.

3 Results and Discussion

In this section we analyze the dynamics resulting from the different configurations listed in Table 1. Notice that all the quantities are always given in code units, that can be rescaled to physical units in terms of $d_{0}$, defined by eq. 1. As a useful quantity for converting between the two we chose the external ram pressure, which in code units for us is $p_{\mathrm{ext}}=\rho_{\mathrm{ISM}}v_{\mathrm{PSR}}^{2}=0.01$. The typical ram pressure for bow shocks in cgs units is $p_{\mathrm{ext}}=10^{-10}\rho_{24}v_{7}^{2}$ g cm*-1* s*-2*, where density is expressed in units of $1\times 10^{-24}$ g cm*-3* and velocity in units of $10^{7}$ cm s*-1*.

We begin by comparing the large scale structure of bow shock in 2D versus 3D.

An overview image comparing different runs from 2D to 3D and from HD to MHD can be found in Fig. 3, where density profiles are shown as 2D slices in the $(x,\,z)$ plane for runs, from the left to the right-side, 2D-H0, H0, I*{0, 0}, A{0, 0}* and A*{π/4, 0}*. Here we can notice a few major differences from case to case: HD runs appear to be more affected by small scales turbulence, while MHD cases (all with $\sigma=1.0$), show a global structure much less influenced by turbulence. As we will see in details in the following, this is a characteristic of high-magnetized cases, in which the properties of the fluid are usually dominated by conditions at injection. The other clear difference is the shape and dimension of the TS, which as expected is very different from the isotropic to the anisotropic model, where it assumes an oblate shape. The outer shape of the bow shock instead appears to be quite similar from case to case. The only exception is the right-most plot, referring the case with $\pi/4$ spin-kick inclination, which shows an evident variation of the FS shape. As we will see in the following it is also the only one characterized by an evident asymmetry around the $z-$axis.

In Fig. 4, the velocity magnitude in the 2D-H0 case, at $t=t_{f}$, is directly compared with 2D slices from runs I*{0,2}* and I*{0,0}* in the $x=0$ plane. These 3D runs have a spin axis aligned with the kick axis, and thus the injection conditions are fully axisymmetric. The shape and size of the TS are essentially the same: the location of the TS in the front region, as well as the distance of the Mach disk on the backward side, show no significative difference. 3D simulations show however a TS having a drop-like shape more than a bullet-shape as in 2D, partly due to the collimating effect of the magnetic field (a toroidal field wounded around the $z-$axis), and partly due to the development of turbulence in the tail. It is interesting to observe that, while in 2D there are several inclined weak shocks in the tail, those tend to disappear in 3D, where a more turbulent flow develops. The behavior of turbulence however is strongly dependent on the magnetization of the wind. High values of $\sigma$ lead to flow structures in the tail that are more laminar, and this has two consequences. On one hand the turbulent ram pressure lead to a sideway expansion of the tail, which is more pronounced in the I*{0,2}* run, while the I*{0,0}* has a tail cross section comparable to the 2D HD result. On the other the development of turbulence slows down the bulk flow in the tail. Typical values of the flow speed in the tail ranges from $0.6c$ in the I*{0,2}* to $0.9c$ in I*{0,0}* run. One also sees that in all cases, the inner slow channel that was found in early 2D simulations (Bucciantini et al., 2004) does not extend more than $5d_{0}$ in the tail, and is completely lost at distance larger than $10d_{0}$.

We have also compared the shape of the forward bow shock among all our runs. This can be seen in Fig. 5, where the position of the FS is given as a contour line of the density for all the configurations with $\sigma=1.0$. Contours are superimposed on a density map from the 2D hydrodynamic case, for direct comparison. Here the forward shock shape and position appear to be very similar from case to case, with the major deviations coming from case A*{π/4,0}*, mostly because the inclination leads to the development of waves that propagate backward modulating the shock shape. A sort of lateral blobs are visible at the FS surface, possibly representing periodic structures. In any case their extension is in general smaller than $d_{0}/2$, or $\sim 1d_{0}$ at most, and should probably be completely invisible to actual instruments, due to resolutions limits. This confirms that the outer shape of a BSPWN appears then to be quite independent on the properties of the pulsar wind, even at high magnetization, but rather dominated by the interaction with the ambient medium, as predicted by previous models (Romani et al., 1997; Vigelius et al., 2007).

We also found that the distance between the bow shock and the contact discontinuity (which is usually named $\Delta$) is in perfect agreement with the predicted value of $5/16d_{0}$ (Chen et al., 1996) for the isotropic model, while in the anisotropic case it is more variable due to the different morphology of the shock, and it goes from 5/16 up to 5/8, depending on the magnetic inclination.

The effect of an anisotropic energy distribution in the pulsar wind is shown in Fig. 6 where the left-half side of each map shows the $v_{z}$ component of the velocity in a cut of the 3D dataset, while in the right-half side of the map the Lorentz factor is drawn, with the purple area indicating the shock surface (at value $\gamma=10$). As expected in the anisotropic model, where the energy flux in the wind is not uniform, the TS assumes an oblate shape, with major extension in the direction in which the energy flux peaks (i.e. the direction orthogonal to the magnetic axis). In this sense the shock tends to reflect more the structure of the energy injection than the motion of the PSR. Here the elongation of the TS into a bullet shape is lost. The deformation of the shock is largely responsible for the formation of turbulence in the anisotropic runs: the bow shock head appears to be more dynamic, with vortexes on smaller scales and less coherent structures. The region in which the high velocity channel streams along the oblate part of the shock appears to be thinner, and randomization of the flow causes this structures to dissolve near to the head region. Moreover the different shape of the shock also causes an important change in the zone occupied by the slow inner channel, which now turns into a “stagnation zone” of the flow behind the TS.

The formation of this “stagnation zone” is caused by the material that is first diverted backward to the shock front in the $+z$ direction, then, due to the mixing of the flow in the region behind the shock, pushed on the back of the shock and diverted again back to the tail, in the $-z$ direction. This can be appreciated in Fig. 7 as the effect of the material flowing along the shock front and then squeezed towards the axis around $z=-6$, giving rise to the formation of a turbulent region behind the shock. The map shows a 2D slice of the velocity magnitude and streamlines from run A0,0, which shows a much extended “stagnation zone” than the isotropic case, visible as a light-blue bubble behind the shock in Fig. 6. This region is more extended in the anisotropic case due to the different shape of the shock front, which causes more mixing of the flow behind the shock and in the tail. This also reflects in a more turbulent motion in the tail, causing the velocity to lower and the tail to broaden. This effect is not seen in 2D simulations, more likely because to the enforced axisymmetry of those models, which prevent the accumulation of material along the axes.

To characterize the differences in the flow structure between the 2D and 3D cases, we have looked at the behavior of various fluid and magnetic quantities along the tail, performing transversal cuts. These have been extracted as 2D slices in the $(x,\,y)$ plane of the 3D data cube, at three distinct positions along the bow shock tail: $z=[-20.0,\,-20.8,\,-21.7]d_{0}$. This choice has been done for two reasons: at these locations the tail has reached a quasi-stationary configuration, and they are far enough from the bow shock head that one can neglect to first order the dynamics of the very head; on the other hand the region very close to the boundary of the numerical grid has been excluded to avoid numerical artifacts. We sample every dynamic variable along different transverse directions, that can be though of as diameters of an arbitrary circumference of fixed radius ($\sim 13d_{0}$), centered on the symmetry axis ($x=y=0$), as illustrated in Fig. 8, where they are drawn as colored lines.

Profiles are given as functions of the position along these directions, and the directions indicated by arrows in Fig. 8. In this representation the symmetry axis is located at $\mathcal{D}=0d_{0}$ and the plane containing the pulsar spin axis is the one selected by the vertical line identified by label $\mathcal{L}$. The profiles obtained along each of these directions are then averaged. We first average at fixed $z$, spanning all the different diameters, and then average over the three different slices (thus varying $z$). The variability is then encoded by the standard deviation, which is also shown (in grey for 2D and different colors for 3D). By the way we perform our averages, the standard deviation has a minimum at the center of the circumference. Profiles obtained with this procedure are then compared with similar profiles obtained from the 2D simulation 2D-H0. In this case we have used six cuts in the same range $z=[-20.0;-21.7]d_{0}$ at three consecutive time-steps, to derive more robust averages. Results are shown in Fig. 10 and Fig. 10.

We found that the mean properties of the tails show very similar trends in all cases, thus we only show few representative runs, focusing on the role of the magnetic field: I*{0,i}* and A*{0,i}*, both with $\phi_{M}=0$. As we will discuss in the following, deviations from these trends are seen only for runs with $\phi_{M}=\pi/4$.

From a first look at Fig. 10 and Fig. 10 no major asymmetry is observed, consistent with the axisymmetric injection of these runs. As discussed before the major difference is the lateral extension of the tail marked by the position of the contact discontinuity. In the 2D hydrodynamic case its diameter is of the order of $\sim 10d_{0}$ on average. In the 3D cases instead the extent of the tail is larger for lower values of the magnetization ($\sim 13d_{0}$), and only in the case $\sigma=1.0$ we recover a size comparable with the 2D-H0 case. What we see here is that the density jump at the CD is shallower al lower magnetization. The extension and variability of the shear layer can be estimated from Figs. 10-10, looking to the $v_{z}$ plots. It can be seen that the shear layer become less extended while increasing magnetization for both the isotropic and anisotropic models, going from $\sim 6.5d_{0}$ for $\sigma=0.01$ to $\sim 3d_{0}$ for $\sigma=1.0$, with a maximum variability of $\sim 0.9d_{0}$. These number must be compared with the 2D case, where the shear layer extension is of order $\sim 4d_{0}$, with a small variability of $\sim 0.2d_{0}$. As already pointed out before, as the magnetization increases the 2D prediction is recovered.

The properties of the shear layer can be equivalently extrapolated from density plots. In the 2D case there is a sharp jump between the PWN material and the shocked ISM. On the contrary this contrast is much less evident in the 3D cases. Here the higher value of the density in the center is a clear signature of the mixing with the denser ambient medium, which is more efficient in 3D than in 2D. The only exception is seen for case I*{0, 1}*, which was already indicated as the most similar to the 2D one. In that case the density within the PWN is compatible with the one found in 2D. In any case the inner density approaches the 2D one with increasing magnetization, meaning again that the importance of the turbulence diminishes with increasing the magnetization level.

This confirms what was seen before in term of turbulence developing at the shear layer between the PSR material and the shocked ISM. Interestingly for models with an anisotropic energy injection the tail appears to be more extended than in the isotropic ones. This can be again explained in terms of turbulence: anisotropic energy injection tends to create stronger shear layers in the very head, and a more turbulent flow is injected from the head of the system.

The third columns of Fig. 10 and Fig. 10 show the thermal pressure inside the bow shock nebulae. A common feature for both the isotropic and anisotropic cases is the decrease of the average pressure with increasing magnetization (as the consequence of the rising of the magnetic pressure). For $\sigma=0.01$ pressure is almost everywhere higher than the 2D case, while for $\sigma=1.0$ the 3D trend is more or less the same than the 2D one. In all cases pressure appears to be highly variable in the tail, with the maximum reached at the FS surface. It is also interesting to notice that in 2D the pressure shows a rapid drop after the CD, at $D=\pm 7d_{0}$, which is not shown by any of the 3D runs, where pressure is maintained quite uniform up to the boundary or it even presents an increase at the same location. Moreover pressure plots do not show evident differences between the isotropic and anisotropic models.

The difference between the isotropic and anisotropic cases is again visible when comparing the velocity profiles. Going back to Fig. 10-10, the first columns displays the $z$ component of the velocity, aligned with the pulsar motion. In the isotropic models the velocity shows a slower channel around the symmetry axis similar to what observed in 2D, even if it is not as marked, where velocity is, in absolute value, of order of $\sim 0.65c$ within a distance of $\sim 0.8d_{0}$, to be compared with the fast outer channel with $\lvert v_{z}\rvert\sim 0.82c$. In the 3D isotropic cases the maximum velocity variation from the inner to the outer channels is seen for case I*{0, 1}*, where velocity goes from $v_{z}\sim 0.7c$ to $v_{z}\sim 0.82c$ and the extension of the slow velocity channel is of order of $d_{0}$ in radius. This behavior is however absent for models with anisotropic energy injection where, on the contrary, velocity is usually maximum along the $z-$axis. Again this is an effect of different level of turbulence that can destroy the coherence of flow structures. Moreover in the isotropic models $v_{z}$ remains higher if compared with the anisotropic runs. In the isotropic cases the average speed within a region of extension $\sim 5d_{0}$ around the symmetry axis is $\lvert v_{z}\rvert\gtrsim 0.75c$, and the flow appears to be more collimated around the $z$ axis. On the contrary in the anisotropic case the velocity is slightly lower (with $0.6c\lesssim\lvert v_{z}\rvert\lesssim 0.8c$) and the high velocity flow is distributed on greater distances from the symmetry axis, following the broadening of the tail.

The different properties of the tail are also influenced by the level of magnetization that survives in the tail and by the structure of the local magnetic field, which play an important role on the development of the turbulence. In Fig. 11 the magnitude and morphology of the magnetic field are shown for some reference runs at $t=t_{f}$. The upper panel compares the field intensity and morphology for the minimum ($\sigma=0.01$, leftmost part of each map) and maximum ($\sigma=1.0$, rightmost part of each map) values of initial magnetization for the aligned isotropic run (i.e. I${}_{\{0,\,i\}}$, left panel) and the anisotropic one (i.e. A${}_{\{0,\,i\}}$, panel on the right). The bottom panel of the same figure shows the whole map of the magnetic field for the same magnetizations for both the isotropic and anisotropic models and magnetic inclinations $\phi_{M}=\pi/4$, $\phi_{M}=\pi/2$. Here it is evident how the survival of highly magnetized regions far from the bow shock head helps in maintaining coherent structures on large spatial scales and that injection with a higher magnetization (left-half part of each map) prevents the formation of small scale turbulence.

The comparison between different values of the magnetization shown in the upper panel of Fig. 11 for runs I*{0,i}* and A*{0,i}*, illustrates that in the case of low magnetization the magnetic field is almost fully dissipated and randomized already at distances of the order of $z\sim 15d_{0}$ from the pulsar. On the contrary for $\sigma=1.0$ the magnetic field survives far further from the BSPWN head: the random motion in the tail is lower and coherent structures are clearly visible on large scales and the tail is maintained more collimated around the symmetry axis. The global maximum of the magnetic field is not so different from case to case and from model to model, with the maximum usually arising near to the front shock at the BSPWN head, where the injected magnetization is still higher (due to the minor dissipation it has suffered) and the field is enhanced by shear instability and compression. Values of the maximum ranges between $2\times 10^{-2}$ and $\sim 10^{-1}$ in code units, corresponding to $\sim 10-50\,(\rho_{24}v_{7}^{2})^{1/2}\,\mu$G. The global average value of the field is of the order of $10^{-3}$ in code units or $0.5\,(\rho_{24}v_{7}^{2})^{1/2}\,\mu$G, as can be seen form the extracted average value in the tail.

Due to turbulence and dissipation the average value of the magnetization is always less than the injected value, with values of a few$\times 10^{-4}$ for $\sigma=0.01-0.1$ up to a few$\times 10^{-3}$ for $\sigma=1.0$. The same behavior is shown from runs with different inclinations of the magnetic field (bottom panel of the same figure).

Compared to other variables as density and velocity, the magnetic field appears also to have the higher level of variability, as can be seen looking at the standard deviations in Fig. 10-10. Interestingly the magnetic field seems to reach the maximum near to the contact discontinuity, possibly reflecting the presence of an effective amplification mediated by the shear instability, even if in the low magnetized regime, early 2D simulations were suggesting also the possibility of local compression into a magnetopause (Bucciantini et al., 2005b). On the other hand for higher magnetization we see a local enhancement of the magnetic field strength close to the axis. This is reminiscent of results in 2D RMHD, and it is due to a residual effect of hoop stresses associated with the toroidal component of the magnetic field, that tends to pinch the flow on the axis. However the effect is reduced in 3D, where there are no symmetry constraints.

A closer look to the magnetic field variations can help in investigating the nature of the turbulence. In Fig. 12, 1D profiles for the $B_{z}$ component of the magnetic field are given, for all the magnetizations, considering cases I*{π/4,i}* and A*{π/4,i}*, for which the effect of turbulence appear stronger. Please beware that the scales of the plots, given in units of $10^{-2}B-\mathrm{code}$, are not the same for the different plots, since these have been chosen in order to highlight the field oscillations. Differences in the amplitude must be then read taking into account the different scales. First thing to notice is that the magnetic field has an intrinsic turbulent nature, clearly visible from the fact that the average value of $B_{z}$ is compatibile with 0, with very high oscillations.

In the case of the anisotropic wind (bottom row) for $\sigma\leq 0.1$ the increase of the injected magnetization does not reflect in a growth of the magnetic field: maximum values and amplitude of the variability do not change, with the maximum approaching 0.002 in code units, i.e. $1\,(\rho_{24}v_{7}^{2})^{1/2}\,\mu$G. This means that in this regimen the magnetic field evolution and dynamics is actually dominated by turbulence rather than injection.

On the contrary, as $\sigma$ exceeds 0.1, the magnetic field is seen to rise with increasing initial magnetization, even if not linearly (juts a factor $2$ for a rise of magnetization by a factor 10), meaning that the system now retains memory of the injection.

In the isotropic wind model (upper row), each step in magnetization corresponds to an increment of a factor of 2 in the maximum value of $B_{z}$, which goes approximately from $0.002$ to $0.1$ in code units, i.e. from $1(\rho_{24}v_{7}^{2})^{1/2}\,\mu$G up to $50(\rho_{24}v_{7}^{2})^{1/2}\,\mu$G. In this case injection is evidently dominant and the development of small scale turbulence is less efficient.

Planar components of the magnetic field can be seen in Fig. 13, where $B_{x}$ and $B_{z}$ are plotted as 2D $(x,y)$ slices (at $z=-20.8d_{0}$) for runs I*{π/4,0}* and A*{π/4,0}*. The nature of turbulence appears very different between the isotropic and anisotropic wind models , with the anisotropic case dominated by small scale structures. Again these differences causes clear diversities in the field morphology: in the isotropic case the field shows coherent structures on large scales, and the injection geometry is still recognizable in the tail far away from the head, as can be seen by the different polarities of the field still well separated (yellow color vs violet). On the contrary in the anisotropic model the high level of mixing causes this global structure to be very poorly recognizable in the tail. The $B_{z}$ component maintains a very faint left-right separation of polarities, while the mixing in the $B_{x}$ component (and $B_{y}$) is almost complete, and no coherent structure is visible any more. This behavior is common for almost all the considered configurations. The exception, as already pointed out in the previous discussion, is represented by the low magnetized cases ($\sigma=0.01$), where the level of turbulence is high enough to dominate on the injection properties and to destroy the large scale structure of the field also in the isotropic case. In the isotropic case the planar component of the magnetic field peaks in a region extending between $1-4d_{0}$, and it shows a shell-like morphology (with thickness of order of $\sim 1.5d_{0}$), with subsequent well separated shells of different polarities. Here the maximum value is of the order of $\pm 5(\rho_{24}v_{7}^{2})^{1/2}\,\mu$G. On the contrary the $z$ component peaks near to the CD in a shell of variable thickness ($\sim 1.3d_{0}$ at the bottom up to $\sim 2.5d_{0}$). The magnetic field is particularly intensified in the bottom region, $y<0$, showing a magnification of $\sim 3.5$, while it shows a minor enhancement (a factor of $\sim 1.5$) for $y>0$. This enhancement is probably the effect of shear instability acting at the CD surface. On the other hand in the anisotropic case both the planar and $z$ components of the magnetic field peak near to the CD, in different regions of the $(x,\,y)$ plane. Moreover the maximum value is more or less the same, still indicating that the magnetic field is mostly dominated by turbulence in all the directions.

A close-up view of the bow shock head is given in Fig. 14. Here the head morphology is compared for cases I*{0,2}* and I*{0,0}* as 2D maps of the module of the current density, given in code units. This zoom-in helps to see that the pulsar wind region is perfectly resolved in our simulations, even for the highly magnetized case (where the shock is expected to be smaller). The injection region, i.e. a sphere of radius $\sim 0.2d_{0}$ centered on the pulsar, appears in fact well far away the TS surface in each direction. From this figure the evolution of the current sheet can be followed, comparing low and high-$\sigma$ cases. The current sheet originates at the pulsar equatorial plane and appears to be evidently different in the two cases. In the low$-\sigma$ case it maintains its structure up to a distance of $\sim 5d_{0}$ from the pulsar. At that point vortexes start to mix the current sheet with the surroundings, causing it to twist and tangles while it is completely mixed up. At a distance of $\sim 12d_{0}$ from the pulsar, the current sheet is no more recognizable from the background. On the other hand, when magnetization is higher, the current sheet survives and it is easily visible up to $\sim 14d_{0}$ from its birth region. It maintains a quasi-laminar structure up to a distance $\sim 10d_{0}$ from the pulsar, then it starts to be diverted.

Again these differences points out that in presence of an higher magnetized flow turbulence is less efficient in mixing up the plasma and quasi-laminar structures can survive far away from the bow shock head. The high$-\sigma$ case also shows a maximum value of the current greater than a factor of $\sim 3$ with respect to the low magnetized case. When magnetization is higher the current also appears to be enhanced at the CD surface. Notice that the high-magnetized case is also the one in which the magnetic field feels the major increment as the consequence of the rise in $\sigma$.

As anticipated major deviations from the common trends shown by the average quantities in the tail arise in cases with $\phi_{M}=\pi/4$. This can be observed in Fig. 15, where profiles for pressure and $v_{z}$ are shown. They present a clear left-right asymmetry, which gives rise to a distorted shape of the BSPWN with respect to the other inclinations. The global structure of the bow shock is visible in Fig. 16, where isocontours of the velocity magnitude are superimposed to a 3D map of density for run A*{0,1}*, to highlight the asymmetry of this configuration with respect to the $z-$axis. As predicted the TS has an oblate shape, with major extension in the direction inclined by $\pi/4$ with respect to the $z-$axis. Velocity is only shown for values $\gtrsim 0.7c$ and only a few spots of $v\sim 0.8c$ are visibile in the tail, concentrated at approximately the half of the tail extension and in the left-side part of the image, with Lorentz factor $\sim 1.6$. This asymmetry, which also reflects in a global asymmetry of the entire bow shock, is a direct effect of the shock shape and inclination. In this configuration the flow coming from the TS has very different evolution depending on the part of the shock at which it is originated: when the elongated part of the TS points to the outer shock (in the right-hand part of the figure, pointing in the $+z$ direction), the flow is directed towards the outer shock, where it is diverted and decelerated. On the contrary, when the elongated part of the shock points to the back of the tail, in the $-z$ direction, the high velocity flow coming from the shock is not decelerated by the interaction with the outer shock and it may survive along the tail. This is in fact what is seen in Fig. 15, where almost all the high velocity flow ($v>0.7c$) is located behind the part of the shock pointing in the $-z$ direction rather than in the opposite. This obviously leads to an asymmetric shape of the entire object and produces a preferential region with high velocity in the direction in which the alignment of the shock front with the pulsar motion is maximized. For obvious geometrical reasons the same effect is not seen when $\phi_{M}=\pi/2$ or $\phi_{M}=0$.

A closer look to the local structure of the flow in the bow shock tails for all the magnetic configurations for $\phi_{M}=\pi/4$ can be seen in Fig. 17 . Here 2D $(y,z)$ maps of the magnetization are shown to compare scales of the turbulence in the isotropic (left column) and anisotropic (column on the right) models.

Maps are normalized to their own maxima and only shown in a region of the tail identified by $z\in[-15,-7]d_{0}$ and $y\in[-4,4]d_{0}$. Even at a local view the different dynamic of the shocked wind between the isotropic and anisotropic models is evident, especially when increasing the initial magnetization. In the anisotropic model small scales turbulence is dominant at all the magnetizations, with the smallest eddies showing a dimension of $\sim 0.3d_{0}$ and being dominant in the flow structure for the lower magnetized cases. On the contrary in the isotropic model small scales turbulence is only poorly dominant and disappears when the initial magnetization is increased. As already discusses the case with $\sigma=1.0$ shows a very ordered flow dynamics, with a quasi-laminar appearance, with elongated structures well recognizable along the pulsar direction of motion.

Among all our runs, case A*{π/4,2}* is the one showing the more developed turbulence in the velocity field. In order to characterize this turbulence we have computed the second order structure function of the velocity. Traditionally the structure function is computed for the longitudinal and trasverse components (with respect to the displacement) of the velocity. Here we instead compute the second order structure function of velocity components $v_{x}$ and $v_{y}$ perpendicular to the axis of the tail, in planes taken at different locations in $z$. We do this for several reasons: first $v_{z}$ in the tail shows a clear pattern that reflects the net bulk motion of the flow, shaped by the tail itself, and defining a turbulent part is non trivial; second, given that the pulsar spin axis belongs to the $z-y$ plane, any memory of the injection condition could translate into a difference in the level of turbulence along the $x$ and $y$ directions; third, by taking planes at different $z$ we can sample how the turbulence varies moving along the tail. Our structure function is defined as:

[TABLE]

where $\bm{\ell}$ is the displacement vector in the $(x,y)$ plane. In general, for fully developed incompressible and compressible MHD turbulence, the structure function is shown to scale as $\propto\ell^{\zeta_{p}}$ ,(She & Leveque, 1994; Grauer et al., 1994; Politano et al., 1995; Müller & Biskamp, 2000; Boldyrev et al., 2002; Padoan et al., 2004), with $\zeta_{p}=p/9+2-2(2/3)^{p/3}$, such that $S_{2}(\ell)\propto\ell^{0.7}$.

In order to increase out statistics we consider different slices in $z$ in the same range used for Fig. 17 (i.e. for $z\in[-15,\,-7]d_{0}$) and square sectors in the orthogonal plane with extension [-4, 4]$d_{0}$ in the $x$ and $y$ directions. The second order structure function is shown in Fig. 18.

Note that moving along the tail away from the pulsar, the values of $S_{2}(\ell)$ drops. This is probably the effect of the sideway expansion of the tail in the orthogonal plane, which obviously increases with distance from the pulsar. We found that our results do not match the analytical prediction: $\zeta_{2}\simeq 0.7$. This should not be surprising at large scales $\gtrsim d_{0}$, that are compatible with the injection scale given by the tail size. On the other hand, at small scales, we found a trend compatible with the value $\zeta_{2}=1.025$ found in relativistic MHD simulations of turbulence in low-magnetized plasmas (Zrake & MacFadyen, 2012). Scales smaller than $0.1d_{o}$ cannot be investigated properly due to the limitations of the numerical grid. We also have verified that these results are not dependent on the resolution or on the AMR resampling.

4 Conclusions

In this paper we present a complete description of the dynamics of pulsar bow shock nebulae arising from full 3D relativistic MHD numerical simulations. In order to discuss how the pulsar parameters influence the tail dynamics we have considered a large set of different configurations, listed in Table 1. In particular we consider three different initial magnetization, from no-magnetization to highly magnetized, with injected values $\sigma=[0.0,\,0.01,\,0.1,\,1.0]$ and different configurations for the magnetic field geometry, with the angle between the magnetic axis and the pulsar velocity inclined by $0,\,\pi/4$ and $\pi/2$. Since inhomogeneities in the ISM density are expected to influence primarily the morphology of the forward shock but not very much the dynamics within the contact discontinuity, we only consider here a uniform ISM density, with standard properties. The pulsar wind is modeled considering both a uniform distribution of the energy flux in the wind (the isotropic model) and a non-uniform distribution, with the maximum of the energy flux in the pulsar equatorial plane (the anisotropic model). This helps in isolating the influence of the wind properties on the overall morphology of the bow shock.

We first compare our findings with results from a 2D simulation with the same model in the hydrodynamic regime. Our comparison begin with the average behavior of the physical quantities in the tail in the plane orthogonal to the velocity direction. Comparison of 3D runs with 2D HD shows a good agreement in the case of the isotropic wind model, with a quasi perfect overlapping for the intermediate value of the magnetization ($\sigma=0.1$). This could appear unexpected, since in principle the major similarity may be expected with the lower magnetized case. But what we found is that for values of the magnetization $\sigma<0.1$ the dynamics is completely dominated by turbulence on small scales. For greater values on the contrary, injection become dominant. At the light of these findings is then not surprising that the 2D hydrodynamic case is more compatible with the intermediate case. The difference in the development of the local turbulence and in the level of mixing of the fluid when moving from 2D to 3D can be better appreciated in Fig. 3, comparing the first two density maps from the leftmost side, that shows exactly the same model from 2D and 3D simulations.

We also found that in the direction aligned with the pulsar motion ($z-$direction) the velocity shows the same structure obtained in 2D for the isotropic wind model, with a lower velocity channel around the $z$ axis ($\sim 0.65c$ in 2D and $\sim 0.7-0.8c$ in 3D) surrounded by an higher velocity flow (with $\sim 0.85c$ in 2D and up to $0.9c$ in 3D). The same similarity is not seen when comparing 2D HD behavior with the anisotropic wind model. In that case velocity is even maximum along the symmetry axis, with the only exception of run A*{0,0}*. Moreover the peak velocity is always lower or equal to the maximum value of the 2D HD case. In the anisotropic cases we also found that the effect of turbulence is even more pronounced: injection becomes important only in for $\sigma>0.1$, even if structures in the flow remain coherent on scales that are much shorter than the same for the corresponding isotropic case. The high level of mixing of the fluid clearly shows the effect of the presence of turbulence.

The same considerations remain valid for the magnetic field. The initial configuration of the magnetic field, and the structured current sheet separating the field polarities, survive in the isotropic cases, while in the anisotropic ones turbulent mixing on small scales tends to destroy the initial configuration, which is only poorly recognizable in the tail, even in the case of higher magnetization. We found that turbulence does not amplify the magnetic field efficiently: when dynamics is dominated by turbulence rather than injection, magnetic field appears to not increase raising the initial value of the magnetization. On the contrary, as soon as turbulence becomes no more dominant (as it happens for $\sigma\gtrsim 0.01$ in the isotropic model and for $\sigma\gtrsim 0.1$ in the anisotropic one), magnetic field starts to increase with the initial magnetization. In 3D magnetic field also appears to be enhanced near to the contact discontinuity, as the effect of an efficient shear instability amplification acting at the CD.

From the morphological point of view the forward shock structure is almost identical from case to case, with the only exception of the anisotropic cases with $\phi_{M}=\pi/4$, which show the larger deviations from the typical shape. Differences in fact appear as small extrusions and blobs, possibly resulting as periodic perturbation of the FS, that will be investigated with major details in following works. Moreover they do not appear to be characterized by high values of the magnetic field and density or pressure, thus it is not clear if they should be visible as emission. Notice also that they in any case arise on very small spatial scales, possibly making it difficult to be revealed by actual instruments due to resolution limits. The observed similarity of the FS between different models should not be surprising, given the forward shock morphology is expected to be primarily influenced by the interaction with the ambient medium and its properties (Romani et al., 1997; Vigelius et al., 2007), taken as fixed in our models. We also confirm that the variation of the magnetization does not reflects in a significant modification of the forward shock, as was previously shown with 2D MHD models (Bucciantini et al., 2005b).

Comparing the bow shocks global morphologies, major differences from case to case are visible mainly in terms of collimation or broadness of the tails, arising as the effect of different magnetic inclination and especially the wind anisotropy. These variations of model also influence the global structure and dynamics of the magnetic field, that may be in principle lead to observed differences in the high-energy emission, especially at X-rays, where emission is dominated by the most relativistic particles near to their injection site. This can be less important for radio emission, where the density distribution of emitting particles should be dominant on the magnetic field configuration (Olmi et al., 2014). Properties of the emission will be addressed with accuracy in a successive papers, to which we postpone this discussion.

We found that the dynamics in 3D is fully dominated by small scale turbulence for values of the magnetization lower than $\sim 0.1$.This is particularly evident for the anisotropic wind model, where turbulence is more efficient in destroying the injection properties. A developed turbulence is well recognizable when looking at the local properties of the flow. For the low magnetized anisotropic case with $\phi_{M}=\pi/4$, which shows the most turbulent dynamics in the tail, we made a more quantitative characterization of the turbulence. We compute velocity structure functions in a sector of the bow shock tail, considering different slices in $z$. We found the velocity structure functions of the two planar components of the velocity, $v_{x}$ and $v_{y}$, have a trend which is compatible with results in the literature for MHD relativistic, low-magnetized, turbulence. For scales comparable with the injection scale ($\sim 1d_{0}$) we indeed found a flatter trend, that is not surprising since we do expect injection, due likely to shear at the contact discontinuity, to involve different scales along the tail.

Overall our results are in line with to those of Barkov et al. (2019). However, due to important differences in the choice of the wind energy and magnetic field distribution, the dynamic of some fluid variables presents some substantial differences. This is not unexpected given that it is well known from standard MHD modeling of PWNe (see for example Del Zanna et al. 2006 vs Komissarov 2004) that these choices affect significantly the post-shock dynamics. For example in our model there is no polar current, which is instead strongly enhanced in the work by Barkov et al. (2019), and this impacts not just the distribution of currents in the nebula but, through its dissipation, also the structure of pressure. This can be seen from a direct comparison of our Fig. 14 vs their figure 13. Moreover we find that for $\phi_{M}=0$ the thickness of the PWN shocked layer in the head, between the TS and the CD, is much smaller than in their equivalent case. However it can be seen that in their simulations the shock is not properly detached from the boundary of the injection region.

Our results also show that the average flow pattern in the tail is not dissimilar between 2D and 3D runs (major differences are mostly in the head due to inclination and anisotropy). However in the low $\sigma$ regime we found a strongly turbulent magnetic field structure, suggesting that perhaps a laminar model is not likely to fully capture the magnetic field, not even at the level of its average strength. For higher values of $\sigma$ instead, a more coherent magnetic field is found, qualitatively in agreement with the prediction of simplified laminar models (Bucciantini, 2017).

We postpone the discussion of the implications of our findings on the observed properties to future works. We in fact intend to present a detailed study of the emission at multi wavelengths, together with the analysis of the time-variability and polarimetric properties of bow shock nebulae from our numerical models, comparing these with available observations. Moreover, it remains to be investigated the escape of high-energy particles from bow shocks. Since pulsars are known to be one of the most efficient antimatter factories in the Galaxy, the study of the trajectories and escape of particles is of great interest. For this problem one must take into account the external magnetic field and additional configurations need to be investigated depending on the different possible inclinations of the ISM magnetic field with respect to the pulsar velocity. This will be particularly relevant in the far tail of the bow shock, where we found a magnetic field comparable with typical ISM values.

Acknowledgements

We acknowledge the “Accordo Quadro INAF-CINECA (2017)” for the availability of high performance computing resources and support. Simulations have been performed as part of the class-A project “Three-dimensional relativistic simulations of bow shock nebulae” (PI B. Olmi). The authors also acknowledge financial support from the PRIN-MIUR project prot. 2015L5EE2Y ”Multi-scale simulations of high-energy astrophysical plasmas”. B. Olmi wishes also to acknowledge Andrea Mignone, from the PLUTO team, for fundamental support, Simone Landi and Luca Del Zanna for fruitful discussions.