PROBING BULK FLOW WITH NEARBY SNe Ia DATA

The use of Type Ia supernovae (SNe) as a tool to probe cosmology has a long and distinguished history, starting with early pioneering works (Riess et al. 1998, 1999; Perlmutter et al. 1999; Garnavich et al. 1998) and continuing to the present (Blakeslee et al. 2003; Tonry et al. 2003; Riess et al. 2004; Barris et al. 2004; Matheson et al. 2005; Hicken et al. 2009, 2009; Hamuy et al. 2006; Folatelli et al. 2010; Astier et al. 2006; Sako et al. 2008; Miknaitis et al. 2007; Wood-Vasey et al. 2007; Holtzman et al. 2008; Kessler et al. 2009; Freedman et al. 2009). Various existing catalogs (Tonry et al. 2003; Hicken et al. 2009; Folatelli et al. 2010; Kessler et al. 2009; Ganeshalingam et al. 2013; Amanullah et al. 2010), containing hundreds of the standardizeable candles, have been used to infer the existence of a late-time accelerating epoch with strong statistical significance (Wood-Vasey et al. 2007; Suzuki et al. 2012; Sullivan et al. 2011). Their complementarity to other cosmological data sets (Tegmark et al. 1998) allows one to set stringent constraints on the now widely accepted ΛCDM model.

To maximize the constraining power of SNe, one must use a combination of both high- and low-redshift data. The low-redshift SNe play an important role, acting essentially as an anchor in the redshift–luminosity distance relation (Hamuy et al. 2006). They are relatively insensitive to cosmological parameters; however, their worth lies in allowing us to determine the relative brightness of the SN. Unfortunately, existing low-redshift samples are sensitive to comparatively large systematic uncertainties, in particular calibration issues and host mass–SN brightness correlations. Future surveys are expected to provide a more homogeneous set of accurately calibrated objects. When using current data, however, it is important that we understand and account for the various measurement and astrophysical uncertainties.

One phenomenon that will affect the low-redshift data is the existence of peculiar velocities. It is well known that the local group has a velocity V_bulk = 627 ± 22 km s⁻¹ relative to the cosmic microwave background rest frame, in the direction (b, ℓ) ∼ (30°, 276°) (Kogut et al. 1993). The origin of this motion is typically attributed to low-redshift superclusters, whose gravitational attraction will produce coherent velocities in nearby galaxies. However, recent claims of a detection of a coherent bulk flow out to very large distances d_c ∼ 800 Mpc (Kashlinsky et al. 2009a, 2009b, 2010, 2011) have led to more exotic explanations such as an inhomogeneous pre-inflationary spacetime (Turner 1991) or higher-dimensional gravity theories (Afshordi et al. 2009). Large-scale cosmological anisotropies can also be generated by magnetic fields produced by Lorentz violating effects during inflation (Campanelli 2009).

Assuming that the bulk flow is due to local attractors, the lack of all-sky and deep peculiar velocity data means that one cannot yet definitively answer the question as to what mass distributions are responsible for the observed flow. Nor can we deduce at what scale the coherent motion ceases to be significant. The issue has been considered for decades, and known superclusters such as Shapley and the so-called Great Attractor are commonly thought to contribute at least partially to the observed bulk motion (Lynden-Bell et al. 1998; Zamorani 1989). More recently, different groups have used various data sets (Watkins et al. 2009; Turnbull et al. 2012; Colin et al. 2011; Nusser et al. 2011; Nusser & Davis 2011; Branchini et al. 2012; Ma & Scott 2013) to estimate the direction and magnitude of the flow at different scales. While consensus appears to have been reached regarding the direction, the magnitude and consistency with ΛCDM remains an open question. It is an issue that has drawn scrutiny from various sources and remains contested (Kashlinsky et al. 2009a; Keisler 2009; Osborne et al. 2011; Mody & Hajian 2012; Lavaux et al. 2013). Some groups have suggested that SNe at redshift z ≳ 0.05 could exhibit evidence of back-infall into the Shapley supercluster, indicating that it might be the source of the observed bulk motion (Colin et al. 2011).

SNe are ideal candidates to study the local bulk flow, being bright distance indicators with precise redshift measurements. The depth of SN surveys is typically large, allowing us to probe the scale at which the bulk flow remains coherent. However, this gain in depth relative to local galaxy catalogs is offset by the comparatively modest quantity of existing data. Future surveys such as the Large Synaptic Survey Telescope (Abate 2012) are expected to yield thousands of SNe and should definitively answer the question as to the magnitude and direction of the bulk flow. The increasing number of detected SNe will allow us to move beyond simply treating the peculiar velocities as systematics in distance measurements. The SNe measurements yield information regarding the peculiar velocity field and open up a new avenue in exploring this dynamical cosmological probe (Johnson et al. 2014).

The question that we address in this work is whether one can observe a statistically significant deviation from isotropy in existing low-redshift SN data, and if so whether this anisotropic signal is consistent with the standard cosmological model. To achieve this aim, we adopt a nonparametric method, which compares the distance modulus residuals with mock realizations drawn from the underlying covariance matrix of the data (Colin et al. 2011). In this way we address the significance of any anisotropic signal within the context of the systematic and statistical uncertainties associated with the data. We make no assumptions regarding the nature of any potential anisotropy. This is an advantage in the sense that we do not need to specify a model and can detect any source of anisotropic behavior. However, a consequence of this model independence is that one cannot measure the magnitude of the bulk flow using our approach.

The paper will proceed as follows. In Section 2 we briefly review the method used to test the isotropic hypothesis. We discuss the catalogs used and how we combine them in Section 3, and we present our results in Section 4. We discuss the effect of inhomogeneous data distributions and the effect of correlations and conclude in Section 5.

The method of residuals has been discussed in detail in Colin et al. (2011), Feindt et al. (2013), and Appleby & Shafieloo (2014), and we begin with a brief review. We note that numerous other approaches to testing isotropy exist in the literature, and we direct the reader to Kolatt & Lahav (2001), Bonvin et al. (2006), Gordon et al. (2007), Schwarz & Weinhorst (2007), Gupta et al. (2008), Appleby & Linder (2013), Cooray et al. (2010), Gupta & Saini (2010), Cooke & Lynden-Bell (2010), Koivisto & Mota (2008), Koivisto et al. (2011), Campanelli et al. (2006, 2007), Appleby et al. (2010), Antoniou & Perivolaropoulos (2010), Blomqvist et al. (2008, 2010), Tsagas (2010), and Appleby & Shafieloo (2014) for searches of (predominantly cosmological) anisotropic signals. The effect of peculiar velocities on SN data has been discussed in Hui & Greene (2006) and Davis et al. (2011).

Our approach entails creating a map on the sky of the sum of distance modulus residuals, smoothed over a particular scale δ. The first step is to calculate a global best-fit cosmology. We fix the background expansion to be flat ΛCDM and fit the parameter Ω_m0 to our chosen SN catalog by minimizing the χ² distribution

$$\begin{eqnarray}&&{{\chi }^{2}}=\delta {{{\boldsymbol{ \mu }} }^{{\rm T}}}{{{\Sigma }}^{-1}}\delta {\boldsymbol{ \mu }} \end{eqnarray} \tag{ 1 }$$

where $\delta {\boldsymbol{ \mu }}$ is the data vector of distance modulus residuals, δμ_i = μ_i(z_i)−μ_th(z_i), and Σ is the covariance matrix of the data. Parameter μ_i(z_i) is the distance modulus of the ith data point, and μ_th(z_i) is the theoretically predicted ΛCDM value at data point redshift z_i. Our approach to calculating the best-fit cosmology differs from more detailed cosmological parameter estimation considered in Amanullah et al. (2010),  Suzuki et al. (2012), and Kowalski et al. (2008), for example. In the Union 2.1 catalog the intrinsic dispersion σ_i is kept as a free parameter to account for the intrinsic scatter of the SN magnitudes, and the cosmology fit such that χ² = 1 (per degree of freedom). However, the best-fit cosmological parameters that we obtain are essentially identical to the values quoted in the various catalogs that we will use. Furthermore, the cosmological dependence of our result is expected to be very weak, as we are focusing solely on the low-redshift z < 0.1 SN sample.

We denote the best-fit distance modulus as $\tilde{\mu }(z,{{H}_{0}},{{{\Omega }}_{{\rm m}0}})$. Using this, we construct the error-normalized residuals of the data from the model (Perivolaropoulos & Shafieloo 2009),

$$\begin{eqnarray}&&{{q}_{i}}({{z}_{i}},{{\theta }_{i}},{{\phi }_{i}})=\frac{{{\mu }_{i}}({{z}_{i}},{{\theta }_{i}},{{\phi }_{i}})-{{\tilde{\mu }}_{i}}({{z}_{i}})}{{{\sigma }_{i}}({{z}_{i}})},\end{eqnarray} \tag{ 2 }$$

Here (θ_i, ϕ_i) are the angular positions of the ith data point on the sphere and σ_i(z_i) is the error associated with the diagonal component of the covariance matrix.

Next, we define a measure Q(θ, ϕ) on the surface of the sphere using these residuals,

$$\begin{eqnarray}&&Q(\theta ,\phi )=\mathop{\sum }\limits_{i=1}^{{{N}_{{\rm SN}}}}{{q}_{i}}({{z}_{i}},{{\theta }_{i}},{{\phi }_{i}})W(\theta ,\phi ,{{\theta }_{i}},{{\phi }_{i}}),\end{eqnarray} \tag{ 3 }$$

where N_SN is the number of SN Ia data points and W(θ, ϕ, θ_i, ϕ_i) is a weight (or window) function that represents a two-dimensional smoothing. We define the weight using a Gaussian distribution

$$\begin{eqnarray}&&W(\theta ,\phi ,{{\theta }_{i}},{{\phi }_{i}})=\frac{1}{\sqrt{2\pi }\delta }{\rm exp} \left[ -\frac{L{{(\theta ,\phi ,{{\theta }_{i}},{{\phi }_{i}})}^{2}}}{2{{\delta }^{2}}} \right],\end{eqnarray} \tag{ 4 }$$

where δ is the width of smoothing and L(θ, ϕ, θ_i, ϕ_i) is the distance on the surface of a sphere of unit radius between two points with spherical coordinates (θ, ϕ) and (θ_i, ϕ_i),

$$\begin{eqnarray}&&L(\theta ,\phi ,{{\theta }_{i}},{{\phi }_{i}})=2\;{\rm arcsin} \frac{R}{2},\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray*}\begin{array}{rcl} R & = & \left( {{\left[ {\rm cos} ({{\theta }_{i}}){\rm cos} ({{\phi }_{i}})-{\rm cos} (\theta ){\rm cos} (\phi ) \right]}^{2}} \right. \\ {} & {} & +{{\left[ {\rm cos} ({{\theta }_{i}}){\rm sin} ({{\phi }_{i}})-{\rm cos} (\theta ){\rm sin} (\phi ) \right]}^{2}} \\ {} & {} & +{{\left. {{\left[ {\rm sin} ({{\theta }_{i}})-{\rm sin} (\theta ) \right]}^{2}} \right)}^{1/2}}. \\ \end{array}\end{eqnarray*}$$

Any anisotropy in the data will translate to Q(θ, ϕ) deviating from zero. Finally, we adopt a value for δ, calculate Q(θ, ϕ) on the whole surface of the sphere, and find the extreme values of this function. In this work we utilize the maximum dipole,

$$\begin{eqnarray}&&{\Delta }{{Q}_{{\rm d}}}=Q(\theta ,\phi )-Q({{\theta }_{{\rm d}}},{{\phi }_{{\rm d}}}),\end{eqnarray} \tag{ 6 }$$

where (θ_d, ϕ_d) are the angles that locate the dipole of (θ, ϕ) on the unit sphere. In what follows we choose the galactic coordinate system (θ, ϕ) = (b, ℓ). In this case, we have (θ_d, ϕ_d) and (−b, ℓ + 180°) if ℓ < 180° and (−b, ℓ−180°) for ℓ > 180°. In principle, one could use any number of different measures of anisotropy—for example, the maximum and minimum of this function on the sky, Q_max and Q_min. However, for practical purposes current SN data are neither sufficiently copious nor sufficiently precise to detect anything other than large-scale flows such as the dipole—hence our use of ΔQ_d.

A large value of ΔQ_d implies a significant anisotropy in the data. However, we expect that this quantity, which denotes the extreme of the Q(θ, ϕ) function on the sphere, will not be zero even if the Hubble residuals are isotropic on the sky. We must test the significance of the magnitude of ΔQ_d by creating mock realizations.

In this work we construct two sets of N_real = 1000 realizations of our SN catalogs. We describe both in turn below, and in what follows we will refer to them as realizations A and B, respectively.

In the first case we simply follow the original method of Colin et al. (2011) and construct our realizations by holding each data point at fixed redshift and position on the sky, and we estimate the distance modulus as the isotropic best-fit value $\tilde{\mu }({{z}_{{\rm i}}},{{{\Omega }}_{{\rm m}0}})$ plus a contribution drawn from a Gaussian of width σ_i, where σ_i is the observational uncertainty on the ith data point quoted in the catalog. Hence, we have ${{\mu }_{{\rm i}}}=\tilde{\mu }({{z}_{{\rm i}}},{{{\Omega }}_{{\rm m}0}})+\delta {{\mu }_{{\rm i},\;{\rm G}}}$, where δμ_{i, G} is the Gaussian component. For each data realization, ΔQ_d is obtained and an empirical probability distribution function (PDF) is constructed for this function. We then ascertain the significance of ΔQ_d obtained from the data by calculating its p-value from the PDF. When comparing the data to this set of realizations, we are addressing the question as to whether any potential anisotropy detected in the data is consistent with observational uncertainties within the data, which we are taking to be Gaussian and uncorrelated. Note that we are implicitly assuming that the error is Gaussian here—for this reason our choice of observable is the distance modulus rather than the velocity (Strauss & Willick 1995). The p-value in this case signifies the frequency of occurrence of a ΔQ_d value obtained from the simulations that is greater than or equal to the data value of ΔQ_d. The simulations are drawn from an isotropic ΛCDM model with no anisotropic signal, so one can interpret the p-value as the probability that the data value of ΔQ_d can be obtained in an isotropic universe.

Our primary aim is to test the consistency of the data with the underlying cosmological model. Inhomogeneous structure present in the late-time universe generates coherent velocity flows, which breaks isotropy at low redshift. One can estimate the expected magnitude of the velocity field on large scales using linear perturbation theory. For the second set of simulated realizations, we use the covariance matrix constructed in Johnson et al. (2014), which accounts for both observational uncertainties and the presence of coherent large-scale peculiar velocities. The catalog is discussed further in Section 3; here we briefly review the contributions to the covariance matrix.

The covariance between the ith and jth data point, C_ij, is given by (Johnson et al. 2014; Hui & Greene 2006)

$$\begin{eqnarray}\begin{array}{rcl} {{C}_{ij}} & = & {{\left( \frac{5}{{\rm ln} 10} \right)}^{2}}\left( 1-\frac{{{(1+{{z}_{i}})}^{2}}}{H({{z}_{i}}){{d}_{{\rm L}}}({{z}_{i}})} \right)\left( 1-\frac{{{(1+{{z}_{j}})}^{2}}}{H({{z}_{j}}){{d}_{{\rm L}}}({{z}_{j}})} \right) \\ {} & {} & \times \int \frac{dk}{2{{\pi }^{2}}}{{k}^{2}}{{P}_{{\rm vv}}}(k,a=1)W(k,{{\alpha }_{ij}},{{r}_{i}},{{r}_{j}}) \\ \end{array}\end{eqnarray} \tag{ 7 }$$

where z_i, z_j are the redshifts of the data points and H(z_{i, j}), d_L(z_i,j) the corresponding Hubble parameter and luminosity distance at said redshifts. P_vv(k, a = 1) is the linear velocity power spectrum obtained for our best-fit cosmological parameters, evaluated at z = 0. The kernel W(k, α_ij, r_i, r_j) is the angular component of the integral over the Fourier modes and is given by

$$\begin{eqnarray}\begin{array}{rcl} W(k,{{\alpha }_{ij}},{{r}_{i}},{{r}_{j}}) & = & \frac{1}{3}({{j}_{0}}(k{{A}_{ij}})-2{{j}_{2}}(k{{A}_{ij}}){{{\hat{r}}}_{i}}.{{{\hat{r}}}_{j}} \\ {} & {} & +\frac{1}{A_{ij}^{2}}{{j}_{2}}(k{{A}_{ij}}){{r}_{i}}{{r}_{j}}{{{\rm sin} }^{2}}({{\alpha }_{ij}}) \\ \end{array}\end{eqnarray} \tag{ 8 }$$

where ${{\alpha }_{ij}}={{{\rm cos} }^{-1}}({{\hat{r}}_{i}}.{{\hat{r}}_{j}})$, ${{A}_{ij}}=|{{r}_{i}}-{{r}_{j}}|$, r_i is the position vector of the ith data point, and j_0,2 are Bessel functions. See (for example) Appendix A of Ma et al. (2011) for a derivation of this function.

The first two redshift-dependent cofactors in Equation (7) transform from correlations between peculiar velocities to correlations between distance modulus fluctuations (Hui & Greene 2006). The integral over the Fourier modes is related to the correlation between data points i and j owing to the fact that they trace the same underlying density field—the trace components of C_ij are cosmic variance and the off-diagonal terms the cross-correlation. The velocity power spectrum P_vv is obtained from linear perturbation theory assuming ΛCDM. Bulk velocities on the scales that are being considered in this work are well modeled using linear theory, and so this contribution to the covariance matrix accounts for large-scale velocity components expected within the standard cosmological model.

To this covariance matrix, the diagonal component of the catalog observational uncertainty is added according to

$$\begin{eqnarray}&&{{{\Sigma }}_{ij}}={{C}_{ij}}+\sigma _{i}^{2}{{\delta }_{ij}}\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&\sigma _{i}^{2}=\sigma _{{\rm obs}}^{2}+{{\left( \frac{5}{{\rm ln} 10} \right)}^{2}}{{\left( 1-\frac{{{(1+{{z}_{i}})}^{2}}}{H({{z}_{i}}){{d}_{{\rm L}}}({{z}_{i}})} \right)}^{2}}\sigma _{{\rm v}}^{2},\end{eqnarray} \tag{ 10 }$$

σ_obs is the observational error quoted in the catalog, and σ_v is the uncertainty associated with nonlinear peculiar velocities of the galaxies. In Johnson et al. (2014) σ_v was treated as a free parameter and fit according to the data. Here we fix its value to σ_v = 400 km s⁻¹, in accordance with Johnson et al. (2014).

To obtain the N_real = 1000 mock realizations, we take the full covariance matrix (9) and diagonalize it by calculating its eigenvectors,

$$\begin{eqnarray}&&{{{\bf \Phi }}^{{\rm T}}}{\bf \Sigma \Phi }={\bf \Lambda }\end{eqnarray} \tag{ 11 }$$

where ${\bf \Phi }$ is the N_SN × N_SN matrix constructed from the eigenvectors of the full covariance matrix ${\bf \Sigma }$ and ${\bf \Lambda }$ is the diagonal matrix constructed from the eigenvalues of ${\bf \Sigma }$. If ${\bf \Sigma }$ is the covariance matrix of the distance modulus residual vector $\delta \underline{\mu }$, then ${\bf \Lambda }$ is the covariance matrix of the decorrelated data vector $\delta {{\underline{\mu }}^{({\rm diag})}}={{{\bf \Phi }}^{{\rm T}}}\delta \underline{\mu }$. To incorporate the correlations into our simulations, we obtain our N_real = 1000 realizations by drawing Gaussian realizations from the diagonal covariance matrix ${\bf \Lambda }$ to obtain a mock data vector $\delta {{\underline{\mu }}^{({\rm diag})}}$ and then transforming back into the original basis according to $\delta \underline{\mu }={\bf \Phi }\delta {{\underline{\mu }}^{({\rm diag})}}$. These δμ residuals are used to construct a PDF of ΔQ_d in our second set of realizations. Once again, the p-value is defined as the fraction of ΔQ_d values from the simulations that are greater than or equal to the data value. The simulations now contain the effect of large-scale velocity perturbations on the data and hence account for correlations between data points owing to the fact that they are tracing the same underlying density and velocity fields.

As a final comment in this section, one should be careful when attempting to infer the direction of maximal anisotropy using this method. It is known (Appleby & Shafieloo 2014) that the method is capable of approximately reconstructing the direction of an anisotropic signal for a homogeneous distribution of data, but the result can be skewed when the data are inhomogeneously scattered on the sky (Feindt et al. 2013). We test the reliability of the method in selecting the "true" anisotropic direction (b, ℓ) in Section 4.

We apply the method outlined in Section 2 to a number of existing data sets in the literature, specifically the Union 2.1 (Amanullah et al. 2010; Suzuki et al. 2012; Kowalski et al. 2008), Constitution (Hicken et al. 2009), and LOSS (Ganeshalingam et al. 2013) samples. Our aim in comparing these different catalogs is to test for consistency among the different light-curve fitting procedures, to ensure that any anisotropic signal is not due to unknown systematics associated with modeling the data. However, the main purpose of this work is to test the consistency of the SN data with the ΛCDM model—for this purpose we utilize a recently developed catalog that combines the majority of low-redshift SNe from a variety of surveys (Johnson et al. 2014).

The Union 2.1 sample contains N_SN = 175 SNe out to redshift z < 0.1, which is the limit to which we perform our analysis. Beyond this, any effect due to bulk flow velocities will be practically negligible. The data set is an amalgamation from numerous sources—and we direct the reader to Amanullah et al. (2010), Suzuki et al. (2012), and Kowalski et al. (2008) for details. As the SN data arise from a number of different telescopes, each with their own systematics and calibrations, each survey's contribution to the overall set is carefully modified to include effects such as photometric zero-point offsets, contamination, Malmquist bias, K-corrections, and gravitational lensing. The light curves of the entire sample are analyzed using the SALT2 light-curve fitting procedure, in which each SN is assigned three parameters: the peak magnitude $m_{{\rm B}}^{{\rm max} }$, the light-curve width x₁, and c, which encodes the effects due to intrinsic color and dust reddening. The distance modulus is then constructed as

$$\begin{eqnarray}\begin{array}{rcl} {{\mu }_{{\rm B}}} & = & m_{{\rm B}}^{{\rm max} }+\alpha {{x}_{1}}-\beta .c \\ {} & {} & +\delta .P(m_{*}^{{\rm true}}\lt m_{*}^{{\rm threshold}})-{{M}_{{\rm B}}} \\ \end{array}\end{eqnarray} \tag{ 12 }$$

where M_B is the absolute B-band magnitude of a Type Ia SN with x₁ = 0, c = 0, and $P(m_{*}^{{\rm true}}\lt m_{*}^{{\rm threshold}})\;=\;0$. Here $P(m_{*}^{{\rm true}}\lt m_{*}^{{\rm threshold}})=0$ denotes the host mass correction that accounts for the correlation between the luminosity of the SN and the host galaxy's mass (Suzuki et al. 2012).

The Constitution catalog contains N_Cons = 256 SNe, also obtained from a variety of surveys. It shares many of the same low-redshift SNe as Union 2.1; however, the method of dealing with systematics and choice of light-curve fitting procedure vary. The dominant systematics are identified as the choice of nearby training set used and the light-curve fitting procedure adopted. Specifically, the treatment of host galaxy reddening is the primary source of systematic uncertainty. To examine these effects, four different light-curve fitting procedures are utilized in Hicken et al. (2009). SALT2 is used in the same manner as in the Union 2.1 sample—where all host reddening effects are incorporated via the color term c, with the empirical relation ${{\mu }_{{\rm B}}}\propto c$.

The MLCS light-curve fitting method utilizes a different approach—fitting the SN distance in conjunction with a shape/luminosity parameter Δ and the host galaxy extinction parameter A_V. The extinction A_V is calculated using a prior on E(B−V) and a reddening law parameterized by R_V. Two different values of R_V are adopted in Hicken et al. (2009): R_V = 3.1 and 1.7. We note that there is some evidence that the MLCS treatment of color introduces a systematic error into the determination of the absolute brightness and can bias cosmological parameter estimation. In particular, it has been argued that the "Hubble bubble" detected in Zehavi et al. (1998), Conley et al. (2007), and Jha et al. (2007) might be due to an incorrect assumption regarding the value of R_V. For the purposes of this work we repeat our analysis using both the MLCS R_V = 1.7 and SALT2 Constitution catalogs.

The recently constructed Lick Observatory Supernova Search (LOSS) sample was developed in Ganeshalingam et al. (2013). The data consist of 165 BVRI light curves of low-redshift SNe, mainly sampled a week before maximum light in the B band. This set of objects was combined with data from the Calan/Tololo sample (Hamuy et al. 1996), CfA1–3 (Riess et al. 1999; Hicken et al. 2009; Jha et al. 2006), and CSP (Folatelli et al. 2010; Contreras et al. 2010). For SNe that are common to multiple surveys, the best sampled light curves are utilized. The SALT2 light-curve fitting procedure was adopted to construct the distance modulus of N_LOSS = 586 data in the range z = 0.01−1.4. Of the 226 low-z SNe, 91 are from LOSS and 45 had distances published for the first time in Ganeshalingam et al. (2013).

Finally, we use the combined catalog of Johnson et al. (2014). In Johnson et al. (2014) the aim was to construct a maximally homogeneous set of data points out to redshift z < 0.07. As such, the sample is drawn from a wide number of catalogs—specifically, the data consist of 40, 128, 135, 58, 33, and 26 SNe from the LOSS, Tonry et al. (2003), MLCS Constitution, Union, Kessler et al. (2009), and Carnegie Supernova Project (Hamuy et al. 2006; Folatelli et al. 2010) samples, respectively. When multiple measurements of the same object are known, the median value of the distance modulus is quoted. One arrives at a catalog containing N = 303 Type Ia SNe, with distance errors of order σ_d ∼ 5%.

A number of small modifications to each catalog are made. Since we are testing for bulk flows of expected order v_bulk ∼ 300–500 km s⁻¹, we must have precise knowledge of the SN redshifts. We therefore update each data point using host galaxy information found in the NASA Extragalactic Database (NED).⁵ We do not include any SNe in our analysis for which no host galaxy information is known. We also neglect any data for which the redshift error is larger than σ_z = 100 km s⁻¹. For SN 2002hu, we do not use the redshift quoted in the NED database. Rather, we adopt a more recent measurement (Blondin et al. 2012). Our use of the host galaxy redshifts is the primary reason for the small difference in results between this work and Colin et al. (2011).

The Constitution set adopts a different value of H_0,const = 65 km s⁻¹ Mpc⁻¹ during the light-curve fit procedure—for the purposes of combining data sets later, we add the term $\delta {{\mu }_{{\rm h}}}=5{\rm log} [{{H}_{0,\;{\rm const}}}/{{H}_{0,{\rm fid}}}]$ to each SN in this catalog, with H_0,fid = 70 km s⁻¹ Mpc⁻¹. All of the catalogs considered in this work introduce an error component in the distance modulus estimation to account for peculiar velocities, which is set at σ_v = 300 km s⁻¹ in the Union 2.1 and LOSS sets and σ_v = 400 km s⁻¹ for the Constitution sample. This contribution to the error is removed for the combined catalog of Johnson et al. (2014) and reintroduced in Equation (9) via the $\sigma _{{\rm v}}^{2}$ contribution to the covariance. In Johnson et al. (2014) σ_v was kept as a free parameter and allowed to vary with the cosmological parameters. It was found that a value of σ_v = 400km s⁻¹ was preferred by the data. We adopt this value whenever we use the combined catalog of Johnson et al. (2014).

We begin by calculating the ΔQ_d values and their associated p-values for each SN catalog, using realization set A (that is, neglecting large-scale velocity perturbations expected within the context of ΛCDM). For each catalog, we decompose the data into five unequally spaced redshift bins; see Table 1 for their limits and number of SNe contained within each bin. We perform our test using both concentric and cumulative redshift shells.

Our test function Q(θ, ϕ, δ) contains a free parameter δ, which determines the width of the smoothing on the unit sphere. Decreasing values of δ will allow us to probe bulk velocity distributions in successively smaller regions of the sky, corresponding to local flows due to small-scale overdensities. However, the significance of such events is expected to be low owing to the small number of participating SNe Ia. Since our statistic ΔQ_d was chosen to search for dipoles, and since a dipole signal would manifest itself as velocities coherent over entire hemispheres on the sky, we fix δ = π/2 in what follows.

We exhibit the ΔQ_d p-value and direction of maximal anisotropy (b_max, ℓ_max) in Table 2 for the Union 2.1, Constitution, LOSS, and combined sets. One can see consistent trends in all of the samples. The direction (b_max, ℓ_max) of maximal anisotropy in the cumulative redshift bins is consistent among the different data sets out to redshift z < 0.06. This is not a particularly surprising outcome, as the majority of low-redshift SNe are present in all four catalogs. However, it serves as a useful check that the previous detection of a bulk flow (Colin et al. 2011) is robust to both treatment of systematics and choice of light-curve fitting procedure. However, the significance of any anisotropic signal is modest, with p-values typically greater than p ∼ 0.1 in all catalogs and all redshift bins. Notable exceptions are the Union 2.1 and Combined samples, in cumulative redshift shell 0.015 < z < 0.045, where the anisotropic signal is largest.

We can state that for each individual catalog, there is no statistical significance of a bulk flow in any concentric redshift bin. The lack of significance is a reflection of the fact that there are an insufficient number of SNe in the individual shells, and coherent motion in a particular region of the sky cannot be strongly distinguished from random velocities when there are only a small number of data points.

It is of interest to compare the results of the Constitution set, analyzed using the two different light-curve fitting procedures. One can see broad consistency in both the directions (b_max, ℓ_max) and significance when applying the SALT II and MLCS light-curve fitting procedures to the same data sets. Both give comparable p-values and (b_max, ℓ_max) directions for the majority of the redshift bins considered. The only interesting deviation occurs in the redshift bin $0.045\leqslant z\lt 0.06$, in which the MLCS procedure assigns considerably higher significance to the anisotropy than the SALT II fit. This bin is of interest, and we discuss it in more detail.

In the $0.045\leqslant z\lt 0.06$ bin one can clearly see some consistency among the anisotropic directions (b_max, ℓ_max) in the Union and Constitution samples—what appears to be a turnaround in direction relative to the z < 0.045 bins. This turnaround is most significant for the Constitution set, using the MLCS R_V = 1.7 light-curve fit. This behavior was first observed in Colin et al. (2011) and was attributed to infall into the Shapley supercluster. However, it was argued in Feindt et al. (2013) that the observed turnaround was primarily due to the SNe SN 1995ac and SN 2003ic, which are located at almost maximal distance from Shapley on the sky. Hence, the effect may not be due to this supercluster, but instead either a random alignment of two large residuals in a particular region of the sky or infall into a different local overdensity. When constructing our residuals, we also find that the change in anisotropic direction in the $0.045\leqslant z\lt 0.06$ bin is due to the SN 1995ac and SN 2003ic data points—when removing them, we find that the direction of maximal isotropy changes to (b_max, ℓ_max) = (−14°, 90°) with insignificant p-value p = 0.492 for the MLCS17 catalog. One can conclude that there is no evidence of back-infall into Shapley in the bin $0.045\leqslant z\lt 0.06$. That the MLCS light-curve fit yields a higher significance than SALT II in this redshift bin is simply due to the fact that the MLCS approach estimates a considerably smaller distance modulus uncertainty for the offending SNe SN 1995ac and SN 2003ic—the observational uncertainty σ_μ is of order ∼30% and ∼15% smaller, respectively, when using MLCS compared to SALT II. Although current results show no evidence of infall into Shapley, more data are required in this redshift bin to definitively answer the question.

There is some evidence of a dipole in the cumulative redshift bins. The direction of the anisotropy remains qualitatively consistent in all cumulative bins and catalogs. The minimum p-value is observed in the $0.015\leqslant z\lt 0.045$ bin and combined and Union catalogs. Beyond z = 0.045, the significance of the detection drops owing to two factors: the small number of SNe in the higher-redshift bins and the effect of peculiar velocities becoming insignificant at high redshift (causing a decrease in the signal-to-noise ratio).

4.1. Consistency with ΛCDM

4.2. Effects due to Inhomogeneous Data

The method adopted in this work yields the p-value for the test function ΔQ_d and the direction of maximal anisotropy, which we denote (b_max, ℓ_max). We now consider how the direction of maximal anisotropy inferred by the smoothed residual method might be biased owing to an inhomogeneous distribution of data points.

A bulk velocity would manifest itself as a directional dependent signal in the distance moduli, specifically a dipole in the luminosity distance. The magnitude of the effect on any given data point will depend on both its position relative to the bulk flow direction and its redshift, and hence the direction of maximal anisotropy detected by our method will be affected by both the anisotropic and inhomogeneous nature of the data distribution on the sky.

To test how well one can reproduce the direction of an underlying anisotropic signal, we create mock realizations of the combined catalog in which we keep the data positions and redshifts fixed and introduce an artificial bulk flow. Specifically, the luminosity distance of the ith data point is estimated as (Bonvin et al. 2006)

$$\begin{eqnarray}&&{{d}_{{\rm L},{\rm i}}}=\frac{(1+{{z}_{{\rm i}}})c}{{{H}_{0}}}\int _{0}^{{{z}_{{\rm i}}}}\frac{d\bar{z}}{h(\bar{z})}+\frac{{{V}_{{\rm bulk}}}{{(1+{{z}_{{\rm i}}})}^{2}}}{H({{z}_{{\rm i}}})}{\rm cos} [{{\theta }_{{\rm i}}}]\end{eqnarray} \tag{ 13 }$$

where V_bulk is the magnitude of the bulk velocity and θ_i is the angle subtending the direction of the ith SN on the sky and the direction of the bulk velocity. The distance modulus of each data point is then constructed using this luminosity distance, and Gaussian noise is added to each point, using the square of the catalog observational uncertainty as the variance. We generate N_real = 15,000 realizations of the data—for each realization we input a bulk velocity into the luminosity distance (Equation (13)) of magnitude V_bulk = 400 km s⁻¹ and with a random direction on the sphere, (b_v, ℓ_v). We then calculate the difference between the actual bulk velocity direction (b_v, ℓ_v) and the reconstructed direction (b_max, ℓ_max) that yields a maximum value of ΔQ_D, for each realization. In the following analysis, we only keep realizations in which there is a "significant" detection of the inputted dipole, that is, those with a p-value p < 0.05 relative to simulations in which V_bulk = 0.

We perform our test using the combined catalog and redshift bin $0.015\leqslant z\lt 0.045$, which constitutes the redshift bin in which an anisotropic signal has the highest significance. For each realization we randomize the direction of the dipole isotropically on the unit sphere, with no priors imposed on (b_v, ℓ_v). In Figure 1 we exhibit the scatter between (b_v, b_max) and (ℓ_v, ℓ_max) for the realizations. For a homogeneous data set, where there would be no bias, the scatter in Figure 1 would be centered on b_v = b_max and ℓ_v = ℓ_max, denoted as black lines in the figures. Indeed, one can observe exactly such behavior in the (ℓ_max, ℓ_v) plot. However, there is a clear bias in the galactic latitude, with the method reconstructing a b_max that is consistently larger (smaller) than b_v in the north (south) galactic hemisphere. This is a reflection of the fact that the data are sparse in the region $|b|\lt 20{}^\circ$, and the method will preferentially select a region of the sky where the data are prevalent.

Zoom In Zoom Out Reset image size

We note two peculiarities in the (ℓ_v, ℓ_max) scatter. The concentrations of points in the top left and bottom right of the plot are due to the periodic nature of the longitude coordinate ℓ→ℓ + 360°, with low ℓ points scattering to ℓ ∼ 360° and vice versa. The small number of points that are otherwise scattered far from the expected ℓ_max = ℓ_v relation are predominantly those in which the underlying bulk velocity in the realization was placed at $|{{b}_{{\rm v}}}|\lt 20{}^\circ$; for these values the uncertainty on ℓ_v increases dramatically owing to the lack of data in this region.

In Figure 2 we exhibit (b_max−b_v) and (ℓ_max−ℓ_v) as empirical probability distributions, where we keep only realizations for which b_v > 20° (top panels) and b_v < −20° (bottom panels). One can see that the galactic longitude yields distributions consistent with ℓ_v = ℓ_max, but the latitude reconstruction is biased. We estimate the bias (Δb_max, Δℓ_max) in the north and south galactic planes by calculating the mean of the distributions in Figure 2, and the uncertainty (δb_max, δℓ_max) as the sample variance. We exhibit the values of (Δb, Δℓ) and (δb, δℓ) in the first row of Table 4.

Table 4.  Bias (Δb, Δℓ) and Uncertainty (δb, δℓ) in the Reconstructed Direction (b_max, ℓ_max) for Realizations with b_v > 20° and b_v < −20°.

	North (bv > 20°)		South (bv < −20°)
Vbulk (km s−1)	(Δb, Δℓ)	(δb, δℓ)	(Δb, Δℓ)	(δb, δℓ)
400	(13°, −°)	(14°, 28°)	(−12°, 2°)	(14°, 29°)
800	(15°, −4°)	(9°, 22°)	(−13°, 2°)	(9°, 21°)

Notes. We adopt two sets of realizations with an inputted bulk velocity magnitude of V_bulk = 400 km s⁻¹ and V_bulk = 800 km s⁻¹. We note that the bias in galactic longitude ℓ is negligible; however, the latitude b_max is systematically shifted away from the galactic plane in both the northern and southern hemispheres. The bias does not appreciably vary with V_bulk, indicating that Δb is due to the inhomogeneous nature of the data. The uncertainty (δb, δℓ) decreases with increasing V_bulk, as a result of increased signal-to-noise ratio.

Download table as:  ASCIITypeset image

We repeat our calculation using realizations in which a bulk velocity of magnitude V_bulk = 800 km s⁻¹ is introduced. In doing so, we obtain similar distributions to Figure 2. We find that the mean values of the distributions do not vary appreciably when we increase V_bulk, indicating that the bias (Δb_max, Δℓ_max) is not a strong function of the underlying bulk flow but rather an intrinsic feature of the data. However, as one might expect, the uncertainty (sample variance, (δb_max, δℓ_max)) of the distributions decreases, reflecting the signal-to-noise ratio increase associated with a larger bulk flow. In Table 4 we exhibit the values of the mean (bias) and sample variance (uncertainty) of the distributions in the north and south galactic planes, for both V_bulk = (400, 800) km s⁻¹ mock data sets.

Zoom In Zoom Out Reset image size

Finally, we use our realizations to estimate the bias and uncertainty associated with the measurement of (b_max, ℓ_max) obtained using the data—specifically the 0.015 < z < 0.045 redshift shell of the combined catalog. To do so, we select the subset of V_bulk = 400 km s⁻¹ realizations that have a reconstructed direction b_max in the range 34° < b_max < 42°, regardless of the actual bulk velocity latitude b_v (recall that the data yielded a direction of maximal anisotropy b_max = 38°; see Table 2). We construct a distribution of (b_max−b_v) of this subset, which we exhibit in Figure 3. We estimate the bias Δb_max and uncertainty δb_max on our data measurement b_max = 38° as the mean and sample variance of this distribution, respectively. We find Δb_max = 18° and δb_max = 12°. As the longitude ℓ_max has been shown to be unbiased, we take Δℓ_max = 0° and δℓ_max = 29° from the whole sample. Hence, for the combined catalog and redshift bin 0.015 < z < 0.045, we estimate the direction of maximal anisotropy to be (b_max,ℓ_max) = (20°, 276°) ± (12°, 29°). We note that the final reconstructed direction is in agreement with various theoretical and observational works (Campanelli 2009; Feindt et al. 2013; Gibelyou & Huterer 2012). Our result is also in qualitative agreement with recent work by the authors (Appleby & Shafieloo 2014), where a nonparametric reconstruction of the bulk flow direction was undertaken using the local galaxy distribution and the associated luminosity function.

Zoom In Zoom Out Reset image size

Although our analysis raises concerns as to the reliability of the b_max values presented in Table 2, one should recall that the method is primarily designed as a null hypothesis test. The principle result of the method are the p-values obtained for each catalog and redshift bin—small values indicate any kind of violation of the isotropic hypothesis. One can make no strong claim as to the direction or underlying cause of anisotropy when performing a null test for the current data sets—we can simply state whether the data are consistent with the null hypothesis. However, by using mock data sets and making assumptions as to the nature of the anisotropy, we are able to make a quantitative statement regarding the accuracy of our recovered (b_max, ℓ_max) directions. Clearly, both the bias and the uncertainty in (b_max, ℓ_max) will depend on the data set used.

In this work we have used a nonparametric method to test the underlying assumption of isotropy in the low-redshift SN data sets. Our first test, using realization set A, is tantamount to the question, given the statistical, systematic, and astrophysical uncertainties associated with the SN measurements, are the data consistent with flat, isotropic ΛCDM expansion with zero bulk flow? Our answer to this question is that there is a hint in the Union and combined catalogs of an anisotropic signal, which is most pronounced in the 0.015 < z < 0.045 redshift bin, with p-value p = 0.07. Owing to the small number of data points, tomography yields no significant findings in any redshift bin. When we compare the data to realization set B, where large-scale velocity perturbations within the context of ΛCDM are accounted for, we find that there is no hint of any anomalous deviation from the null hypothesis, with the smallest p-value p = 0.29. One can conclude that there is some evidence of a dipole in the SN catalogs; however, this bulk flow is consistent with the standard cosmology.

We attempted to use the method of smoothed residuals to estimate the direction of the bulk flow; however, the reconstructed direction will be biased by the inhomogeneous distribution of data on the sky. To estimate this effect, we created mock realizations of the data in which we inserted a bulk flow of magnitude V_bulk = 400 km s⁻¹ and random direction (b_v, ℓ_v). We then used our method to attempt to reconstruct this direction. We found that the reconstructed galactic longitude ℓ_max is consistent with the inputted value ℓ_v; however, the recovered galactic latitude b_max is systematically biased away from the galactic plane. Such behavior is expected given the sparsity of data in the region $|b|\lt 20{}^\circ$, and we expect that such a bias is unavoidable in reconstructions when the data are inhomogeneously distributed. After correcting for a bias of (Δb, Δℓ) = (18°, 0°), we quote the direction of maximal anisotropy as (b_max, ℓ_max) = (20°, 276°) ± (12°, 29°) for the combined catalog and redshift shell 0.015 < z < 0.045.

We note that our result is qualitatively consistent with recent work (Feindt et al. 2013). In Feindt et al. (2013), two approaches are used to conclude that the Union 2 sample contains a dipole contribution. The first of these approaches is an explicit dipole fit to the data—using Equation (13) and fitting for the magnitude of the bulk velocity V_bulk and its direction on the sky (b_v, ℓ_v). Our analysis does not preclude the existence of a detectable dipole—indeed, both our approach and the dipole fit yield a qualitatively consistent direction (b_max, ℓ_max) of maximal anisotropy once we have eliminated the bias in our estimation.

It is clear that there is some evidence of an anisotropic signal in the SN data. However, more data and better control of the observational uncertainties are required before one can pin down the magnitude, direction, and redshift extent of the local bulk flow. However, one can state that the local bulk flow observed in current SN data is consistent with the ΛCDM model.

The authors would like to thank Subir Sakar, Chris Blake, and Jeppe Nielsen for helpful discussions. S.A.A. and A.S. wish to acknowledge support from the Korea Ministry of Education, Science and Technology, Gyeongsangbuk-Do, and Pohang City for Independent Junior Research Groups at the Asia Pacific Center for Theoretical Physics. A.S. would like to acknowledge the support of the National Research Foundation of Korea (NRF-2013R1A1A2013795). A.J. is supported by the Australian Research Council Centre of Excellence for All-Sky Astrophysics (CAASTRO) through project number CE110001020.