Testing the isotropy of the log N - log S slope for the NVSS radio catalogue

TL;DR

This study tests the isotropy of the log N - log S slope in the NVSS radio catalog by analyzing spectral index variations across the sky, finding only weak, statistically insignificant deviations from isotropy.

Contribution

The paper introduces a rigorous likelihood maximization algorithm to assess isotropy of the spectral index in radio source counts, providing a novel method for such analyses.

Findings

Maximum difference in spectral index is 1.3% of the full sky value.

Deviation from isotropy is only at 2σ significance, indicating weak departure.

Dipole fit results are consistent with isotropy.

Abstract

The cumulative number count $N$, of sources above a threshold is known to approximately follow a power law behaviour $N\propto S^{-x}$. We study the variation of spectral index $x$ across the sky in order to look for possible signals of violation of isotropy. We develop a rigorous algorithm of likelihood maximisation to accurately fit for the spectral index. We divide the sky into upper and lower hemispheres for a particular choice of $z$-axis and determine the difference $\Delta x$ between the best fit values of the spectral indices between the two hemispheres. The maximum value of this difference obtained by varying over the $z$-axis provides us with a measure of departure from isotropy. We find that the data support isotropy of the spectral index. The maximum difference is found to be 1.3\% of the full sky best fit value of $x$. The deviation is found to be significant only at 2$\sigma$ level which indicates a weak departure from isotropy. We also perform a dipole fit to the spectral index as a function of the angular coordinates. The result is found to be consistent with isotropy.

Tables

Table 1. Table 1: The results for full sky best fit value of x, average of the best fit value x in each hemisphere, the maximum difference in the value of x in the upper and lower hemisphere and the direction of the upper hemisphere corresponding to the maximum difference. As explained in text, catalogue A only has galactic mask and flux bounds, B has additional 22 sites removed and C has another additional local sources removed. All the errors quoted are 1 σ 𝜎 \sigma errors obtained from the log-likelihood fit. For the hemisphere average values the quoted errors are the standard deviations of the hemisphere best fits.

Data	Full sky $x$	Hemisphere Average $x$	$\text{Maximum }\Delta x$	(l,b)	p-value
A	$0.9939\pm 0.0025$	$0.9938\pm 0.0022$	$0.0134\pm 0.0050$	$({122}^{\circ },{36}^{\circ })$	2.76%
B	$0.9946\pm 0.0025$	$0.9944\pm 0.0023$	$0.0131\pm 0.0050$	$({122}^{\circ },{36}^{\circ })$	3.69%
C	$0.9939\pm 0.0025$	$0.9937\pm 0.0024$	$0.0135\pm 0.0050$	$({122}^{\circ },{36}^{\circ })$	2.76%

Equations

$$N_{>S}\propto S^{-x}.$$

$$n(\theta ,\varphi ,S)=\frac{d^{2}N}{d\Omega dS}=k{S}^{-1-x}.$$

$${\int }_{{\Omega }_r}{\int }_{S_{\text{min}}}^{S_{\text{max}}}n(\Omega ,S)dSd\Omega =\frac{{\Omega }_{r}k}{x}\left[S_{\text{min}}^{-x}-S_{\text{max}}^{-x}\right]=N_r,$$

$$P(S|x)=\frac{1}{N_r}{\int }_{{\Omega }_r}n(\Omega ,S)d\Omega =\frac{x}{\left[S_{\text{min}}^{-x}-S_{\text{max}}^{-x}\right]}{S}^{-1-x}=f(S).$$

$$\ln \mathcal{L}(x)$$

$$=N_r\ln \frac{x}{\left[S_{\text{min}}^{-x}-S_{\text{max}}^{-x}\right]}+(-1-x)\sum _i^{N_r}\ln S_i$$

$$F_i=S_i+e_i,$$

$$P(e_i|{\sigma }_i)=\frac{1}{{\sigma }_i\sqrt{2\pi }}\exp \left(-\frac{e_i^2}{2{\sigma }_i^2}\right)=g(F_i-S_i).$$

$$P(F_i|x,{\sigma }_i)$$

$$=\int \frac{x}{\left[S_{\text{min}}^{-x}-S_{\text{max}}^{-x}\right]}{S}_i^{-1-x}\frac{1}{{\sigma }_i\sqrt{2\pi }}\exp \left(-\frac{(F_i-S_i{)}^2}{2{\sigma }_i^2}\right)d{S}_i$$

$$\ln \mathcal{L}(x)$$

$$=N_r\ln \left(\frac{x}{\sqrt{2\pi }\left[S_{\text{min}}^{-x}-S_{\text{max}}^{-x}\right]}\right)+\sum _i^{N_r}\left[\ln \left\{\int S_i^{-1-x}{e}^{-\frac{(F_i-S_i{)}^2}{2{\sigma }_i^2}}d{S}_i\right\}-\ln {\sigma }_i\right]$$

$${\Delta }_N=\frac{x_{\text{up}}-x_{\text{down}}}{\sqrt{\delta x_{\text{up}}^2+\delta x_{\text{down}}^2}},$$

$$x(\hat{n})=x_0+x_1\hat{\lambda }\cdot \hat{n},$$

$$N_{>S}\propto S^{-x(S)}.$$

Full text

cosmology: observations — large-scale structure of universe — catalogs

††journal: ApJ

1 Introduction

The cumulative number count $N$, of sources above a threshold flux density $S$ per unit solid angle is typically assumed to have a power law behaviour. The slope $x$ of the $\log N-\log S$ curve has played an important role in astronomy and in early cosmology. The variation of the $\log N$ with $\log S$ can be related to the underlying cosmological model (Davidson, 1962; Longair & Scheuer, 1966; Petrosian, 1969). Observationally a power law behaviour is found to provide a reasonably good fit to data, although deviations from a pure power law behaviour have been observed (Condon et al., 1998; Tiwari et al., 2014; Tiwari & Jain, 2015). It is generally assumed that slope is independent of the direction and position of observation. This follows from the cosmological principle. However there currently exist several signals which suggest potential violations of statistical isotropy. These include dipole anisotropy in radio polarization angles (Jain & Ralston, 1999), radio source counts and brightness (Blake & Wall, 2002; Singal, 2011; Gibelyou & Huterer, 2012; Rubart & Schwarz, 2013; Tiwari et al., 2014), radio polarized source counts and brightness (Tiwari & Jain, 2015), alignment of Cosmic Microwave Radiation (CMB) quadrupole and octopole (de Oliveira-Costa et al., 2004), CMB hemispherical anisotropy (Eriksen et al., 2004), alignment of quasar polarizations (Hutsemékers, 1998; Jain et al., 2004) etc. Although not confirmed, these signal are very interesting and are being extensively studied in the literature. Hence it is important to test whether the slope $x$ is consistent with the principle of isotropy. In this work we study the direction dependence of this parameter.

The number counts as well as brightness of radio sources acquire a dipole distribution due to several contributions. These include our local motion relative to the cosmic frame of rest as well as the dipole arising due to local clustering of sources. The slope $x$ provides us with an independent and potentially more robust measure to test the cosmological principle. In particular it is unaffected by the local motion. Furthermore it is a probe of the flux distribution of sources rather than the number density and hence many factors which affect the dipole in number counts may not contribute to $x$ or may contribute differently. These factors include the local cluster sources as well as generalized cosmological models which may assume anisotropic and/or inhomogeneous Universe.

In §2 we discuss the cumulative and differential number count distribution functions and the definition of the likelihood functions used in our data analysis. We discuss the NVSS catalogue and its systematics and their removal in §3. In §4 we discuss our data analysis pipeline and data simulation processes. The results are presented in §5.

2 Theory

The cumulative number count of radio sources per unit solid angle, $N_{>S}$, above a flux density of S can be expressed as (Ellis & Baldwin, 1984):

[TABLE]

The slope of the $\log N-\log S$ relation is equal to the spectral index $x$. It is easy to see that the differential number count per unit flux density per unit solid angle, $n(\theta,\phi,S)$, along a direction of observation $(\theta,\phi)$, assuming isotropy, is given by a power law form

[TABLE]

Here $k$ is a normalisation constant dependent on the total number of sources, the flux range of the data and the solid angle. We are interested in testing the variation in the spectral index $x$ along different directions of the sky.

We use likelihood analysis in order to fit the value of the spectral index. If the differential number count is integrated over solid angle and flux intensity then we should get the total number count. Let us integrate equation (2) over the full sky and a flux density range of $S_{\text{min}}$ and $S_{\text{max}}$. We get,

[TABLE]

where $N_{r}$ is the number of sources in the flux range, $S_{\text{min}}$ to $S_{\text{max}}$, in a region of the sky which subtends a solid angle $\Omega_{r}$. This relation fixes the normalization constant $k$.

We define a probability $P(S|x)$ as the probability that any random radio source observed in the sky has a flux density of $S$. It is given by

[TABLE]

We define our likelihood function as $\mathcal{L}=\prod_{i}^{N_{r}}P(S_{i}|x)$. Using equation (4) we can write the log-likelihood as

[TABLE]

We maximize the log-likelihood given by equation (5) to obtain the best fit value of spectral index $x$. We note here that this process does not incorporate the experimental errors in the measurement of the flux density $\sigma_{i}$.

To account for the measurement error in flux density we write the observed flux density $F_{i}$ as

[TABLE]

where $S_{i}$ is the possible theoretical value drawn from distribution equation (4) and $e_{i}$ is drawn from a standard normal distribution with variance given by $\sigma_{i}^{2}$, i.e.,

[TABLE]

We can then define the probability of an observed value of flux density $F_{i}$ as

[TABLE]

using equations (4), (7) and (6).

The log-likelihood then is given by $\ln\mathcal{L}=\ln\left(\prod_{i}P(F_{i}|x,\sigma_{i})\right)$. Using equation (8) the form of the total log-likelihood is given by:

[TABLE]

We use log-likelihood functions defined in equation (5) and equation (9) to fit for the spectral index $x$.

3 Data

The NRAO VLA Sky Survey (NVSS) (Condon et al., 1998) is a continuum survey at 1.4 GHz covering the northern sky above a declination of $-40^{\circ}$ (J2000) with 82% sky coverage. The NVSS catalogue contains 1773484 radio sources. We remove the galactic plane by removing all sources with galactic longitude $|b|<15^{\circ}$. The galactic plane masking removes an additional 26% of the sky. Hence we have only 56% of the total sky area for analysis. The full catalogue with only the galactic portion masked is termed catalogue A for this work.

3.1 Removing local clustering

Blake & Wall (2002) had identified 22 sites of bright and extended radio galaxies. We follow their procedure of masking these 22 sites. This removes 23386 sources from the full NVSS catalogue. The total area removed is relatively small. As noted by Blake & Wall (2002), the NVSS has poor resolution, so two or more galaxies situated near each other are rarely resolved into the individual galaxies. This is important as this would also increase the flux density for the unresolved object. Having many such unresolved clustered objects can cause deviation in the measured value of the spectral index. This masked catalogue of 1750098 sources is further masked along the galactic plane, as discussed earlier, and the resulting data set is termed catalogue B.

The expectation of isotropy is valid only when looking at cosmological length scales. The NVSS catalogue contains galaxies that are in our local neighbourhood and on these short distance scales local clustering would cause departure from isotropy. Blake & Wall (2002) have shown that the local clustering effects originate from redshifts of $z<0.03$. The Infrared Astronomical Satellite Point Source Redshift (IRAS PSCz) catalogue (Saunders et al., 2000) and the 3rd Reference Catalogue of bright galaxies (RC3) (de Vaucouleurs et al., 1991; Corwin et al., 1994) give a good sample for galaxies in the local supercluster. We prepare another catalogue which removes the NVSS sources which are within 30 arcsec separation from sources listed in IRAS PSCz catalogue and the RC3. This cut removes 13597 sources. We additionally mask the 22 bright and extended radio sources as discussed above. This keeps 1736583 sources. After masking the galactic plane as previously described we prepare catalogue C.

3.2 Systematics

The NVSS catalogue is well known to suffer from systematics. For the survey (Condon et al., 1998; Chen & Schwarz, 2016) between the declinations of $-10^{\circ}$ and $78^{\circ}$ the D configuration of the VLA was used while for the remaining regions of the sky the DnC configuration was used. The D configuration apparently has more brightness sensitivity over the hybrid DnC configuration and hence the portion of the sky sampled with the D configuration has a higher source density. The number density of sources show clear variation with declination as can be seen in Figure 1.

To study the effects of configuration systematics we divide the sky in strips of constant declination with a width of $5^{\circ}$. We look at the ratio of the number density $(n)$ in each strip to the average number density $(\bar{n})$ over the whole sky. We also define a quantity called sky brightness $(S_{\text{sky}})$, which is the sum of the flux densities of all the sources in a given portion of the sky divided by the area of the patch. We look at the ratio of the local sky brightness of a strip to the whole sky average brightness $(\bar{S}_{\text{sky}})$. The variations of these two ratios with declination is shown in Figure 2. From the plot of number density ratios it is very clear that below the declination of $-10^{\circ}$ and above the declination of $78^{\circ}$ there is clear dip in the number density of sources due to the DnC configuration used in surveying as discussed above. It has been shown by Blake and Wall that the declination dependence can be greatly reduced by placing a cut on the minimum flux of the sources to be considered. They showed that over a flux density of 15 mJy the effects due to limited brightness sensitivity do not remain an issue and the number density evens out. For our analysis we choose 20 mJy as the minimum flux density. We also place a cut on the maximum flux density of 1000 mJy to remove sources with large brightness.

The ratio of the sky brightness is also shown in Figure 2. Without any cuts or masking the NVSS data show very large spikes of high sky brightness. This is largely attributable to the bright and clustered sources near the galactic plane as is clearly visible in Figure 1. The masking of the galactic plane with the flux cuts makes the variation of the sky brightness ratio fairly close to unity. For all our analysis in this work we have used the lower flux density cut of 20 mJy and upper flux density cut of 1000 mJy. Following examples in literature we understand that these measures would be sufficient to suppress the effects of systematics arising from the configuration differences of the survey.

4 Data Analysis

We find the best fit value of the spectral index $x$, by maximising the log-likelihood defined in equation (5). This form does not take into account the error in flux density measurements and it is computationally cheap as opposed to the form in relation find (9). For computing equation (5), we need to determine $N_{r}$, $S_{\text{max}}$, $S_{\text{min}}$ and $\sum\ln S_{i}$. Then the calculation of the log-likelihood can be performed which allows us to determine the spectral index $x$ in any selected region of the sky. Since we need to perform the sum $\sum\ln S_{i}$ only once, the maximisation of log-likelihood can be done with good computation speed. We simply scan over all values of $x$ in a chosen range with a very small step size, equal to $0.0001$. However for the form equation (9) we need to compute numerical integrals for each source for every step in the change of spectral index. Given the very large number of sources this process is computationally prohibitive for our requirements. Hence we use this more extensive procedure only for a full sky calculation of the spectral index. This allows us to to test the validity of our procedure based on the likelihood function defined in equation (5). We later discuss the differences in full sky fit of the spectral index obtained with the two procedures.

In our analysis, we are not considering the positional uncertainties in the NVSS catalogue. These uncertainties are not small but for our analysis the absolute position of the source is not relevant. In all our fits to the spectral index we work with a rather large area containing large number of sources. Hence if the source positions have errors in measurement then it would affect sources close to the edges of the area and the net effects should be small for our procedure. Hence we neglect the effects of positional errors for the present analysis.

4.1 Simulated data and method validation

The method of finding the best fit value of spectral index is as described above. Since we have masked the data and imposed flux cuts we would like to test our likelihood maximisation technique against partial sky data and flux cuts. To ensure controlled testing of our analysis method we generate simulated data.

To generate a simulated source catalogue, we used a Mathematica script. We use the probability distribution given by equation (4) with a value of spectral index set from full sky fits, $S_{\text{max}}$ being a randomly chosen upper flux limit from the interval [1200, 1500] mJy and $S_{\text{min}}$ the randomly chosen lower flux limit from the interval [5, 15] mJy. A random number of source positions are generated which are between 1.9 million and 2.1 million and uniformly distributed over the full sky. The random positions are assigned a random flux density drawn from the defined probability distribution. This generates a catalogue of sources uniformly distributed over the sky and obeying the distribution in equation (2).

The catalogue simulated as above is in celestial coordinates. We then proceed to remove all sources south of declination $-40^{\circ}$. This produces a catalogue with 82% sky coverage as the NVSS. We further impose the galactic plane masking and impose flux cuts of 20 mJy and 1000 mJy as discussed previously. We do not mask the 22 sites because such contaminations are not present and the area of these maskings are much smaller. We have followed the most important steps of our data processing pipeline to produce a masked partial sky data within a limited flux range. However, we do accept that this is an idealised simulation as we do not take into account the instrumental limitations of the VLA and the configurations used during the survey. Including such details would require us to incorporate much more information regarding the survey making the simulation computationally prohibitive for the purpose of this work. Testing against these simulations would help us find artefacts that might arise out of the data analysis pipeline. We essentially test the effects of flux density bounds and partial sky data.

We simulated a catalogue of 1971652 objects uniformly distributed over the sky with a minimum flux of 9.87 mJy and maximum of 1471.96 mJy from a flux density distribution with an arbitrarily chosen value of the spectral index of 0.5. For the full sky simulated catalogue and no flux cuts we find the best fit spectral index to be 0.5009. When we mask 44% of the sky to get identical sky coverage as the actual data sets prepared but without any flux cuts, we get spectral index of 0.5015. When analysing the data without any masking but 20-1000 mJy flux density bounds we obtain $x=0.4994$. Working with both the masking and flux density bounds leads to $x=0.4999$. So the deviation we find in the best fit value is very small, less than the error in the spectral index which is found to about $0.2\%$. So our method of fitting is very robust to masking of the sky and the flux density bounds. We do not expect any bias in our results due to our data analysis pipeline.

4.2 Isotropy testing procedure

We divide the sky into a HealPix (Górski et al., 2005) grid of $N_{\text{side}}=32$. We vary the z axis of our system along the centre of each of the pixels. With z axis aligned along this direction we divide the sky in two hemispheres, one centred about $\hat{z}$ and the other along $-\hat{z}$. Then we find the value of the spectral indices $x_{\text{up}}$ and $x_{\text{down}}$ for the upper and lower hemispheres for z axis centred about each of the pixels of the HealPix grid. We define the quantity $\Delta x=x_{\text{up}}-x_{\text{down}}$. This gives us the difference in values of the spectral index in the two halves of sky. We then find the pixel (and its corresponding direction) for which $\Delta x$ is maximum. We have also checked if the full sky value of the spectral index matches to the average of all the half sky spectral indices. This is meant for checking consistency.

The number of sources and hence the error in the best fit value of $x$ depend on the choice of hemisphere. In order to take this into account we define an error normalised quantity $\Delta_{N}$ as:

[TABLE]

where $\delta x_{\text{up}}$ is the error in $x_{\text{up}}$ and $\delta x_{\text{down}}$, the error in $x_{\text{down}}$. The denominator is the two errors added in quadrature. Following the above process of finding the direction of maximum $\Delta x$, we find the direction along which $\Delta_{N}$ maximises.

We also perform a fit to a dipole model for departure from isotropy. We write the spectral index $x$ for a pixel along a direction $\hat{n}$ as:

[TABLE]

where $\hat{\lambda}$ is the direction of maximum anisotropy. We fit for the four parameters of this model: one each for $x_{0}$ and $x_{1}$, and two for $\hat{\lambda}$. The fits are performed on catalogue C. As stated above, in this case our observable is the value of $x(\hat{n})$ in the pixel along the direction $\hat{n}$ using the HealPix pixelization scheme. In this case we use the HealPix map with $N_{\text{side}}=2$ in order have sufficiently large number of sources in any particular pixel. We avoided pixels with less than 4000 sources. We finally ended up with 32 pixels, out of a total of 48, with independent fits of the spectral index. The data set prepared in this way is then fitted to equation (11) by minimising $\chi^{2}$. The goodness of fit is determined by comparison with a null case of fitting to an isotropic full sky spectral index.

5 Results

The results corresponding to maximum difference $\Delta x$ are contained in Table 1. The full sky value of the spectral index was fitted over the entire sky for each of the catalogues A, B and C. We searched for the pixel, for which the upper and lower hemispheres show maximum differences in the values of spectral indices, $\Delta x$. These maximum differences and the corresponding directions are given in Table 1. The means of the spectral indices fitted over all the hemispheres are also given. The plots for the full sky log-likelihood are shown for the three catalogues in first panel of Figure 4. The histogram plots for the best fit spectral indices over the hemispheres over the entire HealPix grid are shown in the second panel of Figure 4. For consistency these values should be very close to the full sky average. The difference between the hemisphere average values and the corresponding full sky values is found to be very small. The quoted error in the values of the full sky fit of $x$ is $1\sigma$ error from the log-likelihood. The uncertainty in best fit values of $x$ is about 0.25%. The standard deviation of the half sky best fit values of $x$ with the z axis along all possible positions of the HealPix grid is quoted as the error is hemisphere average value of $x$.

We also see that the maximum $\Delta x$ is barely a 1.3% variation from the full sky value. Hence, at best, the deviation in spectral index is a tiny effect. The error quoted in the value of $\Delta x$ is the maximum uncertainties in the values of $x_{\text{up}}$ and $x_{\text{down}}$, added in quadratures, for the direction along which $\Delta x$ is maximum. The uncertainty in $\Delta x$ is apparently very large, $\sim 37\%$. The variation of $\Delta x$ with direction is shown in Figure 3. The fit for the case of error normalised difference $\Delta_{N}$ also maximises along the same direction $(l=122^{\circ},b=36^{\circ})$ as that of the maximum $\Delta x$, with a maximum $\Delta_{N}=2.69$ for catalogue C.

To test the results we have simulated 450 catalogues as described previously with the spectral index of simulation set to 0.9939. As explained before all these simulated maps are processed by the entire data processing pipeline. The histogram plot for the values of maximum $\Delta x$ from these simulations is shown in the left graph of Figure 5. We have also shown the observed values of maximum $\Delta x$ for all the catalogues in the figure. The figure shows that all these value lie towards the tail of the histogram and have a relatively small probability. We calculate the p-value as the percentage of simulations for which the value of $\Delta x$ exceeds or equals the observed maximum value. The p-values for each of the results are also shown in Table 1. Hence we find that the deviation from isotropy is significant at 2 $\sigma$ level. This indicates a weak signal of anisotropy which may be probed in future observations. Due to the relatively weak signal we conclude that the current data are consistent with isotropy. The value of $|\Delta x|$ is $0.0025\pm 0.0049$ along the direction of the NVSS dipole (Tiwari et al., 2014). Along the direction of the CMB dipole modulation (Eriksen et al., 2004) $|\Delta x|$ is $0.0023\pm 0.0050$.

We next describe our results for the dipole fit based on equation 11. In this case we determined the spectral index in each pixel using $N_{\text{side}}=2$. Each pixel was found to contain anywhere between 5000 and 9000 sources. When the data were fitted with the isotropic null case we get $x_{0}=0.9929\pm 0.0036$ with a $\chi^{2}=27.5$. When fitted to our model of dipole anisotropy it gives: $x_{0}=0.9910\pm 0.0039$, $x_{1}=0.0105\pm 0.0081$ along ($l=137^{\circ}\pm 37^{\circ},b=17\pm 29^{\circ}$) with $\chi^{2}=24.1$. Hence we do not find a significant signal of dipole anisotropy in the data.

We have thus far used the log-likelihood function defined in equation (5), which does not incorporate the uncertainties in the measurement of the flux density. We have discussed before that the likelihood function defined in equation (9) is more complete because it incorporates the error in flux density measurements. We have computed a value for full sky best fit spectral index by using equation (9). The full sky best fit spectral index is found to be $0.9984\pm 0.0022$ when including the error in flux density in the analysis. This is 0.45% different from the full sky best fit values computed using equation (5). This difference is not significant enough to require the use of computation intensive method of maximising the likelihood function defined in equation (9) for the purpose of this work.

In our paper we have used a power law relationship given in equation 1. Although this is the form most widely used in the literature it is known that a generalization, i.e.,

[TABLE]

with $x(S)=a+b\log S$ provides a much better fit (Condon et al., 1998; Tiwari et al., 2014; Tiwari & Jain, 2015). This is also clear from the second plot of Figure 5 where we have shown the spectral index plotted as a function of the lower flux density bound. We relaxed the upper flux density bound for obtaining this plot. It clearly shows a logarithmic dependence of $x$ on $S$. We obtained a fit for $x(S)=a+b\log S$ with $a=0.2759\pm 0.0124$ and $b=0.2438\pm 0.0019$. It will clearly be of interest to extend our calculation in order to determine the direction dependence of both the parameters $a$ and $b$ rather than restricting ourselves to a single parameter, i.e., the spectral index $x$. We postpone this to future research.

6 Conclusions

In this work we have thoroughly tested the angular dependence of the spectral index $x$ with various flux density cuts and masking. We have searched for deviation of $x$ by dividing the sky in different hemispheres. Our technique of estimating the spectral index of a power law distribution was found to be robust with very small uncertainty in the fitted value. We do not find a significant deviation from isotropy. The deviation is at most a 2$\sigma$ effect which is mildly interesting and should be probed in future observations. We conclude that if there is any departure in the isotropic behaviour of $\log N-\log S$ relation for the distant radio sources in the NVSS catalogue, it is indeed a tiny effect.