# Radius measurement in binary stars: simulations of intensity interferometry

Km Nitu Rai,1 Soumen Basak,1 and Prasenjit Saha2

1School of Physics, Indian Institute of Science Education and Research
Thiruvananthapuram, Maruthamala PO, Vithura,
Thiruvananthapuram 695551, Kerala, India
2Physik-Institut, University of Zurich, Winterthurerstrasse 190, CH-8057
Zurich, Switzerland E-mail<EMAIL_ADDRESS>
###### Abstract

Mass and radius measurements of stars are important inputs for models of
stellar structure. Binary stars are of particular interest in this regard,
because astrometry and spectroscopy of a binary together provide the masses of
both stars as well as the distance to the system, while interferometry can
both improve the astrometry and measure the radii of the stars. In this work
we simulate parameter recovery from intensity interferometry, especially the
challenge of disentangling the radii of two stars from their combined
interferometric signal. Two approaches are considered: separation of the
visibility contributions of each star with the help of differing brightness
ratios at different wavelengths, and direct fitting of the intensity
correlation to a multi-parameter model. Full image reconstructions is not
attempted. Measurement of angular radii, angular separation and first-order
limb-darkening appears readily achievable for bright binary stars with current
instrumentation.

###### keywords:

(stars:)binaries: general, instrumentation: interferometers

##  1 Introduction

The mass and radii of stars are nowadays most often determined from their
asteroseismic frequencies, using scaling relations (e.g., Kjeldsen & Bedding,
1995) that relate stellar parameters to the solar values, or by comparing the
frequencies with stellar-structure models (e.g., Rodrigues et al., 2017). As
the asteroseismic methods are indirect, it is desirable to compare them,
wherever feasible, with dynamical mass measurements and geometrical radius
measurements, which directs attention to binary stars.

It has been known since Russell (1928) that combining astrometric and
photometric observations of a binary stellar system yields the masses of both
stars, as well as the three-dimensional orbit and the distance to the system.
This method has long been used to measure the masses of stars in resolved
binaries (e.g. Stephenson & Sanwal, 1969) and extended to binaries resolved by
interferometry (e.g. Docobo et al., 2018). A further measurement of the radii
of both stars would provide two mass-radius points.

Herbison-Evans et al. (1971) already succeeded in one mass-radius measurement.
They observed the binary $\alpha\rm\,Vir$ (Spica) using optical intensity
interferometry at mas-scale resolution, as well as spectroscopy. Combining
these observations they inferred the distance, the masses of both stars and
the radius of the brighter star, as well as all the orbital parameters (see
their Table III). Only the radius of the second star was lacking.
Unfortunately the progress of intensity interferometry stopped soon after, as
the interest of the comunity shifted to other kinds of interferometry, and an
interferometric measurement of the radius of the second star still remains to
be done. However, a renewed interest in intensity interferometry and new
observations resolving single stars (e.g., Acciari et al., 2020; Rivet et al.,
2020; Abeysekara et al., 2020) suggest that it will soon become possible to
use this technique to do mass-radius determinations in binary stars. The great
advantage of intensity interferometry is that it does not require sub-
wavelength precision on the mirrors or other optical components, and hence can
be extended to very long interferometric baselines. Baselines of $>1\rm\,km$,
corresponding to resolutions of $<1\rm\,mas$ are planned (Dravins et al.,
2013) and still longer baselines have been advocated (Baumgartner et al.,
2020).

Intensity interferometry measures the squared visibility, or equivalently the
2D angular power spectrum of the source brightness, and from the observable it
is possible to reconstruct the image of two stars (Dravins et al., 2013). In
this work we adopt a different strategy, which is to take a parameterised
model for two limb-darkened stars and fit their radii and separation to the
interferometric observable.

This work discusses parameter estimation using simulated data for a binary. A
Spica-like system is taken as an example, but no special properties of Spica
are assumed. We consider two methods to recover the system parameters from the
squared visibility. The first involves comparing the observable at two (or
more) wavelengths to separate the visibility contributions of the two stars.
Essential for this is that the two stars have different effective
temperatures, and hence contribute differently at different wavelengths. The
second method involves direct fitting of the squared visibility to a two-star
model with parameters for the stellar radii, the projected distance, and limb
darkening.

This paper is organized as follows. In Section 2, we describe the methods of
separation in visibilities and direct fitting to estimate parameters, assuming
a uniform-brightness model of the star. However, in reality, stars do not
follow a uniform model of brightness. So in Section 3 we describe a limb
darkening model of the stars. This is followed in Section 4 by simulations of
intensity-correlation data and the recovery of the parameters of interest
using a Markov chain Monte Carlo (MCMC) technique. The results our analysis
are described in Section 5. Finally we conclude in Section 6.

##  2 Methodology

An amplitude interferometer measures the modulus and phase of the Fourier
transform of the brightness distribution of an incoherent source. Hence the
response of the amplitude interferometer, if sufficiently well sampled, can be
inverted to obtain the original image of the source. The processing of an
image obtained in this way, is heavily dependent on the ability of the
measuring device to accurately track the phase of the light signal coming from
sources. The latter is not required in intensity interferometry, which
measures only the modulus of the Fourier transform of a source’s brightness
distribution. Because of the lack of phase information, reconstruction of a
source having a complex structure requires information to be added, as a model
or other prior. In the case where the source’s shape is known, a precise size
estimate can be extracted in a model-dependent manner.

###  2.1 Correlated Signal

We consider first a binary source comprising two stars $A$ and $B$ of uniform
brightness with angular radii $\Theta_{A}$ and $\Theta_{B}$, and effective
temperature $T_{A}$ and $T_{B}$ respectively, at angular distance $\Theta_{d}$
from each other. The total flux of photons coming from the binary source is
the weighted sum of the photon flux coming from the individual stars

$\displaystyle\Phi=I_{A}+I_{B}=|S_{A}(\Omega)|^{2}\pi\Theta_{A}^{2}+|S_{B}(\Omega)|^{2}\pi\Theta_{B}^{2}$
(1)

where $I_{A}=|S_{A}(\Omega)|^{2}$ and $I_{B}=|S_{B}(\Omega)|^{2}$ are the flux
of photons from stars $A$ and $B$ respectively.

For a star of uniform brightness, visibility function at a point $(u,v)$ on
the interferometric plane has the well-known form

$\displaystyle
V_{A}(u,v)\propto\frac{J_{1}(\rho\,\Theta_{A})}{(\rho\,\Theta_{A})}$ (2)

and analogously for $V_{B}(u,v)$, where

$\displaystyle\rho=\sqrt{u^{2}+v^{2}}$ (3)

For both stars together the corresponding visibility function is

$\displaystyle
V(u,v)=\frac{\displaystyle{I_{A}V_{A}(u,v)+I_{B}V_{B}(u,v)\exp\left(i\rho\,\Theta_{d}\cos\psi\right)}}{I_{A}+I_{B}}$
(4)

where the $u$ axis of the baseline is oriented along the line connecting the
sources. The square of the modulus of the visibility function

$|V(u,v)|^{2}=\frac{\Biggr{[}\begin{multlined}I_{A}^{2}V_{A}^{2}(u,v)+I_{B}^{2}V_{B}^{2}(u,v)\\\
+2I_{A}I_{B}V_{A}(u,v)V_{B}(u,v)\cos{\left(\rho\,\Theta_{d}\cos\psi\right)}\end{multlined}I_{A}^{2}V_{A}^{2}(u,v)+I_{B}^{2}V_{B}^{2}(u,v)\\\
+2I_{A}I_{B}V_{A}(u,v)V_{B}(u,v)\cos{\left(\rho\,\Theta_{d}\cos\psi\right)}\Biggr{]}}{\left[I_{A}+I_{B}\right]^{2}}$
(5)

is the main observable in intensity interferometry. By construction
$V(0,0)=1$. In practice, however, the measured quantity is

$\displaystyle g(u,v)=\frac{\Delta\tau}{\Delta t}|V(u,v)|^{2}$ (6)

where $\Delta\tau$ is the coherence time or equivalently the reciprocal of the
frequency bandwidth (Mandel & Wolf, 1965) while $\Delta t$ is the time
resolution.

###  2.2 Methods to fit parameters

The discussion of this paper is based on two possible strategies, which we now
desribe, to estimate the parameters of the stars.

####  2.2.1 Separation in visibilities

In this method, we consider measurements at two different wavelengths.
Provided the stars have different effective temperatures, the squared
visibility will be measured with two different values of the flux ratio
$f=I_{B}/I_{A}$. Then there is a possibility to estimate both radii of the
system.

Along the $v$ axis, i.e., $u=0$ the visibility (4) will be real. We may then
take the square root of the squared visibility (5) to obtain

$\displaystyle V_{1}(0,v)=\frac{V_{A}(0,v)+f_{1}V_{B}(0,v)}{1+f_{1}}$ (7)

and

$\displaystyle V_{2}(0,v)=\frac{V_{A}(0,v)+f_{2}V_{B}(0,v)}{1+f_{2}}$ (8)

Assuming the sign of $V_{1}(0,v)$ and $V_{2}(0,v)$ can be inferred, the
visibility contributions of source $A$ and source $B$ in terms of $V_{1}(0,v)$
and $V_{2}(0,v)$ are

$\displaystyle
V_{A}(0,v)=\frac{\Biggr{[}\begin{multlined}V_{1}(0,v)f_{1}f_{2}-V_{2}(0,v)f_{1}f_{2}\\\
+V_{1}(0,v)f_{2}-V_{2}(0,v)f_{1}\end{multlined}V_{1}(0,v)f_{1}f_{2}-V_{2}(0,v)f_{1}f_{2}\\\
+V_{1}(0,v)f_{2}-V_{2}(0,v)f_{1}\Biggr{]}}{f_{2}-f_{1}}$ (11)

and

$\displaystyle
V_{B}(0,v)=\frac{\Biggr{[}\begin{multlined}f_{1}V_{1}(0,v)-f_{2}V_{2}(0,v)\\\
+V_{1}(0,v)-V_{2}(0,v)\end{multlined}f_{1}V_{1}(0,v)-f_{2}V_{2}(0,v)\\\
+V_{1}(0,v)-V_{2}(0,v)\Biggr{]}}{f_{1}-f_{2}}$ (14)

The expressions (11) and (14) can now be used to fit parameters, specifically
the radii of both stars. But using this method, the distance between stars can
not estimated as the above equation are independent of $\Theta_{d}$.

####  2.2.2 Direct fitting

The second method we consider is directly fitting the squared visibility (5)
to infer the parameters. Since $|V(u,v)|^{2}$ depends explicitly on
$\Theta_{d}$ as well as on $\Theta_{A}$ and $\Theta_{B}$, we can use this
correlation to estimate the angular separation between the sources together
with angular radii of the sources.

###  2.3 Slope of the Signal

In order to better understand the dependence of the squared visibility on each
of the radii in different parts of the $u,v$ plane, we consider derivatives of
the former (equ. 5) with respect to $\Theta_{A}$ and $\Theta_{B}$ for $u=0$,
which are as follows.

$\displaystyle\frac{\partial|V(0,v)|^{2}}{\partial\Theta_{A}}=-\frac{4I_{A}V(0,v)}{(I_{A}+I_{B})^{2}}\frac{J_{2}\left(v\Theta_{A}\right)}{\Theta_{A}}$
(15)

and

$\displaystyle\frac{\partial|V(0,v)|^{2}}{\partial\Theta_{B}}=-\frac{4I_{B}V(0,v)}{(I_{A}+I_{B})^{2}}\frac{J_{2}\left(v\Theta_{B}\right)}{\Theta_{B}}$
(16)

The first factor of eqns. (15) and (16) depends on both sources. However, the
second factor depends on only the source with respect to which equ. (5) is
differentiated. As we can see, the derivative is a product of Bessel functions
$J_{1}$ and $J_{2}$, and can be positive or negative. So, there will be some
regions in the $u,v$ plane where the derivative of the observable will be
positive and some regions where it will be negative. Using these equations, we
can remove the data from the baseline region that reduces the strength of the
signal (i.e., the area of baselines from where signals with negative slope
observed).

##  3 Impact of Limb-Darkening

Limb darkening is the phenomenon that describes the brightness of a star as a
function of its radius, observed relative to the maximum brightness at the
center of the stellar disk. The disk of the star is dimmer towards the limbs
than at the center. As first explained by Schwarzschild (1906) limb darkening
occurs due to the gradient of temperatures through the layers of the star.
Higher temperatures emit more photons at a given wavelength, and lower
temperatures emit fewer. The photons from a star are emitted from one optical
depth within the surface of the star. Because of the star’s geometry, one
optical depth at the center of the stellar disk penetrates through to a deeper
layer of the star than one optical depth near the limb; therefore, we see a
hotter region, and thus middle of the stellar disk appears brighter.

Limb darkening is typically treated as a simple parameterised function of the
angle between a line normal to the stellar surface and the observer’s line of
sight, which makes fitting the stellar intensity profile much simpler and
reduces the number of free parameters. The most common parametrizations are
linear and quadratic relations (Kopal, 1950), but other suggested relations
include a four-parameter relation, a square-root relation as well as
exponential and logarithmic relations (Milne, 1921; Diaz-Cordoves & Gimenez,
1992; Klinglesmith & Sobieski, 1970; Claret & Hauschildt, 2003; Sing et al.,
2009; Claret, 2000). Following the work described in Hanbury Brown et al.
(1974) we have used a limb darkening model of the form

$\displaystyle\left|\frac{S(\Omega,\mu)}{S(\Omega,1)}\right|^{2}=1-l(1-\mu)-m(1-\mu)^{2}-n(1-\mu)^{3}$
(17)

where $l$, $m$ and $n$ are limb coefficient depending on the surface
temperature and frequency, and $\mu$ is the cosine of the angle $\theta$
between the normal to the surface at that point and the line of sight from the
star to the observer. In terms of the radius $R$ of the star, the projected
radial position $R^{\prime}$ of the observed point, and the angle $\theta$ is
expressed as

$\displaystyle\mu=\cos\theta=\sqrt{1-\frac{R^{\prime 2}}{R^{2}}}$ (18)

Given the limb darkening model under consideration and the under the
assumptions that the values of the limb-darkening coefficients are frequency
independent and same for both the stars, the total flux of binary source is
given by

$\displaystyle\Phi=C\left(I_{A}+I_{B}\right)$ (19)

where the constant $C$ in terms of limb coefficients is given by,

$\displaystyle C=\left(1-\frac{l}{3}-\frac{m}{6}-\frac{n}{10}\right)$ (20)

and the visibility functions of star $A$ and $B$ are given by

$\displaystyle
V_{A}(u,v)=\sum_{i}a_{i}\frac{J_{i}(\rho\,\Theta_{A})}{(\rho\,\Theta_{A})^{i}}$
(21) $\displaystyle
V_{B}(u,v)=\sum_{i}a_{i}\frac{J_{i}(\rho\,\Theta_{B})}{(\rho\,\Theta_{B})^{i}}$
(22)

For $i=1,~{}3/2,~{}2,~{}5/2$ and

$\displaystyle\begin{aligned} &a_{1}=1-l-m-n\\\
&a_{3/2}=\sqrt{\frac{\pi}{2}}(l+2m+3n)\\\ &a_{2}=-2(m+3n)\\\
&a_{5/2}=3\sqrt{\frac{\pi}{2}}n\end{aligned}$ (23)

In the absence of star $B$, the result will be according to Hanbury Brown et
al. (1974).

Eqns. (21) and (22) are the generalized form of visibilities of sources $A$
and $B$, according to the coefficient parameters $l,m,n$ the limb darkening
model of star assumed here.

##  4 Simulations and parameter fitting

We now proceed to simulate intensity-correlation data and recover the
parameters using a Markov chain Monte Carlo (MCMC) algorithm.

###  4.1 Signal to noise

We may write the signal as

$\displaystyle g(u,v,\omega)=\frac{\Delta\tau}{\Delta t}\,|V(u,v,\omega)|^{2}$
(24)

which is equ. (6) slightly rewritten with $\omega$ standing for all the
parameters of the system (radii, distance, and limb-darkening coefficients).
For the signal at the $i$-th baseline, we write

$\displaystyle g_{i}(u_{i},v_{i},\omega)=\frac{\Delta\tau}{\Delta
t}\,|V(u_{i},v_{i},\omega)|^{2}$ (25)

and the data collected at that baseline will be

$\displaystyle d_{i}=g_{i}+\sigma_{i}{\cal N}_{i}$ (26)

where $\sigma_{i}$ is the noise level and ${\cal N}_{i}$ is a Gaussian random
number.111The noise level depends on the observing time at each baseline, and
in the simulations we assume all the $\sigma_{i}$ to the same for all
baselines, though this is not essential. The signal-to-noise at the $i$-th
baseline will be

$\displaystyle{\rm SNR}_{i}=\frac{d_{i}}{\sigma_{i}}$ (27)

while

$\displaystyle{\rm SNR}=\Big{(}{\textstyle\sum_{i}}\,{\rm
SNR}_{i}^{2}\Big{)}^{1/2}$ (28)

will be the total signal-to-noise ratio for a full simulation.

For the visibility-separation method, the relation between ${\rm SNR}_{i}$ and
$\sigma_{i}$ is more complicated. First, the observable squared visibility has
to be converted to individual visibilities using eqns. (7) and (8).
Correspondingly, the ${\rm SNR}_{|V|^{2}}$ (say) relating to the observable,
has to be converted to an effective

$\displaystyle{\rm SNR}_{V}=2{\rm SNR}_{|V|^{2}}$ (29)

using the error propagation formula. Further using the latter formula we have

$\displaystyle\sigma_{V_{A}}^{2}=\left(\frac{\partial V_{A}}{\partial
V_{1}}\right)^{2}\sigma_{V_{1}}^{2}+\left(\frac{\partial V_{A}}{\partial
V_{2}}\right)^{2}\sigma_{V_{2}}^{2}$ (30)

and similarly $\sigma_{V_{B}}^{2}$. The values of $\sigma_{V_{A},i}$ and
$\sigma_{V_{B},i}$ are required for parameter-fitting.

To relate the SNR in a simulation to the expected SNR in an observing run, we
recall that for one baseline

$\displaystyle{\rm
SNR}_{i}\approx\left(A_{1}A_{2}\right)^{1/2}\,\Phi\,|V(u_{i},v_{i})|^{2}\left(t^{\rm
obs}_{i}/\Delta t\right)^{1/2}$ (31)

where $A_{1}$ and $A_{2}$ are the effective collecting areas,

$\displaystyle\Phi\approx\frac{\lambda}{1\mu\rm m}\times 10^{-4-{\rm
AB}/2.5}\;\hbox{photons}\;{\rm m^{-2}\,s^{-1}\,Hz^{-1}}$ (32)

is the spectral photon density, and and $t^{\rm obs}_{i}$ is the observing
time at that baseline. The full SNR of the observing run is then obtained
using equ. (28) as before.

###  4.2 Parameter fitting

For a given data set the likelihood in the usual Bayesian formulation will be

$\displaystyle\ln
P\left(\\{d_{i}\\}|\alpha,\omega\right)=-{\textstyle\frac{1}{2}}\sum_{i}\sigma_{i}^{-2}\Big{(}d_{i}-\alpha|V(u_{i},v_{i},\omega)|^{2}\Big{)}^{2}$
(33)

where $\alpha=\Delta\tau/\Delta t$ is an additional parameter, which will not
be accurately known in advance. We can marginalise out the nuisance parameter
$\alpha$ to avoid having to fit it, as follows. Assuming a flat prior for
$\alpha$ we have

$\displaystyle P\left(\\{d_{i}\\}|\omega\right)=\int
P\left(\\{d_{i}\\}|\alpha,\omega\right)d\alpha$ (34)

We now define two scalar products: first of the data with the model

$\displaystyle G=\sum_{i}\sigma_{i}^{-2}\,d_{i}\,|V(u_{i},v_{i},\omega)|^{2}$
(35)

and, the scalar product of model with itself

$\displaystyle W=\sum_{i}\sigma_{i}^{-2}\,|V(u_{i},v_{i},\omega)|^{4}$ (36)

The integral (33) can now be expressed as

$\displaystyle\ln
P\left(\\{d_{i}\\}|\omega\right)=\frac{G^{2}}{2W}-\frac{1}{2}\ln W$ (37)

Eqns. (35–37) are convenient for input to MCMC. An MCMC implementation, such
as the well-known emcee code (Foreman-Mackey et al., 2013) which we used,
generates a sample of points $\omega$ in parameter space sampled according to
some given function of $\omega$. Sampling is desired according to
$P\left(\\{d_{i}\\}|\omega\right)$ times the prior probability of $\omega$,
the latter representing the allowed ranges of the various parameters. The
resulting sample of $\omega$ points cannot in general be visualised directly,
and hence projections to one- and two-parameter subspaces are conventionally
used, as in Fig. 5.

##  5 Results

In this section we first describe the signal of our interest i.e., visibility
and squared visibility. This is followed by the estimation of parameters in
two ways: first by separation of visibilities, and second the direct fitting
of the squared visibility.

The parameter values used for the simulations are
$\Theta_{A}=2.20,\Theta_{B}=1.10$ nanoradians for the angular radii and
$\Theta_{d}=8.30$ nanoradians for the angular separation.222Recall that a mas
is about 5 nanoradians. These values correspond roughly to the Spica system
Herbison-Evans et al. (1971); Tkachenko et al. (2016). Simulated observations
at $\nu=$5\text{\times}{10}^{14}\text{\,}\mathrm{Hz}$$ and
$\nu=$7\text{\times}{10}^{14}\text{\,}\mathrm{Hz}$$ are considered. For
brevity we refer to these as red and blue. The ratio $\Delta\tau/\Delta t$ of
coherence time to resolution time is taken to be $0.11$.

###  5.1 Signal

Fig. 1 shows the main signal for baselines perpendicular to the orientation of
the binary. As expected, the signal attains its maximum value for zero
baseline, and for a single star the larger the source size the narrower is the
visibility pattern. However, the visibility of the total system is weighted
towards the brighter source. Source A dominates the brightness both in blue
and red, but is slightly less dominant in red; this difference is, however,
not discernible from the figure. For certain wavelength (most noticibly around
$175$m at blue wavelength) the visibility $V(0,v)$ is found to be less for the
binary system than those for either of the sources under consideration. This
is due to the fact that the visibilities for sources have opposite sign and
are interfering destructively.

(a) Visibility along the baseline

(b) Squared visibility along baseline

(c) Visibility along unit-less baseline

(d) Squared visibility along unit-less baseline

Figure 1: Analytical plot of the signal at two different wavelengths, for
baselines perpendicular to the line between the stars. The two left panels
show the visibility, while the two right panels show the squared visibility.
For the two upper panels the baseline is in metres, while the two lower panels
use the dimensionless baseline $v$. According to eqns. (4) and (5) the signal
is always unity at zero baseline.

(a) SNR = 32700

(b) SNR =650

(c) SNR =320

Figure 2: Visibility separation from simulated data illustrated at different
SNR levels. The solid red and dotted black line show the combined visibility
with $f_{1}=0.1$ and $f_{2}=0.11$ respectively (see eqns. 7 and 8). The
baseline $v$ is perpendicular to the line between the two sources. The dashed
blue line and the dotted green line show the inferred individual visibilities
of source $A$ and source $B$ respectively.

###  5.2 Measuring Parameters

For visibility separation, the values of the parameters $f_{1}$ and $f_{2}$
were assumed to be $0.10$ and $0.11$ respectively. The separation has been
done for different values of SNR as shown in Fig. 2. The resulting estimates
of the angular radii are shown in Fig. 3. Our results show that the minimum
value of SNR for which this method works is around $300$.

(a) Corner plot for SNR = 32700

(b) Corner plot for SNR=650

(c) Corner plot for SNR=320

Figure 3: Estimation of the angular radii (in nanorad) of stars $A$ and $B$
using separated visibilities shown in Fig. 2.

(a) Model data along baseline

(b) Simulated data along baseline

Figure 4: Observable intensity correlation $g(u,v)$ at differently-oriented
baselines and different observing wavelengths. The labelled
$\cos\psi=u/\sqrt{u^{2}+v^{2}}$ where $u$ is along the line connecting the two
stars. The left panel shows the model according to eqns. (5) and (6). The
right panel shows the signal with added noise corresponding to a total SNR of
370.  Figure 5: Estimation of the angular radii of the two stars and angular
separation between them using the direct fitting method, from the simulated
data shown in Fig. 4(b). Note the positive correlation between angular radii
of the two sources. Figure 6: Derivatives of the squared visibility with
respect to angular radii of the stars (eqns. (15) and (16)), shown along
baselines perpendicular to the line connecting the stars. The regions with a
sharp negative slope in the signal are considered for exclusion in the
following simulations.

(a) Model data along baseline

(b) Simulation data along baseline

Figure 7: Observable HBT correlation $g(u,v)$ after removing the region of
baseline identified in Fig. 6. As Fig. 4 the left panel is for the model
signal and the right panel adds noise. The total SNR here is 208.

For direct fitting of the squared visibility, we simulated the HBT correlation
(Fig. 4(a)) along the baseline for $\cos\psi=0,0.5,1$. It is not necessary to
have the signal along straight lines in the $(u,v)$ plane, but this choice
conveniently illustrates the oscillatory features in the correlation, which
are due to the presence of a cosine factor in the expression of correlation
equ. (5). The added noise is such that the total ${\rm SNR}=370$. Fig. 5
clearly shows good agreement between input and recovered values of the
parameters.

Although Fig. 5 shows no correlation between angular separation and angular
radii of either star and, however, it shows positive correlation between the
angular radii of the stars. In order to explain the origin of the correlation
between the angular radii of the stars, we have plotted in Fig. 6 the change
in squared visibility with angular radii of sources along baseline for the
orientation of the sources in the direction perpendicular to the baseline. The
sharp negative slope in Fig. 6 appears to be the cause of the correlation
between estimated angular radii of the stars. In order to resolve this issue
we tried removing from our analysis the unwanted portion of the baseline and
followed the same methodology to estimate the parameters from the unmasked
portion of the baseline, as shown in Fig. 7. The results in Fig. 8 show no
correlation between the angular radii of the stars.

Figure 8: Estimation of the angular radii of the two stars and angular
separation between them using the direct fitting method, from the simulated
data shown in Fig. 7(b). The difference from Fig. 5 is that some baselines
regions have been excluded, with the result that the positive correlation
between inferred angular radii of both sources has been mostly removed.

(a) SNR = 10346

(b) SNR = 79

(c) SNR = 55

(d) SNR = 41

(e) SNR = 35

(f) SNR = 32

Figure 9: Estimations of the angular radii of both sources and the angular
separation between them, using direct fitting at different total SNR. The
threshold for estimating all three parameters is total SNR of about 30.

###  5.3 Threshold Signal-To-Noise Ratio

Estimation of parameters works requires a minimum value of the signal-to-noise
ratio. Fig. 9 checks this threshold value of the SNR for estimation of angular
radii of both sources and angular separation between them using a direct
fitting with masking unwanted region. The minimum value of SNR is around 30.
If the total SNR is much lower than in Fig. 9(f), parameter estimation becomes
poor, with the uncertainties becoming large and the credible regions excluding
the true value.

(a) limb coefficients $l=0.6,m=0.0,n=0.0$

(b) limb coefficients $l=0.6,m=0.2,n=0.0$

(c) limb coefficients $l=0.6,m=0.2,n=0.1$

Figure 10: Analytical plots of the squared visibility function, comparing a
uniform-brightness model and limb-darkened models. The effect of limb
darkening is measurable only near the peaks of the signal curve. The three
panels show the differences between linear, quadratic and cubic limb-
darkening, which are very small.

###  5.4 Effect of limb-darkening

In order to show the impact of limb-darkening on the squared visibility
function of the stars visibilities for the orientation of the baseline
corresponds to $\psi=0$. Fig. 10 shows that the impact of limb-darkening for
three models: linear, quadratic and cubic. In all cases the difference due to
limb darkening is small, and more prominent around the peaks. Investiging
further we have found that the effect is relatively more noticeable for the
baseline range $(0.5-0.8)\times 10^{9}$ and $(1.1-1.5)\times 10^{9}$. It is to
be noted that for the baseline $1.845\times 10^{9}$, the squared visibility is
always more for the limb darkened model than the uniform model. The reverse is
the situation for the baseline higher than $1.845\times 10^{9}$.

Fig. 11(a) shows the measurement of the parameters of the stars with limb-
darkening taken into account. The input values of the limb-darkening
coefficients used to generate the plots are $l=0.6,m=0.2$, and $n=0.1$ and
signal-to-noise-ratio equal to $370$. For estimating the parameters, angular
radii of both stars, and angular separation between them, limb coefficients
are taken as constant. From Fig. 10 and Fig. 11 it is clear that the effect of
limb-darkening of stars can be measurable and estimation of limb-coefficient
is also possible. Fig. 11 shows the measurement of all the limb-darkening
coefficients together with the parameters of the stars. For the linear model
of the limb-darkening of stars, we have good recovery of the limb darkening
coefficients $l$. On the other hand, for complex limb darkening models the
recovery of the nonlinear limb-darkening coefficients $m$ and $n$ is not
unsuccessful (Fig. 11(c) and Fig. 11(d)). However, in all cases, we have
managed to recover both angular radii and angular separation between stars
also linear limb-darkening coefficient with good accuracy.

(a) SNR = 370

(b) SNR = 1850

(c) SNR = 1850

(d) SNR = 1850

Figure 11: Estimation of various numbers of parameters. The upper left panel
estimates only estimates only the angular radii and the angular separation of
the stars, with limb-darkening coefficients $l=0.6,m=0.2,n=0.1$ given. In the
upper-right panel, $l$ is estimated while $m$ and $n$ are given. The lower-
left corner plot estimates all parameters except $n$ and the lower-right
corner plot is for all parameter being estimated. The estimation of $m$ and
$n$ is not successful because of their small contribution.

##  6 Discussion

Recent years have seen a resurgence of interest in intensity interferometry,
and many science cases have been discussed in the literature (Kieda et al.,
2019). In this paper we argue that radius measurements of binary stars are an
interesting science case for intensity-interferometry observatories now in
development. The basic idea was known during the earlier generation of
intensity interferometry: Herbison-Evans et al. (1971) observed Spica using
intensity interferometry; combining the interferometric observations with
spectroscopy, they reconstructed the binary orbit and inferred the distance to
the system, weighed both the stars, and measured the radius of the brighter
star; the radius of the second star was the main parameter that could not be
determined. It is plausible that the signal-to-noise needed to measure the
radius of the second star was available even in 1971, but interferometric
image reconstruction was still in its infancy. The interferometric pattern of
a binary looks quite similar to a single star modulated with fringes (equ.
(5)) and significant computing is needed to disentangle the contributions of
two stars.

Nowadays, several image-reconstruction algorithms and codes are available
(Berger et al., 2012) including developments specific to intensity
interferometry (Nuñez & Domiciano de Souza, 2015). In particular, the
reconstruction of an image of two stars from simulated data has been
demonstrated (e.g., Figures 4–6 in Dravins et al., 2013). In this work,
however, we suggest that rather than a general-purpose reconstruction of an
image of two stars, parameter-fitting of the observable directed towards
inferring the system parameters may be useful. Parameters that are not of
interest but are unavoidable, such as the normalisation of the visibility, can
be marginalised out. Markov chain Monte Carlo is well suited to this task, and
efficient implementations are readily available. We used the emcee package.

In this paper we carry out a number of simulations, showing that with
sufficient SNR, it is possible to extract both stellar radii, the angular
separation between the stars, and linear limb-darkening. The SNR of the
simulations can be related to source properties, telescope parameters, and
observing time (eqns. 31 and 32) using well-known expressions for signal and
noise in intensity interferometry. For example, suppose the spectral
brightness times the effective area (eqns. 31 and 32) is
$10^{-4}\rm\,s^{-1}\,Hz^{-1}$. This will be the case for zero-magnitude source
and an effective collecting around of about $1\rm\,m^{2}$, or a 5 mag source
and an effective collecting around of about $100\rm\,m^{2}$, and so on. With
nanosec time resolution and an observing time of $10^{5}\rm\,s$ (2–3 nights)
the total SNR will be $\sim 1000$.

We consider two possible approaches to parameter-fitting for a binary-star
system.

If the stars have different effective temperatures and hence differently-
shaped spectra, and are observed at least two wavelengths, it is possible to
separate the visibility contributions of each stars. There are several
complications to the idea. The observable has to be square-rooted while taking
into account possible sign changes. Moreover, the observations would need to
be perpendicular to the line joining the stars, where the visibility is real.
Then, visibility contributions have to be compared at the same values of
dimensionless baseline at two different wavelengths, which will correspond to
different lengths of the physical baseline. Finally, the visibility separation
tends to amplify the noise. Nevertheless, the simulations indicate that
measurement of both stellar radii is possible by this method.

Our second approach is multi-parameter fitting to the observed HBT correlation
(squared visibility). Provided the parametric form being fitted is correct,
one expects higher accuracy in the parameters. We consider a model with two
discs in the sky with limb-darkening up to third order, at some separation
from each other. Rough prior estimates on the radii of the two discs and their
separation are assumed. A limb-darkening model with linear, quadratic, and
cubic terms is used. If the limb-darkening parameters are assumed known, a
total SNR of 350 suffices to estimate the angular radii of the stars and
separation between the stars with reasonable accuracy. A higher total SNR of
about 1800 is required to estimate the linear limb-darkening coefficient. The
higher-order coefficients could not be recovered at this SNR, but that fact
does not degrade the esitmates of the other parameters.

The precision of parameter recovery improves with increasing SNR, but it also
has a strong dependence on the sampling of the $(u,v)$ plane. That is to say,
even if two regions of the $(u,v)$ plane may contribute equally to the total
SNR for given observing time, one of the regions may be much more useful for
parameter estimation. We suggest that the derivatives of the visibility with
respect to the parameters can help identify the most important parts of the
interferometric plane. By concentrating on parts of $(u,v)$ plane so
identified, we found that a total SNR of 30 was sufficient to recover the
angular radii and separation.

## References

  * Abeysekara et al. (2020) Abeysekara A. U., et al., 2020, Nature Astronomy
  * Acciari et al. (2020) Acciari V. A., et al., 2020, MNRAS, 491, 1540
  * Baumgartner et al. (2020) Baumgartner S., et al., 2020, MNRAS, 498, 4577
  * Berger et al. (2012) Berger J. P., et al., 2012, A&ARv, 20, 53
  * Claret (2000) Claret A., 2000, A&A, 363, 1081
  * Claret & Hauschildt (2003) Claret A., Hauschildt P. H., 2003, A&A, 412, 241
  * Diaz-Cordoves & Gimenez (1992) Diaz-Cordoves J., Gimenez A., 1992, A&A, 259, 227
  * Docobo et al. (2018) Docobo J. A., Campo P. P., Gomez J., Horch E. P., 2018, AJ, 156, 185
  * Dravins et al. (2013) Dravins D., LeBohec S., Jensen H., Nuñez P. D., 2013, Astroparticle Physics, 43, 331
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Hanbury Brown et al. (1974) Hanbury Brown R., Davis J., Lake R. J. W., Thompson R. J., 1974, Monthly Notices of the Royal Astronomical Society, 167, 475
  * Herbison-Evans et al. (1971) Herbison-Evans D., Hanbury Brown R., Davis J., Allen L. R., 1971, MNRAS, 151, 161
  * Kieda et al. (2019) Kieda D., et al., 2019, BAAS, 51, 275
  * Kjeldsen & Bedding (1995) Kjeldsen H., Bedding T. R., 1995, A&A, 293, 87
  * Klinglesmith & Sobieski (1970) Klinglesmith D. A., Sobieski S., 1970, AJ, 75, 175
  * Kopal (1950) Kopal Z., 1950, Harvard College Observatory Circular, 454, 1
  * Mandel & Wolf (1965) Mandel L., Wolf E., 1965, Reviews of modern physics, 37, 231
  * Milne (1921) Milne E. A., 1921, Monthly Notices of the Royal Astronomical Society, 81, 361
  * Nuñez & Domiciano de Souza (2015) Nuñez P. D., Domiciano de Souza A., 2015, MNRAS, 453, 1999
  * Rivet et al. (2020) Rivet J. P., et al., 2020, MNRAS, 494, 218
  * Rodrigues et al. (2017) Rodrigues T. S., et al., 2017, MNRAS, 467, 1433
  * Russell (1928) Russell H. N., 1928, AJ, 38, 89
  * Schwarzschild (1906) Schwarzschild K., 1906, Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse, 1906, 41
  * Sing et al. (2009) Sing D. K., Désert J. M., Lecavelier Des Etangs A., Ballester G. E., Vidal-Madjar A., Parmentier V., Hebrard G., Henry G. W., 2009, A&A, 505, 891
  * Stephenson & Sanwal (1969) Stephenson C. B., Sanwal N. B., 1969, AJ, 74, 689
  * Tkachenko et al. (2016) Tkachenko A., et al., 2016, MNRAS, 458, 1964