Balmer Filaments in Tycho's Supernova Remnant: An Interplay between Cosmic-ray and Broad-neutral Precursors

Supernova remnant (SNR) shocks are suspected to be the long-sought Galactic cosmic-ray (CR) sources. Observational evidence for particle acceleration at work in SNRs has been seen in synchrotron emission from electrons, gamma radiation, signatures of amplified magnetic fields, higher compression ratios than predicted by jump conditions, and optical line profiles from shocks propagating in a partially ionized plasma (for a review, see Helder et al. 2012). Since most of the CRs that we detect on Earth consist of protons, finding evidence for accelerated CR protons in SNRs with energies up to the "knee" ($\sim 1\,\mathrm{PeV}$) is thus of particular importance. Insight into CR protons is possible through either the gamma-ray spectrum as a result of neutral pion decay (Ackermann et al. 2013) or the optical wavelength window by carefully analyzing Hα-line profiles (Helder et al. 2009, 2010; Nikolić et al. 2013). The latter shock emission, usually also referred to as optical emission from Balmer-dominated shocks (BDSs), is the subject of the study of this paper.

The spectra of BDSs, typically observed around SNRs that originate from SN Ia explosions, show the presence of strong two-component hydrogen lines (Heng 2010). When a shock wave encounters partly ionized interstellar medium (ISM), the cold pre-shock hydrogen atoms overrun by the shock can either be excited by hot post-shock gas, resulting in the narrow Hα-component emission, or enter a charge exchange (CE) process with the hot post-shock plasma, producing hot neutrals whose collisional excitation then gives rise to the broad Hα component. The two-component Hα-line parameters provide valuable information on the existence of CR precursors in the shocks. A narrow line broadened beyond 10–20 km s⁻¹ gives direct evidence of the presence of non-thermal particles in the shock precursor (Morlino et al. 2013). This is due to the fact that hydrogen atoms are ionized at temperatures larger than ≈10,000 K, but also because the lifetime of the neutral hydrogen in the post-shock region is too short for the collisional interaction to broaden the line profiles (Smith et al. 1994). The CRs will heat the cold neutrals in the ISM, resulting in the broadening of the narrow Hα line, but the CRs will also reduce the broad Hα-line width by removing energy from the protons in the post-shock region.

Several authors proposed that narrow-line broadening can also arise from a broad-neutral (BN) precursor: hot neutrals created in CE processes between hot protons and slow neutrals streaming to the pre-shock region (Hester et al. 1994; Smith et al. 1994). Recent theoretical studies (Blasi et al. 2012; Morlino et al. 2012) show that a BN precursor does not broaden the narrow component, but rather introduces a third intermediate component with the FWHM of around ∼ 150 km s⁻¹ and depends on the shock speed. The reason is that only a small number of incoming neutrals interacts with ions in the BN precursor because its extent, which corresponds to the interaction length of the returning neutrals, is much smaller than the CE interaction length of the incoming neutrals. Therefore, Balmer lines can be used to study the microphysics of collisionless shocks and are currently the only means that give insight into the collisionless shocks.

We observed Tycho's SNR, which has already been well studied across all wavelength ranges (Reynolds & Ellison 1992; Bamba et al. 2005; Lee et al. 2004; Stroman & Pohl 2009; Katsuda et al. 2010; Acciari et al. 2011; Tian & Leahy 2011; Giordano et al. 2012), and evidence for particle acceleration in the shocks of Tycho's SNR, including acceleration up to the knee in the CR spectrum, was found (Warren et al. 2005; Cassam-Chenaï et al. 2008; Eriksen et al. 2011; Slane et al. 2014). In 1572, the star exploded as an SN Ia, leaving a remnant at an estimated heliocentric distance of 2.3 ± 0.5 kpc (Chevalier et al. 1980).¹¹ At that distance, the remnant's diameter of 8' corresponds to ≈5 pc. Density gradients in the medium around the remnant (Williams et al. 2013) modified the evolution of the shock, which in turn resulted in the asymmetric remnant. The lower shock velocity inferred in the northeastern (NE) part suggests that the shock interacted with a dense ambient medium, namely a diffuse cloud (Reynolds & Keohane 1999; Lee et al. 2004).

Previous optical studies of Tycho's SNR have shown indications for CRs (e.g., Ghavamian et al. 2000; Lee et al. 2007, 2010). However, these studies focused on the Hα-bright, but very complex, "knot g," where multiple or distorted shock fronts can contribute to the measured narrow-line (NL) broadening and thus partially mimic the effect that CR acceleration would have. Using the Fabry–Pérot instrument GHαFaS (Galaxy Hα Fabry–Pérot System) on the William Herschel Telescope (WHT), we observed a great portion of the shock front in the NE region of the remnant. The high spatial and spectral resolution, together with the large field of view (FOV) of the instrument, allow us to measure the narrow Hα-line width across individual parts of the shocks simultaneously, and thereby study the indicators of CR presence in a large variety of shock front conditions. In particular, the spatial resolution allows us to distinguish intrinsic line broadening from line broadening originating in geometric distortions and differential kinematics. Moreover, our observational setup provides us with unique insight into the existence of an intermediate component along the entire filament previously only reported for the bright "knot g" (Ghavamian et al. 2000; Lee et al. 2007). Apart from vastly enhancing the amount and quality of the spectroscopic data available for the NE filament, our study also improves the analysis: instead of fitting line models, we employ Bayesian inference to obtain full information and realistic uncertainties on the line parameters, as well as quantitative, reliable evidence for the presence of an intermediate line (IL) originating in a BN precursor and multiple shock fronts. The interpretation of the results is based on predictions of the state-of-the-art shock models that include the effects of BN and CR precursors on the observed Hα profiles (Morlino et al. 2012, 2013). With the GHαFaS spectral coverage of approximately 400 km s⁻¹, we are not able to resolve the broad Hα component that was found to be about ≈2000 km s⁻¹ in previous studies (Chevalier et al. 1980; Ghavamian et al. 2001).

In order to resolve the narrow Hα lines along the rim of Tycho's SNR, we used the instrument GHαFaS mounted on the Nasmyth focus of the 4.2 m WHT (Hernandez et al. 2008), which operates at the Observatory del Roque de Los Muchachos in La Palma, Canary Islands. GHαFaS is a Fabry–Pérot interferometer–spectrometer with an FOV of 34 × 34. Its detector is an Image Photon Counting System (IPCS) for which the absence of readout noise is an advantage for observations of diffuse emission from extended objects. IPCS cameras are almost insensitive to CRs and thus do not require CR rejection. We used a high-resolution mode, acquiring data on 1024 × 1024 pixels² with R ∼ 21,000 resolving power and a pixel scale of nearly 02 pixel⁻¹. The free spectral range (FSR) of the etalon was 8.56 Å or 392 km s⁻¹ centered at 6561 Å and split into 48 channels that differ in their central wavelength, leading to a sampling velocity resolution of 8.16 km s⁻¹. The instrument response function is well approximated by a Gaussian with FWHM of 19 km s⁻¹  (Blasco-Herrera et al. 2010).

The observations were conducted on 2012 November 15–19 under FWHM ≃ 1'' seeing conditions. Successive exposures differ by one channel, and thus 48 successive exposures complete one cycle that covers the full spectral range. The total integration time was ≈9.6 hr, comprising 72 cycles and 3456 $10\,{\rm{s}}$ exposures. Observations conducted over several cycles provide homogeneous airmass and atmospheric conditions for all channels. We reduced the data (see Figure 1) by first applying the phase correction to all exposures individually, where we follow the standard procedure for GHαFaS data described in Hernandez et al. (2008). The phase correction is a process of designating the photons' positions for the interference rings and assigning the corresponding wavelength to each position (x,y) on the image, ${\lambda }_{i}={\lambda }_{i}(x,y)$. As a result, from each exposure ${D}_{i}(x,y)$, we build a data-subcube ${D}_{i}(x,y,\lambda )$ with 48 monochromatic images. In order to use the largest possible GHαFaS FOV, we did not use the optical derotator. Therefore, we have to align and derotate the observed data subcubes before co-adding them. To determine the exposures' relative pointing and orientation, we measure the centroid positions of bright point sources (stars) on a stack that combines each exposure with the 2 × 4 exposures that precede and succeed it. Thus, we are assured that at least three bright point sources are detected with enough flux for a $\lt 1$ pixel centroid precision. This requirement means that we have to discard the first four and the last four D_i of each "run" of consecutive observations, but this concerns only 57 exposures (2%). Along with removing 13 cycles that suffer from reflected light or other defects that we noticed in visual inspection of the data, and also excluding seven cycles that lack reference sources for derotation, we retain 2439 exposures in 52 cycles. The final data cube $D(x,y,\lambda )$ results from summing all data subcubes and consists of 48 calibrated constant-wavelength slices.

Zoom In Zoom Out Reset image size

Besides the data cube, that is, the stack of all aligned data subcubes, we produce a background cube and a flatfield cube. To this end, we model the background flux in individual exposures as well as the flatfield image (the position-dependent throughput of the optical system), and subsequently process the individual background and flatfield frames in exactly the same manner as the corresponding data frames. By constructing the co-added background and flatfield cubes and including them in our parametric models of the observed shock emission, as opposed to subtracting the background from the individual exposures and dividing them by the flatfield image before calibration and co-addition, we preserve the photon (Poisson) statistics in the data cube and simultaneously account for the variable effective exposure time in the stack of aligned subcubes. Details on the modeling of the flatfield and background frames are given in Appendix A.

Our analysis is guided by the following goals and principles:

• 1.  
  To determine the Hα NL width across a maximal area of the observed shock.
• 2.  
  To achieve maximum spatial resolution, implying minimal binning (bin size) and signal per bin.
• 3.  
  To still extract line parameters reliably and, in particular, characterize their uncertainties accurately.
• 4.  
  To do so even when including up to two additional lines (10 model parameters).
• 5.  
  To compare single-NL and multi-line models and to quantify their relative evidence.

In order to achieve these aims, we decide to perform parameter estimation and model comparison using Bayesian inference instead of traditional (maximum-likelihood, minimum-${\chi }^{2}$) fitting routines. We also account for the Poisson statistics of the data, as opposed to the often tacitly applied Gaussian approximation. The details of our method can be found in Section 3.4, and Appendices E and F.

3.1. Motivation for a Multi-line Analysis

3.2. Definition of Models

For each location (bin), we consider several parametrized models (S) to characterize the shock Hα emission. Regardless of the parametrization (type of model) or specific parameter values (θ), we factor in the local flatfield spectrum (F) and add to it the observed background spectrum (B) before comparing the model with the data. This has the advantage of preserving the correct photon statistics, and contrasts with the common approach of subtracting the background from the data and dividing by the flatfield before modeling. Hence, for a given location (bin), the full model is represented by $M(\theta )\,=S(\theta )\times F\,+\,B$, while by model we mostly refer just to S, the intrinsic or "source" component. Note that in this expression, F and B are the result of binning the flatfield cube and the background cube in the same way as the data, and therefore they only depend on wavelength (λ) for the bin that M describes. F and B are constructed separately from the data and inserted into the model without free parameters. More information on how we established F and B can be found in Section 2 and Appendix A.

Every S consists of an NL with a Gaussian profile as well as a constant component (c) that accounts for the sum of the continuum level and the broad Hα line (BL), which is several times wider than our spectral range. Depending on its type, S may further include one or two Gaussian components that represent an IL or an additional NL. Overall, in each bin we therefore have four different models S to compare with the data spectrum:

• 1.  
  NL—constant plus single narrow line.
• 2.  
  NLNL—constant plus two narrow lines.
• 3.  
  NLIL—constant plus one narrow and one IL.
• 4.  
  NLNLIL—constant plus two narrow lines and one IL.

In connection with calculating the models' relative evidence, we also consider a "no-line model" (0L; i.e., with the constant spectrum as the only component), which gives an auxiliary baseline for the relevance of at least the NL being present in the data. A possible choice of model parameters is the component fluxes (${f}_{c}={f}_{\mathrm{constant}}$ and ${f}_{i}={f}_{\{\mathrm{NL},\mathrm{NL}1,\mathrm{NL}2,\mathrm{IL}\}}$), line centroids (${\mu }_{i}$), and the lines' FWHM (W_i). The general form of this model is

$$\begin{eqnarray}&&S(\lambda )={f}_{c}+{(2\pi )}^{-1/2}\displaystyle \sum _{i}{f}_{i}\exp (-{(\lambda -{\mu }_{i})}^{2}/(2{\sigma }_{i}^{2}))/{\sigma }_{i}\,,\end{eqnarray} \tag{ 1 }$$

where ${\sigma }_{i}=\sqrt{{W}_{i}^{2}+{W}_{\mathrm{instr}}^{2}}/\sqrt{4\mathrm{ln}4}$ is the observed (instrumentally broadened) Gaussian dispersion and W_instr the FWHM of the instrumental response.

In practice, our models employ a transformed version of those parameters, which has the advantage of more direct interpretation, and a simpler functional form of the desired parameter priors (see the next paragraph). For example, instead of the two NL centroids, ${\mu }_{\mathrm{NL}1}$ and ${\mu }_{\mathrm{NL}2}$, the NLNL and NLNLIL models use the NL centroid mean ($\langle {\mu }_{\mathrm{NL}}\rangle$) and the separation between the NLs (${\rm{\Delta }}{\mu }_{\mathrm{NL}}$). Similarly, the IL centroid is specified by its offset from the NL centroid or NLNL centroid mean, ${\rm{\Delta }}{\mu }_{\mathrm{IL}}={\mu }_{\mathrm{IL}}-{\mu }_{\mathrm{NL}}$. We replace the component fluxes by the total flux and the components' flux fractions, and use the logarithm of the total flux and line widths as they are strictly positive quantities. A more detailed account of the definition of model parameters can be found in Appendix D.

The parameter priors, $P(\theta )$, are an integral part of the model definition: they encapsulate what we know (or assume) about the relative probabilities of the parameter values a priori, before considering the data. In particular, they can impose parameter boundaries by way of being zero outside of those boundaries. The line centroid parameters are effectively restricted to our spectral window. In our models, the FWHM parameters W_NL and W_IL are limited to [15, 100] km s⁻¹ and [100, 350] km s⁻¹, respectively. The lower boundary of W_NL reflects the lower limit of the pre-shock temperature (≈5000 K), while the upper boundary is based on theoretical models that include the effects of the CR precursor (Morlino et al. 2013). The W_IL range is the theoretical expectation for shock velocities around 2000 km s⁻¹ and a range of shock parameters (see Figure 10 in Morlino et al. 2012).

Inside of the parameter boundaries, we desire to not strongly favor any particular parameter values a priori. However, on physical grounds, we prefer a smooth transition of the prior to zero for parameter values approaching the boundaries (see Figure 2). We therefore employ shifted and scaled Beta distributions with $\alpha =\beta =1.5$ for the centroid parameters and logarithmic line widths, and a Dirichlet distribution with $\alpha =1.5$ for the component flux fractions, which have the constraint of being summed up to unity. For comparison, the special case of "flat" Beta and Dirichlet distributions would be realized by setting $\alpha (=\beta )=1.0$. Thus, our choice slightly favors the center of the allowed parameter range. The logarithm of the total flux has a flat unbounded prior. The model parameters and their priors are summarized in the Table 1, and detailed definitions are presented in Appendix D.

Zoom In Zoom Out Reset image size

Table 1.  Model Parameters and Their Prior PDFs

Parameters	Meaning	Prior
ln $({F}_{\mathrm{tot}})$	Natural log-based total flux	flat prior
${f}_{i},i=[{\rm{c}},\mathrm{NL},\mathrm{NL}1,\mathrm{NL}2,\mathrm{IL}]$	Flux fractions	Dirichlet prior: ${\prod }_{i}{f}_{i}^{\alpha -1}$, $\alpha =1.5$; ${\sum }_{i}{f}_{i}=1,{f}_{i}\in (0,1)$
${\mu }_{{NL}}^{{\prime} },\langle {\mu }_{{NL}}^{{\prime} }\rangle$	NL centroid, NL centroid mean	Beta prior: ${x}^{\alpha -1}{(1-x)}^{\beta -1}$, $\alpha =\beta =1.5;$ x ∈ (0, 1)
${\rm{\Delta }}{\mu }_{\mathrm{NL}}^{{\prime} }$	Separation between the two NLs
${\rm{\Delta }}{\mu }_{\mathrm{IL}}^{{\prime} }$	IL centroid offset from NL centroid (mean)
${w}_{\mathrm{NL}(\mathrm{IL})}^{{\prime} }$	NL (IL) natural log-width
Model	Model Parameters
NL	ln $({F}_{\mathrm{tot}}),\,{f}_{\mathrm{NL}},\,{f}_{{\rm{c}}},\,{\mu }_{\mathrm{NL}}^{{\prime} },\,{w}_{{NL}}^{{\prime} }$
NLNL	ln $({F}_{\mathrm{tot}}),\,{f}_{\mathrm{NL}1},\,{f}_{\mathrm{NL}2},\,{f}_{{\rm{c}}},\,\langle {\mu }_{{NL}}^{{\prime} }\rangle ,\,{\rm{\Delta }}{\mu }_{\mathrm{NL}}^{{\prime} },\,{w}_{\mathrm{NL}1}^{{\prime} },\,{w}_{\mathrm{NL}2}^{{\prime} }$
NLIL	ln $({F}_{\mathrm{tot}}),\,{f}_{\mathrm{NL}},\,{f}_{\mathrm{IL}},\,{f}_{{\rm{c}}},\,{\mu }_{\mathrm{NL}}^{{\prime} },\,{\rm{\Delta }}{\mu }_{\mathrm{IL}}^{{\prime} },\,{w}_{\mathrm{NL}}^{{\prime} },\,{w}_{\mathrm{IL}}^{{\prime} }$
NLNLIL	ln $({F}_{\mathrm{tot}}),\,{f}_{\mathrm{NL}1},\,{f}_{\mathrm{NL}2},\,{f}_{\mathrm{IL}},\,{f}_{{\rm{c}}},\,\langle {\mu }_{{NL}}^{{\prime} }\rangle ,\,{\rm{\Delta }}{\mu }_{\mathrm{NL}}^{{\prime} },\,{\rm{\Delta }}{\mu }_{\mathrm{IL}}^{{\prime} },\,{w}_{\mathrm{NL}1}^{{\prime} },\,{w}_{\mathrm{NL}2}^{{\prime} },\,{w}_{\mathrm{IL}}^{{\prime} }$

Note. All parameters apart from ln $({F}_{\mathrm{tot}})$ are defined in the (0,1) range (see Appendix D for notation).

Download table as:  ASCIITypeset image

3.3. Binning

We analyze two shock filaments, one in the more eastern part of the NE rim which contains "knot g," and the other in the more northern part (Figure 3). We use the Weighted Voronoi Tessellation (Diehl & Statler 2006) with the adaptive bin size to spatially bin the pixels and obtain a signal-to-noise ratio (S/N) of ≳10 (Appendix C) in the wavelength-integrated signal that remains after subtraction of the background. This implies an average minimum S/N of 1.4 per spectral element. Due to the seeing of 1'', we require at least 5 pixels across, so that for a round bin this implies a minimum of 19 pixels. We exclude bins that would require $\gt 400$ pixels for our target S/N, so that unaccounted-for residual background variations, which we estimate to be at most $\sim 2 \%$ of the background level, do not significantly effect our measurements. Following these criteria, we study 73 Voronoi bins in the eastern and 9 Voronoi bins in the northern filament.

Zoom In Zoom Out Reset image size

3.4. Parameter Estimation and Model Comparison

For each Voronoi bin, we want to find which of the four models, M (see Section 3.2), and which vector of model parameters, θ, best explain its data (spectrum). Since the data are noisy, model comparison and parameter estimation are inherently probabilistic. For both tasks and the reasons discussed below, we use Bayesian inference.

The "standard" approach to estimate parameters is the maximum-likelihood (or minimum-${\chi }^{2}$) method. It relies on a well-defined likelihood maximum (mode) and the convergence of the optimizer deteriorates when high noise or multiple modes are present. Both conditions are met in our study: we desire high spatial resolution and therefore small bin size, implying low S/N. Our models are nonlinear and comprise up to 10 parameters, implying generally multiple modes. We wish to characterize all those modes, not just the "main" (global) maximum, toward reliable, high-confidence level parameter uncertainties ("errors"), instead of the minimal but common 68% confidence ("$1-\sigma$") error. Often, and in our models with their non-trivial likelihood function, error propagation over a large parameter range is cumbersome or impossible. Finally, the maximum of the likelihood does not provide a quantitative, well-defined measure for the relative probability of the different models with different parametrizations (model comparison). These circumstances make maximum-likelihood or other fitting methods insufficient for our purposes.

Bayesian inference, by contrast, provides the full, multivariate parameter probability distribution function (PDF), the so-called posterior $P(\theta | D,M)$, as well as its integral (marginalization) over all parameters, $P(D| M)$—the evidence. Evidences of models are the relative model probabilities, without reference to any specific (e.g., best-fit) parameter value. Bayes' theorem states that the posterior as a function of $\theta$ is

$$\begin{eqnarray}&&P(\theta | D,M)=\displaystyle \frac{P(D| \theta ,M)P(\theta | M)}{P(D| M)}.\end{eqnarray} \tag{ 2 }$$

It is proportional to the product of the likelihood, $L=P(D| \theta ,M)$, and prior, $P(\theta | M)$. L reflects the model and is the probability of the data for given model parameters and measurement errors. Our IPCS instrument counts photons, which are described by a Poisson distribution that therefore represents the measurement error, with expectation value and variance equal to the flux predicted by the spectral model. The prior is the parameter PDF that we know or assume before taking into account the data at hand. Apart from making these assumptions or knowledge explicit (fitting methods implicitly assume a flat prior in the chosen parameters), it has the advantage of naturally facilitating self-consistent parameter changes.

We represent the N-dimensional posterior PDF of the N model parameters as a sample, which we obtain using the Markov Chain Monte Carlo (MCMC) method. For details on the sampling algorithm, see Appendix E. The posterior can be summarized in many ways. One is the posterior maximum; we do not emphasize this as it is a relatively noisy estimator and its computation is not unique. Instead, we provide the median and the highest density (shortest) 95% confidence intervals of the one-dimensional marginalized distributions that result from integrating $P(\theta | D)$ over all but one parameter. We deem 95% to be the minimal confidence level worth quoting, and more reliable than the frequently employed 68% ("$1-\sigma$") level, which carries a high probability (32%) of not including the optimal parameter values.

For the parameter estimation, we are interested only in the posterior's shape (relative parameter probabilities), which does not necessitate normalization by the evidence. However, the Bayes factor, i.e., the evidence ratio of two models, is a probabilistically well-defined, quantitative measure for comparing models. For a given M, the evidence is defined as $P(D| M)={\int }_{\theta }P(D| \theta ,M)P(\theta | M)d\theta$. It is the probability of observing the data when assuming that the model is "true" but the parameters are not specified. In practice, the evidence integral is often high dimensional and therefore computationally intensive. In order to approximate it numerically, we use the cross-validation (CV) likelihood (Bailer-Jones 2012), particularly the leave-one-out (LOO) CV likelihood: $P(D| M)\,={L}_{\mathrm{LOO} \mbox{-} \mathrm{CV}}(D| M)$; see also Appendix F. To compute it, samples are drawn from the data partition posteriors instead of the prior as in "standard" evidence integrals, and it has the advantage that it depends on the prior only to second order. Each bin spectrum has 48 elements and the LOO-CV is applied to these elements by predicting each of them from the remaining 47 elements under the model M (marginalized over the parameters).

Since $P(D| M)$ can be a very large (or small) number, its absolute value is meaningless, and only relative values are needed for different models, we express it as the base-10 logarithm of the ratio with some reference model (Figure 4, Tables 2 and 3). As a matter of choice, we consider 0.5 dex log evidence differences as "significant" to clearly prefer one model over another. This choice is somewhat conservative; testing on simulated data reveals that we start to distinguish the correct models at 0.2 dex, while our numerical precision is around 0.05 dex.

With the 0.5 dex criterion per bin, we do not rule out the other respective models, but rather indicate that significant evidence exists that either an IL or double-NL (or both) is present in the data as a "population" in the filament overall, or conversely, that such an additional line emission is most likely not present if the evidence for the NL model relative to other models is larger than 0.5 dex.

The Bayes factors, and in turn the fraction of bins that show significant evidence for a double-NL (e.g. bin in Figure 5) or an additional IL (e.g. bin in Figure 2), depend on the line width used to distinguish an NL from an IL. Our choice of FWHM = 100 km s⁻¹ as the lower IL width limit corresponds to shock speeds as low as 1500 km s⁻¹ while previous results show the shock speed of Tycho's SNR to be at least 2000 km s⁻¹  (Ghavamian et al. 2001). W_IL < 100 km s⁻¹ requires unrealistic full electron−proton temperature equilibrium (Caprioli 2015) or shock speeds lower than 1500 km s⁻¹ and zero equilibration, i.e., ${T}_{e}/{T}_{p}={m}_{e}/{m}_{p}$. The latter case was already debated by Morlino et al. (2012), and furthermore, ${T}_{e}/{T}_{p}\ll 0.01$ was never measured in any of the remnants (Ghavamian et al. 2013). On the other hand, NL widths larger than 60 km s⁻¹ have not been observed before, and $\gt 100$ km s⁻¹ would require a CR acceleration efficiency >40% (Morlino et al. 2013), compared to the more realistic 10%–20% efficiency in SNR shocks (as was also found in Tycho's SNR by Morlino & Caprioli 2012). Therefore, our NL−IL separating line width limit robustly distinguishes the NL and IL that arise from different processes (cold neutral excitation and BN precursor).

Following the analysis described above, we calculated posterior parameter distributions for every model and Voronoi bin, which we then summarized by the median of the 1D-marginalized posteriors as a central parameter estimator, and the boundaries of the shortest (highest-density) 95% confidence interval. We also compared models using the Bayes factors (evidence ratios) calculated using the CV likelihood method. Our focus is on the measured NL width, evidence for and the magnitude of a split in the NL, evidence for an IL and its strength and width, the variation of LOS velocities across the filament, and possible correlations among the line parameters.

In Figure 2, we illustrate the application of our analysis to one of the bins from our data set. The observed spectrum (solid black line) is shown in the top blue panel, and the median NLIL model in solid red, while the model components are shown with dashed lines: the nonparametric background model (dashed black) and the parametrized "source" (SNR emission) model components—a constant component, one NL, and one IL—in dashed red. The background spectrum, which has been derived and fixed independently (see Appendix A) shows geocoronal and Galactic Hα emission, including the Galactic [N ii] line at around −130 km s⁻¹ LOS velocity. The 1D-marginalized posterior parameter distributions (solid black) and prior parameter distributions (dashed black) are shown in the remaining eight panels. The prior parameter distributions are overplotted for comparison with the posterior. The priors are chosen to not strongly prefer any parameter values over the allowed parameter range, but taper off smoothly toward zero at the range limits. Most posteriors are significantly different from the priors, and thus all of the parameters are well-constrained by the data. The parameters of the median model are denoted with red solid vertical lines. In this example, the median NL FWHM width is ${W}_{{NL}}\approx 40$ km s⁻¹. The IL is ${W}_{{IL}}\approx 210$ km s⁻¹ wide and comprises ≈40% of the total flux. With the red dashed vertical lines, we marked the boundaries of the 95% confidence interval. Compared to all other model parameters, the posterior shape of W_IL is more sensitive to the choice of prior (Appendix F). Still, the median and 95% confidence intervals of W_IL of a flat and Beta prior agree within ≈10%. This relatively benign difference occurs because by coincidence the peak of the likelihood is close to the center of the prior range, and because the 95% confidence interval boundaries are similar to the minimum/maximum limits of the prior range.

To summarize the parameter estimation of one bin, we consider the evidence-weighted 1D-posterior (see Figure 6, middle row), to which any of the four models that feature the parameter of interest contribute their marginalized posterior in proportion to their evidence. We then give the median of the evidence-weighted posterior, as well as the boundaries of its 95% confidence interval.

To summarize the results of all bins, we take three different routes.

• 1.  
  We consider the distribution of all 82 evidence-weighted posteriors' medians, as well as their 2 × 82 lower and upper 95% confidence interval boundaries (Figure 6, top panels). This yields information on the variability of line parameters across the filament.
• 2.  
  Second, we combine (average) the evidence-weighted posteriors of all bins, providing a representation of the information that we typically find in one individual bin (Figure 6, middle row of panels, black curves). We also combine model-specific posteriors separately to illustrate the relative contribution of different models and how parameter estimates depend on the model choice (colored curves).
• 3.  
  Finally, we evaluate the parameter constraints as imposed by all data (bins) combined. For that purpose, we sample from all bins' evidence-weighted posteriors, each time computing the median value over all bins. That is, we define one new parameter for each model parameter: the cross-bin median. Its posterior is plotted in the bottom row of Figure 6, and again summarized by its median and 95% confidence interval boundaries. We emphasize that this measure is decidedly distinct from modeling the spectrum of all bins combined, because it is still based on models constrained by all bins individually and independently, and, in particular, allows for local shifts in the line centroids to avoid artificial line broadening.

Apart from estimating model parameters, we compare the probabilities of models (different numbers of emission lines) against each other. In each bin, we define a model as favored if it has the highest evidence and a >0.5 dex logarithmic evidence ratio (>3:1 probability) over the NL model (see Figure 3, right panel). This threshold is a matter of choice; it is twice the 0.2 dex by which the correct model is typically favored in our tests on simulated data and reflects our approximate notion of the minimum for a "significant" probability. If no model satisfies the 0.5 dex requirement, we consider the fiducial single-NL model as favored. In this scheme, the data in one bin may favor an IL (NLIL or NLNLIL model), a double-NL (NLNL or NLNLIL model), both (NLNLIL), or none of them (NL). Using this criterion for each bin separately, we ascertain the fraction of bins in which the evidence indicates an IL, as well as the fraction of double-NL occurrence.

4.1. Narrow-line Width

4.2. Evidence for a Split in the NL

4.3. Intermediate-line Evidence and Parameters

As quantified by the Bayes factors, we find that 34% of the Voronoi bins (28 bins out of 82) are significantly better explained when a line is added to the fiducial NL-only model. In 74% of those, either the NLIL or NLNLIL model is also preferred over an NLNL model, which means 24% of the bins overall. This is illustrated in Figure 4, where the logarithmic (base-10) ratios of all models versus one another are shown in separate panels, and each panel shows the spatial variation of the corresponding ratio. These values are visualized by the histograms and also tabulated in the Tables 2 and 3. Apart from just a few bins, the NLNLIL model is not clearly favored over the NLIL model (panel 6). The situation changes when it comes to the comparison between NLNLIL and NLNL (panel 5), further confirming the necessity of an IL component. One example (bin 17 in Table 2) is illustrated in Figure 2, where the logarithmic evidence ratio of the NLIL model relative to the NL and NLNL models is larger than 1 dex. This means that an IL in addition to a single NL explains our data more than 10 times better than a simple NL model or NLNL model, irrespective of any particular parameter values. Both NLIL and NLNLIL models are equally likely for this bin, which implies no need for the second NL component. For this particular bin, we present posteriors for all other models in Appendix D.

Zoom In Zoom Out Reset image size

Table 2.  Model Comparison for the 73 Spatial Bins in the Eastern Shock Filament of Tycho's SNR

Bin	x ['']	y ['']	Pix	S/N	0L	NL	NLNL	NLIL	NLNLIL
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
1	22.0	59.2	20	14.40	27.073	0.108	0	0.077	0.234
2	21.4	58.4	26	14.21	31.122	0.307	0.273	0	0.098
3	22.4	58.2	21	14.96	40.522	0.245	0	0.236	0.447
4	22.8	59.7	26	12.52	13.405	0	0.318	0.220	0.476
5	21.6	57.4	22	12.33	27.723	0.435	0.343	0	0.110
6	21.0	60.7	46	13.08	31.293	0.459	0.391	0	0.091
7	19.6	64.2	188	12.04	39.630	0.483	0.325	0.270	0
8	21.1	59.5	26	15.19	30.855	0.176	0.113	0	0.160
9	18.3	61.1	378	12.84	37.771	0.376	0.348	0	0.185
10	20.5	67.3	234	11.43	33.046	0.410	0.146	0.473	0
11	20.4	57.1	51	12.56	25.860	0.236	0	0.370	0.295
12	21.3	56.3	32	11.07	25.137	0	0.217	0.074	0.266
13	22.4	57.2	22	11.97	17.165	0	0.161	0.246	0.356
14	19.7	58.5	104	13.81	42.871	0.594	0.366	0.088	0
15	22.5	55.9	48	10.68	14.121	0	0.182	0	0.141
16	21.3	55.0	40	10.47	19.664	0.188	0.178	0.272	0
17	19.5	55.3	137	13.05	30.044	1.456	1.307	0	0.027
18	22.4	54.5	27	9.69	15.533	0	0.203	0.013	0.121
19	21.4	53.8	37	10.44	12.876	1.387	0	0.475	0.203
20	20.2	53.4	51	11.33	17.289	0.094	0.360	0	0.126
21	23.4	55.0	40	8.49	9.724	0.436	0.084	0	0.092
22	21.3	52.6	31	10.51	20.602	0	0.109	0.050	0.207
23	22.5	53.3	30	10.98	19.092	0.518	0.212	0.232	0
24	18.5	52.7	138	11.15	23.991	0.107	0	0.333	0.303
25	22.3	52.1	27	9.17	4.052	0	0.249	0.132	0.333
26	20.4	51.8	38	11.57	10.662	0	0.267	0.281	0.434
27	21.6	51.4	23	9.86	18.144	0.016	0	0.081	0.119
28	19.4	51.0	48	11.10	21.513	0.018	0	0.456	0.444
29	20.8	50.8	19	10.63	20.281	0.166	0.223	0	0.146
30	23.9	53.6	67	9.73	9.985	0.524	0.488	0	0.024
31	19.3	49.6	36	10.09	23.207	0	0.131	0.041	0.254
32	20.3	50.0	28	11.16	18.678	0.459	0.628	0	0.152
33	20.9	49.0	32	10.58	16.697	0.218	0.105	0.010	0
34	19.9	48.7	26	11.08	21.270	0.081	0.098	0	0.099
35	18.7	48.1	47	12.01	23.668	0	0.085	0.084	0.265
36	20.2	47.5	48	11.38	21.533	0	0.008	0.052	0.165
37	17.5	49.6	165	11.17	19.644	0.316	0.352	0	0.112
38	19.0	46.9	37	11.92	25.193	0.730	0.622	0	0.070
39	18.6	46.0	21	12.10	23.762	0.603	0	0.834	0.452
40	20.3	45.5	115	12.72	21.540	1.005	0	0.258	0.033
41	17.3	44.3	59	11.28	23.175	0.503	0.713	0	0.072
42	18.9	44.9	35	11.55	23.385	0.731	0.432	0	0.084
43	18.5	44.2	16	9.10	11.179	1.116	0.790	0.065	0
44	16.4	46.4	279	12.59	27.058	0.125	0.066	0	0.082
45	18.4	43.4	31	11.30	21.451	0.033	0.072	0	0.173
46	16.8	42.6	98	12.74	22.116	1.082	0.673	0	0.020
47	18.7	42.3	52	11.63	20.054	0	0.303	0.172	0.404
48	16.9	40.6	58	12.11	26.149	0.207	0	0.051	0.041
49	17.9	41.3	29	11.62	18.944	0.195	0.282	0	0.118
50	18.0	39.7	42	12.24	20.337	0.324	0.289	0	0.112
51	15.5	35.6	58	11.36	28.214	0.685	0.432	0.459	0
52	17.4	38.1	21	10.47	17.933	0.324	0	0.113	0.045
53	16.9	39.0	38	12.27	22.999	0.265	0	0.404	0.107
54	18.5	38.0	58	10.56	12.943	0	0.022	0.292	0.268
55	16.9	37.3	29	12.04	18.049	0.466	0.152	0.195	0
56	17.9	36.7	29	11.74	18.795	0.198	0.316	0	0.151
57	16.8	36.2	23	12.53	16.550	0.589	0.619	0	0.083
58	19.5	33.6	344	11.89	16.533	0.187	0.395	0	0.117
59	17.8	35.7	34	12.40	30.850	0.695	0	0.361	0.150
60	16.6	35.2	19	11.81	16.934	0.091	0	0.096	0.119
61	17.2	33.5	25	11.15	17.640	0.033	0.203	0	0.118
62	17.3	34.7	28	10.86	19.534	0	0.038	0.172	0.264
63	16.4	34.2	25	11.63	26.490	0.648	0.649	0	0.110
64	17.6	32.1	42	12.86	26.277	2.415	1.670	0	0.098
65	16.5	32.8	29	13.09	22.555	0	0.123	0.128	0.278
66	16.5	31.7	34	12.12	27.519	0.086	0.126	0	0.154
67	16.0	30.7	32	12.02	17.226	0	0.131	0.088	0.203
68	14.4	31.7	155	12.61	32.667	0	0.148	0.165	0.393
69	13.2	22.1	59	9.33	18.851	0.095	0.241	0	0.171
70	12.3	18.4	285	10.04	15.999	0.047	0.152	0	0.112
71	13.8	26.8	47	8.78	8.551	0	0.239	0.030	0.194
72	15.3	33.6	61	11.76	19.728	0.090	0.377	0	0.185
73	13.3	20.6	59	9.43	11.914	0	0.208	0.138	0.300

Note. Columns 1–5: number of the (Voronoi) bin, x and y coordinates of the bin centroid, number of combined pixels, and signal-to-noise ratio. Columns 6–10: relative log CV likelihoods of the favored model (denoted with 0) to other models.

Download table as:  ASCIITypeset images: 1 2

Table 3.  Same as Table 2 just for the 9 bins in the Northern Filament

Bin	x ['']	y ['']	Pix	S/N	0L	NL	NLNL	NLIL	NLNLIL
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
1	13.0	23.7	351	14.30	33.932	0.846	0	0.667	0.357
2	10.4	20.6	62	15.31	43.600	0.483	0.479	0	0.264
3	8.9	22.1	169	15.36	53.524	1.645	0.423	0	0.009
4	12.3	20.8	169	14.77	41.332	0.003	0	0.068	0.205
5	8.6	19.7	129	15.50	49.913	1.194	0	0.644	0.517
6	11.7	17.3	355	16.72	34.586	0.470	0	0.186	0.077
7	9.0	17.8	75	16.31	42.528	0.959	0	0.741	0.269
8	9.1	15.0	190	14.65	40.034	1.090	0.479	0.329	0
9	7.4	16.5	104	14.73	37.961	0.086	0.246	0.025	0

Download table as:  ASCIITypeset image

Prominent ILs are seen across the entire filament (panel 3 in Figure 3), with an IL flux fraction of up to 42%. If we only consider bins that favor the NLIL or NLNLIL model, we get an IL flux fraction of 28% on average, and the intermediate-to-narrow flux fraction ${f}_{\mathrm{IL}}/{f}_{\mathrm{NL}}$ is estimated at 0.61 with a 95% confidence interval of (0.01–1.87).

Zoom In Zoom Out Reset image size

In contrast to the IL flux, W_IL is not well-constrained in individual bins, as indicated by the similarity of the prior and average bin-specific posterior (Figure 6, middle panels). However, several bins with high S/N and/or strong IL emission (as the bin in Figure 2) have good constraints on the IL width. Moreover, the combination of the data in all bins provides even better information: the cross-bin median is ${W}_{\mathrm{IL}}\in [166.77,194.15]$ km s⁻¹ (95% confidence; see the bottom of Figure 6). The global 95% confidence interval of ${f}_{\mathrm{IL}}/{f}_{\mathrm{NL}}$ is [0.34, 0.47]; however, we caution that it varies considerably between bins, with median $\mathrm{log}({f}_{\mathrm{IL}}/{f}_{\mathrm{NL}})\,\in [-0.79,0.05]$.

Zoom In Zoom Out Reset image size

Table 4.  Median and Highest Density 95% Confidence Interval of the Parameter Posteriors in the Favored Model of Voronoi Bins

Bin	Model	${W}_{\mathrm{NL}1}$	${W}_{\mathrm{NL}2}$	WIL	${f}_{\mathrm{IL}}/{f}_{\langle \mathrm{NL}\rangle }$	${\rm{\Delta }}{\mu }_{\mathrm{NL}}$
		[km s−1]	[km s−1]	[km s−1]		[km s−1]
1	NL	60.53	⋯	⋯	⋯	⋯
		[43.84, 80.68]
2	NL	65.12	⋯	⋯	⋯	⋯
		[48.10, 84.66]
3	NL	59.23	⋯	⋯	⋯	⋯
		[44.61, 76.77]
4	NL	72.07	⋯	⋯	⋯	⋯
		[49.95, 94.80]
5	NL	60.90	⋯	⋯	⋯	⋯
		[44.14, 81.94]
6	NLIL	40.63	⋯	191.02	0.51	⋯
		[24.78, 56.31]		[100.15, 308.96]	[0.01, 1.55]
7	NLNLIL	35.23	31.68	195.54	0.26	21.65
		[15.02, 63.35]	[15.11, 79.41]	[102.05, 310.43]	[0, 0.64]	[0.04, 108.74]
8	NL	60.78	⋯	⋯	⋯	⋯
		[45.34, 78.82]
9	NL	46.85	⋯	⋯	⋯	⋯
		[32.28, 61.53]
10	NL	40.87	⋯	⋯	⋯	⋯
		[28.68, 54.37]
11	NL	54.04	⋯	⋯	⋯	⋯
		[35.29, 73.60]
12	NL	78.70	⋯	⋯	⋯	⋯
		[62.54, 95.65]
13	NL	70.71	⋯	⋯	⋯	⋯
		[49.50, 92.40]
14	NLNLIL	40.34	35.18	191.06	0.33	22.74
		[15.76, 81.33]	[15.50, 72.99]	[102.59, 311.14]	[0, 0.80]	[0.16, 101.29]
15	NL	71.88	⋯	⋯	⋯	⋯
		[43.78, 99.09]
16	NL	75.86	⋯	⋯	⋯	⋯
		[56.59, 96.57]
17	NLIL	39.13	⋯	212.36	0.88	⋯
		[22.30, 56.71]		[100.18, 307.88]	[0.01, 6.15]
18	NL	82.43	⋯	⋯	⋯	⋯
		[64.85, 98.87]
19	NLNL	45.04	65.68	⋯	⋯	91.34
		[15.39, 90.16]	[39.58, 90.64]			[41.02, 115.24]
20	NL	74.45	⋯	⋯	⋯	⋯
		[54.71, 96.58]
21	NL	73.00	⋯	⋯	⋯	⋯
		[47.82, 98.64]
22	NL	68.76	⋯	⋯	⋯	⋯
		[50.01, 88.42]
23	NLNLIL	48.04	49.68	168.44	0.42	51.99
		[15.86, 88.08]	[19.25, 84.92]	[100.05, 295.69]	[0, 1.78]	[11.42, 153.20]
24	NL	59.97	⋯	⋯	⋯	⋯
		[43.21, 77.48]
25	NL	66.63	⋯	⋯	⋯	⋯
		[36.35, 96.93]
26	NL	70.97	⋯	⋯	⋯	⋯
		[47.59, 95.35]
27	NL	58.66	⋯	⋯	⋯	⋯
		[39.30, 81.32]
28	NL	56.87	⋯	⋯	⋯	⋯
		[41.50, 74.12]
29	NL	56.07	⋯	⋯	⋯	⋯
		[36.44, 79.31]
30	NLIL	59.03	⋯	181.53	1.44	⋯
		[21.68, 93.56]		[103.90, 288.14]	[0, 7.08]
31	NL	56.38	⋯	⋯	⋯	⋯
		[37.83, 74.36]
32	NLIL	41.53	⋯	150.76	1.02	⋯
		[15.31, 76.64]		[100.17, 278.80]	[0, 3.29]
33	NL	38.17	⋯	⋯	⋯	⋯
		[22.16, 56.19]
34	NL	47.04	⋯	⋯	⋯	⋯
		[30.20, 65.05]
35	NL	72.60	⋯	⋯	⋯	⋯
		[56.45, 89.57]
36	NL	59.26	⋯	⋯	⋯	⋯
		[40.97, 80.03]
37	NL	42.92	⋯	⋯	⋯	⋯
		[23.91, 62.21]
38	NLIL	41.82	⋯	170.06	0.74	⋯
		[21.28, 66.79]		[100.12, 294.73]	[0, 2.26]
39	NLNL	51.41	56.03	⋯	⋯	60.15
		[16.46, 90.57]	[29.12, 80.18]			[9.50, 84.17]
40	NLNL	38.06	63.21	⋯	⋯	39.90
		[15.61, 61.48]	[26.59, 96.32]			[1.32, 144.16]
41	NLIL	48.34	⋯	201.48	0.62	⋯
		[24.99, 70.31]		[103.46, 310.50]	[0.03, 1.67]
42	NLIL	48.09	⋯	162.60	0.61	⋯
		[27.07, 72.44]		[100.39, 291.90]	[0.01, 1.82]
43	NLNLIL	44.57	33.42	171.73	0.75	55.59
		[15.43, 87.20]	[15.03, 77.23]	[100.31, 295.94]	[0.01, 2.23]	[1.57, 97.03]
44	NL	65.49	⋯	⋯	⋯	⋯
		[48.42, 85.13]
45	NL	55.89	⋯	⋯	⋯	⋯
		[38.00, 76.15]
46	NLIL	35.29	⋯	221.61	0.89	⋯
		[17.55, 53.39]		[100.39, 317.73]	[0, 10.90]
47	NL	71.18	⋯	⋯	⋯	⋯
		[51.66, 91.45]
48	NL	49.41	⋯	⋯	⋯	⋯
		[32.16, 66.58]
49	NL	43.21	⋯	⋯	⋯	⋯
		[27.87, 63.95]
50	NL	61.30	⋯	⋯	⋯	⋯
		[42.56, 83.35]
51	NLNLIL	50.84	37.90	201.87	0.31	37.67
		[16.21, 81.76]	[16.10, 66.80]	[104.16, 313.03]	[0, 0.89]	[1.03, 90.84]
52	NL	61.42	⋯	⋯	⋯	⋯
		[39.35, 85.21]
53	NL	54.00	⋯	⋯	⋯	⋯
		[38.21, 71.73]
54	NL	47.73	⋯	⋯	⋯	⋯
		[23.04, 68.53]
55	NLNLIL	33.67	50.83	187.38	0.50	34.31
		[15.24, 72.78]	[17.70, 88.41]	[101.07, 312.58]	[0, 1.49]	[1.00, 82.17]
56	NL	58.08	⋯	⋯	⋯	⋯
		[38.88, 80.69]
57	NLIL	45.30	⋯	199.89	1.12	⋯
		[17.24, 72.11]		[107.28, 312.43]	[0.02, 3.11]
58	NL	74.99	⋯	⋯	⋯	⋯
		[54.48, 97.77]
59	NLNL	57.49	58.12	⋯	⋯	63.14
		[19.54, 93.95]	[36.03, 82.85]			[11.53, 93.23]
60	NL	66.49	⋯	⋯	⋯	⋯
		[46.53, 89.49]
61	NL	73.42	⋯	⋯	⋯	⋯
		[52.04, 96.65]
62	NL	55.28	⋯	⋯	⋯	⋯
		[37.46, 75.62]
63	NLIL	49.93	⋯	202.95	0.49	⋯
		[32.62, 69.52]		[106.47, 316.80]	[0, 1.34]
64	NLIL	38.86	⋯	191.52	1.12	⋯
		[19.78, 60.12]		[107.16, 298.62]	[0.12, 2.53]
65	NL	68.07	⋯	⋯	⋯	⋯
		[47.10, 92.23]
66	NL	57.62	⋯	⋯	⋯	⋯
		[42.10, 77.46]
67	NL	65.35	⋯	⋯	⋯	⋯
		[46.28, 87.53]
68	NL	72.89	⋯	⋯	⋯	⋯
		[57.12, 89.85]
69	NL	60.22	⋯	⋯	⋯	⋯
		[41.25, 82.71]
70	NL	46.09	⋯	⋯	⋯	⋯
		[27.58, 67.10]
71	NL	50.96	⋯	⋯	⋯	⋯
		[29.85, 76.82]
72	NL	78.28	⋯	⋯	⋯	⋯
		[61.10, 96.81]
73	NL	64.24	⋯	⋯	⋯	⋯
		[42.79, 89.21]

Note. Shown here are the spatial bins in the eastern filament of Tycho's SNR. If none of the multi-line models is at least ≈3 times (or 0.5 dex) more likely, the single-line (NL) model is taken to be the favored one.

Download table as:  ASCIITypeset images: 1 2 3

Table 5.  Same as Table 4 Just for Voronoi Bins in Tycho's Northern Filament

Bin	Model	${W}_{\mathrm{NL}1}$	${W}_{\mathrm{NL}2}$	WIL	${f}_{\mathrm{IL}}/{f}_{\langle \mathrm{NL}\rangle }$	${\rm{\Delta }}{\mu }_{\mathrm{NL}}$
		[km s−1]	[km s−1]	[km s−1]		[km s−1]
1	NLNL	33.85	53.16	⋯	⋯	28.04
		[15.04, 61.64]	[21.03, 87.55]			[0.29, 156.41]
2	NLIL	58.50	⋯	195.28	0.39	⋯
		[41.77, 75.74]		[102.92, 316.25]	[0, 1.13]
3	NLIL	54.32	⋯	156.07	0.41	⋯
		[37.83, 69.52]		[100.12, 274.97]	[0.02, 1.08]
4	NL	59.63	⋯	⋯	⋯	⋯
		[37.09, 78.26]
5	NLNL	51.95	69.86	⋯	⋯	39.81
		[25.10, 73.85]	[34.12, 97.51]			[2.97, 173.42]
6	NLNL	42.30	56.58	⋯	⋯	20.37
		[15.65, 83.53]	[20.46, 91.93]			[0.05, 67.38]
7	NLNL	37.89	60.23	⋯	⋯	42.00
		[15.09, 66.94]	[26.97, 93.47]			[6.73, 184.78]
8	NLNLIL	44.66	34.55	149.45	0.36	29.63
		[15.68, 70.78]	[15.22, 80.12]	[100.36, 291.76]	[0, 0.89]	[0.23, 86.09]
9	NL	56.38	⋯	⋯	⋯	⋯
		[37.58, 74.56]

Download table as:  ASCIITypeset image

4.4. Line-of-sight Velocity

We present the observed LOS velocity, i.e., the NL centroid ${\mu }_{\mathrm{NL}}$, in Figure 7. The map of the bin-specific ${\mu }_{\mathrm{NL}}$ is shown in the left panels. In the top-right panel, we show the corresponding histogram (green), as well as the distribution of the 95% confidence interval boundaries. Across the entire NE rim, we find median LOS velocities in the range $[-67.9,-28.6]$ km s⁻¹. Notably, the northern filament moves with respect to the eastern filament, ${\mu }_{\mathrm{NL}}=-32.9$ km s⁻¹ and −40.5 km s⁻¹ on average. The bin-to-bin variations in ${\mu }_{\mathrm{NL}}$ are most noticeable in the eastern filament, indicating inhomogeneities in the shock and, in turn, in the ISM density. The bulk (median) LOS velocity of the filament as a whole is ${\mu }_{\mathrm{NL}}=-38.3\pm 1.5$) km s⁻¹ at the 95% confidence level (see the middle-right panel). Converting to the local standard of rest (LSR), we obtain ${V}_{\mathrm{LSR}}\approx -34$ km s⁻¹. We also check for the correlation between the LOS velocities and surface brightness (bottom-right panel), but find them to be uncorrelated.

Zoom In Zoom Out Reset image size

4.5. "Knot g"

4.6. Correlation in Parameters

We investigate the correlation between various line parameters and surface brightness, as well as correlation among model parameters, using Spearman's rank correlation coefficient (ρ). It is sensitive to any monotonic relationship, even if it is not linear as for Pearson's correlation coefficient, and does not require normally distributed data. The result indicates a strong correlation (or anticorrelation) for values close to ±  1.

In the top row in Figure 8, we plotted the median of evidence-weighted W_NL, ${f}_{\mathrm{IL}}/{f}_{\mathrm{NL}}$, W_IL, and ${\rm{\Delta }}{\mu }_{\mathrm{NL}}$ posteriors against the surface brightness (SB; total flux of the intrinsic model divided by the bin area). We separately test the correlation for bins for which a single-NL model (NL or NLIL, marked with squares) has the highest evidence, and those bins for which a double-NL model (NLNL or NLNLIL, marked with circles) has the highest evidence of all models. We do not find any strong correlation: parameters of single-NL models are uncorrelated with surface brightness, while double-NL models show weak correlations with SB. The highest correlation is with ${\rm{\Delta }}{\mu }_{\mathrm{NL}}$ (ρ = 0.57), meaning that with increasing SB, the two NL centroids tend to be more separated.

Zoom In Zoom Out Reset image size

In the bottom-row panels of the same figure, we estimated the correlation among the parameters. Although we see a hint of a positive correlation between the NL width and NL separation (ρ = 0.47 and ρ = 0.49; see panel 8), we do not find a clear correlation between any of the parameters. We get very similar results when we apply Pearson's correlation: its coefficient is either very close to Spearman's rank coefficient or closer to zero. In addition, the typical uncertainty of all of these points (plotted in the top-right corner of each panel and given as the average shortest 68% confidence interval) propagates into the correlation rank coefficient uncertainty.

Previous Hα observations of Tycho's "knot g" (Lee et al. 2007) were modeled by Wagner et al. (2009) by computing a series of time-dependent numerical simulations of CR-modified shocks. Assuming a distance of 2.1 kpc to the remnant, they found the CR diffusion coefficient of κ = 2 × 10²⁴ cm² s⁻¹ and the lower limit of the injection parameter ${\xi }_{\mathrm{inj}}=4.2\times {10}^{-3}$ to be in good agreement with the observations, suggesting that CR acceleration in the shock is efficient. Diffuse emission 1'' (∼10¹⁶ cm) ahead of the eastern filament was also detected by Lee et al. (2010) and interpreted as emission from the CR precursor with T ∼ 80,000–100,000 K.

In what follows, we will summarize the main results of our paper and look for theoretical explanations for our findings.

Presence of a CR precursor. The NL width is much broader than 20 km s⁻¹  (≈55 km s⁻¹ on average) in the entire NE rim regardless of whether the shock emission is described with a single or double NL. In other words, even when differential velocities (double NL) are present and accounted for, the NL is still significantly broadened. This clearly points toward gas heating in a CR precursor (Morlino et al. 2013). Furthermore, momentum transfer in a CR precursor might result in a split in the NL if we have two inclined shocks projected on our LOS. This is something that we observe—more precisely, we find significant Bayesian evidence for it in 18% of the data. Apart from the separation between the two NLs being 38 km s⁻¹ on average, at the same time we find that their intrinsic widths are around 49 km s⁻¹. If we assume that we have the contribution of two shocks inclined with the same angle to the LOS and that neutrals in the CR precursor acquired 10% of the previously estimated shock velocity of 2500 km s⁻¹  (Ghavamian et al. 2001), we find that the shock normal inclination of 85°–86° explains the centroid separation of 38 km s⁻¹. Obviously, for the same shock speed and smaller acquired bulk velocity, we would need the shocks slightly more inclined, i.e., with angles $\lt 85^\circ$. However, since only few bins show evidence for a split in the NL, we do not have enough information to construct the geometry of the shock.

Presence of a BN precursor. The main signature of the BN precursor is an IL (Morlino et al. 2012); 24% of the bins demand an additional line being specifically IL. We find a median value for the intrinsic IL width of 180 km s⁻¹, comprising on average a 41% intermediate-to-narrow flux ratio. Moreover, the observed high pre-shock neutral fraction of 0.9 (Ghavamian et al. 2001) in combination with the shock velocity of ∼2500 km s⁻¹ supports the BN presence since a large neutral fraction and the specified shock velocity contribute to efficient CE and the larger number of created BNs play an essential role in the formation of and heating in the BN precursor. BNs are ionized inside the BN precursor almost immediately when they cross the shock. The newly created protons move with a bulk speed larger than the Alfvén speed, hence they can trigger the streaming instability (and possibly other kind of instabilities), resulting in an increase of the ion temperature in the BN precursor. In turn, the CE between the pre-shock neutrals and warm ions creates warm neutrals that produce the IL.

For illustrative purposes, in Figure 9 we report the temperature profile of the pre-shock gas for specific values of shock parameters as calculated in Morlino et al. (2013) where the CR and BN precursors can be clearly distinguished. The BN precursor acts on scales of 10¹⁶–10¹⁷ cm in the immediate pre-shock region and can heat the gas to a temperature of ∼10⁶ K. The CR precursor length is much larger, depends on the maximum energy of the accelerated particles, and extends over several 10¹⁷ cm. The expected gas temperature in the CR precursor is ∼10⁵ K.

Zoom In Zoom Out Reset image size

Collisionless shock model prediction in partially ionized medium. Following the results of Morlino et al. (2013) and their Figure 9, 40–50 km s⁻¹ NL widths require efficient turbulent heating ${\eta }_{\mathrm{TH}}$ in the CR precursor, but also accelerated particles with the maximum momentum of p_max = 40 TeV/c or higher. However, this result is dependent on the shock velocity assumption, which is 4000 km s⁻¹ in the cited paper. Furthermore, Figure 10 in the same paper shows the IL width as a function of the CR acceleration efficiency for a fixed shock speed ${V}_{\mathrm{sh}}=4000$ km s⁻¹, ${\eta }_{\mathrm{TH}}=0.5$, and p_max =50 TeV/c, and two values of downstream electron-to-proton temperature ratios ${\beta }_{\mathrm{down}}$ being either 0.01 or 1. Interestingly, our measured IL width of 180 km s⁻¹ on average can be explained with the mentioned shock parameters and the acceleration efficiency of 15%–25%, although any further constraints are difficult since we do not know ${\beta }_{\mathrm{down}}$. The latter value, as well as the shock velocity, is something that can be constrained from the observations of the broad Hα component.

Zoom In Zoom Out Reset image size

In order to obtain the measured ${W}_{\mathrm{NL}}\approx 40$ km s⁻¹, one requires a combination of shock parameters (and CR acceleration included): p_max = 10 TeV/c, V_sh = 2500 km s⁻¹, and ${\beta }_{\mathrm{down}}$ of up to 0.1 (Morlino et al. 2013). At the same time, this configuration predicts an IL width of around 300 km s⁻¹ and intermediate-to-narrow flux ratio of 1.6 on average for various values of ${\eta }_{\mathrm{TH}}$ and . Both of these IL parameters are much higher than the median values we infer from our data and analysis. We speculate that somewhat higher ${\beta }_{\mathrm{down}}$ than 0.1 (see Figure 13 in Morlino et al. 2012) and a shock velocity of around 3000 km s⁻¹ would possibly be able to explain the observed ${W}_{\mathrm{IL}}\approx 180$ km s⁻¹, ${f}_{\mathrm{IL}}/{f}_{\mathrm{NL}}\approx 0.41$, and ${W}_{\mathrm{NL}}\approx 55$ km s⁻¹. Finally, we also notice that keeping ${\eta }_{\mathrm{TH}}$ constant and increasing , the model of Morlino et al. (2013) predicts a simultaneous increase in W_NL and decrease in W_IL and ${f}_{\mathrm{IL}}/{f}_{\mathrm{NL}}$. We find a hint for such an anticorrelation between W_NL and W_IL (panel 7 in our Figure 8).

Average and variation of velocities. Our result for the median LSR-corrected LOS velocity, ${V}_{\mathrm{LSR}}=(-34\pm 2)$ km s⁻¹, is in agreement with the earlier study by Lee et al. (2007), who reported an NL Hα LSR velocity −30.3 km s⁻¹ of "knot g." Furthermore, Ghavamian et al. (2017) reported a value of −45.6 km s⁻¹, which falls within our range [−64, −25] km s⁻¹ of observed LSR velocities in the NE rim. In H i 21 cm observations toward Tycho's NE rim, Reynoso et al. (1999) found V_LSR = −51.5 km s⁻¹ and associated the location of Tycho's SNR and the H i cloud with the Perseus arm. Similarly, ¹²CO emission was found at ${V}_{\mathrm{LSR}}=-62.5$ km s⁻¹ (Lee et al. 2004; see also Zhou et al. 2016) and might be associated with SN 1572's pre-shock gas. However, as pointed out by Tian & Leahy (2011), there is no clear evidence that either the H i cloud or CO cloud is physically associated with Tycho's NE rim. We leave further discussion of the V_LSR result and its possible interpretation to future work, where we will also investigate the shock and CR properties in more detail (see Section 6).

We present Hα spectroscopic observations of Tycho's NE Balmer filaments. This study provides spectroscopic data that, for the first time, are spatially resolved (spectro-imagery), with large coverage that comprises and resolves the entire NE filament. Our analysis is based on Bayesian inference, which enables a quantitative, probabilistic, and well-defined model comparison, and a comprehensive, complete characterization of the parameter probabilities.

We find that the broadening of the NL beyond 20 km s⁻¹ that was noted in previous studies was not an artifact of the spatial integration, and that it extends across the whole filament, not only the previously covered "knot g." The NL width in the NE rim is typically found around 55 km s⁻¹. Such a large width cannot be due to the superposition of multiple lines. In fact, our data analysis allows us to take projection effects into account when interpreting the data. We are able to distinguish between single-NL and double-NL models where we find significant evidence for a split in NL in 18% of the Voronoi bins. The widths of the two NLs are around 49 km s⁻¹ and their centroid separations are 38 km s⁻¹ on average.

An NL width of 55 km s⁻¹ implies a temperature of the upstream gas of ≈68,000 K. If this were the temperature of the unperturbed ISM where the SNR is expanding, no neutral hydrogen would exist in the first place, contradicting the presence of Hα emission. Hence, our finding is the signature of the existence of a mechanism able to heat the upstream plasma in a region ahead of the shock much smaller than the collisional ionization length scale. As shown by the previous study of Morlino et al. (2013), a CR precursor is the best candidate to explain the widening of the NL, opening the possibility of studying particle acceleration in shocks using Hα emission. The fact that the NL width ranges from 35 to 72 km s⁻¹ across the NE rim suggests that the amount of neutrals in the ambient medium varies, which imposes different degrees of damping of magnetic waves excited by CR streaming.

Likewise, we confirm the suspected presence of an IL and show it to be widespread (24% of the bins). Typical IL widths and intermediate-to-narrow flux ratios are 180 km s⁻¹ and 0.41, respectively.

Our model parameters also comprise the LOS velocity centroids. After correction to the local standard of rest, their median is ${V}_{\mathrm{LSR}}=(-34\pm 2)$ km s⁻¹, in agreement with the Lee et al. (2007) investigation of "knot g."

Overall, our results reveal an interplay between two precursors in Tycho's NE rim: broadened NL widths point toward the evidence for the presence of a CR precursor, while detected IL reveals the presence of a BN precursor.

From the knowledge of the NL width only it is not possible to determine the CR acceleration efficiency, because such width depends on many parameters (shock speed, maximum energy of accelerated particles, electron−ion equilibration, and turbulent heating). Nevertheless, we can conclude that, assuming a shock speed between 2500 and 3000 km s⁻¹, our result is compatible with having a maximum CR energy >10 TeV, a turbulent heating >10%, and an acceleration efficiency > few %. The degeneracy between these parameters could be broken using other information coming from the broad-line width and intensity (giving more precise information on shock speed and electron−ion equilibration) and X-ray/gamma-ray observations (determining the maximum energy of accelerated particles). The difficulties in performing such calculations rely on the fact that the broad line is not known with the same accuracy as the NL and IL and that the gamma-ray emission does not have enough spatial and spectral resolution to fix unambiguously the maximum energy (Park 2015; Morlino & Blasi 2016). Improvements in this regards will surely come from the Cherenkov Telescope Array. Furthermore, parallel to the study of the GHαFaS narrow Hα-line profiles, we have conducted an investigation on the same part of the Tycho's remnant using OSIRIS (Optical System for Imaging and low-intermediate Resolution Spectroscopy) on the Gran Telescopio Canarias to observe the broad Hα-line profiles. In a forthcoming paper, we will present results of broad Hα components of the same spatial locations (bins) along the filaments as presented in this paper, which, combined with narrow Hα components and applied shock models, will give a better handle on the overall conditions in the shock and will enable us to quantify CR properties.

We thank the anonymous referee for the constructive report and a valuable contribution toward the completeness of this paper. We also thank René Andrae (MPIA, Heidelberg) and Joonas Nättilä (Tuorla Observatory, University of Turku) for fruitful discussions and helpful suggestions on Bayesian inference. We thank Alex Borlaff (IAC, ULL) for his help with the observations.

In order to compare our models of Tycho's SNR spectra to the data, we need to account for the variable sensitivity ("flatfield") and the background flux. We eschew the standard method of applying background subtraction and flatfield correction to compensate for both effects in the individual exposures. Instead, we reconstruct the flatfield and background in the final product of the data reduction pipeline, which is a cube (position, wavelength) of co-added, "stacked" individual observations. Notably, the flatfield will also have a wavelength dependence in addition to being position dependent. We eventually include the flatfield and background in our models. The advantage of our approach is the preservation of photon statistics and therefore accurate uncertainties of the model parameters and evidence.

Before describing the construction of the individual as well as the co-added flatfield and background, we formalize the process through which the data are generated and how they are propagated by the data reduction pipeline.

A.1. Measurement Process

Each datum (measurement) ${D}_{{xy},i}$ with pixel indices (x,y) and exposure index i is the response of the telescope and instrument to the incoming, seeing-convolved flux ${d}_{0}={d}_{0}(\alpha ,\delta ,\lambda )$, which varies with sky coordinates $(\alpha ,\delta )$ and wavelength λ:

$$\begin{eqnarray}&&{D}_{{xy},i}={\int }_{-\infty }^{\infty }{d}_{0,i}(x,y,\lambda )\,{F}_{{xy}}\,{R}_{{xy},i}(\lambda )\,d\lambda \,,\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray*}&&{d}_{0,i}(x,y,\lambda )\equiv {d}_{0}(\alpha ={\alpha }_{i}(x,y),\delta ={\delta }_{i}(x,y),\lambda )\end{eqnarray*}$$

has the astrometric solution $\alpha ={\alpha }_{i}(x,y),\delta ={\delta }_{i}(x,y)$, and the integral over the area of pixel (x,y) implicitly applied.

F_xy is the spatially varying sensitivity of the detector and optical system—the flatfield image. As expected, and as verified by us (see below), it is the same for all exposures.

${R}_{{xy},i}(\lambda )$ is the pixel- and exposure-specific wavelength filter imposed by the Fabry–Pérot interferometer (etalon). In our case, it is accurately described by a universal line-spread function (LSF), $r(\lambda )$:

$$\begin{eqnarray}&&{R}_{{xy},i}(\lambda )=r({{\rm{\Lambda }}}_{{xy},i}-\lambda )\,,\end{eqnarray} \tag{ 4 }$$

where the wavelength calibration ${{\rm{\Lambda }}}_{{xy},i}$ returns the central, maximum-throughput wavelength for each pixel and exposure (tuning of the etalon). By definition, the LSF is centered on (peaks at) the origin. In our case, it is a Gaussian with dispersion measured from calibration spectra.

Inserting Equation (4) into Equation (3), one sees that the datum ${D}_{{xy},i}$ is the LSF-convolved spectral flux ${d}_{i}(x,y,\lambda$), evaluated at ${{\rm{\Lambda }}}_{{xy},i}\,$:

$$\begin{eqnarray}{D}_{{xy},i} & = & {\displaystyle \int }_{-\infty }^{\infty }{d}_{0,i}(x,y,\lambda )\,{F}_{{xy}}\,r({{\rm{\Lambda }}}_{{xy},i}-\lambda )\,d\lambda \\ & = & {F}_{{xy}}\cdot ({d}_{0,i}{\ast }_{\lambda }r)(x,y,\lambda ={{\rm{\Lambda }}}_{{xy},i})\\ & \equiv & {F}_{{xy}}\cdot {d}_{i}(x,y,{{\rm{\Lambda }}}_{{xy},i}).\end{eqnarray} \tag{ 5 }$$

A.2. Data Processing

By way of converting the observed images into data subcubes, ${D}_{{xy},i}\to {D}_{{xyl},i}$, we make the wavelength information contained in them explicit and obtain a format that allows direct co-addition. Each of the 48 "slices" (third index) of a subcube corresponds to a wavelength, ${\lambda }_{l}$. Each observed ${{\rm{\Lambda }}}_{{xy},i}$ is bracketed by two slices, and the corresponding flux ${D}_{{xy},i}$ is assigned to them via linear interpolation:

$$\begin{eqnarray}{D}_{{xyl},i} & = & {D}_{{xy},i}\cdot {T}_{{xyl},i}\\ {T}_{{xyl},i} & = & {t}_{l}({{\rm{\Lambda }}}_{{xy},i})\\ & \equiv & \,\left\{\begin{array}{rcl}1-| {{\rm{\Lambda }}}_{{xy},i}-{\lambda }_{l}| /{\rm{\Delta }}\lambda & : & | {{\rm{\Lambda }}}_{{xy},i}-{\lambda }_{l}| \leqslant {\rm{\Delta }}\lambda \\ 0 & : & | {{\rm{\Lambda }}}_{{xy},i}-{\lambda }_{l}| \gt {\rm{\Delta }}\lambda \end{array}\right.\,,\end{eqnarray} \tag{ 6 }$$

where ${\rm{\Delta }}\lambda$ is the "size" of each slice, i.e., the distance between the slices' central wavelengths, ${\lambda }_{l}$. Another way to describe this assignment is that each slice imposes a triangle filter ${t}_{l}(\lambda )=\max (0,1-| \lambda -{\lambda }_{l}| /{\rm{\Delta }}\lambda )$. We will use it again to construct the co-added flatfield and background cubes.

As ${D}_{{xy},i}$ itself is a filtered version of the flux, the subcube ${D}_{{xyl},i}$ at a given pixel (x, y) is nearly "empty," except for two slices l. Adequate coverage of the spectrum therefore necessitates multiple tunings of the etalon, or imaging the source in varying locations on the detector, since even for unchanged tuning, Λ depends on x and y. The resulting multiple subcubes are then co-added, albeit after projecting and spatially resampling them onto a common astrometrically calibrated frame $(x^{\prime} ,y^{\prime} )$:

$$\begin{eqnarray}&&{D}_{{xyl},i}\mathop{\longrightarrow }\limits^{\alpha ,\delta }{D}_{x^{\prime} y^{\prime} l,i}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{D}_{x^{\prime} y^{\prime} l}=\displaystyle \sum _{i}{D}_{x^{\prime} y^{\prime} l,i}.\end{eqnarray} \tag{ 8 }$$

${D}_{x^{\prime} y^{\prime} l}$ is the final, co-added data cube which, apart from spatial binning, directly constrains the SNR shock models. We account for the spatial sensitivity variations and any residual non-constant spectral sampling rate by including them in the models, in the form of the co-added flatfield cube, ${F}_{x^{\prime} y^{\prime} l}$, which we derive below along with the background cube, ${B}_{x^{\prime} y^{\prime} l}$.

A.3. Propagation of Flatfield and Background

We represent the incoming flux ${d}_{0,i}$ as the sum of the "source" flux (${s}_{0,i}$) that comprises celestial objects, in particular the SNR filament that we are interested in, and the background flux (${b}_{0,i}$). Because of the linearity of the integral in Equation (5), the same applies to the LSF-convolved flux d_i:

$$\begin{eqnarray*}{d}_{0,i}(x,y,\lambda ) & = & {s}_{0,i}(x,y,\lambda )+{b}_{0,i}(x,y,\lambda )\\ {s}_{i}(x,y,\lambda ) & = & ({s}_{0,i}\,{\ast }_{\lambda }r)(x,y,\lambda )\\ {b}_{i}(x,y,\lambda ) & = & ({b}_{0,i}\,{\ast }_{\lambda }r)(x,y,\lambda )\\ \mathop{\Longrightarrow }\limits^{(3)}\,{d}_{i}(x,y,\lambda ) & = & {s}_{i}(x,y,\lambda )+{b}_{i}(x,y,\lambda )\\ {D}_{{xy},i} & = & {F}_{{xy}}\cdot [{s}_{i}(x,y,{{\rm{\Lambda }}}_{{xy},i})+{b}_{i}(x,y,{{\rm{\Lambda }}}_{{xy},i})]\\ & \equiv & {F}_{{xy}}\cdot {S}_{{xy},i}+{B}_{{xy},i}\,.\end{eqnarray*}$$

Note that in contrast to the source component of the data,

$$\begin{eqnarray}&&{S}_{{xy},i}={s}_{i}(x,y,{{\rm{\Lambda }}}_{{xy},i}),\end{eqnarray} \tag{ 9 }$$

we have absorbed the flatfield in the definition of the data's background component,

$$\begin{eqnarray}&&{B}_{{xy},i}={F}_{{xy}}\cdot {b}_{i}(x,y,{{\rm{\Lambda }}}_{{xy},i}).\end{eqnarray} \tag{ 10 }$$

Due to the linearity of the wavelength assignment (6),

$$\begin{eqnarray*}{D}_{{xyl},i} & = & {T}_{{xyl},i}\,{F}_{{xy}}\,{S}_{{xy},i}+{T}_{{xyl},i}\,{B}_{{xy},i}\\ & \equiv & {F}_{{xyl},i}\,{S}_{{xy},i}+{B}_{{xyl},i}.\end{eqnarray*}$$

Here we have defined the background subcube ${B}_{{xyl},i}={T}_{{xyl},i}\,{B}_{{xy},i}$, and the flatfield subcube ${F}_{{xyl},i}={T}_{{xyl},i}{F}_{{xy}}$, which are the background components of the data and the flatfield image (the response to a flat spectrum), respectively, transformed to a subcube using the same slice-assignment (6) that was used for the data itself. Also, as for the data, background and flatfield subcubes are projected onto the same frame and co-added:

$$\begin{eqnarray*}{F}_{{xyl},i} & \mathop{\to }\limits^{\alpha ,\delta } & {F}_{x^{\prime} y^{\prime} l,i}\\ {S}_{{xy},i} & \mathop{\to }\limits^{\alpha ,\delta } & {S}_{x^{\prime} y^{\prime} ,i}\\ {B}_{{xyl},i} & \mathop{\to }\limits^{\alpha ,\delta } & {B}_{x^{\prime} y^{\prime} l,i}\\ {D}_{x^{\prime} y^{\prime} l} & = & \displaystyle \sum _{i}({F}_{x^{\prime} y^{\prime} l,i}\,{S}_{x^{\prime} y^{\prime} ,i}+{B}_{x^{\prime} y^{\prime} l,i})\end{eqnarray*}$$

${S}_{x^{\prime} y^{\prime} ,i}={s}_{i}(x^{\prime} ,y^{\prime} ,{{\rm{\Lambda }}}_{x^{\prime} y^{\prime} ,i})$ is ${s}_{i}(x,y,{{\rm{\Lambda }}}_{{xy},i})$ with $x={x}_{i}(x^{\prime} ),y\,={y}_{i}(y^{\prime} )$ given by the astrometric solution. The same transformation applies to b. Since the co-added frame's $(x^{\prime} ,y^{\prime} )$ corresponds to a unique sky position, and we assume the filament emission (but not necessarily the background!) to be constant between exposures, ${s}_{i}(x^{\prime} ,y^{\prime} ,\lambda )=s(x^{\prime} ,y^{\prime} ,\lambda )$. Then,

$$\begin{eqnarray}{D}_{x^{\prime} y^{\prime} l} & = & \displaystyle \sum _{i}{F}_{x^{\prime} y^{\prime} l,i}\,s(x^{\prime} ,y^{\prime} ,{{\rm{\Lambda }}}_{x^{\prime} y^{\prime} ,i})\\ & & +{F}_{x^{\prime} y^{\prime} l,i}\,{b}_{i}(x^{\prime} ,y^{\prime} ,{{\rm{\Lambda }}}_{x^{\prime} y^{\prime} ,i})\\ & \approx & \left(\displaystyle \sum _{i}{F}_{x^{\prime} y^{\prime} l,i}\right)\,s(x^{\prime} ,y^{\prime} ,{\lambda }_{l})\\ & & +\displaystyle \sum _{i}{F}_{x^{\prime} y^{\prime} l,i}\,{b}_{i}(x^{\prime} ,y^{\prime} ,{\lambda }_{l})\,\\ & \equiv & {F}_{x^{\prime} y^{\prime} l}\cdot s(x^{\prime} ,y^{\prime} ,{\lambda }_{l})+{B}_{x^{\prime} y^{\prime} l}.\end{eqnarray} \tag{ 11 }$$

In the second line, we took ${F}_{x^{\prime} y^{\prime} l,i}$ to be non-zero only for l nearest to the sampled wavelength, ${\lambda }_{l}\approx {{\rm{\Lambda }}}_{x^{\prime} y^{\prime} ,i}$. Hence, ${F}_{x^{\prime} y^{\prime} l,i}\cdot s(x^{\prime} ,y^{\prime} ,{{\rm{\Lambda }}}_{x^{\prime} y^{\prime} ,i})\approx {F}_{x^{\prime} y^{\prime} l,i}\cdot s(x^{\prime} ,y^{\prime} ,{\lambda }_{l})$ for all $l=1,\ldots ,\,48$, and analogously for b. However, the background does not factor out of the co-addition, as it may change between exposures. In the last step, we have defined the flatfield and background cubes ${F}_{x^{\prime} y^{\prime} l}$ and ${B}_{x^{\prime} y^{\prime} l}$ as the sum of the flatfield and background subcubes.

A.4. Background and Flatfield Model

In order to isolate the source flux component of the data, we still need to model the flatfield and background.

We begin by investigating the background on the co-added data (8), (11). Using Sextractor, we mask sources (stars, galaxies, and the filament itself) on the wavelength-integrated frame ${D}_{x^{\prime} y^{\prime} }={\sum }_{l=1}^{48}{D}_{x^{\prime} y^{\prime} l}$. This frame is the deepest, highest S/N data product we have available, hence the mask is as complete as possible. We obtain the set of unmasked pixels, $({x}_{u}^{{\prime} },{y}_{u}^{{\prime} })$, which we can use to measure the background without "contamination" by celestial sources:

$$\begin{eqnarray*}&&{D}_{{x}_{u}^{{\prime} },{y}_{u}^{{\prime} },l}=\displaystyle \sum _{i}{B}_{{x}_{u}^{\prime} {y}_{u}^{\prime} l,i}={B}_{{x}_{u}^{{\prime} }{y}_{u}^{{\prime} }l}.\end{eqnarray*}$$

Next, in order to ameliorate potential undetected low surface-brightness source flux incursion to measure the variability of B and to average over spectral sensitivity variations in ${D}_{x^{\prime} y^{\prime} l}$, we select 32 boxes of ${(25\mathrm{pix})}^{2}\approx {(5^{\prime\prime} )}^{2}$, in particularly object-poor locations. For each box ${{ \mathcal B }}_{j}=\{{({x}_{u}^{{\prime} },{y}_{u}^{{\prime} })}_{j}\}$, we measure the spatially integrated spectrum:

$$\begin{eqnarray*}{B}_{l}^{(j)} & = & \displaystyle \sum _{({x}_{u}^{{\prime} },{y}_{u}^{{\prime} })\in {{ \mathcal B }}_{j}}{B}_{{x}_{u}^{{\prime} }{y}_{u}^{{\prime} }l}\\ & \mathop{=}\limits^{{\rm{Equation}}(11)} & \displaystyle \sum _{({x}_{u}^{{\prime} },{y}_{u}^{{\prime} })\in {{ \mathcal B }}_{j}}\displaystyle \sum _{i}{F}_{{x}_{u}^{{\prime} }{y}_{u}^{{\prime} }l,i}\,{b}_{i}({x}_{u}^{{\prime} },{y}_{u}^{{\prime} },{\lambda }_{l}).\end{eqnarray*}$$

It turns out that, apart from flux normalization, all ${B}_{l}^{(j)}$ are nearly the same, with differences at the percent level (see Figure 10). This allows us to conclude that

• (i)  
  the spectral part of the co-added flatfield is spatially invariant, as desired and expected after co-adding ∼2500 exposures with different pointing, orientation and wavelength tuning plus averaging over a ∼ 5'' box
• (ii)  
  the spectral shape of the background component is spatially invariant:

  $$\begin{eqnarray}&&{b}_{i}(x^{\prime} ,y^{\prime} ,\lambda )={a}_{i}(x^{\prime} ,y^{\prime} )\cdot b(\lambda )\,,\end{eqnarray} \tag{ 12 }$$

  where a_i is the amplitude of the background.

Proof comes by contradiction: if the spectral shape of the 5'' box-averaged flatfield or background varied with location, ${B}_{l}^{(j)}$ would vary. Notably, we also see no systematic change of the background between the pre- and post-shock regions. Variations of the background flux on smaller (<5'') spatial scales are physically unlikely and in any case must be small as the observed variability can entirely be accounted for by Poisson (measurement) noise (the middle and bottom panels in Figure 10). In the bottom panel of Figure 10, we present differences between the actual background and the modeled background. There is no systematic effect, and the spatial variation is within the measurement uncertainty, which in turn is much smaller than the signal. The variation between background boxes, as measured by the 32 element sample standard deviation normalized by the photon noise, is 0.95 on average across the spectrum (solid purple line), with maximum value of 1.22. The average absolute deviation (dashed purple line) is 0.77 on average, close to the theoretical value of 0.8 expected for the absolute value of a standard normally distributed variable. Therefore, any differences between the actual background ("Data") and model ("Bkg") are fully explained by measurement uncertainties.

We can therefore measure $b(\lambda )$ on the co-added data cube. By combining all 32 background boxes, we additionally minimize noise and residual systematics:

$$\begin{eqnarray*}&&b(\lambda )=\displaystyle \sum _{j}{B}_{l}^{(j)}\equiv \tilde{b}(\lambda ).\end{eqnarray*}$$

The tilde indicates that $\tilde{b}$ is a model of the background spectrum. The normalization of $b(\lambda )$ is absorbed in a_i. We now assume that a_i is spatially invariant, i.e., that the background amplitude changes only between exposures but not across the FOV:

$$\begin{eqnarray*}&&{a}_{i}(x^{\prime} ,y^{\prime} )={a}_{i}.\end{eqnarray*}$$

Therefore, also ${a}_{i}(x,y)={a}_{i}$ is spatially constant, and we have

$$\begin{eqnarray*}&&{D}_{{x}_{u}{y}_{u},i}={B}_{{x}_{u}{y}_{u},i}={F}_{{x}_{u}{y}_{u}}\cdot {a}_{i}\,\tilde{b}({{\rm{\Lambda }}}_{{x}_{u}{y}_{u},i}).\end{eqnarray*}$$

Here, $({x}_{u},{y}_{u})$ are the non-masked pixels of the co-added frame, reprojected onto the individual exposures. We do not know the background amplitudes a_i yet, but can already use the knowledge of $b(\lambda )$ to model F_xy:

$$\begin{eqnarray}&&\left(\displaystyle \sum _{i}\displaystyle \frac{{D}_{{\tilde{x}}_{u}{\tilde{y}}_{u},i}}{\tilde{b}({{\rm{\Lambda }}}_{{\tilde{x}}_{u}{\tilde{y}}_{u},i})}\right)\,\mathop{\longrightarrow }\limits^{\mathrm{fit}+\mathrm{norm}.}\,{\tilde{F}}_{{xy}}.\end{eqnarray} \tag{ 13 }$$

This way, we "divide out" the non-constant background spectrum, which is imprinted on the measurement via ${{\rm{\Lambda }}}_{{xy},i}$. The sum on the left side is only carried out for pixels $({\tilde{x}}_{u},{\tilde{y}}_{u})$ that are not masked in any of the exposures. The result of the sum is a non-normalized flatfield image, which is then modeled by a fourth-order polynomial. The fit eliminates pixel noise and interpolates over masked pixels. It is followed by normalization, such that ${\sum }_{x,y}{\tilde{F}}_{{xy}}={N}_{\mathrm{pix}}$.

Now, we justify that the background amplitude is spatially constant, and even that the flatfield is indeed constant as normally expected. We cannot prove this for individual exposures, as the fluxes are impractically small. However, we can restrict the sum in Equation (13) to different subsets of exposures. This way, we derived ${\tilde{F}}_{{xy}}$ for exposures of only one specific "channel" (tuning) and for subsets of observations at different times. The results are invariant within a few percent: random variations of 3.2% (mean pixel standard deviation) between flatfields of different observing runs, and 1.5% among channel-specific flatfields. Again, the conclusion is by contradiction: given that exposures vary in pointing, field rotation, and tuning, a non-constant flatfield or a spatially variable background amplitude would lead to systematic variations of $\tilde{F}$, but we do not observe such variations.

With ${\tilde{F}}_{{xy}}$ on hand, we measure a_i as the flux scaling required to match the data,

$$\begin{eqnarray*}&&{a}_{i}=\displaystyle \frac{{\displaystyle \sum }_{{x}_{u},{y}_{u}}{D}_{{x}_{u}{y}_{u},i}}{{\displaystyle \sum }_{{x}_{u},{y}_{u}}{\tilde{F}}_{{xy}}\cdot \tilde{b}({{\rm{\Lambda }}}_{{x}_{u}{y}_{u},i})}\,,\end{eqnarray*}$$

and reconstruct the background in each exposure:

$$\begin{eqnarray}&&{\tilde{B}}_{{xy},i}={\tilde{F}}_{{xy}}\cdot {a}_{i}\,\tilde{b}({{\rm{\Lambda }}}_{{xy},i}).\end{eqnarray} \tag{ 14 }$$

We then transform the flatfield and background models ${\tilde{F}}_{{xy}}$ and ${\tilde{B}}_{{xy},i}$ to subcubes in the common astrometric frame and co-add them, in the same way as the data in Equations (6)–(8). As for the "real" flatfield and background, Equation (11) ensures that the resulting co-added cubes correspond to the actual flatfield and background components in the data cube.

We use ${\tilde{B}}_{x^{\prime} y^{\prime} l}$ and ${\tilde{F}}_{x^{\prime} y^{\prime} l}$ to check once more whether our assumption of a uniform background and invariant flatfield are fulfilled: the wavelength-integrated ${(B/F)}_{x^{\prime} y^{\prime} }={\sum }_{l}{B}_{x^{\prime} y^{\prime} l}/{F}_{x^{\prime} y^{\prime} l}$ is indeed flat, and the background residuals ${(D-B)}_{x^{\prime} y^{\prime} }$ are zero apart from Poisson noise and objects (see Figure 11).

Zoom In Zoom Out Reset image size

When gas starts to cool and recombine downstream, the photons produced escape to the pre-shock region and form a PIP (Raymond 1979). Although non-radiative shocks lack recombination zones, PIP can still be created, where the main sources of the photons produced downstream are He i λ584 Å and He II λ304 Å. Ghavamian et al. (2000) reported not only on diffuse Hα, but also [N ii] and [S ii] emission extending over 1' in front of Tycho's NE rim, and suggested that it arises in a PIP. They predicted that the pre-shock gas was heated in the PIP to ∼12,000 K. Subsequently, Lee et al. (2007) measured the PIP spectrum in front of "knot g": narrow Hα with W_NL ≈ 34 km s⁻¹, and [N ii] λ6583 Å with the width of 23 km s⁻¹.

The similarity of the spectra in the pre- and post-shock background boxes (middle panel in Figure 10), where pre-shock boxes partially cover the region of the suspected PIP, unambiguously shows that the PIP emission in our data is negligible. Furthermore, since the signal of a PIP increases toward the shock front, we searched for its signature in nine pre-shock, 13² pixel regions (magenta boxes in Figure 12) that were taken to be closer to the filament than the background boxes (≈7'' away from the filament), but still far enough so that we do not pick up on projected filament emission and emission in the CR precursor. The top-right panel shows the comparison of the putative PIP signal and background model, normalized by measurement noise. The mean normalized PIP level is indicated by the solid -magenta line and is 0.13 on average. Therefore, the possible PIP signal is consistent with zero. We also show the comparison between the mean PIP and background flux, including their 1σ uncertainties (bottom-left panel). Finally, in the bottom-right panel we demonstrate that the PIP signal is negligible in comparison to the filament flux. We chose here the location (bin) with the smallest signal per area; the signal is larger in all other bins, and the putative PIP even smaller in relation to the signal. We therefore conclude that there is no or negligible PIP contribution to our filament flux models.

Zoom In Zoom Out Reset image size

The Voronoi binning (Cappellari & Copin 2003) was performed setting two criteria: the targeted S/N and the minimum bin size. The minimum bin size has to be set in order to account for the seeing of ≃1'', which for the spatial scale of 02 pixel⁻¹ gives the minimal bin size of ≃19 pixels. We set the targeted S/N = 10. The code of Cappellari & Copin (2003) in its standard version tends to create elongated-shaped bins in the direction perpendicular to the shock filament. The HST image of Tycho's SNR seen in Figures 1 and 2 in Lee et al. (2010) shows a complex filamentary structure including several very bright knots as a result of different shock emission projected along the LOS. Small-scale differential gradients are very noticeable in the direction of the shock normal, which is the reason why we used the Weighted Voronoi Tessellation adaptation of Cappellari & Copin (2003; Diehl & Statler 2006) that created a couple of rounder bins in the same direction instead. This way, we created 85 and 15 spatial Voronoi bins in the eastern and northern parts of the rim (Figure 3), respectively. Although we see only 2% level variations between boxes of 25² pixel size, these residual background variations become important relative to the S/N for bin sizes of 400 pixels and larger. This yields the upper size limit of 400 pixels for the bins used in our analysis. Seventy-three bins (out of 85) in the eastern and nine bins (out of 15) in the northern filament fulfill the above criteria and are further considered while the bins with more than 400 pixels are excluded from the analysis.

As already mentioned in Section 3.2, we define the total flux, continuum and line flux fractions, line centroids, and widths as model parameters. In the case of the NLNL and NLNLIL models, we define the NL centroid mean $\langle {\mu }_{\mathrm{NL}}\rangle$ and the separation between the two NLs ${\rm{\Delta }}{\mu }_{\mathrm{NL}}$ as parameters. The IL centroid is introduced with its offset from the NL centroid (mean) ${\rm{\Delta }}{\mu }_{\mathrm{IL}}$.

The parameters that we actually sample from are slightly different from the parameters that we use in Section 3; the difference enables the direct application of the prior PDFs (Dirichlet or Beta distributions) to the parameters. These prior distributions require parameter sets defined in the range (0, 1), except for the naturally based logarithm of the total flux ln $({F}_{\mathrm{tot}})$ for which we use a flat unbound prior. The continuum and line flux fractions are by definition in the (0, 1) range, where we set the continuum flux fraction to be dependent on the flux fractions in the lines ${f}_{{\rm{c}}}=1-\sum {f}_{{\rm{i}}}$. The fact that the flux fractions sum up to 1 and that they are in the range (0, 1) makes the Dirichlet distribution ${\prod }_{{\rm{i}}}{f}_{i}^{{\alpha }_{{\rm{i}}-1}},\,i=[{NL},{NL}1,{NL}2,{IL},c]$ the perfect choice for their prior PDF. Since we do not favor any of the flux components, we use a symmetric Dirichlet distribution with the same index α which we set to $\alpha =1.5$ to disfavor zero fluxes.

The line parameters ${\mu }_{\mathrm{NL}}$, $\langle {\mu }_{\mathrm{NL}}\rangle$, ${\rm{\Delta }}{\mu }_{\mathrm{NL}}$, and ${\rm{\Delta }}{\mu }_{\mathrm{IL}}$ are all defined with the following functional form: ${x}^{{\prime} }=(x-{x}_{\min })/({x}_{\max }-{x}_{\min })$. ${\mu }_{\mathrm{NL}}$ or $\langle {\mu }_{\mathrm{NL}}\rangle$ are defined in the range $[{V}_{\mathrm{cen}}-{V}_{\mathrm{FSR}}/4$, ${V}_{\mathrm{cen}}+{V}_{\mathrm{FSR}}/4$], where V_cen is the center of the free velocity range V_FSR (FSR in velocity units) and ${\rm{\Delta }}{\mu }_{\mathrm{NL}}$ is in the range [0, ${V}_{\mathrm{FSR}}/2]$ so that in the most extreme case the two NL centroids are at the edges of the spectral coverage. ${\rm{\Delta }}{\mu }_{\mathrm{IL}}$ is within $\pm {V}_{\mathrm{FSR}}/4$, having the absolute upper boundary set to the upper (lower) boundary of the NL (IL) width ≈100 km s⁻¹. Instead of line widths W_NL and W_IL defined in the range [15, 100] km s⁻¹ and [100, 350] km s⁻¹, respectively (Morlino et al. 2012, 2013), we used their log-widths (natural logarithm) denoted as ${w}_{\mathrm{NL}(\mathrm{IL})}=\mathrm{ln}{W}_{\mathrm{NL}(\mathrm{IL})}$ in the same functional form as for the centroids and separations. For all of these parameters, we define a symmetric Beta distribution prior ${x}^{\alpha -1}\,{(1-x)}^{\beta -1}$ with the arguments $\alpha =\beta =1.5$—slightly favoring the central values of the defined parameter ranges. The model parameters and their priors are summarized in the Table 1.

In addition to Figure 2, where we plotted posteriors for the favored NLIL model for one of the bins, in Figures 13–15 we show the posteriors for the NL, NLNL, and NLNLIL models, where the latter two figures show the parameters that we actually sample from.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Posterior samples were drawn from an ensemble MCMC sampler (Goodman & Weare 2010), an implementation of which has been popularized as "emcee" (Foreman-Mackey et al. 2013). Among other advantages, this method provides for (near-) optimal tuning at every stage of the sampling, which would otherwise be a substantial challenge and obstacle in the way of efficient sampling considering 82 different data sets, four different models for each of them, up to 10 model parameters, and our intent to test the sampler and results based on hundreds of additional simulated data sets. We draw the initial parameters for 128 parallel chains (walkers) uniformly between the prior boundaries. The unnormalized log-posterior of a model is computed as a sum of unnormalized log-prior and log-likelihood for the proposed walker position. Since our data (fluxes in the spectral bins) result from Poisson processes, the likelihood is the product of each datum's probability under a Poisson distribution with the expectation value equal to the model prediction. After taking the logarithm,

$$\begin{eqnarray}&&\mathrm{ln}L=\displaystyle \sum _{i=1}^{{\rm{N}}=48}{d}_{{\rm{i}}}\mathrm{ln}({m}_{i})-{m}_{i}-\mathrm{ln}({\rm{\Gamma }}({d}_{i}+1))\,,\end{eqnarray} \tag{ 15 }$$

where m_i are the model predictions at each spectral bin i, and d_i are the corresponding data. The last term in Equation (15) (the factorial term) was left out of the posterior sampling (but not the evidence calculation) because it is model independent. We refer to Foreman-Mackey et al. (2013) and Goodman & Weare (2010) for details of the sampling algorithm. Before checking the chains' convergence, we disregard the first 25%, but at least 512, of the samples of each walker ("burn-in"), and further thin the chains until the autocorrelation time of the thinned sample is smaller that 5. In order to achieve low noise and to set a first minimum threshold for the convergence of the chains, we require at least 2¹³ = 8096 total samples (all walkers combined) to be kept after thinning. Once this minimum number of samples is reached, we additionally impose the following convergence (stopping) criterion: we split the sample into subsamples and compute the desired estimators (maximum posterior sample, median, and 95% confidence interval boundaries) for each subsample. The variances of the subsample estimators are then required to be smaller than 5% of the mean parameter value. In order to reduce the probability of coincidentally favorable (small) variances from possible "modes" in the chains, we repeat the process twice, first taking each of the 128 chains as one subsample, and second, taking each of the 128 element walker states as subsamples, and use the arithmetic mean of the resulting eight relative estimator standard deviations toward our 5% criterion.

To ensure that the applied procedure gives the correct results, we performed tests of the posterior sampling routine, using simulated data with known model parameters. We tested models with S/N in the range $[5,50]$ and varying background levels (0%, 50%, and 90% of the total flux). First, we checked if the posterior distribution reproduces the prior in the $S/N=0$ limit, but also if the posterior approached a delta function in the infinity limit, i.e., for very large S/N.

Second, we check if the model parameters of the input model are reproduced statistically. We ran the algorithm for 200 different realizations of the same model, i.e., each time drawing the data from the same model prediction with its specified uncertainty included (each a set of 48 Poisson distributions). The mean of the distribution of the median values is always consistent with the input model parameter and the typical scatter of this mean is 10%–20% for the range of S/N in our real data.

Finally, we vary model parameters by randomly choosing them 200 times from within the prior boundaries. Again, for each of the resulting simulated data, we sample the posterior as described above and evaluate it in the form of the median of the marginalized posteriors. We find that the measured median values scatter symmetrically around the 1:1 relation with the input model parameters, with the scatter being roughly equal to the individual posteriors' standard deviation, as desired. The distribution of the measured values becomes biased as the input parameters approach the parameter range boundaries, as expected for our priors.

We use the CV likelihood, specifically its "LOO" variant, to compute model evidences and to compare models (Bailer-Jones 2012). We prefer it over the standard numeric ("Bayesian") integral, because it draws samples from the posterior instead of the prior. It is hence more efficient, less dependent on the choice of prior, and in some cases numerically more stable. The idea of the CV likelihood is to evaluate the likelihood of part of the data, given the model and the rest of the data. In our case, we have 48 data points and measure how well any 47 data points (the complement, ${D}_{-k}$) under the model M predict the 48th data point (the partition, D_k), as quantified by the complement's posterior $P(\theta | {D}_{-k},M)$ and its prediction for D_k. The process is repeated for all possible (48) partitions. Each time we leave out one datum (D_k), its partition likelihood L_k is given by

$$\begin{eqnarray}{L}_{k} & = & P({D}_{k}| {D}_{-k},M)\\ & = & {\displaystyle \int }_{\theta }P({D}_{k}| \theta ,M)P(\theta | {D}_{-k},M)d\theta .\end{eqnarray} \tag{ 16 }$$

The first term in the integrand is the likelihood of D_k, while the second term is the posterior PDF after considering the information contained in ${D}_{-k}$. We can numerically (Monte Carlo) integrate by drawing a number of samples N from $P(\theta | {D}_{-k},M)$:

$$\begin{eqnarray}&&{L}_{k}\approx \displaystyle \frac{1}{N}\displaystyle \sum _{n=1}^{N}P({D}_{k}| {\theta }_{n},M).\end{eqnarray} \tag{ 17 }$$

Assuming that the data points are independent, the LOO-CV likelihood is the product of all partition likelihoods, ${L}_{\mathrm{LOO} \mbox{-} \mathrm{CV}}={\prod }_{k=1}^{48}{L}_{k}$, or

$$\begin{eqnarray}&&\mathrm{ln}{L}_{\mathrm{LOO} \mbox{-} \mathrm{CV}}=\displaystyle \sum _{k=1}^{48}\mathrm{ln}{L}_{k}.\end{eqnarray} \tag{ 18 }$$

It can be shown that ${L}_{\mathrm{LOO} \mbox{-} \mathrm{CV}}$ for model M is equal to the Bayesian evidence (E(M)), and we henceforth use it to compute the Bayes factors, $E({M}_{1})/E({M}_{2})$, to compare models.

We tested the ability of the Bayes factors to indicate the correct model by employing the simulated data. We find that the typical numerical precision of $\mathrm{ln}{L}_{\mathrm{LOO} \mbox{-} \mathrm{CV}}$ is better than 0.05 dex. At the same time, for data generated from the parameters and with S/Ns that are typical for the actual data, the "right" model's $\mathrm{ln}{L}_{\mathrm{LOO} \mbox{-} \mathrm{CV}}$ is 0.2 dex (≳50%) better than any of the alternative models. Nevertheless, on heuristic grounds (what probability is considered statistically "significant enough"?) and in order to bracket the practically limited scope of such tests, we adopt a more conservative $+0.5\,\mathrm{dex}$ threshold ($3\,:1$ probability) before we consider a model to be preferable over another.

We have also tested the dependence of the evidence ratios on the choice of prior. We used the same functional form for the priors (Table 1), but different α, β distribution parameters. Specifically, we tested ${\alpha }_{{\rm{D}}}$ values for the Dirichlet distribution, and ${\alpha }_{{\rm{B}}}={\beta }_{{\rm{B}}}$ for the Beta distribution with values of $({\alpha }_{{\rm{D}}}$, ${\alpha }_{{\rm{B}}})=\{(1.5,1.5),(1,1.5),(1.5,1),(1,1)\}$, where $(1.5,1.5)$ was used for the results presented in the main part of the paper. We found that for all bins, the mean standard deviation of the appropriate evidence ratios is ≲0.07 dex.

In addition, in Figures 16 and 17 we show the posteriors, medians, and 95% confidence intervals when flat priors, i.e., $({\alpha }_{{\rm{D}}}$, ${\alpha }_{{\rm{B}}})=(1,1)$ instead of the fiducial $(1.5,1.5)$, are applied. We present results for the same two Voronoi bins shown in Figures 2 and 5. The shape of the posteriors and their medians are very similar to those obtained with the adopted priors, demonstrating that our results do not depend strongly on the choice of prior. W_IL posteriors for the adopted and flat prior are different, but similar to the degree that the median and confidence interval boundaries change only by 30 km s⁻¹ ($\approx 10 \%$). Even with the flat prior, 210 km s⁻¹ is clearly more probable than the other W_IL parameter values; in particular, it is preferred over values close to the W_IL limits. This shows that the preference for the central W_IL values is not just borne out of the prior shape or range, but genuinely reflects constraints provided by the data, even if they are not as strong as for other parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size