© ESO 2020

1. Introduction

Supernova remnants (SNRs) are among the most exciting astrophysical objects because they may shed light on a major question of physics: the origin of cosmic rays. An interesting object in this respect is the remnant of the historical Type Ia event SN 1572, first described by Tycho Brahe. Thanks to his observation, the age of Tycho’s SNR (SNR G120.1+1.4, hereafter Tycho) is accurately determined to be ∼440 years. Tycho originated from a Type Ia supernova (SN) (Krause et al. 2008), and is assumed to have a canonical explosion energy of ∼10⁵¹ erg. These details, together with the available broadband spectrum, make Tycho one of the best astrophysical laboratories to study particle acceleration. Nevertheless, several questions still remain unanswered.

Indeed, despite the success of diffusive shock acceleration (DSA) theory, it fails to explain the observed radio spectrum S_ν. The measured radio spectral index σ ≈ 0.65, with S_ν ∝ ν^−σ (Kothes et al. 2006), significantly deviates from the standard DSA prediction σ ≈ 0.5. This discrepancy is generally accounted for using the concept of Alfvénic drift (Bell 1978). Various authors postulate Alfvénic drift only in the upstream region (Völk et al. 2008; Morlino & Caprioli 2012), or in the upstream and downstream (Slane et al. 2014) regions of Tycho’s forward shock. The proper motion of cosmic-ray scattering centers that proceed with Alfvén speed is assumed to decrease the compression ratio felt by the particles, thereby causing a softening of their spectra. To be more precise, the effective compression ratio seen by particles is

$$\begin{array}{c}\hfill r_{\mathrm{eff}}=\frac{u_1+H_{\mathrm{c}1}{v}_{\mathrm{A}1}}{u_2+H_{\mathrm{c}2}{v}_{\mathrm{A}2}},\end{array}$$(1)

where u₁, v_A1 are the plasma velocity and the Alfvén velocity in the upstream and u₂, v_A2 in the downstream regions, respectively. The relative direction of the propagation of the Alfvén waves in the upstream and downstream is reflected by the cross helicities, H_c1, and H_c2, respectively. In the above global models its value is chosen rather freely¹, in order to reduce the effective compression ratio and thereby account for the soft particle spectrum. However, on closer inspection, the concept of Alfvénic drift as an explanation for the spectral softening appears misleading.

First of all, Vainio & Schlickeiser (1999) performed a detailed calculation on self-consistent transmission of the Alfvén waves through the shock and found that the presence of waves results in harder particle spectra than predicted by the standard theory (meaning that Alfvén waves infer exactly the opposite effect to that claimed by Völk et al. 2008; Morlino & Caprioli 2012, and Slane et al. 2014). This is because Alfvén waves that move in the upstream region in the opposite direction to the background plasma (H_c1 = −1) propagate also predominately in the opposite direction in the downstream (H_c2 ≈ −1) (Vainio & Schlickeiser 1999). Therefore, despite the modifications induced by the strong magnetic-field pressure, the effective compression ratio seen by particles (Eq. (1)) exceeds the standard strong-shock value (r_eff >  4). Investigating the impact of the Alfvénic drift within the framework of Vainio & Schlickeiser (1999) we find that for the Alfvénic Mach numbers, M_A ≡ u₁/v_1, A, in the range 10–13 (as presented in Völk et al. 2008 and Morlino & Caprioli 2012) the particle spectral index for a strong shock results in s ≈ 1.9 instead of the s ≈ 2.3 required by radio observations of Tycho. Furthermore, the negative downstream helicity, H_c2 ≈ −1, predicted by Vainio & Schlickeiser (1999) is exactly the opposite to what was assumed by Slane et al. (2014), who used H_c2 = 1.

Secondly, the Alfvénic-drift phenomenon in the global models of Tycho is often referred to as resulting from the non-resonant streaming instability of cosmic rays (Bell 2004). According to Bell (2004), the phase speed of the nonresonant modes is negligible compared to the shock velocity. But in case of Alfvénic drift, the Alfvén velocity required to account for Tycho’s radio spectra has to be enormous.

The following estimation should demonstrate the corresponding discrepancy. Let us assume that Alfvénic drift occurs only in the upstream region (H_c2 = 0 and H_c1 = −1), as in Völk et al. (2008) and Morlino & Caprioli (2012), even though this contradicts the findings of Vainio & Schlickeiser (1999). In this case, we can rearrange Eq. (1) to

$$\begin{array}{c}\hfill M_{\mathrm{A}}={(1-r_{\mathrm{eff}}/r_{\mathrm{sh}})}^{-1},\end{array}$$(2)

where M_A ≡ u₁/v_1, A is the Alfvénic Mach number and r_sh ≡ u₁/u₂ is the gas compression ratio of the shock. The sub-shock compression ratio in Völk et al. (2008) and Morlino & Caprioli (2012) is in the range r_sh = 3.7−3.9, even if the nonlinear effects of the DSA are included. The effective compression ratio seen by particles and required by the radio observations is r_eff ≈ 3.3. Inserting these values into Eq. (2) provides a relatively low Alfvénic Mach number, M_A = 6−9, making it clear that the Alfvén speed exhibits a significant fraction of the shock velocity. Therefore, to explain Tycho’s radio data with Alfvénic drift, the Alfvén phase speed has to attain 11%−16% of the shock velocity. Obviously this value is in conflict with the phase speed of the nonresonant mode, v_ϕ ≈ 0, as described by Bell (2004). Therefore, Alfvénic drift in the models for Tycho cannot be associated with the nonresonant streaming instability.

It is important to note here that Völk et al. (2008) and Morlino & Caprioli (2012) applied in their modeling somewhat smaller magnetic-field values than required to best match the radio data with Alfvénic drift. The post-shock magnetic fields of ∼300 μG (Morlino & Caprioli 2012) and ∼400 μG (Völk et al. 2008) provide sufficient flux in the radio range, but fit the spectral shape of the observed data only moderately well. The Alfvén velocity in these models (with corresponding Alfvénic Mach numbers of M_A ≈ 13 and M_A ≈ 10) is still ∼10% of the shock velocity. The phase speed of the nonresonant mode is much less than that and cannot support Alfvénic drift.

A post-shock magnetic field above ∼300 μG, when combined with the relatively low ambient density of 0.3−0.4 cm⁻³ (Völk et al. 2008; Morlino & Caprioli 2012), becomes dynamically important. The corresponding magnetic-field pressure will affect the shock compression ratio, which results in r_sh <  3.9. Aside from the work of Morlino & Caprioli (2012), this effect has been neglected in global models for Tycho.

The above arguments illustrate that the Alfvénic-drift concept is unable to explain Tycho’s soft radio spectrum consistently. An alternative way to explain the softening of the particle spectra in collisionless shocks is to account for neutral hydrogen in the surrounding medium, first proposed by Blasi et al. (2012). Analytic calculations from Ohira (2012) plus later simulations (Ohira 2016) show that neutrals can leak from the downstream into the upstream region and modify the shock structure. This results in a steeper particle spectrum as produced by standard DSA. Morlino & Blasi (2016) build on that idea to model the rather soft γ-ray spectrum of Tycho. However, the leakage of neutral particles is significant for shock velocities V_sh <  3000 km s⁻¹ (Ohira 2012), which is considerably below the value ascertained for Tycho. Morlino & Blasi (2016) argued that certain regions of Tycho can feature slower shocks that propagate into dense, partially neutral material. The regions with slower shock velocities would have to provide nearly all of the observed overall emission, while contributions from the regions with the fast shock velocities would have to be weak, otherwise the integrated emission would reflect the hard spectra expected for fast shocks in an ionized medium.

In this work we suggest a new approach: besides standard DSA, we consider an additional acceleration process, namely the stochastic re-acceleration of particles in the immediate post-shock region of the SNR. It has been shown that fast-mode waves that survive the transit-time damping by the background plasma are efficient in accelerating charged particles via cyclotron resonance (Yan & Lazarian 2002; Liu et al. 2008). Pohl et al. (2015) demonstrated that particles may be stochastically re-accelerated by fast-mode turbulence, which occurs in the downstream region, after escaping from the forward shock. In this work we build on that idea and incorporate it into a detailed modeling of Tycho.

Fast-mode turbulence may arise behind the shock from velocity fluctuations of the plasma flow, for example via shock rippling (Giacalone & Jokipii 2007), and build a thin turbulent region behind the blast wave. The width of this zone is regulated by the energy transfer from the background plasma into turbulence as well as damping induced by the re-acceleration of cosmic rays. We show that fast-mode turbulence that carries a few percent of the energy density of the background plasma in the downstream region is strong enough to modify the spectrum of particles that have already been accelerated by the shock. Thus, in our treatment, stochastic acceleration and DSA operate together and produce a particle spectrum consistent with the observed radio spectral index. An additional advantage of our approach is that it is fully time dependent². We solve the time-dependent transport equation for cosmic rays that contains a DSA term and diffusion in momentum space, and is coupled to hydrodynamic simulations.

Another interesting question regarding Tycho pertains to the extent to which the magnetic field is amplified inside the remnant. A relatively high post-shock magnetic field, 300 −400 μG, is postulated by several global models (Völk et al. 2008; Morlino & Caprioli 2012). One of the reasons is the aforementioned Alfvénic drift, which demands relatively large magnetic-field values to account for the radio spectrum. Since we do not postulate any Alfvénic drift in our model, our approach of inferring the magnetic-field strength is an alternative to that of previous works on Tycho.

A major argument for a high magnetic field downstream of Tycho are the observed narrow nonthermal X-ray filaments (Hwang et al. 2002; Parizot et al. 2006; Helder et al. 2012). Since electrons can only propagate for a finite distance before they lose their energy due to synchrotron radiation, the rim widths may reflect the magnetic-field strength immediately downstream of the SNR. An alternative scenario is provided by damping of the turbulent magnetic field in the interior of the remnant (Pohl et al. 2005; Ressler et al. 2014; Tran et al. 2015), in which the narrowness of the nonthermal rims can be explained by the damping of the turbulently amplified magnetic field. For Tycho, a distinction between these scenarios by means of the energy dependence of X-ray filaments is difficult (Rettig & Pohl 2012; Tran et al. 2015). Magnetic-field damping is widely considered as the scenario that allows for a weak magnetic-field strength inside SNRs. Nevertheless, Tran et al. (2015) find that in either case the minimum downstream magnetic-field value inferred from Tycho’s nonthermal filaments is at least ∼20 μG. Assuming that the electron acceleration is limited by the age of the remnant, the work from NuSTAR collaboration (Lopez et al. 2015) estimates ∼30 μG for the downstream magnetic field. However, the majority of studies cited above favor the loss-limited interpretation for particle acceleration in Tycho. Furthermore, the most realistic limit is obtained from an analysis of the entire spectral energy distribution (SED), which provides the minimal post-shock magnetic field value of ∼80 μG (Acciari et al. 2011). Attempts to simultaneously fit the radio and the γ-ray data infer that any weaker magnetic field would cause an overproduction of γ-ray photons generated via inverse-Compton scattering. Therefore, the question about the magnetic field value is automatically tied to the question of whether Tycho’s γ-ray emission has a predominately leptonic or hadronic origin. The hadronic scenario has been strongly favored in the literature (Morlino & Caprioli 2012; Zhang et al. 2013; Berezhko et al. 2013; Caragiulo & Di Venere 2014; Slane et al. 2014), as opposed to a leptonic model (Atoyan & Dermer 2012).

For our modeling, we start from the minimal magnetic field compatible with the entire SED. The evolution of the SNR, which is computed using hydrodynamical simulations, occurs in a medium with a constant density. We explicitly model the acceleration of each particle species in the test-particle limit, taking shock acceleration as well as stochastic acceleration in the downstream region into account, with both acceleration processes being time dependent. We explicitly model advection and diffusion of cosmic rays and take synchrotron losses for electrons into account. Furthermore, we consider the nonthermal radio and X-ray filaments and investigate whether they arise from extensive synchrotron losses or magnetic-field damping. For the study of the X-ray filaments, it is especially important to include diffusion of the particles, otherwise the distance that they propagate and thus the rim width would be underestimated.

This paper is organized as follows: Sect. 2 describes our method: hydrodynamical picture, magnetic profiles, and particle acceleration via DSA and Fermi II processes. In Sect. 3.1 we examine the model with the minimal amplified magnetic field compatible with the γ-ray observations. We justify the necessity for magnetic-field damping in Sect. 3.2. Therein, we discuss our favorite model for Tycho, and deduce a new theoretical minimum for the magnetic-field strength based on investigation of the multi-frequency spectrum along with nonthermal filaments. The γ-ray spectrum of the resulting model comprises both hadronic and leptonic components. In Sect. 3.4 we discuss the potential for a purely hadronic scenario with a strongly amplified magnetic field.

2. Modeling

2.1. Hydrodynamics

To accurately calculate the acceleration of particles to relativistic energies and the high-energy emission from the remnant, we start with an evolutionary model of the remnant. This model provides the dynamical and kinematic properties of the remnant as a function of time. We model the hydrodynamic evolution of the SN shock wave in Tycho using the VH-1 code, a 1, 2, and 3 dimensional finite-difference code that solves the hydrodynamic equations using the Piecewise Parabolic Method of Colella & Woodward (1984). Dwarkadas & Chevalier (1998) showed that the ejecta structure of Type-Ia SNRs such as Tycho could be best approximated by an exponential density profile. Using this profile, the parameters needed to model the SNR are the explosion energy of the SN, the ejected mass, and the density of the medium into which the SN is expanding. We use the canonical value of 10⁵¹ erg for the energy of the explosion. Type Ia SNs are presumed to arise from the thermonuclear deflagration and detonation of a white dwarf, and so we assume an ejecta mass of 1.4 M_⊙, appropriate for a Chandrashekhar-mass white dwarf. The remaining parameter necessary to model the evolution is the density of the ambient medium.

Although the density around Tycho varies (Williams et al. 2013), as expected for such an extended structure, the remnant appears to expand in a clean environment without any large inhomogeneities, such as molecular clouds (Tian & Leahy 2011). Therefore, an explosion in a medium with a constant density provides a reasonable description for the dynamics of the remnant. The average ambient density varies in the literature according to the method with which it is measured. X-ray measurements of the expansion rate suggest an ambient density of 0.2–0.6 cm⁻³ (Hughes 2000). Later observations infer an upper limit of 0.2 cm⁻³ (Katsuda et al. 2010). X-ray observations of Cassam-Chenaï et al. (2007) reveal a lack of thermal emission in the post-shock region of Tycho, inferring an ambient density below 0.6 cm⁻³. A low density of 0.2 cm⁻³ around Tycho is obtained by Williams et al. (2013), who determined the post-shock temperature from mid-infrared emission of the remnant. On the other hand, efficient particle acceleration in SNRs can reduce the downstream temperature of the plasma (O’C. Drury et al. 2009), leading to a suppression of thermal emission and hence an underestimation of the ambient density. Higher values for the ambient density are additionally supported by Dwarkadas & Chevalier (1998), who found that densities in the range of 0.6–1.1 cm⁻³ better match the observations of Tycho. Furthermore, Kozlova & Blinnikov (2018) and Badenes et al. (2006) favor a delayed detonation model with ∼1.0 cm⁻³ for the X-ray morphology of Tycho. The density uncertainty of Tycho (0.2–1.0 cm⁻³) is likewise reflected in the uncertainty in the distance to the remnant (see Hayato et al. 2010 for a review and Tian & Leahy 2011), since both quantities are interdependent. For our modeling we choose a value for the hydrogen number density of n_H = 0.6 cm⁻³, which gives a good fit to the observed shock radii and velocities.

Explosion energy, ejecta mass, and ambient density, together with the exponential density profile, form the suite of parameters necessary to model the complete hydrodynamical evolution of Tycho. Our simulations are spherically symmetric, and are similar to those described in Dwarkadas & Chevalier (1998) and Telezhinsky et al. (2012). The evolution of the remnant into the surrounding medium gives rise to a double-shocked structure consisting of a forward shock expanding into the surrounding medium and a reverse shock that propagates back into the ejecta. The two are separated by a contact discontinuity that divides the shocked ejecta from the shocked surrounding medium. The time-evolution of the shock speed, radius, and temperature derived from our simulations is depicted in Fig. 1. At 440 years, our simulations provide a forward shock radius R_sh ≈ 3.5 pc, which implies a distance to the remnant, d ≈ 2.9 kpc. The velocity of the forward shock yields V_sh ≈ 4100 km s⁻¹. The position of the reverse shock, ∼0.69 R_sh, is in good agreement with measurements of Warren et al. (2005). By contrast, the position of the contact discontinuity (CD), R_CD ≈ 0.78 R_sh, falls below the value identified by Warren et al. (2005), who interpreted the closeness of the CD to the blast wave as evidence for the efficient back-reaction of cosmic rays. However, global models that incorporate nonlinear DSA (NLDSA) effects (Morlino & Caprioli 2012; Slane et al. 2014) fail to reproduce the CD position for Tycho. The discrepancy for the CD position can be attributed to the decelerating CD being unstable to the Rayleigh-Taylor instability. Two-dimensional hydrodynamical simulations with the exponential profile show that Rayleigh-Taylor structures can extend almost halfway from the CD to the outer shock (Dwarkadas 2000; Wang & Chevalier 2001). Furthermore, Orlando et al. (2012) showed that Rayleigh-Taylor instabilities and ejecta fingers that extend far beyond the CD can misleadingly suggest that the CD is further out than its actual location. Therefore, we conclude that the position of the CD obtained in our model is quite reasonable.

Fig. 1. Shock speed (blue solid line), radius (green dashed line), and temperature (red dotted line) as functions of time.

According to our simulations, at the age of 440 years the remnant accumulated ∼3.8 M_⊙ of ambient gas, indicating that in our modeling Tycho is in the transition between ejecta-dominated and Sedov-Taylor stages. The total thermal energy in the remnant at 440 years is E_th ≈ 5.6 × 10⁵⁰ erg. The shock profiles and velocity distribution from the simulation are used in the calculation for the particle acceleration described in Sect. 2.3.

2.2. Magnetic field

In this work, the canonical value for the magnetic-field strength of the interstellar medium (ISM) in the far upstream is used: B_ISM = 5 μG. However, its exact value is insignificant for our modeling, because the magnetic field is amplified in the vicinity of the shock, for example due to streaming instabilities (Bell 2004; Lucek & Bell 2000) or turbulent dynamos (Giacalone & Jokipii 2007). Thus, the amplified magnetic field represents a physically relevant parameter. The exact treatment of the magnetic-field amplification is beyond the scope of this work. We adopt instead for the magnitude of the magnetic-field strength in the immediate upstream region B₁ = αB_ISM, where α, a generic amplification factor, is a free parameter in our model. Following the ansatz used by Zirakashvili & Ptuskin (2008) and Brose et al. (2016) the magnetic-field profile in the precursor is assumed to drop exponentially to the interstellar field, B_ISM, at a distance of 5% of the shock:

$$\begin{array}{c}\hfill B(r)=\{\begin{array}{cc}\alpha B_{\mathrm{ISM}}·exp(\frac{-(r/R_{\mathrm{sh}}-1)}{0.05}·ln\left(\alpha \right))\hfill & \text{for}{R}_{\mathrm{sh}}\le r\le 1.05·R_{\mathrm{sh}}\hfill \\ B_{\mathrm{ISM}}\hfill & \text{for}r\ge 1.05·R_{\mathrm{sh}}.\hfill \end{array}\end{array}$$(3)

We stress that the magnetic-precursor length scale is not related to a free-escape boundary, which is usually introduced in the global modeling of SNRs (Völk et al. 2008; Morlino & Caprioli 2012; Slane et al. 2014). The precursor merely reflects typical characteristics of the spatial profile of the amplified magnetic field in the upstream region. In contrast to previous models, our approach does not need an escape boundary, because our system is large enough (ca. 65 shock radii) to retain all injected particles in the simulation. The pre-shock magnetic field is assumed to be isotropic, that its individual components are equal in their magnitudes (for more details see Telezhinsky et al. 2013). Shock compression for the strong unmodified shock yields an immediate downstream value of $B_2=\sqrt{11}{B}_1$.

We consider two scenarios for the magnetic-field distribution inside the SNR. In the first case, the immediate downstream value of the magnetic field is transported inside the SNR with the plasma flow and evolves according to the induction equation for ideal magnetohydrodynamics (MHD; Telezhinsky et al. 2013):

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}=\mathbf{\nabla }\times (\mathit{u}\times \mathit{B}),\end{array}$$(4)

where B is the magnetic-field vector and u the plasma velocity, which are obtained from hydrodynamic simulations. Equation (4) accounts for advection, stretching, and compression of the field, implying that magnetic flux is conserved during the entire evolution of the SNR.

For the second case, we assume that after being amplified in the upstream region and compressed at the shock, the magnetic field decays due to the damping of magnetic turbulence in the downstream region. Magnetic-field damping is one of the key processes in astrophysical plasmas (Kulsrud & Cesarsky 1971; Threlfall et al. 2011), and is expected to operate in SNRs, where it dissipates turbulently amplified magnetic field (Pohl et al. 2005). We do not know the turbulence mode that is most relevant for magnetic-field amplification in SNRs, and hence the particulars of damping are beyond the scope of this paper. Instead, we adopt a simple parametrization for the magnetic profile in the downstream region:

$$\begin{array}{c}\hfill B(r)=B_0+(B_2-B_0)·exp\left(\frac{-(R_{\mathrm{sh}}-r)}{l_{\mathrm{d}}}\right),\end{array}$$(5)

where l_d is the characteristic damping length of magnetic turbulence. The residual field level, B₀, is the limit value to which the magnetic-field tends in the far downstream, behind the CD and the reverse shock of the SNR. In contrast to the previous case, the magnetic flux is not conserved in the scenario described by Eq. (5). Far inside the remnant, the advection of magnetic field has to become important as the magnetic field approaches the constant value B₀. However, this inner field is physically nonrelevant. First of all, we do not include the acceleration of particles at the reverse shock, since we find the radiative contribution of those particles to be negligible. Furthermore, the bulk of the accelerated particles reside in the vicinity of the forward shock and therefore do not experience the constant magnetic field, B₀. Most of the particles cannot cross the CD, which separates the swept-up material from the ejecta because their dynamics are governed by advection. Only a marginal fraction of the very highly energetic particles overcome the CD via diffusive transport. Their radiative contribution to the synchrotron spectrum does not exceed 1% of the overall emission, and is thus negligible for the global SED modeling. Therefore, we can omit the advection of magnetic field for the second scenario that includes magnetic-field damping.

The separate treatment of hydrodynamics and magnetic field does not allow us to include dynamic feedback from magnetic pressure. Therefore, we investigate the range where the impact of magnetic field on the shock structure remains negligible. We analytically solve the classical MHD equations for a plane-perpendicular shock in steady state, assume frozen-in plasma, and derive the dependency of the shock compression ratio on the magnetic-field strength. Apart from the magnetic field, the resulting gas compression ratio depends mainly on shock speed and upstream density, which in our model are V_sh≈ 4100 km s⁻¹ (at 440 years) and n_H = 0.6 cm⁻¹ (see Sect. 2.1). As the extensive magnetic-field pressure reduces the compression ratio of the shock, it consequentially softens the particle spectra. Hence, the corresponding particle spectral index deviates from the canonical DSA solution (s = 2.0) by Δs, resulting in s = 2 + Δs. The resulting radio spectral index, σ = 0.5 + Δσ, is accordingly modified by the value Δσ = Δs/2. We choose Δσ = 0.01 (and thus Δs = 0.02) as the limit within which the magnetic-field pressure can be neglected because the corresponding changes are hardly visible in the radio spectrum when compared to the observed data. We find for the sonic upstream Mach number of M ≈ 3500, which corresponds to the hydrodynamic quantities we use in this work for Tycho, that the magnetic field is dynamically unimportant for Alfvénic Mach numbers above M_A, c ≈ 17. For lower Alfvénic Mach numbers, the deviation from the classical DSA solution exceeds the values of s = 2.02 and σ = 0.51 and hence cannot be neglected. Taking into account that the magnetic-field pressure is exerted only by the tangential components, we find the magnetic-field limits B_1, max ≈ 120 μG for the upstream and B_2, max ≈ 400 μG for the downstream regions, respectively. Any higher magnetic field would efficiently lower the shock compression ratio and produce particle spectra with s >  2.02 and radio spectra with σ >  0.51, and thus demand an MHD treatment.

2.3. Particle acceleration

To determine the evolution of the particle number density, N(r, p, t), we numerically solve the time-dependent transport equation for cosmic rays:

$$\begin{array}{cc}& \mathrm{\nabla }(D_r\mathrm{\nabla }N-\mathit{u}N)-\frac{\partial }{\partial p}(N\dot{p}-\frac{\mathrm{\nabla }\mathit{u}}{3}Np)+\frac{\partial }{\partial p}\left(p^2{D}_p\frac{\partial }{\partial p}\frac{N}{p^2}\right)\hfill \\ \hfill & +Q=\frac{\partial N}{\partial t}.\hfill \end{array}$$(6)

Here, D_r = ξ c r_g/3 is the spatial diffusion coefficient, with the speed of light c and the Larmor radius r_g. The parameter ξ is the ratio between the spatial diffusion coefficient and that for Bohm diffusion, sometimes referred as the gyrofactor. The Larmor radius and hence the spatial diffusion coefficient are calculated using the local magnetic-field strength. Because in our approach the magnetic field is spatially variable (see Sect. 2.2), the spatial diffusion coefficient also varies with position. The synchrotron losses for electrons are included in Eq. (6) via $\dot{p}=(4e^4{B}^2)/(9c^6{m}_e^4)p^2$, with the elementary charge e and the electron rest mass m_e. The plasma velocity, u, is obtained from the hydrodynamic simulations.

The injection of particles is a complex issue that is not fully understood. For simplicity, we use the thermal leakage injection model as presented in Blasi et al. (2005). The source term for electrons and protons is given by

$$\begin{array}{c}\hfill Q_i={\eta }_i{n}_{1,i}{V}_{\mathrm{sh}}\delta (r-R_{\mathrm{sh}})\delta (p-p_{\mathrm{inj},i}),\end{array}$$(7)

where index i denotes the particle species (electrons and protons), η_i the corresponding injection efficiency, and n_1, i the particle number density in the upstream region. The particles in our approach are mono-energetically injected. The associated injection momentum is multiple of their thermal momentum:

$$\begin{array}{c}\hfill p_{\mathrm{inj},i}=4.45p_{\mathrm{th},i}\equiv 4.45\sqrt{2m_i{k}_{\mathrm{B}}T},\end{array}$$(8)

where k_B is the Boltzmann constant and T is the temperature of the plasma.

Electrons and ions are not in equilibrium at collisionless shocks. In fact, the current Rankine-Hugoniot ion temperature is about 30 keV, whereas the observed electron temperature in the post-shock medium is on the order of 1 keV after years of residence in the downstream region. The thermal leakage model is based on the requirement that particles see the shock as a sharp discontinuity and hence are injected only if their Larmor radius exceeds the width of the shock wave. This condition requires that electrons be pre-accelerated to around 100 MeV to participate in DSA. Particle-in-cell simulations provide evidence that thermal electrons can indeed be pre-accelerated by shock-surface and shock-drift acceleration (Matsumoto et al. 2017; Bohdan et al. 2017). Tests demonstrate (Katou & Amano 2019) that this process yields an electron spectra shaped as a power law with a spectral index in the range ∼(1.5−5.5). The exact value of the spectral index depends on the internal structure of the shock transition region. We approximate the two-step energization of electrons by DSA at all momenta and electron injection at very low momentum, which intends to replace an advanced, computationally far more expensive treatment of the pre-acceleration. Thus, the same particle temperature set for Eq. (8) is supposed to mimic the pre-acceleration process of electrons. Although, we do not know the exact value of the spectral index provided by the shock-drift acceleration, it only marginally affects our final spectra, as shown in Pohl et al. (2015, Fig. 4).

Since the main goal of this work is to explain the broad SED of Tycho, we can neglect the thermal and supra-thermal protons. Indeed, as the hadronic γ-rays production starts from the threshold energy ∼0.3 GeV, the protons with momenta below p_inj, p are insignificant for the final photon spectrum. In contrast, the electrons at lower energies play an important role, as the corresponding synchrotron emission is visible in the radio range.

Besides the full time dependence of the particle acceleration and transport, the novel aspect of our work is the stochastic acceleration of particles presented in Eq. (6) by the momentum-diffusion coefficient D_p. Pohl et al. (2015) explicitly derived the momentum-diffusion coefficient for the fast-mode waves and demonstrated that for low-energy particles the acceleration time is independent of momentum and can be of the order of a few years. At higher energy, the process is expected to become less efficient. Stochastic acceleration is an important damping process of fast-mode waves (cf. Thornbury & Drury 2014), and so it is efficient only in a thin region behind the forward shock with the thickness, L_fm. A useful parametrisation of the momentum diffusion coefficient, D_p, is then (Pohl et al. 2015)

$$\begin{array}{c}\hfill D_p(r,p)=\{\begin{array}{cc}0\hfill & \text{for}r<(R_{\mathrm{sh}}-L_{\mathrm{fm}})\hfill \\ \frac{p^2}{{\tau }_{\mathrm{acc}}}f(p)\hfill & \text{for}(R_{\mathrm{sh}}-L_{\mathrm{fm}})\le r\le R_{\mathrm{sh}}.\hfill \end{array}\end{array}$$(9)

Here, τ_acc denotes the acceleration time at small momenta, and f(p) approximates the loss in acceleration efficiency at higher energies as a power law,

$$\begin{array}{c}\hfill f(p)=\{\begin{array}{cc}1\hfill & \text{for}p\le p_0\hfill \\ {\left(\frac{p}{p_0}\right)}^{-m}\hfill & \text{for}p\ge p_0.\hfill \end{array}\end{array}$$(10)

The critical momentum p₀ and the power-low index m are free parameters of our model for now. It is important to note here that the energy in fast-mode waves has mostly a kinetic character. As the magnetic component of the fast-mode turbulence in the typical post-shock plasma is weak, it cannot amplify the magnetic field sufficiently. Therefore, in our approach, other type of waves, such as streaming instabilities (Bell 2004; Lucek & Bell 2000) or turbulent dynamos (Giacalone & Jokipii 2007), are required to provide the magnetic-field amplification at the shock, which is an important scale factor for scattering (see Sect. 2.2). There is not necessarily a simple scaling between momentum and spatial diffusion (e.g., Yan & Lazarian 2004, Eq. (14)). In addition, different types of turbulence may be responsible for the diffusive transport and stochastic re-acceleration (e.g., Shalchi 2009, p. 23). Therefore, the spatial and momentum diffusion coefficients are independent in our model.

The energy in the fast-mode turbulence that occurs in the post-shock region can be primarily provided by the background plasma. To simplify matters, in our approach the energy density in the fast-mode waves scales with the thermal energy density of the post-shock background plasma

$$\begin{array}{c}\hfill U_{\mathrm{fm}}=ϵU_{\mathrm{th}}.\end{array}$$(11)

Here, ϵ is the energy-conversion factor, which is assumed to be of the order of a few percent. The minor value of ϵ provides that the energy transfer from the main plasma flow can be neglected in our hydrodynamic simulations. For a strong shock expanding in a cold plasma, the downstream thermal energy is approximately

$$\begin{array}{c}\hfill U_{\mathrm{th}}\approx \frac{9}{8}{\rho }_{\mathrm{ISM}}{V}_{\mathrm{sh}}^2,\end{array}$$(12)

where ρ_ISM denotes the ambient gas density.

The acceleration time for the fast-mode waves, derived in Pohl et al. (2015), is given by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{acc}}\approx (0.63\mathrm{yr})\left(\frac{U_{\text{th}}}{10U_{\text{fm}}}\right).\end{array}$$(13)

If energy density of the fast-mode waves scales with the thermal post-shock energy density of the plasma, as reflected by Eq. (11), Eq. (13) simplifies to

$$\begin{array}{c}\hfill {\tau }_{\mathrm{acc}}\approx (6.3\mathrm{yr}){\left(\frac{ϵ}{0.01}\right)}^{-1},\end{array}$$(14)

providing that the acceleration time scale depends only on the energy fraction that is transferred to the fast-mode waves.

The re-acceleration of cosmic rays will inevitably damp the fast-mode waves and affect the width of the turbulent region, L_fm. The latter quantity can be estimated by equating the energy density of the fast-mode waves, U_fm, and the total energy per volume that went into cosmic rays. Its value can be obtained from the energy transfer rate from waves to particles, ${\dot{E}}_{\mathrm{tr}}$, and the time period that particles spend in the turbulent region interacting with fast-mode waves, Δt, providing

$$\begin{array}{c}\hfill U_{\mathrm{fm}}={\dot{E}}_{\mathrm{tr}}·\mathrm{\Delta }t=\frac{U_{\mathrm{cr}}}{{\tau }_{\mathrm{acc}}}·\frac{L_{\mathrm{fm}}}{u_2}.\end{array}$$(15)

Here, U_cr is the cosmic-ray energy density in the immediate post-shock region induced by the stochastic re-acceleration. Combining Eqs. (13) and (15) we finally obtain

$$\begin{array}{c}\hfill L_{\mathrm{fm}}\approx (5\times {10}^{13}\mathrm{cm})\left(\frac{V_{\mathrm{sh}}}{1000\mathrm{k}m{s}^{-1}}\right)\left(\frac{U_{\mathrm{th}}}{U_{\mathrm{cr}}}\right).\end{array}$$(16)

Summarizing, the thickness of the re-acceleration region is limited by damping of the turbulence caused by the re-acceleration of particles. Therefore, in our method, U_cr from Eq. (16) is coupled to the intermediate results from Eq. (6) for the immediate post-shock area. Calculating the energy density of cosmic rays, U_cr, we do not take the contribution from electrons into account because in our model their energy is negligible compared to that of protons.

The physically essential quantity for stochastic re-acceleration is the amount of energy available for it, here in the form of fast-mode waves and described by the parameter ϵ. The size of the turbulence region, L_fm, follows from the acceleration timescale, τ_acc, and only their ratio is relevant for the resulting particle spectra. The reason for this is that re-acceleration, once efficient, becomes the main damping mechanism for the waves and hence ceases when the energy supply is exhausted. More details on stochastic re-acceleration in SNRs can be found in the Pohl et al. (2015).

2.4. Radiative processes

To account for the entire SED from Tycho, we calculate synchrotron radiation of electrons and γ-ray radiation resulting from leptonic and hadronic interactions. The synchrotron emission is calculated following Blumenthal & Gould (1970). Care must be exercised to properly account for magnetic field fluctuations in calculating the synchrotron emissivity function. The standard calculation, although not applicable to the turbulent magnetic field, assumes by default a delta function for the probability distribution of local magnetic-field amplitudes. For the turbulent magnetic-field, the amplitude distribution is dispersed compared to the standard case; its exact form is unknown. However, modified emissivities can be found in the literature for exponential (Zirakashvili & Aharonian 2010) and power-law (Kelner et al. 2013) distributions. In this work we use the Gaussian distribution of magnetic-field amplitudes and the corresponding synchrotron emissivity function (Pohl et al. 2015). Given the low gas density, nonthermal bremsstrahlung (Blumenthal & Gould 1970) yields too low a flux for it to be relevant, and so leptonic gamma-ray emission is solely provided by inverse Compton scattering, for which we consider the microwave background as the target photon field (Blumenthal & Gould 1970; Sturner et al. 1997). Tycho shows evidence for infrared emission (Douvion et al. 2001). However, its contribution for the overall inverse-Compton spectrum was found by Acciari et al. (2011) to be negligible. Hadronic γ-rays are the result of decays of neutral pions and other secondaries produced in interactions of cosmic rays with nuclei of the ISM. To calculate its spectrum we follow the method from Huang et al. (2007).

3. Results

The method to solve the transport equation (Eq. (6)) separately for electrons and protons is described in Telezhinsky et al. (2012). We emphasize that, in contrast to previous attempts to model Tycho, we follow the full temporal evolution of the remnant starting at an age of 25 years. However, for conciseness, we show and discuss results for only the current age of Tycho. We find that the reverse-shock contribution to the particle spectrum in Tycho is negligible in the framework of our modeling at an age of 440 years (Telezhinsky et al. 2012), in full agreement with Warren et al. (2005). Therefore we do not discuss it further.

We present two models for the multiwavelength emission of Tycho, which both adequately fit the SED. Model I allows for the weakest magnetic field in the immediate downstream region of the shock, B₂ = 150 μG, which is compatible with the entire γ-ray flux observed with Fermi-LAT in GeV-range and VERITAS in TeV-band (Archambault et al. 2017). In Model I, the resulting γ-ray flux consists of both leptonic and hadronic components. The magnetic profile deeper inside the remnant is determined by advection of frozen-in magnetic field and corresponds to the MHD solution for negligible magnetic pressure and energy density. As we shall demonstrate, this model fails to explain radio and X-ray intensity profiles.

Therefore we present a second model, Model II, involving magnetic field damping and derive constraints on the magnetic field in Tycho. Generic technical details can be found in Sect. 3.1 where we introduce Model I; these also apply to the following sections listing the results for Model II.

3.1. Model I: Moderate advected magnetic field

The efficacy of stochastic acceleration is fully determined by three parameters: m and p₀, and the energy fraction of the plasma transferred to the turbulence, ϵ. The acceleration timescale, τ_acc, and the width of the turbulent region, L_fm, are then self-consistently calculated according to Eqs. (14) and (16). The parameter p₀ provides the critical limit up to which the momentum diffusion coefficient D_p is energy-independent (cf. Eq. (10)). It should be derived from exact scattering theory, which is nevertheless beyond the scope of this work. We tested several values of p₀ and found that for p₀ ≳ 10⁻³ m_pc the resulting curvature of the radio synchrotron spectrum becomes too strong to remain in agreement with radio data. Therefore, in this work we adopt p₀ = 9.3 × 10⁻⁴m_pc, which corresponds to an electron energy of 500 keV. The power-law index, m, may vary up to 25% (to stay in agreement with observations) and can be compensated by a suitable choice of energy-conversion factor, ϵ, implying that these quantities exhibit degeneracy to a certain degree. In the following calculations, we fix m = 0.25 and ϵ = 0.027 (corresponding to 2.7% of the thermal energy density), which adequately reproduce the observed radio spectrum. The acceleration time scale is determined by the energy-conversion factor, which for ϵ = 0.027 provides τ_acc ≈ 2.3 yr.

The thickness of the region in which stochastic re-acceleration is operating is fairly small. For 2.7% of the thermal energy density of the post-shock plasma transferred to the turbulence, it results in L_fm ∼ 10¹⁵−10¹⁶ cm, as illustrated in Fig. 2. The width of the re-acceleration region decays with time as the SNR shock slows down and thus the energy density in the turbulence decreases. At the age of 440 years the width of the turbulent region comprises only ∼3 × 10⁻⁴R_sh. Still, the contribution to the particle spectra is substantial.

Fig. 2. Energy density in fast-mode waves (red solid line) and width of the turbulent re-acceleration region (blue dashed line), both as a function of time for Model I.

Figure 3 shows the differential number density for electrons and protons as obtained from solving Eq. (6). As a result of stochastic re-acceleration, the spectra strongly deviate from the canonical solution, N ∝ p⁻², expected for DSA in the test-particle limit. For comparison, a standard case with stochastic re-acceleration turned off but the same remaining parameters is depicted in Fig. 4. Distinct bumps for both particle species at low energies can be seen in Fig. 3, with concave tails that extend up to the maximum momenta of the spectra, ∼10³m_pc. Deviations between electron and proton spectra at lower momenta result from the different injection criteria for the particles. Since we use the thermal leakage model, electrons and protons feature injection momenta differing by a factor of ∼40 due to $p_{\mathrm{inj},i}\propto \sqrt{m_i}$, as seen from Eq. (8). In other words, we do not explicitly treat electron acceleration, such as for example by stochastic shock drift acceleration below 100 MeV (Katou & Amano 2019), and subsume the entire acceleration from suprathermal to very high energies under DSA, as mentioned above. Electrons and protons at energies between the thermal peak and the injection threshold, which are a factor 4.45 higher in momentum, are considered part of the thermal bulk and are ignored. The electrons between their injection threshold and the proton injection momentum are very numerous and with stochastic re-acceleration can form a larger bump in the particle spectrum than would be observed for protons.

Fig. 3. Number density of protons (red solid) and electrons (green dashed) for Model I involving the weakest advected magnetic field with B2 = 150 μG. The proton number density is multiplied with the factor Ke/p = 1/85.

Fig. 4. Number density of protons (red solid) and electrons (green dashed) for Model I with turned off stochastic re-acceleration. A comparison with Fig. 3 illustrates the quantitative effect of the second-order acceleration in cosmic-ray spectra.

The total energy that went into electrons, E_tot, e ≈ 6.8 × 10⁴⁷ erg, is marginal compared to that of protons, E_tot, p ≈ 2.7 × 10⁴⁹ erg. Reasons for this are the small rest mass of electrons and different injection efficiencies η_e = 9.4 × 10⁻⁶ and η_p = 2.4 × 10⁻⁵, which are chosen to fit the entire SED. The electron injection efficiency simultaneously accounts for the maximum IC peak consistent with the γ-ray data and sufficient radio emission for B_d = 150 μG (amplification factor α ≈ 9), while injection of protons provides an adequate hadronic contribution in the GeV range. The injection efficiencies determine the electron-to-proton flux ratio, K_e/p ≡ N_e/N_p ≈ 1/85 in the range (10−10³)m_p c.

The proton spectrum cuts off at the maximum momentum, p_max ≈ 10⁴m_p c ≃ 10 TeV/c, which is limited by the age of the remnant and particle diffusion in the upstream region. The latter process is determined by two factors: the magnetic precursor which flattens the cosmic-ray precursor, discussed in Sect. 2.2, and the gyrofactor, ξ, which in conjunction with the magnetic-field strength determines the diffusive transport of particles. Since the electrons additionally experience effective synchrotron losses, their maximum momentum, ∼3 TeV/c, is lower than the limit (10 TeV/c) imposed by the spatial diffusion and the age of the remnant. The exact shape of electron and proton cutoffs is shown in Fig. 5, which allows a direct comparison of the cutoff forms. As one can see, the electron density is steeper than that of protons due to effective synchrotron losses. We note that the shape of the spectral cutoff in our approach deviates from the simple exponential and super-exponential forms (e.g., Blasi 2010) due to the full time-dependency of our method. In contrast to the previous works on Tycho, our approach includes the time-dependent transport equation and hydrodynamics. Indeed, the maximum cosmic-ray energy in the age-limited case scales linearly with time, and as the square of the shock velocity, $V_{\text{sh}}^2$. As the shock of the remnant slows down in our simulations, the maximum energy of cosmic rays grows more slowly than one would expect for an approach with a constant shock velocity. As a result, the final particle spectrum in our work, which is qualitatively an average of the instantaneous spectra, gives a sharper cutoff than predicted within the hydrodynamics in steady-state.

Fig. 5. Number density of protons (red solid) and electrons (green dashed) for Model I as in Fig. 3 zoomed to the cutoff region.

Instead of showing the entire SED at one singe plot, we present figures for particular energy bands (radio, X-ray and γ-ray) to provide a detailed comparison of data and models. The synchrotron emission from electrons in the radio and microwave bands for Model I is presented as a red solid line in Fig. 6. The radio spectrum for Tycho is frequently considered to be distorted in response to shock modification by cosmic rays (Reynolds & Ellison 1992; Völk et al. 2008). Figure 6 clearly demonstrates that this is not the only viable interpretation: stochastic re-acceleration in the downstream region of a test-particle shock can reproduce the observed radio spectrum without invoking nonlinear effects.³ It was noted earlier (Planck Collaboration Int. XXXI 2016) that the radiation from Tycho in the microwave band consists of at least two components: synchrotron emission from the nonthermal electrons and thermal dust emission. The latter process is responsible for the sharp rise in flux above ∼30 GHz. We account for the thermal dust emission by calculating the gray-body radiation with a temperature of 25 K and a normalization chosen to fit the flux density measured with Herschel (Gomez et al. 2012). Thus, the red line in Fig. 6 represents the sum of the synchrotron and thermal spectra. While the slightly concave part of the radio spectrum below ∼30 GHz is dominated by the synchrotron emission from electron population shaped by the stochastic re-acceleration, the steep flux increase above ∼30 GHz results from the thermal dust emission.

Fig. 6. Synchrotron plus gray-body emission in the radio and microwave bands for Models I and II. The radio (black error bars) and microwave (black triangles) data are taken from Reynolds & Ellison (1992) and Planck Collaboration Int. XXXI (2016).

The X-ray emission is generated by electrons beyond the cutoff energy, E_max, e ≈ 3 TeV, via synchrotron radiation; it is presented for Model I as a red solid line in Fig. 7. To achieve a good fit for the fixed magnetic field (B₂ = 150 μG) and electron injection efficiency (η_e = 9.4 × 10⁻⁶) we adapt the spatial diffusion coefficient parameter and set ξ = 10. As already discussed in Sect. 2.4, for our calculation of the synchrotron emission we use a Gaussian distribution function for the amplitudes of the turbulent magnetic field. The differences to the standard formula are mostly seen at the cutoff of the synchrotron emission. Indeed, the standard emissivity function produces the steepest slope in the spectral tail while any extended distribution smears the spectral cutoff and causes spectral hardening in the X-ray band. Nevertheless, the predicted X-ray spectrum above 10 keV is softer than that observed. This can be attributed to the spherical symmetry of our model geometry. Lopez et al. (2015) analyzed 66 regions across Tycho and find that their X-ray emission exhibits a varying roll-off energy, which is defined as

$$\begin{array}{c}\hfill E_{\mathrm{rolloff}}\simeq 7\mathrm{eV}\left(\frac{B_2}{100\mu \mathrm{G}}\right){\left(\frac{E_{\max ,\mathrm{e}}}{\mathrm{TeV}}\right)}^2.\end{array}$$(17)

Fig. 7. Synchrotron emission in the X-ray band for Models I and II with Gaussian distribution of magnetic-field amplitudes. Experimental data at 90% confidence level are taken from Tamagawa et al. (2009).

As the magnetic field and the maximum energy of electrons may vary across Tycho’s perimeter, so do corresponding synchrotron spectra. Integration over the individual regions, as in Lopez et al. (2015), results in a harder total spectrum, because the variations in the roll-off energy will invariably smear out the cutoff. This hardening cannot be accounted for in our model due to the 1D geometry, and hence our X-ray spectrum is softer than that observed by Suzaku. We note that since we use a Gaussian distribution function to calculate the synchrotron emission, its cutoff deviates from the usually assumed exponential profile and thus Eq. (17) is no longer valid in our case.

The minimum downstream magnetic field that provides the maximum inverse-Compton contribution compatible with the γ-ray data and sufficient flux contribution in the radio range on account of stochastic re-acceleration is ∼150 μG. A magnetic field any weaker than this would imply an overshooting of γ-ray emission at ∼100 GeV induced by the inverse Compton process. Nevertheless, it is not able to account for the GeV-scale emission and hence the hadronic channel is required. The resulting γ-ray spectra of Tycho and the corresponding γ-ray data are given in the top panel of Fig. 8. The pion bump is represented by the blue dotted and the inverse-Compton peak by the green dashed lines, respectively. The total γ-ray spectral distribution (red solid line) is rather flat with spectral index Γ ≈ 2. The impact from the stochastic re-acceleration of protons is hardly visible in the pion bump. The reason is a relatively small energy fraction in the fast-mode waves. In our model for Tycho, parameters relevant for the stochastic re-acceleration are dictated by the radio data. For other SNRs, their value can differ, potentially resulting in more efficient re-acceleration of protons and consequentially a softer hadronic spectrum. Thus, stochastic acceleration of protons may provide an alternative explanation for the softening of hadronic emission, as opposed to high-density structures in the ambient medium (Berezhko et al. 2013; Morlino & Blasi 2016), if enough energy is available for re-acceleration. Nevertheless, the impact from stochastic re-acceleration is more prominent for electron- than proton spectra. The cutoff in the leptonic and hadronic γ-ray contributions is linked to that in the synchrotron spectrum via the gyrofactor, ξ, and therefore the former is constrained by X-ray data. Therefore, in our approach, we do not specifically adjust the model parameters to fit the hadronic cutoff. Nevertheless, it shows a remarkable match with the observed VERITAS data. The corresponding parameters for Model I are summarized in the first row of Table 1. The amplification efficiency (i.e., the ratio between cosmic-ray pressure and the amplified magnetic-field pressure in the immediate upstream region) is roughly ∼74 for 440 years.

Fig. 8. Calculated γ-ray emission from Tycho in comparison with measurements of Fermi-LAT (black solid) and VERITAS (black circles) from Archambault et al. (2017). The top panel is for Model I involving a moderate transported magnetic field and the bottom panel presents results for Model II, which includes damping of magnetic field.

As already mentioned in Sect 2.4, the contribution from nonthermal bremsstrahlung is negligible for the density of gas at Tycho. Therefore, it is not shown in Fig. 8 and we do not discuss it further.

To further test Model I we consider the spatial profiles of radio and X-ray synchrotron radiation at 1.4 GHz and 1 keV, respectively. The results are depicted in the top panel of Fig. 9. The radio (Slane et al. 2014) and X-ray (Cassam-Chenaï et al. 2007) brightness profiles are extracted from the western rim of the remnant and normalized to their peak values. The measured X-ray brightness is strongly peaked at the shock and rapidly decreases towards the contact discontinuity. The radio profiles are slightly wider but still exhibit a narrow structure close to the shock. The enhancement of the emission towards the interior that is evident in the radio profiles might be attributed to the afore-mentioned Rayleigh-Taylor distortions operating in the vicinity of CD. As seen in the figure, the predictions of Model I involving the smallest possible advected magnetic field are in disagreement with the data. The model can explain neither the narrow X-ray rims nor the structure of the radio emission.

Fig. 9. Radial intensity profiles for Model I (top panel) and Model II (bottom panel). X-ray emission at 1 keV (blue dashed line) was observed with Chandra (Cassam-Chenaï et al. 2007), and radio data at 1.4 GHz (red solid line) were taken from Slane et al. (2014). Following Slane et al. (2014) the radio data were slightly shifted to account for the expansion of the remnant. The gray dotted line in the bottom panel represents a radio profile produced without magnetic-field damping.

There are at least two potential solutions to this problem that we shall explore in the following section. Turbulent magnetic-field damping would affect the synchrotron emissivity and thus create narrow structures close to the shock. Alternatively, a very strong magnetic field at the shock would impose strong synchrotron losses and thus limit the width of the rims.

3.2. Model II: High damped magnetic field

As Model I involving a weak advected magnetic field fails to reproduce the observed intensity profiles of the synchrotron emission in both the radio and the X-ray bands, we now introduce damping of the turbulent magnetic field. Therefore, we set the magnetic-field profile described by Eq. (5). In this scenario the magnetic field strength immediately downstream of the forward shock must be larger than for Model I, because damping suppresses synchrotron emission from the far-downstream region. As inverse Compton radiation is produced wherever high-energy electrons reside, an overproduction of photons in the γ-ray band would arise unless the magnetic field is scaled up. Especially when combined with stochastic re-acceleration as an explanation for the radio data, the magnetic-field damping requires an extensive increase of the overall magnetic field. As reflected by the decreasing width of the turbulent region, the stochastic re-acceleration of particles is more efficient at the early stages of a SNR. As the re-accelerated electrons are advected away from the shock, they mostly experience the damped magnetic field. While the total synchrotron emission depends on the absolute values of the immediate post-shock field strength, B₂, and the residual level B₀, the shape of the radio filaments is determined (apart from the damping scale l_d) by the ratio of B₂ and B₀. We find that the immediate post-shock field strength and the residual level must be at least B₂ = 330 μG (amplification factor α = 20) and B₀ = 70 μG, respectively, in order to maximize the filling factor, and thus to produce sufficient overall synchrotron emission, and to fit the radio filaments simultaneously. Model II is based on these minimum values, and results are displayed as green dashed lines in Figs. 6 and 7, as well as in the bottom panels of Figs. 8 and 9.

Particle acceleration and propagation is modeled following the same procedure as in Model I, but some of the parameters are slightly adjusted. The second row of Table 1 lists the relevant model parameters, which now include the damping length, l_d, and the residual field level, B₀. We note that the damping length, l_d, provides spatial characteristics for the turbulence responsible for the magnetic-field amplification. Likewise, the width of the re-acceleration region, L_fm, represents the damping scale for the fast-mode waves. As both quantities are associated with different types of turbulence, they are independent of each other. At 440 years, the amplification efficiency for the second model is ∼14.

Along the shock surface certain variations in the parameter values are to be expected, as not all filaments are the same (cf. Ressler et al. 2014). Also, Tycho is not a perfectly spherically symmetric SNR and the projection effects may impose a bias. Therefore, our choice for the residual field level and the damping scale rather provide a reasonable value for the order of their magnitudes as in this work we consider only one particular rim of Tycho. Fitting the Suzaku data with an enhanced post-shock magnetic-field B₂ = 330 μG requires a moderate increase of the gyrofactor compared to Model I (ξ = 16).

Damping quenches the synchrotron emissivity in the deep downstream region. While the suppression is achromatic for power-law electron spectra, i.e., at frequencies below the cutoff in the synchrotron spectrum, it becomes progressively stronger beyond the roll-off energy. A competing process for quenching the emissivity at the roll-off frequency is energy losses preventing the propagation of electron from the shock to the deep downstream region. The synchrotron loss time can be expressed in terms of the magnetic-field strength and the energy of photons, E_sy, that the electrons would typically emit

$$\begin{array}{c}\hfill {\tau }_{\mathrm{loss}}\simeq 50\mathrm{yr}{\left(\frac{B_2}{100\mu \mathrm{G}}\right)}^{-3/2}{\left(\frac{E_{\mathrm{sy}}}{\mathrm{keV}}\right)}^{-1/2}.\end{array}$$(18)

The distance electrons travel during the loss time is roughly equal to the thickness of the filaments. The particle transport is governed by diffusion and advection. The latter process dominates at lower energies and its length scale is given by

$$\begin{array}{cc}& l_{\mathrm{adv}}=u_2{\tau }_{\mathrm{loss}}=\frac{V_{\mathrm{sh}}}{r_{\mathrm{sh}}}{\tau }_{\mathrm{loss}}\hfill \\ \hfill & \simeq 2\times {10}^{17}\mathrm{cm}\left(\frac{V_{\mathrm{sh}}}{5000{\mathrm{km}\mathrm{s}}^{-1}}\right){\left(\frac{B_2}{100\mu \mathrm{G}}\right)}^{-3/2}{\left(\frac{E_{\mathrm{sy}}}{\mathrm{keV}}\right)}^{-1/2}.\hfill \end{array}$$(19)

At high energies, diffusion with Bohmian energy scaling becomes the dominant propagation process. The corresponding distance is energy-independent and is given by

$$\begin{array}{cc}& l_{\mathrm{diff}}=\sqrt{D_r{\tau }_{\mathrm{loss}}}\hfill \\ \hfill & \simeq {10}^{17}\mathrm{cm}\sqrt{\xi }{\left(\frac{B_2}{100\mu \mathrm{G}}\right)}^{-3/2}.\hfill \end{array}$$(20)

By setting Eq. (19) equal to Eq. (20) and solving the resulting expression for the photon energy, one can find the critical value, where the transition from advection-dominated transport to diffusion-dominated transport occurs:

$$\begin{array}{c}\hfill E_{\mathrm{sy},\mathrm{c}}\simeq \frac{4\mathrm{keV}}{\xi }{\left(\frac{V_{\mathrm{sh}}}{5000{\mathrm{km}\mathrm{s}}^{-1}}\right)}^2.\end{array}$$(21)

With ξ = 16 and V_sh = 4100 km s⁻¹ we obtain E_sy, c ≈ 0.2 keV for Model II, implying that both advection and diffusion impact the distance covered by electrons that account for the 1 keV rims. Estimating the propagation length for the electrons that radiate 1 keV photons, we find l_adv ≃ 3 × 10¹⁶ cm and l_diff ≃ 7 × 10¹⁶ cm, and therefore diffusive transport is more important, but advective transport is not negligible. Accounting for both advective and diffusive terms, the total transport length that the electrons travel before expending their energy is then given by (Parizot et al. 2006)

$$\begin{array}{c}\hfill l_{\mathrm{loss}}={(\sqrt{\frac{u_2^2}{4D_r^2}+\frac{1}{D_r{\tau }_{\mathrm{loss}}}}-\frac{u_2}{2D_r})}^{-1}\simeq {10}^{17}\mathrm{cm},\end{array}$$(22)

where the last expression applies to an electron that emits 1-keV photons and Model II. The underlying assumption of a constant magnetic-field strength is violated though, as the synchrotron energy losses decrease over a length scale 5 × 10¹⁶ cm. Therefore, the effective loss-length of electrons is much larger than 10¹⁷ cm, and the synchrotron filaments are largely shaped by magnetic-field damping.

Magnetic-field damping is clearly needed to account for the radio profiles. Radio-emitting electrons have energy-loss times far in excess of the age of SNRs, and the radio rims cannot arise from synchrotron losses. The model with magnetic-field damping (Model II) fits the spectral data reasonably well. To be noted from Fig. 9 is that an effective damping length of l_d = 0.01 R_sh ≃ 10¹⁷ cm can indeed reproduce the sharply peaked radio profiles in the shock vicinity, but somewhat underpredicts the radio intensity in the deep downstream region where Rayleigh-Taylors fingers from the contact discontinuity may provide magnetic field amplification (Jun & Norman 1995; Björnsson & Keshavarzi 2017) that is not included in our model. In contrast, a radio emission profile calculated for nondamped magnetic field (gray dotted line in bottom panel of Fig. 9) obviously contradicts the observed data. As a relatively high magnetic-field value is required to maintain the total radio and γ-ray data, the width of the X-ray rims is inevitably affected by synchrotron losses. Still, magnetic-field damping is more important for the formation of the X-ray filaments since l_loss >  l_d. Without magnetic-field damping, the synchrotron losses for the post-shock magnetic field B₂ = 330 μG are able to produce the thin X-ray rims, but fail to form the radio profiles. When introduced as a natural explanation for the radio profiles, magnetic-field damping moreover becomes the dominant process for the production of the X-ray filaments.

We conclude that magnetic-field damping is essential for the radio filaments. Furthermore, for the minimum downstream magnetic field that can explain the complete observed data, B₂ = 330 μG, magnetic-field damping is more important for the X-ray profiles than radiative losses. Nevertheless, for 1-keV emission synchrotron losses provides a subordinate additional process for X-ray filaments formation. Simultaneously the width of the radio rim is solely determined by the magnetic-field damping. In agreement with our results, the joint X-ray and radio analysis of Tran et al. (2015) reveals magnetic-field damping to be the preferable scenario. Both studies find similar damping lengths: 1–2% of the SNR radius.

3.3. Cosmic-ray pressure

Next, we verify that the test-particle approximation is valid for Model II. To this aim, we calculate the cosmic-ray pressure, given by

$$\begin{array}{c}\hfill P_{\mathrm{cr}}(r,t)=\frac{c}{3}{\displaystyle \int }N(r,p,t)\frac{p^2\mathrm{d}p}{\sqrt{p^2+{(m_{p}c)}^2}},\end{array}$$(23)

where m_p is the rest mass of a proton and N(r, p, t) the differential proton number density. As mentioned before, the pressure exerted by the nonthermal electrons is negligible in our model. The relative-velocity change of the plasma in the shock rest-frame in a particular region between r₁ and r₂ caused by the particle pressure is

$$\begin{array}{cc}& \frac{\delta u}{u}=-\frac{1}{u^2\rho }{\displaystyle {\int }_{r_1}^{r_2}}\mathrm{d}x\frac{\partial P_{\mathrm{cr}}(x)}{\partial x}\hfill \\ \hfill & =\frac{P_{\mathrm{cr}}(r_2)-P_{\mathrm{cr}}(r_1)}{P_{\mathrm{flow}}}\equiv \frac{\delta P_{\mathrm{cr}}}{P_{\mathrm{flow}}}.\hfill \end{array}$$(24)

Here ρ denotes the density, u the velocity, and P_flow = ρu² the dynamical pressure of the plasma. Equation (24) is universally valid and holds in any region of the plasma flow.

According to Kang & Ryu (2010), the test-particle approximation is justified as long as the cosmic-ray pressure in the precursor does not exceed 10% of the shock ram pressure. Consequently, according to Eq. (24) the test-particle regime requires δu/u = δP_cr/P_flow ≤ 0.1 in the upstream region. Figure 10 shows the cosmic-ray pressure normalized by the shock ram pressure, $P_{\text{sh}}={\rho }_1{u}_1^2$, as a function of the position, r/R_sh, for different times. To be noted from the figure is that the relative velocity change in the upstream region at 440 years is δu/u = δP_cr/P_sh ≈ 0.021, and hence clearly below the test-particle threshold. Similarly, at the earlier stages, the cosmic-ray pressure at the shock does not exceed δu/u <  0.1, verifying that the dynamic cosmic-ray feedback is indeed negligible in our approach.

Fig. 10. Cosmic-ray pressure, Pcr, normalized by the shock ram pressure, Psh, for Model II (damped magnetic field with B2 = 330 μG) for different times.

Behind the shock wave, the pressure of particles is significantly enhanced due to stochastic re-acceleration. In contrast to standard DSA, where the maximal cosmic-ray pressure occurs at the shock, in our model it reaches its maximum in the immediate downstream region. The corresponding position is determined by the current width of the turbulent zone where the fast-mode waves operate. Taking into account that the dynamic pressure in the post-shock region is a factor four lower than the shock ram pressure in the shock rest-frame, we estimate for the total velocity change at 440 years: δu/u ≈ 0.01. We neglect the dynamical feedback from the re-accelerated particles in this work. Nevertheless, this effect may be of interest for future investigations. Indeed, when strong enough, the stochastic acceleration can have a significant effect on the post-shock structure. The cosmic-ray pressure in the immediate downstream region enhanced by the stochastic acceleration has to press the plasma away from the re-acceleration region and force it to accumulate directly behind the shock. This process increases the plasma density and the associated magnetic-field strength at the shock of the SNR. It is clear that this phenomena fundamentally differs from the standard dynamical reaction of particles described in the NLDSA theory. Summarizing, we conclude that we find a viable self-consistent scenario for Tycho, represented by Model II.

3.4. Hadronic model

In the previous section we present a lepto-hadronic scenario that requires the lowest possible magnetic field in the remnant compatible with radio filaments and the SED. More efficient magnetic-field amplification suppresses the inverse-Compton component until the γ-ray spectrum becomes purely hadronic. In our approach, where we use magnetic-field damping to fit the radio filaments, the magnetic-field values have to be at least B₀ ≈ 120 μG and B₂ ≈ 600 μG. The total magnetic energy in the remnant is still moderate, ∼1.1 × 10⁴⁹ erg, on account of the magnetic-field damping.

As discussed in Sect. 2.2, the separate treatment of magnetic field and hydrodynamics for the ambient density of n_H = 0.6 cm⁻³ is justified for post-shock magnetic field below 400 μG. Thus, for the hadron-dominated case, the dynamical feedback from the magnetic field becomes significant. The corresponding pressure, with the Alfvénic Mach number M_A ≈ 9.6, would lower the shock compression ratio to r_sh ≈ 3.87 and thus soften the particle and radio spectra indices to s ≈ 2.05 and σ ≈ 0.52, respectively. As the electron spectrum becomes slightly softer on account of the impact of the magnetic-field, the contribution from the stochastic re-acceleration becomes less necessary. Hence, in general, as competing explanation for the softening of the spectra, very high magnetic field inevitably decreases the energy fraction converted to turbulence, ϵ.

Due to the very high magnetic field, X-ray filaments at 1 keV would be primarily governed by the radiative losses, with corresponding propagation scale l_loss ≃ 4 × 10¹⁶ cm. Thus, the brightness profile would lack the X-ray flux due to extreme synchrotron losses. Nevertheless, as pointed out above, the radio as well the X-ray rims can strongly vary with position along the perimeter due to the natural asymmetry of the remnant. As the electrons emitting the synchrotron radiation in the radio range do not experience significant losses, the radio profiles do require magnetic-field damping.

In summary, we conclude that a purely hadronic model is possible for Tycho, but requires an elaborate and cautious treatment which, among other effects, includes the dynamical feedback from magnetic-field pressure. Nevertheless, the lepto-hadronic scenario, referred to as Model II, is able to explain the broadband observations of Tycho in a satisfactory manner.

4. Summary and conclusion

Here, we conduct extensive multiwavelength modeling of Tycho. For the very first time, we account for stochastic re-acceleration in the downstream region of the forward shock, which provides a consistent explanation for the soft particle spectra without resorting to Alfvénic drift (Völk et al. 2008; Morlino & Caprioli 2012; Slane et al. 2014) and its inherent problems, or deducing the particle spectral index from the radio observations (Atoyan & Dermer 2012; Zhang et al. 2013; Caragiulo & Di Venere 2014). As discussed in the introduction, although it has been widely used in the global models for various SNRs, we find the concept of Alfvénic drift contradictory. We present instead a new approach that adds diffusion in momentum space to the standard DSA approach. Furthermore, in our method both hydrodynamics and particle acceleration are fully time-dependent.

The stochastic acceleration of particles in the immediate downstream region is assumed to arise from the fast-mode turbulence, which is supplied by the energy of the background plasma. We find that 2.7% of the thermal energy density of the downstream background plasma is sufficient to explain Tycho’s soft radio spectra. Simultaneously, the magnetic field is assumed to be amplified by streaming instabilities (Bell 2004; Lucek & Bell 2000) or turbulent dynamos (Giacalone & Jokipii 2007). Based on this, we present a self-consistent global model (Model II) that is able to accurately reproduce the observed radio, X-ray, and γ-ray emission, and simultaneously account for the nonthermal filaments in the radio and X-ray ranges. The radio filaments are generated due to magnetic-field damping, which is widely considered to allow a relatively low magnetic-field value inside a remnant. Combining this scenario with stochastic re-acceleration, we find that the minimum magnetic-field required to explain the entire observed dataset of Tycho is B₂ ≈ 330 μG. This value is similar to the results of Morlino & Caprioli (2012), Völk et al. (2008), and Zhang et al. (2013) although the underlying physical assumptions are quite different. We find that for this minimum magnetic-field strength of B₂ ≈ 330 μG, the X-ray filaments at 1 keV are primarily produced by magnetic-field damping, while the synchrotron losses play a secondary role. This finding is inconsistent with the work of Morlino & Caprioli (2012), who concluded that X-ray filaments shaped by magnetic-field damping are not possible for Tycho. An important criterion here is that the propagation length, and therefore the synchrotron loss scale for electrons radiating at 1 keV, is dominated by diffusion. For this reason, the diffusive transport of particles, which has been previously neglected in all global models of Tycho, must be taken into account to adequately model the X-ray filaments. We stress that an accurate modeling of the filaments is seamlessly tied to the determination of the post-shock magnetic field. Magnetic-field damping is additionally needed for its unique capability to explain the radio rims. In our model for Tycho, the magnetic-decay length is of the order of l_d ∼ 0.01 R_sh, which is consistent with the estimate of Tran et al. (2015). The total magnetic energy in the remnant for Model II is 3.4 × 10⁴⁸ erg. This value is moderate on account of magnetic-field damping, in spite of the efficient magnetic-field amplification at the shock.

In the framework of our model, we predict relatively inefficient Bohm diffusion, reflected in the value of the gyrofactor ξ ≈ 16. This value is required to consistently fit the synchrotron cutoff observed in the X-ray range. In line with most previous works on Tycho, we also conclude that acceleration of protons is required to explain the γ-ray flux observed by VERITAS and Fermi-LAT. Our research cannot account for a purely leptonic origin for the γ-ray emission as in Atoyan & Dermer (2012). We instead favor a mixed model with a composite flat spectrum (Γ ≈ 2). The electron-to-proton ratio for Model II is K_e/p ≈ 1/80 and the maximum energy for protons is E_max, p ≈ 10 TeV, since it is linked to the roll-off energy of the synchrotron emission via gyrofactor, ξ. This result falls below previously presented ∼500 TeV (Morlino & Caprioli 2012) and ∼50 TeV (Slane et al. 2014). The maximum energy of electrons is ∼3 TeV, which is similar to the value 5–6 TeV suggested by Zhang et al. (2013) and Caragiulo & Di Venere (2014). The total energy in protons in our model is E_tot, p ≈ 2.7 × 10⁴⁹ erg, implying that a few percent of the explosion energy went into cosmic rays. This value is significantly below the 16% claimed by Slane et al. (2014). We can use the test-particle approximation, as a relatively marginal energy fraction is transferred to particles and the cosmic-ray pressure at the shock does not exceed 2.1% of the shock ram pressure during the entire evolution of the remnant.

Furthermore, we explicitly show that NLDSA effects are not required, either to explain the hydrodynamic structure of Tycho or to produce its slightly concave radio spectrum. First, the hydrodynamic simulations with an ambient gas density of n_H = 0.6 cm⁻³ and canonical explosion energy of 10⁵¹ erg provide a decent fit to the observed radii and a reasonable remnant distance ∼2.9 kpc. Second, the radio spectrum is produced by synchrotron emission generated by electrons that are re-accelerated in the immediate downstream region. In general, the imprint left by stochastic re-acceleration is more prominent in electrons than in protons. Thus, future γ-ray observations that can successfully discriminate between leptonic and hadronic models for various SNRs may also be able to distinguish between NLDSA and stochastic re-acceleration scenarios.

We find that a purely hadronic model may also be possible for Tycho for the immediate post-shock magnetic-field B₂ ≈ 600 μG and far-downstream field B₀ ≈ 120 μG. However, we do not explicitly model the hadronic scenario in this work because the dynamical reaction of the magnetic field has to be taken into account, an effect that is not yet included in our method. Nevertheless, we favor a lepto-hadronic scenario (Model II) as it might better produce the X-ray filaments. Additionally, extremely efficient magnetic-field amplification is not required, as in the case of a purely hadronic model.

In our model for Tycho, the contribution of the stochastic re-acceleration in the proton spectrum is marginal due to the small energy fraction converted into downstream turbulence. For a larger amount of energy in fast-mode waves, the corresponding γ-ray spectrum has to be naturally softer than predicted by the standard DSA, which could provide a viable explanation for the observed soft γ-ray spectra of some SNRs. This subject is beyond the scope of this work, but could be of interest for future investigations.

Finally, we emphasize that the dynamical feedback on the background plasma from the cosmic rays re-accelerated in the immediate downstream region of the SNR is of great interest for future studies because it clearly differs from the classical nonlinear effects discussed in the literature.

In summary, we find that stochastic re-acceleration is indeed a promising and natural alternative to explain the soft spectra of Tycho, and potentially other SNRs.

Acknowledgments

We are very grateful to Asami Hayato for kindly providing us the Suzaku data. Furthermore, we acknowledge NASA and Suzaku grants, and support by the Helmholtz Alliance for Astroparticle Physics HAP funded by the Initiative and Networking Fund of the Helmholtz Association. VVD is supported by NSF grant 1911061.

References

All Tables

All Figures

	Fig. 1. Shock speed (blue solid line), radius (green dashed line), and temperature (red dotted line) as functions of time.
In the text

	Fig. 2. Energy density in fast-mode waves (red solid line) and width of the turbulent re-acceleration region (blue dashed line), both as a function of time for Model I.
In the text

	Fig. 3. Number density of protons (red solid) and electrons (green dashed) for Model I involving the weakest advected magnetic field with B2 = 150 μG. The proton number density is multiplied with the factor Ke/p = 1/85.
In the text

	Fig. 4. Number density of protons (red solid) and electrons (green dashed) for Model I with turned off stochastic re-acceleration. A comparison with Fig. 3 illustrates the quantitative effect of the second-order acceleration in cosmic-ray spectra.
In the text

	Fig. 5. Number density of protons (red solid) and electrons (green dashed) for Model I as in Fig. 3 zoomed to the cutoff region.
In the text

	Fig. 6. Synchrotron plus gray-body emission in the radio and microwave bands for Models I and II. The radio (black error bars) and microwave (black triangles) data are taken from Reynolds & Ellison (1992) and Planck Collaboration Int. XXXI (2016).
In the text

	Fig. 7. Synchrotron emission in the X-ray band for Models I and II with Gaussian distribution of magnetic-field amplitudes. Experimental data at 90% confidence level are taken from Tamagawa et al. (2009).
In the text

	Fig. 8. Calculated γ-ray emission from Tycho in comparison with measurements of Fermi-LAT (black solid) and VERITAS (black circles) from Archambault et al. (2017). The top panel is for Model I involving a moderate transported magnetic field and the bottom panel presents results for Model II, which includes damping of magnetic field.
In the text

Fig. 9. Radial intensity profiles for Model I (top panel) and Model II (bottom panel). X-ray emission at 1 keV (blue dashed line) was observed with Chandra (Cassam-Chenaï et al. 2007), and radio data at 1.4 GHz (red solid line) were taken from Slane et al. (2014). Following Slane et al. (2014) the radio data were slightly shifted to account for the expansion of the remnant. The gray dotted line in the bottom panel represents a radio profile produced without magnetic-field damping.

In the text

	Fig. 10. Cosmic-ray pressure, Pcr, normalized by the shock ram pressure, Psh, for Model II (damped magnetic field with B2 = 330 μG) for different times.
In the text