{'fraction_non_alphanumeric': 0.03743288291727104, 'fraction_numerical': 0.023258772934320687, 'mean_word_length': 4.234847475689357, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Molecular environment, reverberation, and radiation from the pulsar wind nebula in CTA 1

We estimate the molecular mass around CTA 1 using data from Planck and the Harvard CO survey. We observe that the molecular mass in the vicinity of the complex is not enough to explain the TeV emission observed by VERITAS, even under favorable assumptions for the cosmic-ray acceleration properties of the supernova remnant. This supports the idea that the TeV emission comes from the PWN. Here, we model the spectrum of the PWN at possible different stages of its evolution, including both the dynamics of the PWN and the SNR and their interaction via the reverse shock. We have included in the model the energy lost via radiation by particles and the particles escape when computing the pressure produced by the gas. This leads to an evolving energy partition, since for the same instantaneous sharing of the injection of energy provided by the rotational power, the field and the particles are affected differently by radiation and losses. We present the model, and study in detail how the spectrum of a canonical isolated PWN is affected during compression and re-expansion and how this may impact on the CTA 1 case. By exploring the phase-space of parameters that lead to radii in agreement with those observed, we then analyze different situations that might represent the current stage of the CTA 1 PWN, and discuss caveats and requirements of each one.


INTRODUCTION
CTA 1 (G119.5+10.2) was proposed as a supernova remnant (SNR) by Harris & Roberts (1960). Based on HI observations, Pineault et al. (1993) derived that its distance is 1.4±0.3 kpc. In X-rays, CTA 1 was detected for the first time by ROSAT (Seward, Schmidt & Slane 1995). The latter authors also reported the presence of a faint compact source in the central region, RX J0007.0+7302. Slane et al. (1997) detected non-thermal emission coming from it using ASCA, indicating the possible presence of a synchrotron nebula. Chandra observations done by Halpern et al. (2004) confirmed it, resolving a point-like source embedded in a compact nebula of 3" in radius, with a jet-like feature. Lin et al. (2012) performed a 105 ks observation with Suzaku and detected a ∼10' extended feature (∼ 4.1d1.4 pc), which was interpreted as the bow shock of the pulsar wind nebula (PWN). They fit the spectrum using a power-law with an index of 1.8 suggesting a slow cooling synchrotron scenario.
Based on X-ray thermal emission, Slane et al. (1997) suggested CTA 1 has an age of ∼20 kyr. Later, Slane et al. (2004) improved on this estimate based on a Sedov expansion model, and proposed that the age is ∼13 kyr. The latter value is in agreement with the characteristic age of the pulsar (∼ 13.9 kyr).
At high energies, Mattox et al. (1996) proposed that the EGRET source 3EG J0010+7309 (which lies in spatial coincidence with RX J0007.0+7302), was a potential candidate for a radio-quiet γ-ray pulsar. Brazier et al. (1998) also pointed out that this source was pulsar-like, but a search for γ-ray pulsations using EGRET data failed (Ziegler et al. 2008). Flux variability and correlation studies with the position of the SNR using EGRET data gave support to this interpretation (Romero, Benaglia & Torres 1999;Torres et al. 2003).
A radio-quiet pulsar in CTA 1 was finally discovered during the commissioning phase of the Fermi satellite (Abdo et al. 2008). X-ray pulsations from this source were detected shortly after (Caraveo et al. 2010;Lin et al. 2010). The pulsar in CTA 1 has a period of ∼316 ms and a spindown power of ∼ 4.5 × 10 35 erg s −1 (Abdo et al. 2012). No radio counterpart to RX J0007.0+7302 was identified, most likely due to beaming. No optical counterpart is known either (Mignani et al. 2013).
A few models of the spectral energy distribution (SED) of the nebula in CTA 1 were presented. Zhang, Jiang & Lin (2009) proposed a model using the SNR radio flux as if it was coming entirely from the PWN, an assumption that may be considered as hardly tenable. Nevertheless, they expected that the PWN in CTA 1 could be detected with the VERITAS array.
The TeV detection of a source at the position of the CTA 1 was confirmed by Aliu et al. (2013), who reported an extended source with dimensions of 0.3 • × 0.24 • (7.3d1.4 × 5.9d1.4 pc) at 5' from the Fermi γ-ray pulsar PSR J0007.0+7303. A hint for a detection of an extended GeV source in the off-pulse emission of PSR J0007.0+7303 was reported by Abdo et al. (2012), using 2 years of Fermi-LAT data. However, only upper limits were imposed. The use of the newest response functions of the Fermi satellite, together with several years of additional data, apparently does not confirm this earlier hint and rather verifies the upper limits to the GeV emission of the CTA 1 PWN (Li et al. (Fermi-LAT Collaboration) 2016, in preparation).
More recently, Aliu et al. (2013) produced a 10 kyr model based on the PWN setup described in Gelfand, Slane & Zhang (2009), using a low ejected mass of 6.1M⊙, a radius of ∼7 pc, and a magnetic fraction (see Section 3 below) of 0.2 (implying a magnetic field of 6.3 µG). Torres et al. (2014) performed a population study of young (tage ∼ 10 kyr or less) TeV-detected PWNe using a time-dependent model. The fit for the PWN in CTA 1 led to an age of 9 kyr and an ejected mass of 10M⊙, a radius of 8 pc, and a magnetic fraction of 0.4. Such high magnetic fraction, as the one obtained by Aliu et al. (2013), contrast with all others magnetic fractions found for young PWNe, which are in the range 0.01-0.05. The nature of this difference is unclear. If the CTA 1 PWN is in, or has passed a reverberation phase, the synchrotron flux could perhaps be fitted with a lower magnetic fraction.
In this work, we include a detailed accounting of the dynamical evolution of the PWN into the radiative model described by Martin et al. (2014) to study the effects of the reverberation phase in PWNe. We then perform a parameter space exploration for the PWN in CTA 1, varying the ejected mass, the age, and the interstellar medium (ISM) density to see whether it is possible to explain the features of the SED with a lower magnetic fraction. We shall also present a detailed study of the molecular environment of the CTA 1 nebula, in order to consider (and discard) possible alternative origins of the TeV radiation. We shall particularly analyze the possibility that the interaction of cosmic-rays with molecular clouds could give rise to the TeV source. With the latter aim, we shall present CO data from the Planck satellite and the Harvard survey of Dame, Hartmann & Thaddeus (2001).

THE MOLECULAR ENVIRONMENT IN
THE VICINITY OF CTA 1 The molecular content of a region can be estimated through the intensity of the 12 CO(1→0) line at 115 GHz (2.6 mm, see e.g., Dame et al. 1987). A survey of the intensity of the CO line, WCO, is provided in Dame et al. (1987) and Dame, Hartmann & Thaddeus (2001), hereafter referred to as DHT dataset. The data are released 1 in cubes of radial velocity, galactic latitude, and longitude. The radial velocity connects to the estimation of the distance. From WCO one can estimate the mass of the molecular material as follows: where µ is the mean weight, mH is the mass of the H nucleon, D is the distance, ∆Ω is the solid angle. The latter is taken out of the summation as it is the same for each pixel (squares of 0.25 • per side). XCO = NH2/WCO is a conversion factor defined from the assumption of proportionality of the CO content to the molecular content. Hereafter, recognizing that there is some uncertainty in this value, we will use µ = 2 and XCO = 2.0 × 10 20 cm −2 (K km/s) −1 . This is the average galactic value from the review of Bolatto, Wolfire & Leroy (2013). We shall find that our conclusions below do not depend on the value of XCO. Fig.  1 presents the data of Dame, Hartmann & Thaddeus (2001), dubbed DHT18 in the archive, cut at the location of CTA 1 and integrated in appropriate velocity ranges. We select the maps that were derived with the moment masking technique for noise suppression (Dame 2011). The left panel shows the WCO integrated along the entire velocity range. The moment masking technique is evidently suppressing as noise most of the region. The green countours overlaid in Fig. 1 are taken from the emission at 4850 MHz (GB6 survey, Condon et al. 1994 2 ). The region coincident with the VERITAS source (Aliu et al. 2013) is at the limit of the emission recognized as a real signal. The middle panel of Fig. 1 shows the same data set, but clipped in the velocity range −26 < v < −6, where v is given in units of km/s. This velocity range corresponds to the distance of CTA 1 and of the HI shell seen by Pineault et al. 1993 (v = −16.9 km/s). However, most of the data in this velocity range is tagged as noise and clipped away. Therefore, considering a circle centered at the VERITAS VHE emission centroid (l, b) = (119.6 • ,10.4 • ), with an extension of 0.25 • (recalling that the VERITAS ellipse is 0.3 • x 0.24 • ), we derive a very conservative upper limit on the mass content of the region using the WCO integrated over all the velocity range: the total mass amounts to MCO ∼ 1000M⊙.
We have checked these results using Planck satellite data. We use the Type-2 CO foregorund map of the second release of the Planck satellite data 3 . We shown in Fig.  1 the cut on the region, rebinned in pixels of 2" from an original resolution of 15". The Planck satellite data on CO can only be compared to the DHT dataset when integrated along the whole range of velocity. The Planck Collaboration warns that a factor 1.4 in normalization can be expected in their maps with respect to DHT, due to uncertainties in the Planck satellite observations and the calibration in the  Pineault et al. (1993). Right: Planck dataset. Images are on saturated scale to highlight the faint emission in the region of CTA 1. Second row: same as row above, but smoothed with a 3 pixel Gaussian kernel (1px=15" for DHT and 1px=2" for Planck).
DHT survey. We find only a 15% difference in the mass estimation from the DHT and Planck data, which given the caveats above and the differences in binning is to be considered acceptable. This difference is in any case irrelevant for the conclusion that follows. The predicted flux from a molecular cloud embedded in the sea of CR is Aharonian (1991); Gabici (2008); Pedaletti et al. (2013): where Eγ is expressed in TeV, the distance D in kpc, the mass M5 = 10 5 M⊙, and κ is the enhancement factor of CRs, assumed to be unity for passive clouds (i.e. a case where only the CR background is included) and larger in the presence of a nearby accelerator. We are using the PAMELA results by Adriani et al. (2011) for the Galactic cosmic-ray sea at energies > 200 GeV. We note that the factor k parametrizes our ignorance of the putative source, hence includes both normalization factor and could be considered as an enhancement produced by a hardening of the cosmic-ray spectrum too, at least for relevant CR energies of about 10 TeV. Considering the MCO calculated above, one would have F (E > 1 TeV) ∼ k 5 × 10 −16 cm −2 s −1 , requiring a very large en-hancement factor (k ∼ 1.5 × 10 3 ) to reach the VHE flux of F (E > 1 TeV) = 8.5 × 10 −13 cm −2 s −1 as reported in Aliu et al. (2013). This number can be compared with the case of the interacting SNR/molecular cloud complex W28, which also produces GeV and TeV sources, where the value k ranges from 13 to 32 (Aharonian et al. 2008 Thus, for the TeV emission to originate in a molecular cloud illumination (e.g., see Torres, Rodriguez Marrero & de Cea Del Pozo (2008)) the whole molecular mass along the line of sight should be concentrated right at the neighbourhood of CTA 1, and the latter should be acting as a very powerful accelerator. This is not supported by data (see middle panels of Fig. 1) nor by the arguments in the discussion above, and thus we believe that this possibility is untenable.

PWN MODEL DESCRIPTION
The evolution of the pair distribution in the PWN, N , as a function of the Lorentz factor γ at time t is described by where the term on the left hand side of the equation is the variation in time of the pair distribution. The second and the third term on the right hand side take into account the energy losses due to synchrotron, inverse Compton (IC), and Bremsstrahlung interactions, the adiabatic losses, and the escaping particles (we assume Bohm diffusion), respectively. A detailed description of these terms and the associated formulae can be found in Martín, Torres & Rea (2012); . Q(γ, t), the injection function, is to be assumed as a broken power law and the normalization term Q0(t) is computed using the spin-down luminosity of the pulsar L(t) where the parameter η is is usually called the magnetic fraction of the nebula, i.e. the fraction of the energy injected by the pulsar that powers the magnetic field of the PWN at each moment of time, see e.g., Gelfand, Slane & Zhang (2009) (5) is actually used to define the parameter η, effectively assigning a fraction of the total rotational power available at each time, to particle's energization. However, η is perhaps more properly named as done in Table 1, as the injection sharing parameter, since it describes the instantaneous distribution of the spin-down power into the nebula components. Since particles radiate along the evolution and are subject to losses and escape, the total amount of energy contained by them at any moment is not simply the summation of the power injected at all times. The actual energy partition (how much energy is stored in the magnetic field as opposed to how much energy is stored by the particles) will be relatively enhanced towards the field content in comparison to the case in which such radiation and losses are neglected.
Studies of the nebula evolution have taken into account particles' radiative losses for determining the spectral energy distribution, but have not, to our knowledge, consider them when computing the evolutionary stage. Particles' energy and pressure would be smaller in such cases, for a given value of η. Put otherwise, if we have equipartition at injection at each time (η = 0.5), the nebula would be magnetically dominated in reality; since the particles will loose energy more efficiently than the field. Simple as it sounds, taking into account this effect produces some complexities and requires a special algorithm to deal with them, described below.
Additionally, referring to the formulae above, me and c are the electron mass and the speed of light, respectively.
The minimum energy at injection is a free parameter. If not otherwise stated, we assume it as γmin = 1. The maximum energy at injection is determined by the most restrictive of two different criteria (see e.g., de Jager & Djannati-Ataï 2009). On the one hand, we consider the synchrotron limit where the maximum energy is established by the balance between the energy loss by particles due to synchrotron radiation and the energy gain during acceleration, which yields where e is the electron charge and B the mean magnetic field of the PWN. On the other hand, demanding confinement of the particles inside the termination shock of the PWN where particles are accelerated implies being κ the magnetic compression ratio (we assume it as 3) and ǫ is the containment fraction (ǫ < 1), which is the ratio between the Larmor radius of the pairs and the termination shock radius. The spin-down luminosity of the pulsar evolves in time as being n the braking index of the pulsar and τ0 the initial spin-down age defined as where τc = P/2Ṗ is the characteristic age of the pulsar. The magnetic field evolution is calculated as in Torres, Cillis & Martín Rodriguez (2013). The variation of the magnetic energy in the PWN is balanced by the fraction of the rotational energy that powers the magnetic field ηL(t) and the adiabatic losses due to the expansion of the PWN. This is described by the equation being WB = B 2 R 3 /6 the total magnetic energy. This equation can be solved analytically resulting the following formulae Here again, take into account that in the equations above the instantaneous sharing η is fixed by Eq. (5), and is not varying along the PWN evolution. The radius of the PWN plays an important role in the evolution of the energy losses and the spectrum. In order to take into account the effects of the reverberation of the PWN after interacting with the reverse shock produced by the ejecta, we model as well the evolution of the SNR. Specifically, in order to compute the radius of the PWN during the free expansion phase and compression, we solve the equations given by Chevalier (2005), using a similar prescription as in Gelfand, Slane & Zhang (2009). Specifically, we solve the following system of equations We note that Eq. (14) only applies if vej (R, t) < v(t). Otherwise, dM (t)/dt = 0. In the latter set of equations, v, R, M and P are the velocity, radius, mass and pressure of the PWN shell. vej , ρej and Pej correspond to the values of the velocity, density, and pressure of the SNR ejecta at the position of the PWN shell. γpwn is the adiabatic coefficient of the PWN relativistic gas, which is fixed as 4/3. We consider now the pressure of the PWN. In order to take into account the radiative (as well as other) losses by the particles inside the PWN in the computation of the pressure P (t), we solve Eq. (3) in each time-step and compute the total energy contained in pairs, Ep, by doing the integral, Since Eq.
(3) contains a full accounting of all losses when determining N (γ, t), the equation above is accurately representing the particle's energy share as a function of time.
Then, the pressure contributed by particles is The magnetic field energy also contributes to the total pressure. The magnetic pressure is given by Finally, the total pressure is obtained by summing both contributions As said above, vej , ρej and Pej correspond to the values of the velocity, density, and pressure of the SNR ejecta at the position of the PWN shell. These profiles change if the PWN shell is surrounded by unshocked ejecta (thus the radius of the PWN is smaller than the radius of the reverse shock of the SNR, R < Rrs), or by shocked ejecta (where Rrs < R < Rsnr, being Rsnr the radius of the SNR). The initial profiles for the unshocked medium are assumed, following Truelove & McKee (1999); Blondin, Chevalier & Frierson (2001), as vej (r, t) = r/t for r < Rsnr, 0 for r > Rsnr, where The parameters Esn and Mej are the energy of the SN and the total ejected mass during the explosion, respectively. In our model, for simplicity, we assume ω = 9 as in Chevalier & Fransson (1992 (2009) for a type II SN. When the PWN shell is surrounded by the shocked medium, then we use the prescription given by Bandiera (1984) to obtain the vej , ρej and Pej profiles. For the shock trajectories, we use the semianalytic model by Truelove & McKee (1999) for a non-radiative SNR.
After the compression, the PWN bounces and starts the Sedov phase when its pressure reaches the pressure of the SNR's Sedov solution. At this time, the evolution of the radius of the PWN follows the relation given in where P (t) = ρismv 2 f s /(γsnr + 1) is the pressure in the SNR forward shock. The term v 2 f s is the velocity of the forward shock given by the Truelove & McKee equations.
The integration of Eq. (3), critical in all this modelling, is done following the subsequent prescription: (ii) We proceed to integrate Eqs. (12 -14) with a fourth order Runge Kutta method.
(iii) We compute the new values for the spin-down luminosity and magnetic field using Eqs. (8) and (11), respectively.
(v) We numerically integrate the energy loss, Eq. (3), to get the electron spectrum N (t).
(vi) We obtain the energy and the pressure done by particles by using Eq. (15) and (16). In this way, we take into account the energy losses of the system in order to compute the pressure.
(vii) We obtain the pressure done by the magnetic field by using Eq. (17) and then obtain the total pressure.
(viii) We compute the radius and the velocity of the reverse and forward shock of the SNR using the prescription of Truelove & McKee (1999).
(ix) If the radius of the reverse shock Rrs is greater than the radius of the PWN Rpwn, we approximate the pressure done by the SNR ejecta as 0, and the velocity and the density of the ejecta are described by Eqs. (19) and (20), respectively. If not, we compute these profiles using the prescription given in Bandiera (1984).
(x) We recurse and go back to the first step.  After step viii, when the first compression has already happened, we use Eqs. (12) to compute the radius until we reach the pressure of the Sedov solution for the SNR. At this moment, the evolution of the radius is described by Eq. (24).

SPECTRAL IMPACT OF REVERBERATION
During the reverberation phase, when there are evident effects on the dynamics (compression and bouncing of the nebula radius), there are also important differences in the radiated PWN spectrum when compared to a freely expanding case. To briefly explore this for a PWN having a typical (and very low) instantaneous sharing η, in Table 1 we set the input parameters for a 40 kyr PWN. The evolution of the radius and the magnetic field for the reverberation and the free expansion models of this nebula is shown in Fig. 2. We can observe that both quantities evolve in the same way during the first ∼4 kyr (see in Fig. 2). Also the spectrum is similar. The age at which the transition between different stages of the evolution occurs varies with the energy of the SN explosion, the ejected mass, and the initial profiles of the SNR ejecta. At an age of ∼4 kyr, the PWN shell goes into the shocked medium of the remnant and starts the compression. During this phase, the magnetic field and the internal pressure increases as they evolve approximately as B(t) ∼ R −3/2 , P ∼ R −3 . In order to see this first deviation, in Fig. 3 we show the comparison of the pairs and the PWN spectrum at the age of 6 kyr. The increase of the magnetic field enhance the synchrotron flux and cools down the high energy pairs (γ ∼ 10 8 − 10 10 ) while increasing slightly the population at γ ∼ 10 5 − 10 8 . The pairs at this range of energy interact with the IC target photon fields (CMB, FIR and NIR), thus also increasing slightly the IC radiation flux between 1 GeV and 10 TeV. These effects become more evident as the compression goes on, as we can see at 10 kyr (second row of panels in Fig. 3).
At the moment of maximum compression, the magnetic field is high enough to burn off the majority of high energy pairs and move them to lower energies in the pair population spectrum. When the pressure is high enough too, the PWN bounces and starts re-expanding (∼16 kyr). Then, the magnetic field in the nebula decreases again and the new high energy particles that are injected by the pulsar start accumulating. In Fig. 3, we observe these effects at 18 kyr, where we see how the old high energy pairs are accumulated at the lower energy part of the pair spectrum and the quantity of high energy particles is orders of magnitude lower in comparison to the free expansion model. The PWN spectrum shows that the IC flux is two orders of magnitude lower after compression. The synchrotron flux is still higher for the reverberation case, where the radius of the PWN is smaller and the magnetic field is high enough to produce a synchrotron flux two orders of magnitude higher than in the free expansion model. Differences between models considering reverberation and not, are large, as expected (see Gelfand, Slane & Zhang (2009)).
After the bounce (see the 40 kyr case in Fig. 3), the radius of the PWN keeps growing freely until reaching the pressure of the Sedov solution. At this stage, the PWN ex-  pansion is in equilibrium with the ambient medium and follows the Taylor-Sedov law. The pulsar injects newly accelerated pairs that contributes to reestablish the high energy pair population. However, we still note differences due to the extinction of high energy particles during the compression phase and the synchrotron losses, which are higher than in the free expansion model, where the magnetic field is very low due to the unrealistic size of the PWN at this age if a free-expansion model is considered. Note that at the same age after the bounce, considering the same instantaneous sharing, the synchrotron flux in the free expansion model may be up to 3 to 4 (depending on the frequency) orders of magnitude lower than what we get if we include reverberation for the case we studied.

An evolving energy partition
Since the equation for the pressure evolution of the particle content in the nebula considers the particles's adiabatic and radiative losses, while the magnetic energy follows the adiabatic law, the ratio of particle's energy to magnetic field energy, is no longer constant in time, neither it is equal to 1 − η. However, the modelling provides B(t), as well as the particles's energy via Eq. (15), and thus the energy ratio can be computed directly. This is shown in Fig. 4, where we plot the ratio of EB with the total energy EB + Ep = E total and with the particle's energy (the averaged σ parameter). Due to the high-energy losses at the beginning of the evolution, the ratio between the magnetic energy and the total energy content of the nebula is high with respect the parameter η (the free parameter for the instantaneous sharing, as explained in the text above) used during injection. During the free expansion evolution of the nebula, the energy partition approaches to η, but it is still larger than it (what can be understood due to the fact that particles loose energy via radiation, and the field does not). In compression, then, the increase of the magnetic field makes the energy partition to significantly deviate from the free expansion case.

MODELLING OF THE PWN IN CTA 1
We shall first consider the constraints that the SNR and PWN radius impose on models. The CTA 1 radio SNR has a radius of ∼ 20.4d1.4 pc (Pineault et al. 1993); therefore, taking into account the uncertainty on the distance, we will accept models able to cope with a constraint in the SNR radius of 19.5 < Rsnr < 21.5 pc. On the other hand, our model (being 1D) assumes a spherical system. We take as a reference the dimensions of CTA 1 PWN described by VERITAS (7.3d1.4 × 5.9d1.4 pc, Aliu et al. 2013), and we will consider acceptable any value for the PWN radius within this range.
The radius of the SNR and the PWN are determined by three parameters: the ejected mass Mej , the age the system tage, and the interstellar medium (ISM) density ρism. The rest of the assumed parameters (including those taken directly from observations) related with the evolution are We show them separated in three different plots to avoid superposition of the points. The appearance of models with the same ejected mass and ISM density but with different ages is due to the uncertainty on the CTA 1 distance and the tolerance that we allow to fit the radius of the PWN. The right plot on the bottom correspond to models in compression phase.
quoted in Table 2. We consistently computed the initial spindown luminosity and the spin-down age. We explore models that fit the radii assuming a range of ejected masses between 9 and 13M⊙ (in the usual range for a SN explosion with a typically massive progenitor, note instead that Aliu et al. (2013) have assumed a lower ejected mass of 6.5 M⊙) and ages that range from 9 kyr to 13 kyr (Slane et al. 2004). We have computed ∼650 spectra varying the assumed average density from 0.05 cm −3 to 0.1 cm −3 . The latter value is the upper limit for the ISM density of the Galaxy at the CTA 1 position (Dickey & Lockman 1990;Slane et al. 1997). We have considered too the observational indications for a density of 0.017 Slane et al. (2004) and 0.037 cm −3 Slane et al. (1997) in particular regions. In any case, all these values are still much lower than the upper limits we have found above from our analysis of molecular material in the vicinity.
Depending on the values of Mej , tage, and ρism, we have different compatible models fitting the constraints for the radii, given the observational uncertainties. Their spectral features are different, since the PWN would be in a different stage of evolution (free expansion or compression), and we discuss them next. Our first result is that no models in the late Sedov expansion phase, within the explored parameter ranges mentioned above, fit the observed radii.
In Fig. 5, we show those models that do fit the observational values for the radii of the SNR and the PWN. Note that some models have, e.g., similar ejected mass and ISM density, but the age of the system is different. In order to see the variety of these models, Fig. 5 shows three plots for those models which are in free expansion (see the caption for a detailed explanation). The models in free expansion fit the spectrum with instantaneous sharings between 0.25 and 0.4. An upper limit for the FIR and NIR densities is established at 0.037 and 0.075 eV cm −3 , respectively. These values are lower than those obtained from GALPROP 4 (Porter, Moskalenko & Strong 2006), which predicts w f ir = 0.3 eV cm −3 and wnir = 0.6 eV cm −3 for the Galactic position of CTA 1, and can be understood in the context of local fluctuations.
For models in compression and those cases in which the constraints in the radii are fulfilled, as shown in Fig. 5, we do not get an acceptable fit for the TeV flux. In all these cases, the IC flux due to the CMB target field would already exceed the observational constraints. This can be understood by the accumulation of cooled electrons with Lorentz factors between γ ∼ 10 5 and 10 8 , as shown in Fig. 3. Fig. 6 shows a free expansion model for CTA 1 that fits the spectrum, and a compression stage model that does not, as an example. The parameters for these models are summarized in Table 3. The fitted parameters are mainly those related with the particle injection, photon environment, and instantaneous sharing fraction. Some other parameters in the table are derived, using model equations (see Table 2). Note that the free expansion model considers an ISM density of 0.037 cm −3 . At ρISM = 0.017 cm −3 , all models have a SNR radius larger than 21.5 pc. The free expansion model fitting that we show comes from the refinement of the parameters in order to fit an "average" radius for CTA 1 PWN of 6.7 pc. Taking this into account, we get only one model Energy (GeV) Figure 6. Two models for CTA 1. In free expansion (left), a high instantaneous sharing parameter of ∼0.25 is needed to achieve the X-ray fluxes and the energy densities for the FIR and NIR target fields need to be low (see Table 3). In compression ( Data points in radio come from Giacani et al. (2014), in X-rays from Lin et al. (2012), in VHE from Aliu et al. (2013). The Fermi upper limits are given in Acero et al. (2013).
as viable, from the two shown in Fig. 5 for ρISM = 0.037 cm −3 , with an ejected mass of 12M⊙ and an age of 9.2 kyr. Regarding the evolution of the energy partition in the CTA nebula, we show it in Fig. 7 for the two models considered above, as an example.

Discussion
We see that the observational data for the CTA 1 PWN is problematic, and cannot be satisfactorily fitted by our relatively-detailed model without considering important caveats.
If we consider a fit during the free expansion phase, the instantaneous sharing (and the energy partition of the nebula) needs to be relatively high (in comparison to all other young PWNe, see Torres et al. 2014). This is needed in order to obtain the appropriate magnetic field to fit the X-ray data with synchrotron radiation. If the PWN is in a compression state, the problem arise from over-predicting the high energy fluxes, as shown above.
We recall that additional considerations done in our simulations are the fixing of the braking index as n = 3, the energy break as γ b = 10 6 , and the low energy index as α l = 1.5. The latter two parameters are unconstrained, since at low energies in the spectrum, we only have an upper limit on the radio flux from the estimation done by Giacani et al. (2014). The values we have assumed for these parameters are typical (Bucciantini, Arons & Amato 2011;Torres et al. 2014). However, in Fig. 8 we show the results of variations in these parameters taking as a reference model the fit for the PWN in free expansion that we show in Table 3. In the upper panels of Fig. 8, we show the variation of the spectrum with the braking index going from 2.1 to 3. Note that the spectrum is almost independent on this parameter. In the middle panels, we show the variation with the energy break from γ b = 10 5 to 10 7 . Note that models with γ b > 10 6 exceeds the observational synchrotron and IC flux. The synchrotron flux can be corrected by decreasing the in-stantaneous sharing (which plays in favour to the tendency that we observe in Torres et al. 2014), but the IC flux generated considering only the CMB contribution is too high. Note also that those models with γ b < 10 6 could fit better the IC flux using the FIR and NIR energy densities given by GALPROP, but the synchrotron X-ray flux is too low and we would need to increase the instantaneous sharing even more, up to the point to exceed the upper limit in the radio band. Regarding the lower panel, we see a similar behaviour when we change the value of the low energy index α l . We vary this parameter from 1 to 2. We have also varied the slope of the density describing the ejecta, from ω=7 to 13, finding that the differences with our reported results are relatively minor.
From a modelling perspective, the age of the PWN as well as its morphology does not favor it being in a pure free expansion phase, despite the important caveats that we have seen for the compression phase in our models. A way out could be that the cooling of pairs during compression is different from what we have been able to model, due perhaps to a radial dependence of the magnetic field (e.g., Gaensler & Slane 2006). If so happens, particles in different regions of the PWN would not be affected by synchrotron energy losses in the same way, changing the resulting pair population. An asymmetric evolution of the PWN due to e.g., different ISM densities, could explain this scenario. In this situation, one can imagine that some of the pairs would effectively evolve as in a freely expanding nebula, and some others would be affected by the compression of the SNR reverse shock. An intermediate situation where the overall IC yield is not increased beyond the measurements while an increased magnetic field reduce the need of a larger magnetization can be imagined. However, this exploration is beyond the ability of our numerical scheme and would probably require a morphological/radiative approach to be consistently solved at once. Finally, in terms of our model, the nebula cannot be re-expanding (Sedov phase) after compression.  Figure 7. Ratio between the magnetic energy content to the total energy content (particle + field) for the CTA 1 PWN model. We include the σ-parameter in a dashed line (hardly visible in the right plot). Left panel: free expansion model with an age of 9.2 kyr; right panel: (compression model at an age of 11.4 kyr; the change in these ratios due to compression (as shown above in Fig. 4 is still not visible, but would show up similarly, and well before 20 kyr, if we let this model continue run in time. For more details see text. Energy (GeV) Figure 8. Variation of the electron (left) and spectral (right) model for CTA 1 with the braking index n (upper panels), energy break γ b (middle panels) and low energy spectral index α l (lower panels). The color legend is: for upper panels, the red, orange, green, blue, purple and black lines correspond to n = 2.1, 2.3, 2.5, 2.7, 2.9, 3, respectively; for middle panels, the red, orange, yellow, green, cyan, blue, purple, grey and black lines correspond to γ b = 10 5 , 3 × 10 5 , 5 × 10 5 , 7 × 10 5 , 10 6 , 3 × 10 6 , 5 × 10 6 , 7 × 10 6 , 10 7 ; and for lower panels, the red, orange, green, blue, purple and black lines correspond to α k = 1, 1.2, 1.4, 1.6, 1.8, 2.

CONCLUDING REMARKS
We have found that the TeV emission from the direction of CTA 1 is very unlikely related to hadronic illumination of a molecular cloud in the vicinity of the SNR. Planck and CO Harvard survey data are inconsistent with this interpretation. Thus, understanding the multifrequency emission with a PWN model is important.
We have included the prescription given by Chevalier (2005); Gelfand, Slane & Zhang (2009) to solve the Euler equations of the PWN during the reverberation phase into the radiative model described in Torres et al. (2014), and use it to study the effects of the compression and re-expansion of a PWN due to the interaction with the shocked medium of the SNR. We have included the energy radiated by particles and their escape when computing the pressure produce by the gas. This leads to an evolving energy partition, since for the same instantaneous sharing (magnetic fraction) of the injection of energy provided by the rotational power, the field and the particles are affected differently by radiation and losses.
We observed that there are significant effects in the spectrum compared with a free expansion model that induce differences in flux predictions at all energies of one or two orders of magnitude. We note that with the same instantaneous sharing, models that take into account reverberation have more synchrotron flux during the bounce than models considering the nebula to be in free expansion at the same age. This may severely affect the determination of the instantaneous sharing η, and ultimately, of the energy partition.
We have applied this more complete radiative/dynamical model to CTA 1, which is a critical case due to its age (in the high end of those considered to be young PWNe) and atypical because of its requirement of a relatively larger magnetization when free expansion models are considered to fit the spectrum. We have made a parameter space exploration for the ejected mass, age, and ISM density that would give SNR and PWN radii as measured, and consider those in order to fit the CTA1 PWN spectrum. We could not find any completely acceptable representation at low values of magnetization, in any phase. Free expansion models (the only ones that can actually fit the data in the explored phase space) need instantaneous sharing parameters that are more than one order of magnitude larger than the ones typically obtained for young PWNe. Compression models exceed the IC flux data even only considering the CMB target photon field. And there are no models in a late Sedov phase after reverberation that is compatible with the radii of the CTA 1 PWN and SNR inside the parameter space that we have considered. A better low and high-energy determination of the SED, in radio and GeV, in particular, would help distinguishing models better. CTA 1 thus remains as a difficult case, and a testbed for new models incorporating more complicated dynamics.