Implementation details

This section describes the technical details of the simulation flow for a standard model component and for the conversion of a Chandra observation. In order to take advantage of the efficient array manipulation capabilities provided by numpy, the entire implementation is vectorized, i.e. we don’t have an explicit event loop in python.

Standard model component

The basic flow of the simulation for a single standard model component is coded in ixpeobssim.srcmodel.roi.xModelComponentBase.rvs_event_list(). Mathematically speaking, the simulation algorithm can be spelled out in the form of the following basic sequence:

  1. Given the source spectrum \mathcal{S}(E, t) and the effective area A_{\rm eff}(E), we calculate the count spectrum as a function of the energy and time:

    \mathcal{C}(E, t) = \mathcal{S}(E, t) \times A_{\rm eff}(E)
\quad [\text{s}^{-1}~\text{keV}^{-1}].

    In case of a source with cosmological redshift Z or hydrogen column density n_H different from zero, these two quantities are used to correctly re-scale \mathcal{C}(E, t):

    \mathcal{C}'(E, t) = \mathcal{C}(E(1 + E), t) \exp{(-n_H \sigma(E))}

  2. We calculate the light curve of the model component (in counts space) by integrating over the energy:

    \mathcal{L}(t) = \int_{E_{\rm min}}^{E_{\rm max}} \mathcal{C}'(E, t) dE
\quad [\text{Hz}].

  3. We calculate the total number of expected events N_{\rm exp} by integrating the count rate over the observation time:

    N_{\rm exp} = \int_{t_{\rm min}}^{t_{\rm max}} \mathcal{L}(t) dt

    and we extract the number of observed events N_{\rm obs} according to a Poisson distribution with mean N_{\rm exp}.

  4. We treat the count rate as a one-dimensional probability density function in the random variable t, we extract a vector \hat{t} of N_{\rm obs} values of t according to this pdf—and we sort the vector itself. (Here and in the following we shall use the hat to indicate vectors of lenght N_{\rm obs}.)

  5. We treat the array of count spectra \hat{\mathcal{C}}(E, \hat{t}), evaluated at the time array \hat{t}, as an array of one-dimensional pdf objects, from which we extract a corresponding array \hat{E} of (true) energy values. (In an event-driven formulation this would be equivalent to loop over the values t_i of the array \hat{t}, calculate the corresponding count spectrum

    \mathcal{C}_i(E, t_i)

    and treat that as a one-dimensional pdf from which we extract a (true) energy value E_i, but the vectorized description is more germane to what the code is actually doing internally.)

  6. We treat the energy dispersion \hat{\mathcal{D}}(\epsilon; \hat{E}) as an array of one-dimensional pdf objects that we use to extract the measured energies \hat{\epsilon} and the corresponding PHA values.

  7. We extract suitable arrays of (true) \hat{\text{RA}}, \hat{\text{Dec}} values and, similarly to what we do with the energy dispersion, we smear them with the PSF in order to get the correponding measured quantities.

  8. We use the polarization degree P and angle \alpha of the model component to calculate the visibility \hat{\xi} and the phase \hat{\phi}_0 of the azimuthal distribution modulation, given the modulation factor \mu(E) of the polarimeter

    \hat{\xi} =
\hat{P}(\hat{E}, \hat{t}, \hat{\text{RA}}, \hat{\text{Dec}}) \times

\hat{\phi}_0 =
\hat{\alpha}(\hat{E}, \hat{t}, \hat{\text{RA}}, \hat{\text{Dec}}),

    and we use these values to extract the directions of emission of the photoelectron.

(For periodic sources all of the above is done in phase, rather than in time, and the latter is recovered at the very end using the source ephemeris, but other than that there is no real difference.)

For source models involving more than one component, this is done for each component separately, and the different resulting phothon lists are then merged and ordered in time at the end of the process.

Conversion of a Chandra observation

In this case the starting point is a Chandra photon list, publicly available on the on-line database, which includes temporal, energetic and spatial information for each detected event. Being the Chandra angular and energetic resolution (respectively HPD \sim 1 arcsec and FWHM \sim 5% at 5.9 keV) much better than those expected for the IXPE mission), energies E_C and positions RA_C, Dec_C measured by Chandra can be assumed as Monte Carlo truth.

From the implementation standpoint, the conversion is coded in ixpeobssim.srcmodel.roi.xChandraObservation.rvs_event_list() and can be summarized in the following steps:

  1. We compute the ratio T_R between the Chandra observation time T_C and the provided IXPE simulation time T_X:

    T_R = T_C / T_X

  2. We scale Chandra detected events by first concatenating the measured arrays (E_C, RA_C and Dec_C) as many times as the integer part of T_R (denoted \lfloor T_R \rfloor) and then appending to the end of the sequence a portion whose relative length is equal to the fractional part of T_R (called \lbrace T_R \rbrace). For instance, in case of the energy:

    \hat{E} = \underbrace{E_C \: E_C \cdots E_C}_{\lfloor T_R \rfloor}
\underbrace{E_C}_{\lbrace T_R \rbrace}.

  3. We calculate the effective area ratio A_R between one IXPE telescope and Chandra as a function of energy and off-axis angle \theta:

    A_R(E, \theta) = \dfrac{A_{\mathrm{eff},X} (E,\theta)}

  4. This ratio allows to conveniently down-sample the Chandra time-scaled photon list. In practice, this is done by throwing an array of random numbers \hat{r} between 0 and 1, and accepting as IXPE events only those that meet the condition:

    \hat{r} < A_R(\hat{E},\hat{\theta})

  5. We randomly extract the event times \hat{t} uniformly between the starting and ending observation time.

  6. From this point, once the arrays \hat{E}, \hat{t}, \hat{RA} and \hat{Dec} of true values are selected, the simulation goes ahead as in the standard model component case. In particular, we smear these quantities with the IXPE response functions and we extract the emission angle \hat{\phi} based on the input polarization model.

In case of definition of multiple sub-regions of the ROI, the conversion is done separately for each of them (skipping those flagged as to be removed), and the resulting photon lists are merged and ordered in time at the end of the process.