From electrons to Janskys: Full stokes polarized radiative transfer in 3D relativistic particle-in-cell jet simulations

The underlying plasma composition of relativistic extragalactic jets remains largely unknown. Relativistic magnetohydrodynamic (RMHD) models are able to reproduce many of the observed macroscopic features of these outflows. The nonthermal synchrotron emission detected by very long baseline interferometric (VLBI) arrays, however, is a by-product of the kinetic-scale physics occurring within the jet, physics that is not modeled directly in most RMHD codes. This paper attempts to discern the radiative differences between distinct plasma compositions within relativistic jets using small-scale 3D relativistic particle-in-cell (PIC) simulations. We generate full Stokes imaging of two PIC jet simulations, one in which the jet is composed of an electron-proton ($e^{-}$-$p^{+}$) plasma (i.e., a normal plasma jet), and the other in which the jet is composed of an electron-positron ($e^{-}$-$e^{+}$) plasma (i.e., a pair plasma jet). We examined the differences in the morphology and intensity of the linear polarization (LP) and circular polarization (CP) emanating from these two jet simulations. We find that the fractional level of CP emanating from the $e^{-}$-$p^{+}$ plasma jet is orders of magnitude larger than the level emanating from an $e^{-}$-$e^{+}$ plasma jet of a similar speed and magnetic field strength. In addition, we find that the morphology of both the linearly and circularly polarized synchrotron emission is distinct between the two jet compositions. We also demonstrate the importance of slow-light interpolation and we highlight the effect that a finite light-crossing time has on the resultant polarization when ray-tracing through relativistic plasma.


Introduction
Relativistic extragalactic jets are among the most persistent energetic objects in the universe. They are composed of collimated beams of magnetized relativistic plasma that can extend up to thousands (and in some cases millions) of parsecs from their host galaxies. The current theoretical paradigm postulates that the ultimate physical mechanism powering these relativistic outflows is the energy released from matter accreting onto spinning supermassive black holes (Blandford & Znajek 1977;Blandford & Payne 1982). The plasma content of these jets (and their central engines) remains an active area of research (see, e.g., Croston et al. 2018;Fan et al. 2018;Thum et al. 2018;Myserlis et al. 2018a,b;Enßlin et al. 2019;Anantua 2020;Sikora et al. 2020;Emami et al. 2021;Moscibrodzka et al. 2021;Ricarte, Qiu, & Narayan 2021).
With the advent of global (and recently space-based) millimeter-wave very long baseline interferometry (VLBI), we are able to probe the polarized emission emanating from the innermost regions of a number of jets. In particular, the linearly and circularly polarized synchrotron emission from these jets carry imprints of both the strength and orientation of the collimating magnetic fields as well as the plasma content of each jet. Studying the nature of this synchrotron emission can be used to infer the physical conditions both within the jet and in the surrounding environment into which the jet propagates.
In parallel to this observational advance, modern computational resources have allowed for increasingly sophisticated numerical plasma simulations. In particular, 3D particle-incell (PIC) simulations (e.g., Nishikawa et al. 2014Nishikawa et al. , 2016aAlves 2018;Guo et al. 2020;Nishikawa et al. 2020a) have enabled, for the first time, a self-consistent treatment of the kinetic effects occurring within relativistic plasma outflows, such as the following: plasma instabilities, jet shear and entrainment, and magnetic reconnection (see Birdsall &Langdon 1995 andNishikawa et al. 2020b for a summary of PIC methods). These PIC simulations, however, are numerically intensive, and this kinetic precision comes at the cost of small (relative to relativistic magnetohydrodynamic -RMHD) simulation sizes. There exists, therefore, a synergy between PIC and RMHD jet modeling that should be exploited in order to gain a more Article number, page 1 of 21 arXiv:2106.04915v1 [astro-ph.HE] 9 Jun 2021 holistic understanding of both the micro and macro physics of relativistic outflows. While RMHD simulations can effectively model the large-scale fluid motions of the jet (e.g., Mizuno et al. 2015;Tchekhovskoy et al. 2016;Fuentes et al. 2018;Fromm et al. 2019;Mukherjee et al. 2020Mukherjee et al. , 2021, PIC can be used to model the microphysics and radiative processes occurring within the jet plasma (e.g., Zhang et al. 2018;Petropoulou et al. 2019;Zhang et al. 2020;Hosking & Sironi 2020;Davelaar et al. 2020;Sironi, Rowan, & Narayan 2021). These kinetic-scale processes form a direct link to VLBI observations of the polarized synchrotron emission.
In this paper, we compare the radiative differences between two PIC jet simulations (computed using the TRISTAN-MPI code 1 ; see Niemiec et al. 2008& Buneman 1993, one in which the jet is composed of an electron-proton (e − -p + ) plasma (i.e., a normal plasma jet), and the other in which the jet is composed of an electron-positron (e − -e + ) plasma (i.e., a pair plasma jet). We make use of a polarized radiative transfer scheme (see MacDonald & Marscher 2018) that has been embedded into a ray-tracing code for post-process imaging of each numerical jet simulation. We use this ray-tracing code to create synthetic full Stokes (I, Q, U, & V) images of each numerical PIC jet simulation. This paper is organized as follows: In §2 we summarize the scaling relations used in our PIC jet simulations. In §3 & 4 we outline the radiative transfer theory adopted in our study. In §5 we present the results of our ray-tracing calculations through the e − -p + plasma and e − -e + plasma jet simulations. Finally, in §6 we present our summary and conclusions. We adopt the following cosmological parameters: H o = 71 km s −1 Mpc −1 , Ω m = 0.27, and Ω Λ = 0.73.

Scaling
PIC calculations are carried out in dimensionless grid units which must be scaled (via scaling factors) into physical units (i.e., cgs) as a post-process step for use in our ray-tracing calculations. The PIC simulations presented in this paper are small (120 × 120 × 240 cells) segments of the e − -p + and e − -e + PIC jet simulations published in Nishikawa et al. (2017).
The simulation cell size (∆l ), the time step (∆t ), and the speed of light (c) are computed in the following manner: where by definition the grid values: ∆l grid = ∆t grid = c grid = 1. The quantities: ∆l scale , ∆t scale , and v scale represent scaling factors that set the physical dimensions of each jet simulation. In both jet simulations presented in this paper, we set v scale = 3 × 10 10 cm s −1 . The plasma skin depth (≡ c/ω pe ) parameterizes the length scale across which low-frequency electromagnetic radiation can propagate within a collisionless plasma. The quantity ω pe denotes the relativistic electron plasma frequency and is given by: where n e , e, and m e denote the number density [cm −3 ], charge [cm 3/2 g 1/2 s −1 ], and mass [g] of the electrons within the 1 https://ascl.net/1908.008 plasma and γ e is the electron Lorentz factor. The plasma frequency is the fundamental frequency of oscillation of the electrons and together with the plasma skin depth sets a fundamental length scale within our PIC jet simulations. In order to model the kinetic scales within our jet simulations, the plasma skin depth of the jet plasma is resolved across ten computational grid cells by the TRISTAN code: Equation 5, therefore, defines the characteristic length scale of the numerical cells within our PIC jet simulations. We define a fiducial jet electron number density as follows: to which we scale our PIC values. By combining Equations 1, 3, 4, 5, and 6 we obtain the following PIC cell size scale factor: 4π n e e 2 6.5 × 10 4 cm , having set m e = 9.1094 × 10 −28 g, e = 4.8066 × 10 −10 cm 3/2 g 1/2 s −1 , and γ e 15 (the limiting value at the jet base in our models). This scaling implies that the jet axis of both of our PIC simulations (i.e., 240 cells in length) corresponds to a physical size of ∼ 150 km. This physical size demonstrates the vastly different spatial scales probed by PIC simulations in comparison to RMHD and represents an inherent challenge PIC codes face when attempting to model the plasma of parsec scale astrophysical jets.
As in Equations 1, 2, and 3, we scale the magnetic field within our PIC simulations in the following manner: where the dimensionless grid value at the base of our jet models is B grid 2.7 (and varies throughout each simulation box), and B scale is a magnetic field scaling factor [Gauss]. We define a fiducial electron plasma beta (β e ) at the jet base in each of our PIC simulations to which we scale the magnetic field strength: The value of β e has been chosen to be indicative of a highlymagnetized synchrotron emitting relativistic plasma. By combining Equations 3, 8, and 9 we arrive at the following magnetic field scale factor: again having set n e 10 1 cm −3 and γ e 15. The above scale factor results in magnetic field strengths of B 10 0 Gauss at the base of each jet (see Figure C.1), which in combination with our fiducial electron number density of n e 10 1 cm −3 produces numerically tractable levels of synchrotron emission at radio/mm wavelengths within our PIC simulations. This choice of physical scaling also helps to ensure that the electron Larmor radii do not vastly exceed the scaled PIC cell sizes (see Appendix A).
The dimensionless electron grid value at the base of our jet models is n e grid 64 (and varies throughout each simulation box). Combining Equations 1, 7, and 11 results in a scaled electron number density of n e 2.3 × 10 −13 cm −3 at the jet base in each of our models. This value is many orders of magnitude below the levels inferred from observational and theoretical modeling of relativistic jets (i.e., our fiducial value of n e 10 1 cm −3 ). This difference in particle number again highlights the challenge faced by PIC codes when attempting to model the plasma of parsec scale astrophysical jets.
To help mitigate the numerical limitation in our ability to simulate astrophysically plausible numbers of electrons in our PIC jet calculations, we introduce a PIC "super particle" parameter ( f p ) which we apply to each computational cell within our jet models: n e grid = f p × n e grid , where n e grid is a proxy for the number of 'real' electrons represented by our 'simulation' electrons. As in Equations 1, 2, 3, and 8, we formulate the following electron number density scaling relation which we apply throughout our simulations: n e = n e scale × n e grid , where: Since ∆l 6.5 × 10 4 cm (Equations 1 and 7) and n e grid 64, in order to scale the dimensionless grid electron number densities at the jet base to our fiducial jet value of n e 10 1 cm −3 , we need to tune our PIC 'super particle' parameter to f p 4.6×10 13 . This immense value highlights a numerical limit PIC codes face in simulating plasmas on astrophysical scales and also sets a numerical benchmark for future simulations with larger particle populations. With our dimensionless PIC grid values scaled into physical units (i.e., cgs), we are now in a position to apply polarized radiative transfer via ray-tracing, in order to infer both the level and the morphology of the polarized synchrotron emission emanating from the normal plasma jet and the pair plasma jet (both of which are illustrated in Figure 1). Jones & O'Dell 1977a,b;Jones 1988 present solutions to the full Stokes equations of polarized radiative transfer for synchrotron emission emanating from a homogeneous and an inhomogeneous magnetized plasma containing isotropic distributions of electrons (see Appendix B for further discussion). We use these solutions to compute the levels of LP and CP emanating from our PIC jet simulations at radio frequencies. In particular, we solve the matrix presented in Equation 14 along individual rays passing through our PIC simulations.

Polarized radiative transfer
Here, I ν , Q ν , U ν , and V ν denote the frequency-dependent Stokes parameters, while (κ I , κ Q , κ U , κ V ) and (η I ν , η Q ν , η U ν , η V ν ) represent, respectively, the synchrotron absorption and emission coefficients for each Stokes parameter. The terms (κ * Q , κ * U ) and (κ * V ) account for the effects of Faraday conversion and rotation (see Appendix C), respectively, within our PIC plasma. The term l denotes the path length of each ray through the computational cells of our PIC calculations. Jones & O'Dell (1977a) present an analytic solution to this matrix which we apply along each ray. The radiative transfer is carried out in the 'co-moving' frame of the plasma with: (i) a relativistic abberation correction being applied, cell-by-cell, to obtain the angle between each cell's local magnetic field vector and the observer's inclination to the jet axis , and (ii) a rotation correction being applied, cell-by-cell, to transform the linear polarization ellipse from the local co-moving frame onto the plane of the sky (see Appendix D for further discussion). Once the Stokes parameters have been generated for each sightline (i.e., for each pixel/ray in our synthetic maps) a Doppler factor is applied (which incorporates the velocity/angle dependence of the larger scale jet) to obtain the flux levels in the observer's frame. This Doppler boosting/de-boosting is apparent upon comparison of the relative flux levels between im-Article number, page 3 of 21 A&A proofs: manuscript no. paper Schematic representation of our PIC slow-light interpolation scheme. A ray (I 0 ν ) enters the jet and encounters successive upstream plasma cells as the ray propagates across the jet (similar to stepping through a fast moving stream). These successive encounters contribute to the emission stored in each ray and result in the final total intensity value (I ν ) that exits the jet and produces one pixel in our synthetic radio maps. ages generated when viewing each simulation edge-on and at right angles to the jet axis. We have embedded this polarized radiative transfer scheme into the ray-tracing code RADMC-3D 2 . For the images presented in this paper, RADMC-3D casts 640, 000 individual rays through our cartesian PIC grids forming 800 × 800 square-pixel images of the polarized emission produced by our jet models. The resultant emission maps from these ray-tracing calculations for the normal plasma and pair plasma jet simulations are presented in Figures 3, 4, and 6.

Slow-Light Interpolation
In order to resolve (both spatially and temporally) the kineticscale processes occurring within a relativistic jet plasma, the time step (∆t ) of the TRISTAN code is set such that: Therefore, it takes ten numerical time steps for a light ray to traverse the length of an individual plasma cell within each jet calculation. Equation 15 has profound implications for our ray-tracing calculations which typically assume/invoke the fast-light approximation. Fast-light ray-tracing makes the assumption that the light-crossing time of the jet is far smaller than the dynamical time step of the jet simulation being imaged (i.e., that we are taking a snapshot of static jet plasma). This assumption is not in general valid for PIC simulations. The fact that the plasma is evolving in time as each ray propagates through the computational grid must be taken into account within our ray-tracing calculations. In particular, we have constructed a slow-light interpolation scheme that accounts for this effect and is illustrated in Figure 2. In essence, our interpolation scheme builds up a hybrid computational grid composed of the stratified components of successive time steps within each jet simulation. This scheme ensures that the plane parallel rays of our ray-tracing calculations encounter the correct upstream plasma values as each ray propagates through the jet. Within our PIC jet simulations: 2 http://ascl.net/1202.015 • The time scale for light crossing is ∼ 500 simulation time steps. • The time scale for jet plasma evolution is ∼ 10 simulation time steps. • The time scale for particle acceleration is < 10 simulation time steps.
Since the time scale for light crossing than the time scale for jet plasma evolution, slow light interpolation has a noticeable effect on the resultant synchrotron emission (see Figures 3 and 6). Our radiative transfer scheme, however, is unable to incorporate the smallest particle acceleration time scales within our emission calculations, since we implicitly assume that the plasma properties of a given cell remain constant during the 10 simulation time steps it takes a light ray to traverse each cell. We also emphasize that when scaled into physical units, a simulation time step within our jet models corresponds to ∼ 2.2 −7 s. This infinitesimal value is again a reminder of the drastically smaller spatial and temporal scales probed by PIC models in comparison to RMHD jet calculations.

Results
To gain a better understanding of the effects of slow-light interpolation within our ray-traced images, we first construct fast-light images for the purposes of comparison between the two ray-tracing methods. Both the normal plasma jet and the pair plasma jet simulations were run for a sufficient time span to allow both jets to propagate across the full extent of both computational grids (shown in Figure 1). At each time step within our PIC simulations we apply the scaling relations presented in §2 to our PIC grid values in order to create ray-tracing output files. These output files consist of three dimensional arrays containing: • magnetic field strength (and orientation): where in general E = γ e m e c 2 , and the indices [i , j , k] denote the [x, y, z] directions within our cartesian computational grids, respectively. In contrast to magnetohydrodynamic calculations where prescriptions for the electron values must be applied in order to infer these quantities from the thermal fluid variables (e.g., Porth et al. 2011;Fromm et al. 2016), in our PIC calculations, we are able to compute these values directly from our jet simulations. The output files from the TRISTAN particle-in-cell code then become input files for the RADMC-3D ray-tracing code. Given the fact that the entire computational grid in each PIC calculation (when scaled into physical units) only spans roughly ∼ 150 km (as discussed in §2) we have arbitrarily placed each jet at a distance of 1 Astronomical Unit (AU) from the 'observer' (i.e., the distance to the Sun). This distance scale, in combination with our simulation box sizes, generates radio maps with an angular extent of ∼ 200 milli-arcseconds (mas). Clearly, 150 cubic kilometers of relativistic jet does not exist at such a close distance to the Earth, but, in the spirit of theoretical study, we proceed with these calculations to infer what: (i) the morphology and (ii) the fractional level of polarized synchrotron emission would be if we could image each of these plasma jets with an idealized interferometric array at such close proximity.

Fast-light Images
The results of our fast-light ray-tracing calculations for both simulations are presented in Figures 3 and 4. We explore the nature of the polarized synchrotron emission when viewing each jet at right angles to the jet axis (i.e., θ obs = 90 • , shown in Figure 3) and when viewing each jet edge-on to the jet axis (i.e., θ obs = 0 • , shown in Figure 4). For both sets of images we have performed our ray-tracing calculations through the last time step of each simulation (illustrated in Figure 1).
We point out that in the case of a 'pure' pair plasma, Stokes V would be exactly zero, since the equal numbers of electrons and positrons would cancel out their respective contributions to the circularly polarized synchrotron emission. Similar to the methods presented in Wardle & Homan 2003;Homan et al. 2009, and more recently Anantua (2020), we parametrize within each plasma cell the proton-to-position ratio (r + ≡ n p + /n e + ). This ratio is incorporated into the radiative transfer coefficients, including the CP specific terms: η V ν , κ V , and κ * V (see Appendix C; see also MacDonald & Marscher 2018). We have initially fixed this ratio to r + = 100 (i.e., 100 protons for every positron) in the normal plasma jet and r + = 0.01 (i.e., 100 positrons for every proton) in the pair plasma jet. We plan on running a larger set of PIC plasma simulations in the future in which we explore a wider range of jet plasma compositions than the two cases investigated here.

Radio Jet Orientation
In Figure 3, we present ray-traced images of each jet simulation at θ obs = 90 • for an observing frequency of ν obs = 1 GHz (see Appendix E for further discussion). In the upper panels we present Stokes I maps of both jets (normal plasma on the left and pair plasma on the right), in the middle panels we present the corresponding maps of linearly polarized intensity (P ≡ Q 2 +U 2 ) for each jet with electric vector position angles ( EVPAs; χ ≡ 1 2 arctan U Q ) overlaid in white. In the lower panels of Figure 3 we present the corresponding Stokes V images for both simulations (again normal plasma on the left and pair plasma on the right). To mimic the effects of interferomet-ric resolution, all images have been convolved with a circular Gaussian beam of FWHM 7.5 mas. Upon inspection of Figure  3, one can see that the pair plasma jet exhibits a slightly more filamentary emission morphology in contrast to the normal plasma jet. This morphological difference in emission is a reflection of the distinct jet dynamics occurring within the two plasmas: the e − -e + plasma jet is prone to larger plasma instabilities (such as the Kelvin-Helmholtz and Weibel instabilities) and is, as a result, less stable than an e − -p + plasma jet of similar speed and magnetic field strength (see, e.g., Nishikawa et al. 2016bNishikawa et al. , 2017Nishikawa et al. , 2019.

Blazar Jet Orientation
In Figure 4, we present ray-traced images of each jet simulation at θ obs = 0 • for an observing frequency of ν obs = 230 GHz (see Appendix E for further discussion). In the upper panels we present Stokes I images of both jets (normal plasma on the left and pair plasma on the right), in the middle panels we present LP images, and in the lower panels we present the corresponding CP images for both simulations. We have convolved these higher frequency images with a beam size similar to the lower frequency 1 GHz images for ease of comparison with Figure 3. The effect of Doppler beaming is apparent upon comparison of the Stokes I flux levels in Figures 3 and 4. The fractional circular polarization (m c ≡ −V /I ) is minimal in each simulation ( 1%) but is many orders of magnitude larger in the e − -p + plasma jet in comparison to the e − -e + plasma jet (integrated values are listed to the lower right in each Stokes V image). We present spectropolarimetry of the integrated levels of fractional polarization for each jet in Appendix F.
In addition to mimicking the resolution of an interferometric array (via beam convolution), our ray-tracing algorithm can also mimic the sensitivity of an interferometric array by including a synthetic Gaussian noise floor within our ray-traced images. This is discussed further in relation to the detectability of our PIC jets in Appendix G.

Individual Ray-Properties
In Figure 5, we illustrate the variations of the Stokes parameters (upper panels) and the Faraday rotation (τ F ) and conversion (τ C ) depths (lower panels) of both jet plasmas along the sightlines indicated by the red dots in the lower panels of Figure 4. The Faraday depths (which are written out explicitly in Appendix C and are themselves functions of B , n e , and γ min ) parameterize the ability of the plasma cells along each sightline to attenuate both the linearly and circularly polarized synchrotron emission within the jet. It is evident upon comparison of the lower panels of Figure 5 that, for these particular sightlines, the Faraday rotation depth is larger in the e − -p + plasma jet. Both jets, however, are optically (and Faraday) thin and the Stokes V images presented in Figure 4 are intrinsic in origin (see the lower panels of Figure F.1 in Appendix F).

Slow-light Images
The results of our slow-light interpolated ray-tracing calculations are presented in Figure 6. The slow-light interpolation results in an averaging of the emission along each sightline. In contrast to the fast-light images, the morphology of the resultant emission becomes blurred and the distinct morphologies present between the normal plasma and pair plasma jet com-  positions (evident in the fast-light images shown in Figure 3) are less pronounced in the corresponding slow-light images illustrated in Figure 6. In particular, the hybrid computational grids (constructed via the algorithm illustrated in Figure 2), through which we ray-trace, are composed of ∼ 50 distinct/stratified jet epochs (each separated in time by 10 code time steps). This number of epochs corresponds to the number of distinct plasma cells that a plane-parallel ray encounters when each jet is imaged at right angles to the jet axis. As discussed in §4, as each ray propagates across the jet, it encounters newer upstream values of jet plasma due to the finite light crossing time across each PIC plasma cell (i.e., ten time steps -see Equation 15). Further computational time (i.e., more epochs) is required in order to produce slow-light images of these simulations when each jet is viewed edge-on (i.e., with the orientation of a blazar).
We finally point out that features/blobs in our slow-light images, in contrast to the fast-light images, do not necessarily map to individual plasma structures within the jet flow and are instead a mixture of emission from multiple plasma components along various sightlines through the jet.

Summary & Conclusions
We have carried out full Stokes polarized radiative transfer calculations (via ray-tracing) through 3D relativistic PIC simulations of a normal plasma (e − -p + ) and of a pair plasma (e − -e + ) jet. We generate two sets of images of each jet simulation, one in which we invoke the fast-light approximation (see Figure 3), and the other in which we implement slow-light interpolation (see Figure 6). It is clear, upon comparison of Figures 3 and  6, that slow-light interpolation has a discernible effect on the emission emanating from within each jet. In particular: • The finite light-crossing times through our relativistic jet simulations results in 'blending' of various plasma emission features along each sightline.
It is also clear that there are differences both in the morphology and in the fractional level of polarization emanating from the two jet plasma compositions. Specifically: • The pair plasma jet exhibits a more filamentary synchrotron emission morphology in comparison with the normal plasma jet. • The normal plasma jet has a much larger value of integrated fractional circular polarization in comparison with the pair plasma jet.
These differences are a reflection of the distinct plasma dynamics occurring within each jet: the pair plasma jet is prone to larger kinetic instabilities (e.g., Kelvin Helmholtz and Weibel instabilities) in comparison with the normal plasma jet and is, as a result, less stable, resulting in a more filamentary emission structure. The lower levels of circular polarization in the pair plasma jet are in keeping with the synchrotron theory we have incorporated into our ray-tracing calculations. We emphasize that the applicability of our calculations to jets on parsec scales is quite limited. In particular: • Both jet simulations, when scaled into physical units, are only hundreds of kilometers in extent and are extremely tenuous in nature.
We present these calculations, despite this limitation in scale, for two main purposes: (i) to provide a point of comparison for future polarimetric imaging of larger-scale PIC simulations, and (ii) to emphasize that kinetic scale jet dynamics can produce distinct morphologies in the jet's synchrotron radiation. These calculations reveal how the micro physics of the jet affect the macro emission we detect in the radio. This relationship between the micro and the macro is not commonly addressed in most relativistic jet simulations. Refinement of: (i) the jet injection scheme, (ii) the grid size, and (iii) the jet particle content are planned for future analysis. This work clearly emphasizes the dire need for vastly larger computational grid sizes and particle populations when attempting to model kinetic scale effects in relativistic plasmas on astrophysical length scales. Hybrid techniques (e.g., Mignone et al. 2018;Vaidya et al. 2018;Davelaar 2019;Parfrey, Philippov, & Cerutti 2019;Bacchini et al. 2020) will be crucial for future jet simulations. We also highlight the care that must be taken when interpreting features in interferometric radio maps of relativistic jets, especially if the light crossing time of the jet exceeds the relevant dynamical timescales of the emitting plasma. radii of our simulated PIC electrons are largely contained within our numerical cells. The regions where the electron Larmor radii exceed the plasma cell size occur where the magnetic field drops precipitously (see the top panels of Figure C.1). These zones contribute minimally to the integrated synchrotron emission along each ray since the synchrotron emissivity is ∝ B 2 .

Appendix B: Electron Phase-Space & Energy Distributions
The synchrotron emission/absorption coefficients (contained in Equation 14) assume that the underlying electron distribution within the jet is isotropic. While we have attempted to scale our PIC simulations in a manner that ensures the individual electron Larmor radii are largely contained within our grid cells (see Appendix A), this scaling does not guarantee that the particle distribution is indeed isotropic in nature. To explore this point further, we have generated 2D phase-space distribution plots of the electrons contained within the jet regions we generate synchrotron emission from (see the left and middle panels of Figure B.1). Clearly, anisotropy is present within the underlying particle distribution used in our emission calculations (especially along the jet axis - Figure B.1 left panel). Despite this anisotropy, we compute the synchrotron emission using the fully isotropic emission coefficients of Equation 14 in order to gain a 'first approximation' of the synchrotron emission emanating from our PIC simulations.
The synchrotron emission is calculated in a post-process fashion. We first obtain (within each cell) the local 3D fluid velocity and magnetic field components which are averaged over 27 neighboring cells with weighting to help ensure good statistical properties within each cell. These values are then used to compute the cell's 'local' magnetic field strength and 'local' electron Lorentz factor. In addition, the number of particles are counted and scaled (via Equation 12) to obtain a 'local' electron number density.
As a further simplification in our emission calculations, we have arbitrarily fixed a constant power-law index (s = 2.3) within each plasma cell. Future refinement of our algorithm will be needed in order to connect the spectral indices used in our synchrotron emission calculations more closely to the 'local' particle energy spectrum within each cell. To explore this point further, we have generated a particle energy spectrum of the electrons contained within the jet regions we compute synchrotron emission from (see the right panel of Figure B.1). Clearly, the particle spectrum (especially the high-energy tail) is much steeper than the power-law index of s = 2.3 used in our emission calculations. This remains a limitation/discrepancy between our PIC models and our synchrotron emission calculations.
The angle ϑ denotes the angle that each sightline within our ray-tracing calculations makes with respect to the local magnetic field vector in the co-moving frame of each plasma cell. Under the assumption of charge neutrality (n e − n e + + n p + ), the term ( n e − − n e + )/( n e − + n e + ) in Equation C.2 may be re-written as ( 1 + 2/r + ) −1 , where r + ≡ n p + /n e + and parameterizes the plasma composition of each cell (see §5.1). In all of the computations presented in this paper, we have assumed a constant optically thin spectral index of α = 0.65. In future, we plan on refining our algorithm to vary α cell-to-cell and compute it directly from the local electron power-law index (s): α = (s − 1)/2. The synchrotron opacity (τ) of an individual plasma cell may be written as: where l is the path length of the ray through the plasma cell (computed using RADMC-3D) and where n e is the electron number density within the plasma cell. The length scale of the PIC plasma cells in the simulations presented in this paper is 6.5×10 4 cm (see Equation 7). The parameter κ α is a physical constant (of order unity) that is tabulated in Jones & O'Dell (1977a) and depends on the value of the spectral index α. The term g (ϑ), represents the electron pitch angle distribution. It is evident, upon inspection of equations: C.1, C.2, C.3, and C.4, that τ F , τ C , and τ are together themselves functions of each plasma cell's: B , n e , and γ min . The ray profiles (similar to those presented in Figure 5) of these three variables (for the sightlines indicated by the red dots in the lower panels of Figure 4) are illustrated in Figure C.1. Upon inspection of these panels one can see that our choice of physical scaling (i.e., §2) results in magnetic field strengths that range from B 10 −2 − 10 0 Gauss and electron number densities that range from n e 10 0 − 10 2 cm −3 . These values are in rough agreement with the plasma conditions that have been inferred in blazar jets from theoretical modeling of shock acceleration and turbulence (see, e.g., Marscher 2014).
We present here an alternate physical scaling method (commonly used in magnetohydrodynamic simulations) which generates a similar magnetic field scale factor for comparison to our PIC approach ( §2). In an RMHD simulation one typically defines the following three scale factors: (see the PLUTO 3 code manual for further discussion). All other physical scaling quantities can be derived from these three unit values. In particular, the magnetic field strength scale factor (in cgs) may be written as: As discussed in §2, we set v o = 3×10 10 cm s −1 . Again, under the assumption of charge neutrality (n e − n e + +n p + ), and invoking the normal plasma case, in which n p + n e + , one may make the assumption that: n p + n e − . While the electrons dominate the nonthermal synchrotron emission, the protons (in contrast) dominate the thermal fluid dynamics. It follows that: ρ o = n o m p , where n o 10 1 cm −3 (i.e., setting the proton number density equal to the electron value from §2). This results in ρ o 1.7 × 10 −23 g cm −3 , having used m p = 1.6726231×10 −23 g cm −3 . Inserting v o and ρ o into Equation C.8 yields a magnetic field strength scaling factor of B o 0.4 Gauss which is in agreement with the scaling value we obtain in §2 (Equation 10) after invoking a fiducial jet plasma beta of β e 10 −3 . (ii) With ϑ determined (cell-by-cell) the radiative transfer is then computed across each cell. These calculations are done in a cell specific 'Stokes U = 0' linear polarization basis. The analytic solution used in our computations is presented in Jones & O'Dell (1977a) (Appendix B) and summarized in MacDonald & Marscher (2018) (Appendix A). (iii) Once the radiative transfer across a cell is complete, we then rotate the cell specific linear polarization basis onto a generalized observer's plane (illustrated in Figure D.1). In particular, the angle of this rotation (φ rot ) is: x + B y − sin(θ obs ) B y sin(θ obs ) − B z cos(θ obs ) 2 + B z − cos(θ obs ) B z cos(θ obs ) − B y sin(θ obs ) and is equal to the angle the projected magnetic field B proj =B − (B ·k )k makes with respect to the x-axis (i.e., west) on the observer's plane. This rotation is then applied to the Stokes parameters in the following manner: (D.5) (iv) These 'generalized' Stokes parameters are then recorded and passed on to the next cell via our slow-light interpolation scheme at which point we rotate back into the jet plasma frame and repeat steps (i-iii) for the next plasma cell along each sightline.  Figure 4. The red dashed line indicates a numerical limit below which we enforce the criterion I ν ,Q ν ,U ν ,V ν → I 0 ν ,Q 0 ν ,U 0 ν ,V 0 ν directly in our PIC ray-tracing calculations. Right panel: Variation of the minimum frequency (ν min ) within the normal plasma (e − -p + ) jet along the sightline indicated by the red dot in the lower left panel of Figure 4. The red dashed line indicates the emitted frequency ν em in the plasma frame corresponding to ν obs = 230 GHz. In both panels, the radiative transfer progresses from left to right: starting on the far side of the jet relative to the observer (cell 0), and then advancing through the plasma toward the near side of the jet (cell 240).
algorithm for comparison to the calculations presented here. We also plan on implementing an integration scheme (similar to the methods presented in Ruszkowski & Begelman 2002;Dexter 2016;Moscibrodzka & Gammie 2018) to solve Equation 14 numerically, thus providing a further comparison to the analytic expressions of Jones & O'Dell (1977a) used in this manuscript.

Appendix F: Particle-in-cell Spectropolarimetry
As a further comparison between the radiative properties of the normal plasma jet and the pair plasma jet we generate spectra of each jet's: (i) integrated total intensity -I , (ii) integrated fractional linear polarizationm l , (iii) integrated EVPA -χ , and (iv) integrated fractional circular polarizationm c , over the frequency range ν obs = 86 − 230 GHz (see Figure F.1), when each jet is oriented edge-on to the observer as shown in Figure 4 (i.e., like blazars). In particular, we follow the formalism of Kim et al. (2019) and compute integrated Stokes values for each simulation in the following manner: where I i , j , Q i , j , U i , j , and V i , j are summations of the map pixel values contained within the outermost Stokes I contour of each ray-traced image (conservatively set at 20% of each map's peak value). The term A pixel (= 0.4 × 0.4 square mas) denotes the angular extent of each pixel in our ray-traced images and is computed from the ratio of the RADMC-3D pixel size (∼ 3 × 10 4 cm) to the source distance (∼ 1.5 × 10 13 cm; 1 AU). The term A beam (= π ψ maj × ψ min /4 ln 2) denotes the angular extent of the convolving beam, where ψ maj and ψ min are the FWHM of the major and minor beam axes, respectively. For this analysis we have used a fixed circular beam of ψ maj = ψ min = 3.5 mas (illustrated in the lower left of each panel in Figure 4). From these integrated Stokes values we compute the following polarimetric quantities: The above quantities are then recorded for each ray-tracing calculation at individual frequencies ranging from ν obs = 86 − 230 GHz (see Figure F.1).
The normal plasma jet and the pair plasma jet are both optically thin in this frequency range. Both spectra exhibit a general trend of decreasing intensity (top panels) with increasing fractional linear polarization (upper middle panels) accompanied by decreasing fractional circular polarization (bottom panels) in agreement with synchrotron theory. An anticorrelation between linear and circular polarization across frequency has been detected in the quasar PKS B2126-158 with high-precision CP measurements made using the Australian Compact Telescope Array (ATCA) (O'Sullivan et al. 2013). The solid red and dashed red lines in the lower panels highlight theoretical predictions for the frequency dependence of m c (ν) in the limits of intrinsic emission and emission dominated by Faraday conversion in the high-rotation limit, respectively (see Pacholczyk 1970;Jones & O'Dell 1977a). Clearly, the emission produced in our PIC models is intrinsic in origin due to the very small Faraday depths through the jet plasma along each sightline (see Figure 5). Fast-light Stokes I image of the normal plasma jet, but in which an artificial Gaussian noise floor (set at ∼ 10 −6 Jy/beam) has been included to illustrate the extreme interferometric array sensitivity needed to detect our PIC jets even at such close proximity.