Reliable computation of a multiple integral involved in the neutron star theory

Abstract

The following multiple integral is involved in the neutron star theory:

$$\tau (\epsilon \text{,}v)=\frac{1}{\omega (\epsilon )}{\int }_0^{\pi /2}\mathrm{d}\theta \sin (\theta ){\int }_0^{\infty }\mathrm{d}n{n}^2{\int }_0^{\infty }\mathrm{d}ph(n\text{,}p\text{,}\theta \text{,}\epsilon \text{,}v)$$

where

$$h(n\text{,}p\text{,}\theta \text{,}\epsilon \text{,}v)=\psi (z)\varphi (n-\epsilon -z)+\psi (-z)\varphi (n-\epsilon +z)-\psi (z)\varphi (n+\epsilon -z)-\psi (z)\varphi (n+\epsilon +z)$$

and

$$z=\sqrt{p^2+{(v\sin (\theta ))}^2}\text{,}\psi (x)=\frac{1}{\exp x+1}\text{,}\varphi (x)=\frac{x}{\exp x-1}\text{.}$$

$\omega (\epsilon )$ is a normalization function.

The aim is to get a table for $\tau (\epsilon \text{,}v)$ for some values of $(\epsilon \text{,}v)$ in $[10^{-4}\text{,}10^4]\times [10^{-4}\text{,}10^3]$ and then to interpolate for the others. We present a new strategy, using the Gauss–Legendre quadrature formula, which allows to have one code available whatever the values of $v$ and $\epsilon$ are. We guarantee the accuracy of the final result including both the truncation error and the round-off error using Discrete Stochastic Arithmetic.

Introduction

This paper describes how to dynamically control the computation of a multiple improper integral. This integral, which is involved in the neutron star theory, has been established by Villain and Haensel [23], [24]. First we briefly present in which physical context this integral arises.

For any isolated star, the final stages of life depend mostly on its initial mass. The current idea is that massive stars end up within the collapse of their sterile iron core. This collapse is followed by a final explosion which ejects all the matter, except a central residue. Depending on the details of the scenario, this residue is a neutron star (NS) or a black hole. As the global electric charge of a NS has to be zero, there are as many electrons as protons inside. Due to the high density of a NS, most of these electrons interact with protons to give neutrons. Thus, NS are mainly done of neutrons with a typical proportion of protons smaller than 10%. With time, equilibrium between protons and neutrons, called beta equilibrium, takes place. Time to reach beta equilibrium is obtained by calculating high order integrals involving several variables for any particle appearing in the reactions. These variables are roughly the direction of their motion, their velocity and their energy. Moreover, NS involve a lot of physical phenomena that will make this calculation more complicate. For instance, at such densities, nucleons may be super-fluid. The main effect of super-fluidity in these calculations is that it excludes some values of the possible energy of nucleons, and sometimes only for some directions. There is then a gap in the spectrum of energy and a possible breaking of symmetries; the integrals involved become much more difficult to calculate and numerical calculation is unavoidable. Super-fluidity has been taken into account to establish the integral studied in this paper. It is a three-dimensional improper integral which depends on two parameters: the amplitude of the gap and another variable describing how far the matter is from beta equilibrium.

There are other fields in physics where the computation of an integral is essential. For instance, in 1981, Harrison [9] evaluated by numerical integration an improper integral arising in a study of the total electronic energy of crystals using the tight-binding approximation. The evaluation of this integral was also a problem posed in [15]. More recently, Benkaci and Maugin [1] presented a study of piezo-ceramic materials, where integrals were computed using a Gaussian quadrature method. Verified quadrature methods have been used for integrals involved in some physical problems, such as the determination of Newton’s constant of gravitation [10], [13] or the evaluation of the transport coefficients of polar gases [17].

This paper is dedicated to the computation of the following integral:

$$\tau (\epsilon \text{,}v)=\frac{1}{\omega (\epsilon )}{\int }_0^{\pi /2}\mathrm{d}\theta \sin (\theta ){\int }_0^{\infty }\mathrm{d}n{n}^2{\int }_0^{\infty }\mathrm{d}ph(n\text{,}p\text{,}\theta \text{,}\epsilon \text{,}v)\text{,}$$

where

$$(\epsilon \text{,}v)\in [10^{-4}\text{,}10^4]\times [10^{-4}\text{,}10^3]$$

$$h(n\text{,}p\text{,}\theta \text{,}\epsilon \text{,}v)=\psi (z)\varphi (n-\epsilon -z)+\psi (-z)\varphi (n-\epsilon +z)-\psi (z)\varphi (n+\epsilon -z)-\psi (z)\varphi (n+\epsilon +z)$$

with

$$z=\sqrt{p^2+{(v\sin (\theta ))}^2}\text{,}\psi (x)=\frac{1}{\exp x+1}\text{,}\varphi (x)=\frac{x}{\exp x-1}$$

and

$$\omega (\epsilon )=\frac{17{\pi }^4}{60}\epsilon \left(1+\frac{10{\epsilon }^2}{17{\pi }^2}+\frac{{\epsilon }^4}{17{\pi }^4}\right)\text{.}$$

$\omega$ is a normalization function, which has been defined such that $\tau (\epsilon \text{,}0)=1$. Our aim is to get a table for $\tau (\epsilon \text{,}v)$ for some values of $(\epsilon \text{,}v)$ and to interpolate for the others, in order to determine the surface $\tau (\epsilon \text{,}v)$. We present in this paper how to compute a guaranteed value of the integral $\tau (\epsilon \text{,}v)$.

This paper is organized as follows. In Section 2, we present the expression of the improper integral to compute and the numerical problems which must be taken into account. As this multiple integral can be expressed as several one-dimensional integrals, it can be computed using a quadrature method. In Section 3, we show how the Gauss–Legendre method can be dynamically controlled. In Section 4, we briefly review methods and concepts which enable one to estimate round-off error propagation with a probabilistic approach. In Section 5, we describe a strategy to dynamically control the computation of the integral presented in Section 2. Furthermore, we show that we can determine in the result obtained the significant digits common with the exact value of the integral. The last section presents numerical experiments carried out using the strategy described in Section 5.