Jet base synchrotron radiation

Synchrotron radiation of self-collimating relativistic MHD jets

Abstract

The goal of this paper is to derive signatures of synchrotron radiation from state-of-the-art simulation models of collimating relativistic magnetohydrodynamic (MHD) jets featuring a large-scale helical magnetic field. We perform axisymmetric special relativistic MHD simulations of the jet acceleration region using the PLUTO code. The computational domain extends from the slow magnetosonic launching surface of the disk up to Schwarzschild radii allowing to reach highly relativistic Lorentz factors. The Poynting dominated disk wind develops into a jet with Lorentz factors of and is collimated to . In addition to the disk jet, we evolve a thermally driven spine jet, emanating from a hypothetical black hole corona. Solving the linearly polarized synchrotron radiation transport within the jet, we derive VLBI radio and (sub-) mm diagnostics such as core shift, polarization structure, intensity maps, spectra and Faraday rotation measure (RM), directly from the Stokes parameters. We also investigate depolarization and the detectability of a -law RM depending on beam resolution and observing frequency. We find non-monotonic intrinsic RM profiles which could be detected at a resolution of 100 Schwarzschild radii. In our collimating jet geometry, the strict bi-modality in polarization direction (as predicted by Pariev et al.) can be circumvented. Due to relativistic aberration, asymmetries in the polarization vectors across the jet can hint to the spin-direction of the central engine.

Subject headings:
galaxies: active - galaxies: jets - ISM: jets and outflows - magnetohydrodynamics (MHD) - radiation mechanisms: non-thermal - relativistic processes

1. Introduction

Relativistic jets are launched from accretion disks around compact objects and accelerated and collimated by magnetohydrodynamic forces (Blandford & Payne, 1982). The disk jet flow may be accompanied by a highly relativistic central spine jet resulting from electrodynamic effects within the black hole magnetosphere (Blandford & Znajek, 1977).

The formation of relativistic MHD jets has been investigated by a number of authors since the seminal paper of Blandford & Znajek. The early attempts were to look for stationary state MHD solutions for the asymptotic jet structure (Chiueh et al., 1991; Appl & Camenzind, 1993; Fendt, 1997b), or the jet formation domain of collimation and acceleration in case of self-similarity (Contopoulos, 1994, 1995), or taking into account the 2.5D force-balance in case of special relativity (Camenzind, 1987; Fendt & Memola, 2001), or general relativity (Fendt, 1997a; Ghosh, 2000; Fendt & Greiner, 2001).

Time-dependent simulations of jet formation from the disk surface have first been investigated in the non-relativistic approximation (Ustyugova et al., 1995; Ouyed & Pudritz, 1997). General relativistic MHD simulations of accretion disks launching outflows (Nishikawa et al., 2005; McKinney & Narayan, 2007; Beckwith et al., 2008; McKinney & Blandford, 2009), indicate that highly relativistic jets may not be launched by the disk itself, but from the black hole magnetosphere. These jets - so-called funnel flows - may reach Lorentz factors up to 10 however, the mass loading is put in by hand (maybe corresponding to a pair-creation process) and is not related to the accretion process itself. Komissarov et al. (2007) proposed special relativistic MHD simulations of AGN jets over a large spatial scale and find asymptotic Lorentz factors of about 10. These jets were pressure-confined by the outer boundary condition.

In a recent paper we have treated another setup of relativistic jet formation (Porth & Fendt, 2010). We have applied special relativistic MHD to launch outflows from a fixed-in-time hot surface of an accretion disk which is rotating in centrifugal equilibrium with the central compact object. This allowed us to pass through the slow-magnetosonic point, thus obtaining consistent mass and Poynting fluxes, and to investigate the subsequent acceleration and collimation degree of a variety of jets. Essentially, we find these relativistic MHD jets to be truly self-collimating - similarly to their non-relativistic counterparts. In this paper, we will prescribe certain mass fluxes and explore the influence of the poloidal current distribution on the jet formation process. In contrast to the previous work, the jet energy and its partitioning is now a true parameter of the models, allowing us to investigate the regime of highly relativistic flows. More importantly, observational signatures of the jet formation models are discussed.

The existence of an ordered large-scale G- mG magnetic field in extragalactic jets is well established by the detection of radio synchrotron emission, however the exact geometry of the field structure cannot be derived. The recent observational literature strongly suggests helical magnetic fields (e.g. O’Sullivan & Gabuzda (2009)), a scenario which is consistent with theoretical models and numerical simulations of MHD jet formation (see e.g. (Blandford & Payne, 1982; Ouyed & Pudritz, 1997; Krasnopolsky et al., 1999; Fendt & Čemeljić, 2002; Fendt, 2006; Porth & Fendt, 2010)). The wavelength dependent rotation of the polarization plane known as Faraday rotation provides a valuable diagnostic of magnetic field structure in astrophysical jets. Consistent detections of are found in resolved jets as well as in unresolved radio cores (Zavala & Taylor, 2003). Helical magnetic fields are generally perceived to promote transversal Faraday rotation measure (RM) gradients owing to the toroidal field component. Observationally, such gradients were first detected by Asada et al. (2002) and Zavala & Taylor (2005) in the jet of . The RMs are generally found to follow a monotonic profile across jets (Gabuzda et al., 2004; O’Sullivan & Gabuzda, 2009; Croke et al., 2010).

While the magnetic field structure of pc and Kpc-scale jets can in principle derived from radio observations, not much is known about the field structure within the jet forming region very close to black hole - mainly because of two reasons. Firstly, this very core region cannot be spatially resolved and, thus, cannot readily be compared with expected RM profiles of helical jets. This is somewhat unfortunate as close to the jet origin we expect the magnetic field helix to be well preserved and not much affected by environmental effects. Secondly, the observed rotation measures are so high that it is impossible to draw firm conclusions about the intrinsic field geometry from the polarization vector with mere radio observations.

Only very few cases exist where this core of jet formation could be resolved observationally. Among them is the close-by galaxy M87, where the VLBI/VLBA resolution of  mas is sufficient to resolve about 0.01 pc within the central region, and allows to trace the jet origin down to Schwarzschild radii when the recent mass estimate of by Gebhardt et al. (2011) is adopted. This pinpoints the launching area within (Junor et al., 1999; Kovalev et al., 2007). The radio maps clearly show limb brightening and indicate an initial jet opening angle of about .

Despite a vast amount of observational data spanning over a huge frequency scale, and also time series of these multifrequency observations, rather little is known about the dynamical status of relativistic jets. There are no direct unambiguous observational tracers of jet velocity or density as only (if at all) the pattern speed of radio knots is detected. Kinematic modelling of knot ejection suggests pc-scale Lorentz factors of typically , while Kpc-jet velocities are believed to be definitely lower and of the order of c. Kinematic modeling of jet propagation has been combined with synchrotron emission models of nuclear flares, resulting in near perfect fits of the observed, time-dependent radio pattern of jet sources such as 3C 279 (Lindfors et al., 2006).

However, what was missing until recently is a consistent combination of dynamical models with radiation models of synchrotron emission resulting in theoretical radiation maps which can then be compared with observations. Zakamska et al. (2008) and Gracia et al. (2009) have taken a step into this direction by providing optically thin synchrotron and polarization maps from self-similar MHD solutions. Broderick & McKinney (2010) presented synchrotron ray-castings from 3D general relativistic jet formation simulations that also include the evolution of a turbulent accretion disk performed by McKinney & Blandford (2009). In their approach, the MHD solution is extrapolated by means of an essentially self-similar scheme in order to reach distances up to . Their study focussed on the rotation measure provided by Faraday rotation in the disk wind external to the emitting region in the Blandford-Znajek jet.

In comparison, using axisymmetric large scale simulations, we do not rely on an extrapolation within the AGN core (up to ) and treat the Faraday rotation also internal to the emitting region in the fast jet that gradually transforms into a sub-relativistic disk wind. The observational signatures obtained in our work are derived entirely from the (beam convolved) Stokes parameters which allows us to also investigate the breaking of the law due to opacity effects.
As the above studies, we also rely on post-hoc prescriptions for the relativistic particle content. Towards a more consistent modeling of the non-thermal particles, Mimica & Aloy (2010) have presented a method to follow the spectral evolution of a seed particle distribution due to synchrotron losses within a propagating relativistic hydrodynamic jet. The question of particle acceleration and cooling is in fact essential to close the loop for a fully self-consistent treatment of jet dynamics, jet internal heating, and jet radiation.

To derive signatures of synchrotron emission from relativistic MHD jet formation is the main goal of the present paper. We apply the numerically derived dynamical variables such as velocity, density, temperature, and the magnetic field strength & configuration to calculate the synchrotron emission from these jets, taking into account proper beaming and boosting effects for different inclination angles. In particular, we apply the relativistically correct polarized radiation transfer along the line of sight throughout the jet.

2. The relativistic MHD jet

We perform axisymmetric jet acceleration simulations with the PLUTO 3.01 code (Mignone et al., 2007) solving the special relativistic magnetohydrodynamic equations. As in Porth & Fendt (2010), the simulations are of the disk-as-boundary type, however in the present study the jet starts out with the slow magnetosonic velocity opening the freedom to assign all energy channels as input parameters. The energy of the jet base is dominated by Poynting flux, driving a large-scale poloidal current circuit. This current distribution is prescribed as a boundary condition in the toroidal magnetic field component. We investigate three cases for with , where is the cylindrical radius, resulting in asymptotic jets which are either in a current-carrying (one of them) or current-free (two others) configuration. Injected initially with slow magnetosonic velocity, the jet material accelerates through Alfvèn and fast magnetosonic surfaces within the simulated domain.

2.1. Numerical grid setup

To eradicate any artificial collimation effect from the outflow boundaries, we decided to move the domain boundaries to such large distance that they are out of causal contact with the region of interest. An inner equidistant grid of cells extends to the scale radius corresponding to the inner radius of the accretion disk. Beyond the grid is linearly stretched to with a scaling factor of , and for with a larger scaling factor of . So far, we apply a square box of scale radii corresponding to grid cells. Magnetic fields are advanced on a staggered grid using the method of constrained transport (Balsara & Spicer, 1999) supplied by PLUTO.

For the calculation of the radiation maps, we consider only a subdomain of grid cells corresponding to scale radii.

The first possible contamination by boundary effects takes place when the bow shock reaches the upper outflow boundary after . A typical simulation is terminated after maximally time units when a signal traveling at the speed of light could have returned from the boundary to the subdomain. This conservative treatment ensures that no spurious boundary effects can occur.

In order to test this we have re-run one of our simulations on a five times larger computational domain with comparable resolution. This simulation provided the same result as the lower grid size simulation and acquired a nearly stationary state when terminated. Thus, the region of interest is not at all affected by any outflow boundary effects on account of an increased computational overhead.

Inflow boundary conditions

In a well posed MHD boundary, the number of outgoing waves (i.e. seven minus the downstream critical points) must equal the number of boundary constraints provided. Thus, in addition to the and profiles, we choose to prescribe the thermal pressure and density as boundary conditions for the jet injection (the jet inlet). The fifth condition sets and thus constrains the toroidal electric field . This is a necessary condition for a stationary state to be reached by the axisymmetric simulation.

The remaining primitive MHD variables and are extrapolated linearly from the computational domain, while the component (which determines the magnetic flux) follows from the solenoidal condition.

Specifically, the fixed profiles read

(1)
(2)
(3)
(4)

where

(5)

is the step function and and denote the cylindrical and spherical radius, respectively. Here and in the following, a subscript 1 indicates the quantity to be evaluated at .

The slow magnetosonic velocity profile along the disk of relation (4) is updated every time step to account for the variables which are extrapolated from the domain. Within , an inner “black hole corona” with relativistic plasma temperatures is modeled. The physical processes responsible for the formation of an inner hot corona could be an accretion shock or a so-called CENBOL shock (CENtrifugal pressure supported BOundary Layer shock) (see also Kazanas & Ellison, 1986; Das & Chakrabarti, 1999). Also Blandford (1994) has proposed a mechanism of dissipation near the ergosphere as a consequence of the Lense-Thirring effect. Due to the large enthalpy and decreasing Poynting flux of the inner heat bath, a thermally driven outflow is anticipated from this region (e.g. Meliani et al., 2010). We apply the causal equation of state introduced by Mignone et al. (2005) to smoothly join the relativistically hot central region to the comparatively cold disk jet. This choice permits a physical solution for both regions of the flow (see also Mignone & McKinney, 2007).

Outflow boundary conditions

At the outflow boundaries we apply power-law extrapolation for and the parallel magnetic field component, while we apply the solenoidal condition to determine the normal magnetic field vector respectively . For the velocities and the toroidal field, zero-gradient conditions are applied. This choice is particularly suited to preserve the initial condition that is well approximated by power laws. A more sophisticated treatment such as force-free (introduced by Romanova et al., 1997) or zero-current boundaries as in Porth & Fendt (2010) is rendered unnecessary by the increased computational grid as detailed before.

2.2. Initial conditions

As initial setup we prescribe a non-rotating hydrostatic corona threaded by a force-free magnetic field. For the initial poloidal field distribution we adopt and (Ouyed & Pudritz (1997), see also Porth & Fendt (2010)). The hydrostatic corona is balanced by a point-mass gravity in a Newtonian approximation. In order to avoid singularities in the density or pressure distribution, equations (1), (2), we slightly offset the computational domain from the origin by .

2.3. Parametrization

In the present setup we focus on effects of the poloidal current distribution and parametrize accordingly. Thus, we set the Kepler speed at to for convenience, yielding a sound speed from the hydrostatic condition. To further minimize the number of free parameters, we tie the toroidal field strength given by to the poloidal field strength via . This Ansatz is consistent with the sub-Alfvèninc nature of the flow, since at the Alfvèn point (e.g. Krasnopolsky et al., 1999) is valid. The two remaining parameters are the poloidal magnetic field strength measured by the plasma-beta and the toroidal field profile power law index .

Note that with the choice of a fixed in time toroidal field, the injected Poynting flux is controlled via the parameter, implying the toroidal field being induced in a non-ideal MHD disk below the domain. The conserved total jet energy flux and its partitioning between Poynting and kinetic energy given by the -parameter essentially become boundary conditions and determine the terminal Lorentz factor.

2.4. Physical scaling

In order to convert code units to physical units we need to define the two scaling parameters of length and density, while the velocity is naturally normalized to the speed of light . To obtain an approximate radial scale, we assume that the transition between the inner corona to the disk-driven jet at corresponds to the innermost stable circular orbit at (for a Schwarzschild black hole). Hence the physical length scale is given by

(6)

The physical density is then obtained by assuming a total jet power in units of ,

(7)

where is the corresponding power in code units for a certain simulation run (of order ). With this, the magnetic field strength follows to

(8)

and the physical time scale becomes

(9)

Unless stated otherwise, we adopt a black hole mass of and a total jet power of . In order to calculate the observable radiation fluxes, we assume a photometric distance of . The angular scale of the Schwarzschild radius then becomes . For the case of M87’s supermassive black hole with and we would have yielding an increase in resolution by a factor of compared to our fiducial scaling.

3. Jet dynamics: Acceleration and Collimation

The large intrinsic scales of the relativistic MHD jet acceleration process require substantial numerical effort when simulated with a dynamical code. Codes optimized for such tasks were developed by Komissarov et al. (2007) or Tchekhovskoy et al. (2008), involving particular grid-extension techniques which can speed up the simulations of a causally de-coupled flow. In both seminal papers the flow acceleration could be followed substantially beyond the equipartition regime to establish tight links to analytical calculations. However, it can be argued that by using a rigid, reflecting boundary of certain shape as done by Komissarov et al. (2007), or a force-free approach as applied by Tchekhovskoy et al. (2008), the rate of jet collimation resulting from those simulations could be altered.

Here we aim at studying MHD self-collimation including inertial forces. We therefore solve the full MHD equations omitting the outer fixed funnel around the jet and replace it with a stratified atmosphere which may dynamically evolve due to the interaction with the outflow. By placing the outer boundaries out of causal contact with the solution of interest, we can be certain to observe the intricate balance between jet self-collimation and acceleration that is inaccessible otherwise. We like to emphasize that the density has to be considered in the flow equations for two reasons - one is to take into account the inertial forces which are important for collimation and de-collimation, the other is our aim to consistently treat the Faraday rotation which is given by cold electrons. The latter could in principle be taken into account in magnetodynamic simulations by introducing electron tracer particles as additional degree of freedom, but is not immediately satisfied by applying the force-free limit of ultrarelativistic MHD.

We summarize our parameter runs of different jet models in table 1, indicating the maximum Lorentz factor attained, and other dynamical quantities of interest on which our results discussed in the following are based.

run ID
1h 1 0.005 8.5 25 0.21
2h 1.25 0.005 7.9 24 0.16
3h 1.5 0.005 7.9 23 0.17
1m 1 0.01 5.9 13 0.14
2m 1.25 0.01 5.6 13 0.16
3m 1.5 0.01 5.9 13 0.24

Note. – Columns denote: Simulation ID; Radial power law slope of the toroidal field; Plasma beta at (r,z)=(1,0); Maximal Lorentz factor obtained in the Raycasting domain; Maximal attainable Lorentz factor assuming complete conversion into kinetic energy . Collimation angle of the field line rooted at the inner disk radius evaluated at z=1000.

Table 1Simulation runs

Although the different electric current distributions applied in the inflow boundary condition promote a distinct jet dynamics as seen for example in the position of the light cylinder shown in figures 1 to 3, the geometry of the field lines turns out to be quite similar. At a height (corresponsing to ), the fast jet component is collimated into an opening angle less than in all our models. More obvious differences are found in the radial distribution of the Lorentz factor which is peaked at the maximum of vertical current density . In the case of closed-current models, the fast jet component becomes narrower as the integral electric current levels off more steeply. We find an acceleration efficiency in terms of the total energy per rest mass energy, of for the axial spine and acceleration efficiencies varying between for the outer parts of the jet (see the outer field lines in the middle panels). Since the flow has not reached equipartition within the considered domain, the acceleration efficiencies we obtain represent only a lower limit to the total efficiency .

Figure 1.— Jet model 1h allowing no outgoing current in the disk. Left: Current lines (dashed) shown on Lorentz factor color-contours in the (r,z) plane (note the extreme aspect ratio). Field lines are given in solid white and the light-cylinder is indicated by the solid black line. Center: Cuts through for Lorentz-factor , integral current , lab-frame density , Poynting-to-kinetic energy flux ratio (solid) and total normalized energy (dashed), acceleration efficiency , field line collimation angle and the pitch angles of the co-moving system (dashed), respectively the lab-frame (solid). Right: Acceleration along selected field-lines against the cylindrical radius r showing thermal acceleration for a field line in the spine (footpoint , above) and magnetic acceleration in the jet (, below). Vertical lines indicate the crossing of the Alfvèn (A) and fast (F) critical point as well as the light cylinder (lc).

Energy conversion is depicted in the right panels of figures 1 to 3 showing the individual energy channels normalized to the conserved rest mass energy flux along selected field lines. The terms are defined: , , and the Poynting flux . 3 In the latter relation, we have introduced Ferraro’s iso-rotation parameter defining the “angular velocity of the field line” given by

(10)
Figure 2.— As figure 1 for model 2h.
Figure 3.— As figure 1 for model 3h.

3.1. Poynting dominated flow

We first consider the MHD acceleration of the disk component of the jet flow. Here, the bulk of the acceleration takes place in the relativistic regime beyond the light surface , and can therefore be approximated asymptotically , . For a cold wind initially dominated by Poynting flux, the total energy flux per rest mass energy flux can be expressed as

(11)

Following the asymptotic relations by Camenzind (1986) we have , and hence can be used to eliminate the toroidal field from equation 11. With we can write

(12)

where , , and are conserved quantities along the stationary streamline (for details see e.g. Porth & Fendt (2010)). Thus, the asymptotic flow acceleration depends solely on the decrease of

(13)

along the flow line by differential fanning out of the field lines.4 We show the evolution of the function along selected field lines for the intermediate model 2h in figure 4 (left panel).

Figure 4.— Acceleration along field lines shown respectively as {solid,dotted,dashed,dash-dotted} lines in model 2h. Left: The function of the expanding flux tube normalized by footpoint value against radius in terms of light cylinder radii . Alfvèn () and fast () critical point transitions are marked accordingly. Right: Total energy flux ratio (top) and Lorentz factor at the fast point compared to the expected value of for various field line footpoints (bottom). We find Michel’s scaling to be satisfied within .

In the non-asymptotic regime, increases until the surface, while it is decreasing for as expected.

The second term of (12) corresponds to the Michel magnetization parameter (Michel, 1969) 5 For a critical solution in a monopole field geometry where the fast magnetosonic velocity is reached at infinity (), the terminal Lorentz-factor becomes

(14)

Different derivations of this fundamental result are given by Camenzind (1986); Tchekhovskoy et al. (2009). For small perturbations from the monopole field geometry Begelman & Li (1994); Beskin et al. (1998) could show that can be crossed at a finite distance, where again (14) is satisfied. This general scaling was also found by Fendt & Camenzind (1996); Fendt & Greiner (2001); Fendt & Ouyed (2004) for collimating relativistic jets. Our jet solutions quickly accelerate to the fast magnetosonic point and, despite the departure from the monopolar shape, follow Michel’s scaling there remarkably well. As illustrated in figure 4 (right panel), the deviation from the expectation of is less than .

Figure 5.— Characterization of the acceleration in model 2h. Left: Comparison of the field strengths for along the flux tube. We find for the most part of the domain yielding the linear acceleration regime. Also shown are power-law fits to the super-fast regime (thin solid lines). Right: The quantity along the same field lines of figure 4. Efficient acceleration in the power-law regime would yield . Alfvèn () and fast () critical point transitions are marked accordingly.

Insight into the ongoing acceleration process can be gained by an analysis of the trans-field force equilibrium as performed for example by Chiueh et al. (1991); Vlahakis (2004). The asymptotic relativistic force balance can conveniently be decomposed into “curvature”, “electromagnetic” and “centrifugal” contributions. Depending on the dominating terms, at least two regimes are possible (see also the discussion by Komissarov et al., 2009): When the curvature term is negligible, the equilibrium is maintained by balancing of the centrifugal force with the electromagnetic contribution. This constitutes the first or linear acceleration regime. The transition to the second regime occurs when field line tension begins to dominate over the centrifugal force, maintaining the equilibrium between purely electro-magnetic forces. The occurrence of curvature in the force equilibrium leads to a tight correlation between collimation and acceleration since the tension force also becomes the governing accelerating force6.

As far as a stationary state is reached, we find that the acceleration is well described by the linear acceleration regime , or

(15)

as suggested for the initial acceleration of rotating flows by various authors (e.g. Beskin et al. (1998), Narayan et al. (2007), Tchekhovskoy et al. (2008) and Komissarov et al. (2009)). This corresponds to

(16)

which is satisfied for the most part of the flow in our simulation domain. Figure 5 (left panel) shows and for a sample field line in the fast jet. We find that the critical field strengths are fairly well approximated by power-laws in the asymptotic super fast-magnetosonic regime. For the particular case shown, we have and such that the flow experiences a transition to the second, or power-law acceleration regime where the inverse of relation 16 becomes true.

For the power-law regime, a correlation between Lorentz factor and half-opening angle of the jet ,

(17)

was discovered by Komissarov et al. (2009) in the context of ultra-relativistic gamma-ray bursts. Figure 5 (right panel) illustrates the run of along a set of field lines in our fiducial model. Compared to the suggestion of equation 17, our simulation setup shows efficient MHD self-collimation, but appears less efficient in terms of acceleration.

We note that only when a substantial part of flow acceleration takes place in the power-law regime, relation 17 will hold. Our AGN jet models are however collimated to and accelerated with efficiencies of () already in the linear regime. Even if the flow acceleration is followed indefinitely, can not be recovered as this would require terminal Lorentz factors of and thus violate energy conservation.

It could be argued that the low acceleration efficiency is due to the loss of causal connection for the relativistic flow. In this case, the bunching of field-line can not be communicated across the jet anymore, thus stalling the acceleration process. This should in fact occur when the fast Mach-cone half opening-angle does not comprise the jet axis, hence (see also Zakamska et al., 2008; Komissarov et al., 2009). We have checked this conjecture by comparing both angles and found our still moderate Lorentz factor, highly collimated, jet models to be in causal connection throughout the whole acceleration domain.

3.2. Thermal spine acceleration

In this work, the very inner jet spine is modeled as a thermal wind. An alternative approach would be to prescribe a Poynting dominated flow originating in the Blandford & Znajek (1977) process. In this case, the toroidal field would be generated by the frame dragging in the black hole ergosphere below our computational domain similar to the induction in the disk. However, our attempts to increase the central magnetization by further decreasing the coronal density failed at the inability of the numerical scheme to handle the steep density gradients emerging at the boundary. Due to the vanishing toroidal field at the axis, also the Blandford & Znajek (1977) mechanism is not able to provide acceleration of the axial region (see also McKinney, 2006).

In principle, it would be possible to convert the thermal enthalpy first into Poynting flux when the jet is expanding, and then back into kinetic energy via the Lorentz force as observed by Komissarov et al. (2009). However, as we see in figures 1 to 3 (right top panels), this does in fact not occur in our simulations since the Poynting flux is approximately conserved along the inner flux lines that show little expansion. It is the magnetic field distribution and the collimated structure of the outer (disk-jet) component which merely provide the shape of the trans-sonic nozzle for the thermal wind. We can thus understand the acceleration in the jet spine by using the relativistic Bernoulli equation, which we cast in the form

(18)

An order of magnitude estimate sufficiently far from the compact object yields . Using mass conservation and a polytropic equation of state with the enthalpy , we obtain a scaling relation . For a relativistic polytropic index of this results in .

Applying a non-relativistic index of , the latter relation would yield , however, the non-relativistic limit also implies , and thus and the acceleration ceases.

Our simulations are performed employing the causal equation of state (obeying the Taub (1948) inequality) introduced by Mignone et al. (2005). Thus we obtain a variable effective polytropic index

(19)

between and . Figure 6 (top panel) shows the effective polytropic index along a stream line / flux surface.

Figure 6.— Thermal acceleration along the field line in model 2h. Top: Effective polytropic index along the flow. Bottom: Individual terms of equation 18 showing thermal energy conversion and the conservation of . In this plot, we normalized to c=1. Alfvèn () and fast () critical point transitions are marked accordingly.

In the sample stream line we find to vary between as the plasma adiabatically cools from relativistic to non-relativistic temperatures. Thermal acceleration saturates for as the enthalpy approaches the specific rest mass energy (see also Fig. 6, bottom panel). The maximum attainable Lorentz factor is given by the footpoint values at the sonic point to and depends on the detailed modeling of the inner corona. In our approach the jet spine Lorentz factor is thus limited to values of .

4. Synchrotron radiation and Faraday rotation

The numerical MHD simulations discussed above provide an intrinsic dynamical model for the parsec-scale AGN core. In the following we will use this information - kinematics, magnetic field distribution, plasma density and temperature - to calculate consistent synchrotron emission maps. What is still missing for a fully self-consistent approach is the acceleration model for the highly relativistic particles which actually produce the synchrotron radiation. However, we have compared a few acceleration models and discuss differences in the ideal resolution synchrotron maps (see below).

Radio observations of nearby AGN-cores show optically thick and thin emission features with a high degree of Faraday rotation (e.g. Zavala & Taylor, 2002). The nature of the Faraday rotation could either be internal, thus directly produced in the emitting volume, or due to an external Faraday sheet, possibly comprised of a magnetized disk-wind as ventured e.g. by Broderick & McKinney (2010), or an ambient jet cocoon. On these scales, even with global VLBI experiments, the radio emission is barely resolved for most of the known sources. In order to confront the existing observations, we perform linearly polarized synchrotron radiation transport in the relativistically moving gas, taking into account self-absorption and internal Faraday rotation. We apply beam averaging to examine the resolution dependence of the results. An illustration of our ray-tracing procedure with a rendering of an exemplary MHD solution of a collimating jet is shown in Fig. 7.

Figure 7.— Illustration of the ray-casting geometry on an exemplary solution. Color-coding in slices ( and ) and field-lines represents the bulk Lorentz-factor. An inclination angle corresponds to looking directly into the jet.

For a grid of lines of sight, each corresponding to one pixel in the final image, we solve for the parameters of linear polarization as defined e.g. by Pacholczyk (1970). This treatment provides the equivalent information as the Stokes parameters . Within the aforementioned notation, the transport equation is a linear system of equations

(20)

where denotes the emissivity vector and the opacity matrix, taking into account relativistic beaming, boosting and swing of the polarization as defined in Appendix A. Faraday rotation of the relativistically moving plasma has first been considered by Broderick & Loeb (2009b) and is directly incorporated into the previous relation via the observer system Faraday rotation angle

(21)

in cgs units, where denote the electron charge, mass and number density, the observed photon frequency and direction, the Doppler factor and the co-moving field. The dimensionless function takes into account that for high plasma temperatures the natural wave modes do not remain circular (Melrose, 1997), suppressing Faraday rotation in favor of conversion between linear and circular polarization. We follow Huang et al. (2008) and Shcherbakov (2008) in defining

(22)

in terms of the thermal electron Lorentz-factor to interpolate between the cold and relativistic limits. Especially for the hot axial flow, Faraday rotation is thus substantially suppressed. Assuming an electron-proton plasma, the electron number density follows from the mass density of the simulations. We assume further that a small subset of these “thermal” electrons is accelerated to a power-law distribution and thus responsible for the non-thermal emission of synchrotron radiation. The modeling of particle acceleration is detailed further in section 4.1.

A fraction of the Faraday rotation thus takes place already in the emitting region of the relativistic jet, such that the radiation undergoes depolarization due to internal Faraday rotation. In this case, the angular difference between the observed polarization angle and the () case can depart from the integral

(23)

which is customarily used in the diagnostics of jet observations. For example, in a uniform optically thin medium with internal Faraday rotation, the value of is just half of relation 23. Non-uniform optically thin media will break the -law and exhibit depolarization once exceeds (e.g. Burn, 1966). For optically thin media with , -law rotation measures can be recovered also in the non-uniform case, however the observed rotation angle is always less than .

We show the effect of internal Faraday rotation along an individual ray compared to the case with no Faraday rotation in Fig.8.

Figure 8.— Raytracing for an individual line of sight. Thick blue lines including Faraday rotation compared to a case where the latter was neglected - illustrated by thin red lines. The upper panels show intensity () and the Stokes parameters (dotted), (dashed) and (second panel). Faraday depth , polarization degree and the optical depth is shown in the subsequent panels. Internal Faraday rotation and the accompanying depolarization is observed in the emitting region near the plane ().

In the emitting volume the polarization degree oscillates as expected for an optically thin medium with Faraday rotation. As a consequence, the observed polarization degree is lowered. Following the density and magnetic field strength, the differential Faraday depth decreases fast enough towards the observer, so we can be confident not to miss a substantial part of the Faraday screen in the ray-casting domain.

To speed up the computation in cases of high optical or Faraday depths, we (i) limit the integration to , and (ii) solve the polarized transport only for the last internal Faraday rotations . Both optimizations do not at all affect the resulting emission maps, as the observed radiation typically originates in the photosphere of , and only a few internal Faraday rotations suffice to depolarize the radiation in the models under consideration.

Once is recovered, we obtain beam-averaged quantities via the convolution

(24)

with a Gaussian beam . The beam-averaged Stokes parameters are then used for mock observations providing spectral indices, polarization maps, rotation measure maps, and spectra to be compared to the model parameters.

4.1. Particle acceleration recipes

Within the MHD description of the jet plasma, knowledge about the relativistic particle distribution, which is needed as input for the synchrotron emission model, is not available. To recover the information from the velocity-space averaged quantities of MHD, we have to rely on further assumptions. To mention other approaches, Mimica et al. (2009) were able to follow the spectral evolution of an ensemble of relativistic particles embedded in a hydrodynamic jet simulation. Their treatment includes synchrotron losses, assuming a power-law seed distribution derived from the gas thermal pressure and density at the jet inlet.

Alas, for our purposes a consistent prescription for in-situ acceleration and cooling would be required - which seems unfeasible at the time. We therefore take a step back and assume that relativistic electrons are distributed following a global power law with index as for where signifies the overall normalization of the distribution and , denote the lower and upper cutoffs. The optically thin flux density for the synchrotron process then reads with . Optically thick regions radiate according to the source function .

This choice of particle distribution is justified by observations as well as theoretical expectations for the particle acceleration. The major physical mechanisms capable of producing non-thermal relativistic electrons are (internal) shock acceleration of relativistic seed electrons (e.g. Kirk et al., 2000) and MHD processes like magnetic reconnection (Lyubarsky, 2005) or hydromagnetic turbulence (Kulsrud & Ferrari, 1971). Considering differential rotation in relativistic jets Rieger & Mannheim (2002); Rieger & Duffy (2004); Aloy & Mimica (2008) suggested particle acceleration by shear or centrifugal effects.

In addition to and , the normalization depends highly on the mechanism under consideration. A straight-forward recipe is to connect the particle energy to the overall mass density (similar to Gracia et al., 2009)

(25)

where an ionic plasma consisting of equal amounts of protons and relativistic electrons is assumed. Thus, all available electrons are distributed following this relation and the Faraday effect is maintained by the “equivalent density of cold electrons” (e.g. Jones & Odell, 1977) in contrast to relation (21) where we assumed that the most part of electrons is non-relativistic.

An alternative to (25) is to specify the integral particle energy density and thus the first moment of the distribution function. For their leptonic jet models, Zakamska et al. (2008) have assumed that the internal energy is carried by relativistic particles, hence the relation

(26)

can be used to provide from the gas pressure resulting from the simulations.

In contrast to relativistic shock acceleration where the energy reservoir for the particles is the bulk kinetic energy of the flow, MHD processes directly tap into the co-moving magnetic energy density, and can effectively accelerate the particles up to equipartition. Accordingly, for the equipartition fraction ,

(27)

is customarily used to estimate jet magnetic field strength from the observed emission (e.g. Blandford & Königl, 1979) or vice versa (Lyutikov et al., 2005; Broderick & McKinney, 2010). To obtain peak fluxes in the Jy range, we have applied for our fiducial model. For the following discussion we have adopted (), and specified for application with relation (25). For and applying the recipes (26,27), only the cutoff energy ratio enters logarithmically into the determination of . Within these assumptions, the influence of the cutoff values on resulting jet radiation is marginal. The magnitude of then merely determines the number density of relativistic particles, to be chosen consistent with the number of particles available for acceleration. In this first study, we neglect the spectral changes introduced by the cutoff energies as this would require a more detailed modeling of the particle content which is beyond the scope of the current paper.

The observed morphology of the intensity maps is mainly given by the various prescriptions of mentioned above. In the following we briefly compare the resulting radio maps.

4.2. Radio maps for different particle acceleration models

Figure 9 shows ideal resolution maps for the aforementioned particle acceleration tracers. For the sake of comparison, Faraday rotation is neglected and with we choose a high radio frequency to penetrate through the opaque jet base.

Figure 9.— Ideal resolution logarithmic maps for model 2h at , using various tracers for the relativistic particles: Density (left), thermal pressure (center) and magnetic energy density (right). Linear () polarization vectors are overlaid as white sticks. The (x-) scale is given in terms of ray-casting footpoint and corresponds to a physical extent of .

All tracers show an almost identical polarization structure, and highlight a thin “needle” owing to the cylindrically collimated axial flow with high density, and high magnetic and thermal pressure. The axial flow is slower than the Poynting dominated disk wind which is de-beamed and, thus, not visible at this inclination. Since the axial flow features (cf. A17), the resulting polarization vector reduces to the classical case, and points in direction perpendicular to the projected vertical field of the axial spine. In the case of the density tracer, the emission becomes optically thick, as indicated by the contour. Correspondingly, the polarization degree is lowered and the direction of the spine polarization turns inside the surface. Depending on the radial density and pressure profiles , and the pressure distribution in the disk corona , the emission at the base of the jet is more or less extended and dominates the flux in all three cases. Relativistic beaming cannot overcome the energy density which is present in the disc corona, and therefore necessitates a more elaborate modeling of the accelerated particles in the jet. This will be provided in section 5.

4.3. Relativistic swing and beaming

In optically thin, non-relativistic synchrotron sources, the observed polarization vector directly corresponds to the projected magnetic field direction of the emitting region and thus carries geometric information about the jet. This allows us to interpret parallel vectors in terms of toroidal fields, while perpendicular vectors indicate a poloidal field (e.g. Laing, 1980). Similar to all realistic models of MHD jet formation, our simulations feature a helical field structure that is tightly wound within the fast jet, but increasingly poloidal further out. Hence, the resulting polarization structure is that of the telltale spine and sheath geometry - across the jet, the polarization direction flips from being perpendicular to parallel and eventually returns to a perpendicular orientation.

Due to aberration and the accompanying swing of the polarization (e.g. Blandford & Königl, 1979), an interpretation in terms of pure geometrical effects is not longer applicable in flows with relativistic velocities, instead a kinematic jet model is required. For cylindrical (i.e. -symmetric) relativistic jets, Pariev et al. (2003) have demonstrated how the optically thin polarization follows a strictly bimodal distribution, since the inclined polarization vector from the front of each annulus cancels with the corresponding polarization vector from the back side. This remains also valid for differentially rotating jets.

For the case of a collimating and accelerating jet as shown here, we loosen the constraint of cylindrical symmetry to mere axisymmetry in the direction. Additionally, our simulations feature a non-constant pitch of the magnetic field and a small degree of rotation. Together, this results in inclined polarization vectors which deviate from the strict bimodality observed in -symmetry. Figure 10 shows the optically thin polarization structure in the presence of relativistic effects (left panel), and in absence thereof (right panel). To produce the non-relativistic map, we had simply set before conducting the radiation transport.

Figure 10.— Relativistic effect on optically thin polarization (run 2h) at . The polarization degree is indicated by the background coloring, black contours show total intensity levels spaced by factors of two. Left: Including relativistic aberration. At this high inclination, the jet is de-beamed and rotation is apparent in the asymmetry of the intensity contours. Right: In the absence of relativistic effects, the polarization pattern is point symmetric about the origin and the intensity clearly promotes the jet.

At the base of the outflow where the velocities are only mildly relativistic, the polarization vectors are found predominantly perpendicular to the collimating poloidal magnetic field in both cases. Further downstream, the bimodal spine-and-sheath polarization structure prevails as the jet dynamics becomes increasingly cylindrical. The figure clearly demonstrates how the relativistic swing skews the spine towards the approaching side of the jet. Note that the jet rotation is also apparent in the beamed asymmetric intensity contours. At high inclination , the main emission from the high-speed jet is de-beamed and only the low-velocity “needle” of the thermal spine along the axis can be recognized. It is worthwhile to note that both intensity and polarization of axisymmetric non-relativistic synchrotron sources exhibit a point symmetry about the origin as illustrated in Fig. 10.

4.4. Pitch-angle dependence

Several Authors (Marscher et al., 2002; Lister & Homan, 2005) found indications for a bimodal distribution of the electric vector position angles (EVPA) of quasars and BL-Lac objects either aligned or perpendicular to the jet direction. It was also supposed that BL-Lacs tend to aligned vectors and overall higher degree of polarization. In a recent polarimetric survey however, Agudo et al. (2010) found no such correlation in their flat-radio-spectrum AGN sample, rather are their data consistent with an isotropic (mis-)alignment. Alignment is customarily attributed to oblique shocks, while perpendicular EVPAs are then interpreted in terms of a shearing of the magnetic field with the surrounding. Another possible explanation for the bimodality is due to large-scale helical fields as shown in the previous section. By varying the emitting region within the collimating jet volume, we investigate to what degree the polarization still conveys the geometric information of the emitting region. As before, Faraday rotation is neglected and we introduce the co-moving pitch angle in analogy to Lyutikov et al. (2005). In a non-rotating, cylindrical jet, the pitch-angle transforms as

(28)

and hence the co-moving fields appear much less twisted than their laboratory frame counterparts. Relation 28 together with equation 15 for the linear acceleration regime yields the co-moving fields

(29)

Substantially higher pitches are not realized within the simulations. An impression on the pitch angle distribution throughout the jet can be obtained with the cuts shown in the middle panels of figures 1 to 3. We restrict the emission to originate from , corresponding to where equipartition particle energy density (eq. 27) is assumed. This way, only the regions of the current driven jet contribute to the emission and no radiation is observed from the spurious axis where must vanish. Figure 11 shows the polarization for the two cases and various viewing angles.

Figure 11.— Polarizations for (from left to right) emitted from regions with co-moving pitches (above) and (below). The polarization degree is color-coded and contours are shown. Contours are spaced by a factor of out to where the image is cropped. Spatial scale is given in milli arcseconds and a restoring beam with FWHM= was used. The right-hand panel shows polarization angles along the cuts along core (black) and jet (gray).

In the high pitch-angle case , the resulting EVPA become parallel for viewing angles while for the standard case , this happens only at viewing angles . Most structure is observed at moderate viewing angles where we can find the core polarization perpendicular to the ridge line polarization for the case as observed in some sources (e.g. Pushkarev et al., 2005). Here we also see a spine and sheath polarization profile across the jet. In the adopted parametrization of the emission region, the spine and sheath structure is only observed for as higher pitches tend to alignment and lower pitches tend to counter align.

Beam depolarization is apparent in regions where the polarization turns and, consistent with most observations, the degree of polarization increases towards the boundary of the jet. In the high-pitch case, a left-right asymmetry is most significant with parallel vectors on the approaching side of the rotating jet and perpendicular ones at the receding one. Lyutikov et al. (2005) have proposed that based on the asymmetric polarization signal, the handedness of the magnetic field and thus the spin direction of the black hole / accretion disk can be inferred. Our results support this finding. Clearly, when increasing the pitch, the EVPAs turns from the perpendicular to the parallel direction as a general trend. However, pitches were not realized by our simulations, such that the emitting region would vanish. Thus for the models under consideration, when approaching the Blazar case, the intrinsic polarization always appears perpendicular. In the following we will mostly characterize the emission region by as this best selects the relativistic jet contribution.

5. Radiation in the jet models

We have now introduced all parameters to build models for the synthetic observations. A radiation model comprises a MHD simulation run, a particular physical scaling given by and the specific parameters of the radiation transport and . To normalize the flux, a photometric distance is assumed in all models. Table 2 summarizes the parameters adopted in this work. We will mostly report results for the fiducial model A and consult the other models only for comparison with the standard case.

model ID run ID
A 2h 0.1 1 20
B 2h 0.001 0.5 20
C 1h 0.1 1 20
D 3h 0.1 1 20

Note. – Columns denote: Model ID; Simulation ID; Black hole mass; Total energy flux rate; Minimal co-moving pitch angle of emission; Standard inclination; Resulting jet mass loss rate in terms of the Eddington accretion rate (with a radiative efficiency of ).

Table 2Jet radiation models

For the interpretation of the results, an understanding of the qualitative dependence of the observables on our parametrization is helpful and therefor discussed in the following. In terms of the physical scaling and equipartition fraction, we have

(30)
(31)
(32)

The observable quantities then become

(33)
(34)

and the opacities follow to

(35)
(36)

In reverse, based on the latter relations, spectrum and Faraday rotation measurements will allow us to constrain the physical parameters.7

We deliberately chose model B to feature times higher Faraday depth compared to model A by increasing the jet energy , while maintaining a similar spectrum with the choice of (although the total flux is thus decreased by a factor of ).

5.1. Core shift

Lobanov (1998) showed how opacity effects in optically thick jet cores provide valuable information that can help to constrain the dynamical jet quantities. As the bulk of the radiation originates in the photosphere , the projected distance of the surface will result in a specific core offset. Due to the frequency dependent opacity of a synchrotron self absorbing radio source, the measured core position varies systematically with the observing frequency. In simple model jet where the magnetic field and relativistic particle density is modeled as and (Königl, 1981), this projected distance becomes

(37)

where is a combination of the parameters and the spectral index

(38)

In a conical jet, the (predominantly toroidal) magnetic field strength will follow due to flux conservation. Conservation of relativistic particles sets for the number density . It is noteworthy that in this most reasonable case, the exponent becomes independent of the spectral index . Although our jet models are collimating and thus not exactly conical, qualitatively we expect close to unity when equipartition particle energy (eq. 27) and thus is assumed. The observed core shift is illustrated in figure 12 for the fiducial model A. Fits of the exponent for all runs are shown in the right panel.

Figure 12.— Left: Half maximum intensity contours and peak positions illustrating the core shift in mas in model A for different frequencies: 344GHz (blue) to 15 GHz (maroon). The images are aligned with respect to the imaginary black hole line of sight located at the origin. Right: Fit of the core distances relative to the 344 GHz peak for models C,A,D (top to bottom). The fit function reads with free parameters and .

The best fits are slightly steeper than the conical expectation of . Within the fitting error, models C and D are however still consistent with while the deviation in model A becomes significant. It is tempting to interpret this behavior as a consequence of jet collimation which increases the shift for larger distances from the origin. As the far regions are probed by lower frequencies, the distance law of a collimating jet is expected to systematically steepen towards the low frequency side. The quality of the fits for our model jets suggests that core shifts can provide a robust diagnostic of the AGN jet acceleration region. Unfortunately, core shift largely complicates the interpretation of the Faraday rotation maps, as we will show in section 5.3.

5.2. Depolarization

To understand the polarization signal of the simulated jets, let us briefly review some considerations on the polarization degree in general (e.g. Pacholczyk, 1970), in the presence of Faraday rotation (e.g. Burn, 1966) and when an observing beam is used.

In the simplified case of a uniform magnetic field and for power-law electron distributions, the polarization degree reads

(39)

for the optically thin regime and

(40)

for optically thick radiation. Such high polarization degrees are however rarely observed in AGN cores which go down to the percentage level. Clearly a mechanism for depolarization is required.

Already Burn (1966) considered the admixture of a isotropic random field component to an otherwise ordered field in the emitting region. He found the simple relation for depolarization relating the energies of the fields where is now the reduced polarization due to the additional random component. In result, to obtain significant depolarization the random component needs to be comparable to the ordered field. The accompanying dissipation of energy would notably decrease the efficiency of the jet acceleration process, increasing the scales of jet acceleration possibly beyond the parsec scale. The occurrence of turbulence and thus randomly oriented fields in the AGN core can serve as an explanation for various multifrequency observations as ventured by Marscher (2011). Under additional compression, even a completely randomized field structure can account for the high and low observed polarization degrees as demonstrated by Laing (1980). However, the complex physics of turbulence within the jet can not be incorporated into MHD simulations of relativistic jet formation at this time. 8 Our current MHD simulation models thus provide highly ordered near force-free fields around which relativistic electrons following an isotropic distribution are assumed to gyrate. We propose that the most promising site for finding such ordered (helical) fields is in fact the Poynting dominated regime of the jet acceleration region that is simulated here.

Under the influence of Faraday rotation within the emitting volume, the polarization degree will also depend on the Faraday opacity .
Considering first optically thick radiation where the polarization degree is governed by the ratio of ordinary to Faraday opacity . Specifically, with and it becomes . Once the photon mean free path is smaller than the correlation length of the field and the mean rotation length in the low frequency limit (), we expect to approach the uniform magnetic field case . Accordingly, for a small value of , the radiation is depolarized .

For optically thin radiation, the impact of Faraday rotation can be parametrized by describing the angular change during the emission length . The polarization degree will then decrease and oscillate according to . This is known as differential Faraday rotation and its influence on is shown in figure 8 along a line of sight in the simulation. Since , internal Faraday rotation has vanishing influence also in the high frequency limit and approaches the maximal polarization degree . In a non-uniform field, changes in the emitting geometry lower the maximal polarization degree and can only be assumed at the edges of the emission region. The direction of preferred emission perpendicular to the projected magnetic field is also the direction of dominant absorption, such that in the uniform field case the optically thick polarization direction is flipped by with respect to the optically thin polarization.

Additional depolarizing effects occur when an extended source is observed with finite resolution. The radiation is depolarized when the beam encompasses (1) intrinsic changes in the emitting geometry, (2) optical depth transitions leading to flips and (3) varying Faraday depths, known as (internal or external) Faraday dispersion (see also Gardner & Whiteoak, 1966; Sokoloff et al., 1998).

For completeness, we should also mention depolarization via “blending” - contamination with an unpolarized (thermal) component - for example radiation from the torus in the infrared. This occurs when the intensity of the contaminant becomes a notable fraction of the total intensity and thus imprints on the spectral energy distribution as well. Naturally, our results are only valid as long as the jet-synchrotron radiation is dominating the total flux.

Low Faraday rotation case

Spectrum and core polarization degree for model A is shown in figure 13. The observable core polarization degree is defined as

(41)

where the quantities under the integral are themselves subject to beam convolution.

Figure 13.— Spectrum, core polarization degree and core polarization direction in model A showing depolarization due to beam- and Faraday effects. The curves in the middle panel are shown for increasing beam-FWHM in mas as indicated and converge towards the unresolved case. The unresolved polarization direction is shown in the lower panel. Thick blue lines take Faraday rotation into account and thin red lines are calculated in its absence.

In absence of the latter, the averaged polarization degree increases monotonically from the expected value of in the low frequency range to in the optically thin case. Due to the varying emitting geometry, the theoretical maximum of is not realized. The influence of differential Faraday rotation seen in the ideal resolution case is in fact small. Depolarization occurs through ordinary beam depolarization and through Faraday dispersion, once the angle between the (unresolved) emitted polarization direction and the Faraday rotated polarization becomes larger than (compare with lower panel of Fig. 13). We note that also for the unrotated polarization direction, a flip of between the thick and thin case is not observed. This is also expected, since the photosphere probes various pitch angles as the frequency is decreased and so the simple uniform field case does not apply.

With increasing beam size, the polarization degree ultimately converges to the unresolved case. The convergence is faster for the optically thin regime where the intrinsic emission is less extended and thus quickly masked by the beam.

High Faraday rotation case

To observe the effect of depolarization due to internal Faraday rotation, we perform the same analysis as before for the high Faraday rotation model B. Figure 14 shows the resulting quantities.

Figure 14.— As figure 13 but for model B with high Faraday rotation.

As anticipated, a similar spectrum is obtained but the polarization degree and direction behave differently. Even when observed with ideal resolution, the polarization degrees of the two cases (Faraday active vs. neglected Faraday rotation) separate clearly as a result of internal depolarization. Also the unresolved polarization angles separate at higher frequencies. Due to multiple rotations, the polarization angle appears to fluctuate below observing frequencies of . In principle, the ideal resolution polarization degree is expected to rise again for lower frequencies as . However, this did not yet occur at frequencies above that were under investigation.

5.3. Rotation measure

A helical field geometry is generally perceived to promote transversal rotation measure (RM) gradients owing to its toroidal field component. First evidence for RM gradients was found by Asada et al. (2002) and Zavala & Taylor (2005) in the jet of . In several unresolved radio cores, law RMs have been detected and are found to follow a monotonic profile for example by Gabuzda et al. (2004); O’Sullivan & Gabuzda (2009); Croke et al. (2010).
Taylor & Zavala (2010) point out the observational requirements of a RM gradient detection as follows: 1. At least three resolution elements across the jet. 2. A change in the RM by at least three times the typical error. 3. An optically thin synchrotron spectrum at the location of the gradient. 4. A monotonically smooth change in the RM from side to side (within the errors).

In the following paragraphs we will touch up on each of the aforementioned points. Due the flip between the optically thick and thin polarization direction, the measurement of -law RM around the spectral peak require extra caution. Using sufficiently small spacings in to recover rotations, we can fit the rotation measure law

(42)

to the optically thick and thin cases, where now denotes the effective angle of emission. The fits are shown for a particular line of sight in figure 15.

Figure 15.— Determination of the rotation measure for an ideal resolution line of sight. The value of RM (indicated in ) necessitates multiple rotations by for the high case. In the inlay, a different fit for the optically thin regime is shown. Between the two cases, the effective angle of emission (indicated in ) is rotated by .

Here, differs by almost , whereas RM is of comparable size. However, the latter two findings are not necessarily true for all lines of sight, since the optically thin photons can originate in higher Faraday depths and different emitting geometry, as mentioned previously. We stress that with ideal resolution, consistent -laws are found both for optically thick and thin photons. Alas, if taken together, we would not be able to fit a linear function for the whole range.

The two-dimensional RM maps shown in figure 16 demonstrate the effect of beam averaging on the low- and high frequency regime. In the high frequency case, the intensity is strongly peaked close to the central object where the Faraday depth is highest, leading to steep radial gradients also in the rotation measure. As the emission is more extended, the RM is lower and smoother in the low frequency case. We find an interesting relation between the core shift and the rotation measure:

Figure 16.— Rotation measure maps () in the optically thick (top) and in the optically thin regime (bottom). White sticks indicate the direction of effective emission, . In the top plot inlay, only the region is fitted while the additional points illustrate the -law breaking due to the core shift (see text).

As the photosphere moves outwards with decreasing frequency, the flux tends to originate further away from the central object leaving systematically less Faraday active material between the observer and the source. 9 After beam convolution, the polarization angle is biased towards the outward moving photosphere, resulting in a shallower rotation measure for lower frequencies which ultimately breaks the -law that is valid for each ideal resolution line of sight, (e.g. figure 15).

Once the core shift distance is large compared to the scale of typical changes in the Faraday depth or emission angle, a -law will not be observed when the observations are aligned according to the core position. Even absolute positioning as done here10 can successfully reestablish the -law only when the beam size is smaller than the core shift. In practice, unresolved core shifts could well be the origin of optically thick -law breakers.

5.4. Resolution, mm-VLBI

In order to observe the helical fields of the jet acceleration region via the associated rotation measure, beam sizes able to resolve the dynamics are essential. In addition, the frequency must be sufficiently high to peer through the self-absorption barrier. With the advance of global mm-VLBI experiments, substantial progress will be made on both of these fronts. The theoretical resolution of a mm-observatory (at ) evaluates to , similar to the resolution of space-VLBI at . Corresponding to a physical scale of , our reference object is thus entering the regime of interest for RM studies.

At present, the record holders in terms of physical and angular resolution are the observations of the radio source Sgr A* near the galactic center black hole that were reported by Doeleman (2008) and Fish et al. (2011). Coherent structures on scales less than or could already be detected11 and sub-mm rotation measure magnitudes in excess of were discovered and confirmed by Macquart et al. (2006) and Marrone et al. (2007).

For the prominent case of M87, Broderick & Loeb (2009a) elaborate that the inner disk and black hole silhouette at could be observed with (sub) mm-VLBI. In M87, high frequency radio observations are already pushing towards the horizon scale (e.g. Krichbaum et al., 2006; Ly et al., 2007), only the important core polarization signal is still inconclusive (e.g. Walker et al., 2008). Rotation measures at core distances of vary between and depending on the location in the jet (e.g. Zavala & Taylor, 2002). It is tempting to extrapolate these values to the scale, assuming the observed RM values are non-local enhancements due to cold electron over densities. Following this argument that was put forward by Broderick & Loeb (2009a), the resulting core RM’s would be on the order of , enough to account for the low observed polarization degree via Faraday depolarization.

In the mm wavelength range, large rotation measures are needed to produce detectable deviations in the polarization angle due to the decreasing coverage of space. Typical calibration errors of require for a detection between and GHz and when GHz observations are added. The steep spectrum synchrotron flux of the jet rapidly declines when higher frequencies are considered, however, opacity and Faraday depth also decrease to yield a higher contribution of polarized flux from the core (e.g. Figure 13). Also the deviations introduced by the core-shift as described in the previous section will pose lesser problems in (sub-) mm observations.

sub-mm Rotation measure maps

To obtain detectable deflections of the polarization angle also at short wavelengths, we now focus on the high Faraday depth model B. As discussed in section 5.2.2, at frequencies beyond , the unresolved polarization vector exhibits changes below and internal depolarization is not observed (see also Figure 14), such that we find consistent -law rotation measures in the mm wavelength range. For an increasing beam FWHM from to , sub-mm RM maps for the optically thin radiation are shown in figure 17. The corresponding physical resolution is .

RM []

Figure 17.— RM maps in the optically thin wavelength range () of model B when observed with decreasing resolution from to doubling the beam size for each image. The effective emission angle is indicated by white sticks.

Steep gradients of RM across the jet axis and “spine and sheath” polarization structures are observable down to a resolution of , below which most information is destroyed by the beam. With increasing beam size, only the central Faraday pit remains detectable. Here the core rotation measure reaches values as high as .

We show the transversal cuts along the core and jet in figure 18.

Figure 18.— RM cuts along the paths indicated in figure 17 for various beam sizes. The physical resolution of the cuts is as indicated and detection limits for and observations are shaded grey in the top panel. Curves are cropped where the intensity falls below as in figure 17.

At a beam size of and observing frequencies between the transversal RM gradients of the jet fall below the assumed detection limit of and become consistent with a constant. Interestingly, beam convolution decreases the magnitude of RM not only in the cuts exhibiting a sign change, but also for the cuts across the intensity peak. We note that the intrinsic jet RM profiles are non-monotonic. Only when under-resolved, the RM features the monotonic profiles “from side to side” that are typically observed.

5.5. Viewing angle

To investigate the viewing angle dependence on the observations we show radio observables for different inclinations from to in figure 19. With a resolution of in the 43-86 GHz frequency range these images preview the next generation space VLBI experiments. We observe asymmetric features in the spectral index as proposed by Clausen-Brown et al. (2011). Due to the core shift, the maximum of the spectral index appears “behind” the low frequency intensity peak. When seen right down the jet, the polarization vectors are radially symmetric and show an inclination about the radial direction. With increasing viewing angle, the predominating polarization direction with respect to the jet flips from perpendicular to parallel orientation. Also the transversal polarization structure exhibits asymmetries as mentioned in section 4.4. To produce rotation measure maps, we fitted the law to observations at and GHz. In order to avoid the problems due to optical depth effects (see section 5.3) we exclude regions of spectral index between in the maps. The strong feature in the RM map at is most likely an artifact from the finite ray-casting domain and corresponds to the region where the axis pierces through the “lid” of the domain. For this region is excluded from the RM map and for the problem does not arise. At high inclinations, we observe steep RM gradients that coincide with the spine-sheath flip of polarization.

Figure 19.— Mock observations based on model A at and for viewing angles with resolution of . Contours according to intensity and spaced by factors of two. From left to right: Spectral index , polarization degree and rotation measure in optically thin regions.

As such high resolution data are not readily available, we also quantify the integral values that should reflect unresolved core properties. The viewing angle dependence of flux, spectral index, unresolved polarization degree, polarization direction and rotation measure are shown in figure 20.

Figure 20.— The dependence on viewing angle of unresolved quantities for 86 GHz (red) and 43 GHz (blue). Top to bottom: Beamed flux, spectral index, integral polarization degree, observed polarization vector and apparent rotation measure derived from .

The flux peaks at viewing angles of roughly and not when looking directly down the jet as would be the case for simplified cylindrical flows with toroidal fields. This reflects the fact that our jet solutions are not perfectly collimated and exhibit a small degree of rotation such that the Doppler factor attains its maximum value of at . The dominating polarization direction flips from perpendicular () to parallel () for viewing angles . At this flip, the polarization degree shows a local minimum through beam depolarization. In the framework of the unified model proposed by Urry & Padovani (1995), this suggests that the core polarization direction in the radio galaxy case (viewed at high inclination) should be clearly distinct from the blazar case (looking down the jet). For , the polarization degree approaches zero and its direction fluctuates between and due to the axial symmetry.

5.6. Towards modeling actual observations

Individual sources are modeled by constraining all free parameters (tables 1 and 2) with observations of spectra, core shift, polarization and rotation measures. Taken together, we obtain seven interrelated parameters in our radiation models.
Two of these are related to jet dynamics, parametrizing the energy partitioning with and the current distribution with the parameter . We find that the current distribution has little influence on the observables and hence this dependence can eventually be dropped, which reduces the number of free parameters but at the same time looses predictive power.
Two parameters are related to the physical scaling of the simulations, namely the total energy flux and the black hole mass . With the scale free nature of the underlying MHD, a single dynamical simulation can be used to construct a multitude of objects for arbitrary total energy flux and black hole mass. Thus in principle a library of physical jet models can be constructed from a few dynamically distinct simulations (see also: Physical scaling, section 2.4).
Compared to the simulations, the raytracing is fast and can be used to vary the remaining parameters and .

An alternative to modeling individual sources is to compare the statistical properties of a set of models with a sample of AGN cores as observed e.g. by Agudo et al. (2010).

Such an effort is beyond the scope of the current paper and we leave this open for future investigation. However, already with the acquired data (e.g. figure 20), our simulations strongly suggest a bimodal distribution of the polarization angle.

6. Summary and Conclusions

We have performed axisymmetric, special relativistic MHD simulations of jet acceleration and collimation. The resulting dynamical variables are applied to calculate polarized synchrotron radiation transport in postprocessing, providing emission maps consistent with the jet dynamical structure.

6.1. Jet acceleration

Our jets are realistically modeled to consist of two components: an inner thermal spine assumed to originate in the a black hole corona, and a surrounding self-collimating disk jet driven by Poynting flux.

We follow the flow acceleration for more than 3000 Schwarzschild radii reaching Lorentz factors in the disk jet of within the AGN “blazar zone” where we calculate the synchrotron emission maps. In application to a black hole this translates to a distance of .

Although the Poynting dominated jet flow becomes super fast-magnetosonic within the domain, it has not yet reached equipartition between Poynting and kinetic energy - jet acceleration is still ongoing. According to the available energy budget, in the case of high energy disk jets terminal Lorentz factors of would be aquired asymptotically. At the fast magnetosonic point we find the Michel scaling to be satisfied within . We find that the jet acceleration up to a distance is well described by the linear relation as proposed by Tchekhovskoy et al. (2008) and Komissarov et al. (2009). We do however not reproduce the tight coupling of acceleration and collimation observed in the latter communications but instead find to monotonically decrease along the flow. The fast jet component in all models considered collimates to half-opening angles of . We find the causal connection within the flow - its ability to communicate with the axis via fast magneto-sonic waves - to be well maintained in our simulations.

Also the thermal spine acceleration is shown to be efficient with and limited only by the amount of enthalpy available at the sonic point, in our case to .

We have placed the location of the outflow boundaries out of causal contact with the propagating jet beam of interest. Thus, we can be sure that the calculated jet structure is purely self-collimated, and does not suffer from spurious boundary effects leading to an artificial collimation. We have investigated jets with a variety of poloidal electric current distributions. We find - somewhat surprisingly - that the topology of the current distribution, e.g. closed current circuits in comparison to current-carrying models, has little influence on the jet collimation.

6.2. Jet radiation

The benefit of having performed relativistic MHD simulations of jet formation is that we could apply them to produce dynamically consistent emission maps to predict VLBI radio and (sub-) mm observations of nearby AGN cores. For this purpose, we have developed a special relativistic synchrotron transport code fully taking into account self-absorption and internal Faraday rotation. Since the acceleration of non-thermal particles can not be followed self-consistently within the framework of pure MHD, it remains necessary to resort the particle energy distribution to simple recipes. We have compared three prescriptions of the non-thermal particle energy distribution. We found good agreement in the alignment of the polarization structure, but considerable differences in the intensity maps. Thus, the polarization maps derived in this work can be considered as robust, while the intensities distribution should be regarded with caution.

The strict bi-modality of the polarization direction suggested by Pariev et al. (2003); Lyutikov et al. (2005) and others can be circumvented when the structure of a collimating jet is considered. However, the efficient collimation to near-cylindrical jet flows in general confirms these results obtained for optically thin cylindrical flows when the fast jet is considered. Thus, depending on the pitch angles of the emission region, also a spine-and-sheath polarization structure could be observed. The relativistic swing effect skews the polarization compared to the non-relativistic case. Our radiation models affirm the finding of Lyutikov et al. (2005) and Clausen-Brown et al. (2011) that relativistic aberration promotes asymmetries in the polarization (half spine-sheaths) and also in the spectral index. The observational detection of such features would allow to determine the spin direction of the jet driver, be it the accretion disk or the central black hole.

The frequency-dependent core shift in the radiation maps following our jet simulations is consistent with analytical estimates of conical jets by Lobanov (1998) in two jet models and slightly steeper in the third case considered. We attribute this discrepancy to fact of jet collimation. The overall good agreement with the analytical estimate suggests that the standard diagnostics should provide robust results capable of determining the jet parameters. With our radiation models we have confirmed the intuition that unresolved core shifts should lead to a breaking of the rotation measure law. Further, we have demonstrated that that law can be restored again as soon as the resolution is increased. Opacity effects do not allow to obtain a consistent -law across the spectral peak. Once the regimes are separated however, we obtain two valid relations for which the optically thin rotation measure is substantially increased over the optically thick case as it peers deeper into the Faraday pit.

The interpretation of observations featuring both internal Faraday rotation and changes in opacity is one of the most challenging aspects in polarimetric imaging of jets. With our detailed modeling, we were able to disentangle the depolarizing effects of opacity transition, differential Faraday rotation, and also beam effects such as ordinary beam depolarization or Faraday dispersion for two exemplary jet models. We find that the unresolved, optically thin mm-wavelength radiation is depolarized due to both the changing emission geometry (down to ), and the additional beam depolarization (down to ). Increasing the Faraday opacity by observing at lower frequencies would lead to depolarization below the level due to 1. Faraday dispersion and 2. differential Faraday rotation.

We have also investigated the influence of resolution on the detectability of rotation measure (RM) gradients in the optically thin parsec-scale jet-core previewing mm-VLBI observations. To detect such gradients across the jet, a resolution of would be required. Increasing the beam size leads to more monotonic transversal RM profiles. We find the peak magnitude of the RM to increase with resolution. High RMs beyond are required to obtain a noticeable deflection in the mm-wavelength range. From the sources where high resolution data is already available, namely Sgr A* and M87, such high rotation measures are in fact observed, and we predict that many more objects in this class will be found at the advent of ALMA and global mm-VLBI.

Finally, we have presented mock observations of spectral index, polarization degree and rotation measure for various inclinations. Asymmetries in the spectral index and polarization degree can be observed most clearly at high inclinations . The necessary resolution for this detection in a fiducial low Faraday rotation case (our model A) amounts to , which could be reached with the next generation space-VLBI. At , the predominant polarization vector flips from perpendicular alignment (with respect to the projected jet direction) for the blazar case to parallel alignment for the radio galaxy case at high inclinations. The flip in polarization is clearly detectable also from unresolved quantities. In summary, these findings suggest a bimodal distribution of the observed polarization direction of AGN core jets. However, by adding a substantial amount of Faraday rotation (our model B), this signature will be scrambled unless the observing frequency is chosen high enough - confirming the popular intuition.

In this paper we have focussed on general signatures of the synchrotron radiation in the large-scale helical fields in the acceleration region of relativistic MHD jets. With the developed tool set in hand, further progress can be made when calibrating the observational diagnostics with the mock observations that are detailed here. We expect a substantial improvement from a more consistent treatment of the non-thermal particles, taking into account particle acceleration and cooling. Potentially, also modeling of individual sources, or the cumulative statistics of AGN surveys applying dynamical simulations could be undertaken in the future.

This work was partly carried out under the HPC-EUROPA2 project (project number: 228398), with the support of the European Community - Research Infrastructure Action of the FP7. Post-processing of the simulations was performed on the THEO cluster of MPIA and fits were obtained with the flexible MPFIT routines provided by Markwardt (2009). The authors thank an anonymous referee for comments and suggestions that have helped to improve the presentation of this work. O.P. likes to thank Christophe Sauty for discussions and kind hosting during an interesting research visit.

Appendix A Stokes transport

In the basis, the linear transfer equation (20) has the coefficients

(A1)
(A2)
(A3)