∎

11institutetext: P. Mendes 22institutetext: Center for Quantitative Medicine
and Department of Cell Biology, University of Connecticut School of Medicine,
Farmington, CT 06030, USA
Tel.: +1-860-6797632
22email<EMAIL_ADDRESS>
# Reproducible research using biomodels

Pedro Mendes

(Received: date / Accepted: date)

###### Abstract

Like other types of computational research, modeling and simulation of
biological processes (biomodels) is still largely communicated without
sufficient detail to allow independent reproduction of results. But
reproducibility in this area of research could easily be achieved by making
use of existing resources, such as supplying models in standard formats and
depositing code, models, and results in public repositories.

###### Keywords:

Reproducibility Biomodels

††journal: Bulletin of Mathematical Biology

##  1 Introduction

Reproducibility is at the core of the scientific process. A basic aspect of
science is the ability to establish reproducible effects. Karl Popper’s notion
of falsifiability depends on reproducibility: a theory is falsified when
reproducible effects refute it popper1959 . This means that results of
scientific investigations must be reproducible or else support for
“discoveries” becomes discredited (e.g. see maddox88 ). Therefore scientific
reports must describe experiments in sufficient detail to allow other
researchers to reproduce them. A recent survey suggests that there is a
problem with lack of reproducibility in a large proportion of scientific
articles baker2016 .

This reproducibility crisis has been highlighted mainly for experimental
research. At first glance it could appear that computational research would
not suffer from such problems since, after all, computers follow specific sets
of instructions (programs) and thus can be run in a reproducible manner. This
idea was formalized in the 1990s by Claerbout and Karrenbach claerbout1992
who described how electronic publications could easily be made into
reproducible publications. For example they could allow the reader to re-
execute analyses and plots directly from the data through an appropriately
associated program. These authors already highlight, though, that a critical
issue underlying computational reproducibility is that the programs used
should be open source in order to be available to all readers. Claerbout and
Karrenbach were essentially optimistic in how computational resources were
going to make publications more reproducible, writing that “With workstations
becoming widespread and software available, the burdens imposed on the author
to create reproducible results are little more than the task of filing
everything systematically” claerbout1992 . Unfortunately this optimism did not
materialize.

It is now widely accepted that the problem with reproducibility extends also
to computational research mesirov2010 ; peng2011 ; stodden2016 . And while
some journals have adopted a few measures to minimize this problem
greenbaum2017 ; guerreiro2017 ; loew2015 ; peng2009 , a large proportion of
articles describing computational research are still hard or impossible to
reproduce hothorn2009 ; hothorn2011 ; hubner2011 ; stodden2018 .

##  2 Which reproducibility?

It is rather unhelpful that the word ‘reproducibility’ has been used with
different meanings in this context, as summarized by Plesser plesser2018 . The
terminology proposed by Goodman et al. goodman2016  seems the most appropriate
and I will follow it here. These authors distinguish three different types of
reproducibility:

  * •

reproducibility of methods requires one to be able to exactly reproduce the
results using the same methods on the same data;

  * •

reproducibility of results requires one to obtain similar results in an
independent study applying similar procedures;

  * •

reproducibility of inferences requires the same conclusions to be reached in
an independent replication potentially following a different methodology.

All of these types of reproducibility are desirable but should be addressed
differently. The reproducibility of inferences requires that similar
conclusions be reached when an independent approach is applied to a problem.
Thus this aspect is perhaps the least problematic since it can be satisfied by
clearly describing the problem and the conclusions. On the other end of the
spectrum, to achieve reproducibility of methods in computational research
requires that the exact same software and input data be made available. To
achieve reproducibility of results the methods and algorithms must be well
defined and the input data is still required. Claerbout and Karrenbach were
right in that to achieve all of these levels of reproducibility in
computational research “merely” requires to describe all steps in detail — but
therein lies the problem.

##  3 Biomodels

A subset of computational research concerns the use of dynamic models of
biological systems (‘biomodels’). This includes a wide range of models, from
those representing basic physico-chemical phenomena, such as ligand-receptor
binding, all the way to models of disease transmission in populations or of
entire ecosystems. Models of biochemical reaction networks or pathways are
perhaps the most widely represented in this class. Computational biomodels
have been reported since the dawn of computers (e.g. chance1960 ) and are now
increasingly used to help understand phenomena and make predictions, as can be
witnessed by reading the pages of this journal and many others.

Research with biomodels may be better equipped to deal with reproducibility
than other types of computational research. A number of researchers have
agreed on various standards and principles that promote reproducibility. Most
notorious of all is the systems biology markup language (SBML, hucka2003 ) a
specification for biomodels that many software packages can read and write,
thus promoting reproducibility (see also CellML lloyd2004 ), NeuroML
gleeson2010 , and PharmML swat2015 ). Standards have also been proposed for
model diagrams (SBGN, le_novere2009 ), simulation specifications (SED-ML,
waltemath2011 ), and data (SBRML, dada2010 ). Finally there are also minimal
information recommendations for publishing models (MIRIAM, miriam2005 ) and
simulations (MIASE, miase2011 ) that, if followed, assure a good level of
reproducibility — essentially formalizing the process of “filing everything
systematically” advocated by Claerbout and Karrenbach claerbout1992 .

SBML allows biomodels to be specified in a manner that describes the biology
and mathematics without prescribing the algorithms to be applied. An SBML file
describes the transformations (“reactions”) that variables (“species”) can
undergo; the kinetics of the transformations are well specified with all
necessary constants, and the initial state of the system is also included. In
essence an SBML file contains all that is needed for software to construct and
solve the equations. Because the algorithms are not prescribed, this allows
SBML models to be used in different contexts and analyzed with different
formalisms in addition to those used in the original research.

With all the standards mentioned above computational results obtained with
biomodels can easily be made reproducible. Publishing the biomodels as SBML
files, either as attached supplementary material or by inclusion in a database
like BioModels chelliah2015 , satisfies reproducibility of results because it
makes the model immediately available for simulation with a range of software
(e.g. COPASI hoops2006 , VCell moraru2008 , and many others). Then, by using
the same software as the authors used originally, even reproducibility of
methods is achieved. Because SBML does not prescribe the mathematical
formalism, in some circumstances it can also facilitate reproducibility of
inferences. For example, conclusions may have been obtained from an analysis
using the linear noise approximation (e.g. pahle2012 ), and others may
reproduce the same conclusions using the Gillespie stochastic simulation
algorithm gillespie1976 . This is possible because both methods are
implemented in software packages that read SBML.

##  4 Publication of electronic materials

The most basic aspect to promote reproducibility of research using biomodels
is to publish the electronic materials used. This includes any programs used,
the input data (usually constants and initial conditions), and results. Where
should these materials be published? There are a number of options in current
use: public repositories; journal website as supplementary materials; author’s
website; supplied by author upon reader’s request — listed in decreasing order
of utility to the research community.

The choice of supplying materials only when requested is not practical and
there is plenty of evidence that too frequently it is not honored by authors
stodden2018 ; additionally those materials become inevitably lost as authors
retire or die. Any results relying on materials “supplied upon request” should
be considered non-reproducible and journal editors should simply not allow
this practice.

Publication of materials in authors’ web sites is only minimally better: for a
short while the materials are indeed immediately available to all, but dead
links appear at a fast pace. This is due to frequent website redesigns,
authors moving to other institutions, retirement, etc.

Publication of materials in public repositories and in journal websites (as
supplementary materials) are much better because the materials are immediately
available and likely to be findable for a longer time. Both options have some
advantages over the other and it is unclear which of the two may be better.
Thus it is recommended that authors follow both whenever possible.

Several public repositories are conveniently available and are generally being
used by a growing number of authors. For code the most popular are GitHub
(https://www.github.com) and CRAN (for R programs,
https://cran.r-project.org/), though there are several are other options. For
models the most widely used repository is the BioModels database
(https://www.ebi.ac.uk/biomodels-main/), which has the advantage of having
curators that ensure the models do indeed produce the results that are
described in the publications le_novere2006 . In the many cases when this is
not true, they contact authors and correct the issues. Thus submission of
models to BioModels already ensures a major verification of reproducibility of
methods chelliah2015 . For data of any kind (which could include programs and
models) other repositories could be used: Zenodo (https://zenodo.org/), Dryad
(https://datadryad.org//), and FigShare (https://figshare.com/), all of which
issue digital object identifiers (DOI) for data sets.

##  5 Some recommendations

Hübner et al. surveyed a sample of some 400 articles reporting research with
biomodels and found that only a minority of them properly described the
computational research performed in the study such that it could be reproduced
hubner2011 . Thus it seems that despite the readily available tools to promote
reproducibility described above, authors and journals are not applying them
widely. In order to improve the present situation, the list below includes
actions that authors should take to make their biomodel research more
reproducible. This short list is partly based on the MIRIAM proposal
le_novere2009 . In addition it is also important to consult recommendations
made for computational research in general piccolo2016 ; sandve2013 ;
stodden2016 .

  1. 1.

Whenever possible use existing peer-reviewed, actively maintained and open-
source software to create and analyze models. That way the algorithms and
their implementation have already been reviewed, are available to all readers,
and most likely will make step 3 trivial. Remember to cite the software and
mention the version number used.

  2. 2.

If the research required specially written software, deposit the code in a
public repository, or at a minimum include it as supplementary material in the
manuscript. This includes code that may require proprietary software (such as
Matlab, Comsol, etc.); your programs need to be published!

  3. 3.

Whenever possible, encode the model in an accepted standard (SBML, CellML,
etc.), include it as supplementary material, and submit it to a repository.

  4. 4.

If it is not possible to encode the model in a standard, then include the full
set of equations, parameter values, and initial conditions in the manuscript
(at least in a supplement). Make sure that all algorithms used are specified
unequivocally, as well as the software used (including version number).

  5. 5.

Publish the numerical results as data files, either in a repository, or as
supplemental files. (Note that the generic repositories mentioned above allow
very large data sets to be deposited.)

Because some authors may disregard these recommendations, either by ignorance,
for convenience of publishing quickly, or to make their research harder to
reproduce, journal editors should enforce them. In particular 2, 4 and 5 are
essential when 1 and 3 are not possible. As mentioned earlier, point 3 is the
most comprehensive and enables all three types of reproducibility. Item 2 will
only fulfill reproducibility of methods. Items 4 and 5, without any of the
others, do not guarantee reproducibility but at least describe the model and
results in detail.

##  6 Conclusion

Reproducibility is clearly an important aspect of science both for
experimental as well as computational research. Research using biomodels
should be communicated in ways that make it reproducible too. A few actions
can be taken that will greatly facilitate this objective. Scientific journals
should not publish non-reproducible research and thus should promote, or even
enforce, such actions.

## References

  * (1) Baker, M.: 1,500 scientists lift the lid on reproducibility.  Nature 533(7604), 452–454 (2016).  DOI 10.1038/533452a
  * (2) Chance, B., Garfinkel, D., Higgins, J., Hess, B.: Metabolic control mechanisms. V. A solution for the equations representing interaction between glycolysis and respiration in ascites tumor cells.  Journal of Biological Chemistry 235(8), 2426–2439 (1960)
  * (3) Chelliah, V., Juty, N., Ajmera, I., Ali, R., Dumousseau, M., Glont, M., Hucka, M., Jalowicki, G., Keating, S., Knight-Schrijver, V., Lloret-Villas, A., Natarajan, K.N., Pettit, J.B., Rodriguez, N., Schubert, M., Wimalaratne, S.M., Zhao, Y., Hermjakob, H., Le Novère, N., Laibe, C.: BioModels: ten-year anniversary.  Nucleic Acids Research 43, D542–548 (2015).  DOI 10.1093/nar/gku1181
  * (4) Claerbout, J.F., Karrenbach, M.: Electronic documents give reproducible research a new meaning.  pp. 601–604. Society of Exploration Geophysicists (1992).  DOI 10.1190/1.1822162
  * (5) Dada, J.O., Spasić, I., Paton, N.W., Mendes, P.: SBRML: a markup language for associating systems biology data with models.  Bioinformatics 26(7), 932–938 (2010).  DOI 10.1093/bioinformatics/btq069
  * (6) Gillespie, D.T.: A general method for numerically simulating the stochastic time evolution of coupled chemical reactions.  Journal of Computational Physics 22, 403–434 (1976)
  * (7) Gleeson, P., Crook, S., Cannon, R.C., Hines, M.L., Billings, G.O., Farinella, M., Morse, T.M., Davison, A.P., Ray, S., Bhalla, U.S., Barnes, S.R., Dimitrova, Y.D., Silver, R.A.: NeuroML: A Language for Describing Data Driven Models of Neurons and Networks with a High Degree of Biological Detail.  PLoS Computational Biology 6(6), e1000815 (2010).  DOI 10.1371/journal.pcbi.1000815
  * (8) Goodman, S.N., Fanelli, D., Ioannidis, J.P.A.: What does research reproducibility mean?  Science Translational Medicine 8(341), 341ps12 (2016).  DOI 10.1126/scitranslmed.aaf5027
  * (9) Greenbaum, D., Rozowsky, J., Stodden, V., Gerstein, M.: Structuring supplemental materials in support of reproducibility.  Genome Biology 18(1) (2017).  DOI 10.1186/s13059-017-1205-3
  * (10) Guerreiro, M.: Forking software used in eLife papers to GitHub (2017).  URL https://elifesciences.org/inside-elife/dbcb6949/forking-software-used-in-elife-papers-to-github
  * (11) Hoops, S., Sahle, S., Gauges, R., Lee, C., Pahle, J., Simus, N., Singhal, M., Xu, L., Mendes, P., Kummer, U.: COPASI — a COmplex PAthway SImulator.  Bioinformatics 22(24), 3067–3074 (2006).  DOI 10.1093/bioinformatics/btl485
  * (12) Hothorn, T., Held, L., Friede, T.: Biometrical Journal and Reproducible Research.  Biometrical Journal 51(4), 553–555 (2009).  DOI 10.1002/bimj.200900154
  * (13) Hothorn, T., Leisch, F.: Case studies in reproducibility.  Briefings in Bioinformatics 12(3), 288–300 (2011).  DOI 10.1093/bib/bbq084
  * (14) Hübner, K., Sahle, S., Kummer, U.: Applications and trends in systems biology in biochemistry: Systems biology in biochemical research.  FEBS Journal 278(16), 2767–2857 (2011).  DOI 10.1111/j.1742-4658.2011.08217.x
  * (15) Hucka, M., Finney, A., Sauro, H.M., Bolouri, H., Doyle, J.C., Kitano, H., Arkin, A.P., Bornstein, B.J., Bray, D., Cornish-Bowden, A., Cuellar, A.A., Dronov, S., Gilles, E.D., Ginkel, M., Gor, V., Goryanin, I.I., Hedley, W.J., Hodgman, T.C., Hofmeyr, J.H., Hunter, P.J., Juty, N.S., Kasberger, J.L., Kremling, A., Kummer, U., Le Novère, N., Loew, L.M., Lucio, D., Mendes, P., Minch, E., Mjolsness, E.D., Nakayama, Y., Nelson, M.R., Nielsen, P.F., Sakurada, T., Schaff, J.C., Shapiro, B.E., Shimizu, T.S., Spence, H.D., Stelling, J., Takahashi, K., Tomita, M., Wagner, J., Wang, J.: The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models.  Bioinformatics 19(4), 524–531 (2003).  DOI 10.1093/bioinformatics/btg015
  * (16) Le Novère, N., Bornstein, B., Broicher, A., Courtot, M., Donizelli, M., Dharuri, H., Li, L., Sauro, H., Schilstra, M., Shapiro, B., Snoep, J.L., Hucka, M.: BioModels Database: a free, centralized database of curated, published, quantitative kinetic models of biochemical and cellular systems.  Nucleic Acids Research 34(Database issue), D689–91 (2006).  DOI 10.1093/nar/gkj092
  * (17) Le Novère, N., Finney, A., Hucka, M., Bhalla, U.S., Campagne, F., Collado-Vides, J., Crampin, E.J., Halstead, M., Klipp, E., Mendes, P., Nielsen, P., Sauro, H., Shapiro, B., Snoep, J.L., Spence, H.D., Wanner, B.L.: Minimum information requested in the annotation of biochemical models (MIRIAM).  Natute Biotechnology 23(12), 1509–15 (2005).  DOI 10.1038/nbt1156
  * (18) Le Novère, N., Hucka, M., Mi, H., Moodie, S., Schreiber, F., Sorokin, A., Demir, E., Wegner, K., Aladjem, M.I., Wimalaratne, S.M., Bergman, F.T., Gauges, R., Ghazal, P., Kawaji, H., Li, L., Matsuoka, Y., Villéger, A., Boyd, S.E., Calzone, L., Courtot, M., Dogrusoz, U., Freeman, T.C., Funahashi, A., Ghosh, S., Jouraku, A., Kim, S., Kolpakov, F., Luna, A., Sahle, S., Schmidt, E., Watterson, S., Wu, G., Goryanin, I., Kell, D.B., Sander, C., Sauro, H., Snoep, J.L., Kohn, K., Kitano, H.: The Systems Biology Graphical Notation.  Nature Biotechnology 27(8), 735–741 (2009).  DOI 10.1038/nbt.1558
  * (19) Lloyd, C.M., Halstead, M.D., Nielsen, P.F.: CellML: its future, present and past.  Progress in Biophysics and Molecular Biology 85(2-3), 433–450 (2004).  DOI 10.1016/j.pbiomolbio.2004.01.004
  * (20) Loew, L., Beckett, D., Egelman, E.H., Scarlata, S.: Reproducibility of Research in Biophysics.  Biophysical Journal 108(7), E1 (2015).  DOI 10.1016/j.bpj.2015.03.002
  * (21) Maddox, J., Randi, J., Stewart, W.W.: ”High-dilution” experiments a delusion.  Nature 334(6180), 287–290 (1988).  DOI 10.1038/334287a0
  * (22) Mesirov, J.P.: Accessible Reproducible Research.  Science 327(5964), 415–416 (2010).  DOI 10.1126/science.1179653
  * (23) Moraru, I., Morgan, F., Li, Y., Loew, L., Schaff, J., Lakshminarayana, A., Slepchenko, B., Gao, F., Blinov, M.: Virtual Cell modelling and simulation software environment.  IET Systems Biology 2(5), 352–362 (2008).  DOI 10.1049/iet-syb:20080102
  * (24) Pahle, J., Challenger, J.D., Mendes, P., McKane, A.J.: Biochemical fluctuations, optimisation and the linear noise approximation.  BMC Systems Biology 6(1), 86 (2012).  DOI 10.1186/1752-0509-6-86
  * (25) Peng, R.D.: Reproducible research and Biostatistics.  Biostatistics 10(3), 405–408 (2009).  DOI 10.1093/biostatistics/kxp014
  * (26) Peng, R.D.: Reproducible Research in Computational Science.  Science 334(6060), 1226–1227 (2011).  DOI 10.1126/science.1213847
  * (27) Piccolo, S.R., Frampton, M.B.: Tools and techniques for computational reproducibility.  GigaScience 5(1) (2016).  DOI 10.1186/s13742-016-0135-4
  * (28) Plesser, H.E.: Reproducibility vs. replicability: a brief history of a confused terminology.  Frontiers in Neuroinformatics 11, 76 (2018).  DOI 10.3389/fninf.2017.00076
  * (29) Popper, K.: The logic of scientific discovery.  Hutchinson, London (1959)
  * (30) Sandve, G.K., Nekrutenko, A., Taylor, J., Hovig, E.: Ten simple rules for reproducible computational research.  PLoS Computational Biology 9(10), e1003285 (2013).  DOI 10.1371/journal.pcbi.1003285
  * (31) Stodden, V., McNutt, M., Bailey, D.H., Deelman, E., Gil, Y., Hanson, B., Heroux, M.A., Ioannidis, J.P.A., Taufer, M.: Enhancing reproducibility for computational methods.  Science 354(6317), 1240–1241 (2016).  DOI 10.1126/science.aah6168
  * (32) Stodden, V., Seiler, J., Ma, Z.: An empirical analysis of journal policy effectiveness for computational reproducibility.  Proceedings of the National Academy of Sciences 115(11), 2584–2589 (2018).  DOI 10.1073/pnas.1708290115
  * (33) Swat, M., Moodie, S., Wimalaratne, S., Kristensen, N., Lavielle, M., Mari, A., Magni, P., Smith, M., Bizzotto, R., Pasotti, L., Mezzalana, E., Comets, E., Sarr, C., Terranova, N., Blaudez, E., Chan, P., Chard, J., Chatel, K., Chenel, M., Edwards, D., Franklin, C., Giorgino, T., Glont, M., Girard, P., Grenon, P., Harling, K., Hooker, A., Kaye, R., Keizer, R., Kloft, C., Kok, J., Kokash, N., Laibe, C., Laveille, C., Lestini, G., Mentré, F., Munafo, A., Nordgren, R., Nyberg, H., Parra-Guillen, Z., Plan, E., Ribba, B., Smith, G., Trocóniz, I., Yvon, F., Milligan, P., Harnisch, L., Karlsson, M., Hermjakob, H., Le Novère, N.: Pharmacometrics Markup Language (PharmML): Opening New Perspectives for Model Exchange in Drug Development: PharmML - Pharmacometrics Markup Language.  CPT: Pharmacometrics & Systems Pharmacology 4(6), 316–319 (2015).  DOI 10.1002/psp4.57
  * (34) Waltemath, D., Adams, R., Beard, D.A., Bergmann, F.T., Bhalla, U.S., Britten, R., Chelliah, V., Cooling, M.T., Cooper, J., Crampin, E.J., Garny, A., Hoops, S., Hucka, M., Hunter, P., Klipp, E., Laibe, C., Miller, A.K., Moraru, I., Nickerson, D., Nielsen, P., Nikolski, M., Sahle, S., Sauro, H.M., Schmidt, H., Snoep, J.L., Tolle, D., Wolkenhauer, O., Le Novère, N.: Minimum information about a simulation experiment (MIASE).  PLoS Computational Biology 7(4), e1001122 (2011).  DOI 10.1371/journal.pcbi.1001122
  * (35) Waltemath, D., Adams, R., Bergmann, F.T., Hucka, M., Kolpakov, F., Miller, A.K., Moraru, I.I., Nickerson, D., Sahle, S., Snoep, J.L., Le Novère, N.: Reproducible computational biology experiments with SED-ML - The Simulation Experiment Description Markup Language.  BMC Systems Biology 5(1), 198 (2011).  DOI 10.1186/1752-0509-5-198


quadratically with time. The necessary high-precision astrometric observations
could be provided by Gaia itself. As an example, if the initial position and
proper motion of the contaminated star were known exactly and it lies 100 pc
away, an acceleration within the sky plane of $1000\,a_{{}_{0}}$ over 5 years
would cause it to appear $100\,\mu$as off from where it would be if it was
unaccelerated. In reality, the initial conditions would not be known exactly,
so astronomers would try to fit the positions with a linear trend in
accordance with Occam’s Razor. This would reduce the maximum deviation from
the unaccelerated model by a factor of 6, but ultimately a parabola can’t be
fit by a straight line. Given that DR2 of Gaia shows that it can already
measure the positions of some 15th magnitude stars to within $20\,\mu$as, it
seems likely that observers could remove a significant fraction of stars
contaminated by close binary companions.

$g_{{}_{close}}$ would fall below the Gaia detection threshold if $d$ were
much larger. Although this would reduce the resulting reflex velocity of the
contaminated star, the undetected star would still affect the total mass of
the system. For the typical parameters just discussed, this would be a rather
small effect as the system mass would be ${2.1\,M_{\odot}}$, 5% more than
expected based on the two visible Solar mass stars. This raises the Newtonian
expectation for the WB orbital velocity by 2.5%, reducing $\widetilde{v}$ by
essentially the same proportion (Equation 3). The MOND effect is very likely
much larger $-$ our calculations show that $\widetilde{v}$ can exceed the
Newtonian limit of $\sqrt{2}$ by ${\approx 20\%}$ (Section 5) for plausible
assumptions about the MOND interpolating function (Section 7.1).

A WB system could also contain planetary mass objects. If $\beta$ was reduced
ten-fold to 0.01, then the perturber would be a ${\approx 10}$ Jupiter mass
object. In this case, Equation 54 shows that the orbital separation must be
only 0.5 AU for the perturber to significantly affect $\widetilde{v}$. At this
distance, the orbital period would be short enough that the radial velocity
would change significantly over a few months. Detecting this would not require
us to know the absolute radial velocity of the star, thereby avoiding
uncertainty from the convective blueshift correction. Any time variations in
this correction would still be relevant, but these are expected to be rather
small (e.g. Kürster et al., 2003). The WBs of interest for the WBT have
orbital velocities of ${\sim 0.2}$ km/s (Figure 7), much larger than the
stellar reflex motion used to discover one of the first known exoplanets, 51
Pegasi b (Mayor & Queloz, 1995). As technology has improved greatly in the
subsequent generation, it should be straightforward to discover systems in
such a configuration. Rather than removing them from the sample, it may be
preferable to average over several orbital periods to deduce the velocity of
the star-planet system’s barycentre.

Secular accelerations decrease for companions at larger $d$. At these
distances, the companion must have a higher mass to have the same effect on
$\widetilde{v}$ (Equation 54). Thus, some initially undetected companions
might be easier to discover directly using deeper exposures of the system.
Even objects as faint as white dwarfs are routinely discovered nowadays, with
the discovery rate accelerating in the Gaia era (Kilic et al., 2018). Of
particular relevance is the recent discovery of a 5 Gyr old
${0.034\,M_{\odot}}$ object $0.535\arcsec$ from a ${1.05\,M_{\odot}}$ star
that lies 42 pc from the Sun (Cheetham et al., 2018). Their magnitude 18.5
detection in the near-infrared $J$-band suggests that, once some follow-up
observations are taken of high-$\widetilde{v}$ systems, only companions even
fainter than this could evade direct detection. Such an object would generally
be less massive, reducing its effect on $\widetilde{v}$ unless it is much
closer in (${\la 50}$ AU for a ${0.1\,M_{\odot}}$ companion). This would cause
a larger orbital acceleration (lower $\beta$ in Equation 56), likely allowing
for indirect detection in a manner analogous to exoplanets.

Our estimates suggest that a companion marginally too faint for direct
detection yet close enough to affect the WBT could be detected via secular
astrometric acceleration with ${\approx 5}$ years of Gaia observations.999We
assume a future Gaia data release will provide sky positions at different
epochs rather than just quantities derived from them. Alternatively, the data
release could include the fitted secular acceleration, which would be useful
for studies of distant companions more generally. This may be preferable to
searches for radial velocity trends as the latter only measures motion along
one dimension and requires additional follow-up. Therefore, the combination of
direct detection and secular acceleration (in proper motion or radial
velocity) should be able to definitively confirm or rule out the undetected
companion hypothesis for any high-$\widetilde{v}$ systems that may be
discovered.

####  8.2.1 The wide binary orthogonality test

Figure 11: Distribution of the angle $\psi$ between the sky-projected
separation and relative velocity of WB systems governed by MOND (red curve).
Our result is based on systems with $\widetilde{v}>0.97$ and
$r_{p}=\left(3-20\right)$ kAU, the optimal range for the WBT (Section 5).
Different orbital eccentricities are weighted according to ${\gamma=1.2}$
(Equation 22). If the relative velocity does not arise from the WB orbit but
due to e.g. an undetected close companion, we expect $\psi$ to have a uniform
distribution (black line).

If many high-$\widetilde{v}$ systems are found but these systems are
contaminated by an additional low-mass object (to evade direct detection),
then this object would have to be quite close to one of the stars forming the
WB. Thus, the extra object could be located in any direction with respect to
the star whose velocity it is contaminating. This would lead to an isotropic
relative velocity between WB stars with high $\widetilde{v}$ as
$\widetilde{v}$ would not be measuring the WB orbital motion.

However, this is not true if the high measured $\widetilde{v}$ is a genuine
consequence of the WB orbit. The systems with the highest $\widetilde{v}$
would be those observed close to pericentre, with an orbital plane aligned
closely with the sky plane. Consequently, we expect these systems to have
${\psi\approx\frac{\pi}{2}}$, where
$\psi\equiv\cos^{-1}\left(\widehat{\bm{r}}\cdot\widehat{\bm{v}}\right)$ is the
angle between the sky-projected separation and relative velocity of the stars
in a WB.101010${\psi\approx\frac{\pi}{2}}$ near apocentre as well, when
$\widetilde{v}$ would be very low. This is unexpected if the
high-$\widetilde{v}$ systems arise due to contamination. Consequently, the
distribution of $\psi$ could serve as an important check on any future claim
of detecting modified gravity using the WBT.

To see how this might work, we determine the $\psi$ distribution using our
MOND orbit library, assuming only systems with ${r_{p}\geq 3}$ kAU will be
used for the WBT and restricting ourselves to systems with sky-projected
${\widetilde{v}>0.97}$. This is because the important thing for the WBT is
verifying the orbital nature of the relative velocity in the controversial
high-$\widetilde{v}$ systems. Our results in Section 5 suggest that the
optimal $\widetilde{v}$ range for the WBT is ${0.97-1.68}$ in the more
conservative scenario where only sky-projected information is available. After
applying these restrictions on $r_{p}$ and $\widetilde{v}$, we obtain the
result shown in Figure 11.

Figure 12: The probability of detecting the difference in $\psi$ distributions
shown in Figure 11. Our results are shown against the number of systems with
$r_{p}=\left(1-20\right)$ kAU (solid red) or $r_{p}=\left(1-40\right)$ kAU
(dashed red). In both cases, we only use systems with $r_{p}>3$ kAU in the
comparison to avoid nearly Newtonian systems. Importantly, we also restrict to
systems with $\widetilde{v}>0.97$ based on our results in Section 5. Because
only a small fraction of WB systems have such a high $\widetilde{v}$ (even in
MOND), it is rather difficult to detect the seemingly obvious anisotropy in
their $\psi$ distribution. Our results show that the best way to do so is to
focus on the fraction of systems with $\psi$ between ${64\pm 2^{\circ}}$ and
${116\pm 2^{\circ}}$.

To quantify how many WB systems would be needed to detect this anisotropic
distribution in $\psi$, we use the method described in Section 4. Our results
are shown in Figure 12 for two different upper limits on $r_{p}$. The
probability distributions we compare have been normalised to the number of
systems with $r_{p}$ between 1 kAU and the indicated upper limit, though we
assume only systems with ${r_{p}>3}$ kAU are used in the comparison to avoid
nearly Newtonian systems (Section 5). Importantly, we do not restrict the
range in sky-projected $\widetilde{v}$ when normalising. This captures the
fact that only a rather small fraction of WB systems are expected to have
${\widetilde{v}>0.97}$. As a result, accurate observations of ${\approx 1000}$
systems would be required to detect the anisotropy in $\psi$ _for the high-
velocity systems_. This rather high number arises from the fact that binning
WB systems in both $\widetilde{v}$ and $\psi$ requires more systems than
binning in $\widetilde{v}$ alone, which is sufficient for the WBT.

Our results in Figure 12 assume that $\psi$ would be distributed isotropically
if the relative velocity between the stars in a WB mostly arose from
contamination by a close undetected companion to one of the stars. This is
because we expect only a weak tidal influence from the distant star in the WB.
If we suppose it is ${200\times}$ further away, then its gravity on the close
binary is $\approx 200^{2}\times$ weaker than the gravity binding the close
binary. However, the close binary is expected to be in the Newtonian regime
(Equation 56) such that a uniform external gravitational field merely moves it
without affecting its internal dynamics. Thus, only tides raised by the
distant star can be relevant.

This tide is another factor of $\frac{d}{r}$ weaker, making it $8\times
10^{6}$ times weaker than the gravity binding the close binary. It therefore
needs to complete about that many orbits to ‘notice’ the distant star. As the
orbital period would be ${\approx 350}$ years at a 50 AU separation, this
would take 2.8 Gyr. However, the WB orbit is much shorter than this. Even the
Galactic orbit is much shorter. We will see later that the WB orbital pole
precesses to maintain its angle with the Galactic EF (Figure 17). It also
precesses on shorter timescales (Figure 14). Thus, the WB orbital motion and
precession do not give the close binary enough time to significantly adjust to
the presence of its distant companion.

##  9 Proxima Centauri as a wide binary

For the WBT to work, it is important that WB systems are fairly common.
Fortunately, this seems very likely if we consider our nearest external star,
Proxima Cen. This orbits the close binary $\alpha$ Cen A and B at a current
distance of 13 kAU (Kervella et al., 2017). As first pointed out by Beech
(2009), this puts the Proxima Cen orbit well within the regime where MOND
would have a significant effect.

###  9.1 Short-term evolution

Due to its proximity, even this single system could allow a direct test of
MOND with the proposed Theia mission (Theia Collaboration, 2017). To see how
this might work, we use the method described in Section 2.1 to find the
orbital acceleration of Proxima Cen. We treat it as a test particle orbiting
the much more massive $\alpha$ Cen A and B, which we consider as a single
point mass of ${2.043\,M_{\odot}}$ at their barycentre given that they are in
a tight orbit separated by only ${\approx 18}$ AU (Kervella et al., 2016).
This is $<\frac{1}{700}$ of their distance to Proxima Cen, so the force on it
should be affected by $\approx\frac{1}{700^{2}}$ due to the finite separation
of $\alpha$ Cen A and B. This would only significantly affect the WB orbit of
Proxima Cen after $\approx 5\times 10^{5}$ orbits. However, its long orbital
period means that it has only completed $\approx 1.3\times 10^{4}$ orbits in
the past 5 Gyr, the estimated age of the system (Bazot et al., 2016).

Based on the mass-luminosity relation of Mann et al. (2015), Proxima Cen has a
mass of only ${0.122\,M_{\odot}}$, justifying our assumption that it can be
treated as a test particle. To account for the fact that the MOND gravity
between two bodies is weakened if their total mass is split more equally
(Section 7.3), we neglect the mass of Proxima Cen altogether rather than add
it to the mass of $\alpha$ Cen AB. However, when we require the Galactocentric
velocity of the system as a whole (Section 9.3), we take a weighted mean of
the values for Proxima Cen and the $\alpha$ Cen AB barycentre, thus using the
barycentre of all three stars.

Parameter |  $\alpha$ Cen A+B | Proxima Cen
---|---|---
Right ascension | $14^{h}39^{m}40.2068\arcsec$ | $14^{h}29^{m}47.7474\arcsec$
Declination | $-60^{\circ}50\arcmin 13.673\arcsec$ | $-62^{\circ}40\arcmin 52.868\arcsec$
Parallax | 747.17 mas | 768.77 mas
$\mu_{\alpha,*}$ | -3619.9 mas/yr | -3773.8 mas/yr
$\mu_{\delta}$ | 693.8 mas/yr | 770.5 mas/yr
Heliocentric | -22.332 km/s | -22.204 km/s
radial velocity
Mass | ${2.0429\,M_{\odot}}$ | ${0.1221\,M_{\odot}}$
Table 5: Observed parameters of the barycentre of the $\alpha$ Cen A+B close
binary (Kervella et al., 2016, table 1) and their companion Proxima Cen
(Kervella et al., 2017, table 2). $\mu_{\alpha,*}$ is the time derivative of
the right ascension multiplied by the cosine of the declination, making it a
true (East-West) angular velocity.

We determine the gravitational field of $\alpha$ Cen AB using the method
described in Section 2.1, assuming it feels the same $\bm{g}_{ext}$ as the Sun
because both are in much the same part of the Galaxy. We then integrate the
binary orbit forwards, starting with the initial conditions given in Table 5.
We also find the Proxima Cen trajectory in Newtonian gravity. In both cases,
any near-term observations would span a negligibly short fraction of the
${\approx 400}$ kyr orbital period, allowing us to assume a parabolic
trajectory during the observing campaign. The acceleration would be
${0.60\,a_{{}_{0}}}$ in Newtonian gravity but ${0.87\,a_{{}_{0}}}$ in MOND
(though the MOND force points ${3.1^{\circ}}$ away from the radial direction).

To see if this might be detectable, we show the angular difference between the
trajectories on our sky (Figure 13). By the time Theia is flown, we assume the
LSR parameters would be known very accurately, allowing a reliable adjustment
for the Galactic acceleration of Proxima Cen. This also has a vertical
component towards the MW disk, but this is expected to be very small because
Proxima Cen lies very close to the disk mid-plane (Ferguson et al., 2017).
These uncertainties could be reduced with accurate astrometric observations of
both Proxima Cen and $\alpha$ Cen AB.

Figure 13: Difference in the sky position of Proxima Cen depending on whether
Newtonian gravity or MOND governs its orbit about $\alpha$ Cen A and B. The
same initial conditions are used for both trajectories (Table 5). $RA_{*}$ is
the difference in right ascension multiplied by the cosine of the declination.
The total angular difference grows quadratically with time and is
$7.18\,\mu$as after a decade. Astronomers might try to fit the data by varying
the initial conditions, in which case the angular differences would be
${\approx\frac{1}{6}}$ that shown here.

Unless the initial conditions were known exactly, the difference in sky
position would actually be $\frac{1}{6}$ that shown because astronomers would
try to fit the data using different initial conditions.111111The exact ratio
will depend on spacecraft performance and other factors. We assume the fit to
data is designed to minimise its $\chi^{2}$ with respect to observations taken
at regular intervals with equal accuracy. In this case, the best linear fit to
the parabola $y=t^{2}$ over the range ${t=0-1}$ is given by $y=t-\frac{1}{6}$.
Even so, a parabola can only be fit with a straight line for so long. Thus, if
Theia is flown and achieves $\mu$as astrometric precision over a few years, it
should be able to directly measure how much Proxima Cen accelerates towards
$\alpha$ Cen AB. This would yield a much-needed strong yet direct constraint
on gravity at low accelerations.

In principle, the radial velocity $v_{r}$ of Proxima Cen could also be used to
distinguish these theories. However, a constant acceleration causes $v_{r}$ to
change linearly with time, whereas the position would respond quadratically.
Thus, $v_{r}$ would only differ by 0.5 cm/s between the models after a decade
of observations. This would be very challenging to detect, making it a much
less plausible test of MOND than using precise astrometry of Proxima Cen.

One possible complication with such tests is that an undetected exoplanet
could also cause an extra acceleration. However, as perceived at Proxima Cen,
the exoplanet is quite likely to be in a different direction than $\alpha$
Cen. If the acceleration of Proxima Cen is directed to within a few degrees of
$\alpha$ Cen but is stronger than expected in Newtonian gravity, then this
would be strong evidence against it. Moreover, a short period exoplanet would
show up in multi-epoch observations. This would not be the case for a
sufficiently long period, but in this case the greater distance implies the
exoplanet must be more massive and so more likely to be detected in other
ways. This is especially true given our proximity to the system enlarging the
angles involved, thus making it easier to efficiently suppress the light from
a star already at the faint end of the main sequence. If an anomalous
acceleration was detected, then intensive observations could be taken in its
direction from Proxima Cen.

###  9.2 Long-term evolution and orbital stability

In this section, we consider the evolution of the $\alpha$ Cen system over a
longer period. We begin by generalising the calculations of Section 9.1 to
cover a little over one full revolution, a period still short enough that the
Galactic orbit of the system is expected to have a negligible effect on the EF
it experiences. In Figure 14, we show the trajectory within the plane
orthogonal to the initial $\bm{h}$ of the system. The orbit is not closed,
with a small amount of precession evident.

Figure 14: Part of the Proxima Cen orbit relative to $\alpha$ Cen, treating
its A and B components as a single point mass at their barycentre (red * at
the origin). The $x$ and $y$ axes used here are orthogonal to each other and
to the initial orbital angular momentum $\bm{h}$. We show the orbit over the
next 909 kyr, starting with the initial conditions in Table 5. By the end of
the simulation, $\widehat{\bm{h}}$ precesses ${8.0^{\circ}}$ from its initial
orientation (not shown). We indicate the position of each pericentre (pink
squares) and apocentre (green squares). The radial period is $\approx$ 210
kyr, though the azimuthal period is slightly longer because the system rotates
${<360^{\circ}}$ between successive pericentres/apocentres.

Over much longer periods, the system rotates around the Galaxy, changing the
direction towards the Galactic Centre and thus $\widehat{\bm{g}}_{ext}$. The
magnitude of the EF might also change, something we consider in Section 9.3.
Gravity from the Galactic disk would have some effect, but this is not
expected to be significant in the Solar neighbourhood (Section 9.3.1).

To explore these issues, we integrate the WB orbit of Proxima Cen backwards
for 5 Gyr, allowing $\widehat{\bm{g}}_{ext}$ to rotate at the angular rate
$-\frac{v_{c,\odot}}{R_{\odot}}$ appropriate for the LSR. We make the
approximation that $\bm{g}_{ext}$ rotates at a constant rate and does not
change in magnitude, as appropriate for a purely circular Galactic orbit. A
more rigorously calculated Galactic orbit is considered in Section 9.3.

Our calculations show that the $\alpha$ Cen system ought to be stable over its
lifetime. This allows us to address whether just 20 revolutions is enough to
accurately determine the probability distribution of observing different
$\left(r_{p},\widetilde{v}\right)$ combinations. Figure 15 suggests that this
should be sufficient as a similar distribution is obtained when integrating
backwards over 5 Gyr instead. Moreover, when discussing the WBT, we consider a
wide range of WB orbital parameters (Table 1) rather than just the values
currently appropriate for the $\alpha$ Cen system. Even if 20 revolutions is
too short to consider how the shape and size of an individual orbit changes,
this should still get accounted for statistically when considering a
sufficiently broad range of initial conditions. Our results also indicate that
the pericentre and apocentre of the $\alpha$ Cen system typically cycle around
a closed loop in only ${\approx 8}$ radial periods (blue curve in Figure 20).

Figure 15: The probability distribution for a randomly located observer to see
different combinations of $r_{p}$ and $\widetilde{v}$ if they observe the
motion of Proxima Cen relative to $\alpha$ Cen. _Top_ : Based on integrating
over 20 revolutions. _Bottom_ : Integrating backwards for 5 Gyr, assuming the
whole system is on a circular orbit about the Galactic Centre at the speed of
the LSR. The different simulation durations give rather similar results,
suggesting that 20 orbits is long enough to accurately estimate the
statistical properties of the system.

Although the $\alpha$ Cen system is stable in our model, this may not be the
case for WB systems with different parameters. To investigate this, we explore
rotated versions of the system that preserve $r$, $v$ and $\bm{r}\cdot\bm{v}$.
The tangential velocity is rotated to explore a range of $\widehat{\bm{h}}$
along the great circle orthogonal to the present separation $\bm{r}$. To
investigate other orbital poles, we first rotate both $\bm{r}$ and $\bm{v}$
about the axis $\bm{r}\times\bm{h}$ before rotating the tangential velocity.
This lets us explore the full range of $\widehat{\bm{h}}$.

Figure 16: Blue dots show orbital poles for which a rotated version of the
$\alpha$ Cen WB orbit would have undergone a close passage ($<50$ AU) when
integrated backwards for 5 Gyr, assuming a purely circular Galactic orbit. The
observed orbital pole (large red dot) yields a stable orbit that lacks such a
close encounter. We also show red vertical gridlines at Galactic longitudes of
${90^{\circ}}$ and ${270^{\circ}}$ as these correspond to orbits for which
$\bm{h}$ is initially orthogonal to the EF direction. These orbits are
expected to be most vulnerable to a collision. As no system ever reaches a
separation exceeding 19 kAU, we infer that instability rarely arises in the
sense of a system becoming unbound by gaining energy from the time-dependent
potential.

We use Figure 16 to show the values of $\widehat{\bm{h}}$ which cause the
orbit to reach $r<50$ AU at some point in the last 5 Gyr. Although some orbits
‘crash’ in this sense, we find no orbital poles for which $r$ rises to very
large values such that the system could be considered unbound. In fact, $r$
never exceeds 19 kAU.

If full 6D phase space information becomes available for more WB systems, then
we could find whether $\widehat{\bm{h}}$ always falls in the region allowed by
MOND. Discovering systems where this is not the case would challenge the
theory. Equally, as not all directions are allowed by it, finding no such
systems would lend it some credence. This is especially true as Newtonian
gravity does not readily provide an explanation for why some WB orbital poles
should be preferred over others. The main reason is that Newtonian gravity
lacks an EFE, so only the Galactic tide can affect the internal dynamics of a
WB. Because both the WB and Galactic orbits have accelerations ${\sim
a_{{}_{0}}}$, the Galactic tide across a WB is smaller than its internal
acceleration by $\sim\frac{r}{R_{\odot}}\approx 10^{-5}$. Thus, ${\approx
10^{5}}$ WB orbits would be needed for Galactic tides to have a significant
effect. As a single orbit takes ${\approx 1}$ Myr at 10 kAU separation, there
is probably insufficient time since the Big Bang for Galactic tides to become
important. Even so, some anisotropic tidal effects may arise from the
flattened distribution of MCs and of the Galactic disk in general.

Figure 17: _Top_ : Frequency distribution of the mean angle between
$\widehat{\bm{h}}$ and $\widehat{\bm{g}}_{ext}$ as a function of its initial
value for rotated versions of the $\alpha$ Cen system integrated backwards 5
Gyr. Different initial spin vectors $\widehat{\bm{h}}$ are weighted to have an
isotropic distribution. The line of equality is shown as a black dotted line.
We use a pink square to show the result for the system whose
$\widehat{\bm{h}}$ most closely matches observations of the $\alpha$ Cen
system (Table 5). For each system, the mean is taken by considering the angle
at each pericentre and apocentre as it would be too intensive to consider the
value at every timestep. _Bottom_ : The standard deviation of these angles is
shown for each system. The blue * represents the system whose initial
conditions most closely resemble $\alpha$ Cen.

Because the Galactic orbit of $\alpha$ Cen is much slower than its WB orbit,
we expect approximate conservation of the component of $\bm{h}$ in the
$\widehat{\bm{g}}_{ext}$ direction. We might also expect the magnitude of
$\bm{h}$ to be roughly conserved, as often happens in Galactic dynamics. Thus,
we show the mean and dispersion in
$\widehat{\bm{h}}\cdot\widehat{\bm{g}}_{ext}$ for rotated versions of the
$\alpha$ Cen system (Figure 17). This reveals that
$\widehat{\bm{h}}\cdot\widehat{\bm{g}}_{ext}$ is indeed approximately
conserved despite each system rotating around the Galaxy many times. As a
result, systems with high
$\left|\widehat{\bm{h}}\cdot\widehat{\bm{g}}_{ext}\right|$ should remain
stable over the long term whereas systems with low
$\left|\widehat{\bm{h}}\cdot\widehat{\bm{g}}_{ext}\right|$ would be more
vulnerable to crashing.

It is interesting to consider whether such $\widehat{\bm{h}}$-dependent
effects could be detected using WB systems without full 6D phase space
information. One possibility is to focus on ${\bm{h}}_{{}_{LOS}}$, the line of
sight (LOS) component of ${\bm{h}}$. This should be known fairly accurately as
it requires only sky-projected positions and velocities. For a WB observed
towards the Galactic anti-centre, ${\bm{h}}_{{}_{LOS}}$ should be
approximately conserved as the LOS is aligned with $\widehat{\bm{g}}_{ext}$,
the direction about which the gravitational field should be nearly
axisymmetric. Thus, if a WB observed in this direction has a large
${\bm{h}}_{{}_{LOS}}$, it is unlikely to ever crash. This is much more
plausible if it has a very low ${\bm{h}}_{{}_{LOS}}$, perhaps leading to a
paucity of systems with low ${\bm{h}}_{{}_{LOS}}$ in the directions towards or
away from the Galactic Centre.

For WBs in the orthogonal direction on our sky but still within the Galactic
disk, the LOS is nearly orthogonal to $\widehat{\bm{g}}_{ext}$. Thus, the
${\bm{h}}_{{}_{LOS}}$ of a WB system would not be so strongly correlated with
its long-term orbital stability. This should lead to the frequency
distribution of ${\bm{h}}_{{}_{LOS}}$ varying with sky position. Such an
effect would arise even though we expect the intrinsic physical
characteristics of nearby (within ${\,{\sim}\,100}$ pc) WB systems to not
depend on viewing direction as the disk scale length is much longer (Bovy &
Rix, 2013).

Figure 18: The minimum orbital separation of the $\alpha$ Cen system during
the last 5 Gyr for rotated versions of it integrated backwards. The blue *
shows the system with $\widehat{\bm{h}}$ most closely matching the observed
system. The overdensity at very low separations arises from orbits which were
terminated once they came within 50 AU. This only occurs when
$\widehat{\bm{h}}$ has a sufficiently small component along the present EF
direction $\widehat{\bm{x}}$, a consequence of approximate conservation of
${\bm{h}\cdot\widehat{\bm{g}}_{ext}}$ due to the potential being
instantaneously axisymmetric about the very slowly varying EF direction
(Figure 17). For systems whose $\widehat{\bm{h}}_{x}$ is currently high, this
provides an angular momentum barrier that precludes very small separations.

Of course, there are numerous other complications that could arise with this
test. The main problem is that $\bm{h}_{{}_{LOS}}$ is also correlated with the
total $h$ of the system, very low values for which necessarily render a system
liable to a destructive close approach. However, this is true regardless of
the observing direction. The benefit of this test is that, when observing
along $\pm\widehat{\bm{g}}_{ext}$, a system with low $\bm{h}_{{}_{LOS}}$ is
vulnerable to a crash regardless of how large $h$ is. This is because other
components of $\bm{h}$ are not conserved, making it quite possible that they
reached very low values at some time in the history of the system.

Though we do not consider this in much detail, we use Figure 18 to show the
frequency distribution of different closest approach distances $r_{min}$ as a
function of
$\widehat{\bm{h}}\cdot\widehat{\bm{g}}_{ext}\equiv\widehat{\bm{h}}_{x}$ for
the same grid of models as are shown in Figures 16 and 17. There is a strong
correlation between $r_{min}$ and the initial value of $\widehat{\bm{h}}_{x}$
despite each system having undergone ${\ga 10^{4}}$ orbits and initially
having the same $h$. This is evident from the fact that systems which reach
very low $r_{min}$ (blue dots in Figure 16) initially had low values for
$\widehat{\bm{h}}_{x}$ (near $y$-axis in Figure 18), a consequence of the
binary orbit adiabatically adjusting to the much slower Galactic orbit.

###  9.3 Towards a full Galactic orbit

To consider the $\alpha$ Cen Galactic orbit in more detail, we need to find
its Galactic velocity by adding the heliocentric velocity of the system to the
Galactocentric velocity of the Sun, $\bm{v_{\odot}}$. In the usual Galactic
Cartesian co-ordinates (Section 2.3.1), this Solar velocity is

$\displaystyle\bm{v}_{\odot}=\begin{bmatrix}U_{\odot}\\\
V_{\odot}+v_{c,\odot}\\\ W_{\odot}\end{bmatrix}$ (57)

Here, $\left(U_{\odot},V_{\odot},W\odot\right)$ is the velocity of the Sun
with respect to the LSR, which has a speed of $v_{c,\odot}$. Our adopted
values for these parameters are given in Table 6.

Parameter | Meaning | Value
---|---|---
$R_{\odot}$ | Galactocentric radius of Sun | 8.2 kpc
$v_{c,\odot}$ | Local Standard of Rest velocity | 232.8 km/s
$U_{\odot}$ | See Equation 57 | 14.1 km/s
$V_{\odot}$ | See Equation 57 | 14.6 km/s
$W_{\odot}$ | See Equation 57 | 6.9 km/s
Table 6: Parameters governing the position and velocity of the Sun with
respect to the Galaxy. The position and velocity of the Local Standard of Rest
are obtained from McMillan (2017), which yields similar results to more recent
estimates that use Gaia data (Kawata et al., 2019). The non-circular velocity
of the Sun is obtained from Francis & Anderson (2014).

The $\alpha$ Cen system is assumed to orbit the Galactic centre in a fixed
plane, though the motion within this plane need not be circular. Although
vertical forces due to the disk must matter somewhat, we neglect these as the
Sun is ${\approx 4}$ disk scale lengths from the Galactic Centre (Bovy & Rix,
2013). We discuss the accuracy of our spherical potential approximation in
Section 9.3.1.

The velocity of $\alpha$ Cen out of the MW disk means that its orbit is tilted
with respect to the Galactic plane. We find that the tilt angle is
${4.61^{\circ}}$, causing the $\alpha$ Cen orbit to rise ${\approx 700}$ pc
out of the MW. Given that the Sun is only ${\approx 15}$ pc from the MW disk
plane (Ferguson et al., 2017, figure 6) and that $\alpha$ Cen is only 1.3 pc
from the Sun (Kervella et al., 2016), we assume that $\alpha$ Cen currently
lies within the line of nodes between its Galactic orbital plane and the MW
disk.

To estimate how the Galactocentric distance $R$ of $\alpha$ Cen varies with
time, we rotate its velocity into the Galactic disk plane. We then integrate
its orbit using the potential of the nominal MOND Galactic model from Banik &
Zhao (2018b). This shows that the mean of the perigalacticon and apogalacticon
distances is $R_{0}=8.94$ kpc, with oscillations in $R$ of amplitude $dR=0.84$
kpc. The fact that ${dR\ll R_{0}}$ means that the orbit is nearly circular,
allowing us to make an epicyclic approximation.

$\displaystyle
R~{}\approx~{}R_{0}+dR\,\sin\left(\frac{2\mathrm{\pi}t}{P_{R}}+\phi_{0}\right)$
(58)

To find the radial period $P_{R}$, we numerically advance the orbit until $R$
and $\dot{R}$ get back to their initial value, with $\dot{R}$ denoting a time
derivative. The presently observed values of $R$ and $\dot{R}$ fix the current
radial phase $\phi_{0}$ of the $\alpha$ Cen Galactic orbit. We demonstrate the
accuracy of the epicyclic approximation in Figure 19.

Figure 19: The Galactocentric radius $R$ of the $\alpha$ Cen system over one
radial period using orbit integration (solid black line) and in our epicyclic
approximation (Equation 58), shown with a dotted red line.

As $R$ varies, we also vary the tangential speed of $\alpha$ Cen around the
Galaxy to maintain angular momentum conservation. The Galactic latitude of
$\alpha$ Cen is always rather small, allowing us to make the approximation
that it varies sinusoidally with the Galactic longitude. In a spherical
potential, the maximum Galactic latitude of $\alpha$ Cen is directly
determined by its present position and velocity because the orbit is confined
to a plane.

An important objective of considering the Galactic orbit of $\alpha$ Cen in
more detail is to include the time variation of the EF strength. This is done
self-consistently with our numerical integration of the $\alpha$ Cen Galactic
orbit. As a result, the relation between the orbital separations at each
pericentre and the successive apocentre becomes less tight than for a purely
circular Galactic orbit (Figure 20). However, the result remains broadly
similar. In both cases, significant evolution of the orbit is apparent due to
the non-spherical nature of the potential. Although its time dependence must
also have some effect, the $\alpha$ Cen system takes only ${\approx 8}$ radial
periods to cycle around the curve apparent in Figure 20. This is too short for
the Galactic orbit of the system to have much effect.

Less extreme effects would be expected in the Solar System, where the stronger
$g$ makes for more nearly Newtonian behaviour. Even so, MOND combined with the
Galactic EF could conceivably explain objects like Sedna with perihelion well
outside the orbit of Neptune, the most distant Solar System planet (Brown et
al., 2004). Significant MOND effects are only possible if Sedna had a fairly
large apocentre, most likely due to interaction with a large planet. Although
such an interaction is possible in Newtonian gravity too, this would cause
Sedna to return to the heliocentric distance where it strongly interacted with
the planet. This requirement is relaxed in MOND as its Keplerian orbital
parameters are no longer fixed, perhaps providing a way to explain Sedna’s
orbit. In fact, its perihelion would take only ${\approx 10}$ Myr to double
from Neptune’s orbital radius to its observed perihelion (Paučo & Klačka,
2016). However, we note that there are alternative Newtonian explanations for
Sedna such as capture from the planetesimal disk of a passing star (Jílková et
al., 2015).

As well as increasing the perihelion distance, MOND effects can reduce it,
explaining why some orbits crash in Figure 16. This could provide an
additional mechanism by which Oort Cloud comets get perturbed onto orbits that
reach the inner Solar System. It is unclear how this would affect observations
as there are also conventional mechanisms for generating such perturbations.

Figure 20: The separation of $\alpha$ and Proxima Cen at each pericentre
against their distance at the successive apocentre, with the current orbit
shown as a large green dot. The blue dots show results for a simple model
where the whole system rotates around the Galaxy for 5 Gyr on a purely
circular orbit at the LSR speed. The red dots are based on a more realistic
Galactic orbit (Section 9.3). The points fall on a well-defined curve which is
thicker for the latter case. In both cases, the system cycles around this
curve in ${\approx 8}$ radial periods.

####  9.3.1 The effect of the Galactic disk

So far, we have not considered the fact that the MW is a disk rather than a
point mass. This is because the Sun lies several disk scale lengths out. To
test how the extended nature of the MW might affect our results, we improve
our algebraic MOND approximation in Equation 7 by including the vertical
force. Specifically, we assume that

$\displaystyle
g_{{}_{N,r}}\,\nu\left(\frac{\sqrt{{g_{{}_{N,r}}}^{2}+{g_{{}_{N,z}}}^{2}}}{a_{{}_{0}}}\right)~{}=~{}g_{r}~{}=~{}\frac{{v_{c,\odot}}^{2}}{R_{\odot}}$
(59)

This is known to work rather well in disk galaxies, at least once the vertical
Newtonian gravity $g_{{}_{N,z}}$ is taken into account when determining $\nu$
so as to more rigorously implement the algebraic MOND approximation (Jones-
Smith et al., 2018). For a thin disk, $g_{{}_{N,z}}$ may be estimated as

$\displaystyle
g_{{}_{N,z}}~{}=~{}-2\mathrm{\pi}G\Sigma\operatorname{sign}{\left(z\right)}$
(60)

We obtain the local Galactic surface density $\Sigma$ using MW parameters from
table 1 of Banik & Zhao (2018b). Using our previous estimate of
$g_{{}_{N,r}}$, we see that $g_{{}_{N,z}}\approx\frac{1}{3}g_{{}_{N,r}}$. This
means that including the vertical component of $\bm{g}_{{}_{N}}$ raises its
magnitude by $\approx\frac{1}{2\times 3^{2}}$ fractionally. As
$K_{0}\approx-0.26$, we expect that this reduces $\nu$ by only ${\approx
1.4\%}$. The lower value of $\nu$ makes $\nu g_{{}_{N,r}}$ fall below $g_{r}$,
even though these are equal when neglecting $g_{{}_{N,z}}$ (Equation 7). As a
result, $g_{{}_{N,r}}$ would have to increase ${\approx 1.4\%}$, thereby
reducing $\nu$ by ${\approx 0.35\%}$.

To check this, we solve Equation 59 using the Newton-Raphson algorithm. This
shows that $g_{{}_{N,r}}$ increases by 1.9% compared to our previous
expectation. As a result, $\nu_{ext}$ in the Solar neighbourhood falls from
1.5603 to 1.5525. Perhaps more meaningful is
$\nu_{ext}\left(1+\frac{K_{0}}{3}\right)$, the angle-averaged ratio between
Newtonian and MOND gravity (Equation 16). This decreases from 1.423 to 1.417,
implying that our results are hardly affected by the vertical component of the
Galactic gravitational field.

##  10 Conclusions

We investigated the feasibility of testing MOND using the orbital velocities
of wide binary stars whose very low orbital accelerations fall below the MOND
$a_{{}_{0}}$ threshold. Several thousand such systems have recently been
identified with a low contamination rate (Andrews et al., 2018) using Gaia
observations (Gaia Collaboration, 2018a). To keep the complexity manageable,
we solved the MOND equations assuming all the mass in each binary was in one
of the stars, though this assumption should not much affect forces in the
Solar neighbourhood (Section 7.3). We then self-consistently included the
external field (EF) from the rest of the Galaxy, leading to an axisymmetric
gravitational field. By using this to integrate systems with a range of total
masses, semi-major axes, eccentricities and orbital planes over 20
revolutions, we showed that wide binary systems have relative velocities which
can exceed Newtonian expectations by a significant amount. The excess is
${\approx 20\%}$ if we use the simple MOND interpolating function to
transition between the Newtonian and modified regimes (Famaey & Binney, 2005).
This rather gradual transition works very well at explaining a variety of
Galactic and extragalactic observations, which distinctly prefer this form
over the sharper transition of the so-called ‘standard’ function (Section
7.1). Once the interpolating function is fixed, the two major formulations of
MOND yield numerically very similar results (Table 3).

Consequently, Newtonian gravity should be distinguishable from MOND with
${\approx 500}$ well-observed systems with sky-projected separation
$r_{p}=\left(1-20\right)$ kAU, minimising contamination of the sample (Section
5). Our results show that the test is best performed by considering the
proportion of systems with $r_{p}>3$ kAU and whose scaled velocity
$\widetilde{v}$ (Equation 3) lies in the range 0.97$-$1.68.

The number of systems needed for the test could be halved if full 3D relative
velocities were available, requiring follow-up measurements of the radial
velocity accurate to within ${\approx 0.02}$ km/s. In this case, the best
$\widetilde{v}$ range is 1.05$-$1.68. Additional improvements might come from
allowing wider binaries, where MOND effects would be somewhat larger (e.g. an
upper limit of 40 kAU was used by Andrews et al., 2018). However, wider
systems would rotate slower and lead to a more contaminated sample.

It is very likely that a sufficient number of systems can be found within 150
pc (Andrews et al., 2018). At this distance, the astrometric accuracy required
is comparable to that of the most recent Gaia data release (Section 6). Using
stellar mass-luminosity relations (e.g. Mann et al., 2015), the observed
magnitudes of the stars in each binary should be sufficient to accurately
constrain their masses as their distances would be known very well (Section
6.3). It should also be feasible to obtain radial velocities of ${\approx
1000}$ nearby stars at the required accuracy (Chubak et al., 2012).

Our results include how the EF from the rest of our Galaxy weakens the self-
gravity of wide binary systems. However, some versions of MOND lack this
external field effect (EFE). Such theories predict much larger deviations from
Newtonian gravity that should become detectable with data from only ${\approx
100}$ systems, even without radial velocity measurements (Section 7.4).
Although we do not consider it likely that MOND applies without the EFE, this
would be strongly favoured by the discovery of systems with $\widetilde{v}>2$
as the EFE limits the maximum value of $\widetilde{v}$ to ${\approx 1.7}$.
Without the EFE, our results suggest that $\widetilde{v}$ could reach values
as large as 3.2. Even larger values might arise if WB systems with apocentres
beyond 100 kAU remain stable against encounters with passing stars.

We considered a number of systematic errors that could hamper our ability to
test gravity using wide binary stars (Section 8). The most serious issue would
probably be a low-mass undetected companion to at least one of the stars
(Section 8.2). There are several ways in which the effect of such companions
could be minimised. Most straightforwardly, it can cause $\widetilde{v}$ to
exceed the maximum possible in MOND. This allows such high-$\widetilde{v}$
systems to create a statistical model for the properties of low-mass binary
companions. This model could be extended down to the critical $\widetilde{v}$
range of 0.9$-$1.7 that is most sensitive to the underlying law of gravity. In
this way, we might be able to estimate the true fraction of systems whose
$\widetilde{v}$ lies in this range.

To evade direct detection, any such companion must have a rather low mass. For
such an object to noticeably affect $\widetilde{v}$, it would need to be quite
close to one of the detected stars. This would lead to a much shorter orbital
period (though still several centuries), creating a much larger acceleration
than for the wide binary orbit (Equation 56). Astrometric observations over a
few years should be able to detect this by finding a parabolic sky trajectory,
with a deviation from a linear trend of $\approx 100\,\mu$as after 5 years.
This might also be detectable as a radial velocity trend of ${\approx 4}$ m/s
per year.

Wide binary stars are not very rare $-$ the nearest star to the Sun is in just
such a system. Motivated by this, we considered the orbit of Proxima Cen about
$\alpha$ Cen in some detail (Section 9). The orbital acceleration is ${\approx
40\%}$ higher in MOND, something that could be tested with the proposed Theia
mission (Theia Collaboration, 2017) if it achieves an astrometric precision of
$\approx 1\,\mu$as over 5 years (Figure 13).

To explore the long-term stability of the orbit, we integrated it backwards
for 5 Gyr in an EF that varies with time due to the Galactic orbit of the
system. We also considered rotated versions of the Proxima Cen-$\alpha$ Cen
orbit that trace all possible orbital poles $\widehat{\bm{h}}$. Some
$\widehat{\bm{h}}$ turn out to be unstable in the sense that there was a very
close encounter at some time in the past 5 Gyr. These $\widehat{\bm{h}}$ are
nearly orthogonal to the EF direction at the present time. This arises because
the component of $\widehat{\bm{h}}$ along the EF direction is conserved in the
limit that the EF is constant with time, as the potential is axisymmetric
about the EF direction. Because changes in the EF arise only from the rather
slow Galactic orbit of the system, we expect that systems with a lot of
angular momentum along the EF direction would benefit from an angular momentum
barrier that prevents them from undergoing an extremely close encounter.

Systems like this are indeed generally stable, with the minimum separation
well correlated with the fraction of the angular momentum that presently lies
along the EF direction (Figure 18). Moreover, the angle between
$\widehat{\bm{h}}$ and the EF is conserved rather well despite the systems
completing many Galactic orbits that completely reorient the EF (Figure 17).
This suggests that we may gain a deeper analytic understanding of wide
binaries by assuming that the binary orbit adjusts adiabatically to the much
slower Galactic orbit. The dependence of orbital stability on
$\widehat{\bm{h}}$ could also provide a basis for testing MOND using novel
direction-dependent effects (Section 9.2).

Wide binary stars promise to be an exciting new frontier in the quest to
understand the law of gravity relevant to low-acceleration astrophysical
systems like galaxies. Such understanding almost certainly requires
independent non-galactic tests of the gravity law, an area that has received
little attention so far. Wide binary systems promise to provide just such a
test, perhaps using the catalogue of such systems identified by Andrews et al.
(2018) based on the Gaia mission (Perryman et al., 2001). Although Gaia has
already released data on more than a billion stars (Gaia Collaboration,
2018a), precise data on only a minute fraction of them should be sufficient to
answer fundamental questions about the laws governing our Universe.

##  Appendix A Derivation of Equation 11

In this section, we determine the contribution to the potential at the point
$\left(r,\theta\right)$ arising from phantom dark matter exterior to the
radius $r_{{}_{out}}$ (itself $\geq r$) relative to a point mass $M$ embedded
in a uniform external field $\bm{g}_{ext}$. Combining equations 24 and 28 of
Banik & Zhao (2018a), we get that

$\displaystyle\nabla^{2}\Phi~{}=~{}-\frac{GM\nu_{{}_{ext}}K_{0}}{r^{3}}\left(3\cos^{2}\theta-1\right)$
(61)

We now consider the phantom dark matter within a thin shell of radius
$r_{{}_{sh}}$ and thickness $dr_{{}_{sh}}$. Because $\nabla^{2}\Phi=0$
everywhere outside the thin shell and its surface density has the same angular
dependence as the second Legendre polynomial, the solution must be
$\Phi\propto r^{n}\left(3\cos^{2}\theta-1\right)$. The Laplacian of this is 0
if and only if $n=2$ or $-3$. The choice is fixed by the need to avoid the
potential diverging at small and large radii.

$\displaystyle\Phi=\left\\{\begin{array}[]{ll}u\left(3\cos^{2}\theta-1\right)\left(\frac{r}{r_{{}_{sh}}}\right)^{2}&\mbox{if
}r\leq r_{{}_{sh}}\\\
u\left(3\cos^{2}\theta-1\right)\left(\frac{r}{r_{{}_{sh}}}\right)^{-3}&\mbox{if
}r\geq r_{{}_{sh}}\end{array}\right.$ (64)

To find the multiplicative constant $u$, we need to consider the discontinuity
in $\frac{\partial\Phi}{\partial r}$ across the shell at $r=r_{{}_{sh}}$. This
must equal the surface density of phantom dark matter in the shell (Equation
61).

$\displaystyle-5u\left(3\cos^{2}\theta-1\right)~{}$ $\displaystyle=$
$\displaystyle~{}-\frac{GM\nu_{{}_{ext}}K_{0}}{{r_{{}_{sh}}}^{2}}\left(3\cos^{2}\theta-1\right)dr_{{}_{sh}}$
$\displaystyle u~{}$ $\displaystyle=$
$\displaystyle~{}\frac{GM\nu_{{}_{ext}}K_{0}}{5{r_{{}_{sh}}}^{2}}dr_{{}_{sh}}$
(65)

This lets us find the contribution to $\Phi$ from shells with $r_{{}_{sh}}\geq
r_{{}_{out}}$.

$\displaystyle\Delta\Phi~{}$ $\displaystyle=$
$\displaystyle~{}\int_{r_{{}_{out}}}^{\infty}\frac{GM\nu_{{}_{ext}}K_{0}\left(3\cos^{2}\theta-1\right)}{5{r_{{}_{sh}}}^{2}}\left(\frac{r}{r_{{}_{sh}}}\right)^{2}dr_{{}_{sh}}$
$\displaystyle=$
$\displaystyle~{}\frac{GM\nu_{{}_{ext}}K_{0}r^{2}\left(3\cos^{2}\theta-1\right)}{15{r_{{}_{out}}}^{3}}$
(66)

Differentiating this reproduces our result in Equation 11, something that is
most easily verified using Cartesian co-ordinates. We can perform a similar
calculation to get the contribution to $\Phi$ from all shells with
$r_{{}_{sh}}\leq r_{{}_{in}}$ for some $r_{{}_{in}}\leq r$.

$\displaystyle\Delta\Phi~{}$ $\displaystyle=$
$\displaystyle~{}\int_{0}^{r_{{}_{in}}}\frac{GM\nu_{{}_{ext}}K_{0}}{5{r_{{}_{sh}}}^{2}}\left(\frac{r}{r_{{}_{sh}}}\right)^{-3}dr_{{}_{sh}}$
(67) $\displaystyle=$
$\displaystyle~{}\frac{GM\nu_{{}_{ext}}K_{0}\left(3\cos^{2}\theta-1\right){r_{{}_{in}}}^{2}}{10r^{3}}$
(68)

To find the total potential at some point $\left(r,\theta\right)$, we can
combine the results of Equations 66 and 68 evaluated for
$r_{{}_{in}}=r_{{}_{out}}=r$.

$\displaystyle\Phi~{}$ $\displaystyle=$
$\displaystyle~{}\frac{GM\nu_{{}_{ext}}K_{0}\left(3\cos^{2}\theta-1\right)}{r}\left(\frac{1}{15}+\frac{1}{10}\right)$
(69) $\displaystyle=$
$\displaystyle~{}\frac{GM\nu_{{}_{ext}}K_{0}\left(3\cos^{2}\theta-1\right)}{6r}$
(70)

This reproduces equation 29 of Banik & Zhao (2018a), which was obtained using
equation 2.95 of Binney & Tremaine (2008). Another consistency check is that
our solution for $\Phi$ satisfies Equation 61.

## Acknowledgements

IB was partially supported by Science and Technology Facilities Council
studentship 1506672. He is grateful to the Weizmann Institute of Science for
hosting him during a week-long visit. During this visit, he had some very
helpful discussions with Mordehai Milgrom, Boaz Katz and Doron Kushnir that
greatly advanced this contribution. He also thanks Will Sutherland for
providing an observational perspective and Duncan Forgan for comments on
tides. The authors wish to thank the referee for comments which improved this
manuscript.

## References

  * Alcock et al. (2000) Alcock C., et al., 2000, ApJ, 542, 281
  * Andrews et al. (2017) Andrews J. J., Chanamé J., Agüeros M. A., 2017, MNRAS, 472, 675
  * Andrews et al. (2018) Andrews J. J., Chanamé J., Agüeros M. A., 2018, Research Notes of the American Astronomical Society, 2, 29
  * Bahcall et al. (1985) Bahcall J. N., Hut P., Tremaine S., 1985, ApJ, 290, 15
  * Banik & Zhao (2016) Banik I., Zhao H., 2016, MNRAS, 459, 2237
  * Banik & Zhao (2017) Banik I., Zhao H., 2017, MNRAS, 467, 2180
  * Banik & Zhao (2018a) Banik I., Zhao H., 2018a, SciFed Journal of Astrophysics, 1, 1000008
  * Banik & Zhao (2018b) Banik I., Zhao H., 2018b, MNRAS, 473, 419
  * Banik & Zhao (2018c) Banik I., Zhao H., 2018c, MNRAS, 473, 4033
  * Banik et al. (2018) Banik I., O’Ryan D., Zhao H., 2018, MNRAS, 477, 4768
  * Barnes & Hernquist (1992) Barnes J. E., Hernquist L., 1992, Nature, 360, 715
  * Bazot et al. (2016) Bazot M., Christensen-Dalsgaard J., Gizon L., Benomar O., 2016, MNRAS, 460, 1254
  * Beech (2009) Beech M., 2009, MNRAS, 399, L21
  * Beech (2011) Beech M., 2011, ApSS, 333, 419
  * Bekenstein & Milgrom (1984) Bekenstein J., Milgrom M., 1984, ApJ, 286, 7
  * Bell & de Jong (2001) Bell E. F., de Jong R. S., 2001, ApJ, 550, 212
  * Belloni et al. (2017) Belloni D., Askar A., Giersz M., Kroupa P., Rocha-Pinto H. J., 2017, MNRAS, 471, 2812
  * Bellstedt et al. (2018) Bellstedt S., et al., 2018, MNRAS, 476, 4543
  * Belokurov et al. (2018) Belokurov V., Erkal D., Evans N. W., Koposov S. E., Deason A. J., 2018, MNRAS, 478, 611
  * Bershady et al. (2010) Bershady M. A., Verheijen M. A. W., Swaters R. A., Andersen D. R., Westfall K. B., Martinsson T., 2010, ApJ, 716, 198
  * Bílek et al. (2018) Bílek M., Thies I., Kroupa P., Famaey B., 2018, A&A, 614, A59
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press
  * Bovy (2017) Bovy J., 2017, MNRAS, 470, 1360
  * Bovy & Rix (2013) Bovy J., Rix H.-W., 2013, ApJ, 779, 115
  * Brada & Milgrom (1999) Brada R., Milgrom M., 1999, ApJ, 519, 590
  * Brown et al. (2004) Brown M. E., Trujillo C., Rabinowitz D., 2004, ApJ, 617, 645
  * Caldwell et al. (2017) Caldwell N., et al., 2017, ApJ, 839, 20
  * Chae et al. (2017) Chae K.-H., Bernardi M., Sheth R. K., 2017, preprint, Arxiv (arXiv:1707.08280)
  * Cheetham et al. (2018) Cheetham A., et al., 2018, A&A, 614, A16
  * Chubak et al. (2012) Chubak C., Marcy G., Fischer D. A., Howard A. W., Isaacson H., Johnson J. A., Wright J. T., 2012, preprint, Arxiv (arXiv:1207.6212)
  * Claydon et al. (2017) Claydon I., Gieles M., Zocchi A., 2017, MNRAS, 466, 3937
  * D’Onghia (2015) D’Onghia E., 2015, ApJL, 808, L8
  * Desmond (2017a) Desmond H., 2017a, MNRAS, 464, 4160
  * Desmond (2017b) Desmond H., 2017b, MNRAS, 472, L35
  * Dubinski et al. (2008) Dubinski J., Gauthier J.-R., Widrow L., Nickerson S., 2008, in Funes J. G., Corsini E. M., eds, Astronomical Society of the Pacific Conference Series Vol. 396, Formation and Evolution of Galaxy Disks. p. 321 (arXiv:0802.3997)
  * Durazo et al. (2017) Durazo R., Hernandez X., Cervantes Sodi B., Sánchez S. F., 2017, ApJ, 837, 179
  * Eker et al. (2015) Eker Z., et al., 2015, AJ, 149, 131
  * Fall & Lynden-Bell (1981) Fall S. M., Lynden-Bell D., 1981, The structure and evolution of normal galaxies. NATO Advanced Study Institute held at the Institute of Astronomy and Clare College, Cambridge, UK, 3 - 15 August 1980.. Cambridge University Press
  * Famaey & Binney (2005) Famaey B., Binney J., 2005, MNRAS, 363, 603
  * Famaey & McGaugh (2012) Famaey B., McGaugh S. S., 2012, Living Reviews in Relativity, 15, 10
  * Famaey et al. (2018) Famaey B., McGaugh S., Milgrom M., 2018, MNRAS, forthcoming, forthcoming
  * Farr et al. (2011) Farr W. M., Sravan N., Cantrell A., Kreidberg L., Bailyn C. D., Mandel I., Kalogera V., 2011, ApJ, 741, 103
  * Feng & Jones (2018) Feng F., Jones H. R. A., 2018, MNRAS, 473, 3185
  * Ferguson et al. (2017) Ferguson D., Gardner S., Yanny B., 2017, ApJ, 843, 141
  * Forero-Romero & Arias (2018) Forero-Romero J. E., Arias V., 2018, MNRAS, 478, 5538
  * Francis & Anderson (2014) Francis C., Anderson E., 2014, Celestial Mechanics and Dynamical Astronomy, 118, 399
  * Fuchs (2003) Fuchs B., 2003, ApSS, 284, 719
  * Gaia Collaboration (2018a) Gaia Collaboration 2018a, A&A, 616, A1
  * Gaia Collaboration (2018b) Gaia Collaboration 2018b, A&A, 616, A11
  * Gentile et al. (2011) Gentile G., Famaey B., de Blok W. J. G., 2011, A&A, 527, A76
  * Gilmore & Reid (1983) Gilmore G., Reid N., 1983, MNRAS, 202, 1025
  * Güdel et al. (2007) Güdel M., et al., 2007, A&A, 468, 353
  * Haghi et al. (2016) Haghi H., Bazkiaei A. E., Zonoozi A. H., Kroupa P., 2016, MNRAS, 458, 4172
  * Hees et al. (2016) Hees A., Famaey B., Angus G. W., Gentile G., 2016, MNRAS, 455, 449
  * Hernandez et al. (2012) Hernandez X., Jiménez M. A., Allen C., 2012, European Physical Journal C, 72, 1884
  * Hernandez et al. (2013) Hernandez X., Jiménez M. A., Allen C., 2013, MNRAS, 428, 3196
  * Hohl (1971) Hohl F., 1971, ApJ, 168, 343
  * Hu & Sijacki (2018) Hu S., Sijacki D., 2018, MNRAS, 478, 1576
  * Ibata et al. (2013) Ibata R. A., et al., 2013, Nature, 493, 62
  * Iocco et al. (2015) Iocco F., Pato M., Bertone G., 2015, Physical Review D, 92, 084046
  * Jiang & Tremaine (2010) Jiang Y.-F., Tremaine S., 2010, MNRAS, 401, 977
  * Jílková et al. (2015) Jílková L., Portegies Zwart S., Pijloo T., Hammer M., 2015, MNRAS, 453, 3157
  * Jones-Smith et al. (2018) Jones-Smith K., Abraham R., Kell L., Mathur H., 2018, New Journal of Physics, 20, 063042
  * Kahn & Woltjer (1959) Kahn F. D., Woltjer L., 1959, ApJ, 130, 705
  * Kaib & Raymond (2014) Kaib N. A., Raymond S. N., 2014, ApJ, 782, 60
  * Kawata et al. (2019) Kawata D., Bovy J., Matsunaga N., Baba J., 2019, MNRAS, 482, 40
  * Kent (1987) Kent S. M., 1987, AJ, 93, 816
  * Kervella et al. (2016) Kervella P., Mignard F., Mérand A., Thévenin F., 2016, A&A, 594, A107
  * Kervella et al. (2017) Kervella P., Thévenin F., Lovis C., 2017, A&A, 598, L7
  * Kilic et al. (2018) Kilic M., Hambly N. C., Bergeron P., Genest-Beaulieu C., Rowell N., 2018, MNRAS, 479, L113
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kroupa et al. (2005) Kroupa P., Theis C., Boily C. M., 2005, A&A, 431, 517
  * Küpper et al. (2010) Küpper A. H. W., Kroupa P., Baumgardt H., Heggie D. C., 2010, MNRAS, 407, 2241
  * Kürster et al. (2003) Kürster M., et al., 2003, A&A, 403, 1077
  * Lee (2018) Lee J., 2018, ApJ, 856, 57
  * Lee & Li (2017) Lee J., Li B., 2017, ApJ, 842, 2
  * Lee et al. (2015) Lee J., Kim S., Rey S.-C., 2015, ApJ, 815, 43
  * Lelli et al. (2016) Lelli F., McGaugh S. S., Schombert J. M., 2016, AJ, 152, 157
  * Lelli et al. (2017) Lelli F., McGaugh S. S., Schombert J. M., Pawlowski M. S., 2017, ApJ, 836, 152
  * Lépine & Bongiorno (2007) Lépine S., Bongiorno B., 2007, AJ, 133, 889
  * Li et al. (2018) Li P., Lelli F., McGaugh S., Schombert J., 2018, A&A, 615, A3
  * Lin & Shu (1964) Lin C. C., Shu F. H., 1964, ApJ, 140, 646
  * Mann et al. (2015) Mann A. W., Feiden G. A., Gaidos E., Boyajian T., von Braun K., 2015, ApJ, 804, 64
  * Martin et al. (2018) Martin N. F., Collins M. L. M., Longeard N., Tollerud E., 2018, ApJL, 859, L5
  * Mayor & Queloz (1995) Mayor M., Queloz D., 1995, Nature, 378, 355
  * McGaugh (2011) McGaugh S. S., 2011, Physical Review Letters, 106, 121303
  * McGaugh (2016a) McGaugh S. S., 2016a, ApJ, 816, 42
  * McGaugh (2016b) McGaugh S. S., 2016b, ApJL, 832, L8
  * McGaugh & Milgrom (2013a) McGaugh S., Milgrom M., 2013a, ApJ, 775, 139
  * McGaugh & Milgrom (2013b) McGaugh S., Milgrom M., 2013b, ApJ, 775, 139
  * McGaugh & Wolf (2010) McGaugh S. S., Wolf J., 2010, ApJ, 722, 248
  * McGaugh et al. (1995) McGaugh S. S., Schombert J. M., Bothun G. D., 1995, AJ, 109, 2019
  * McGaugh et al. (2016) McGaugh S., Lelli F., Schombert J., 2016, Phys. Rev. Lett., 117, 201101
  * McMillan (2017) McMillan P. J., 2017, MNRAS, 465, 76
  * Milgrom (1983) Milgrom M., 1983, ApJ, 270, 365
  * Milgrom (1986) Milgrom M., 1986, ApJ, 302, 617
  * Milgrom (1989) Milgrom M., 1989, ApJ, 338, 121
  * Milgrom (1999) Milgrom M., 1999, Phys. Lett. A, 253, 273
  * Milgrom (2010) Milgrom M., 2010, MNRAS, 403, 886
  * Milgrom (2011) Milgrom M., 2011, Acta Physica Polonica B, 42, 2175
  * Mirabel et al. (1992) Mirabel I. F., Dottori H., Lutz D., 1992, A&A, 256, L19
  * Müller et al. (2018) Müller O., Pawlowski M. S., Jerjen H., Lelli F., 2018, Science, 359, 534
  * Nakanishi & Sofue (2006) Nakanishi H., Sofue Y., 2006, PASJ, 58, 847
  * Norris et al. (2016) Norris M. A., et al., 2016, ApJ, 832, 198
  * Ostriker & Peebles (1973) Ostriker J. P., Peebles P. J. E., 1973, ApJ, 186, 467
  * Ostriker & Steinhardt (1995) Ostriker J. P., Steinhardt P. J., 1995, Nature, 377, 600
  * Paučo & Klačka (2017) Paučo R., Klačka J., 2017, preprint, Arxiv (arXiv:1705.09273)
  * Paučo & Klačka (2016) Paučo R., Klačka J., 2016, A&A, 589
  * Pawlowski (2018) Pawlowski M. S., 2018, Modern Physics Letters A, 33, 1830004
  * Pazy (2013) Pazy E., 2013, Phys. Rev. D, 87, 084063
  * Peebles (2017a) Peebles P. J. E., 2017a, Nature Astronomy, 1, 0057
  * Peebles (2017b) Peebles P. J. E., 2017b, preprint, Arxiv (arXiv:1705.10683)
  * Perryman et al. (2001) Perryman M. A. C., et al., 2001, A&A, 369, 339
  * Peters (1981) Peters P. C., 1981, American Journal of Physics, 49, 564
  * Peters & Kuzio de Naray (2018) Peters W., Kuzio de Naray R., 2018, MNRAS, 476, 2938
  * Pittordis & Sutherland (2018) Pittordis C., Sutherland W., 2018, MNRAS, 480, 1778
  * Reid et al. (2006) Reid I. N., Lewitus E., Allen P. R., Cruz K. L., Burgasser A. J., 2006, AJ, 132, 891
  * Riess et al. (1998) Riess A. G., et al., 1998, AJ, 116, 1009
  * Rogstad & Shostak (1972) Rogstad D. H., Shostak G. S., 1972, ApJ, 176, 315
  * Rowland (2015) Rowland D. R., 2015, International Journal of Modern Physics D, 24, 1550065
  * Saburova (2011) Saburova A. S., 2011, Astronomy Reports, 55, 409
  * Salucci & Turini (2017) Salucci P., Turini N., 2017, preprint, Arxiv (arXiv:1707.01059)
  * Scarpa et al. (2017) Scarpa R., Ottolina R., Falomo R., Treves A., 2017, International Journal of Modern Physics D, 26, 1750067
  * Schutz et al. (2018) Schutz K., Lin T., Safdi B. R., Wu C.-L., 2018, Physical Review Letters, 121, 081101
  * Shao et al. (2018) Shao S., Cautun M., Frenk C. S., Grand R. J. J., Gómez F. A., Marinacci F., Simpson C. M., 2018, MNRAS, 476, 1796
  * Shaya & Olling (2011) Shaya E. J., Olling R. P., 2011, ApJS, 192, 2
  * Smolin (2017) Smolin L., 2017, Physical Review D, 96, 083523
  * Spada et al. (2013) Spada F., Demarque P., Kim Y.-C., Sills A., 2013, ApJ, 776, 87
  * Theia Collaboration (2017) Theia Collaboration 2017, preprint, Arxiv (arXiv:1707.01348)
  * Tisserand et al. (2007) Tisserand P., et al., 2007, A&A, 469, 387
  * Tokovinin (2017) Tokovinin A., 2017, MNRAS, 468, 3461
  * Tokovinin & Kiyaeva (2016) Tokovinin A., Kiyaeva O., 2016, MNRAS, 456, 2070
  * Toomre (1964) Toomre A., 1964, ApJ, 139, 1217
  * Trifonov et al. (2018) Trifonov T., et al., 2018, A&A, 609, A117
  * Vallée (2014) Vallée J. P., 2014, AJ, 148, 5
  * Vallée (2017) Vallée J. P., 2017, ApSS, 362, 79
  * Verlinde (2016) Verlinde E. P., 2016, SciPost Physics, 2, 16
  * Wetzstein et al. (2007) Wetzstein M., Naab T., Burkert A., 2007, MNRAS, 375, 805
  * Williams et al. (2017) Williams A. A., Belokurov V., Casey A. R., Evans N. W., 2017, MNRAS, 468, 2359
  * Yu & Liu (2018) Yu J., Liu C., 2018, MNRAS, 475, 1093
  * Zhao et al. (2010) Zhao H., Li B., Bienaymé O., 2010, Physical Review D, 82, 103001
  * Zhao et al. (2013) Zhao H., Famaey B., Lüghausen F., Kroupa P., 2013, A&A, 557, L3
  * Zwicky (1937) Zwicky F., 1937, ApJ, 86, 217
  * Zwitter (2003) Zwitter T., 2003, in Munari U., ed., Astronomical Society of the Pacific Conference Series Vol. 298, GAIA Spectroscopy: Science and Technology. p. 329 (arXiv:astro-ph/0301213)
  * de Almeida et al. (2016) de Almeida Á. O. F., Piattella O. F., Rodrigues D. C., 2016, MNRAS, 462, 2706
  * den Heijer et al. (2015) den Heijer M., et al., 2015, A&A, 581, A98
  * van Dokkum et al. (2018a) van Dokkum P., et al., 2018a, RNAAS, 2, 54
  * van Dokkum et al. (2018b) van Dokkum P., et al., 2018b, Nature, 555, 629
  * van Leeuwen (2007) van Leeuwen F., 2007, A&A, 474, 653