# Revisiting the proposed circumbinary multi-planet system NSVS14256825

## Abstract

In this work we carry out an analysis of the observed times of primary and secondary eclipses of the post-common envelope binary NSVS14256825. Recently, Almeida et al. (2013) proposed that two circumbinary companions orbit this short-period eclipsing binary, in order to explain observed variations in the timing of mutual eclipses between the two binary components. Using a standard weighted least-squares minimisation technique, we have extensively explored the topology of parameter space of a single planet model. We find the data set to be insufficient to reliably constrain a one-companion model. Various models, each with similar statistical significance, result in substantially different orbital architectures for the additional companion. No evidence is seen for a second companion of planetary nature. We suspect insufficient coverage (baseline) of timing data causing the best-fit parameters to be unconstrained.

## 1 Introduction

The discovery of planets within binary star systems has recently sparked an increased interest in their formation, occurrence frequency and dynamical evolution (Portegies Zwart, 2013). Several techniques exist to detect additional bodies accompanying binary stars. In addition to the traditional radial velocity technique, Han (2008) outlines the possibility to infer such planets from microlensing observations. Recently, transiting circumbinary planets have been detected using the Kepler space telescope (Doyle et al., 2011; Welsh et al., 2012; Orosz et al., 2012a, b; Schwamb et al., 2013; Kostov et al., 2013). Furthermore, companions can be detected from pulsar timing measurements (Wolszczan & Frail, 1992). The formation and dynamical evolution of planets around binary star systems have been the subject of recent theoretical studies (Quintana & Lissauer, 2006; Haghighipour & Raymond, 2007; Marzari et al., 2009; Shi et al., 2012).

Utilising ground-based observations, a number of multi-planet systems around short-period eclipsing binary stars have been proposed in recent years (Lee et al., 2009; Beuermann et al., 2010; Potter et al., 2011; Qian et al., 2011). From measuring the times at minimum light (either primary and/or secondary eclipse) one can use the light-travel time (LTT) effect to detect additional companions by measuring periodic changes in the binary period (Irwin, 1952; Hinse et al., 2012a; Horner et al., 2012a). In contrast to other detection methods (radial velocity, microlensing and transit) the LTT technique is sensitive to massive companions on a long-period orbit: the semi-amplitude of the LTT signal scales with the companions mass and period as and , respectively. In addition, low-mass binary components will favour the detection of low-mass companions on short-period orbits (Pribulla et al., 2012).

From ground-based photometric observations, the first two-planet circumbinary system (HW Virginis, a.k.a HW Vir) was proposed by Lee et al. (2009). Additional multi-body systems of planetary nature were subsequently proposed by Beuermann et al. (2010); Marsh et al. (2013) (NN Serpentis, a.k.a NN Ser), Potter et al. (2011) (UZ Fornacis, a.k.a. UZ For) and Qian et al. (2011) (HU Aquarii, a.k.a. HU Aqr). Recently, Lee et al. (2012) proposed a quadruple system with two circumbinary sub-stellar companions orbiting the Algol-type binary SZ Hercules (a.k.a. SZ Her).

For a secure detection of a multi-planet circumbinary system, at least two criteria need to be satisfied. First, any period variation, due to additional companions, must be recurring and periodic in time. The data should extend over at least two complete cycles of the longest period. Second, the proposed system should be dynamically stable on time scales comparable to the age of the binary components. Horner et al. (2011) first studied the dynamical stability of the two planets in HU Aqr. Their study allowed them to conclude that the system is highly unstable with disruption times of a few hundred years. Subsequent studies of the same system were carried out by Hinse et al. (2012a), Wittenmyer et al. (2012) and Goździewski et al. (2012). The overall conclusion of these studies is that the planets, as proposed in the discovery work, are simply not feasible. More observational data is necessary before any further constraints can be imposed on the orbital parameters of any companions in the HU Aqr system.

The proposed planets orbiting the close binary system HW Vir (Lee et al., 2009) is another case where the proposed planets do not stand up to dynamical scrutiny (Horner et al., 2012b). In that case, the dynamical character of the HW Vir system was studied, and the planets proposed were found to follow highly unstable orbits most likely due to their crossing orbit architecture and relatively high masses. However, Beuermann et al. (2012b) presented new timing measurements of HW Vir allowing them to conclude stable orbits under the assumption of fixing some of the orbital elements in their least-square analysis.

The NN Ser system was also recently studied by Horner et al. (2012a). These authors found stable orbits for the Beuermann et al. (2010) solutions, if the planets are locked in a mean-motion resonant (MMR) configuration. However, an in-depth remodeling of the timing data renders the system unstable when all parameters are allowed to vary freely. Very recently, Beuermann et al. (2013) published additional timing data of NN Ser. Their re-analysis allowed them to conclude the existence of two companions orbiting the binary pair involved in a 2:1 MMR.

Unstable orbits in proposed multi-body circumbinary systems have not only been found among companions of planetary nature. The SZ Her system with two sub-stellar mass companions was recently investigated within a dynamical analysis (Hinse et al., 2012b). Here, the authors also found that the proposed companions followed highly unstable orbits.

In a recent work, Almeida et al. (2013) interpreted observed eclipse timing variations of the post-common envelope binary NSVS14256825 as being the result of a pair of light-travel time effect introduced by two unseen circumbinary companions. The proposed companions are of planetary nature, with orbital periods and years, and masses of and , respectively. Once again, however, a recent study (Wittenmyer et al., 2013) reveals that the proposed planetary system would be dynamically unstable on very short timescales - with most plausible orbital architectures being unstable on timescales of just a few hundred years, and only a small fraction of systems surviving on timescales of years.

The aim of this paper is as follows. In section 2 we present the available timing data of NSVS14256825, which forms the basis of our analysis. In particular we augment the timing measurements presented in Almeida et al. (2013) with three additional data points presented in Beuermann et al. (2012a). We also introduce the light-travel time model using Jacobian coordinates and outline the derivation of the minimum mass and projected semi-major axis for a single circumbinary companion along with a short description of our least-squares minimisation methodology. In section 3, we carry out a data analysis and perform a period analysis based on Fourier techniques and present our results describing the main properties of our best-fit solutions. In particular, we present results from finding a best-fit linear, best-fit quadratic and best-fit one-companion model. Finally, we summarise our results and discuss our conclusion in section 4.

## 2 Data acquisition and Jacobian light-travel time model

As the basis of this work we consider the same timing data set as published in Almeida et al. (2013). However, we noticed that three timing measurements published in Beuermann et al. (2012a) were not included in Almeida et al. (2013). We have therefore carried out two independent analysis based on the following data sets. Dataset I: Data as presented in Table 3 in Almeida et al. (2013). This data set spans the period from June 22, 2007 to August 13, 2012, corresponding to an observing baseline of around 5 years. Dataset II: Data as presented in Table 3 in Almeida et al. (2013) plus three data points (primary eclipse) from Beuermann et al. (2012a). The additional points are as follows. BJD days, BJD days and BJD days. The second data set spans the period from June 10, 1999 to August 13, 2012, corresponding to an observing baseline of around 13 years (i.e., doubling the time window).

The aim of considering the second data set (Dataset II) is to investigate the effect of the additional timing data on the overall best-fit solution and compare the results obtained from considering the first data set (Dataset I) since it covers a longer observing baseline.

The time stamps in Beuermann et al. (2012a) are stated using the terrestrial time (TT) standard while the times in Almeida et al. (2013) states timing measurements in the barycentric dynamical time (TDB) standard. However, the difference between these time standards (TT vs TDB) introduces timing differences on a milli-second (approx. 0.002s) level due to relativistic effects (Eastman et al., 2010). In light of the quoted measurement uncertainties (from the literature) of the eclipse timings in the two data sets, the two time stamps (TT and TDB) can be combined and no further transformation of one time standard to the other is necessary.

Considering the binary as an isolated two-body system and in the absence of mechanisms that cause period variations, the linear ephemeris of future (or past) eclipses is given by (Hilditch, 2001)

(1) |

where denotes the cycle number, is the reference epoch and is the nominal binary period. Additional effects that cause variations of the binary period would be observed as a systematic residual about this best-fit line.

We use the formulation of the light-travel time effect based on Jacobian coordinates (Goździewski et al., 2012). In the general case a circumbinary -body system is a hierarchical system and employing Jacobian elements therefore seems natural. This is particularly true for the case of a single companion (the first object in a hierarchical multi-body ensemble), where the Jacobian coordinate is equivalent to astrocentric coordinates and readily returns the geometric osculating orbital elements of the companion relative to the binary. Here we assume the binary to be a single massive object with mass equivalent to the sum of the two component masses. For a single circumbinary companion the LTT signal can be expressed as (Goździewski et al., 2012)

(2) |

where is the speed of light and is given as

(3) |

Here denotes the orbital eccentricity and measures the argument of pericentre of the companion relative to the combined binary representing the dynamical centre. The eccentric anomaly is given as . Following Goździewski et al. (2012) the semi-amplitude of the LTT signal is given as

(4) |

with measuring the speed of light, the semi-major axis, the inclination of the orbit relative to the skyplane. The quantities and denote the masses of the combined binary and companion, respectively.

In summary, the set represent the five free osculating orbital parameters for the companion with and denoting the orbital period and time of pericentre passage, respectively. These latter two quantities are introduced implicitly via Kepler’s equation and the eccentric anomaly (Goździewski et al., 2012; Hinse et al., 2012a)

### 2.1 Deriving minimum mass and projected semi-major axis

Once a weighted least-squares best-fit model has been found the minimum mass of the companion is obtained from solving the following transcendental function

(5) |

where

(6) |

The projected minimum (with ) semi-major axis () is then found from Kepler’s third law

(7) |

where the gravitational parameter is given by with denoting Gauss’ gravitational constant. The combined mass of the two binary components is assumed to be (Almeida et al., 2013).

Considering only the case of a single circumbinary companion, the timings of minimum light for primary eclipses is given as

(8) |

We therefore have a total of seven model parameters describing the light-travel time effect caused by a single circumbinary companion. For a description of two companions we refer to Goździewski et al. (2012). The LTT signal is a one-dimensional problem similar to the radial velocity
technique. We therefore only derive the minimum mass and minimum (projected) semi-major axis of the companion. For simplicity, we henceforth write and for the minimum masses and minimum semi-major axis of the companion^{1}

It is worth pointing out that no gravitational interactions have been taken into account in the above formulation of the LTT signal. Only Keplerian motion is considered. It is possible to include additional effects (such as mutual gravitational interactions) that can cause period variations and we refer to Goździewski et al. (2012) for more details.

Finally, we stress that the case of a single companion the Jacobian-based description of the one-companion LTT effect is equivalent to the formulation given in Irwin (1952, 1959). Hence, the should be identical to those parameters obtained using the Irwin (1952) LTT model along with the derived minimum mass. The only parameter which is different is the semi-major axis of the binary due to the different reference systems used and we refer the reader to Goździewski et al. (2012) for details. For consistency, we tested our results for the presently (Jacobian) derived LTT formulation using the procedure detailed in Irwin (1952), and obtained identical results. However, one complication could arise in the argument of pericentre which can differ depending on the defined direction of the line-of-sight axis. Either this axis can point towards or away from the observer. The difference will affect the argument of pericentre and can be rectified using the relation , where is the argument of pericentre defined in a reference system with opposite line-of-sight direction compared to the formulation outlined in Goździewski et al. (2012). Hence the difference is only a matter of convention and does not affect the quantitative results obtained from the two formulations.

### 2.2 Weighted least-squares fitting

We have implemented the Jacobian-based Kepler-kinematic LTT model in IDL^{2}

(9) |

where is the number of data points, measures the vertical difference between the observed data and the computed model at the th cycle, and measures the 1-sigma timing uncertainty (usually obtained formally). However, in this work we will quote the reduced defined as with denoting the degree of freedom. The MPFIT routine attempts to minimise iteratively using free parameters.

In the search for a global minimum of the underlying space we utilise a Monte Carlo approach by generating a large number of random initial guesses. Two approaches can be used to explore the space for a global minimum. The first involves generating random initial guesses in a relatively narrow region of a given parameter and may be applied when information about the periodicity and amplitude of the LTT signal is inferred from other means (e.g. Fourier analysis). For example, if a Fourier analysis reveals a given frequency within the data one can then generate random initial guesses from a normal distribution centred at that period with some (more or less narrow) standard deviation for the variance. In the second approach, random initial guesses are generated from a uniform distribution defined over a broad interval for a given parameter. However, in both approaches we randomly choose the eccentricity from a uniform distribution within with the argument of pericentre chosen from . In all our searches we recorded the initial guess and final parameters along with the goodness-of-fit value, the corresponding root-mean-square (RMS) statistic and formal 1-sigma uncertainties. A single LM iteration sequence is terminated following default values of accuracy parameters within MPFIT or after a maximum of 3000 iterations (rarely encountered with the average number of iterations required being just 11).

## 3 Data analysis & results

### 3.1 Period analysis and linear ephemeris

As a starting point for our analysis, we first determined the parameters of the linear ephemeris () by calculating a linear least-squares regression line to the same data (Dataset I) as considered by Almeida et al. (2013). A best fit line resulted in a with free parameters and data points. The corresponding value was found to be 853 and the (rounded) linear ephemeris was determined to be

(10) | |||||

(11) | |||||

(12) |

For Dataset II we obtained the slightly different ephemeris, with little improvement in the precision of the binary period

(13) | |||||

(14) | |||||

(15) |

We applied the PERIOD04^{3}

### 3.2 Quadratic ephemeris model - Dataset I

In some cases a change of the binary period can be caused by non-gravitational interaction between the two components of a short-period eclipsing binary. Often the period change is described by a quadratic ephemeris (linear plus secular) with the times of primary eclipses given by Hilditch (2001)

(16) |

where is a period damping factor (Goździewski et al., 2012) which can account for mass-transfer, magnetic braking, gravitational radiation and/or the influence of a distant companion on a long-period orbit. Following Brinkworth et al. (2006) the rate of period change, in the case of mass-transfer, is then given by

(17) |

with denoting the currently measured binary period. In Fig. 3 we show the best-fit quadratic ephemeris given as

(18) | |||||

(19) |

with unreduced for (67-3) degrees of freedom. In Fig. 4 we show the location of the best-fit surrounded by the (68.3%, ), (95.4%, ) and (99.7%, ) joint-confidence contours (Press et al., 2002; Bevington & Robinson, 2003; Hughes & Hase, 2010) for the parameter space. Similar results were obtained for the remaining two parameter combinations. Considering Dataset I we found the average period change to be . This value is about one order of magnitude smaller than the period decrease reported in Almeida et al. (2013).

### 3.3 Single companion model - Dataset I

To reliably assess the validity of a two-companion model we first considered a one-companion model. Our period analysis yielded a shortest period of around years (2557 days) with a semi-amplitude of days (20 seconds). We therefore searched for a best-fit solution in a narrow interval around these values by seeding 523,110 initial guesses. The best-fit solution with is shown in Fig. 5. In Table 1 we show the corresponding best-fit parameters and derived quantities for the companion along with their formal (derived from the covariance matrix) uncertainties as obtained from MPFIT. Formal errors in the derived quantities were obtained from numerical error propagation, as described in Bevington & Robinson (2003). The residual plot in Fig. 5 (middle panel) shows no obvious trend above the 5 seconds level. The average timing uncertainty in the Almeida et al. (2013) data set is 5.5 seconds. An additional signal associated with a light-travel time effect should be detected on a level equivalent to a timing semi-amplitude of seconds. Usually timing measurement are assumed to distribute normally around the expected model. We have therefore also plotted the normalised residuals (Hughes & Hase, 2010) as shown in the bottom panel of Fig. 5. The corresponding histogram is shown in Fig. 6. Whether the timing residuals follow a Gaussian distribution is unclear at the moment.

Again, we have explored the function in the vicinity of the best-fit parameters and determined two-dimensional joint-confidence intervals. We show all 21 two-parameter combinations in Fig. 7 and Fig. 8. While the two considered parameters in a given panel were kept fixed, we allowed all the remaining parameters to re-optimise (with an initial guess given by the best-fit values listed in Table 1) during a LM iteration (Bevington & Robinson, 2003).

We note that several of the parameters correlate with each other. This is especially true for the pair shown in Fig. 8. Choosing our reference epoch to be close to the middle of the data set results in almost no correlation between and (see top left panel in Fig. 7). In addition, we note that the topology around the best-fit parameters deviates from its expected parabolic form. This is most readily apparent in the panel in Fig. 8.

Finally, we note that the confidence level in the (bottom-right) panel of Fig. 7 appears open, and stretches toward longer periods and larger semi-amplitudes (). With this in mind, we then recalculated the space of considering a larger interval in the two parameters. The result is shown in Fig. 9, demonstrating that the 3-sigma joint-confidence contour remains open for orbital periods larger than around 22 years. We therefore suspect that our best-fit model resides within a local minimum.

To test whether we are dealing with a local minimum we explore the parameter space on a wider search grid by following the approach as outlined previously. Surprisingly, we found a marginally improved solution with a smaller best-fit value of 1.96, a reduction by compared to the first best-fit solution of 1.98. Computing the space around the new best-fit solution over a large interval in the parameters and resulted in Fig. 10.

In each panel, our (new) improved best-fit solution is marked by a cross-hair. The corresponding model parameters are shown in Table 2. We have omitted quoting the formal uncertainties for reasons that will become apparent shortly. In Fig. 10 we also show the 1-sigma (68.3%) joint-confidence contour of (black line) encompassing our best-fit model. Our results suggest that a plethora of models, with within the 1-sigma confidence level, are equally capable of explaining the timing data. Statistically, within the 1-sigma uncertainty region, essentially no differences in the exist between the various solutions.

For this reason, the considered parameters (semi-amplitude, eccentricity and period) span a vast range, making it impossible to place firm confidence intervals on them. From Fig. 10 possible periods span from days (6.8 years) to at least 80,000 days (219 years) chosen as our upper cut-off limit in the search procedure. We have tested this result by selecting three significantly different pairs of in Fig. 10a. We label them as follows: Example 1): . Example 2): . Example 3): . We re-calculated a best-fit model with the parameters held fixed, and remaining parameters allowed to vary freely (starting from the best-fit solution given by the cross-hair in Fig. 10a) to find new optimum values. We show the results of this experiment in Fig. 11. All models have within the 1-sigma confidence level (1.993) but differ significantly in their orbital periods, eccentricities and semi-amplitudes. Our best-fit model (cross-hair) is shown in Fig. 11d and Table 2. We calculated the following values for the companion’s minimum mass and semi-major axis for our three examples. Example 1): , au. Example 2): , au. Example 3): , au. In light of the large range of possible parameters we omit quoting parameter uncertainties. Minimum mass and semi-major axis for our improved best-fit solution (Fig. 11d) are given in Table 2.

Up to this point our analysis allows us to conclude that the data is not spanning a sufficiently long observing baseline to firmly constrain the parameters of a single companion model. We stress that the model itself could still be valid. With the data currently at hand it is impossible to establish firm confidence intervals on the parameters. Our first solution (comparable with the solution presented in Almeida et al. (2013)) likely represents a local minimum in the parameter space, or appears to be a solution within the confidence interval characterised by a shallow topology of space. In such a case we cannot distinguish isolated models in the continuum of possible solutions. All three panels in Fig. 10 indicate the existence of local minima with statistics close to our first best-fit solution with (Table 1). In fact, Fig. 10 suggests the existence of multiple local minima in the space. Since in Fig. 10 we have not found the 1-sigma confidence level to render as a closed-loop contour line, we suspect that the data can be fit to an infinite number of models each having the same statistical significance, but exhibiting significant differences in their orbital architectures. In light of this result any efforts to search for a second companion in Dataset I seems unfruitful.

### 3.4 Single companion model - Dataset II

We have noted that three datapoints from Beuermann et al. (2012a) were not included in the analysis presented in Almeida et al. (2013). Although they are accurate (placing them well on the linear ephemeris) their timing precision is lower. However, the large timing uncertainty for these points should not disqualify them from being included in the analysis. In principle, the precision of the eclipsing period should increase for a dataset of increased baseline, and could eventually help to constrain any long-period trend. We have repeated our search procedure as outlined previously to find a best-fit model based on dataset II. We show our best-fit solution in Fig. 12 and state the best-fit parameters within the figure area.

For dataset II, the main characteristics of the Keplerian orbit for the companion are similar to the parameters shown in Table 2. The period, minimum semi-major axis and eccentricity are comparable in both cases. We also explored the topology of space for a large region around the best-fit solution and found similar results as discussed previously by generating two-dimensional joint-confidence interval maps. The 1-sigma confidence contour around the best-fit solution extends over a large interval in the period, eccentricity and semi-amplitude.

From examining the residual plot in Fig. 12 we are not convienced about any additional light-travel time periodicity above the RMS level of about six seconds. A light-travel time signal with amplitude of around six seconds would require a dataset with RMS of about one second or less. Hence, from a qualitative judgment, the data in Dataset II does not currently support the inclusion of five additional parameters corresponding to a second companion. The results from examining Dataset II reinforces insufficient coverage of the orbit as presented in (Almeida et al., 2013). Because Dataset II covers two-times the best-fit period found for Dataset I, one would expect Dataset II to constrain the orbital period to a higher degree than for Dataset I. However, this is not the case for the present situation.

## 4 Summary and conclusions

In this work we have carried out a detailed data analysis of timing measurements of the short-period eclipsing binary NSVS14256825. In particular we have examined the one-companion model bearing in mind that additional valid companions should be readily visible in the resulting residuals. On the basis of Dataset I, we first carried out an initial local search for a weighted least-squares best-fit solution. A best-fit model (Table 1) with resulted in an inner circumbinary companion with orbital characteristics comparable to the short-period companion presented in Almeida et al. (2013). Extending our search grid of parameter space resulted in a similar best-fit statistic with significantly different orbital characteristics (Table 2). We were able to show quantitatively that the present timing data does not allow us to firmly constrain a particular model with well-established parameter confidence limits. In light of this, quoting formal errors for the model parameters seems meaningless. We concluded that the best-fit solution found by Almeida et al. (2013) most likely represents a local minimum. We explain the lack of constraint in the parameters by the limited monitoring baseline over which timing data was acquired. Dataset I represented a baseline of about 5 years. If a periodicity is present, the principle of recurrance should apply, requiring two full orbital periods to be covered in order to establish firm evidence for the presence of a companion. This would correspond to a light-travel time period of at most 2.5 years for Dataset I and 6 years for Dataset II (spanning about 13 years). However, for Dataset I, the data did not allow models with periods shorter than days. Simultaneously Dataset II does not constrain the period any better than Dataset I.

Our analysis did not allow us to find convincing evidence of a second light-travel time signal. The RMS scatter of timing data around the best-fit model was found to be around 5 seconds. Signals with a semi-amplitude comparable with the measurement uncertainties seem unlikely to be supported by the present data. The claimed second companion in Almeida et al. (2013) has a semi-amplitude of seconds. It is likely that noise was wrongly interpreted as a light-travel time signal. We recommend that a secure detection requires a signal semi-amplitude of at least three times above the noise level (i.e ). Future timing data (Pribulla et al., 2012; Park et al., 2012) of this system will be important to help constraining the parameters significantly.

#### Acknowledgements

Research by T. C. H is carried out at the Korea Astronomy and Space Science Institute (KASI) under the 2012 KRCF (Korea Research Council for Science and Technology) Young Scientist Research Fellowship Program. K. G. is supported by Polish NSC, grant N/N203/402739. Numerical computations were partly carried out using the SFI/HEA Irish Centre for High-End Computing (ICHEC) and the PLUTO computing cluster at KASI. Astronomical research at Armagh Observatory is funded by the Department of Culture, Arts and Leisure (DCAL). T. C. H and J. W. L acknowledges support from KASI grant 2013-9-400-00. J. H. gratefully acknowledges financial support of the Australian government through ARC Grant DP0774000. R.A.W is supported by a UNSW Vice-Chancellor’s Fellowship.

Dataset I | ||

, , , | ||

RMS | 5.4 | seconds |

2,455,408.74450(36) | BJD | |

0.11037415(5) | days | |

AU | ||

- | ||

radians | ||

2,455,197(67) | BJD | |

days | ||

AU | ||

- | ||

radians | ||

days |

Dataset I | ||

, , , | ||

RMS | 5.3 | seconds |

2,455,408.74455(41) | BJD | |

0.11037411(6) | days | |

AU | ||

- | ||

radians | ||

2,455,330 | BJD | |

days | ||

AU | ||

- | ||

radians | ||

days |

### Footnotes

- Technically, the values obtained represent the minimum possible values of and - but in standard papers dealing with eclipse timing or radial velocity studies authors use the shortened versions, for brevity.
- http://www.exelisvis.com/ProductsServices/IDL.aspx
- http://www.univie.ac.at/tops/Period04/

### References

- Almeida, L. A., Jablonski, F., Rodrigues, C. V., 2013, ApJ, 766, 11
- Beuermann, K., Hessman, F. V., Dreizler, S. et al. 2010, A&A, 521, L60
- Beuermann, K. et al., 2012, A&A, 540, 8
- Beuermann, K., Dreizler, S., Hessman, F. V., Deller, J., 2012, A&A, 543, 138
- Beuermann, K., Dreizler, S., Hessman, F. V., 2013, arXiv:1305.6494
- Bevington, P. R., Robinson, D. K., ”Data reduction and error analysis for the physical sciences”, 2003, 3rd edition, McGraw-Hill, USA
- Brinkworth, C. S., Marsh, T. R., Dhillon, V. S., Knigge, C., 2006, MNRAS, 365, 287
- Doyle, L. R., et al., 2011, Science, 333, 1602
- Eastman, J., Siverd, R., Gaudi, B. S., 2010, PASP, 122, 935
- Goździewski, K., Nasiroglu, I., Słowikowska, A., Beuermann, K., Kanbach, G., Gauza, B. et al., 2012, MNRAS, 425, 930
- Han, C., 2008, ApJ, 676, 53
- Haghighipour, N., Raymond, S. N., 2007, ApJ, 666, 436
- Hilditch, R. W., ”An introduction to close binary stars”, 2001, Cambridge University Press, UK
- Hinse, T. C., Lee, J. W., Goździewski, K., Haghighipour, N., Lee, C.-U., Scullion, E. M., 2012, MNRAS, 420, 3609
- Hinse, T. C., Goździewski, K., Lee, J. W., Haghighipour, N., Lee, C.-U., 2012, AJ, 144, 34
- Horner J., Marschall J. P., Wittenmyer R., Tinney C. G., 2011, MNRAS, 416, L11
- Horner, J., Wittenmyer, R. A., Hinse, T. C., Tinney, C. G., 2012, MNRAS, 425, 749
- Horner, J., Hinse, T. C., Wittenmyer, R. A., Marshall, J. P., Tinney, C. G., 2012, MNRAS, 427, 2812
- Hughes, I. G., Hase, T. P. A., ”Measurements and their Uncertainties - A Practical Guide to Modern Error Analysis”, 2010, Oxford University Press, UK
- Irwin, J. B., 1952, ApJ, 116, 211
- Irwin, J. B., 1959, AJ, 64, 149
- Kostov, V. B., McCullough, P., Hinse, T. C. et al. 2013, ApJ, accepted and in press
- Lee J. W., Kim S.-L., Kim C.-H., Koch R. H., Lee C.-U., Kim H. I., Park J.-H., 2009, AJ, 137, 3181
- Lee, J. W., Lee, C.-U., Kim, S.-L., Kim, H.-I., Park, J.-H. 2012, AJ, 143, 34
- Lenz, P.; Breger, M.; 2005, CoAst, 146, 53
- Markwardt, C. B., 2009, ASPCS, ”Non-linear Least-squares Fitting in IDL with MPFIT”, Astronomical Data Analysis Software and Systems XVIII, eds. Bohlender, D. A., Durand, D., Dowler, P.
- Marsh, T. R., Parsons, S. G., Bours, M. C. P. et al. 2013, pre-print, arXiv:1310.1391
- Marzari, F., Scholl, H., Thébault, P., Baruteau, C., 2009, A&A, 508, 1493
- Orosz, J. A., et al., 2012, ApJ, 758, 87
- Orosz, J. A., et al., 2012, Science, 337, 1511
- Title:”Korea Microlensing Telescope Network: science cases” Park, B.-G., Kim, S.-L., Lee, J. W., Lee, B.-C., Lee, C.-U., Han, C., et al., 2012, SPIE, 8444, 47, (2012SPIE.8444E..47P)
- Potter, S. B., Romero-Colmenero, E., Ramsay, G. et al. 2011, MNRAS, 416, 2202
- Press, W. H. et al., ”Numerical Recipies in Fortran90: The Art of Scientific Computing”, 2002, Cambridge University Press, UK
- Pribulla, T. et al., 2012, AN, 333, 754
- Qian, S.-B., Liu, L., Liao, W.-P. et al. 2011, MNRAS, 414, L16
- Quintana, E. V., Lissauer, J. J., 2006, Icarus, 185, 1
- Schwamb, M. E., Orosz, J. A., Carter, J. A., Welsh, W. F. et al., 2013, ApJ, 768, 127
- Shi, J.-M., Krolik, J. H., Lubow, S. H., Hawley, J. F., 2012, ApJ, 749, 118
- Welsh, W. F., et al., 2012, Nature, 481, 475
- Wittenmyer, R. A., Horner, J., Marshall, J. P., Butters, O. W., Tinney, C. G., 2012, MNRAS, 419, 3258
- Wittenmyer, R. A., Horner, J., Marshall, J. P., 2013, MNRAS, 431, 2150
- Wolszczan, A., Frail, D. A., 1992, Nature, 355, 145
- Portegies Zwart, S., 2013, MNRAS, 429, 45