Analytical expressions and numerical evaluation of the luminosity distance in a flat cosmology

TL;DR

This paper introduces a new hypergeometric function-based expression for calculating luminosity distances in flat mbda CDM cosmology, offering improved computational efficiency especially at high redshifts, and compares it with existing methods.

Contribution

It presents a novel, simple hypergeometric function expression for luminosity distance in flat mbda CDM cosmology, enhancing computational efficiency over existing elliptic integral methods.

Findings

The new expression is the most efficient for redshifts z  1.

Combining the new formula with Carlson's elliptic integral approach covers the full redshift range.

Approximate fitting formulae provide acceptable accuracy for practical purposes.

Abstract

Accurate and efficient methods to evaluate cosmological distances are an important tool in modern precision cosmology. In a flat $\Lambda$CDM cosmology, the luminosity distance can be expressed in terms of elliptic integrals. We derive an alternative and simple expression for the luminosity distance in a flat $\Lambda$CDM based on hypergeometric functions. Using a timing experiment we compare the computation time for the numerical evaluation of the various exact formulae, as well as for two approximate fitting formulae available in the literature. We find that our novel expression is the most efficient exact expression in the redshift range $z\gtrsim1$. Ideally, it can be combined with the expression based on Carlson's elliptic integrals in the range $z\lesssim1$ for high precision cosmology distance calculations over the entire redshift range. On the other hand, for practical work where relative errors of about 0.1% are acceptable, the analytical approximation proposed by Adachi & Kasai (2012) is a suitable alternative.

Equations

$$d_{\text{L}}(z)=\frac{c(1+z)}{H_0}{\int }_0^z\frac{\text{d}t}{\sqrt{{\Omega }_{\text{m}}(1+t{)}^3+(1-{\Omega }_{\text{m}})}}$$

$$s=\sqrt[3]{\frac{1-{\Omega }_{\text{m}}}{{\Omega }_{\text{m}}}}$$

$$d_{\text{L}}(z)=\frac{c(1+z)}{H_0\sqrt{s{\Omega }_{\text{m}}}}\left[T(s)-T\left(\frac{s}{1+z}\right)\right]$$

$$T(x)={\int }_0^x\frac{\text{d}u}{\sqrt{u(1+u^3)}}$$

$$T_{\infty }={\int }_0^{\infty }\frac{\text{d}u}{\sqrt{u(1+u^3)}}=\frac{\Gamma (\frac{1}{6})\Gamma (\frac{1}{3})}{3\sqrt{\pi }}$$

$$T(x)\sim 2\sqrt{x}\left(1-\frac{1}{14}{x}^3+\cdots \right)$$

$$T(x)\sim T_{\infty }-\frac{1}{x}+\frac{1}{8x^4}+\cdots$$

$$T(x)=\frac{1}{\sqrt[4]{3}}F\left(\arccos \left(\frac{1+(1-\sqrt{3})x}{1+(1+\sqrt{3})x}\right),\cos \left(\frac{\pi }{12}\right)\right)$$

$$T(x)=4R_{\text{F}}\left(m,m+3-2\sqrt{3},m+3+2\sqrt{3}\right)$$

$$m=\frac{2\sqrt{x^2-x+1}}{x}+\frac{2}{x}-1$$

$$t={\left(\frac{u}{x}\right)}^3$$

$$T(x)=\frac{1}{3}\sqrt{x}{\int }_0^1\frac{\text{d}t}{t^{5/6}{\left(1+x^{3}t\right)}^{1/2}}$$

$$T(x)=2\sqrt{x}{}_2{F}_1\left(\frac{1}{6},\frac{1}{2};\frac{7}{6};-x^3\right)$$

$$T(x)=\{\begin{array}{ll}2\sqrt{\frac{x}{1+x^3}}{}_2{F}_1\left(\frac{1}{2},1;\frac{7}{6};\frac{x^3}{1+x^3}\right)& x⩽1\\ T_{\infty }-\sqrt{\frac{x}{1+x^3}}{}_2{F}_1\left(\frac{1}{2},1;\frac{4}{3};\frac{1}{1+x^3}\right)& x>1\end{array}$$

$$T(x)\approx \frac{2\sqrt{x}}{\sqrt[8]{1+a_{1}x+a_2{x}^2+a_3{x}^4+a_4{x}^4}}$$

$$a_1$$

$$a_2$$

$$a_3$$

$$a_4$$

$$T(x)\approx \sqrt{x}\left(\frac{2+b_1{x}^3+b_2{x}^6+b_3{x}^9}{1+c_1{x}^3+c_2{x}^6+c_3{x}^9}\right)$$

$$b_1$$

$$b_2$$

$$b_3$$

$$c_1$$

$$c_2$$

$$c_3$$

Full text

††pubyear: 2017††pagerange: Analytical expressions and numerical evaluation of the luminosity distance in a flat cosmology–References

1 Introduction

The calculation of cosmological distances is one of the most fundamental tasks in cosmological studies. As the conversion between distance and redshift depends on the parameters of the underlying cosmological model, distance measurements are one of the key ingredients for cosmological tests. Accurate and efficient methods to evaluate these distances are an important tool in modern precision cosmology.

The relation between cosmological distance and redshift can be derived from the solution of the Friedmann equation, and involves an integral over the expansion history that depends on the cosmological model (e.g., Carroll et al., 1992). In a general $\Lambda$CDM cosmological model, the luminosity distance, possibly the most important distance scale from an observational point of view111Other distance scales such as angular diameter distance or proper distance are easily calculated from the luminosity distance (e.g., Weinberg, 1972)., cannot be expressed as a simple analytical formula of redshift and the cosmological parameters. Even in a cosmological model with zero curvature, as the current observations convincingly suggest (Hinshaw et al., 2013; Planck Collaboration, 2016), the luminosity distance–redshift relation can only be written as a nontrivial integral that cannot be evaluated using elementary functions.

Obviously, this integral can be evaluated numerically using standard quadrature algorithms, but these can become computationally demanding when high numerical accuracy is required.

One option to avoid numerical quadrature is to use analytical approximations for the luminosity distance. Several analytical recipes to approximate the luminosity distance in a flat cosmology have been put forward. Pen (1999) presented an algebraic fitting formula that has a relative error of less that 0.4% for $0.2<\Omega_{\text{m}}<1$ for any redshift, and a global relative error of less than 4% over the entire parameter space. Wickramasinghe & Ukwatta (2010) and Liu et al. (2011) used a similar approach and presented alternative analytical approximations that are show smaller error variations than Pen (1999). Finally, Adachi & Kasai (2012) used the technique of Padé approximants to come to an analytical formula with even smaller error variations (see also Wei et al., 2014).

An alternative option, particularly when higher accuracy is required, is to make use of exact analytical expressions for the luminosity distance. These will necessarily involve transcendental or special functions, which are more demanding to evaluate numerically than elementary functions. However, numerical software libraries such as Boost222http://www.boost.org (Schäling, 2014), NAG333http://www.nag.com (Phillips, 1987), GSL444http://www.gnu.org/software/gsl/ (Galassi et al., 2011), or SciPy555http://www.scipy.org (Oliphant, 2007) contain specialised algorithms to evaluate such special functions in a very efficient way, up to arbitrary precision. In the past few years, several authors have presented analytical expressions for the luminosity distance in a flat cosmology, all involving elliptic integrals (Eisenstein, 1997; Mészáros & Řípa, 2013; Liu et al., 2011; Zaninetti, 2016).

In this paper we derive and present an alternative and simple expression for the luminosity distance in a flat $\Lambda$CDM cosmology, based on hypergeometric functions. This is done in Section 2. In Section 3 we test the numerical performance of this new formula against the existing exact formulae, and against two of the approximate recipes that have been proposed. Section 4 sums up.

2 Analytical expressions

The general expression for the luminosity distance in a flat $\Lambda$CDM cosmology is given by

[TABLE]

where $\Omega_{\text{m}}=1-\Omega_{\Lambda}$ is the matter density of the Universe. If we introduce

[TABLE]

and we apply the substitution $u=s/(1+t)$, this expression can conveniently be written as

[TABLE]

with the function $T(x)$ defined as

[TABLE]

In spite of its apparent simplicity, this integral cannot be expressed in terms of elementary functions. The function $T(x)$ is a smooth function that continuously rises from $T_{0}=0$ at $x=0$ to a finite value

[TABLE]

in the limit $x\rightarrow\infty$ (see also Pen, 1999). The asymptotic behaviour at small and large values of $x$ is

[TABLE]

2.1 Expressions in terms of elliptic integrals

Recently, Mészáros & Řípa (2013) derived an analytical expression for the function $T(x)$,

[TABLE]

where $F(\phi,k)$ is the Legendre incomplete elliptical integral of the first kind.666Mészáros & Řípa (2013) mentioned that they found nothing in the literature about the non-numerical integration of this equation, and hence that their result is new and original. In fact, exactly the same equation had already been presented by Eisenstein (1997), but, although well cited, this contribution apparently never appeared in the refereed literature. Legendre elliptic integrals are typically evaluated numerically using infinite series or polynomial expansions, or using Newton-Raphson integration schemes (Byrd & Friedman, 1971; Lemczyk & Yovanovich, 1988).

Two years before Mészáros & Řípa (2013), Liu et al. (2011) presented an alternative analytical expression for the function $T(x)$:

[TABLE]

with

[TABLE]

Rather than the most common form of elliptical integrals, the standard Legendre format, this expression uses the Carlson elliptic integral of the first kind $R_{F}(x,y,z)$ (Carlson, 1977). The advantage of using Carlson’s form of elliptical integrals compared to the standard Legendre form is that the former have a number of interesting symmetry properties. In particular the so-called duplication theorem is extremely useful, and guarantees that the integrals can be calculated in a fast and robust way (Carlson, 1979, 1995).

For the sake of completeness, a third formulation in terms of the elliptic integrals was recently presented by Zaninetti (2016). This expression was substantially more complicated than the expressions given above, and involves a Legendre incomplete elliptical integral of the first kind with complex arguments. It therefore seems of little practical use.

2.2 Expressions in terms of hypergeometric functions

An alternative, and to the best of our knowledge, novel expression for the function $T(x)$, and hence for the luminosity distance in a flat cosmology, can be obtained by applying the transformation

[TABLE]

to expression (4). This yields

[TABLE]

This integral can be recognised as an Euler-type representation of a hypergeometric function. It results in the simple expression

[TABLE]

In general, the numerical evaluation of the hypergeometric function is notoriously difficult, in particular when one has the ambition of finding an expression for all complex values for the parameters and the argument (e.g., Michel & Stoitsov, 2008; Doornik, 2015). The main reasons are that the power series expansion of the hypergeometric function is thwarted by instabilities induced by cancellations between very large terms, and that certain regions in the complex plane are hard to include into the convergence domain (Bühring, 1986; López & Temme, 2013).

For the case of purely real parameters that avoid certain combinations, such as equation (13), the numerical evaluation is much more straightforward. The main disadvantage of this expression is that the argument of the hypergeometric function is negative and that it can fall outside the convergence domain of the corresponding hypergeometric series. However, it is easy to apply hypergeometric function transformation formulae to obtain equivalent expressions in which the argument falls within the convergence domain. Following the recommendations from Forrey (1997), one obtains

[TABLE]

For all positive $x$, the argument of the hypergeometric functions in (14) falls in the range between zero and one-half. This means that the power expansion series not only converges, but that it converges rapidly.

3 Numerical tests

We implemented the three expressions (8), (9) and (14) in C++, in order to test their relative performance. We used the implementations of the Legendre and Carlson elliptic integrals and the hypergeometric function from the GNU Scientific Library (GSL).

As a comparison, we also implemented two of the analytical approximations that have appeared in the literature. First, we used the original approximation put forward by Pen (1999). Converting his notation to the one used in this paper, the approximation can be written as

[TABLE]

Secondly, we use the approximation based on Padé approximants proposed by Adachi & Kasai (2012). In our notation,

[TABLE]

3.1 Accuracy

The accuracy of the three different analytical expressions was readily verified: over the entire range, the three expressions yielded values that are compatible up to the last significant digit. Their accuracy was also tested by evaluating the same expressions using Mathematica with much higher precision. For what concerns the approximate solutions (15) and (16), we confirm the findings already discussed in the literature: the approximations are generally accurate to less than half a percent at the lowest redshifts, and typically an order of magnitude more accurate at $z>1$ (see e.g. Fig. 1 of Adachi & Kasai, 2012).

3.2 Speed

The main goal of our numerical tests is a comparison of the efficiency of the different exact formulas to evaluate the function $T(x)$. For this goal, we set up a timing experiment, where we accurately timed the execution time of each of the three exact formulae and the two approximations. For this experiment, we fixed the value of the matter density to $\Omega_{\text{m}}=0.308$ (Planck Collaboration, 2016). We calculated each implementation 10 million times in 101 points distributed logarithmically in redshift space, between $z_{\text{min}}=0.01$ and $z_{\text{max}}=1000$. We used a setup of the timing exercise similar to the one described in Saftly et al. (2013). The actual run times were determined using the chrono functionality available in the C++11 standard library (Gregoire et al., 2011, Ch. 16).

The results of this timing experiment are shown in Figure 1. The evolution of the computation time as a function of redshift for both expressions based on the elliptic integrals has a similar behaviour: the computations are relatively expensive at low redshifts, and gradually become more efficient at larger $z$. The pattern is characterised by a number of discrete jumps between different plateaus, corresponding to a decrease or increase of the number of iterations in the calculation. The corresponding line for the hypergeometric formula (14) shows quite a different pattern. The computation time first increases as $z$ increases, until it reaches a maximum value, and then it continuously decreases sharply with increasing $z$. The maximum corresponds to $x=1$, and hence to the point where the calculation swaps from the second to the first branch in expression (14). At this maximum, the argument of the hypergeometric functions is equal to one-half, and this corresponds to the poorest convergence of the power series. For $x$ increasingly further away from this value on either side, the argument of the hypergeometric function gets increasingly smaller, and hence fewer terms in the power series need to be calculated to reach convergence. Finally, the computation time for each of the approximations is obviously independent of the redshift, as a fixed number of operations needs to be performed, without a convergence criterion.

Comparing the actual values of the timing, it is clear that the expression using Legendre’s elliptic integral derived by Eisenstein (1997) and Mészáros & Řípa (2013) has the poorest efficiency overall. It is about two times less efficient than the version of Liu et al. (2011) based on Carlson’s symmetric elliptic integral. This is not very surprising given the efficiency with which Carlson’s elliptic integrals can be evaluated (Carlson, 1979, 1995). More surprising is that our novel expression using the hypergeometric function is numerically the fastest one for $z\gtrsim 1$. For redshifts between about 3 and 100, it is another factor two faster than the implementation based on Carlson’s elliptic integral. Interestingly, the new analytical formula is as fast as the simple approximation by Pen (1999) for $z\gtrsim 10$ and even more efficient for $z\gtrsim 60$. The rational function approximation by Adachi & Kasai (2012), however, is always the fastest method to evaluate the luminosity distance (but it remains an approximation of course). In the low redshift regime ($z\lesssim 1$) it beats the exact methods by an order of magnitude, and also at high redshift it remains a factor three faster than the new method based on hypergeometric functions.

4 Discussion and conclusions

We have derived a novel and simple analytical expression for the luminosity distance in a flat $\Lambda$CDM cosmology that makes use of hypergeometric functions. Given that these functions are readily available from different numerical software libraries, our expression forms a useful alternative for quadrature formulae (which can be numerically expensive) or for analytical approximations (which remain approximations).

Apart from our novel formula, we have implemented two other exact analytical formulae for the luminosity distance, both based on elliptic integrals, and two analytical approximate fitting formulae. We have set up a timing experiment to determine the computation speed of the various possible exact and approximate formulae. The results of this experiment, and the corresponding implications, are as follows:

• •

The expression using Legendre’s elliptic integral derived by derived by Eisenstein (1997) and Mészáros & Řípa (2013) is the least efficient method of the three exact methods. Its use is hence not encouraged for numerical evaluation.

• •

The expression using Carlson’s symmetric elliptic integral, first derived by Liu et al. (2011), is a factor two more efficient than the one using Legendre’s elliptic integral. It is the most efficient exact formula at $z\lesssim 1$. On the other hand, our novel expression involving the hypergeometric function is about a factor two faster at $z\gtrsim 1$. Ideally, both functions can be combined for high precision cosmology distance calculations.

• •

The rational function approximation of the luminosity distance proposed by Adachi & Kasai (2012) remains a significantly faster alternative. It is a factor three faster than the relatively popular fitting function proposed by Pen (1999). The approximation is roughly an order of magnitude faster than any exact method at $z\lesssim 1$ and still a factor of three at the highest redshifts. For practical work where relative errors of $\sim 0.1\%$ are allowed, this seems the most intelligent option.

As a final remark, we note that the results obtained in this paper, which strictly deal with a flat cosmology, can also be useful for a more general $\Lambda$CDM cosmology. For example, Mészáros & Řípa (2015) showed that the expression for the luminosity distance in cosmological models with a small curvature can expanded in a power series, and that the dominant term is equivalent to the expression studied in this paper.

Acknowledgements

D.V.D.P. acknowledges the support from the Special Research Fund (BOF) of Ghent University. The authors thank the anonymous referee for useful comments.