On the population of remnant FRII radio galaxies and implications for radio source dynamics

TL;DR

This study quantifies the rarity of ultra-steep spectrum remnant FRII radio galaxies in a flux-limited sample and uses simulations to explore their evolution, suggesting most remnants are young and over-pressured, with implications for future observations.

Contribution

It provides the first quantitative estimate of ultra-steep spectrum remnant FRII radio galaxies and tests models of their evolution through Monte-Carlo simulations.

Findings

Less than 2% of FRII galaxies are ultra-steep spectrum remnants.

Models with Sedov-like expansion fit the observed upper limit.

Most remnants are predicted to be young and over-pressured.

Abstract

The purpose of this work is two-fold: (1) to quantify the occurrence of ultra-steep spectrum remnant FRII radio galaxies in a 74 MHz flux limited sample, and (2) perform Monte-Carlo simulations of the population of active and remnant FRII radio galaxies to confront models of remnant lobe evolution, and provide guidance for further investigation of remnant radio galaxies. We find that fewer than 2$\%$ of FRII radio galaxies with S$_{ \rm74~MHz} > 1.5$ Jy are candidate ultra-steep spectrum remnants, where we define ultra-steep spectrum as $\alpha_{\rm 74~MHz}^{\rm 1400~MHz} > 1.2$. Our Monte-Carlo simulations demonstrate that models involving Sedov-like expansion in the remnant phase, resulting in rapid adiabatic energy losses, are consistent with this upper limit, and predict the existence of nearly twice as many remnants with normal (not ultra-steep) spectra in the observed frequency range as there are ultra-steep spectrum remnants. This model also predicts an ultra-steep remnant fraction approaching 10$\%$ at redshifts $z < 0.5$. Importantly, this model implies the lobes remain over-pressured with respect to the ambient medium well after their active lifetime, in contrast with existing observational evidence that many FRII radio galaxy lobes reach pressure equilibrium with the external medium whilst still in the active phase. The predicted age distribution of remnants is a steeply decreasing function of age. In other words young remnants are expected to be much more common than old remnants in flux limited samples. For this reason, incorporating higher frequency data $\gtrsim 5$ GHz will be of great benefit to future studies of the remnant population.

Tables

Table 1. Table 1: Results of cross-matching our VLSSr selected sample with NVSS.

Number of NVSS matches	Number of VLSSr sources
0	4
1	3466
2	382
3	9
$>3$	0
Total	3861

Table 2. Table 2: Results of cross-matching our VLSSr sample with the FRII candidate sample of van Velzen et al. ( 2015 ) .

	Entire VLSSr selected Sample	Ultra-steep Spectrum Sample
		(${\alpha }_{\mathrm{VLSSr}}^{\mathrm{NVSS}}>1.2$)
Sample Size	3861	57
Fraction matched with FRII candidates	65$\%$ (2498/3861)	45$\%$ (26/57)
Fraction of FRII candidates with detected core	14$\%$ (348/2498)	23$\%$ (6/26)

Table 3. Table 3: Results of Monte Carlo Simulations: All redshifts.

	Maximal (Sedov-like) expansion			No expansion
	in remnant phase			in remnant phase
Quantity	Number of sources	Fraction of sample		Number of sources	Fraction of sample
Active sources	9499	94.6$\%$		5372	82$\%$
All Remnant sources	539	5.4$\%$		1199	18$\%$
Low $\nu$ ultra-steep spectrum remnants (${\alpha }_{74\mathrm{MHz}}^{1400MHz}>1.2$)	187	1.8$\%$		751	11.4$\%$
Mid $\nu$ ultra-steep spectrum remnants (${\alpha }_{1400MHz}^{5GHz}>1.2$)	378	3.8$\%$		1036	15.8$\%$
High $\nu$ ultra-steep spectrum remnants (${\alpha }_{5\mathrm{GHz}}^{10GHz}>1.2$)	441	4.4$\%$		1110	16.9$\%$
Curved spectrum remnants (${\alpha }_{1400\mathrm{MHz}}^{5000\mathrm{MHz}}-{\alpha }_{74\mathrm{MHz}}^{1400\mathrm{MHz}}>0.5$)	292	2.9$\%$		976	14.9$\%$
Curved spectrum remnants (${\alpha }_{5000\mathrm{MHz}}^{10000\mathrm{MHz}}-{\alpha }_{74\mathrm{MHz}}^{1400\mathrm{MHz}}>0.5$)	385	3.8$\%$		1068	16.3$\%$

Table 4. Table 4: Results of Monte Carlo Simulations: Low redshifts only ( z < 0.5 𝑧 0.5 z<0.5 ).

	Maximal (Sedov-like) expansion			No expansion
	in remnant phase			in remnant phase
Quantity	Number of sources	Fraction of sample		Number of sources	Fraction of sample
Active sources	496	80$\%$		292	37$\%$
All Remnant sources	123	20$\%$		505	63$\%$
Low $\nu$ ultra-steep spectrum remnants (${\alpha }_{74\mathrm{MHz}}^{1400MHz}>1.2$)	20	3.2$\%$		265	33$\%$
Mid $\nu$ ultra-steep spectrum remnants (${\alpha }_{1400\mathrm{MHz}}^{5GHz}>1.2$)	57	9.2$\%$		419	53$\%$
High $\nu$ ultra-steep spectrum remnants (${\alpha }_{5\mathrm{GHz}}^{10GHz}>1.2$)	79	12.8$\%$		456	57$\%$
Curved spectrum remnants (${\alpha }_{1400\mathrm{MHz}}^{5000\mathrm{MHz}}-{\alpha }_{74\mathrm{MHz}}^{1400\mathrm{MHz}}>0.5$)	47	7.6$\%$		397	50$\%$
Curved spectrum remnants (${\alpha }_{5000\mathrm{MHz}}^{10000\mathrm{MHz}}-{\alpha }_{74\mathrm{MHz}}^{1400\mathrm{MHz}}>0.5$)	69	11.1$\%$		448	56$\%$

Equations

$$S_{\nu }=\frac{(1+z)}{D_L^2}\frac{\sqrt{3}{e}^3}{16{\pi }^2{ϵ}_{0}c{m}_e}{\int }_{{\gamma }_{\min }}^{{\gamma }_{\max }}N(\gamma )\left[{\int }_0^{\infty }B{p}_B\bar{F}(y)dB\right]d\gamma$$

$$p_B=\sqrt{\frac{2}{\pi }}\frac{B^2\exp (-B^2/2a^2)}{a^3}$$

$$a=\frac{B_0}{\sqrt{3}}$$

$$\int B^2{p}_{B}dB=B_0^2$$

$$\rho \propto r^{-\beta }$$

$$V(t)\propto \{\begin{array}{ll}{t}^{9/(5-\beta )},& \text{if }t<t_{\mathrm{on}}.\\ t^{6/(6-\beta )},& \text{if }t>t_{\mathrm{on}}.\end{array}$$

$$B(t)=\{\begin{array}{ll}{B}_0{\left(\frac{t}{t_0}\right)}^{\frac{-4-\beta }{2(5-\beta )}}{\left(\frac{Q}{Q_0}\right)}^{\frac{2-\beta }{2(5-\beta )}},& \text{if }t<t_{\mathrm{on}}.\\ B(t_{\mathrm{on}}){\left(\frac{t}{t_{\mathrm{on}}}\right)}^{\frac{-4}{(6-\beta )}},& \text{if }t>t_{\mathrm{on}}.\end{array}$$

$$\frac{dN}{d{\gamma }_{i}d{t}_i}=\{\begin{array}{ll}{q}_0{\gamma }_i^{-a},& \text{if }t<t_{\mathrm{on}}\text{ and }{\gamma }_{\mathrm{i},\min }<{\gamma }_i<{\gamma }_{\mathrm{i},\max }.\\ 0,& \text{otherwise }\end{array}$$

$$V(t)\propto \{\begin{array}{ll}{t}^{9/(5-\beta )},& \text{if }t<t_{\mathrm{on}}.\\ \text{mbox}constant,& \text{if }t>t_{\mathrm{on}}.\end{array}$$

$$B(t)=\{\begin{array}{ll}{B}_0{\left(\frac{t}{t_0}\right)}^{\frac{-4-\beta }{2(5-\beta )}}{\left(\frac{Q}{Q_0}\right)}^{\frac{2-\beta }{2(5-\beta )}},& \text{if }t<t_{\mathrm{on}}.\\ \text{mbox}constant,& \text{if }t>t_{\mathrm{on}}.\end{array}$$

$$\frac{dN}{d\gamma }\left(\gamma ,t\right)=q_0{\gamma }^{-a}{\int }_{t_{\mathrm{i},\min }}^{t_{\mathrm{i},\max }}{\left(\frac{V(t_i)}{V(t)}\right)}^{(a-1)/3}{\left(1-\frac{\gamma }{{\gamma }_{\ast }(t_i,t)}\right)}^{a-2}d{t}_i$$

$$\frac{1}{{\gamma }_{\ast }(t_i,t)}={\int }_{t_i}^t{a}_0{\left(\frac{V(\tau )}{V(t)}\right)}^{-1/3}\left(\frac{B^2(\tau )+B_{\mathrm{CMB}}^2}{2{\mu }_0}\right)d\tau ,$$

$$t_{\mathrm{i},\min }=\text{mbox}MAX(0,t_{\mathrm{i},\min }^{\ast }).$$

$$\frac{1}{\gamma }=\frac{1}{{\gamma }_{\mathrm{i},\max }{\left(\frac{V(t_{\mathrm{i},\min }^{\ast })}{V(t)}\right)}^{1/3}}+\frac{1}{{\gamma }_{\ast }(t_{\mathrm{i},\min }^{\ast },t)}$$

$$p(z)\propto \rho (z)\frac{dV}{dz}$$

$$\rho (z)\propto \exp \left[-\frac{1}{2}\left(\frac{z-z_{h0}}{z_{h1}}\right)\right]$$

$$\frac{dV}{dz}\propto \frac{(1+z{)}^2{D}_A^2}{\sqrt{{\Omega }_M(1+z{)}^3+{\Omega }_{\Lambda }}}$$

$$p(Q_{\mathrm{jet}})\propto \{\begin{array}{l}{Q}^{-n_Q}\text{ if }5\times 10^{36}W<Q_{\mathrm{jet}}<2\times 10^{42}W\\ =0\text{ otherwise}\end{array}$$

$$\frac{dN}{d\gamma }\left(\gamma ,t\right)=q_0{\gamma }^{-a}{\int }_{t_{\mathrm{i},\min }}^{t_{\mathrm{i},\max }}{\left(\frac{V(t_i)}{V(t)}\right)}^{(a-1)/3}{\left(1-\frac{\gamma }{{\gamma }_{\ast }(t_i,t)}\right)}^{a-2}d{t}_i$$

$$\frac{dN}{d\gamma }\left(\gamma ,t\right)<q_0{\gamma }^{-a}{\int }_{t_{\mathrm{i},\min }}^{t_{\mathrm{i},\max }}{\left(\frac{V(t_i)}{V(t)}\right)}^{(a-1)/3}d{t}_i$$

$$L_{\nu }$$

$$\frac{dN}{dtd{L}_{\mathrm{peak}}}$$

$$\frac{dN}{dtd{L}_{\nu }}\propto t^{a_2(1-p)}$$

$$V(t)\propto \{\begin{array}{ll}{t}^{9/(5-\beta )},& \text{if }t<t_{\mathrm{on}}.\\ V(t_{\mathrm{on}}),& \text{if }t>t_{\mathrm{on}}.\end{array}$$

$$B(t)=\{\begin{array}{ll}{B}_0{\left(\frac{t}{t_0}\right)}^{\frac{-4-\beta }{2(5-\beta )}}{\left(\frac{Q}{Q_0}\right)}^{\frac{2-\beta }{2(5-\beta )}},& \text{if }t<t_{\mathrm{on}}.\\ B(t_{\mathrm{on}}),& \text{if }t>t_{\mathrm{on}}.\end{array}$$

$$j_{\nu }={\int }_0^{\infty }{\int }_0^{\pi }{\int }_{{\gamma }_{\min }}^{{\gamma }_{\max }}\frac{\sqrt{3}B{e}^3\sin \Theta }{8{\pi }^2{ϵ}_{0}c{m}_e}F(x)n(\gamma )p_{\Theta }{p}_{B}d\gamma d\Theta dB$$

$$j_{\nu }=\frac{\sqrt{3}{e}^3}{8{\pi }^2{ϵ}_{0}c{m}_e}{\int }_{{\gamma }_{\min }}^{{\gamma }_{\max }}n(\gamma )\left[{\int }_0^{\infty }B{p}_B\left[{\int }_0^{\pi }\sin \Theta F(x)p_{\Theta }d\Theta \right]dB\right]d\gamma$$

$$y=\frac{\nu \sin \Theta }{{\nu }_c}=\frac{4\pi \nu }{3{\Omega }_0{\gamma }^2}.$$

$$\bar{F}(y)={\int }_0^{\pi }\sin \Theta F(x)p_{\Theta }d\Theta .$$

$$j_{\nu }=\frac{\sqrt{3}{e}^3}{8{\pi }^2{ϵ}_{0}c{m}_e}{\int }_{{\gamma }_{\min }}^{{\gamma }_{\max }}n(\gamma )\left[{\int }_0^{\infty }B{p}_B\bar{F}(y)dB\right]d\gamma$$

Full text

††pubyear: 2017††pagerange: On the population of remnant FRII radio galaxies and implications for radio source dynamics.–C

1 Introduction

Do Fanaroff-Riley type II (FRII; Fanaroff & Riley, 1974) radio galaxies remain over-pressured with respect to the external medium throughout their entire active life, or do they at some point reach pressure equilibrium with their environment? The answer to this question has potentially far-reaching implications for the understanding of radio galaxy dynamics and energetics. If the lobes remain over-pressured throughout the entire life, the source is always surrounded by an elliptical bow shock and the lobes undergo supersonic self-similar expansion (Kaiser & Alexander, 1997). This self-similar scenario is typically assumed when modelling populations of FRII radio galaxies (for example Blundell et al., 1999; Wang & Kaiser, 2008; Kapińska et al., 2012, to name a few). However, some studies now point to a more complex situation, in which the lobes start out being highly over-pressured, but reach pressure equilibrium with their environment well before the jet activity comes to an end. Estimates of internal lobe pressures based on the assumption of minimum energy (Hardcastle & Worrall, 2000) as well as those based on inverse-Compton measurements (Hardcastle et al., 2002; Croston et al., 2004) are observed to be comparable to the external pressures, assuming no large contribution to the lobe pressure from thermal material within the lobes. Furthermore, Mullin et al. (2008) find a trend of increasing axial ratio with sources size in a flux limited sample of FRIIs with $z<1$. If the lobes remained over-pressured throughout their lifetime, undergoing self-similar expansion, then the axial ratio is expected to remain constant. The observed size-dependent axial ratio distribution demonstrated by Mullin et al. (2008) is therefore inconsistent with models in which the expansion remains self-similar. Hardcastle & Worrall (2000) and Mullin et al. (2008) suggest that FRII radio galaxies may only grow self-similarly early on in their lifetime, and reach pressure equilibrium with their surroundings in middle-age.

One way to discriminate between these two scenarios is to consider the remnant phase of radio galaxy evolution. If the lobes remain highly over pressured throughout their lifetime, then the remnant phase will be governed by supersonic Sedov-like expansion (Kaiser & Cotter, 2002). This continued supersonic lobe expansion in the remnant phase will cause rapid dimming of the lobe radio emission due to a decrease in the magnetic field strength and decrease in particle energies due to adiabatic expansion losses. In contrast, if the lobes are already in pressure equilibrium with the external medium at the end of the active phase, the luminosity evolution in the remnant phase is expected to be much more sedate.

Kaiser & Cotter (2002) showed that models of remnant radio galaxy spectra are highly degenerate, and modelling of individual remnant radio galaxy spectra cannot constrain the history of lobe evolution. The only way to constrain the lobe evolution in the remnant phase is via a statistical approach. That is the approach taken in this paper. We compose a flux-limited sample radio sources that is dominated by FRII radio galaxies (Section 2). With this sample, we obtain an upper limit on the number of ultra-steep spectrum remnant radio galaxies in a low frequency selected, flux limited sample of FRII radio galaxies (Section 3). We then perform Monte-Carlo simulations to assess whether models of remnant phase lobe evolution are consistent with the observed limit on ultra-steep spectrum remnants in our flux limited sample (Sections 4 and 5).

We make a strong distinction between the population of ultra-steep spectrum remnants, which we define as remnant radio galaxies with spectral index $\alpha>1.2$ (defined such that $S_{\nu}\propto\nu^{-\alpha}$) between our chosen frequencies, and the entire population of remnant radio galaxies, which includes any remnant radio galaxy regardless of spectral characteristics. This distinction is necessary because not all remnant radio galaxies have ultra-steep spectra in the observed frequency range. A good example of this is the remnant radio galaxy discovered by Brienza et al. (2016), which was identified purely based on morphological characteristics, and only shows an ultra-steep spectrum above 1.4 GHz.

In this work, we restrict our analysis to FRII radio galaxies. A follow-up paper (Brienza et al., 2017 submitted) will focus on the study of lower luminosity, FRI class remnant radio galaxies.

2 Sample selection

2.1 Sample Definition

Our sample selection is based on the 74 MHz VLA Low-freqeuncy Sky Survey Redux catalog (VLSSr; Lane et al., 2014). We calculate the flux density of sources in the VLSSr catalogue from the catalogued peak intensity, and fitted major and minor axes, based on the expressions given in Condon (1997); Cohen et al. (2007). We then restrict the VLSSr catalog as follows:

9 hrs $<$ RA $<$ 16 hrs 2. 2.

$0^{\circ}<$ DEC $<60^{\circ}$ 3. 3.

Distance to nearest neighbour D $>$ 4 arcminutes. 4. 4.

Fitted major axis $<$ 120 arcseconds. 5. 5.

Flux density $S_{\rm 74~{}MHz}>1.5$ Jy

The reason for each of the restrictions is as follows:

(i and ii) The restrictions on RA and DEC are imposed in order to match the sky area of the FIRST survey, so that we can assess the morphology of the selected objects.

(iii) The median distance to the nearest neighbour in the VLSSr catalogue is approximately 15 arcminutes. However, a histogram of distance to nearest neighbour shows a clear peak in the range 1 - 4 arcminutes. The narrow peak at 1-4 arcminutes is a result of radio galaxies with large angular size and complex morphology, in which the individual radio galaxies are fitted by more than one Gaussian component. To simplify the cross-matching with NVSS, we restrict our sample to “isolated" VLSSr sources, for which the distance to nearest neighbour is greater than 4 arcminutes. This selection criterion reduces the catalogue size by 8$\%$. However, for sources with nearest neighbour $<$ 4 arcminutes, several catalogued sources are often related to the same radio galaxy, and so the reduction in the number of radio galaxies is likely to be significantly less than 8$\%$.

(iv) We next remove all objects for which the fitted major axis size is equal to 120 arcseconds, which is the upper limit on fitted major axis (Lane et al., 2014). For these objects, we are unable to accurately calculate the flux density based on the catalogued values of peak intensity, major and minor axes. This selection criterion reduces the catalogue size by a further 8$\%$ (1237 objects). These very large angular size sources would provide an interesting sample for searches for remnant radio galaxies, as they most likely represent very large, relatively nearby radio galaxies. However, the need to obtain accurate flux density estimates from the catalogued values of peak intensity, major axis and minor axis, necessitates the removal of these large sources from our analysis.

(v) Finally, we impose a flux density limit of $S_{\rm 74~{}MHz}>1.5$ Jy, which corresponds to the knee in the source-counts distribution (eg. Massardi et al., 2010). With this flux limit, the sample is expected to be highly dominated by high excitation (predominantly FRII) radio galaxies, with an FRII fraction of up to 80$\%$ (Willott et al., 2001). We note that the S-cubed simulation of Wilman et al. (2008) predicts that our sample will comprise 55$\%$ FRII, 5$\%$ Gigahertz Peaked Spectrum (GPS), and 40$\%$ FRI111The lowest frequency in the S-cubed database is 151 MHz, so we converted our 74 MHz flux limit to 151 MHz assuming a spectral index of 0.75.. The morphologies of our sample agree with the prediction of the S-cubed simulation: $\gtrsim 50\%$ of our sample have a double-lobed FRII-like appearance. Whilst the sample contains only $\gtrsim 50\%$ FRII radio galaxies, we account for this factor in our analysis of the FRII remnant population, as described in section 3.1.

The flux limit further reduces the sample size by 72$\%$, giving a final sample size of 3861 objects.

2.2 Cross-matching with NVSS

We cross-match our 74MHz flux limited sample described in Section 2.1 with the NVSS catalogue at 1.4 GHz (Condon et al., 1998). Due to the higher frequency and higher resolution of the NVSS relative to the VLSSr, individual catalogue entries in the VLSSr may be associated with multiple catalogue entries in the NVSS. For this reason, we must use a large matching radius. Given that the largest fitted major axis in the VLSSr sample is 120 arcseconds, we use a matching radius of 60 arcseconds, and sum the flux densities of all NVSS catalogue entries within the matching radius. The results of our cross-matching is summarised in Table 1.

The probability of finding an NVSS source within 60 arcseconds of a random position is approximately 4$\%$ (Condon et al., 1998). We therefore expect that approximately 4$\%$ of our sample ($\sim 150$ objects) will contain an unrelated NVSS source within the search radius of 60 arcseconds. However, these unrelated sources are clustered near the NVSS flux limit (2.5 mJy), while the majority of our sample have NVSS flux densities more than 10 times higher than the flux limit (see Figure 1). Therefore, only a very small fraction of our sample will be significantly affected by the detection of an unrelated NVSS object within the search radius.

The differing survey resolution is not likely to cause significant errors in spectral index calculation. Our sample is restricted to sources with fitted major axes less than 120 arcseconds in the VLSSr catalogue, a factor of only 2.7 times the NVSS angular resolution of 45 arcseconds. Furthermore, our sample has a high signal to noise in the NVSS survey, and therefore the catalogued flux densities are unlikely to have missed flux.

3 Empirical Results

In Figure 1 we present the histogram of spectral index between 74 MHz and 1400 MHz for our sample, as well as the histogram of 74 MHz and 1400 MHz flux densities. The spectral index distribution shows a broad, symmetric central peak, along with a flat spectrum tail and an ultra-steep spectrum tail. We wish to place an upper limit on the fraction of our sample that are ultra-steep spectrum remnant radio galaxies. To do so, we split our sample into 3 spectral categories: flat ($\alpha<0.45$); normal ($0.45<\alpha<1.2$); and ultra-steep ($\alpha>1.2$). The choice of the dividing line between flat, normal and ultra-steep spectrum is somewhat arbitrary, but is guided by the shape of the histogram in Figure 1, and models of radio source emission. The ultra-steep spectrum limit at $\alpha=1.2$ corresponds to the maximum spectral index for active radio galaxy models with particle injection corresponding to $\alpha<0.7$ and a cooling break of $\Delta\alpha=0.5$. We stress that our choice of dividing line at $\alpha=1.2$ is not intended to capture all remnant radio galaxies. Indeed, we expect that most remnant radio galaxies above our flux limit will not have ultra-steep spectra according to our definition. However, the key results are not strongly affected by our choice of dividing line, since our modelling approach is designed specifically to accommodate arbitrary selection criteria.

Flat spectrum objects comprise $2\%$ of the sample (78 objects), “normal" spectrum objects comprise 96.5$\%$ of our sample (3726 objects), and ultra-steep spectrum objects comprise only $1.5\%$ of the sample (57 objects). Note, however, that our sample is comprised of both FRI and FRII radio galaxies. In this work, we seek an upper limit on the fraction of FRII radio galaxies that are ultra-steep spectrum remnants. We therefore need to account for the fraction of the sample that is FRII. We consider this question in the following sections.

3.1 Morphology of the S${}_{\rm 74~{}MHz}>$ 1.5 Jy sample in FIRST

Our flux density limit of 1.5 Jy at 74 MHz corresponds to the knee in the source counts distribution. This flux limit is chosen to provide a large sample that is dominated by high excitation (predominantly FRII) radio galaxies.

To assess the dominant morphology of the sample, we compare our sample to that of van Velzen et al. (2015), who compiled a sample of 59,192 candidate FRII radio sources from the FIRST survey catalogue, using an algorithm that selected double, triple, or multiple sources with lobe-lobe separation of less than 60 arcseconds. This sample of “dubbeltjes" (small doubles) is considered to be a relatively clean sample of FRII radio galaxies on these angular scales. We cross-match our VLSSr flux limited sample with the dubbeltjes sample of van Velzen et al. (2015), using a matching radius of 30 arcseconds ($\lesssim 5$ expected random associations), and find 65$\%$ of our sample (2498 objects) are matched to candidate FRII radio galaxies in the catalogue of van Velzen et al. (2015). Not all of the van Velzen et al. (2015) sample are confirmed FRIIs. However, the selection of FRIIs in van Velzen et al. (2015) only includes radio galaxies up to 1 arcminute in angular size. Our sample includes sources with de-convolved major axes up to approximately 100 arcseconds. We therefore visually inspected a random sample of FIRST images centred on the VLSSr coordinates, and find that more than 50$\%$ have a double-lobed FRII-like appearance, while 10$\%$ are unresolved at the FIRST resolution (see Figure 2). This random sample is drawn from our VLSSr flux limited sample, and includes sources that are within the van Velzen sample, as well as some that are not in that sample.

Based on the preceding discussion, we can be confident that our $S_{\rm 74~{}MHz}>1.5$ Jy sample is dominated by FRII radio galaxies, as expected. We therefore place a robust lower limit on the number of FRII radio galaxies in our sample: $N_{\rm FRII}>0.5\times 3861>1930$ objects.

3.2 A limit on the fraction of FRIIs in our sample that are ultra-steep spectrum remnants

The total number of ultra-steep spectrum ($\alpha_{\rm 74~{}MHz}^{\rm 1400~{}MHz}>1.2$) sources in our sample is 57 (53 detections + 4 non-detections in NVSS). Therefore, as an absolute upper limit, fewer than 3$\%$ (57/1930) of FRII radio galaxies in our flux limited sample are ultra-steep spectrum remnants.

Of course, not all ultra-steep spectrum radio sources are remnant radio galaxies, and therefore we can place a more stringent upper limit on the fraction of remnant FRIIs in our sample by considering what fraction of the ultra-steep spectrum objects in our sample are active. We do this initially by considering the ultra-steep spectrum sample object’s morphology in FIRST. Of the 53 ultra-steep spectrum sources detected in NVSS, 10 are point-like in FIRST and may be high redshift radio galaxies, or core dominated systems in which the spectral index has been influenced by variable flux between the dates of observation in VLSS and NVSS surveys. A further 6 objects are triple systems, exhibiting two lobes plus an AGN core, indicating that these systems are still active (although some may be re-started systems, with remnant lobes and a re-born AGN). After removing the point-like and triple systems from the ultra-steep spectrum sample, we find that less than 41 objects are remnant radio galaxies, indicating that less than 2$\%$ of the FRIIs in our sample are candidate remnant radio galaxies.

The large core-detection fraction in the ultra-steep spectrum sample suggests that many of the ultra-steep spectrum objects without a core detection would have a detected core if we were to follow-up with sensitive high resolution observations. Therefore, the true number of ultra-steep spectrum FRII remnants in our sample is likely to be less than 41, because we detect such a large fraction of triples in our ultra-steep spectrum sample. Out of the 2498 FRII candidates from van Velzen et al. (2015) in our VLSSr sample, only 14$\%$ have a detected core. In comparison, of the 26 FRII candidates from van Velzen et al. (2015) in our ultra-steep spectrum sample, 23$\%$ have a detected core. Extrapolating from the 23$\%$ of ultra-steep sources with detected cores, it is conceivable that most of the ultra-steep spectrum objects in this flux limited sample are active, and would be found to have core emission if they were followed up with more sensitive high resolution observations. We note that the higher core fraction in the ultra-steep spectrum dubbeltjes from (van Velzen et al., 2015) could arise if the ultra-steep spectrum dubbeltjes tend to be larger angular size objects, enabling core detection in a larger number of cases.

Some of the ultra-steep spectrum triple sources may be restarted objects (eg. double-doubles) or high redshift radio galaxies, and a follow-up study of this sample of ultra-steep spectrum triples is warranted, particularly given that they are significantly lower angular size objects than most known double-doubles.

In summary, fewer than 2$\%$ of the FRII radio galaxies in our 74 MHz selected flux limited sample are candidate ultra-steep spectrum remnants. Due to the large core-detection fraction in our ultra-steep spectrum sample, we expect that many of the canidate remnants are in fact active, so that the true ultra-steep remnant fraction is likely to be significantly less than 2$\%$. Follow-up radio observations at higher resolution are required to confirm the nature of the ultra-steep spectrum objects, and which, if any, are indeed remnant FRII radio galaxies.

4 Simulating the population of active and remnant FRII radio galaxies

The second part of our study involves simulating the flux limited sample of FRII radio galaxies, including both the active and remnant phase, using an accurate model of spectral evolution. This Monte-Carlo approach allows us to compare the models with the empirical results using the exact same selection criteria. It also allows us, for the benefit of future studies, to investigate the efficiency of alternative selection criteria including such things as spectral curvature, high frequency spectra, and redshift limits. The degeneracy involved in modelling the radio spectrum of an individual remnant radio galaxy at one point in it’s life (Kaiser & Cotter, 2002) can be broken by observing many sources at different stages of their life, and modelling their spectral distributions.

In this section we present the relevant model equations used to construct the time dependent radio galaxy spectra, and we describe the procedure used to generate the simulated catalogues.

4.1 Flux density calculation assuming non-uniform magnetic field strength

In order to accurately simulate a population of radio galaxies including remnants, we must accurately trace the spectral evolution in both the active and remnant phases. Tribble (1991, 1993) showed that a non-uniform distribution of magnetic field strengths gives rise to aged synchrotron spectra that differ significantly from those obtained under the assumption of a uniform magnetic field strength. We therefore follow Tribble (1991, 1993) and Hardcastle (2013), and calculate the synchrotron spectra assuming that within each volume element of the lobes the magnetic field is distributed according to a Guassian-random field. That is, at each point within a volume element, the cartesian components of the field are each drawn independently from a Gaussian distribution with mean of zero. This magnetic field configuration is expected to arise from homogenous, isotropic turbulence (Tribble, 1991; Hardcastle, 2013). We further assume that each volume element within the lobes has the same magnetic field distribution, that the electron energy distribution is independent of the local magnetic field strength, and that the particle pitch angle distribution is isotropic. With these assumptions, we can calculate the synchrotron flux density as

[TABLE]

where $N(\gamma)=\frac{dN}{d\gamma}$ is the volume integrated electron energy distribution, $\bar{F}(y)$ is the angle-averaged synchrotron function (see Section C), $p_{B}$ is the probability distribution for the magnetic field strength, $\epsilon_{0}$ is the permittivity of free space, $c$ is the speed of light, $m_{e}$ is the electron mass, and $e$ is the electron charge.

For a Guassian-random field, the probability distribution for the magnetic field strength $p_{B}$ is the Maxwell-Boltzmann distribution (Hardcastle, 2013):

[TABLE]

where

[TABLE]

and $B_{0}$ specifies the mean magnetic energy density, defined such that

[TABLE]

4.2 Radio source evolution

In the current work, we are focused on FRII radio galaxies, and therefore we assume that during the active phase the radio sources evolve according to the self-similar dynamical model of Kaiser & Alexander (1997). During the remnant phase, the lobe evolution is not well understood, however, the true lobe expansion is likely to be bounded by two extreme cases of (i) maximal energy driven expansion and (ii) no expansion. For the maximal expansion rate of case (i) we consider the Sedov-like expansion described by Kaiser & Cotter (2002). Sedov-like expansion refers to the situation in which the bow shock surrounding the lobes is driven by the adiabatic expansion of the lobes. This provides a solution that is similar to, but not entirely analogous to the Sedov solution for a point explosion, and is therefore referred to by Kaiser and Cotter as "Sedov-like". In order to apply the models of Kaiser & Alexander (1997) and Kaiser & Cotter (2002), we must assume that the radio lobes expand into a power-law radial density profile, in which the ambient density, $\rho$, scales with the radial distance, $r$, according to

[TABLE]

4.2.1 Case (i), Sedov-like expansion in the remnant phase

In our model, for case (i), the evolution of the radio source volume is a piece-wise power law, with:

[TABLE]

where $t_{\rm on}$ is the length of the active phase. Our assumption of piece-wise power-law volume evolution means that the integrated electron energy distribution $N(\gamma)$ does not depend on the absolute value of the volume, only the exponents of the time evolution. For this reason, we have purposefully left the above expression as a proportionality. We assume that the magnetic field is isotropically distributed (“tangled") on all scales, and therefore can be treated as a magnetic “fluid" with adiabatic index $\Gamma_{B}=4/3$ (Leahy, 1991). In this case, the magnetic field evolves according to

[TABLE]

(Kaiser & Alexander, 1997; Kaiser & Cotter, 2002), where $Q$ is the jet power and $Q_{0}=10^{39}$ W is a normalisation constant. We note that for $\beta\approx 1.5-2$, the dependence of magnetic field strength on jet power is extremely weak, and therefore $Q_{0}$ does not strongly affect the results. We assume that the particle energy distribution injected into the lobes is a power-law, such that

[TABLE]

where $q_{0}$ is proportional to the jet power, and a number of other assumed model parameters, as described in Section A.3.

4.2.2 Case (ii), no expansion in the remnant phase

The only difference in our model for case (i) and case (ii) is in equations 6 and 7 at $t>t_{\rm on}$. For completeness, they are defined below:

[TABLE]

and

[TABLE]

4.3 Volume integrated electron energy distribution

Equations 6, 7 and 8 completely define our radio source model. The volume integrated electron energy distribution can be obtained from the completely general solution to the continuity equation (see Appendix A.2)

[TABLE]

where $t_{\rm i,max}=t$ for active sources, and $t_{\rm i,max}=t_{\rm on}$ for remnant sources.

[TABLE]

where $B_{\rm CMB}=0.325(1+z)^{2}$ nT is the equivalent magnetic field strength of the CMB, and the integration limit $t_{\rm i,min}$ is given by

[TABLE]

The parameter $t_{\rm i,min}^{*}$ corresponds to the injection time at which a particle injected with Lorentz factor $\gamma_{\rm i,max}$ will have cooled to Lorentz factor $\gamma$ at time $t$, and is given by the solution to the equation

[TABLE]

where $\gamma_{\rm i,max}$ is the maximum electron Lorentz factor of the particle distribution injected into the lobes (see Appendix A.2).

4.4 Simulation Approach

To generate a mock catalogue of radio galaxies, several of the model parameters are sampled from probability distributions, while others remain fixed. Each of the free parameters and the relevant distributions are discussed below.

4.4.1 Redshift

Redshifts are sampled from a probability distribution given by

[TABLE]

where $\rho(z)$ is the volume density of radio galaxies as a function of redshift and $\frac{dV}{dz}$ is the differential comoving volume element for a spherical shell. For high power FRII radio galaxies, the volume density $\rho(z)$ is typically taken to be a Gaussian (Blundell et al., 1999; Willott et al., 2001; Grimes et al., 2004). We assume

[TABLE]

with $z_{h0}=1.95$ and $z_{h1}=0.55$ (Grimes et al., 2004). The comoving volume element for a flat Universe ($\Omega_{k}=0$) is

[TABLE]

(Hogg, 1999).

4.4.2 Radial density profile exponent, $\beta$

As discussed, we assume that our mock radio galaxies expand into a power-law density profile with $\rho_{\rm ext}\propto r^{-\beta}$. We assume that the exponent $\beta=1.9$ for all of our mock radio galaxies. Lower values of $\beta$ increase the fraction of remnants, but not by a large factor, due to the less rapid evolution of volume and magnetic field.

4.4.3 Active lifetime $t_{\rm on}$

The typical active lifetime of FRII radio galaxies remains poorly constrained. Estimates based on self-similar dynamical models range from $\gtrsim$ 10 Myr (Bird et al., 2008; Kapińska et al., 2012), to 200 Myr (Antognini et al., 2012).

Analysis of the length asymmetry of the most powerful double lobed radio sources indicates that the lobe advance speeds are typically a few percent of the speed of light, and not more than 0.15c (Scheuer, 1995), implying a typical active lifetime of a few tens of Myr.

Estimates based on spectral ageing (Alexander & Leahy, 1987; Liu et al., 1992) systematically underestimate dynamical ages by a significant factor, and both dynamical and spectral age estimates have an uncertain relationship to the true source age (Eilek et al., 1997; Blundell & Rawlings, 2000; Kaiser, 2005; Hardcastle, 2013).

For the purposes of our simulated catalogues, we assume that the active lifetimes follow a truncated log-normal distribution, with mean $\langle\log(t_{\rm on})\rangle=7.5$ and standard deviation $\sigma_{\log(t_{\rm on})}=0.1$, truncated such that $7.3<\log(t_{\rm on})<8.3$, where $t_{\rm on}$ is specified in years. The active lifetime $t_{\rm on}$ refers to the time period during which the lobes are fed with fresh electrons. Increasing the mean active lifetime by a factor $f$ will cause a decrease in the simulated catalogue remnant fraction by a factor $\lesssim f$, while decreasing the mean active lifetime by a factor $f$ will increase the remnant fraction by a factor $\lesssim f$.

4.4.4 Age at which the source is observed, $t_{\rm obs}$

Source ages are sampled from a uniform distribution between $t_{\rm obs,min}=0.1$ Myr and $t_{\rm obs,max}=200$ Myr, which is several times the active lifetime. The results are insensitive to the assumed value for $t_{\rm obs,max}$, because the age distribution of remnants declines very steeply, as shown in Figure 4 (b). The results are insensitive to the assumed value for $t_{\rm obs,min}$, because $t_{\rm obs,min}\ll t_{\rm obs,max}$ and source age $t_{\rm obs}$ is sampled uniformly.

4.4.5 Magnetic field normalisation $B_{0}$

One of the most important and influential parameters of our model is the normalisation of the magnetic field strength, $B_{0}$ in Equation 7. Croston et al. (2005) measured lobe magnetic field strengths in a sample of 33 powerful FRII radio galaxies, based on inverse Compton modelling of the observed X-ray emission. They find that the magnetic field strength in their sample is a strong function of source size. For sources greater than 300 kpc, Croston et al. (2005) find the median magnetic field strength is 0.6 nT, and reaches 0.2 nT for the largest sources in that sample. We therefore fix $B_{0}=0.6$ nT at time $\log(t_{\rm years})=7.3$ corresponding to the lower limit on the active lifetime for sources in our mock sample. In this way, we ensure that the magnetic field strength at the end of each source’s active lifetime is $\lesssim 0.6$ nT, consistent with the results of Croston et al. (2005). In equation 7, we have set $Q_{0}=10^{39}$ W. However, we note that the dependence of magnetic field strength on jet power is extremely weak, and $Q_{0}$ is therefore not of great importance.

4.4.6 Energy injection index, $a$

The injection index is represented by the parameter $a$ in equation 8. We sample the injection index for each source from a truncated Gaussian distribution with $2.0<a<2.4$, mean $\bar{a}=2.2$ and standard deviation $\sigma_{a}=0.2$. This distribution of injection indices results in a relatively broad distribution of “observed" spectral index in our mock catalogues, ranging between $0.5<\alpha_{74}^{1400}<1.2$ with a peak at $\alpha_{74}^{1400}\sim 0.8$, providing a good match to the spectral index distribution observed in our VLSSr sample.

Note that with this injection index distribution, even the oldest active sources cannot have spectral index $\alpha>1.2$ between any pair of frequencies, unless those frequencies are near to the frequency corresponding to the cutoff in the electron distribution. Therefore, active sources do not contaminate our ultra-steep spectrum sample in our mock catalogue. To explain this in more detail, consider a source with injection index $a=2.4$ – the maximum possible value that could be drawn from the truncated Gaussian described above. Imagine this source is observed when it is very old, so that the radiative cooling break of $\Delta\alpha=0.5$ occurs well below the lowest frequency we have observed. In that case, the spectral index will be $\alpha=1.2$, but cannot get any steeper with age.

4.4.7 Minimum/maximum injected Lorentz factor, $\gamma_{\rm min,max}$

We assume $\gamma_{\rm min}=100$ for all mock radio galaxies, and $\gamma_{\rm max}=5\times 10^{6}$. The maximum and minimum injected Lorentz factors do not strongly affect the results presented here, since the electrons at $\gamma_{\rm min}$ and those at $\gamma_{\rm max}$ emit at frequencies well outside the observed range. They do however affect the scaling between jet power and particle injection rate $q_{0}$ (see Section A.3. ).

Our choice for low energy electron cutoff is based on the observation in several sources of a low-frequency flattening in the hotspot spectra, consistent with a low energy cutoff in the electron distribution at a Lorentz factor of $\gamma_{\rm min}\sim 100-700$ (Godfrey et al., 2009, and references therein). Such values of $\gamma_{\rm min}$ are likely to be the result of dissipation of jet bulk kinetic energy (Godfrey et al., 2009). In the case of Cygnus A, absorption is at least partially responsible for the low frequency turnover in the hotspot spectra (McKean et al., 2016), however, this does not rule out the possible involvement of a low energy cutoff in the hotspots of Cygnus A. The interpretation of the hotspot turnover in Cygnus A remains problematic (McKean et al., 2016).

4.4.8 Jet power, $Q_{\rm jet}$

The jet power is sampled from a power-law probability distribution, with

[TABLE]

We assume $n_{Q}=2.3$, which is an average of the values derived by Blundell et al. (1999) (2.6) and Wang & Kaiser (2008) (2.0). We note that the results presented here are not strongly dependent on the assumed value of $n_{Q}$, for reasonable departures from our assumed value.

4.4.9 Electron energy fraction, $\epsilon_{e}$

To calculate the particle injection rate $q_{0}$, we must specify the fraction of jet power that is converted to the internal energy of the relativistic electron population, $\epsilon_{e}$ (see Equation 53). Here we assume that the jet power is equally distributed between magnetic energy, relativistic electron energy, and the energy of non-radiating particles, so that $\epsilon_{e}=1/3$. Our conclusions are not sensitive to the value of $\epsilon_{e}$.

4.4.10 Ratio of hotspot pressure to lobe pressure, $\frac{p_{\rm HS}}{p_{\rm lobe}}$

Another parameter that is necessary to calculate the particle injection rate into the lobes, is the ratio of hotspot pressure to lobe pressure (see Equation 53). Here we assume $\frac{p_{\rm HS}}{p_{\rm lobe}}=10$, however we note the very weak dependence of $q_{0}$ on the value of $\frac{p_{\rm HS}}{p_{\rm lobe}}$. Again, our conclusions are not sensitive to the assumed value of $\frac{p_{\rm HS}}{p_{\rm lobe}}$.

4.5 Contribution of the jets and hotspots in active sources

After the jets stop supplying energy to the hotspots, the bright emission from them will disappear in the order of a sound-crossing time ($\lesssim 1$ Myr for a hotspot with a diameter of a few kpc). If the hotspots provide a significant fraction of the total source flux at 74 MHz, then the rapid disappearance of the hotspots at the start of the remnant phase could help to explain the rapid disappearance of remnant FRII radio galaxies from our flux limited samples.

Jenkins & McEllin (1977) defined the hotspots as regions 15 kpc in diameter, and found that the ratio of hotspot to total luminosity is strongly correlated with the total luminosity, and furthermore, that the hotspots can often contribute more than 90$\%$ of the source flux density at 178 MHz, particularly at high luminosities.

In direct contrast to this result, is the more recent study by Mullin et al. (2008) of a sample of 100 FRII 3CRR radio galaxies with $z<1$. Mullin et al. (2008) are able to more accurately identify the hotspot regions than Jenkins & McEllin (1977), and find no correlation between hotspot prominence and the total radio luminosity in this sample (see their figure 32). Furthermore, Mullin et al. (2008) find that the hotspots contribute typically a few percent to a few tens of percent to the total radio luminosity at 178 MHz222Mullin et al. (2008) define hotspot prominence as the ratio of hotspot luminosity at 8.4 GHz to the total source luminosity at 178MHz. We have converted the hotspot prominence of Mullin et al. (2008) to a hotspot emission fraction, or compactness in the terminology of Jenkins & McEllin (1977) by multiplying the hotspot prominence values by a factor of $(8400/178)^{-\alpha}$ assuming $\alpha=0.8$..

The median hotspot prominence (summed for both the north and south hotspot) for the sample of Mullin et al. (2008) is 0.008. We convert this to a hotspot fraction (the ratio of total source luminosity contributed by the hotspots at 74 MHz) by multiplying the hotspot prominence by a factor of $(74/8400)^{-\alpha}$ assuming $\alpha=0.8$. This gives a median hotspot fraction of 35$\%$. However, this does not account for the fact that hotspot spectra often become much flatter towards low frequency, particularly in high luminosity radio galaxies (Leahy et al., 1989; Godfrey et al., 2009; McKean et al., 2016). We therefore do not attempt to account for the contribution to the source flux from hotspot related emission. We simply note that if the source flux density drops by a factor $f_{\rm drop}$ at the end of the active phase due to the rapid disappearance of the hotspots, then we will over-estimate the fraction of remnant to active sources by a factor of $f_{\rm drop}^{p-1}$, where $p\approx 2.3$ is the slope of the luminosity function.

The jet prominence is typically much lower than the hotspot prominence, by an order of magnitude or more (Mullin et al., 2008). While it is true that in a few percent of cases, the jet prominence is comparable or greater than the hotspot prominence, in the vast majority of cases, the jet prominence is a negligible component of the total source flux density.

4.6 Procedure for generating mock catalogues

The mock catalogues are created in several steps, as described below.

First, we generate several million radio galaxies, each with the model parameters fixed or sampled from their corresponding probability distributions as described in the preceding sections. 2. 2.

For each source, we calculate an upper limit to the flux density at the sample selection frequency (74 MHz) using a fast, analytic expression. 3. 3.

We apply the flux cut to the sample, using the upper limits calculated in the previous step: All sources for which the calculated flux upper limit at the selection frequency is below the flux limit are removed from the mock catalogue. 4. 4.

For the remaining sources, we calculate the model radio galaxy spectrum accurately using numerical integration of equations 11 and 1. 5. 5.

We again apply a flux cut to the sample, this time using the accurate model flux densities as the basis of the flux cut. Only those sources whose flux densities lie above the flux limit remain in the sample.

The flux density upper limit in step (ii) described above, is obtained by using an analytic approximation to equation 11, along with the $\delta-$function approximation to the synchrotron emission spectrum. The analytic approximation to equation 11 is derived by neglecting the term $(1-\gamma/\gamma^{*})^{a-2}$. That is, we replace the following equation (Equation 11)

[TABLE]

with the following integral,

[TABLE]

which in the case of power law volume evolution, has an analytic solution. Since the term $(1-\gamma/\gamma^{*})^{a-2}<1$, Equation 19 provides an upper limit on $N(\gamma)$, and therefore an upper limit on the flux density at any frequency.

5 Modeling Results

We carried out two Monte-Carlo simulations, each one representing a bound on the possible evolution scenarios in the remnant phase, and each one differing only in the prescribed evolution of lobe volume and magnetic field strength in the remnant phase (see the description of these scenarios in section 4.2). In the remnant phase, for case (i), we assume that the lobes remain over pressured with respect to the ambient medium, and evolve in a Sedov-like manner as described by Kaiser & Cotter (2002). This is the maximal expansion rate that can be achieved by inactive radio galaxies. In case (ii) we assume that there is no lobe expansion in the remnant phase.

In the active phase, for both of our simulations, we assume the lobes evolve according to self-similar expansion models of Kaiser & Alexander (1997).

We note that our model is only applicable to radio galaxies with a single active phase: re-started radio galaxies will show different characteristic luminosity and spectral evolution.

5.1 Results of Monte-Carlo simulation with Maximal (sedov-like) expansion in the remnant phase (Case (i))

In Figures 4 and 5 and Tables 3 and 4 we summarise the results of our simulations of active and remnant radio galaxy populations with maximal (sedov-like) expansion in the remnant phase.

Our simulated catalogue contains a total of 10,038 sources333We purposefully simulated a larger catalogue of sources than was obtained in our VLSSr selected sample in order to increase the fidelity of the simulated distributions., of which 9499 are active (94.6$\%$), and 539 are remnant (5.4$\%$). However, of the 539 remnants, only 187 (1.8$\%$ of the entire sample) have ultra-steep spectra between 74 MHz and 1400 MHz with $\alpha_{\rm 74MHz}^{\rm 1400MHz}>1.2$. The fraction of ultra-steep spectrum remnants in our mock sample is therefore consistent with our upper limit of $2\%$ on the ultra-steep spectrum FRII remnant fraction obtained in section 3. However, this model predicts almost two times as many remnants would have been missed by our remnant selection criterion. The 74 MHz to 1400 MHz spectral index is clearly an inefficient selection criterion, with a selection efficiency of only $35\%$ (187/539). Spectral indices measured at higher frequency are significantly more efficient at remnant selection. If remnant selection incorporated higher frequency data using the criterion $\alpha_{\rm 1.4GHz}^{\rm 5GHz}>1.2$, our Monte-Carlo simulation predicts that we would more than double the number of remnants selected, and achieve $70\%$ remnant selection efficiency. If we could have incorporated both 5GHz and 10 GHz data into our analysis, and selected remnants based on the criterion $\alpha_{\rm 5GHz}^{\rm 10GHz}>1.2$, our Monte-Carlo simulation predicts a selection efficiency of 82$\%$.

One of the problems with selecting remnants based on ultra-steep spectra alone is that there are other types of radio source that can satisfy the ultra-steep spectrum criterion, such as high-redshift radio galaxies (HzRGs). To overcome this problem, some authors have used spectral curvature as a more robust indicator of remnant radio galaxies (eg. Murgia et al., 2011). Spectral curvature is calculated as the difference between a high frequency spectral index and low frequency spectral index, and requires flux density measurements at at least three different frequencies. We have considered the spectral curvature selection in our mock sample using measurements at three frequencies, up to 5 GHz ($\alpha_{\rm 1400~{}MHz}^{\rm 5000~{}MHz}-\alpha_{\rm 74~{}MHz}^{\rm 1400~{}MHz}>0.5$), and also using four frequencies up to 10 GHz ($\alpha_{\rm 1400~{}MHz}^{\rm 5000~{}MHz}-\alpha_{\rm 74~{}MHz}^{\rm 1400~{}MHz}>0.5$), achieving remnant selection efficiency of 54$\%$ and $71\%$, respectively.

Figure 4 (a) demonstrates that the flux density distribution for the remnant population has the same slope as that of the active population, despite their different luminosity evolution. This implies that the overall remnant fraction is independent of the flux limit of the sample.

Figure 4 (b) demonstrates the very sharp decline in the number of remnants as a function of source age. This implies that remnants in flux limited samples tend to be relatively “new" remnants, as expected due to the rapid luminosity evolution that arises from a combination of decreasing magnetic field strength along with the adiabatic and radiative losses during the remnant phase. This tendency for remnant radio galaxies in flux limited samples to be young is the reason that many remnants in our Monte-Carlo simulations do not have ultra-steep spectra in the observed frequency range. Most importantly, unlike the luminosity distribution, the age distribution of a flux limited sample of remnants can provide meaningful constraints on the luminosity evolution during the remnant phase. To illustrate this, we consider the following simplified example. Suppose that sources are “born" with some peak luminosity $L_{\rm peak}$, and evolve according to a piecewise power law with

[TABLE]

Now assume that sources are “born" with luminosity $L_{\rm peak}$ at a rate given by the birth function

[TABLE]

Then the age distribution of remnant radio galaxies at a given luminosity $L_{\nu}$ is given by

[TABLE]

Thus, given an estimate of the “birth" luminosity function (i.e. $p$ in the above model), the age distribution can be used to constrain the luminosity evolution in the remnant phase.

Figures 4 (c) and (d) demonstrate the existence of many remnant radio galaxies with “normal" spectral index values in our mock catalogue, as well as the broad tail of ultra-steep spectrum remnants.

Figure 4 (e) demonstrates that high redshift remnant radio galaxies are extremely rare, and low redshift remnant radio galaxies are predicted to be much more common. This is due to the combination of two factors: (1) the increased rest-frame frequency at higher redshifts corresponds to higher energy electrons, and consequently corresponds to a faster radiative cooling rate; and (2) the increased energy density of the CMB at higher redshifts causes a faster radiative cooling rate from inverse Compton scattering of the CMB. The number of remnants per redshift bin is approximately constant up to $z\gtrsim 1$, but decreases dramatically for $z\gtrsim 1.5$. In Table 4 we list results for the low redshift ($z<0.5$) subsample of our simulated catalogues. In the simulated catalogue with Sedov-like expansion in the remnant phase, 20$\%$ of the flux limited FRII sample is predicted to be remnant, and of those remnants, approximately half will show ultra-steep spectra between 1.4 GHz and 5 GHz. It is clear from Table 4 and Figures 4 (e) and 5 that there is a significant advantage to be gained by studying the remnant population at low redshifts, particularly when high frequencies (5 and 10 GHz) are included in the analysis.

5.2 Results of Monte-Carlo simulation with no expansion in the remnant phase

We repeated the Monte-Carlo simulations with no expansion in the remnant phase, so that

[TABLE]

and

[TABLE]

This model is not physically realistic: it is unlikely that all sources will reach pressure equilibrium right at the moment the central engine shuts off. However, it is a useful exercise to demonstrate the strength of the effect of expansion in the remnant phase. It also clearly demonstrates that models with no expansion in the remnant phase are in great conflict with observations. This is important, given the evidence that FRII radio galaxies might reach pressure equilibrium with the surroundings before the end of their active life.

The results for this model with no expansion in the remnant phase are presented in figures 8 and 6 and Tables 3 and 4. It is immediately clear, as expected, that the remnant fraction is greater in this simulated catalogue, particularly at low redshift ($z<0.5$), where the model predicts that nearly two-thirds of the flux limited sample are remnants, and that a third of the flux limited sample are ultra-steep spectrum remnants with $\alpha_{\rm 74~{}MHz}^{\rm 1400~{}MHz}>1.2$ (see table 4 and Figure 6).

6 Discussion and Conclusions

We have carried out an empirical study based on a flux limited sample from the VLSSr radio catalogue, in order to place a limit on the occurrence of remnant FRII radio galaxies in a 74 MHz flux limited sample. We have also performed Monte-Carlo simulations of the population of active and remnant FRII radio galaxies to assess whether models of remnant lobe evolution are consistent with the observed remnant fraction. Our main conclusions may be summarised as follows:

In our VLSSr selected sample, fewer than 2$\%$ of FRII radio galaxies with 74 MHz flux density greater than 1.5 Jy are ultra-steep spectrum remnants with $\alpha_{\rm 74~{}MHz}^{\rm 1400~{}MHz}>1.2$. 2. 2.

Our Monte-Carlo simulation with Maximal (sedov-like) expansion in the remnant phase produced a remnant fraction of $2\%$, marginally consistent with the upper limit described above. 3. 3.

Our Monte-Carlo simulation with Maximal (sedov-like) expansion in the remnant phase predicts the existence of nearly twice as many remnants with “normal" spectra ($\alpha_{\rm 74~{}MHz}^{\rm 1400~{}MHz}<1.2$) as there are ultra-steep spectrum remnants in our sample. 4. 4.

The above conclusion can be phrased another way: the ultra-steep selection criterion $\alpha_{\rm 74~{}MHz}^{\rm 1400~{}MHz}>1.2$ is not efficient at selecting remnant radio galaxies, with a selection efficiency of only $\sim 35\%$. Higher frequency spectral indices are significantly more efficient at remnant selection. Remnant selection based on the criterion $\alpha_{\rm 1.4GHz}^{\rm 5GHz}>1.2$ increases the selection efficiency to $70\%$, and remnant selection based on the criterion $\alpha_{\rm 5GHz}^{\rm 10GHz}>1.2$ increases the selection efficiency to 82$\%$. For redshifts less than 0.5, the number of identified ultra-steep spectrum remnants increases by a factor of 4 when using $\alpha_{\rm 5~{}GHz}^{\rm 10~{}GHz}$ as opposed to $\alpha_{\rm 74~{}MHz}^{\rm 1400~{}MHz}$ (Table 4). 5. 5.

The remnant fraction increases rapidly towards low redshift. This is the result of (1) the increased rest-frame frequency at higher redshifts and (2) the increased energy density of the CMB at higher redshifts, resulting in and increase in the radiative cooling rate from inverse Compton scattering of the CMB. 6. 6.

The model predicts an ultra-steep remnant fraction approaching 10$\%$ at redshifts $z<0.5$, when considering ultra-steep selection based on higher frequency data ($\alpha_{\rm 1.4~{}GHz}^{\rm 5~{}GHz}>1.2$). 7. 7.

The age distribution of remnant radio galaxies in flux limited samples is a steeply decreasing function of sources age, indicating that most remnant radio galaxies in flux limited samples are young remnants (Figures 4 (b) and 8 (b)). Due to the steep remnant age distribution, incorporating higher frequency data into the analysis will be more important in future studies than incorporating lower frequencies. 8. 8.

The age distribution of remnant radio galaxies in a flux limited sample can constrain the luminosity evolution of the remnant phase. The luminosity distribution of remnants cannot. 9. 9.

The spectral index distribution of remnant radio galaxies peaks at “normal" spectral index values. This is a direct result of the age distribution: most remnants are young in flux limited samples (Figure 4 (c)). 10. 10.

In the idealised situation that we have modelled, and with the selection criteria we have used to identify candidate remnants, the high frequency spectral index is more efficient at selecting remnants than the spectral curvature. However, when considering heterogeneous samples as obtained for example from flux limited samples, the spectral curvature is likely to be more robust and may result in fewer false remnant candidates. 11. 11.

The remnant fraction is independent of flux limit (Figure 4 (a)). Therefore, going to fainter flux limits will not increase the remnant fraction.

As discussed in Section 1, several studies suggest that FRII radio lobes may reach pressure equilibrium before the end of their lifetime. However, we have shown that models without rapid remnant phase expansion significantly over-predict the FRII remnant fraction. Rapid luminosity evolution in the remnant phase resulting from Sedov-like expansion is required to match the low observed remnant fraction in our flux limited sample. Our results imply that either the previous evidence for internal/external pressure equilibrium in FRII radio galaxy lobes is flawed, or alternative mechanisms other than adiabatic expansion are required explain the low remnant fraction in our flux limited sample.

Acknowledgements

The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Advanced Grant RADIOLIFE-320745.

Appendix A The Generalised Continuous Injection Model for Lobed Radio Galaxies

We describe a mathematical framework for modelling the integrated spectra of radio galaxies in various phases of their life, under the following set of assumptions.

The electron energy distribution is independent of the local magnetic energy density. 2. 2.

The cooling rate is a function of time and particle energy only. Specifically, this assumption implies that the cooling rate is independent of position within the lobes. This is a common assumption used in modelling integrated radio spectra, and is valid provided that the average magnetic energy density (averaged over the cooling length scale) is independent of position within the lobes.

This mathematical framework allows for arbitrary evolution of lobe volume and magnetic field strength, it allows for arbitrary, and time-dependent injection spectra (important for modelling objects in which jet power evolves with time), and it allows for a distribution of magnetic field strengths within each volume element of the lobes.

However, this type of model is only valid for lobed radio galaxies, in which the emission is dominated by the lobes, and the mean magnetic energy density is approximately constant throughout the lobe volume. This model is therefore not suited to the modelling of tailed or naked jet FRI radio galaxies.

We present a simplified version of the generalised model, which is relevant to modelling objects with multiple power-law phases of lobe evolution, such as remnant radio galaxies, as well as cluster radio-relics (Enßlin & Gopal-Krishna, 2001).

A.1 Integrated flux density

Consider a small volume of lobe plasma $\Delta V$ in which the number density of relativistic particles per unit volume and per unit Lorentz factor is given by $\frac{dN}{d\gamma dV}=n(\gamma)$, and in which the magnetic field is described by $p_{\Theta}$ and $p_{B}$, where $p_{\Theta}$ and $p_{B}$ are the probability distributions of the pitch angle $\Theta$ and the magnetic field strength $B$, respectively. Hardcastle (2013) give an expression for the volume-averaged emissivity from such a volume element (see Hardcastle, 2013, for details):

[TABLE]

where $e$ and $m_{e}$ are the charge and mass of the electron, respectively, $c$ is the speed of light in vacuum, $\epsilon_{0}$ is the permttivity of free space, $x=\frac{\nu}{\nu_{c}}$, $\nu$ is the rest-frame frequency, $\nu_{c}=\frac{3}{4\pi}\Omega_{0}\gamma^{2}\sin\Theta$ is the characteristic synchrotron frequency, $\Omega_{0}$ is the non-relativistic gyrofrequency which in SI units is given by $\Omega_{0}=eB/m_{e}$, and the synchrotron function $F(x)=x\int_{x}^{\infty}K_{5/3}(z)dz$, where $K_{5/3}(z)$ is the modified Bessel function of order 5/3.

If we assume that $n(\gamma)$, $B$ and $\Theta$ are all independent (assumption (i)), then we can separate the integrals, and the equation becomes

[TABLE]

Let us define the parameter

[TABLE]

We can write the inner integral over pitch-angle as (see Appendix C)

[TABLE]

Then

[TABLE]

The flux density of the lobe is

[TABLE]

[TABLE]

If we assume that the magnetic field distribution $p_{B}$ and $p_{\Theta}$ is the same within each volume element $\Delta V$ throughout the lobes (assumption (ii)), then we can switch the order of integration, so that

[TABLE]

Let us define

[TABLE]

Then

[TABLE]

For a Guassian-random field, the probability distribution for the magnetic field strength $p_{B}$ is the Maxwell-Boltzmann distribution (Hardcastle, 2013):

[TABLE]

where

[TABLE]

and $B_{0}$ specifies the mean magnetic energy density, and is defined such that

[TABLE]

Then, for a Gaussian-random magnetic field distribution,

[TABLE]

[TABLE]

It therefore remains to specify the volume integrated electron energy distribution as a function of time, $N(\gamma,t)$. We do so by solving the continuity equation for the electron population within the lobes, as described in the following section.

A.2 Volume integrated electron distribution $N(\gamma)$

Consider a radio lobe of volume V(t) with mean magnetic energy density $U_{B}(t)=B_{0}(t)^{2}/2\mu_{0}$, into which relativistic particles are injected at a rate of $\frac{dN}{d\gamma_{i}dt_{i}}=\sigma_{\rm lobe}(\gamma_{i},t_{i})$. Provided that the cooling rate $\frac{d\gamma}{dt}$ is a function of time and particle energy only (assumption (ii)), the evolution of the volume integrated particle energy distribution is described by the continuity equation

[TABLE]

The general solution to the continuity equation is

[TABLE]

and is valid for any functional form of the cooling rate $d\gamma/dt$. Here we consider only radiative and adiabatic losses (it is straightforward to modify the following to include acceleration and/or additional loss terms). In this case we can write the cooling rate as

[TABLE]

where, in both c.g.s. and S.I. units,

[TABLE]

and $U(t)$ is the sum of the magnetic energy density and the CMB energy density. For arbitrary volume and magnetic field evolution, the injected Lorentz factor, $\gamma_{i}$, injected at time $t_{i}$, is related to the particle’s Lorentz factor $\gamma$ at time t, via

[TABLE]

and the derivative is

[TABLE]

The completely general solution for arbitrary V(t), $B_{0}(t)$ and $\sigma_{\rm lobe}(\gamma_{i},t_{i})$ is

[TABLE]

where

[TABLE]

and

[TABLE]

and the integration limit $t_{\rm i,min}$ is given by

[TABLE]

and $t_{\rm i,min}^{*}$ is given by the solution to the equation

[TABLE]

where $\gamma_{\rm i,max}$ is the maximum electron Lorentz factor of the particle distribution injected into the lobes. In the case of an infinite power-law, $\gamma_{\rm i,max}\rightarrow\infty$, in which case the first term on the RHS of equation 48 is equal to zero. For injection times earlier than $t_{\rm i,min}^{*}$, all particles will have cooled to below the Lorentz factor $\gamma$ at time t.

Equations 44 - 48 represent the generalised continuous injection model. The standard $CI_{\rm off}$ model of Komissarov & Gubanov (1994), as well as the equations of Murgia et al. (1999) and the model of Enßlin & Gopal-Krishna (2001) can all be derived from the above expression under the relevant simplifying assumptions.

The above expression is valid for an arbitrary particle injection spectrum $\sigma_{\rm lobe}(\gamma_{i},t_{i})$, and arbitrary evolution of the cooling function $d\gamma/dt$. We note that $\sigma_{\rm lobe}(\gamma_{i},t_{i})$ is the electron energy distribution that is “injected" into the lobe. Typically the particles are accelerated (injected) in regions of higher pressure than the lobes (i.e. the hotspots). The relationship between the particle distribution that is injected into the lobes and the particle distribution that is accelerated in the hotspots is discussed below in Section A.3.

A.3 The relation between $\sigma_{\rm lobe}(\gamma_{i},t_{i})$ and Jet Power

The particle injection rate into the lobe, per unit energy, is defined by

[TABLE]

We assume that the electron energy spectrum injected into the lobes can be approximated by a power law

[TABLE]

We can relate the proportionality constant $q_{0}$ to the jet power, if we make a number of further model assumptions. We assume that the particles are injected into a high pressure region (the hotspot) with pressure $p_{\rm hs}$ and flow into the lobe at some lower pressure $p_{\rm lobe}$. The particles lose energy as they flow from the high pressure injection site into the lobes via adiabatic expansion and radiative cooling, and may gain energy due to re-acceleration processes (eg. Orienti et al., 2012). Due to the complexity and uncertainty surrounding the radiative cooling and re-acceleration in the head of the source, we consider only adiabatic expansion losses in this region, but note that alternative specific models may be considered (Blundell et al., 1999; Manolakou & Kirk, 2002).

Suppose that some (constant) fraction of the jet power, $\epsilon_{e}$, is converted to the internal energy of the relativistic particle population, and that the energy per particle is reduced due to adiabatic expansion losses by a factor

[TABLE]

The energy injection rate (into the lobes) is

[TABLE]

where the integration is over the particle distribution in the lobes, not the particle distribution injected into the hotspots. We can then write

[TABLE]

where

[TABLE]

and

[TABLE]

Appendix B Limiting cases of the generalised Continuous Injection Model

B.1 Power-law injection spectrum

It is typically assumed that a power-law energy spectrum is injected, with

[TABLE]

In this case,

[TABLE]

The effects of cooling and source evolution enter via $\gamma_{*}$ as defined in Equation 46. The different limiting cases are obtained by evaluating the integral 46 analytically under different sets of assumptions.

B.2 Single power law phase

Let the source volume and magnetic field strength evolve according to power-laws

[TABLE]

The energy density of the cosmic microwave background is $U_{\rm CMB}=4.2\times 10^{-12}~{}(1+z)^{4}$ J m*-3*. From Equation 46, it can be shown that:

[TABLE]

This is equivalent to equation 7 of Enßlin & Gopal-Krishna (2001), and equation 7 of Kaiser et al. (1997).

B.3 Multiple power law phases

Suppose we have K power-law phases of evolution, bounded by the times $(t_{0},t_{1},t_{2},...,t_{K})$ for phases (1, 2, …, K) and each $j^{\rm th}$ phase has its own $n_{V_{j}}$ and $n_{B_{j}}$. Then for each $j^{\rm th}$ phase we can define $\gamma_{*,j}$ as

[TABLE]

The total effect of all of the phases is given by

[TABLE]

The above is a slightly different, but equivalent form to that given by Enßlin & Gopal-Krishna (2001).

Suppose that injection occurs only within the first of several phases, i.e. prior to $t_{1}$, then

[TABLE]

where $\gamma_{*}(\tau,t)$ is given by the above expression, with $t_{0}=\tau$ and $t_{K}=t$.

To test the formalism described in this section, we compared the output with a numerical solution of the continuity equation using a finite volume method with piecewise linear reconstruction and a van Leer flux limiter. The tests demonstrate that the above formalism is correct, and provides much faster evaluation than is possible with the finite volume method while avoiding the effects of numerical diffusion.

Appendix C The angle-averaged Synchrotron function $\bar{F}(y)$ and it’s evaluation

For an isotropic distribution of magnetic field orientation, the distribution of angle $p_{\Theta}=\frac{1}{2}\sin\Theta$. We can then write the inner integral

[TABLE]

where $F(x)$ is defined in equation 23. Let

[TABLE]

Then we can write

[TABLE]

where $K_{5/3}(z)$ is the modified Bessel function of order 5/3. Crusius & Schlickeiser (1986) obtain an expression of $\bar{F}(y)$ in terms of Whittaker functions:

[TABLE]

where $W_{\lambda,\mu}(y)$ denotes Whittaker’s function (Abramowitz & Stegun, 1970, pg. 505).

Enßlin et al. (1999) give an approximation to the function $\bar{F}(y)$ to within a few percent

[TABLE]

However, such an approximation is not necessary when spline interpolation can provide an efficient and accurate computational solution. For small values of $y$ ($y\lesssim 0.01$), an asymptotic form is available (G. Bicknell, private communication):

[TABLE]

where

• •

$c_{1}=1.808418021$

• •

$c_{2}=-1.813799364$

• •

$c_{3}=0.8476959474$

• •

$c_{4}=-0.510131$