3D Relativistic MHD numerical simulations of X-shaped radio sources

TL;DR

This paper uses three-dimensional relativistic MHD simulations to explore how low-power jets in elliptical galaxies can naturally form X-shaped radio sources, revealing the influence of jet power and orientation on morphology.

Contribution

It extends previous two-dimensional studies by performing 3D simulations, demonstrating the conditions under which X-shaped radio sources form, including jet power and misalignment effects.

Findings

Jets with power $ extless 10^{44}$ erg s$^{-1}$ produce X-shaped morphologies.

A misalignment of $30^ ext{o}$ still favors X-shaped formation.

Synthetic maps suggest different observed shapes result from viewing angles.

Abstract

A significant fraction of extended radio sources presents a peculiar X-shaped radio morphology: in addition to the classical double lobed structure, radio emission is also observed along a second axis of simmetry in the form of diffuse wings or tails. In a previous investigation we showed the existence of a connection between the radio morphology and the properties of the host galaxies. Motivated by this connection we performed two-dimensional numerical simulations showing that X-shaped radio sources may naturally form as a jet propagates along the major axis a highly elliptical density distribution, because of the fast expansion of the cocoon along the minor axis of the distribution. We intend to extend our analysis by performing three-dimensional numerical simulations and investigating the role of different parameters is determining the formation of the X-shaped morphology. The problem is addressed by numerical means, carrying out three-dimensional relativistic magnetohydrodynamic simulations of bidirectional jets propagating in a triaxial density distribution. We show that only jets with power $\lesssim 10^{44}$ erg s$^{-1}$ can give origin to an X-shaped morphology and that a misalignment of $30^o$ between the jet axis and the major axis of the density distribution is still favourable to the formation of this kind of morphology. In addition we compute synthetic radio emission maps and polarization maps. In our scenario for the formation of X-shaped radio sources only low power FRII can give origin to such kind of morphology. Our synthetic emission maps show that the different observed morphologies of X-shaped sources can be the result of similar structures viewed under different perspectives.

Tables

Table 1. Table 1: Parameters used in the simulation models. In the first column we have the case identifier, in the second column the angle between the jet axis and the major axis of the triaxial density distribution, in the third column the jet Lorentz factor, in the fourth column the density ratio, in the fifth column the magnetization parameter σ 𝜎 \sigma , in the sixth column the jet power, in the seventh column the size of the computational domain in units of r c subscript 𝑟 𝑐 r_{c} and in the eighth column the grid size. L j subscript 𝐿 𝑗 L_{j} is measured in erg s -1 and is computed for r c = 4 subscript 𝑟 𝑐 4 r_{c}=4 kpc and ρ 0 = 1 subscript 𝜌 0 1 \rho_{0}=1 cm -3 , the dependence of L j subscript 𝐿 𝑗 L_{j} on r c subscript 𝑟 𝑐 r_{c} and ρ 0 subscript 𝜌 0 \rho_{0} is given by Eq. 9

Case	$\theta$	${\gamma }_j$	$\nu$	$\sigma$	$L_j$	Domain	Grid
A	0	$5$	$1\times {10}^{-4}$	${10}^{-3}$	$1.6\times {10}^{46}$	$2\times 4\times 2$	$384\times 768\times 384$
B	0	$2.5$	$4\times {10}^{-4}$	${10}^{-3}$	$1.6\times {10}^{46}$	$2\times 4\times 2$	$384\times 768\times 384$
C	0	$1.25$	$2.66\times {10}^{-3}$	${10}^{-3}$	$1.6\times {10}^{46}$	$4\times 4\times 4$	$768\times 768\times 768$
D	0	$5$	$1\times {10}^{-5}$	${10}^{-2}$	$1.6\times {10}^{45}$	$4\times 8\times 4$	$768\times 1024\times 768$
E	0	$5$	$1\times {10}^{-6}$	${10}^{-2}$	$1.6\times {10}^{44}$	$4\times 8\times 4$	$768\times 1024\times 768$
F	${30}^o$	$5$	$1\times {10}^{-6}$	${10}^{-2}$	$1.6\times {10}^{44}$	$4\times 8\times 4$	$768\times 1024\times 768$

Equations

$${\partial }_{\mu }\left(\rho u_{\mu }\right)=0,$$

$${\partial }_{\mu }\left(w{u}^{\mu }{u}^{\nu }-b^{\mu }{b}^{\nu }+p{\eta }^{\mu \nu }\right)=0,$$

$${\partial }_{\mu }\left(u^{\mu }{b}^{\nu }-u^{\nu }{b}^{\mu }\right)=0,$$

$$\rho =\frac{{\rho }_0}{[1+(x^{\prime }/a_x{)}^2+(y^{\prime }/a_y{)}^2+(z/a_z{)}^2{]}^{3/2\beta }}$$

$$x^{\prime }=x\cos (\theta )-y\sin (\theta )$$

$$y^{\prime }=x\sin (\theta )+y\cos (\theta )$$

$$B_x=-B_{t}r\sin (\varphi )$$

$$B_y=B_{t}r\cos \varphi$$

$$L_j=\pi r_j^2{\rho }_j{\gamma }_j^2{c}^3=2.5\times 10^{-3}\pi {\rho }_0{r}_c^2\nu {\gamma }_j^2{c}^3$$

$$v_h=\frac{{\gamma }_j\sqrt{\eta }}{1+{\gamma }_j\sqrt{\eta }}{v}_j,$$

$$\eta =\frac{{\rho }_j{h}_j}{{\rho }_a{h}_a}$$

$$v_w=\sqrt{\frac{p_c}{{\rho }_a}}$$

$$j^{\prime }({\nu }^{\prime })=K|{\mathbf{B}}^{\prime }\times {\mathbf{n}}^{\prime }{|}^{\alpha +1}{\nu }^{\prime -\alpha }$$

$$j(\nu )=K|{\mathbf{B}}^{\prime }\times {\mathbf{n}}^{\prime }{|}^{\alpha +1}{D}^{\alpha +2}{\nu }^{-\alpha }$$

$$D=\frac{1}{\gamma (1-\mathbf{\beta }\cdot \mathbf{n})}$$

$$I(\nu ,Y,Z)={\int }_{-\infty }^{\infty }j(\nu ,X,Y,Z)dX$$

$$Q(\nu ,Y,Z)=\frac{\alpha +1}{\alpha +5/3}{\int }_{-\infty }^{\infty }j(\nu ,X,Y,Z)\cos 2\chi dX$$

$$U(\nu ,Y,Z)=\frac{\alpha +1}{\alpha +5/3}{\int }_{-\infty }^{\infty }j(\nu ,X,Y,Z)\sin 2\chi dX$$

Full text

*Aims. * We intend to extend our analysis by performing three-dimensional numerical simulations and investigating the role of different parameters is determining the formation of the X-shaped morphology.

*Methods. * The problem is addressed by numerical means, carrying out three-dimensional relativistic magnetohydrodynamic simulations of bidirectional jets propagating in a triaxial density distribution.

Results. * We show that only jets with power $\lesssim 10^{44}$ erg s-1* can give origin to an X-shaped morphology and that a misalignment of $30^{o}$ between the jet axis and the major axis of the density distribution is still favourable to the formation of this kind of morphology. In addition we compute synthetic radio emission maps and polarization maps.

*Conclusions. * In our scenario for the formation of X-shaped radio sources only low power FRII can give origin to such kind of morphology. Our synthetic emission maps show that the different observed morphologies of X-shaped sources can be the result of similar structures viewed under different perspectives.

1 Introduction

Fanaroff & Riley (1974) classified the extended radio sources on the basis of their morphologies and identified two classes, since then called Fanaroff-Riley I and II. The characteristic structure of FRII source is dominated by two hot spots located at the edges of the radio lobes that, in most cases, show bridges of emission linking the core to the hot spots. A significant fraction of FRII, about 10%, show a strong distortion in the bridges, that can be mirror symmetric (or C-shaped) when the bridges bend away from the galaxy in the same direction, or center symmetric, when they bend in opposite direction forming an X-shaped. In many X-shaped sources the radio emission along the secondary axis, although more diffuse, is still quite well collimated and can be even more extended than the main lobed structure. In the past 40 years, many authors suggested different interpretations of X-shaped radio galaxies, Ekers et al. (1978) suggested that the tails of radio emission in one of these sources, NGC 326, are the result of the trail caused by a secular jet precession (see also Rees 1978). A similar model accounts for the morphology of 4C 32.25 (Klein et al. 1995). In a similar line, Wirth et al. (1982) noted that a change in the jet direction can be caused by gravitational interaction with a companion galaxy. Dennett-Thorpe et al. (2002), from the analysis of spectral variations along the lobes, proposed that the jet reorientation occurs over short time scales, a few Myr, and are possibly associated to instabilities in the accretion disk that cause a rapid change in the jet axis. In all these models, the secondary axis of radio emission represents a relic of the past activity of the radio source. An alternative interpretation was suggested by Leahy & Williams (1984) and Worrall et al. (1995). They emphasize the role of the external medium in shaping radio sources, suggesting that buoyancy forces can bend the back-flowing material away from the jet axis into the direction of decreasing external gas pressure. Capetti et al. (2002) (Paper I) proposed a different scenario based on evidence for a strong connection between the properties of the radio emission and those of the host galaxy, a comparison that has been overlooked in the past, but that provides crucial new insights on the origin of X-shaped radio-sources. A more complete summary of the models proposed for the origin of X-shaped radiogalaxies can be found in Gopal-Krishna et al. (2012). In Paper I we compared the radio and host galaxy orientation of X-shaped radio sources showing that there is a close alignment between the radio wings and the galaxy minor axis. We also showed that X-shaped sources occur exclusively in galaxies of high ellipticity. These results have been confirmed by the analysis of a larger sample by Gillone et al. (2016). They selected, from an initial list of 100 X-shaped RGs from Cheung (2007), the 53 galaxies with the better defined wings and with available SDSS images; in 22 sources it was possible to measure the optical position angle. The orientation of the secondary radio structures shows a strong connection with the optical axis, with all (but one) wing forming an angle larger than 40o with the host major axis. Motivated by the observational evidences of the relation between the radio morphology and the optical galactic structure, in Paper I we presented the results of 2D numerical simulations of the propagation of a jet in a stratified medium with an elliptical symmetry. The simulations show that the radio source evolution and morphology is strongly influenced by the distribution of the external density and by the orientation of the jet propagation with respect to this distribution. If the jet propagates along the major axis of the density distribution, initially it inflates an almost spherical cocoon that becomes soon elongated in the jet direction, but also substantially expands along the minor axis due to the faster decrease of the external density in this direction. In this way we can observe the formation of a typical X-shape morphology. In Paper I we showed that, on the contrary, when the jet propagates along the minor axis we don’t observe the wing formation and we have a classical double radio source.

A limitation of the investigation performed in Paper I is the dimensionality. More recently Hodges-Kluck & Reynolds (2011) performed three-dimensional simulations and investigated the effect of a 3D elliptical stratification on the propagation of jets by using Newtonian hydrodynamical simulations. They emphasized the importance of ellipticity of the host galaxy and they stressed how the hot atmosphere may play a dominant role in shaping the morphology of radiogalaxies. In this paper, we intend to reconsider our previous work proposing a much more detailed numerical approach and extend the results obtained by Hodges-Kluck & Reynolds (2011) by including magnetic and relativistic effects, assuming a triaxial equilibrium density distribution of the external medium. Our approach differs from that of Hodges-Kluck & Reynolds (2011) also because they consider the long term evolution of the radio sources looking at the effect of switching-off the jet and of buoyancy forces, while we study the first part of the evolution and consequently also our spatial scales are smaller. Our purpose is to investigate how the formation of the wings depends on the jet parameters like their Lorentz factor $\gamma_{j}$ and their density ratio $\nu$ with respect to the external medium. Moreover we want also to explore the consequence of misalignment between jet axis and galaxy major axis.

In the next section we will describe the numerical setup and the set of simulations that we performed, in Section 3 we present our results focusing in particular on constraining the jet parameters for the formation of the wings and finally our conclusions are outlined in the last section.

2 Problem Description

Numerical simulations are carried out by solving the equations of relativistic magnetohydrodynamic

[TABLE]

where $\rho$ is the rest mass density, $u^{\mu}\equiv\gamma(1,\mathbf{v})$ and $b^{\mu}=(b^{0},\mathbf{B}/\gamma+b^{0}\mathbf{v})$ are respectively the four-velocity and covariant magnetic field written in terms of the three-velocity $\mathbf{v}$ and laboratory magnetic field $\mathbf{B}$. A flat metric $\eta^{\mu\nu}=\rm{diag}(-1,1,1,1)$ is considered. In the previous equations $w=w_{g}+b^{2}$ and $p=p_{g}+b^{2}/2$ express the total enthalpy and total pressure in terms of their thermal ($w_{g}$ and $p_{g}$) and magnetic field contributions ($b^{2}=b^{\mu}b_{\mu}$), respectively. We assume a single-species relativistic perfect fluid (the Synge gas) described by the approximated equation of state proposed by Mignone et al. (2005).

We study the propagation of bidirectional jets flowing into a triaxial ellipsoidal atmosphere initially isothermal and in hydrostatic equilibrium. The simulations are performed on a three-dimensional cartesian grid with the $y$ axis aligned with the jet direction. The computational domain is initially filled with a density distribution chosen as a $\beta-$model (Cavaliere & Fusco-Femiano 1976) deformed into a triaxial shape by taking different effective core radii along the three axes. The effective core radii $a_{x}$, $a_{y}$, $a_{z}$ and all the lengths are measured in unit of a fictitious core radius $r_{c}$, whose value can be chosen to be consistent with observational data. The initial density distribution is given by

[TABLE]

where $\rho_{0}$ is the central density,

[TABLE]

[TABLE]

are the coordinate along the axes of the ellipsoidal distribution and $\theta$ is the angle between the jet and the major axis of the atmosphere. A sketch of this setup is shown in Fig. 1. In all the simulations we used $a_{x}=1/3$, $a_{y}=2/3$ and $a_{z}=1/4$. The pressure distribution follows the density distribution and equilibrium is maintained by an external gravitational potential. Two jets, with initial radius $r_{j}=0.05r_{c}$, flowing in opposite directions are injected at the center of the domain. No magnetic field is present in the initial configuration at $t=0$, and a toroidal magnetic field is injected along with the jets. The injection region is a cylinder of radius $r_{j}$ and length $0.04r_{c}$ aligned with the $y$ axis in which the variables are not evolved but they are kept fixed at the injection values. As shown in Fig. 2 this region is divided in three part, a central part (shown in black in the figure) in which velocity and magnetic field are zero, and two side regions (shown in red and blue in the figure) in which the jet velocities are opposite and the magnetic field is defined as

[TABLE]

[TABLE]

in the blue region and is reversed in the red region. Here $r$ and $\phi$ are polar coordinates in a plane normal to the $y$ axis and the maximum field strength $B_{t}$ is fixed by the value of the jet magnetization, $\sigma$, which represents the ratio between the Poynting flux and the matter energy flux, in the jets. A jet with a purely toroidal field is potentially unstable to pinch and kink magnetohydrodynamic instabilities, however the effect of these instabilities starts to show up only at late times where some jet wiggling can be observed.

The two jets are slightly perturbed at their base, in order to destroy all symmetries, the following evolution is independent from the details of the perturbation. The parameters that characterize this configuration are the jet Lorentz factor $\gamma_{j}$, the density ratio $\nu$, between the jet density and the central density, the jet temperature and the magnetization $\sigma$, The temperature of the external medium is fixed by the jet Mach number $M_{j}$, which is 400 for all simulations. In Table 1 we show all the parameters for the simulations that we present in this paper. In the sixth column of Table 1 we give the values of the jet energy flux defined by

[TABLE]

and computed for the representative values $r_{c}=4{\rm kpc}$ and $\rho_{0}=1\hbox{cm}^{-3}$, which are consistent with observational data. The details of the injection setup have of course some influence on the characteristics of jet propagation. However, we can expect limited differences on the large scale structure, also in view of the fact that the perturbations imposed at the base will quite likely cause the jet to lose memory of the details of the initial state, while propagating into the ambient medium.

The first three simulations (A, B and C) have constant energy fluxes, but different and decreasing Lorentz factors $\gamma_{j}$. Cases A, D and E have fixed Lorentz factor $\gamma_{j}$, but decreasing density ratio and consequently decreasing energy flux. We have also to notice that cases D and E have a ten time higher magnetization than case A, the Poynting flux of case D is then equal to that of case A. Finally, cases $E$ and $F$ differ only for the angle $\theta$ between the jet direction and the major axis of the atmosphere. In Table 1 we give, in the last two columns, for each case, the size of the computational domain and the grid size. For all the simulations the grid is uniform in the central region ($1\times 2\times 1$) and is stretched in the outer parts in order to achieve larger distances. The simulations have been performed with the PLUTO code (Mignone et al. 2007), with a second order accurate linear reconstruction and HLL Riemann solver.

3 Results

In the presentation of the results, all quantities are expressed in physical units, their dimensional values depend however on the choice of our unit of length $r_{c}$ and density $\rho_{0}$, which we have assumed $r_{c}=4\rm{kpc}$ and $\rho_{0}=1\rm{cm}^{-3}$. A different choice for these units will give different values for all the relevant physical quantities, more precisely lengths scale as $r_{c}$, times scale also as $r_{c}$, pressure scales as $\rho_{0}$ and magnetic field strength scales as $\sqrt{\rho_{0}}$.

In the first three simulations (A, B and C) we kept a constant jet power of $1.6\times 10^{46}$ erg s${}^{-1},$ but we varied the Lorentz factor and the jet density, keeping $\gamma_{j}^{2}\rho_{j}$, and hence the energy flux, constant. In all cases, at the injection point, the jet results to be overpressured with respect to the ambient medium, with increasing pressure going form case A to case C. In cases D and E, instead, we kept a constant Lorentz factor but we decreased the jet density and therefore the jet power which results to be $1.6\times 10^{45}$ erg s*-1* in case D and $1.6\times 10^{44}$ erg s*-1* in case E.

In Fig. 3, the panels in the central column show, for cases A to E, 2D cuts in the $yz$ plane of the logarithmic density distribution at the time when the jets have reached a length of approximately 8 kpc. We recall that the $z$ direction corresponds to the minor axis of the triaxial distribution. The figures give a first impression of the relative importance of lateral wings in the different cases. The formation of a winged structure depends on the lateral expansion velocity of the cocoon relative to the propagation velocity of the jet head. A theoretical estimate of the latter can be derived by equating the momentum flux of the jet and of the ambient medium in the frame of reference of the jet head. In the relativistic case one can derive the following expression (Martí et al. 1997):

[TABLE]

where

[TABLE]

and $\rho_{j}$, $\rho_{a}$, $h_{j}$, $h_{a}$ are respectively the jet and ambient densities and the jet and ambient specific enthalpies. The ambient density decreases as the jet propagates outward according to Eq. (4) and, according to Eq. (10), the jet head should accelerate. In the right column of Fig. 3, we can compare the actual jet head position as a function of time, represented by the black curves, with the theoretical estimates given by Eq. (10), represented by the green curves. In the same figures we also show the size reached by the cocoon as a function of time, in the two transverse direction, the size along the $z$ direction is represented in blue, while the size along the $x$ direction is represented in red. Considering cases A, B and C, the initial predicted jet head velocities are the same, since $\gamma_{j}^{2}\rho_{j}$ is kept constant, and equal to about $2.1\times 10^{9}$ cm s*-1*. The actual values are in very good agreement with this estimate. These three cases, according to Eq. (10) should behave in the same way, however, in reality, they show marked differences. In case A the jet head accelerates reaching a velocity of about $4.3\times 10^{9}$ cm s*-1* when the jet extends up to 8 kpc. The agreement between the actual and estimated velocities is very good. In case B we still observe acceleration and the velocity at 8 kpc is about $3.15\times 10^{9}$ cm s*-1*, the comparison with the theoretical estimate shows an actual lower velocity, in fact the estimated velocity at 8 kpc should be as before $4.3\times 10^{9}$ cm s*-1*. This effect is much more pronounced in case C, where, instead of accelerating, the jet head decelerates, reaching a velocity of about $8.7\times 10^{8}$ cm s*-1* at 8 kpc. According to Eq. (10), the jet head velocity, in our conditions, scales as $\sqrt{\nu}$ and therefore the initial values reduce to $6.7\times 10^{8}$ cm s*-1* in case D and to $2.1\times 10^{8}$ cm s*-1* in case E. In case D, we still observe an acceleration of the jet head, with, however, a velocity lower than expected (compare the black and green curves), with an actual value of about $1.6\times 10^{9}$ cm s*-1*. In case E the deviation from the theoretical estimate is more evident and the acceleration barely visible (compare the black and green curves), in fact the velocity at 8 kpc is only $3.13\times 10^{8}$ cm s*-1*. Clearly the jet reaches a size of 8 kpc at progressively later times in cases D and E. The discrepancy between the actual and estimated velocity can be understood by considering that, in deriving Eq. (10), it is implicit the assumption that the two areas used for balancing the momentum fluxes are equal (Lind et al. 1989; Massaglia et al. 1996), in reality the area on which the jet deposits its momentum can increase leading to a velocity lower than expected from Eq. (10). The increase in the area can occur either because of an expansion of the working surface, particularly evident in case C, or because of instabilities that can fragment or induce wiggling of the jet head. In Figure 4 we compare the structures of the jet head for case A, in which the actual head velocity perfectly agrees with the theoretical estimate, and for case E, in which the actual head velocity is well below the theoretical estimate. The top panels refer to case A, while the bottom panels refer to case E, left panels show a cut of the pressure distributions in the plane $yz$, i.e. the plane defined by the major and minor axes, while right panels show a cut of the longitudinal velocity distributions in the same plane. We can see that while in case A all the jet thrust is concentrated on a section which is comparable to the jet width, in case E the jet close to the head appears to have widened and the thrust is distributed on a larger area. This is consistent with the discussion made above and the effect appear to become more pronounced as we decrease the density ratio.

Concerning the wing expansion, their velocity can be estimated as (Begelman & Cioffi 1989):

[TABLE]

where $p_{c}$ is the average cocoon pressure, whose evolution as a function of the jet length is represented in Fig. 5, for all cases. We see that the cocoon pressure decreases as the system evolves and its size increases. The wings expansion velocity which depends on the ratio between $p_{c}$ and $\rho_{a}$, will increase or decrease depending on which between pressure and density decreases faster. Our results displayed in the right column of Fig. 3 show that the lateral expansion velocities for all cases remain approximately constant in time (blue and red curves), the decrease in the cocoon pressure is then almost balanced by the decrease in the external density.

In cases A, B and C, the cocoon pressure has approximately the same value and in fact the lateral expansion velocities are almost the same for all these cases, approximately $10^{9}$cm s*-1*, with the velocity along the $x$ direction about $10\%$ lower than the velocity along the $z$ direction. Notice then that the higher wing prominence in case C is due to a lower jet head velocity and not to a higher lateral velocity. The results for cases A, B and C show that for getting a winged source at high power the jets have to be subrelativistic. For relativistic jets we can get X-shape morphology only decreasing their power.

The decrease of the jet power leads to a decrease of the cocoon pressure, as shown in Fig. 5, and therefore to a slowing down of the wing expansion. The lateral velocities reduce to about $5\times 10^{8}$ cm s*-1* in case D and to about $2.5\times 10^{8}$ cm s*-1* in case E. However, as seen above, also the jet head velocity reduces with the jet power and the wing prominence depends on which reduces faster. In Fig. 6 we compare, for cases A, D and E, the average expansion velocities along the three axes as a function of jet power. Circles represent the jet head velocity, while diamonds and squares represent the lateral expansion velocities respectively along the $z$ and $x$ directions. The figure shows that the jet head velocity decreases with jet power faster than the lateral velocities and moreover the expansion velocity along the minor axis $z$ decreases slower than the corresponding velocity along the intermediate $x$ axis. Therefore, in conclusion, only low power relativistic jets can give origin to X-shaped morphologies, while high power jets have to be subrelativistic, as shown by cases A, B and C discussed above. Based on the comparison between the velocities along the $z$ and $x$ axes, shown in Fig 6, we expect that, increasing the ellipticity, the decrease of the lateral expansion velocity along the minor axis will be lower and the critical power, under which a winged morphology can be obtained, will depend on the ellipticity of the gas distribution, with a higher critical power for higher ellipticity.

3.1 Case F: The effects of misalignment between jet and major axis

In this subsection we will describe the evolution of Case F, which has the same parameters as case E, but has a misalignment of $30^{o}$ between the jet axis and the major axis of the density distribution and which we followed for a somewhat longer time, until it reached a jet length of about 12 kpc at $t\sim 3.6\times 10^{6}$yrs. The most favourable condition for the formation of the wings is when the jet is aligned with the major axis of the galaxy, with case F we want then explore whether the X-shape morphology can develop also in a less favourable condition. Fig. 7 compares the evolution of case E and case F by showing three-dimensional volume renderings of the distribution of the $y$ component of the velocity, with blue colors representing negative values and red-yellow colors representing positive values. The renderings show in particular the structure of the backflows, i.e. of the jet material that after crossing the jet terminal shock flows back towards the jet origin. Superimposed, we also show a cut of the density distribution in the $x-y$ plane, while the blue lines are density iso-contours at $\rho=0.2$cm*-3* on the same plane. The density distributions and the iso-contour lines allow to follow how the initial ellipsoidal density structure is progressively modified by the inflating cocoon. The panels on the left are for case E, while the panels on the right are for case F. The three panels on each column are for three different times showing the cocoon evolution. Finally in Fig. 8 we show only the structure of case F at a later time, for which we don’t have the corresponding image for case E, since, for this case we stopped the simulation at an earlier time. The figure shows how the cocoon progressively inflates expanding in the direction of the minor axis of the density distribution, which in case F is not perpendicular to the jet axis. As shown in the figure, the times at which case F reaches the same length as case E are somewhat earlier, indicating a slightly larger jet head velocity. In fact the average velocity for case F is $3.26\times 10^{8}$ cm s*-1* to be compared to $2.8\times 10^{8}$ cm s*-1*. The lateral expansion velocity is, on the contrary, slightly lower in case F than in case E. The initial density distribution, that at the beginning strongly influences the cocoon evolution, is completely modified by the jet propagation, which inflates a cavity and forms a compressed shell at the border of the cavity. In addition we can follow the interaction of the two backflows, with one of them always prevailing on the other and breaking, in case E, the initial symmetry. In case F the jets form non symmetric backflows that bend towards the minor axis of the initial density distribution flowing up in the wings, as it can be seen in Fig. 8. In case F the wing axis form an angle of about $70^{\circ}$ with the initial jet axis, almost constant with time, and slightly larger than the angle of $60^{\circ}$ formed by the jet axis with the minor axis of the density distribution.

Fig. 9 shows the distribution of pressure in the two planes $(xy)$ and $(yz)$, when the radio source has reached the maximal size at $t\sim 3.2\times 10^{6}$ yrs, with superimposed the velocity field. The red ellipses represent an equicontour of the initial density distribution. In the two panels we show respectively cuts in the $x-y$ plane (left panel) and in the $y-z$ plane (right panel). The two panels allow to observe the effect of the misalignment on the final geometrical structure, already evident in Fig. 7. In addition we can see that the cocoon is still strongly overpressured with respect to the ambient. Looking at the velocity field, whose streamlines are visible in the same figure, we can observe the presence of flows rising in the wings.

The conclusion of the analysis of case F is that the prominence of wings is lower than in case E, as could be expected, since the jet direction has been moved away from the most favourable condition of perfect alignment with the major axis of the initial density distribution. This inclination of $30^{o}$, however, is still favourable for the formation of a winged source.

Summarizing, we showed that jets with a power larger than $10^{45}$ erg s*-1* have to be non relativistic in order to form a winged structure, while relativistic jets can form a winged structure only if they have $L_{j}\lesssim 10^{44}$ erg s*-1*. Bîrzan et al. (2004) estimated the jet power from the observations of the X-ray cavities inflated by radio AGN. The ratio between $L_{j}$ and the radio power at 1.4 GHz ($P_{1.4}$) ranges from $\sim 3$ to $\sim 300$ and this leads to an estimate of $P_{1.4}\lesssim 3\times 10^{43}$ erg s*-1* for the maximum radio luminosity of winged sources; this value compares favourably with the power distribution of X-shaped radio sources, that reaches $P_{1.4}\sim 2\times 10^{43}$ erg s*-1* (Gillone et al. 2016).

3.2 Magnetic field structure

An important question to be analysed is whether the flows rising in the wings are able to carry the magnetic field necessary for sustaining the synchrotron radiation. In order to answer this question, we first examine in Fig. 10 the distribution of the magnetic field strength, for case E at time $t=2.6\times 10^{6}$ yrs, in the $xy$ and $yz$ planes, corresponding, respectively, to the planes of the major and minor axes ($xy$, left panel) and of the major and intermediate axes ($yz$, right panel). The field strength is expressed in microgauss and is computed by assuming for the basic physical units $r_{c}$ and $\rho_{0}$ their reference values, recalling that the field strength scales as $\rho_{0}^{1/2}$. From the figure we see that the field strength in the jets reaches maximum values between $30\mu$G and $60\mu$G, similar values can be also found in the wings in the $xy$ plane and in the lower part of the wings in the $yz$ plane, while the upper part of the wings in the $yz$ plane shows lower values. One question that immediately arises is whether the magnetic field becomes dynamically important, we have then examined the distribution of the local values of the plasma $\beta$ representing the ratio between thermal and magnetic pressure and we have seen that there are only a few points where it reaches values of the order of unity, while for the most part it has values much larger than unity and therefore it is generally dynamically unimportant. A more quantitative comparison between the magnetic field strength in the jets and in the wings is provided by Fig. 11 where we plot the distribution functions of the magnetic field strength in the jet and in the wings. The jet is defined as a cylinder of radius $0.3r_{c}$ with its axis aligned with the $y$ axis, while the wings are defined as the region external to this cylinder. From the figure we see that the maximum of the distribution is about $15\mu$G in the jets and about $5\mu$G in the wings, the medians of the distributions are however much closer, since they are again $15\mu$G in the jets, but about $10\mu$G in the wings. This means that half of the volume in the jets has field strength larger than $15\mu$G, while in the wings, in order to get half of the volume, we have to lower the threshold to $10\mu$G. Regions with field strength larger than $30\mu$G have very small filling factors, of about $7\%$ in the jets and $1.5\%$ in the wings. The distribution of magnetic field described above represents a snapshot of the evolution at $t=2.6\times 10^{6}$ yrs, at the end of the simulations, if we look at the complete temporal evolution, we can see that the spatial distribution and the distribution function of the field strength do not change much in time, however the average field strength decreases with time as represented in Fig. 12. The value of the field depends of course on the value of $\sigma$, representing the ratio of the Poynting flux to the matter energy flux, in this case $\sigma=10^{-2}$. Note that the values of the resulting field amplitude are quite consistent with observations. Increasing $\sigma$ we expect that the magnetic field scale as $\sqrt{\sigma}$ until it reaches values for which it becomes dynamically important. The regime in which the magnetic field becomes dynamically important is not examined in this paper, however it could be of interest for further investigations.

3.3 Synthetic radiomaps

In Fig. 13 we display synthetic surface brightness and polarization maps for different viewing angles. The calculations of the surface brightness and of the polarization are described in the appendix A. In the left panels we show the brightness maps, while in the right panels we show the polarization maps, plotting the magnetic field direction. The line of sight is in the $xy$ plane and from top to bottom makes with the jet direction angles respectively of $90^{\circ},45^{\circ}$ and $20^{\circ}$. Decreasing the viewing angle, the size of the wings becomes more prominent as the projected jet length decreases, until, at the smallest viewing angles, the main lobes are projected onto the wings.

We can compare the synthetic maps of the surface brightness distribution in Fig. 13 with radio maps of X-shaped sources. The synthetic map in the top-left panel of Fig. 12, viewing angle of $90^{\circ}$ is reminiscent, for instance, of the radio map 3C 136.1 (Fig. 14, left panel, from Lal & Rao 2007). Considering the middle-left panel of Fig. 12, viewing angle of $45^{\circ}$, the synthetic map can be compared to the morphology displayed by the source J1434+5906 (from Roberts et al. 2015) and shown in Fig. 14 (middle panel). The bottom-left panel of Fig. 12, viewing angle of $20^{\circ}$, compares nicely with the source J1043+3131 (again from Roberts et al. (2015)) and shown in Fig. 14 (right panel). Within our model, the variety of morphologies displayed by X-shaped sources can the be interpreted, beside the general effects of the the main physical parameters, by the combination of different viewing angles and the orientation of the jet direction with respect to the major axis of the host galaxy. In particular, sources with a morphology similar to J1043+3131 are not usually classified as classical X-shaped, however our 3D results show that X-shaped sources can indeed take this morphology if the line of sight form a small angle with the jet direction.

From our polarization maps we see that the magnetic field is directed along the jets and along the wings, becoming transverse only at the edges. This kind of morphology can be observed also in actual sources, as we can see in Fig. 15, where we show the polarization map of 3C 223.1 (from Dennett-Thorpe et al. (2002), left panel) and the one of 3C 136.1 (from Rottmann (2001), right panel). These maps have to be compared with the right panels of Fig. 13.

The emission here is treated in a simplified way and, since we do not model the evolution of the emitting non-thermal particles, we have to give a prescription for the their density and distribution. Different prescriptions may lead to some differences in the synthetic maps. We have then considered an alternative in which we prescribe local equipartition between the magnetic energy and the energy of non-thermal particles. In Fig. 16 we show, for this choice, the synthetic map of the surface brightness distribution for line of sight at $90^{\circ}$ with the jet direction, corresponding to the top left panel in Fig.13. Comparing the two maps, we can observe an overall similarity, the main differences are in the jets and in the central part of the wings which show a somewhat higher brightness.

A final remark concerns the distribution of the offset between the host major axis and the wings. Capetti et al. (2002) and Gillone et al. (2016) found that in only one source (out of 31) the wings form an angle smaller than 40o and this is the key observing result that motivated the simulations presented in this paper. However, this finding carries further information on the wings structure. In fact, in case the wings were linear and essentially one dimensional features, a significant fraction of X-shaped sources should show, due to projection effects, small misalignments; this would occur when the plane defined by the wings axis and the major axis of the host forms a large angle with the plane of the sky. Our simulations show instead that the wings form a flattened emitting region aligned with a plane including the host minor axis; the projection of such essentially two-dimensional structure leads to a narrower range of offsets, as observed.

4 Summary

We presented 3D relativistic magnetohydrodynamic simulations of bidirectional jets propagating in a triaxial density distribution in order to investigate the mechanism of formation of X-shaped radiogalaxies. In Paper I we proposed a hydrodynamic mechanism for the formation of this class of sources and we performed 2D simulations of the propagation of a jet in an ellipsoidal density distribution. We showed that, when the jet propagates along the major axis of the distribution a winged structure can form. Here we extended our previous results, relaxing the axial symmetry, introducing magnetic field and relativistic effects. We investigated, in particular, under which conditions of jet power, Lorentz gamma factor, and density ratio, an X-shaped structure arises. We showed that jets with a power larger than $10^{45}$ erg s*-1* have to be non relativistic in order to form a winged structure, while relativistic jets can form a winged structure only if they have $L_{j}\lesssim 10^{44}$ erg s*-1*. By applying a scaling relation between jet and radio power, we estimated that the maximum radio luminosity of a winged source is $\sim 3\times 10^{43}$ erg s*-1*, a value compatible with the observations of X-shaped radio sources. Moreover we considered a case in which the jet is misaligned with respect to the major axis of the density distribution and showed that an angle of $30^{o}$ is still favourable to the formation of an X-shaped structure.

Hodges-Kluck & Reynolds (2011) also performed three-dimensional simulations for a configuration similar to ours, they did not include magnetic field and relativistic effects. An important difference between their analysis and ours is the range of spatial and temporal scales considered. We focus on the first phases of the evolution limiting to time scales of the order of $10^{6}$ yrs and spatial scales of the order of $10$ kpc, when the jet still propagates within the galactic atmospheres, while they consider the full life of the radiosource and study its evolution when the jet is switched-off, the spatial scales are also larger, and the jet propagates in the cluster environment. In this sense the two approaches can be considered as complementary: the first phases of the evolution, that we consider, are essential for the formation of an X-shaped radiosource and the conditions that we find are necessary conditions, while the results obtained by Hodges-Kluck & Reynolds (2011) give important insights on the subsequent long term evolution.

We have also calculated the brightness distribution producing synthetic emission maps. Being the simulated sources intrinsically three-dimensional, i.e. displaying a triaxial structure, the same objects may appear differently as seen from varying viewing angles. We show that different observed morphologies of X-shaped sources can be the result of similar structures viewed under different perspectives. By considering both the structure of the wings and the triaxiality of their hosts, the relative orientation of radio and optical structures has a complex dependence on the projection; this leads not only to favour large offsets between the wings and the host major axis, but also to the observed paucity of offset smaller than $\sim 40^{o}$.

Finally we remark that we carried out a simplified treatment of the emission. In reality, one should consider the evolution of a distribution of relativistic electrons, carried along the jets and streaming in the wings, that lose energy via synchrotron losses and adiabatic expansion, and gain energy by shock acceleration and/or via turbulent acceleration by IInd order stochastic Fermi-like process (see Vaidya et al. 2017 to be submitted). Our results show that shocks in the wings are quite weak and therefore they will most likely contribute in a limited way to particle reacceleration, however reacceleration by turbulence may be more effective. Further studies are then required to examine the behavior of the spectral index in the wing, compared to observations.

Acknowledgments

We acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support. This work has been partly funded by a PRIN-INAF 2014 grant.

Appendix A Surface brightness and polarization maps

The surface brightness radio maps and polarization maps have been obtained by integrating along the line of sight respectively the synchrotron emissivity $j$ and the Stokes parameters $Q$ and $U$. The synchrotron emissivity in the comoving frame, for a power law distribution of relativistic electrons is (Del Zanna et al. 2006)

[TABLE]

where all primed quantities are computed in the proper frame, $\mathbf{n}^{\prime}$ represents the versor corresponding to the line of sight direction, $K$ is a constant proportional to the number density of relativistic electrons and $-(2\alpha+1)$ is the slope of the energy distribution of relativistic particles. In order to obtain the emissivity in the observer frame, we have to take into account relativistic corrections and we get

[TABLE]

where $D$ is the Doppler factor

[TABLE]

Considering a cartesian reference frame in which $X$ is along the line of sight and $Y$ and $Z$ are in the plane of the sky, the surface brightness can be obtained as

[TABLE]

and the Stokes parameters $Q$ and $U$ are given by

[TABLE]

[TABLE]

where $\chi$ is the local polarization position angle, measured clockwise from the $Z$ axis. From $Q$ and $U$ we can then get the observed polarization direction. In our calculations we assumed the density of relativistic electrons to be proportional to the local fluid density and the spectral index $\alpha=0.5$.