FM stars II: a Fourier view of pulsating binary stars -- determining binary orbital parameters photometrically for highly eccentric cases

Continuous and precise space-based photometry has made it possible to measure the orbital frequency modulation of pulsating stars in binary systems with extremely high precision over long time spans. Frequency modulation caused by binary orbital motion manifests itself as a multiplet with equal spacing of the orbital frequency in the Fourier transform. The amplitudes and phases of the peaks in these multiplets reflect the orbital properties, hence the orbital parameters can be extracted by analysing such precise photometric data alone. We derive analytically the theoretical relations between the multiplet properties and the orbital parameters, and present a method for determining these parameters, including the eccentricity and the argument of periapsis, from a quintuplet or a higher order multiplet. This is achievable with the photometry alone, without spectroscopic radial velocity measurements. We apply this method to Kepler mission data of KIC8264492, KIC9651065, and KIC10990452, each of which is shown to have an eccentricity exceeding 0.5. Radial velocity curves are also derived from the Kepler photometric data. We demonstrate that the results are in good agreement with those obtained by another technique based on the analysis of the pulsation phases.


INTRODUCTION
The study of binary and multiple stars provides the best measurements of stellar parameters, giving direct, nearly modelindependent knowledge of stellar masses, and radii as well from eclipsing binaries. These are fundamental to our understanding of stellar structure and evolution, the bedrock upon which much of stellar astrophysics rests. Asteroseismology can now also provide stellar masses and radii, and even ages for individual stars, but our confidence in asteroseismology rests on calibration of asteroseismic masses and radii with those determined from binary stars.
The primary observational data for binary star studies are the light curve and radial velocity curve. Additional constraints are given by measurement of stellar effective temperature, T eff , from spectroscopy, and, for a few stars, stellar radius from interferometry. In the case of T eff and interferometric radius, single measurements are sufficient. The light curve and radial velocity curve, however, are by definition time series, and it is these that demand so much time and effort from observers.
Traditionally, light curves have been measured with groundbased photometry from a single observing site in one or more pho-tometric bandpasses. Radial velocity curves have been measured by obtaining spectra from a single site repeatedly until the orbital phase is well covered. So, a large amount of observing time is required, often on large telescopes for fainter stars. This is the traditional bottleneck to binary star studies, and until recently the result has been that high precision masses and radii for stars were only known for a matter of hundreds of stars. That all of stellar structure and evolution theory, and beyond that galactic astronomy, should be based on such a sparse base has been a significant concern. Now the study of stellar light curves has been revolutionised. Whereas until recently light curves were studied for no more than a few stars from a few sites, for a limited number of nights per year with great gaps in the data, we now have nearly continuous spacebased photometry for hundreds of thousands of stars. Whereas from the ground photometric measurements were precise to parts per thousand, or in the best cases, about 10 times better than that, the space data give precision as good as parts per million. The satellite telescope photometry revolution has a pedigree from the Canadian MOST mission 1 , through the ESA CoRoT mission 2 to the NASA Kepler mission 3 , and it anticipates the upcoming NASA mission TESS 4 and ESA mission PLATO 5 .
It is the Kepler mission that has led to the greatest results, with -as of this writing -over 4600 exoplanet candidate discoveries 6,7,8 , over 1800 confirmed exoplanets, asteroseismic results for thousands of pulsating stars, and with light curves of more than 2700 eclipsing binary stars 9 . The main Kepler mission data set provides 4 yr of nearly continuous photometric white light data with 30-min time resolution for 150 000 stars, with further data for shorter time spans for another 50 000 stars. The follow-up studies for exoplanets involve many hundreds of astronomers and require large amounts of ground-based telescope time -often on some of the world's largest telescopes -for radial velocity measurements necessary to confirm that masses are indeed planetary. The thousands of binary stars also require large amounts of ground-based telescope time to acquire the radial velocities needed for full astrophysical study of the multiple systems.
The impediment in these studies is the acquisition of the spectra that have traditionally been used to measure radial velocity from the Doppler shift of wavelength of spectral lines. However, Doppler shift can be measured for any astrophysical source of a stable frequency. We have developed two methods to do this directly from Kepler mission photometric light curves for pulsating stars, without recourse to ground-based observations, where the pulsation frequencies are the standard, rather than spectroscopic wavelength. These methods obviate the need for any ground-based observations to determine radial velocity, and they provide nearly continuous radial velocity curves that also cover the entire Kepler 4-yr observational time span -an unreachable goal with ground-based radial velocity measurements. Our methods are limited to stars with stable pulsation frequencies, but there are thousands of these in the Kepler data with amplitudes high enough to apply our techniques.
Our two methods are the frequency modulation (FM) method (Shibahashi & Kurtz 2012) and the phase modulation method (Murphy et al. 2014; see also Balona 2014;Koen 2014). Our techniques are related to older observed-minus-calculated (O − C) methods of studying stellar frequency variability, but the FM and PM techniques both have the advantage that they use the entire data set for maximum signal-to-noise ratio and that they are particularly suited to multi-mode pulsators. The FM technique has the highest possible frequency resolution, and the PM technique is both easily automated and can combine the results from many pulsation frequencies easily. In addition, many studies of pulsating stars for asteroseismic inference begin with the frequency spectrum, and the FM technique shows explicitly the patterns expected in an amplitude spectrum, or power spectrum, for pulsation frequencies undergoing periodic Doppler shifts from binary orbital motion.
We have demonstrated the validity of the FM method by showing the consistent results obtained from it when compared to a traditional eclipsing binary light curve analysis (Kurtz et al. 2015).
Furthermore, the FM technique is also applicable to non-eclipsing systems and is a powerful tool to determine the orbital parameters.
In practice, the information desired from the radial velocity curve includes the orbital period, the mass function, f (m1, m2, sin i), the orbital semi-major axes, a1 sin i and a2 sin i, the eccentricity, e, the angle between the nodal point and the periapsis, ̟, and the time of periapsis passage, tp. Shibahashi & Kurtz (2012) showed how to derive all of these except e and ̟ from photometric data using the FM method, while Murphy et al. (2014) showed how to derive all of them from the PM method. In this paper we show how to extract e and ̟ with the FM method. We also show how to derive a radial velocity curve from the FM method, because many investigators have developed software to analyse binary star orbits using light curves combined with radial velocity curves. For example, the pre-eminent program for binary star analysis is now PHOEBE 10 , which is designed for radial velocity curve input (Prša & Zwitter 2005). We emphasise that the FM technique can provide directly all the information that is extracted from the radial velocity curve, but we also provide the method to generate the radial velocity curve itself, both for input to programs such as PHOEBE, and because the radial velocity curve is visually helpful for thinking about the binary star orbit.
We have discussed Doppler shift to help visualise what our techniques do, but the FM and PM methods do not rely directly on Doppler shift; they are built on the equivalent Rømer time delay. In Sections 2-4 of this paper we derive the relations to determine the mass function, orbital semi-major axis for a pulsating component, eccentricity, ̟ and tp, as well as the radial velocity curve, from the frequencies, amplitudes and phases of the components of a multiplet in the Fourier spectrum, which are split by the orbital frequency of the binary star.
In Section 5, we provide some examples using Kepler data for stars with high eccentricity. We compare the results with those obtained with the PM method, and discuss the circumstances for which FM or PM is the preferable technique. We show that the method is robust, and that it is capable of detecting companions with masses down to brown dwarf (m 0.08 M⊙) and even super-Jupiter masses (m 0.03 M⊙).

THE LIGHT TRAVEL TIME EFFECT ON PULSATION IN A BINARY STAR
Let us consider a star sinusoidally pulsating with a single angular frequency ω0 in a binary. We name the stars '1' and '2', and suppose that star 1 is pulsating. As the pulsating star moves in its binary orbit, the path length of the light between us and the star varies, leading to a periodic variation in the arrival time of the signal at Earth. The difference in the light arrival time, τ (t), compared to the case of a signal arriving from the barycentre of the binary system is given by where c is the speed of the light and v rad (t) denotes the radial velocity, due to the orbital motion, of the star 1 at the time t, where the epoch is the time at which the star passes the nodal point ('N' in Fig. 1) directed toward us. We follow the convention that the sign of the radial velocity is positive when the star is receding from us.
10 http://phoebe-project.org/1.0 Hereafter we call τ (t) the "time delay". With the help of equation (1) for the time delay, the observed luminosity variation at time t is written as where A is the amplitude of pulsational luminosity variation and φ denotes the pulsation phase at t = 0. In the following we derive the Fourier transform of the luminosity variation with the form of equation (2).

The radial velocity as a function of time
First, we have to derive the radial velocity as a function of time.
Let us suppose a plane that is tangential to the celestial sphere on which the barycentre of the binary is located, and let the z-axis that is perpendicular to this plane and passing through the barycentre of the binary be along the line-of-sight toward us (see Fig. 1). The orbital plane of the binary motion is assumed to be inclined to the celestial sphere by the angle i, which ranges from 0 to π. The orbital motion of the star is in the direction of increasing position angle, if 0 i < π/2, and in the direction of decreasing position angle, if π/2 < i π. We write the semi-major axis and the eccentricity of the orbit as a1 and e, respectively. Let ̟ be the angle measured from the nodal point N, where the motion of the star is directed toward us, to the periapsis in the direction of the orbital motion of the star. Also let f be the 'true anomaly', which is the instantaneous angle measured from the periapsis to the star in the direction of the orbital motion of the star, and let r be the distance between the barycentre and the star 11 . To make some complicated formulae derived later in this paper more easily traceable, we must repeat the fundamental derivation of the radial velocity as a function of time, although it was already given in Shibahashi & Kurtz (2012). The z-coordinate of the position of the star is written as and the radial velocity, v rad : With the help of the known laws of motion in an ellipse (e.g., Brouwer & Clemence 1961), and the radial velocity of star 1 along the line-of-sight is expressed as where Ω denotes the orbital angular frequency. The time dependence of radial velocity is implicitly expressed by the true anomaly 11 The range of i, the definition of ̟, and the sign of z are often differently given by different authors. It should also be noted here that we cannot distinguish a motion illustrated in Fig. 1 and a motion in the inverted image of the bottom panel of Fig. 1 from the radial velocity measurement alone. The former is a prograde motion, while the latter is a retrograde motion with the same value of ̟, but with the inclination angle being the supplementary angle of i shown in the top panel of Fig. 1. Schematic side view of the orbital plane, seen from a faraway point along the intersection of the orbital plane and the celestial sphere, NFN ′ , where the points N and N ′ are the nodal points, respectively, and the point F is the barycentre of the binary system; that is, a focus of the orbital ellipses. The orbital plane is inclined to the celestial sphere by the angle i, which ranges from 0 to π. In the case of 0 i < π/2, the orbital motion is in the direction of increasing position angle of the star, while in the case of π/2 < i π, the motion is the opposite. The z-axis is the line-of-sight toward us, and z = 0 is the plane tangential to the celestial sphere. Bottom: Schematic top view of the orbital plane along the normal to that plane. The periapsis of the elliptical orbit is P. The angle measured from the nodal point N, where the motion of the star is directed toward us, to the periapsis in the direction of the orbital motion of the star is denoted as ̟. The star is located, at this moment, at S on the orbital ellipse, for which the focus is F. The semi-major axis is a 1 and the eccentricity is e. Then OF is a 1 e. The distance between the focus, F, and the star, S, is r. The angle PFS is 'the true anomaly', f , measured from the periapsis to the star at the moment in the direction of the orbital motion of the star. 'The eccentric anomaly', u, also measured in the direction of the orbital motion of the star, is defined through the circumscribed circle that is concentric with the orbital ellipse.
f , which can be written in terms of 'the eccentric anomaly', u (see Fig. 1), defined through the circumscribed circle that is concentric with the orbital ellipse, satisfying and r sin f = a1 1 − e 2 sin u.
The eccentric anomaly u is written as where Jn(x) denotes the Bessel function of the first kind of integer order n, and tp denotes the time of the periapsis passage of the star. The trigonometric functions of the true anomaly f are expressed in terms of a series expansion of trigonometric functions of  . e(1 − e 2 ) −1/2 J ′ n (ne)/Jn(ne) for n = 1 (red), n = 2 (green), n = 3 (blue), and n = 4 (magenta) as functions of eccentricity e. This shows that e(1 − e 2 ) −1/2 J ′ n (ne)/Jn(ne) for n = 1, ..., 4 is ∼ 1 except for cases of large e, which means ϑ(e, ̟) ≃ ̟ except for cases of large e. This characteristic is used later in Sect. 3.3. the time after the star passed the periapsis with Bessel coefficients (see Appendix 1): where Jn ′ (x) := dJn(x)/dx. The radial velocity is then written explicitly as a function of time: where ξn(e, ̟) := 2 √ 1 − e 2 ne Jn(ne) and where arctan(x) returns the principal value of the inverse tangent of x. Since J ′ n (ne) Jn(ne) = 1 e − ne/2 (n + 1) − n 2 e 2 /2 (n + 2) − n 3 e 3 /2 (n + 3) − . . . , where the series in square brackets denotes continued fractions, the ̟-dependence of ξn(e, ̟) is weak for e ≪ 1, but it becomes substantial with the increase of e (see Fig. 3). Also, since the lowest order of the series expansion of ξn(e, ̟) in terms of power of e is e n−1 ; ξn(e, ̟) ∼ O(e n−1 ) for e ≪ 1. The band of each curve shows the range of ̟ from 0 to 2π. The ̟-dependence of ξ 2 , ξ 3 , and ξ 4 is substantially weaker than ξ 1 in a wide range of e. This figure is used later in Sect. 3.3.

The time delay as a function of time
Integrating equation (13), we obtain an expression for the time delay as a function of time where τ0(e, ̟) := − a1 sin i c ∞ n=1 ξn sin(−nΩtp + ϑn). Fig. 3 shows ξ1, ξ2, ξ3, and ξ4(e, ̟) as functions of eccentricity e. The band of each curve shows the range of ̟ from 0 to 2π. It is seen that the ̟-dependence of ξ2, ξ3, and ξ4 is substantially weaker in a wide range of e than is the case of ξ1. This is because the terms ξ2, ξ3, ξ4 and the higher terms represent a departure from a pure sinusoidal variation with the orbital angular frequency Ω, which is essentially caused by the eccentricity of the orbit. On the other hand, ξ1(e, ̟) shows substantial ̟-dependence, in particular with increase in e. This is because the projected light travel time is dependent not only on the eccentricity of the orbit, but also on the direction of the apparent elongation of the projected orbit, particularly in cases of high eccentricity.

Phase modulation
Shibahashi & Kurtz (2012) define a parameter α that measures the amplitude of the phase modulation when the pulsation frequency is treated as fixed: which is the ratio between the light travel time across the projected semi-major axis and the pulsation period of the mode in consideration. It should be noted that the value of α is dependent on the mode frequency. The larger the orbit size and the shorter the pulsation period, the larger the resulting phase modulation. Then the observed luminosity variation, whose amplitude is assumed to be unity, is rewritten as where and θn := ϑn − nΩtp.
If we regard equation (21) as luminosity variation with the intrinsic angular frequency ω0, the pulsation phase is regarded as being modulated by the orbital motion.

METHODOLOGY FOR DETERMINING BINARY PARAMETERS FROM PHOTOMETRY
The problem we address here is how to reproduce the radial velocity v rad (t) from the observed luminosity variation given by equa- or equivalently, how to deduce the orbital parameters, that is, the orbital angular frequency Ω, the semi-major axis a1, the eccentricity e, the angle between the periapsis and the nodal point ̟, and the time of periastron passage of the star tp.

Fourier transform of the luminosity variation
Here we describe the Fourier analysis of the luminosity variation, given by equation (21), showing periodic phase modulation. This is done with the help of Bessel functions, written with arbitrary real x and θ: Applying this, we rewrite equation (21) as where ℜ(x) means the real part of x and N denotes a large number. This means that a frequency multiplet appears around the intrinsic frequency ω0 in the frequency spectrum, and that each component of the multiplet is separated from its neighbouring components by the orbital frequency Ω: where A0 and A±m are the complex amplitudes of the central peak and the sidelobes separated from the central component by ±mΩ, respectively.

Truncation in the case of α < 1
It is straightforward to compute numerically the terms in equation (26) for a given set of binary parameters. The purpose of the present paper is, however, the reverse: we deduce the binary parameters from the Fourier transform of the light curve. An important point is that the complex amplitudes, A0 and A±m, are determined by the quantities that characterize the binary system, α, ξn(e, ̟), θn(e, ̟), allowing us to inversely determine those binary parameters from the amplitude spectrum of the oscillations.
In the following, we derive analytical expressions of the Fourier transform of the light curve in terms of the binary parameters. In practice, we truncate the infinite expansion series in equation (26) with a finite number of terms, N . There appear (2N +1) N cross-terms, and it is neither practical nor necessary to take account of all the terms. Rather, it is better to keep only the leading terms. Indeed, in the case of α < 1, the Bessel functions with the argument of αξn become negligibly small with the increase of n and the order kn. Hereafter, we limit ourselves to consider the cases of α < 1; that is, the cases for which the light travel time across the projected orbit is longer than the pulsation period.
Later, in Section 7, we will demonstrate some observational examples that show up to septuplets, nonuplets, and even undecuplets. So, here we truncate the infinite series of equation (26) with N = 5. Then, equation (26) produces (2N + 1) N = 11 5 = 161051 cross-terms. Among them, we discard terms separated from the central peak at ω0 further than ±6Ω and keep only the undecuplet. For each frequency component of mΩ, we keep only the terms up to O(α 2 ).
Then, equation (27) shows that the light-time effect in a binary star leads to a frequency undecuplet, ω = ω0 ±mΩ (m = 1, ..., 5), for which complex amplitudes of the first and second sidelobes are given as and respectively. The complex amplitudes of the higher order sidelobes are given in Appendix 2.

Amplitude ratios and phase differences between the sidelobes and the central component
Let A0 and φ0 be the amplitude and the phase of the central component of the multiplet at the angular frequency ω0, and let A±m and φ±m be the amplitudes and the phases of the ±m-th sidelobes, respectively. In the following, we derive the amplitude ratios and the phase differences between the sidelobes and the central component.
Also, the phase differences between the first sidelobes and the central component are given as Let us choose the zero point for the phases so that the phases of the first sidelobes are equal. This condition is fulfilled twice during one orbital period, at t = t0 being where k =    0 and 1 if 0 θ1 π/2 1 and 2 if π/2 θ1 3π/2 2 and 3 if 3π/2 θ1 < 2π.
At each of these phases, the Fourier component proportional to cos Ωt of the radial velocity becomes zero. In the case of a circular orbit, this corresponds to the moment at which the orbital motion of the star is perpendicular to the line-of-sight. At t = t0, which means that the phases of the first sidelobes apparently differ from that of the central peak approximately by either +π/2 or −π/2. This can be used to confirm that FM is caused by the orbital motion. Equation (35) also indicates that the deviation of the phase difference from +π/2 or −π/2 is dependent on the sign of sin(2ϑ2 − ϑ1). Since ϑi ≃ ̟, this means that the phase difference is slightly larger (smaller) if the periapsis is in the far (near) side of the orbit, with respect to us.

The case of the second sidelobes
Similarly, up to the order of O(αξ2), and Equation (37) indicates whether A+2 or A−2 has the higher amplitude is determined by the sign of cos(2ϑ1 − ϑ2). Note that if A+1 is higher (lower) than A−1, then A+2 is lower (higher) than A−2.
Also the phase differences are Combining equations (33) and (39), we obtain equivalently, It should be noted that there remains an uncertainty of 2nπ in the measured phase difference, φ+2(t0) − φ−2(t0), where n is an integer, and then we cannot uniquely determine (2ϑ1 − ϑ2) from this equation alone. This uncertainty can be eliminated, however, by taking account of the inequality of A±1 and the deviation of φ±1(t0) − φ0 from ±π/2.

General description: the higher-order sidelobes
Similarly, up to O(αξm), for the m-th sidelobes (m 1), Since the value of the left-hand side is observationally obtained, this equation determines the value of αξm. A practical numerical method of solving equation (42) is given in Section 4.1. The graphic solutions are shown in Fig. 4. Also, This relation allows us to deduce

Formal reconstruction of radial velocity from photometry
Equation (42) indicates that the value of αξm is observationally determined from the amplitude ratio between the sum of the sidelobes and the central peak. The graphical solution is shown in Fig. 4. The solution is numerically obtained with sufficient accuracy by expanding the Bessel functions Jn(x) up to the fifth power of x.
On the other hand, the phase differences between the sidelobes determine θm from equation (44), where θ1 is given by equation (33).
By substituting the values of αξm and θm into equation (24), For the whole range of 0 x 1, it is well approximated by a polynomial of the fifth order, which is derived from equation (17). From the value of (A +m + A −m )/A 0 , the value of αξm satisfying equation (42) is obtained.
we can formally reconstruct the radial velocity as a continuous function of time from photometry alone: where M is the number of components in the multiplet. In the limiting case of M → ∞, this is the time-varying radial velocity component, added to the systemic velocity, that would be measured by spectroscopy. The number of terms, M, is, however, very limited and not large enough in practical cases, hence we need to devise a more practical way of solving the problem. Nevertheless, equation (45) is instructive in the sense that it proves that the radial velocity can be reproduced from the frequency spectrum of the photometric variations.

FM METHOD
The amplitudes and phases of the multiplet components are determined by the binary orbit. Hence, inversely, once we have observed the frequency spectrum of a pulsating star in a binary system, we can extract information of the binary orbit. Essentially, this is the the FM method (Shibahashi & Kurtz 2012). While the strict values of the orbital elements can be obtained by analysing the continuous radial velocity, [equation (45)], if the number of components of the multiplet M is large enough, in practice M is limited to 3 or 4. With such a limited number of terms, equation (45) may still substantially deviate from the real radial velocity. In order to improve this, we derive the orbital parameters first and reconstruct the radial velocity in the case of M → ∞. It is then instructive to derive some 'easy solutions' of the orbital parameters to illustrate the technique. More approximate solutions can be obtained directly from the Fourier transform of the light curve without carrying out further numerical computations. In this section, we treat such cases.

The eccentricity, e
Equation (42) leads to the following relation between the amplitudes of the m + 1-th and m-th sidelobes: =: Rm(α, e, ̟).
The RHS of equation (46), defined as R(α, e, ̟), is a function of α, e and ̟, while the LHS is observable. Fig. 5 shows R(α, e, ̟) as a function of e, for m = 1, 2, and 3, for a fixed value of α = 0.217. The width of each curve shows the range of ̟ from 0 to 2π.
In the case of αξm ≪ 1, J1(αξm) ≃ αξm/2 and J0(αξm) ≃ 1, hence and then R is regarded as being no longer dependent on α. It should be noted that, while {ξm} are dependent on e and ̟, the ratios {ξm+1/ξm}, and then Rm, are dependent mainly on e, but relatively weakly on ̟ (see Fig. 5). A good initial guess for the eccentricity is then derived from the amplitude ratio of the (m + 1)-th sidelobes and m-th sidelobes. Indeed, as seen in equation (14), the ratio ξm+1/ξm is, in a good approximation, given by Whether or not αξm ≪ 1 can be checked by seeing the amplitude ratio of the sidelobe to the central component: =: Sm(αξm).
If αξm ≪ 1, Sm(αξm) ≃ αξm ≪ 1. So, if the amplitude ratio is substantially smaller than unity, the assumption of αξm ≪ 1 is justified. In such a case, the eccentricity is estimated as (Shibahashi & Kurtz 2012) Use of equation (49) for more than two values of m determines multiple values of e, each with its own uncertainty. The final value of e should be determined as a weighted average of those, which, given that Am decreases as m increases, will lead to the eccentricity being dominated by the ξ1 and ξ2 terms.

The angle between the periapsis and the nodal point, ̟
The angles ϑ1 and ϑ2 are regarded as functions of ̟ alone, since the eccentricity e has already been obtained. The angle ̟ between the periapsis and the nodal point is then derived from equation (41).
In the case of e ≪ 1, as seen in Fig. 2, the e-dependence of both ϑ1 and ϑ2 is weak, and 2ϑ1 − ϑ2 ≃ ̟. Hence, in such a case, The uncertainty in φ+2(t0) − φ−2(t0) by 2nπ can be eliminated by taking the inequality of A±1 and the deviation of φ±1(t0) − φ0 from π/2 into account.
It should be noted that the angle ̟ is obtained from information of the second sidelobes. If the amplitudes of A±2 are not statistically significant, we have to conclude that the orbit is close to a circular one.

The projected semi-major axis, a1 sin i
By substituting the values of e and ̟ thus obtained into the first line of equation (14), we estimate ξ1(e, ̟). The value of α is then obtained by dividing αξm (equation (51)) by ξm. Then, the projected semi-major axis is derived by equation (20):

4.4
The mass function, f (m1, m2, sin i) Once the value of a1 sin i is derived, with the help of the pulsation frequency and the orbital period, the mass function can be derived: where G is the gravitational constant. With a suitable asteroseismic estimate (or a reasonable assumption) of the primary mass, the minimum mass of the secondary star, m2 sin i, can be deduced.

The radial velocity curve deduced from photometry
Once e and ̟ are determined, we can deduce cos f and sin f from equations (11) and (12), respectively, as functions of the mean anomaly, Ω(t − tp), all the terms of ξn(e, ̟) and ϑ(e, ̟). Then substituting them, along with a sin i and Ω, into equation (7), we obtain the radial velocity as a function of the orbital phase.

Iteration for an improvement of parameters
Once the value of ̟ is estimated, we may iterate the above mentioned processes, by taking account of the ̟-dependences of ξm.
The value of ̟ newly derived in this way may be slightly different from that obtained previously. The iteration should be repeated until consistent solutions are obtained. An initial guess for ξ is obtained by ignoring ̟-dependence: ξ (0) m = ξ(e (0) , ̟ = 0). The first iteration gives ξ (1) = ξ(e (0) , ̟ (0) ). We improve the value of e as e (1) by adopting ξ (1) , and then ̟ (1) . We repeat these processes until we get consistent solutions. Fig. 3 is not used in this process, but it is illustrative.

Kepler mission data
The Kepler mission data available for this star are Q0-17 LC data, and Q5.1-5.3, Q7.1-7.3 short cadence (SC) data with 1-min integration time. An examination of the Q5 SC data shows that there are no significant peaks with amplitudes greater than 10 µmag at frequencies higher than 60 d −1 . In the frequency range 0 − 50 d −1 the four highest amplitude peaks were studied for FM. Fig. 6 shows the amplitude spectrum for the Q5 SC data in the frequency range 0 − 50 d −1 where the highest peaks chosen for analysis can be seen. They are in order of amplitude: ν1 = 19.47768 d −1 , ν2 = 30.80189 d −1 , ν3 = 21.71214 d −1 and ν4 = 36.14643 d −1 . There is no clear separation of p mode and g mode frequency ranges; there are peaks in all frequency ranges, suggesting p modes, g modes and mixed modes.
The Q0, Q1 and Q17 LC data were of shorter duration than full quarters, and they had a significant slope from the pipeline reductions. We therefore chose to analyse the Q2-16 LC data covering Table 1. A least-squares fit of the frequency septuplet for the highest amplitude mode to the Q2-Q16 Kepler data for KIC 9651065. The frequencies of the multiplet are separated by the orbital frequency, ν orb = 0.0036524 ± 0.0000036 d −1 (P orb = 273.8 ± 0.3 d). Column 4 shows that the first sidelobe phases are close to π/2 = 1.57 rad out of phase with the central peak, as required by theory for FM. The zero point for the phases has been chosen to be a time when the phases of the first sidelobes are equal, t 0 = BJD 2455783.05262, for the first multiplet. Column 5 shows that the phases of the first sidelobes of the other multiplets are equal within the errors at this time.  Fig. 7 shows the prewhitening process. The first, second and third FM sidelobes are significant and detected. There is some amplitude and/or frequency modulation of the ν1 mode, hence there is amplitude left near to that frequency after prewhitening. This is astrophysical, since amplitude and frequency variations on the time scale of this data set are common in δ Sct stars (see, e.g., Bowman & Kurtz 2014). This process was carried out for the other three of the highest four peaks. Figures are not shown for these other peaks. The results are given in Table 1. In the small sections of the amplitude spectrum between significant peaks the highest noise peaks have amplitudes of 8 µmag. The error on the amplitude determination is taken to be one quarter of that, with phase and frequency errors scaled accordingly (Montgomery & O' Donoghue 1999). The formal least-squares errors contain variance from all the peaks, so must be scaled down by this estimate.

FM analysis
For all four modes presented in Table 1 all sidelobes are equally split from the central frequency, being consistent with theoretical expectation of FM stars. This splitting gives the orbital frequency, ν orb = 0.0036524 ± 0.0000036 d −1 (P orb = 273.8 ± 0.3 d), where we have taken the best value here to be from the splitting of the first sidelobes of the highest amplitude mode. The zero point for the phases has been chosen to be a time when the phases of the first sidelobes of the highest amplitude peak are equal. The phases of the first sidelobes are close to π/2 rad out of phase with the central peak. That they are slightly less than this is a consequence of the values of the eccentricity and argument of periapsis. The same is true for the other modes within the error bars. Furthermore, the phase relation (φ+n − φ−n) for n = 2, 3, 4 is the same within the errors for the ν1 and ν2 modes. As for the ν3 and ν4 modes, the uncertainty is too large to confirm this. These facts are consistent with the theory for FM discussed in Section 3.3.
For each set of sidelobes, the values of mαξm and θm − mθ1 that were derived by equations (42) and (44), respectively, are summarized in Table 2. The eccentricity e is then estimated from the first and second sidelobes around ν1 by equation (52) to be e = 0.569 ± 0.030. The technique explained in Section 4.1 and inspection of Fig. 3 gives the coefficient ξ1(e, ̟) = 0.840±0.015. Using this value and the amplitude ratio for ν1, αξ1 = 0.2247 ± 0.0015, we obtain α = 0.2675 ± 0.0018. The mass function is then estimated to be 0.0992 ± 0.0057, which gives a minimum secondary mass of m2 = 0.87 ± 0.02 M⊙ based on an assumption of m1 = 1.7 M⊙, so the secondary is probably a main sequence G star. We derive the semi-major axis of the primary star about the barycentre from equation (20) to be a1 sin i = 0.378 ± 0.007 au from the data for ν1. The top, right panel shows the amplitude spectrum after prewhitening ν 1 where two equally spaced sidelobes can be seen easily. These are the first FM sidelobes. A careful examination shows the second and third FM sidelobes also. The bottom left panel shows the amplitude spectrum after prewhitening ν 1 and its first FM sidelobes, where the second and third sidelobes can be seen. The remaining variance close to the mode peak after pre-whitening the second and third sidelobes is seen in the bottom right panel; this is caused by amplitude modulation of the mode frequency over the time span of the data set, a common characteristic of δ Sct stars. The highest remaining peak may be that of an independent pulsation mode. Table 2. The estimated quantities of the terms for KIC 9651065, mαξm and θm, appearing in the expression of the radial velocity. Column 1 shows the frequency of the central peak, and Column 2 shows its ratio to the orbital frequency. Here c is the speed of light. Column 4 shows mαξm estimated from equation (42). Column 5 shows θm −mθ 1 estimated by equation (44). The facts that A+2 − A−2 < 0 and φ±1 − φ0 > −π/2 indicate 0 < 2ϑ1 − ϑ2 < π. The phase difference between φ+2 and φ−2, along with this constraint, leads to 2ϑ1 − ϑ2 = 2.17 ± 0.03 rad. With use of the derived value of e, this gives ̟ = 2.22 ± 0.04 rad. Fig. 8 shows the radial velocity curve of KIC 9651065 derived from the Kepler photometric data alone. This demonstrates again that the present FM method extracts binary information from the Kepler light curve without spectroscopic observations. Murphy et al. (2014) developed a different method, by which the phase modulation (PM) in the time domain of intrinsic pulsation frequencies of the star is tracked. In this method, first of all, the frequency of the central component of the multiplet is measured in the Fourier transform. Then the light curve is divided into short segments and, with the frequency fixed, the pulsation phase in each segment is measured by a least-squares method. This provides us with "time delays" τ (t) as a function of time, where τ (t) is defined as in equation (1). Then by carrying out a Fourier analysis, we obtain the amplitudes αω −1 0 ξn(e, ̟) and phases ϑn(e, ̟) given in equation (18) for n = 1, 2, 3, · · ·. The amplitude ratios between ξn and ξn+1 allow us to estimate the eccentricity e. On the other hand, the radial velocity is deduced as a function of time by taking the time derivative of the time delay, v rad (t) = c dτ dt .

Comparison with the result obtained by PM
Then, from the maximum and minimum radial velocity, along with the eccentricity deduced from the amplitude ratios, the projected semi-major axis is deduced;   Murphy et al. (2014) applied this PM method to several pulsating stars observed by the Kepler Mission, including KIC 9651065. They derived the orbital elements of this star, as well as the radial velocity curve. The results obtained in this paper by the FM method and their PM results are in good agreement (see Table 3).
It should be noted here that the PM method has been further developed by . The latest method derives directly the orbital elements without converting time delays to radial velocities, offering higher precision.

KIC 10990452
KIC 10990452 is a Kp = 12.4 binary system showing first, second, third and one component of the fourth FM sidelobes, indicating that the system is highly eccentric. According to Huber et al. (2014), its effective temperature and gravity are T eff = 7585 +244 −287 K and log g = 4.022 +0.126 −0.316 (cgs units), and its mass is m = 1.691 +0.332 −0.219 M⊙, which we round to T eff = 7600±250 K, log g = 4.02 ± 0.13 and m1 = 1.69 ± 0.33 M⊙. These are consistent with the original KIC parameters. We used Q1-16 data in our analysis, which have gaps because the star fell on the failed Module 3. Light curves in quarters Q7, Q11 and Q15 are therefore missing. The short 9.7-d engineering Q0 and the short 31-d final Q17 were not included. Fig. 9 shows the LC amplitude spectrum where it can be seen that there is a dominant highest peak. We examine that and the second highest amplitude peak next to it for FM. The top left panel of Fig. 10 shows the highest peak at ν1 = 17.72374 d −1 and the second highest peak close to it, ν2 = 17.85683 d −1 . Because of the complexity of the spectral window caused by the gaps in the data, the window patterns for the two frequencies interfere, so they must be modelled together. The top right panel shows the amplitude spectrum after the highest two mode peaks have been prewhitened; the FM sidelobes are apparent for both modes. The middle left panel is one step further with prewhitening of the first sidelobes, too. The remaining central amplitude of ν2 represents amplitude modulation over the time span of the data set -either astrophysical or instrumental. Given that ν1 does not show this as much, it is probably astrophysical. The middle right panel shows the third sidelobes and the low-frequency component of the fourth sidelobes. The bottom panel shows the amplitude spectrum of the residuals after the two nonuplets given in Table 4 have been prewhitened.
When searching for significant peaks associated with pulsation modes among thousands of independent frequencies in the amplitude spectrum -as we do in this case -it is found that the highest noise peaks have amplitudes about 4 times that of the rms noise in amplitude. The least-squares fitting that we do to find frequencies, amplitudes and phases estimates the errors based on the total variance in the data. Since there are so many other pulsation frequencies present that have not been modelled, the noise is not white and errors are significantly overestimated. We therefore looked at a section of the amplitude spectrum that is pure noise and found that the highest noise peaks are around 12 µmag. The amplitude error was therefore scaled to one quarter of that, or 3 µmag. The phase and frequency errors were scaled by the same factor, since the errors on those quantities are proportional to the amplitude signal-to-noise ratio (Montgomery & O' Donoghue 1999).
The results of the least-squares fits of the two pulsation mode peaks and their first, second, third, and fourth FM sidelobes are given in Table 4. All sidelobes are equally split from the central frequency, and is consistent with theoretical expectation for FM stars. This splitting gives the orbital frequency, ν orb = 0.008190 ± 0.000024 d −1 (P orb = 122.11 ± 0.36 d). The zero point for the phases has been chosen to be a time when the phases of the first sidelobes of the highest amplitude peak are equal. The phases of the first sidelobes are π/2 rad out of phase with the central peak. The same is true for the second highest amplitude mode with ν2. Furthermore, the phase relation (φ+m − φ−m) for m = 2, 3, 4 is the same within the errors for the two modes. These facts are consistent with theory for FM discussed in Section 3.3.
For each set of sidelobes, the values of mαξm and θm − mθ1 were derived by equations (42) and (44), respectively, and are sum- Figure 10. Top left: Amplitude spectrum for KIC 10990452 in the frequency range of the two highest amplitude p-mode peaks. The spectral window has significant sidelobes because of the missing quarters of data. Top right: This panel shows the first FM sidelobes after prewhitening by ν 1 and ν 2 ; further sidelobes can be seen. Bottom left: The this panel shows the data after prewhitening by the first sidelobes, too, where the second and third sidelobes are easily visible. Bottom right: In this panel the third sidelobes are easily seen for both main frequencies, and the low frequency fourth sidelobe is also seen for the ν 1 group. Table 4. A least-squares fit of the frequency nonuplets for the two highest amplitude modes to the Q1 to Q16 Kepler data for KIC 10990452. The frequencies of the multiplet are separated by the orbital frequency, ν orb = 0.008190 ± 0.000024 d −1 (P orb = 122.11 ± 0.36 d). The zero point for the phases has been chosen to be a time when the phases of the first sidelobes for the highest amplitude peak are equal, t 0 = BJD 2455673.01863. Column 4 shows that the first sidelobe phases are π/2 = 1.57 rad out of phase with the central peak, as required by theory for FM. Column 5 shows that the phases of the first sidelobes are equal within the errors at this time.   Table 5. The values of these quantities derived from ν1 and ν2 are in good agreement, and this is consistent with theoretical expectation of FM stars. The eccentricity e is then estimated from the first and second sidelobes around ν1 by equation (52) to be e = 0.569 ± 0.030.
The technique explained in Section 4.1 and inspection of Fig. 3 gives the coefficient ξ1(e, ̟) = 0.815 ± 0.015. Using this value and the amplitude ratio for ν1, αξ1 = 0.0640 ± 0.0015, we obtain α = 0.0785 ± 0.002. The mass function is determined from these values to be f (m1, m2, sin i) = 0.0163 ± 0.0011 M⊙. Assuming m1 = 1.7 M⊙ (Huber et al. 2014) gives a minimum secondary mass of m2 = 0.41 ± 0.01 M⊙, so the secondary is probably a main sequence M or K star. We also derive the semi-major axis of the primary star about the barycentre. From equation (20), we find a1 sin i = 0.122 ± 0.003 au from the data for ν1.

PM analysis
KIC 8264492 is a binary with a δ Sct star component. Huber et al. (2014) give T eff = 7992 +231 −315 K, log g = 3.947 +0.199 −0.224 (cgs units), and M = 1.87 M⊙. The binary nature of KIC 8264492 was discovered using the PM method (Murphy et al. 2014). The Fourier transform of its light curve (Fig. 12a) shows a dense frequency spectrum, in which the highest amplitude peak occurs at 31.29 d −1 . This peak is distinguished from its Nyquist alias at 17.7 d −1 by the larger amplitude of the former (Murphy, Shibahashi & Kurtz 2013). The non-sinusoidal nature of the time delays indicated the orbit is highly eccentric, which is confirmed by the high amplitude of the harmonics (see figure 21 of . From the amplitude ratios of those harmonics, a graphic solution for the eccentricity, e ≃ 0.73 ± 0.05, is obtained with the help of Fig. 5. Table 5. The estimated quantities of the terms for KIC 10990452, mαξm and θm, appearing in the expression of the radial velocity. They are estimated from the frequency nonuplets for the two highest amplitude modes, ν 1 and ν 2 , shown in Table 4. Column 1 shows the frequency of the central peak, and Column 2 shows its ratio to the orbital frequency. Here c is the speed of light. Column 4 shows mαξm estimated from equation (42). Column 5 shows θm estimated by equation (44).

FM analysis
The following FM analysis is based on the highest amplitude peak from Fig. 12a. We removed just over 5000 points from the start of the Q0-Q16 LC data set because of some trends and flux discontinuities in the early quarters. The data set we used has 64 031 LC points, and is 1318.9-d long (BJD 2 455 072.035 to 2 456 390.959). Fig. 13a zooms in on the highest peak at ν1 = 31.2920323 d −1 . The first sidelobes are clearly visible, separated from the central peak by 0.004 d −1 . Fig. 13b shows the amplitude spectrum after the central peak (at the dashed red line) has been prewhitened. The FM sidelobes are apparent, as shown by arrows, at exact multiples of the orbital frequency distant from the central peak. Table 7 gives the frequencies, amplitudes, phases and their uncertainties for the three highest amplitude modes and their detectable pairs of sidelobes. We determined the orbital frequency splitting from the splitting of the first pair of sidelobes with respect to the central component, ν orb = 0.0039621 ± 0.0000087 d −1 , giving P orb = 252.39 ± 0.56 d.
For each set of sidelobes, the values of mαξm and θm − mθ1 were derived by equations (42) and (44), respectively, and are summarized in Table 8. The eccentricity e is then estimated from the first and second sidelobes around ν1 by equation (52) to be e = 0.761 ± 0.045.
The technique explained in Section 4.1 and inspection of Fig. 3 gives the coefficient ξ1(e, ̟) = 0.630 ± 0.030. Using this value and the amplitude ratio for ν1, αξ1 = 0.369 ± 0.004, we obtain α = 0.577 ± 0.028. The mass function is determined from these values to be f (m1, m2, sin i) = 0.274 ± 0.040 M⊙. Assuming m1 = 1.87 M⊙ gives a minimum secondary mass of m2 = 1.44 ± 0.10 M⊙. We also derive the semi-major axis of the primary star about the barycentre. From equation (20), we find a1 sin i = 0.508 ± 0.024 au from the data for ν1.  31.29, 23.92, 33.76, 34.21, 23.54, 21.14, 24.47, 20.99 and 13.53 d −1 . The dashed red line represents the Nyquist frequency of the Kepler LC data. Table 7. A least-squares fit of the frequency undecuplet for the three highest amplitude modes in the KIC 8264492 data set. The zero-point in time was set at BJD = 2 455 754.1266, so that the first sidelobes of the highest amplitude mode had exactly equal phases. From the first pair of sidelobes, we obtain a frequency splitting of 0.0039621 ± 0.0000087 d −1 , giving P orb = 252.39 ± 0.56 d. Frequency uncertainties were calculated from a separate non-linear least-squares calculation where sidelobes were not forced to be equally split by the orbital frequency, and should therefore be taken as representative, only. Since only three modes were investigated, there is substantial variance left in the data. Treating the remaining variance would reduce the errors. The fact that the differences between the phases of the first sidelobes and the central peak for each mode are so close to π/2 confirms the binarity within the FM theory.   and θm, appearing in the expression of the radial velocity. They are estimated from the frequency nonuplets for the three highest amplitude modes, ν 1 , ν 2 and ν 3 , shown in Table 7. Column 1 shows the frequency of the central peak, and Column 2 shows its ratio to the orbital frequency. Here c is the speed of light. Column 4 shows mαξm estimated from equation (42). Column 5 shows θm estimated by equation (44). The minimum mass derived above puts the star in the mid-F range, or more massive. That suggests that some of the peaks seen in Fig. 12a may come from pulsations in the secondary. However, no evidence of anti-phase time delay variations in p modes are found, and g modes are low amplitude and not PM sensitive.
The upper panel of Fig. 14 shows the radial velocity curve of KIC 8264492 derived with the FM method from the Kepler photometric data alone. The time delay is also derived with the FM method. The lower panel of Fig. 14 shows the variation in the light arrival time with respect to the barycentre. This is in good agreement with time delays derived by the PM method (see ).

DISCUSSION AND CONCLUSIONS
We have shown in this paper that the Fourier transform of the light curve of a pulsating star in a binary system leads to frequency multiplets in the amplitude spectra and that all the binary orbital information traditionally derived from spectroscopic radial veloc- ity measurements can be derived from the frequency splitting and the amplitudes and phases of the components of the frequency multiplet. This is an extension of the FM method described in Shibahashi & Kurtz (2012), in which the case of e ≪ 1 was mainly discussed. We have improved the applicability to highly eccentric cases, and presented a method for determining binary orbital parameters, including the eccentricity and the argument of periapsis, from photometry alone. The orbital elements to be determined are (i) the orbital angular frequency, Ω, (ii) the projected semi-major axis, a sin i, (iii) the eccentricity, e, and (iv) the argument of the periapsis, ̟. Hence, we need at least four independent relations among the components of the frequency multiplet. The orbital angular frequency is directly measured from the frequency spacing of the adjacent components of the multiplet. The combination of the amplitude ratio of the sum of a pair of ±1-components to the central component, (A+1 + A−1)/A0, and the ratio for the case of ±2-components, (A+2 + A−2)/A0, gives a1 sin i and e. The phase information of the first and the second sidelobes leads to the argument of the periapsis. Careful examination of the phase differences also indicates whether the periapsis is located at the far side or the near side of the orbit, with respect to us. Hence, even if both components of the binary are pulsating, we can separate the pulsation spectra of two pulsating stars.
Binary information can be extracted by the PM method (Murphy et al. 2014;. The present FM method and the PM analysis are complementary. The PM analysis offers a clear visualization of the binary orbit in the time domain by dividing the light curve into segments, while the present FM method utilises the full frequency resolution of the data in the Fourier domain, and leads to high precision, particularly in the case of short-period binaries.
For pulsating stars, such as δ Sct stars, with many significant peaks in the amplitude spectrum, some care is needed to test for possible unresolved contamination of FM multiplet components by other, independent mode or combination frequencies. This is easily done for multi-periodic pulsators, because the amplitude and phase relationships for each FM multiplet must agree within the relations given in this paper. In the PM method, unresolved contamination results in obvious outliers in the time delay diagram, again making it easy recognise and discard contaminated mode frequencies.
A promising application of our method is in the search for exoplanets orbiting pulsating δ Scuti stars, or more generally, exoplanets around upper main-sequence stars. Both the transit method and the ground-based Doppler method, which are widely used in the search for exoplanets around solar-like stars and have discovered most of the known exoplanets to date, are more difficult to apply for upper main-sequence stars. Planetary transits are smaller in the case of upper main-sequence stars because of the much larger stellar disks. On the other hand, the typical rotationally broadened spectral lines of A-type and earlier type stars make it harder to detect the Doppler shift of spectral lines. Our new technique of pho-tometric measurement of binary orbital parameters opens a possibility of exoplanet detection around upper main sequence pulsating stars, such as δ Sct stars, and around compact pulsating stars.
Another application of the present method is at the opposite extreme, that is, in the search for invisible massive companions in binary systems. If a pulsating δ Sct star forms a binary system with a star born initially with mass less than about 8 M⊙ but more than the pulsating star itself, the massive component must have already evolved into a white dwarf while the less massive component is still in the main-sequence stage as a pulsating δ Sct star. The mass of the white dwarf must be less than the Chandrasekhar mass limit, but still be of ∼ 1 M⊙, and the star must be too faint to be detected. On the other hand, if a binary system is composed of a δ Sct star and a star born initially more massive than ∼ 8 M⊙ but less than ∼ 30 M⊙, the massive star must have already been turned to a neutron star, probably with mass ∼ 1 M⊙, and the star must be again undetectable except perhaps as a pulsar. Although the survival of the companion after the explosion of the massive star is uncertain from the theoretical viewpoint, we know of binary systems composed of a neutron star and either an early-type massive star, as a high-mass X-ray binary (HMXB), or a cool star, as a low-mass X-ray binary (LMXB). Hence it seems natural to expect a binary system of a neutron star and an A-type δ Sct star. The further extreme case is a star born with an initial mass 30 M⊙. Such a star becomes a black hole. A well-known binary system composed of a black hole candidate is Cygnus X-1, whose companion is a blue supergiant variable star. If a black hole forms a binary system with a δ Sct star, the present FM technique is a unique, promising method of finding such an exotic system.