subtleties

Hammurabi's integrator (class_Integrator.cpp) has two main running modes. Its default behavior is quite straightforwardly to compute the emission in each cell numerically by taking a single magnetic field vector (returned by the magnetic field module function class_B_field2::return_B_cart()) at each point and computing the Stokes parameters that depend on that vector and the given line of sight vector, and then summing each of those cells along the LOS.

Alternatively, it has an "analytic" mode that follows the method used in Jansson & Farrar (2012, ApJ 757,14, aka JF12a). This method takes advantage of an observational distinction made between magnetic field components as: a coherent field, an isotropic random field, and an anisotropic random field (which JF12 called "striated"). See Fig. 1 of Jaffe et al. (2010, MNRAS 401, 1013, aka Jaffe10). First, an explanation of how the magnetic field module handles these components:

Hammurabi computes the coherent magnetic field components analytically based on models implemented in the class_B_field2.cpp module. In other words, there is a simple mathematical description for the field vector based on a given position in the Galaxy. Alternatively, the user can specify an external file that contains a Cartesian grid where the vector field is defined in each grid cell. It will then interpolate from that grid onto the (non-Cartesian) integration grid.

The class_B_field2.cpp also includes a module to compute an isotropic Gaussian random field. It then scales this GRF following any of several (analytic) models for the distribution of the average strength of the isotropic random component at a given position. For example, this Gaussian random field is usually modeled as having a roughly exponential profile in Galacto-centric radius and height.

The magnetic field module then returns to the integrator the vector sum of the coherent and random components at each position requested. The third component, the anisotropic random component, can be handled several ways. For the moment, we will assume it is zero and return to it later.

The "analytic" versus "numeric" modes of the hammurabi integrator are identical if only a coherent magnetic field component is included. The difference is in how the random field components are handled. In the numeric case, the vector sum of the coherent field and the single realization of the Gaussian random field (scaled according to the average strength of the random component at the given position) is given to the integrator, which computes the Stokes component from that position and adds them all along the LOS. But in the analytic case, the magnetic field module does not compute a Gaussian random field. It simply returns to the integrator the average strength of the random field at the given position according to whichever model has been specified. (In practice, it returns the vector with x,y,z components each equal to that average strength/sqrt(3).) The integrator then determines the average amount of synchrotron emission from this cell, in total intensity following the expressions above Eq. 1 of Jansson & Farrar (ApJL, 2012, 761, L11, aka JF12b) for B_perp, the average strength of the component of B perpendicular to the LOS. The important part of that equation is the factor of 2/3, which is the average increase in B_perp due to the addition of the random component. This fudge is only used for the total intensity. The polarized intensity is computed solely from the regular field component, so the result therefore represents the ensemble average of the Stokes Q and U. (The isotropic random component would cancel out in each of Q and U separately in the ensemble average.) Note that the RM returned is also the ensemble average RM, which is identically the RM from the coherent field alone.

But it's important when thinking about these issues to recall that the random component turns the LOS integration of the polarization signal into a 2D random walk. The starting point of the random walk is basically the polarization (Q,U) one would see if only the coherent component of the magnetic field were present. The random component then causes the polarization vector to random walk in the Q-U plane. On average, the polarization vector would be the same as in the presence of the coherent field alone. But with the additional random component, the polarization vector random walks a distance away from that average which is proportional to the strength of the random component and the square root of the number of steps in the random walk, i.e. the number of turbulent cells along the LOS. (Though in reality or in a GRF, it's not discrete, this is a useful simplification for visualizing.) Therefore, even if there is no coherent field, the ensemble average polarization signal would not be zero; the average Q and U would be, but not the average P=sqrt(Q^2+U^2). (This is analogous to noise bias in any quantity dependent on the square of a noisy quantity.) So the products of hammurabi in the analytic case are only representative of the ensemble average Q and U, but do no reflect the average PI.

(Note that the first version of this note got this precisely backwards. My apologies.)

Note also that the RM is a 1D random walk. The average RM is indeed that assuming only a coherent field, but adding a random field means that the RM deviates by that average by an amount proportional again to the strength of the random field and the sqrt of the number of turbulent cells along the LOS. For the analytic integration, RM is computed using only the coherent field component. For the numerical integration, RM is computed numerically and will include the variance due to the random component; if you generate an ensemble of realizations, their average should be equal to the analytic result.

Note that hammurabi's parameter use_B_analytic=T can be specified
for any of it's models, not just the JF12 model (which is implemented as
B_field_type=7.) Note also that the JF12 model can be computed
numerically, and will be by default since the use_B_analytic=F
is the default behavior. This means that the coherent field of JF12a is
added to a numerically generated GRF realization scaled following the
model described in JF12b. In this case, the I,Q, and U are each correct
for a single realization of the coherent field plus the isotropic random component.

Now we come to the additional complication of the anisotropic random field component, called "striated" by JF12 and "ordered random" by Jaffe10. This component of the magnetic field retains an ordered orientation but changes sign stochastically. The JF12 model implemented this component as a fudge factor. Because this component contributes to total and polarized intensity but not to RM (see Jaffe10 Fig. 1 again), their analysis simply scaled up the analytically computed emission by a factor of beta; see again Eq. 1 in JF12b or section 5.2 of the earlier JF12a. This was presumably done in the likelihood evaluation outside of hammurabi. But now the public version of hammurabi implements this internally so that the resulting I, Q, and U maps reflect this increased component. If the parameter B_analytic_beta is not zero, the coherent field strength is multiplied by the sqrt() of it. Note that this then messes up the RM by the corresponding factor. In order to avoid undue complexity in the hammurabi integrator, the code must be run twice to produce all observables: once for the synchrotron emission using B_analytic_beta=1.36 and once for the RMs with B_analytic_beta=0. This applies regardless of whether the computation is done analytically or numerically, because it will be simply a multiplicative factor applied to the coherent field component in the magnetic field module. The physical interpretation of this fudge factor is that there is additional field strength that is morphologically identical to the coherent field component but changes sign stochastically. The hammurabi implementation of the striated field using B_analytic_beta>0 is always computed analytically, not numerically.

The Jaffe10 analysis also includes an anisotropic random component. This model is B_field_type=10 in hammurabi. But the anisotropic component is generated from the random component instead of from the coherent component, as in the JF12 case. The reality is likely to be somewhere between these two approaches, and both can indeed be used together. As in the case of the JF12 fudge, the Jaffe10 fudge is also now a generic parameter of hammurabi: B_f_ord. See Eq 7 of Jaffe10. This adds a field component parallel to the coherent field but proportional (by f_ord) to the random field component at that position. BUT note that this does not work in analytic mode; the average of N numeric realizations does not match the analytic result, presumably due to a bug I have not found. So don't use B_f_ord in analytic mode.

Note also that the analytic integration mode depends on the assumption that the emissivity is proportional to the strength of the magnetic fields squared. (Or more generically, a power law with non-zero index B^n. See Longair, High_Energy_Astrophysics Eq 8.87. But this general case is not implemented in hammurabi, which assumes the square.) Since dust emission is not proportional to the magnetic field strength, this analytic method cannot be used for dust polarization. If anybody can derive an analytic expression for the ensemble average polarization from dust, please do let me know. (A theorist collaborator had a go, and though the idea seems simple enough at the start, the math quickly becomes a mess that we did not know how to solve.)