3 The Lucy-Richardson inversion algorithm applied to Abel's integral equation

The Lucy-Richardson algorithm provides another robust way to invert Eq. (4) (see also Cerf & Boffin 1994). The method starts from the Bayes theorem on conditional probability in the form

$$\Psi(M_2) \; \Pi(Y\vert M_2) = \Phi(Y) \; R(M_2 \vert Y),$$ (9)

where R(M₂ | Y) is the reciprocal kernel corresponding to the integral equation inverse to the one that needs to be solved (Eq. (1)):

$$\Psi(M_2) = \int_0^{M_2} \Phi(Y)\; R( M_2\vert Y)\; {\rm d}Y.$$ (10)

The reciprocal kernel represents the conditional probability that the binary system has a companion mass M₂ when the observed  $M_2 \sin i$ value amounts to Y. Thus, one has:

$$\frac{\Psi(M_2) \; \Pi(Y\vert M_2)}{\Phi(Y)}$$ R(M₂|Y) = (11)

$$\frac{\Psi(M_2) \; \Pi(Y\vert M_2)}{\int_0^{\infty} \Psi(M_2)\; \Pi(Y\vert M_2)\; {\rm d}M_2},$$ = (12)

which obviously satisfies the normalization condition $\int_0^{\infty} R(M_2\vert Y)\; {\rm d}M_2 = 1$. The problem in solving Eq. (10) arises because R(M₂|Y) also depends on $\Psi (M_2)$, so that an iterative procedure must be used. If $\Psi_r(M_2)$ is the rth estimate of $\Psi (M_2)$, it can be used to obtain the (r+1)th estimate in the following way:

$$\Psi_{r+1}(M_2) = \int_0^{M_2} \Phi(Y)\; R_r(M_2\vert Y)\; {\rm d}Y$$ (13)

with

$$R_r(M_2\vert Y) = \frac{\Psi_r(M_2)\; \Pi(Y\vert M_2)}{\Phi_r(Y)}$$ (14)

and

$$\Phi_r(Y) = \int_0^{\infty} \Psi_r(M_2)\; \Pi(Y\vert M_2)\; {\rm d}M_2.$$ (15)

Thus, $\Phi_r(Y)$ represents the corresponding rth estimate of the observed distribution $\Phi (Y)$. Equations (13) and (14) together yield the recurrence relation for the $\Psi_r$'s,

$$\Psi_{r+1}(M_2) = \Psi_r(M_2) \; \int_0^{M_2} \frac{\Phi(Y)}{\Phi_r(Y)}\; \Pi(Y\vert M_2)\; {\rm d}Y,$$ (16)

with $\Pi(Y \vert M_2)$ given by Eq. (3) for the problem under consideration. The conditions for convergence of this recurrence relation are discussed by Lucy (1974) and Cerf & Boffin (1994). It needs only be remarked here that (i) the iterative scheme converges if $\Phi_r(Y)$ tends to $\Phi (Y)$, given the normalization of the probability $\Pi(Y\vert M_2) {\rm d}Y$, and (ii) the full convergence of the method is not necessarily desirable, as the successive estimates $\Phi_r(Y)$ will tend to match $\Phi (Y)$ on increasingly smaller scales, but the small-scale structure in $\Phi (Y)$ is likely to be dominated by the noise in the input data. This is well illustrated in Fig. 2 below.

When the number of data points is small (typically N <100; Cerf & Boffin 1994), it is advantageous to express $\Phi (Y)$ as

$$\Phi(Y) = \frac{1}{N}\sum_{n = 1}^{N} \delta(Y - y_n)$$ (17)

where the $y_n \; (n=1,\dots~ N)$ are the N individual measured  $M_2 \sin i$ values and $\delta(x)$ is the Dirac "function'' such that $x_0 = \int \delta (x - x_0)\; {\rm d}x$. Substitution in Eq. (13) then yields

$$\Psi_{r+1}(M_2) = \frac{1}{N}\sum_{n=1}^N R_r(M_2\vert y_n)$$ (18)

where R_r(M₂|y_n) is defined as in Eq. (14). The sample size should nevertheless be large enough for the functions R_r(M₂|y_n) to have sufficient overlap so as to produce a smooth $\Psi_{r+1}(M_2)$ function.

In the application of the method described in Sect. 4, the initial mass distribution $\Psi_0(M_2)$was taken as a uniform distribution, but it has been verified that this choice has no influence on the final solution.